How to use ALE for moving boundary problems?

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.

Figure_1
Figure_2

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()

Consider the following minimal example:

from dolfin import * 
mesh = UnitSquareMesh(10, 10)

V = VectorFunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)
disp = project(as_vector((x[1],0)), V)

disp_bc = DirichletBC(V, disp, "near(x[0], 1)")
b_disp = Function(V)
disp_bc.apply(b_disp.vector())

outfile = File("mesh.pvd")
outfile << mesh
ALE.move(mesh, b_disp)
outfile << mesh
for i in range(10):
    mesh.smooth()
    outfile << mesh
1 Like

Thank you so much for giving me the example. :slightly_smiling_face:

May I ask you one more favor, if don’t mind?
I am also wondering if I can move the circle boundary by using ALE with causing minimal change of the mesh. If it is possible, could please give me a hint of it? (I want to study based on your hint to be familiar with FEniCS :slight_smile: )

The code above illustrates such a progress. Project the displacement you would like to apply on the circle as above. Mark the surface of the circle with a mesh function, similar to

Here i have an empty CG 1 vector function, Which I apply the displacement vector (x[1],0) to along the boundary where x[0] is close to 1. then I use the smoothing operation to improve the mesh quality.

2 Likes

Thank you so much. I found what I really want from one of your solutions for other users.
https://fenicsproject.discourse.group/t/moving-mesh-causes-entangled-cells/4159

I really appreciate your help. :+1:

Dear Mr.Dokken

https://fenicsproject.discourse.group/t/moving-mesh-causes-entangled-cells/4159

Can we do the same thing in FunctionSpace as well?(not vector function space like the example above) What I want to do is to compute the numerical solution u(x,y) with original mesh(star in center) at time=0, and then move the star a bit to the right side and compute the numerical solution u(x,y), like this.
How can I move the star from the center to anywhere?

Please spend some time thinking about your problem, and consider the replies that I’ve already given you. In particular, consider my previous answer: How to use ALE for moving boundary problems? - #2 by dokken

2 Likes

Thank you for your advice. As you said, I tried to analyze the code you gave me and tried to make my own code. My goal is to move circle from center point (0.5,0.5) to center point (0.57,0.57).
When I execute my code, I found that the circle shape does not preserve after moving.

Before ALE

After ALE

I think I set all boundary conditions correctly, and their numbering as well.
What I found is that if I tried to move it very tiny, from (0.5,0.5) to (0.51,0.51), for example, the shape of the circle does preserve quite well.

Could you please give me some clue or hint what the problem is?
My code is below.

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

# Generate mesh
rect = Rectangle( Point(0,0), Point(2,2) )
circ = Circle( Point(0.5,0.5),0.1,30 )
domain = rect - circ
mesh = generate_mesh(domain,30)

# Define function space
V = VectorFunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)


plot(mesh)
plt.show()

# Define boundary and set the boundaries number.
class on_inclusion(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and not(near(x[1], 0) or near(x[1], 2) or near(x[0], 0) or near(x[0], 2)))

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

class Circle_boundary(SubDomain):
    tol = 1e-3
    def inside(self, x, on_boundary):
        return on_boundary and near( (x[0]-0.5)*(x[0]-0.5) + (x[1]-0.5)*(x[1]-0.5), 0.01, self.tol )


boundaries = MeshFunction('size_t', mesh, 1, mesh.domains())
boundaries.set_all(0)
all = all_boundaries()
all.mark(boundaries,1)
inc = on_inclusion()
inc.mark(boundaries,2)
circle = Circle_boundary()
circle.mark(boundaries,3)


# Set each boundary condition.
bcs = [DirichletBC(V, Constant((0.0, 0.0)), boundaries, 1),
   DirichletBC(V, Constant((0.0, 0.0)), boundaries, 2),
   DirichletBC(V, Constant((0.07, 0.07)), boundaries, 3)]

disp_bc1 = bcs[0]
disp_bc2 = bcs[1]
disp_bc3 = bcs[2]
b_disp = Function(V)
disp_bc1.apply(b_disp.vector())
disp_bc2.apply(b_disp.vector())
disp_bc3.apply(b_disp.vector())


# Define new mesh
new_mesh = Mesh(mesh)
new_boundaries = MeshFunction("size_t", new_mesh, 1)
new_boundaries.set_values(boundaries.array())
ALE.move(mesh, b_disp)


plot(mesh)
plt.show()

You need to smooth the movement after displacement, using
mesh.smooth() as I suggested in my post. However, you are moving your circle quite alot, and mesh.smooth() will not be enough to preserve mesh quality.

You could either add in a linear elasticity equation to displace your mesh (see for instance: Drag minimization over an obstacle in Stokes-flow), or you could use the MultiMesh module in FEniCS, which is for instance described in this paper: An Error Occurred Setting Your User Cookie (Figure 1 is for instance very close to what you would like to do).

Code for this paper is available at: MultiMeshShapeOpt_code/StokesSolver.py at master · jorgensd/MultiMeshShapeOpt_code · GitHub

You could try using gradual loading as shown below:

# Set each boundary condition.
disps = [0.01, 0.01, 0.01, 0.01,0.01, 0.01, 0.01]
i = 0
for disp in disps:
    bcs = [DirichletBC(V, Constant((0.0, 0.0)), boundaries, 1),
    DirichletBC(V, Constant((0.0, 0.0)), boundaries, 2),
    DirichletBC(V, Constant((disp, disp)), boundaries, 3)]

    disp_bc1 = bcs[0]
    disp_bc2 = bcs[1]
    disp_bc3 = bcs[2]
    b_disp = Function(V)
    disp_bc1.apply(b_disp.vector())
    disp_bc2.apply(b_disp.vector())
    disp_bc3.apply(b_disp.vector())


    # Define new mesh
    new_mesh = Mesh(mesh)
    new_boundaries = MeshFunction("size_t", new_mesh, 1)
    new_boundaries.set_values(boundaries.array())
    ALE.move(mesh, b_disp)

    for j in range(50):
        mesh.smooth()
    plt.figure()
    plot(mesh,  linewidth=0.1)
    plt.savefig(f"mesh_post_{i}.png", dpi=500)
    i+=1

However, as you observe, you cannot move the mesh as much as you would like without it degenerating.

2 Likes

Thank you so much. This code helped me a lot.