Applying mesh.Mesh() on a mesh domain makes it incompatible

mesh
1

Hello, I would appreciate it if someone could explain why applying the mesh.Mesh() function to the domain makes it incompatible with the fem.form() function?

Here is my MWE:

from dolfinx import mesh
from dolfinx import fem
import ufl
from mpi4py import MPI
import numpy as np
comm = MPI.COMM_WORLD

no_elem = 5
domain = mesh.create_unit_square(comm, no_elem, no_elem, mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("CG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

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

mesh.Mesh(domain.comm, domain.topology, domain.geometry, domain.ufl_domain())

V = fem.FunctionSpace(domain, ("CG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

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

which produces the following error:

RuntimeError                              Traceback (most recent call last)
/home/hosseini/Programming/Fenicsx/010411/packeshkon.ipynb Cell 1 in <cell line: 24>()
     21 u = ufl.TrialFunction(V)
     22 v = ufl.TestFunction(V)
---> 24 a = fem.form(
     25     ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
     26 )

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:168, in form(form, dtype, form_compiler_params, jit_params)
    165         return list(map(lambda sub_form: _create_form(sub_form), form))
    166     return form
--> 168 return _create_form(form)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:163, in form.<locals>._create_form(form)
    160 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    161 return form argument"""
    162 if isinstance(form, ufl.Form):
--> 163     return _form(form)
    164 elif isinstance(form, collections.abc.Iterable):
    165     return list(map(lambda sub_form: _create_form(sub_form), form))

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:157, in form.<locals>._form(form)
    151 # Subdomain markers (possibly None for some dimensions)
    152 subdomains = {_cpp.fem.IntegralType.cell: subdomains.get("cell"),
...
     52 ffi = cffi.FFI()
---> 53 super().__init__(ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
     54                  V, coeffs, constants, subdomains, mesh)

RuntimeError: Incompatible mesh
2

Keep in mind that dolfinx.mesh.Mesh is a class and not a function. It returns an instance of itself.

As for the error: domain.ufl_domain() actually has a reference to the mesh it was created with in domain.ufl_domain().ufl_cargo(). Consider:

print(domain is domain.ufl_domain().ufl_cargo())

which yields True. As such you should create a ufl domain explicitly (or copy the other without the cargo), e.g.,

domain_ufl = ufl.Mesh(ufl.VectorElement("Lagrange", domain.ufl_cell(), 1))
mesh.Mesh(domain.comm, domain.topology, domain.geometry, domain_ufl)

“Under the hood” this ufl_cargo() is used implicitly by ufl.dx to be able to interpret the problem integration domain.

Edit: Thanks to @conpierce8 for the following suggestion which will copy the previous mesh’s UFL domain without having to define the vector element explicitly:

domain_ufl = ufl.Mesh(domain.ufl_domain().ufl_coordinate_element())
mesh.Mesh(domain.comm, domain.topology, domain.geometry, domain_ufl)
2 Likes