ALE formulation with prescribed boundary motion

Question

Dear all,
I just started working with FEniCS for my master thesis. The project is about the flow simulation inside the left ventricle of the heart. Therefore, I have a moving geometry which I collected from Magnetic Resonance Imaging. Now, I am trying to use FEniCS to simulate the flow in the moving domain.
At each time step I would like to load the new mesh (which is already interpolated from the first one, all the meshes are topologically consistent, same connectivity elements), find the velocity from the mesh displacement, use this velocity inside the convective term of the weak formulation, and use it as DirichletBC on the left ventricular wall. I tried something like that, but it does not work and I do not know what I have to change:

mesh_1 = Mesh("LV_Fenics/LV_00_vol.xml")

V = VectorFunctionSpace(mesh_1,"CG",3)
# Extract dofs
dms = [V.sub(i).dofmap().dofs() for i in range(3)]

u_mesh = Function(V)

while t < T + DOLFIN_EPS:
   # Load next geometry
   mesh_2 = Mesh("LV_Fenics/LV_0"+str(i)+"_vol.xml")
   # Compute mesh velocity
   u_mesh_new = Function(V)
   # Get the coordinates in dof order
   dof_coor_1 = V.dofmap().tabulate_all_coordinates(mesh_1)
   # Coordinates of dofs 
   coor_x_1 = dof_coor_1[dms[0]]
   coor_y_1 = dof_coor_1[dms[1]]
   coor_z_1 = dof_coor_1[dms[2]]

   dof_coor_2 = V.dofmap().tabulate_all_coordinates(mesh_2)
   coor_x_2 = dof_coor_2[dms[0]]
   coor_y_2 = dof_coor_2[dms[1]]
   coor_z_2 = dof_coor_2[dms[2]]

   u_mesh_new.vector()[dms[0]] = (coor_x_2-coor_x_1)/dt
   u_mesh_new.vector()[dms[1]] = (coor_y_2-coor_y_1)/dt
   u_mesh_new.vector()[dms[2]] = (coor_z_2-coor_z_1)/dt

   u_mesh.assign(u_mesh_new)

I perform this point at each iteration in time. The mesh velocity is present inside the weak formulation in the convective term:

    # Tentative velocity
U = 0.5*(u0 + u)
F1 = rho*(1/k)*inner(u-u0, v)*dx \
   + rho*inner(grad(u0)*(u0), v)*dx \
   + inner(sigma(U, p0), epsilon(v))*dx \
   + inner(p0*n, v)*ds \
   - mu*inner(grad(U).T*n, v)*ds \
   - inner(f, v)*dx

a1 = lhs(F1)
L1 = rhs(F1)

As last point, I have to use the mesh velocity u_mesh as dirichlet boundary condition on the boundary mesh. Do you have any suggetions on this?
In advance many thanks for your help!

Comments

It is a bit hard to understand what you are trying to do with the DirichletBC and what is not working. As far as I can see it should be possible to directly use a Function as the DirichletBC value.

If you remove all the ALE parts of the problem and reduce it to a minimum it will be easier to help. You can then post a minimal complete code example. Using a known Function as a DirichletBC is not really related to ALE and can be used also in other applications, so you can remove all the mesh handling if it is not related to the question you are asking.

Regarding ALE, one thing to have in mind is that unless you redefine your equation F1 every time it will still be referring to the old mesh (through the function space of u). It may be better to keep the original mesh and deform it each time step according to the mesh velocities that you calculate. This may be what you are already doing.

---

Thank you for your reply. I wrote also the ALE part because the velocity u_mesh computed in this part is the one which I will use as Dirichlet boundary condition on the wall.
Therefore I define the BC as follow:

V = VectorFunctionSpace(mesh_1,"CG",1)
bcu_wall = DirichletBC(V, u_mesh, subdomain_wall, 1)

where subdomain_wall is a MeshFunction to mark all the boundary facets with 1.
Now, I would like to update the BC in each time step after moving the mesh and compute the mesh velocity u_mesh:

while t < T + DOLFIN_EPS:
   bcu_wall.apply(u_mesh.vector())
   ....
   u_mesh.assign(u_mesh_new)

where u_mesh_new is computed as I wrote before in the ALE description. This formulation stops working when I try to apply the new vector of the mesh velocity to the BCs.

Another point about the movement of the mesh. I tried to use the class move() provided by dolfin in the following way:

mesh_1.move(mesh_2)

And I had the following error:

*** Error: Unable to access mesh data.
*** Reason: Mesh data array named "parent_vertex_indices" does not exist.
*** Where: This error was encountered inside MeshData.cpp.

I understand that my mesh.data() object is empty but I have no idea why it has no data.
Many thanks for any help you can give me on this topic.

Answer 1

When you run bcu_wall.apply() you are applying the boundary conditions, either to the left hand side (if you pass a matrix) or to the right hand side (if you pass a vector). You seem to run apply in order to update the boundary condition value, previously passed to the constructor of the DirichletBC. This is not needed and is indeed not correct as this will modify u_mesh in this case, not the BC itself.

The BC will ask the object you passed, in this case u_mesh, for the boundary values. So, as long as u_mesh is updated by your code then the BC will be updated as well and you do not need to tell the BC about updates to u_mesh.

For the mesh moving problem, I have never tried to move the mesh according to some other mesh. Try something like this (not tested, may contain minor errors)

V = VectorFunctionSpace(mesh_1,"CG",1)
displacement = Function(V) # same function space as u_mesh

displacement.vector().zero()
displacement.vector().axpy(dt, u_mesh)

mesh.move(displacement)

Comments

Thank you very much for the suggestion on the boundary conditions, if I update u_mesh without calling bc.apply() it works fine.
However, using mesh.move() I always have the error about the mesh data "parent_vertex_indices" which does not exist for my mesh.
I tried to create this array like this:

mesh_1.data().create_array("parent_vertex_indices", len(mesh_1.coordinates()))

and also for the second mesh:

mesh_2.data().create_array("parent_vertex_indices", len(mesh_1.coordinates()))

But also with this I have the same error:

*** Error: Unable to access mesh data.
*** Reason: Mesh data array named "parent_vertex_indices" does not exist.

using mesh_1.move(displacement).
Do you have an idea about how I could create an array of indices called "parent_vertex_indices" inside the object meshdata?
Many thanks.