Hi, all
I am just a beginner of the fenics user. I tried to make 64 spheres in cubic and each of sphere has their own boundary conditions, like du/dn = q - s*u (Robin boundary).
When I made my own code and run it, I found that it always failed to recognize all 64 spheres. Instead, they only recognize a few of spheres like the image below.
Could you please help me what the problem is? How can I define circle or sphere boundary?
For example, if I want to set sphere boundary like Sphere( Point(1,1,1),0.1 ), are there any way to define
Below is the code that I made.
from mshr import *
from fenics import *
from dolfin import *
import matplotlib.pyplot as plt
# Set Cubic length, cubic margin, cell_radius, and cell to cell distance(the closest distance, not center to center.)
cube_l = 90 # Cubic length
cube_m = 10 # Cubic margin
c_r = 5 # Cell radius
c_d = 10 # Cell to cell distance
n = 0
while n * (2 * (c_r) + (c_d)) + 2 * (c_r) <= (cube_l) - 2 * (cube_m):
if n * (2 * (c_r) + (c_d)) + 2 * (c_r) == (cube_l) - 2 * (cube_m):
break
elif (n + 1) * (2 * (c_r) + (c_d)) + 2 * (c_r) > (cube_l) - 2 * (cube_m):
break
n = n + 1
else:
print("Cell cannot be larger than Cubic.")
# Redefine the cubic length so that it becomes compact to cell size and margin.
cube_l = n * (2 * c_r + c_d) + 2 * c_r + 2 * cube_m
total_cell = (n + 1) ** 3
# Sphere numbering for boundary marker
# List structure is [ [ [0,1,2,3...] ex> numbering[i][j][0] = 0 , numbering[i][j][1] = 1
numbering = []
tt = 0
for i in range(n + 1):
numbering.append([])
for j in range(n + 1):
numbering[i].append([])
for k in range(n + 1):
numbering[i][j].append(tt)
tt = tt + 1
# Coordinate setting
# Each cell's center coordinates information. [ [ [ [x0,y0,y0] , [x1,y0,z0] , ... <- List structure. the coordinate information is in center_coordinates[i][j][k] = [x,y,z]
center_coordinates = []
for i in range(n + 1):
center_coordinates.append([])
for j in range(n + 1):
center_coordinates[i].append([])
for k in range(n + 1):
center_coordinates[i][j].append([cube_m + c_r + k * (2 * c_r + c_d), cube_m + c_r + j * (2 * c_r + c_d),
cube_m + c_r + i * (2 * c_r + c_d)])
####################################################
### 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
cube = Box(Point(0.0, 0.0, 0.0), Point(cube_l, cube_l, cube_l))
# Collect sphere volume
sphere_collection = []
# Draw each cells in domain
for i in range(n + 1):
for j in range(n + 1):
for k in range(n + 1):
sphere = Sphere(Point(center_coordinates[i][j][k][0], center_coordinates[i][j][k][1], center_coordinates[i][j][k][2]), c_r)
sphere_collection.append(sphere)
# Add all cells area
# mm is sum of all spheres(cells) inside the cubic.
mm = Sphere(Point(center_coordinates[0][0][0][0], center_coordinates[0][0][0][1], center_coordinates[0][0][0][2]), c_r)
for i in range((n + 1) ** 3):
mm = mm + sphere_collection[i]
# Subtract sum of all cells area from global domain cubic.
domain = cube - mm
mesh = generate_mesh(domain, 30) # The more you have cells in cubic, the larger number you need to set.
V = FunctionSpace(mesh, 'CG', degree)
tol = 1 # Keep in mind that it won't work if this tolerance is too harsh(like 1e-7 or -8). Set it as 1e-1 or 1e-2 in 3D case.
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(n + 1):
for j in range(n + 1):
for k in range(n + 1):
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near((x[0] - center_coordinates[i][j][k][0]) ** 2 + (x[1] - center_coordinates[i][j][k][1]) ** 2 + (x[2] - center_coordinates[i][j][k][2]) ** 2,c_r ** 2, tol)
bc = Boundary()
bc.mark(boundary_markers, numbering[i][j][k])
# 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)
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs)
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()
return u
####################################################################
### Setting Part
####################################################################
import sympy as sym
import csv
x, y, z = sym.symbols("x[0],x[1],x[2]")
f = 0.0
k_m = 2 # Constant k_m
R_th = 25 # # of receptor of Th1 cell
R_nk = 1.3 # # of receptor of NK cell
s_0 = k_m * R_th
s_1 = k_m * R_nk
# Import paracrine secreting rate from .txt file. It indicates each cell's secreting rate.
secreting_list = []
with open('PDE_t', 'r') as ff:
for line in ff:
secreting_list.append(line.strip())
# Each entry's type in secreting_list is string. So we need to change them to int or float.
aa = []
for i in range((n + 1) ** 3):
aa.append(float(secreting_list[i]))
# Now, the list aa's all entries are float
# Distinguish Th secreting cells from NK cell.
Th_cell_nr = []
for i in range(len(aa)):
if aa[i] > 3000:
Th_cell_nr.append(i)
aa1 = [f, s_0, s_1]
# Merge two list aa and aa1 as a name of data_set. Now we have all boundary condition's constants in here.
data_set = aa1 + aa
data_set = [sym.printing.ccode(var) for var in data_set]
data_set = [Expression(var, degree=2) for var in data_set]
# Redefine values of f, s_0 (Th1), and s_1 (NK).
f = data_set[0]
s_0 = data_set[1]
s_1 = data_set[2]
# Create dictionary for paracrine secreting constants of each cell. It shapes like 'q_0':1, 'q_1':7 , ... ,'q_((n+1)**3-1)'
constant_dictionary = {}
for i in range((n + 1) ** 3):
constant_dictionary["q_{0}".format(i)] = data_set[i + 3]
# Define each cell's boundary conditions.
boundary_conditions = {}
for i in range((n + 1) ** 3):
if i not in Th_cell_nr:
boundary_conditions[i] = {'Robin': (constant_dictionary["q_{0}".format(i)], s_1)}
if i in Th_cell_nr:
boundary_conditions[i] = {'Robin': (constant_dictionary["q_{0}".format(i)], s_0)}
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/ThreeD_Two_Cells_Simulation_Th1_NK.pvd")
vtkfile << u
I think this part has some problems…
mm = Sphere(Point(center_coordinates[0][0][0][0], center_coordinates[0][0][0][1], center_coordinates[0][0][0][2]), c_r)
for i in range((n + 1) ** 3):
mm = mm + sphere_collection[i]
# Subtract sum of all cells area from global domain cubic.
domain = cube - mm
mesh = generate_mesh(domain, 30) # The more you have cells in cubic, the larger number you need to set.
V = FunctionSpace(mesh, 'CG', degree)
tol = 1 # Keep in mind that it won't work if this tolerance is too harsh(like 1e-7 or -8). Set it as 1e-1 or 1e-2 in 3D case.
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(n + 1):
for j in range(n + 1):
for k in range(n + 1):
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near((x[0] - center_coordinates[i][j][k][0]) ** 2 + (x[1] - center_coordinates[i][j][k][1]) ** 2 + (x[2] - center_coordinates[i][j][k][2]) ** 2,c_r ** 2, tol)
bc = Boundary()
bc.mark(boundary_markers, numbering[i][j][k])
