Mixed formulation with different trial and test mesh

dolfinx
1

Hello everyone,

I am interested in a mixed formulation where the test and trial spaces are defined on different (but nested) meshes. Is this somehow possible in FEniCSx?

Thanks!

2

I am adding a MWE that tries to implement it.

from dolfinx import mesh, fem
import ufl
from mpi4py import MPI

msh_trial = mesh.create_unit_square(MPI.COMM_WORLD, 5, 5, mesh.CellType.quadrilateral)
msh_test = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.quadrilateral)

U = fem.functionspace(msh_trial, ("Lagrange", 2))
V = fem.functionspace(msh_test, ("Lagrange", 1))

u = ufl.TrialFunction(U)
v = ufl.TestFunction(V)
dx = ufl.Measure("dx", domain=msh_trial)

a = ufl.inner(u, v) * dx
A = fem.assemble_matrix(fem.form(a))

As it could be expected, the problem is in the meshes, as I get the following error.

RuntimeError: Incompatible mesh. entity_maps must be provided.

Note that I still get the same error with msh_trial = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.quadrilateral).

3

Hello Ivan,

currently dolfinx does not allow for this use case. Internally it is assumed that a variational form lives on a single mesh and thus integration only supports assembly with different function spaces defined on the same mesh.

However, there is some recent work regarding the interaction of meshes that are produced by refinement operations - in the context of geometric multigrid. This should (in the future - see this PR) allow for the transfer of functions between the two spaces you proposed. Note however that there is currently no quadrilateral refinement routine available, so this requires the usage of triangles/tetrahedrons.

In general the code of how such a matrix between two different spaces (with different meshes) may be constructed can be found here. Note: currently this only supports P^1 elements and needs to be extended to general dof maps in the future.

Another thought: Maybe introducing an auxiliary space P^2 on the finer mesh as an ‘integration space’ (for at least a qualitative analysis) could be worth considering.

Hope that helps,
Paul

1 Like
4

Good morning folks,

Just in sake of curiosity, what are the physical problem that require a FE method where the trial and test function are not defined on the same space/mesh ?

Regards,
Pierre.

5

I know some mixed dimensional methods with non-conforming grids where this could be required, ref: https://arxiv.org/pdf/2312.16565

6

Another possible application I’m aware of is with partition of unity methods such as XFEM or GFEM. Here the approximation space is constructed by connecting local approximation spaces (FEM spaces themselves) \{V_i \}_{i\in I} by a partition of unity \{ \phi_i \}_{i \in I} (most commonly this is just given by P^1 elements) to form an approximation space \{ \phi_i V_i \}_{i \in I}. In this setting it is often advantageous to be able to control if the PU is considered in the testing and trial arguments. This results in more control over runtime (not using the PU is faster to integrate/assemble) and allows for great preconditioning opportunities which make use of the localized approximation space structure.

Equal ideas can be employed by changing the quadrature rules used for the assembly of forms with classic FEM, which may be facilitated by using a different space/mesh.

7

Thank you for your reply and your efforts. Unfortunately, I would need to use higher degree than p=1 (something like p=5,6). At this stage of my project, I am only doing some preliminary experiments, and an ugly/inefficient (but working) solution would be more than enough.

Say that my coarse mesh has quadrilateral elements with mesh-size h and I use P^k trial functions. Then I would like to use P^1 test functions on a finer mesh with mesh-size h/k (so that the dimension is the same). I thought that one solution could be to use P^k “auxiliary” trial functions on the finer mesh, and then represent the P^k trial functions on the coarse mesh as linear combinations of the auxiliary ones. I am perfectly aware that this is not efficient at all, but is fine for the moment.

To do so, I would need to assemble the matrix that expresses the P^k basis functions on the coarse mesh as a linear combination of the ones on the fine mesh. It should be enough to evaluate the basis function of the coarse mesh on the nodes of the fine one. Is there any way to do so in FenicsX, possibly assembling a sparse and not a dense matrix? Maybe also @dokken can help with this.

8

I might have found a way to do it, but it is extremely slow. Any advice for making it more efficient?

from dolfinx import fem, mesh
from mpi4py import MPI
from tqdm import tqdm 
import numpy as np

p = 5
nelems = 5

msh_to = mesh.create_unit_square(
    comm=MPI.COMM_WORLD,
    nx=nelems*p,
    ny=nelems*p,
    cell_type=mesh.CellType.quadrilateral,
)
V_to = fem.functionspace(msh_to, ("Lagrange", p))  # Trial space
ndofs_to = V_to.dofmap.index_map.size_global

msh_from = mesh.create_unit_square(
    comm=MPI.COMM_WORLD,
    nx=nelems,
    ny=nelems,
    cell_type=mesh.CellType.quadrilateral,
)
V_from = fem.functionspace(msh_from, ("Lagrange", p))  # Trial space
ndofs_from = V_from.dofmap.index_map.size_global

interp_data = fem.create_nonmatching_meshes_interpolation_data(V_to.mesh, V_to.element, V_from.mesh)

interpolation_matrix = np.zeros((ndofs_to, ndofs_from))
for i in tqdm(range(ndofs_from)):
    u_from = fem.Function(V_from)
    u_from_array = np.zeros_like(u_from.x.array)
    u_from_array[i] = 1
    u_from.x.array[:] = u_from_array
    u_to = fem.Function(V_to)
    u_to.interpolate(u_from, nmm_interpolation_data=interp_data)
    interpolation_matrix[:, i] = u_to.x.array