Hi guys,
As a very beginner of fenics, I tried to make this code where there are 4 circles inside the box. The right bottom one has a intense boundary condition, while the other three circles have very little value of this, as you can see the image below.
Here is my question.
My professor asked me to compute moving boundary condition FEM. He asked me that if it is possible for us to move the cell (it means boundaries are moving) without changing mesh grid. As far as I know, there’s ALE which can be used for moving boundary problems. Can I do this with ALE?
If it is possible, how can I make this code? Does anyone show me a very simple example code?
I really appreciate this website because there are so many experts on fenics!
This is the code that I made.
from mshr import *
import numpy as np
from dolfin import *
import matplotlib.pyplot as plt
# Number of cells for one line
# Users can change this number
# Margin, cell distance, cell radius are decided by this.
number_of_cells_in_line = 2
# numbering list for boundary marker
numbering = []
tt=0
for i in range(number_of_cells_in_line):
numbering.append([])
for j in range(number_of_cells_in_line):
numbering[i].append(tt)
tt=tt+1
# Coordinate setting
# Each cell's center coordinates information
center_coordinates = []
for i in range(number_of_cells_in_line):
center_coordinates.append([])
for j in range(number_of_cells_in_line):
center_coordinates[i].append( [(3*j+2) / (3*(number_of_cells_in_line)+1) , (3*i+2) / (3*(number_of_cells_in_line)+1) ] )
####################################################
### Solving part
####################################################
def solver_bcs(kappa,f,boundary_conditions,degree=1,subdomains=[],linear_solver="Krylov",abs_tol=1E-5,rel_tol=1E-3,max_iter=1000):
# Global domain
rect = Rectangle( Point(0.0 , 0.0) , Point(1.0 , 1.0) )
# Collect cells area
circ_collection=[]
# Draw each cells in domain
for i in range(number_of_cells_in_line):
for j in range(number_of_cells_in_line):
circ = Circle( Point(center_coordinates[i][j][0] , center_coordinates[i][j][1]) , (1 / ( 3*(number_of_cells_in_line)+1) ) , 30 )
circ_collection.append(circ)
# Add all cells area
mm = Circle( Point(center_coordinates[0][0][0], center_coordinates[0][0][1]) , (1 / (3 * (number_of_cells_in_line) + 1)), 30)
for i in range(number_of_cells_in_line**2):
mm = mm + circ_collection[i]
# Subduct sum of all cells area from global domain rectangle
domain = rect - mm
mesh = generate_mesh(domain,30)
V = FunctionSpace(mesh, 'P', degree)
tol = 1e-3 # Keep in mind that it won't work if this tolerance is too harsh(like 1e-7 or -8). Set it as 1e-4 or 1e-5.
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(9999)
# Boundary numbering for all the cells in domain. For example, cell_1's boundary is ds_1, cell_2's boundary is ds_2, and so on.
for i in range(number_of_cells_in_line):
for j in range(number_of_cells_in_line):
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near( (x[0]-center_coordinates[i][j][0])**2 + (x[1]-center_coordinates[i][j][1])**2 , 1 / ( 3*(number_of_cells_in_line)+1 )**2 , tol )
bc = Boundary()
bc.mark(boundary_markers,numbering[i][j])
# Redefine boundary integration in terms of our boundary markers
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Collect Dirichlet conditions and debugging
# In our case, since there's no Dirichlet conditions, we can just ignore this part.
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
bcs.append(bc)
debug1 = False
if debug1:
# Print all vertices that belong to the boundary parts
for x in mesh.coordinates():
if bx0.inside(x, True): print('%s is on x=0' % x)
if bx1.inside(x, True): print('%s is on x=1' % x)
if by0.inside(x, True): print('%s is on y=0' % x)
if by1.inside(x, True): print('%s is on y=1' % x)
print("Number of Dirichlet conditions:", len(bcs))
if V.ufl_element().degree() == 1: # P1 elements
d2v = dof_to_vertex_map(V)
coor = mesh.coordinates()
for i, bc in enumerate(bcs):
print('Dirichlet condition %d' % i)
boundary_values = bc.get_boundary_values()
for dof in boundary_values:
print(' dof %2d: u = %g' % (dof, boundary_values[dof]))
if V.ufl_element().degree() == 1:
print(' at point %s' %
(str(tuple(coor[d2v[dof]].tolist()))))
# From the line 105 to here, we can ignore them.
u = TrialFunction(V)
v = TestFunction(V)
# 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)
# Set Linear solver parameters
prm = LinearVariationalSolver.default_parameters()
if linear_solver == 'Krylov':
prm["linear_solver"] = 'gmres'
prm["preconditioner"] = 'ilu'
prm["krylov_solver.absolute_tolerance"] = abs_tol
prm["krylov_solver.relative_tolerance"] = rel_tol
prm["krylov_solver.maximum_iterations"] = max_iter
else:
prm["linear_solver"] = 'lu'
# Computation
u = Function(V)
solve(a == L, u, bcs, solver_parameters=prm)
plot(mesh)
plt.title("Mesh Grid")
plt.show()
return u
####################################################################
### Setting Part
####################################################################
import sympy as sym
import csv
x,y = sym.symbols("x[0],x[1]")
f = 0.0
k_m = 2 # Constant k_m
R_th = 25 # # of receptor of Th cell
R_nk = 1.3 # # of receptor of NK cell
s_0 = k_m * R_th
s_1 = k_m * R_nk
q_0 = 0
q_1 = 300
q_2 = 0
q_3 = 0
data_set = [f,s_0,s_1,q_0,q_1,q_2,q_3]
data_set = [sym.printing.ccode(var) for var in data_set]
data_set = [Expression(var,degree=2) for var in data_set]
f,s_0,s_1,q_0,q_1,q_2,q_3 = data_set
# Create dictionary for paracrine secreting constants of each cell. It shapes like 'q_0':1, 'q_1':7 , ... ,'q_n':
constant_dictionary = {}
for i in range(number_of_cells_in_line**2):
constant_dictionary["q_{0}".format(i)]=data_set[i+3]
boundary_conditions = { 0: {'Robin':(constant_dictionary["q_0"],s_1)} , 1: {'Robin':(constant_dictionary["q_1"],s_0)} , 2: {'Robin':(constant_dictionary["q_2"],s_1)} , 3: {'Robin':(constant_dictionary["q_3"],s_1)} }
kappa = Constant(1) # This is corresponding to the diffusion coefficient.
u= solver_bcs(kappa,f,boundary_conditions,degree=1,linear_solver="direct")
# Save the image in Paraview with defined path
vtkfile = File("paraview_ex/2D_Two_Cells_Simulation.pvd")
vtkfile << u
plot(u)
plt.show()
)
