Mapping between meshes in dolfinx

[chennachaos]

Hi,

I am adapting FEniCSx for FSI problems using the ALE formulation. I am having difficulties transferring data from the fluid mesh to the solid mesh and vice versa.

Say, I know fluid velocity (fluid_velo) and pressure (fluid_pres) and want to calculate the force vector on the fluid-solid boundary.

fluid_stress = mu*grad(fluid_velo) - fluid_pres*Id

fluid_traction = fluid_stress*normals_fluid

#force_on_solid = inner( solid_disp_test, fluid_traction )  * ds_solid(2)
force_on_solid = inner(solid_disp_test, fluid_traction)  * ds_fluid(5)

solid_force_temp = dolfinx.fem.assemble_vector(dolfinx.fem.form(force_on_solid))

When run, I get the following error.

Traceback (most recent call last):
  File "/home/chenna/Documents/myCode/fenics/FSI/channel-thickbeam/channel-thickbeam-FSI.py", line 622, in <module>
    solid_force_temp = dolfinx.fem.assemble_vector(dolfinx.fem.form(force_on_solid))
                                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 229, in _form
    f = ftype(
        ^^^^^^
RuntimeError: Incompatible mesh. entity_maps must be provided.

With ds_solid(2), I get a different error.

I would appreciate any help with this. Thanks in advance!

[dokken]

Without a definition of the parent mesh, and how you have extracted the fluid and solid components, it is hard to give you any further guidance. You can for instance have a look at:

Discontinuity at interface using mixed domains

I think this is the variational form that you would like: # SPDX-License-Identifier: MIT from mpi4py import MPI import dolfinx import dolfinx.fem.petsc import ufl import numpy as np from petsc4py import PETSc class NewtonSolver: max_iterations: int bcs: list[dolfinx.fem.DirichletBC] A: PETSc.Mat b: PETSc.Vec J: dolfinx.fem.Form b: dolfinx.fem.Form dx: PETSc.Vec def __init__( self, F: list[dolfinx.fem.form], J: list[list[dolfinx.fem.for…

for inspiration on how to couple data between parts of the mesh.

[chennachaos]

Thanks for your prompt response! The code is about 800 lines. So, I did not want to post it all here.

The meshes for the fluid and solid domains are separate, generated using GMSH and read and stored separately.

domain_fluid, domain_fluid_markers, domain_fluid_facet_tags = io.gmshio.read_from_msh("channel-thickbeam-2D-fluid-P2.msh", MPI.COMM_WORLD)

dx_fluid = ufl.Measure('dx', domain=domain_fluid, metadata={'quadrature_degree': 4})
ds_fluid = ufl.Measure('ds', domain=domain_fluid, subdomain_data=domain_fluid_facet_tags, metadata={'quadrature_degree': 4})



domain_solid, domain_solid_markers, domain_solid_facet_tags = io.gmshio.read_from_msh("channel-thickbeam-2D-solid-P2.msh", MPI.COMM_WORLD)

dx_solid = ufl.Measure('dx', domain=domain_solid, metadata={'quadrature_degree': 4})
ds_solid = ufl.Measure('ds', domain=domain_solid, subdomain_data=domain_solid_facet_tags, metadata={'quadrature_degree': 4})

Each mesh has its own tags for the fluid-solid interface.
The fluid mesh is shown in the figure below. The respective FE spaces and variables are defined separately for each domain.

I can solve for the fluid velocity and pressure. Next, I want to calculate the force on the interface in terms of solid DOFs so I can add the array to the RHS of the solid problem.

Let me know if you need additional information.

[dokken]

If the meshes are meshed separately you need to map the function or derived quantity from one mesh to the other, see for instance

Interpolate on boundary

Consider the minimal working example: from mpi4py import MPI import dolfinx import numpy as np L = 1 N = 5 def interface(x): return np.isclose(x[0], L) mesh_1 = dolfinx.mesh.create_box( comm=MPI.COMM_WORLD, points=[(0, 0, 0), (L, L, L)], n=[N, N, N] ) V = dolfinx.fem.functionspace(mesh_1, ("Lagrange", 2)) u = dolfinx.fem.Function(V) u.interpolate(lambda x: x[0] + np.cos(x[1]) + np.sin(np.pi * x[2])) Nx, Ny, Nz = 12, 13, 15 mesh_2 = dolfinx.mesh.create_box( comm=MPI.COMM_WORLD…

which illustrates this.

[chennachaos]

Thank you! I will go through this.

[chennachaos]

Continuing on, I realised that I needed to calculate the force on the fluid mesh first and then interpolate it onto the boundaries of the solid domain.

I am calculating the force using the following code once fluid_velo and fluid_pres are available.

fluid_stress = mu*grad(fluid_velo) - fluid_pres*Id

fluid_traction = fluid_stress*normals_fluid

force_on_solid = inner(fluid_velo_test, fluid_traction)  * ds_fluid(5)

fluid_force_temp = assemble_vector(form(force_on_solid))

I get the fluid_force_temp. I would like to know
(i) How do I export this in VTK format? I use 6-noded triangles for the velocity field and 3-noded for the pressure as a mixed element.

(ii) The correct type (like fem.Function) I should use for fluid_force_temp in order to map this force data to the solid domain using the code suggested for interpolating on the boundary.

[dokken]

@chennachaos

fluid_force_temp = assemble_vector(form(force_on_solid))

Rewrite this as

fluid_force = dolfinx.fem.Function(Velocity_space)
dolfinx.fem.assemble_vector(fluid_force.x.array, form(force_on_solid))
fluid_force.x.scatter_reverse(la.InsertMode.add)

This will assemble the force into a dolfinx function that you can then output, and in a way that you can transfer it between meshes with interpolation.

[chennachaos]

Thank you!

With your help, I calculated the force and exported it in VTK format.

While the mapping from the fluid mesh to the solid mesh also seems to work fine, there is an issue with the forces. The distribution of forces should be smoother but the values of force components are much higher at mid-nodes (see the attached figure). I checked velocity and pressure profiles. They are smoother.

I tried various options, interpolating pressure from P1 to P2, using pressure and viscous terms separately, etc., but the behaviour persists. I suspect the issue might be in evaluating the integral, but I don’t know how to check further and fix it.

This is the code block.

    normals_fluid = FacetNormal(domain_fluid)

    fluid_pres_P2 = Function(functionspace(domain_fluid,P2))
    fluid_pres_P2.interpolate(w.sub(1))

    fluid_velo_test2 = TestFunction(V2)

    force_on_solid = inner(fluid_velo_test2, -fluid_pres_P2*normals_fluid)  * ds_fluid(5)

    fluid_force_temp = Function(V2)

    dolfinx.fem.assemble_vector(fluid_force_temp.x.array, form(force_on_solid))
    fluid_force_temp.x.scatter_reverse(dolfinx.cpp.la.InsertMode.add)

[dokken]

@chennachaos

While the mapping from the fluid mesh to the solid mesh also seems to work fine, there is an issue with the forces. The distribution of forces should be smoother but the values of force components are much higher at mid-nodes (see the attached figure). I checked velocity and pressure profiles. They are smoother.

Could you clarify one thing:

1. Is the forces wrong when you compute them on the initial mesh, or when you transfer them from the fluid to solid mesh?
   Then, for my own understanding:

• Is what you want to visualize the fluid_traction, or the assembled vector force_on_solid?
  Fluid traction can be visualized (at quadrature points) with for instance:
  Visualizing quadrature functions as point clouds — scifem
  which even works for boundary integrals. See for instance: presentations/fenics_ecole_des_ponts/code/expressions.py at dokken/edp · jorgensd/presentations · GitHub
  with corresponding plot:
  
  

[chennachaos]

The force vector I got earlier was correct but my interpretation was wrong and there was a bug in the boundary conditions. I fixed everything. The FSI code is working now.

Thank you very much for your assistance!

[dokken]

It would be good if you could post your solutions, to the benefit of the open source community:)

[chennachaos]

Certainly! I will share via GitHub once I test it further with more examples.

Scripts for CFD problems using a coupled velocity-pressure formulation with Taylor-Hood elements are shared at GitHub - chennachaos/CFDwithFEniCSx: Computational Fluid Dynamics with FEniCSx library.