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.
For your information, I uploaded full code below.
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[1], 0) or near(x[1], box_l) or near(x[0], 0) or near(x[0], 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[0] - 15) ** 2 + (x[1] - 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[0] - 75) ** 2 + (x[1] - 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[0] - 15) ** 2 + (x[1] - 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[0] - 75) ** 2 + (x[1] - 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
a = kappa * dot(grad(u), grad(v)) * dx + sum(integrals_L)
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[0]
disp_bc2 = bcs[1]
disp_bc3 = bcs[2]
disp_bc4 = bcs[3]
disp_bc5 = bcs[4]
disp_bc6 = bcs[5]
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()