 # How to compute solution with ALE moving boundaries?

Hi everyone,

I tried to draw four cells inside the box like the image, and tried to move it by using ALE. Two right side cells are moving +1 for x-axis, and the other two cells are moving -1 for x-axis. When I try to do this, I found that four cells are successfully moved. The numerical solution u, however, doesn’t look correct because there are hardly value differences (I mean, all u values are bewteen 0.49 and 0.51) and near the lower left cell’s values are not very big (It supposed to have greater values than the other cells since q_0 is much greater than q_1.)

This is initial mesh.

First step’s numerical solution

Second step’s solution

Third step

Forth step

Fifth step

So I am wondering if there’s any problem with my code.
What am I missing? Can you guys please tell me how to compute numerical solution with new moved boundary by ALE?

In my opinion, some suspicious parts are

``````boundary_marker = MeshFunction('size_t', mesh, 1, mesh.domains())
boundary_marker.set_all(0)
``````

this part,

``````all = all_boundaries()
all.mark(boundary_marker,1)
inc = on_inclusion()
inc.mark(boundary_marker,2)
circle1 = Circles1()
circle1.mark(boundary_marker,3)
circle2 = Circles2()
circle2.mark(boundary_marker,4)
circle3 = Circles3()
circle3.mark(boundary_marker,5)
circle4 = Circles4()
circle4.mark(boundary_marker,6)
``````

this part

``````for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
q, s = boundary_conditions[i]['Robin']  # q represents q_0, q_1, ... , s represents s_0 and s_1, which are the multiplication of k_m and R_th(R_nk).
integrals_R.append(q * v * ds(i+3))
integrals_L.append(s * u * v * ds(i+3))
``````

and this part as well.

Thanks a lot for your help. Please teach me how to compute numerical solutions with moving boundary by ALE.

``````from mshr import *
import numpy as np
from dolfin import *
import matplotlib.pyplot as plt

box_l = 90  # Box length
box_m = 10  # Box margin
c_r = 5  # Cell radius
c_d = 10  # Cell to cell distance

# Global domain
rect = Rectangle( Point(0.0 , 0.0) , Point(box_l , box_l) )
circ1 = Circle( Point(15,15), c_r, 30)
circ2 = Circle( Point(75,15), c_r, 30)
circ3 = Circle( Point(15,75), c_r, 30)
circ4 = Circle( Point(75,75), c_r, 30)
domain = rect - circ1 - circ2 - circ3 - circ4
mesh = generate_mesh(domain,30)

plot(mesh)
plt.show()

V = VectorFunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)

boundary_marker = MeshFunction('size_t', mesh, 1, mesh.domains())
boundary_marker.set_all(0)

# Define boundary conditions and get boundary marker
class on_inclusion(SubDomain):
def inside(self, x, on_boundary):
return (on_boundary and not(near(x, 0) or near(x, box_l) or near(x, 0) or near(x, box_l)))

class all_boundaries(SubDomain):
def inside(self, x, on_boundary):
return (on_boundary)

class Circles1(SubDomain):
tol = 0.5
def inside(self, x, on_boundary):
return on_boundary and near( (x - 15) ** 2 + (x - 15) ** 2, (c_r) ** 2,self.tol)

class Circles2(SubDomain):
tol = 0.5
def inside(self, x, on_boundary):
return on_boundary and near( (x - 75) ** 2 + (x - 15) ** 2, (c_r) ** 2,self.tol)

class Circles3(SubDomain):
tol = 0.5
def inside(self, x, on_boundary):
return on_boundary and near( (x - 15) ** 2 + (x - 75) ** 2, (c_r) ** 2,self.tol)

class Circles4(SubDomain):
tol = 0.5
def inside(self, x, on_boundary):
return on_boundary and near( (x - 75) ** 2 + (x - 75) ** 2, (c_r) ** 2,self.tol)

all = all_boundaries()
all.mark(boundary_marker,1)
inc = on_inclusion()
inc.mark(boundary_marker,2)
circle1 = Circles1()
circle1.mark(boundary_marker,3)
circle2 = Circles2()
circle2.mark(boundary_marker,4)
circle3 = Circles3()
circle3.mark(boundary_marker,5)
circle4 = Circles4()
circle4.mark(boundary_marker,6)

center_coordinate = [ [15,15] , [75,15] , [15,75] , [75,75] ]

####################################################
### Solving part
####################################################

def linear_solver(kappa,f,boundary_conditions,degree=1,subdomains=[],linear_solver="Krylov",abs_tol=1E-5,rel_tol=1E-3,max_iter=1000):
W = FunctionSpace(mesh, 'P', degree)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_marker)

bcs_new = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_marker, i)
bcs_new.append(bc)

u = TrialFunction(W)
v = TestFunction(W)

# Collect Robin integrals
integrals_R = []
integrals_L = []

for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
q, s = boundary_conditions[i]['Robin']  # q represents q_0, q_1, ... , s represents s_0 and s_1, which are the multiplication of k_m and R_th(R_nk).
integrals_R.append(q * v * ds(i))
integrals_L.append(s * u * v * ds(i))

# Define LHS and RHS of weak form
L = sum(integrals_R)
u = Function(W)

problem = LinearVariationalProblem(a, L, u, bcs_new)
solver = LinearVariationalSolver(problem)

solver.parameters["linear_solver"] = 'gmres'
solver.parameters["preconditioner"] = 'ilu'
solver.parameters["krylov_solver"]["absolute_tolerance"] = abs_tol
solver.parameters["krylov_solver"]["relative_tolerance"] = rel_tol
solver.parameters["krylov_solver"]["maximum_iterations"] = max_iter

solver.solve()

plot(u)
plt.title("Numerical Solution")
plt.show()

return u

####################################################################
###  Setting Part
####################################################################

import sympy as sym

f = 0.0
q_0 = 25
q_1 = 1.3
s_0 = 2 * q_0
s_1 = 2 * q_1

data_set = [f,q_0,s_0,q_1,s_1]
data_set = [sym.printing.ccode(var) for var in data_set]
data_set = [Expression(var,degree=2) for var in data_set]

f,q_0,s_0,q_1,s_1 = data_set

boundary_conditions = { 0: {'Robin':(q_0,s_0)} , 1: {'Robin':(q_1,s_1)} , 2: {'Robin':(q_1,s_1)} , 3: {'Robin':(q_1,s_1)} }

kappa = Constant(1)  # This is corresponding to the diffusion coefficient.

# Five times moving loop
for mm in range(5):

# Figure out solution u before using ALE
u = linear_solver(kappa, f, boundary_conditions, degree=1, linear_solver="direct")
vtkfile = File("paraview_ex/ttttt.pvd")
vtkfile << u

# Prepare for cell's moving
bcs = [DirichletBC(V, Constant((0.0, 0.0)), boundary_marker, 1), DirichletBC(V, Constant((0.0, 0.0)), boundary_marker, 2),
DirichletBC(V, Constant(( 1 , 0 )), boundary_marker, 3), DirichletBC(V, Constant(( -1 , 0 )), boundary_marker, 4),
DirichletBC(V, Constant(( 1 , 0 )), boundary_marker, 5), DirichletBC(V, Constant(( -1 , 0 )), boundary_marker, 6) ]

disp_bc1 = bcs
disp_bc2 = bcs
disp_bc3 = bcs
disp_bc4 = bcs
disp_bc5 = bcs
disp_bc6 = bcs

b_disp = Function(V)

disp_bc1.apply(b_disp.vector())
disp_bc2.apply(b_disp.vector())
disp_bc3.apply(b_disp.vector())
disp_bc4.apply(b_disp.vector())
disp_bc5.apply(b_disp.vector())
disp_bc6.apply(b_disp.vector())

# Move the cells
new_mesh = Mesh(mesh)
new_boundaries = MeshFunction("size_t", new_mesh, 1)
new_boundaries.set_values(boundary_marker.array())
ALE.move(mesh, b_disp)

# mesh optimization
for j in range(380):
mesh.smooth()
plot(mesh)
plt.show()``````