# Solve equation in only one phase of a two-phase problem

**URL:** <https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610>\
**Category:** General\
**Created:** [December 23, 2024, 7:18am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610 "2024-12-23T07:18:02Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [December 23, 2024, 7:18am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/1 "2024-12-23T07:18:02Z")

</div>

I have a two-phase system for FSI simulation. With everything built on a single big mesh, I need to compute variables on the solid phase **only** as shown in the attached figure.

 ![two_phase_FSI](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/09ceb98a26fc7c7f8f6e8edc7a5a12a69cba94ff.jpeg)  
(a) mask and velocity vector arrows, (b) velocity magnitude. (c) result of \omega=0.5\nabla \times u.

The test in the figure shows that when the velocity field `u` is given, the rotation velocity `curlU` of the solid is computed via projection of \omega=0.5\nabla \times u. Interpolation cannot do this as the velocity is continuous but not differentiable on the interface.

To do so, I have 2 scenarios:

1. Do projection on the big mesh with the help of a function `mask` for solid and apply Dirichlet BC for fluid-only DOFs as

```auto
a = fem.form(ufl.inner(Pv, w) * mask * ufl.dx)
L = fem.form(ufl.inner( v, w) * mask * ufl.dx)

problem = LinearProblem(a, L, bcs, u=target_func,
        petsc_options={"ksp_type": "preonly", 
                       "pc_type": "lu", 
                       "pc_factor_mat_solver_type": "mumps"})
problem.solve()

```

The full MWE is attached in the end.

1. Create a submesh for solid with `submesh = mesh.create_submesh(big_mesh, tdim, solid_cells)[0]`. First, I interpolate `u` from the big mesh to the submesh. Then, solve the projection on the submesh to get `curlU`. Finally, interpolate `curlU` to the function on the big mesh. The whole process involves 2 non-matching meshes interpolation.

In my application, I do two-phase FSI simulations in 3D with 2nd order iso-parametric meshes. Scenario 1 gets really slow as it involves too many unnecessary DOFs. Scenario 2 is fast but it often crashes with interpolation related GJK or Newton method failed to converge runtime errors. I sense the scenario 1 could be tuned for this problem, but still I am looking for a hand on this. Thanks.

MWE with DOLFINx-v0.8.0:

```auto
from dolfinx import mesh, io, fem
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import numpy as np
import ufl

comm = MPI.COMM_WORLD
eps = 1e-8

def project(v, target_func, mask, bcs=[]) -> None:

    V = target_func.function_space
    w = ufl.TestFunction(V)
    Pv = ufl.TrialFunction(V)
    a = fem.form(ufl.inner(Pv, w) * mask * ufl.dx)
    L = fem.form(ufl.inner( v, w) * mask * ufl.dx)

    problem = LinearProblem(a, L, bcs, u=target_func,
            petsc_options={"ksp_type": "preonly", 
                           "pc_type": "lu", 
                           "pc_factor_mat_solver_type": "mumps"})
    problem.solve()
    target_func.x.scatter_forward()
    
def solid_domain(x):
    return np.logical_and(
        np.logical_and(x[0] > 0.25 - eps, x[0] < 0.75 + eps),
        np.logical_and(x[1] > 0.25 - eps, x[1] < 0.75 + eps))
    
_mesh = mesh.create_unit_square(comm, 32, 32)
tdim = _mesh.topology.dim

# Mask on a DG-0 function space
DG0 = fem.functionspace(_mesh, ("DG", 0))
mask = fem.Function(DG0)
_mesh.topology.create_connectivity(tdim, tdim)

solid_cells = mesh.locate_entities(_mesh, tdim, solid_domain)
solid_dofs = fem.locate_dofs_topological(DG0, tdim, solid_cells)
mask.x.array[solid_dofs] = 1.0

with io.XDMFFile(comm, "output/mask.xdmf", "w") as file:
    file.write_mesh(_mesh)
    file.write_function(mask)

# Build velocity field for 2 phases
omega = 1.0 # rotation velocity
center = ufl.as_vector([0.5, 0.5])
x = ufl.SpatialCoordinate(_mesh) - center
r = ufl.sqrt((x[0] - 0.5)**2 + (x[1] - 0.5)**2)
tri = ufl.as_vector([- x[1] / (r + eps), x[0] / (r + eps)])
V = fem.functionspace(_mesh, ("Lagrange", 1, (tdim,)))
Q = fem.functionspace(_mesh, ("Lagrange", 1))
u_fluid = fem.Function(V)
u_solid = fem.Function(V)
u = fem.Function(V)
mask_V = fem.Function(V)

solid_dofs_V = fem.locate_dofs_topological(V, tdim, solid_cells)
for d in range(tdim):
    mask_V.x.array[solid_dofs_V * tdim + d] = 1.0
mask_V.x.scatter_forward()

expr = fem.Expression(omega * r * tri, V.element.interpolation_points())
u_solid.interpolate(expr)
u_solid.x.scatter_forward()

u_fluid.interpolate(lambda x: ((x[0] - 0.5), (x[1] - 0.5)))
u_fluid.x.scatter_forward()

u.x.array[:] = (u_solid.x.array[:] * mask_V.x.array[:] +
                u_fluid.x.array[:] * (1 - mask_V.x.array[:]))

with io.XDMFFile(comm, "output/u.xdmf", "w") as file:
    file.write_mesh(_mesh)
    file.write_function(u)
    
curlU = fem.Function(Q)
total_dofs_Q = (Q.dofmap.index_map.size_local + 
                Q.dofmap.index_map.num_ghosts)
solid_dofs_Q = fem.locate_dofs_topological(Q, tdim, solid_cells)
fluid_only_dofs_Q = np.setdiff1d(np.arange(total_dofs_Q), solid_dofs_Q)

expr = 1. / 2. * ufl.curl(u)
f = fem.Function(Q)
bc_f = fem.dirichletbc(f, fluid_only_dofs_Q)
project(expr, curlU, mask, bcs=[bc_f])

with io.XDMFFile(comm, "output/omega.xdmf", "w") as file:
    file.write_mesh(_mesh)
    file.write_function(curlU)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 23, 2024, 7:29am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/2 "2024-12-23T07:29:54Z")

</div>

> [@chenyongxin](#):
>
> Create a submesh for solid with `submesh = mesh.create_submesh(big_mesh, tdim, solid_cells)[0]`. First, I interpolate `u` from the big mesh to the submesh. Then, solve the projection on the submesh to get `curlU`. Finally, interpolate `curlU` to the function on the big mesh. The whole process involves 2 non-matching meshes interpolation.

You don’t need to interpolate u from the parent to the sub mesh. At least if you upgrade to 0.9, you can have u in a variational form on the submesh.

See for instance: [Multiphysics: Solving PDEs on subdomains — FEniCS Workshop](http://jsdokken.com/FEniCS-workshop/src/multiphysics/submeshes.html)

---

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [December 24, 2024, 10:09am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/3 "2024-12-24T10:09:46Z")

</div>

I had a go with the matched section `Integration with function from parent and sub-mesh` in the tutorial and the result looks correct with my case. Please correct me if I am wrong with the following idea:

`u` is defined on the parent mesh, while `omega` function, trail function and test function are defined on the submesh. **It needs extra `entity_maps` info to explicitly formulate and assemble to A and b**. Then solve projection on the submesh to get `omega`.

Several follow-up questions:  
(1) I noticed `entity_map` argument of `fem.form` also exists in v0.8.0 with the same description, but I failed to call `a_solid, L_solid = dolfinx.fem.form(ufl.system(F), entity_maps=entity_maps)` as that in v0.9.0. Is it a v0.9.0-only related feature?  
(2) With the use of `entity_maps` of `fem.form`, the `Ax=b` problem can only be solved explicitly rather than to use `LinearProblem` wrapper?  
(3) My practical application is FSI + ALE, so the mesh moves in every step. If with a submesh, any idea to associate the movement of a submesh to the parent mesh?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 24, 2024, 10:15am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/4 "2024-12-24T10:15:28Z")

</div>

> [@chenyongxin](#):
>
> `u` is defined on the parent mesh, while `omega` function, trail function and test function are defined on the submesh. **It needs extra `entity_maps` info to explicitly formulate and assemble to A and b**. Then solve projection on the submesh to get `omega`.

You don’t need to project anything. You can leave u in the variational form, and it will behave as any function defined in a standard variational form.

> [@chenyongxin](#):
>
> Several follow-up questions:  
> (1) I noticed `entity_map` argument of `fem.form` also exists in v0.8.0 with the same description, but I failed to call `a_solid, L_solid = dolfinx.fem.form(ufl.system(F), entity_maps=entity_maps)` as that in v0.9.0. Is it a v0.9.0-only related feature?

Various fixes related to submeshes, as well as a slight api change happened from 0.8 to 0.9

> [@chenyongxin](#):
>
> (2) With the use of `entity_maps` of `fem.form`, the `Ax=b` problem can only be solved explicitly rather than to use `LinearProblem` wrapper?

You can pass in compiled forms to LinearProblem.

> [@chenyongxin](#):
>
> 1. My practical application is FSI + ALE, so the mesh moves in every step. If with a submesh, any idea to associate the movement of a submesh to the parent mesh?

Then you should interpolate any displacement in the parent mesh onto the submesh and move the mesh afterwards

---

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [January 6, 2025, 5:56am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/5 "2025-01-06T05:56:54Z")

</div>

Assume I have two functions `f` and `fs` defined on parent mesh and submesh respectively. Is it able to find a DOF map from submesh to parent mesh? Say, the `dof_sub_to_parent_map` array in the following code:

```auto
tdim = cell_tags.dim
submesh, cell_map, vertex_map, _ = dolfinx.mesh.create_submesh(mesh, tdim, cell_tags.find(marker))

# Function space and functions
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (tdim,)))
Vs = dolfinx.fem.functionspace(submesh, ("Lagrange", 1, (tdim,)))
f = dolfinx.fem.Function(V)
fs = dolfinx.fem.Function(Vs)

# Wanted: dof_sub_to_parent_map
dof_sub_to_parent_map = ...
f.x.array[dof_sub_to_parent_map] = fs.x.array[:]

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 9, 2025, 8:23am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/6 "2025-01-09T08:23:27Z")

</div>

You can make this function if you want to.  
However, it is much quicker and easier call

```python
submesh, parent_cells, _, _ = create_submesh(mesh, tdim, cells)
f.interpolate(fs, cells0=np.arange(len(parent_cells)), cells1=parent_cells)

```

as done in:

> <https://github.com/FEniCS/dolfinx/blob/main/python/test/unit/fem/test_interpolation.py#L1057-L1080>

---

<div class="post-metadata">

**Author:** ![chenyongxin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/chenyongxin/32/6501_2.png) [@chenyongxin](https://fenicsproject.discourse.group/u/chenyongxin)\
**Post date:** [January 10, 2025, 3:49am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/7 "2025-01-10T03:49:25Z")

</div>

It works with interpolation. My understanding is the interpolation between submesh and parent mesh is matching meshes interpolation. An unmatured question: for the quoted interpolation, does DOLFINx use 1-to-1 assigning on specified cells or interpolate through compute?

> [@dokken](#):
>
> `f.interpolate(fs, cells0=np.arange(len(parent_cells)), cells1=parent_cells)`

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 10, 2025, 5:59am UTC](https://fenicsproject.discourse.group/t/solve-equation-in-only-one-phase-of-a-two-phase-problem/16610/8 "2025-01-10T05:59:25Z")

</div>

> [@chenyongxin](#):
>
> It works with interpolation. My understanding is the interpolation between submesh and parent mesh is matching meshes interpolation

Yes it is matching interpolation.

> [@chenyongxin](#):
>
> An unmatured question: for the quoted interpolation, does DOLFINx use 1-to-1 assigning on specified cells or interpolate through compute?

As the spaces for both parent and submesh are Lagrange 1, it becomes a task of assigning degrees of freedom (not tabulate/compute basis functions on the input space)