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.

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