Spatial derivative on coupled problem with different dimension

UFL
1

Hello everyone,

I would like to implement a problem that could be written like this (\Omega is a volumic domain and \Gamma is a surface of \Omega):
(p, q)\in\Omega\times\Gamma
(v, u)\in\Omega\times\Gamma
\begin{cases}\int_\Omega p \cdot vd\Omega + \int_\Gamma q\cdot v d\Gamma = 0 \\ \int_\Gamma p \cdot \frac{du}{dx}d\Omega + \int_\Gamma q\cdot u d\Gamma = 0\end{cases}
The mathematical meaning of this formulation is not relevant, but this formulation illustrates my issue, and more important, the code that comes below.
The following code :

from mpi4py import MPI
import gmsh

from dolfinx.io import gmshio

from ufl import (TestFunction, TrialFunction, TestFunctions, TrialFunctions,
                 Measure,
                 FacetNormal, CellNormal,
                 inner, grad, div,
                 Dn,
                 MixedFunctionSpace,
                 extract_blocks)
from basix.ufl import element
from dolfinx.fem import functionspace, form, petsc

import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio

import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio
import dolfinx.mesh as dmesh
import numpy as np

def generate_cube_gmshio(side_length=1.0, filename="cube.msh"):
    """
    Generate a cube with side length `side_length`,
    export the mesh to `filename`, and return:
      - final_mesh (full mesh)
      - cell_tags (tags for volume)
      - facet_tags (tags for surfaces)
      - submesh (mesh of one face)
      - entity_map (relation between submesh entities and parent mesh)
    """
    # Initialize the Gmsh API
    gmsh.initialize()

    comm = MPI.COMM_WORLD
    model_rank = 0
    model = gmsh.model
    model.add("cube")

    # Create the cube
    box_tag = model.occ.addBox(0, 0, 0, side_length, side_length, side_length)
    model.occ.synchronize()

    # Retrieve surfaces (dim=2)
    surfaces = model.occ.getEntities(dim=2)

    # Identify the face at z=0
    face_z0 = None
    for s in surfaces:
        com = model.occ.getCenterOfMass(s[0], s[1])
        if abs(com[2]) < 1e-8:  # z=0 face
            face_z0 = s[1]
            break
    if face_z0 is None:
        raise RuntimeError("Face at z=0 not found!")

    # Add physical groups
    model.addPhysicalGroup(3, [box_tag], tag=1)         # volume
    model.addPhysicalGroup(2, [face_z0], tag=10)        # surface z=0

    # Mesh options
    gmsh.option.setNumber("Mesh.ElementOrder", 2)
    gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
    model.mesh.generate(3)

    # Save mesh
    gmsh.write(filename)

    # Convert to FEniCSx
    final_mesh, cell_tags, facet_tags = gmshio.model_to_mesh(
        model, comm, model_rank
    )

    gmsh.finalize()

    # --- Create submesh of face z=0 ---
    # Extract entities of the facet with tag=10
    facets = facet_tags.find(10)

    # Create the submesh (dim=2 since we want a surface)
    submesh, entity_map = dmesh.create_submesh(
        final_mesh, dim=2, entities=facets
    )[0:2]

    tdim = final_mesh.topology.dim
    fdim = tdim - 1
    submesh_tdim = submesh.topology.dim
    submesh_fdim = submesh_tdim - 1

    # Create the entity maps and an integration measure
    mesh_facet_imap = final_mesh.topology.index_map(fdim)
    mesh_num_facets = mesh_facet_imap.size_local + mesh_facet_imap.num_ghosts
    extract = np.array([entity_map.tolist().index(entity) if entity in entity_map else -1 for entity in range(mesh_num_facets)])
    entity_maps_mesh = {submesh: extract}

    return final_mesh, cell_tags, facet_tags, submesh, entity_maps_mesh



# Example usage:
if __name__ == "__main__":
    # Generate a cube mesh and retrieve mesh data
    final_mesh, cell_tags, facet_tags, submesh, entity_maps_mesh = generate_cube_gmshio(1.0, "cube.msh")

    # Define volume (dx) and surface (ds) measures
    dx = Measure("dx", domain=final_mesh)
    ds = Measure("ds", domain=final_mesh, subdomain_data=facet_tags)
    dx1 = Measure("dx", domain=submesh)

    # Define a Lagrange element of degree 3
    P1 = element("Lagrange", final_mesh.basix_cell(), 2)
    P = functionspace(final_mesh, P1)
    Q1 = element("Lagrange", submesh.basix_cell(), 2)
    Q = functionspace(submesh, Q1)


    # Define trial and test functions
    p, v = TrialFunction(P), TestFunction(P)
    q, u = TrialFunction(Q), TestFunction(Q)

    # Define the weak form
    a_00 = inner(p, v) * dx
    a_01 = inner(q, v) * ds(10)
    a_10 = inner(p, u.dx(0)) * ds(10)
    a_11 = inner(q, u) * dx1


    # Assemble the total form
    A_00 = form(a_00)
    A_01 = form(a_01, entity_maps=entity_maps_mesh)
    A_10 = form(a_10, entity_maps=entity_maps_mesh)
    A_11 = form(a_11)
    print("end")

provides this error :

libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2489:10: error: redefinition of 'J_c3'
 2489 |   double J_c3 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2487:10: note: previous definition is here
 2487 |   double J_c3 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2490:10: error: redefinition of 'J_c4'
 2490 |   double J_c4 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2486:10: note: previous definition is here
 2486 |   double J_c4 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2491:10: error: redefinition of 'J_c5'
 2491 |   double J_c5 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2488:10: note: previous definition is here
 2488 |   double J_c5 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2495:10: error: redefinition of 'J_c0'
 2495 |   double J_c0 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2484:10: note: previous definition is here
 2484 |   double J_c0 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2496:10: error: redefinition of 'J_c1'
 2496 |   double J_c1 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2483:10: note: previous definition is here
 2483 |   double J_c1 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2497:10: error: redefinition of 'J_c2'
 2497 |   double J_c2 = 0.0;
      |          ^
libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c:2485:10: note: previous definition is here
 2485 |   double J_c2 = 0.0;
      |          ^
6 errors generated.
Traceback (most recent call last):
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/spawn.py", line 87, in spawn
    subprocess.check_call(cmd, env=_inject_macos_ver(env))
  File "/opt/anaconda3/envs/clean777/lib/python3.10/subprocess.py", line 369, in check_call
    raise CalledProcessError(retcode, cmd)
subprocess.CalledProcessError: Command '['/opt/anaconda3/envs/clean777/bin/arm64-apple-darwin20.0.0-clang', '-ftree-vectorize', '-fPIC', '-fstack-protector-strong', '-O2', '-pipe', '-isystem', '/opt/anaconda3/envs/clean777/include', '-D_FORTIFY_SOURCE=2', '-isystem', '/opt/anaconda3/envs/clean777/include', '-I/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/ffcx/codegeneration', '-I/opt/anaconda3/envs/clean777/include/python3.10', '-c', 'libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.c', '-o', './libffcx_forms_9b64614a007c2a27dc174c07011f9d15d75b6f69.o', '-std=c17', '-O2', '-g0']' returned non-zero exit status 1.

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/compilers/C/unix.py", line 221, in _compile
    self.spawn(compiler_so + cc_args + [src, '-o', obj] + extra_postargs)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/compilers/C/base.py", line 1158, in spawn
    spawn(cmd, dry_run=self.dry_run, **kwargs)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/spawn.py", line 93, in spawn
    raise DistutilsExecError(
distutils.errors.DistutilsExecError: command '/opt/anaconda3/envs/clean777/bin/arm64-apple-darwin20.0.0-clang' failed with exit code 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/cffi/ffiplatform.py", line 48, in _build
    dist.run_command('build_ext')
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/dist.py", line 1104, in run_command
    super().run_command(command)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/dist.py", line 1021, in run_command
    cmd_obj.run()
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/command/build_ext.py", line 99, in run
    _build_ext.run(self)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/command/build_ext.py", line 368, in run
    self.build_extensions()
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/command/build_ext.py", line 484, in build_extensions
    self._build_extensions_serial()
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/command/build_ext.py", line 510, in _build_extensions_serial
    self.build_extension(ext)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/command/build_ext.py", line 264, in build_extension
    _build_ext.build_extension(self, ext)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/command/build_ext.py", line 565, in build_extension
    objects = self.compiler.compile(
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/compilers/C/base.py", line 655, in compile
    self._compile(obj, src, ext, cc_args, extra_postargs, pp_opts)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/setuptools/_distutils/compilers/C/unix.py", line 223, in _compile
    raise CompileError(msg)
distutils.compilers.C.errors.CompileError: command '/opt/anaconda3/envs/clean777/bin/arm64-apple-darwin20.0.0-clang' failed with exit code 1

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/pierremariotti/Documents/PhD/FEniCSx/new_try2/root/ABC_FEniCSx_classical/mwe_FEniCSx.py", line 135, in <module>
    A_10 = form(a_10, entity_maps=entity_maps_mesh)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 337, in form
    return _create_form(form)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 331, in _create_form
    return _form(form)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 254, in _form
    ufcx_form, module, code = jit.ffcx_jit(
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/dolfinx/jit.py", line 62, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/dolfinx/jit.py", line 212, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 225, in compile_forms
    raise e
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 205, in compile_forms
    impl = _compile_objects(
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 380, in _compile_objects
    ffibuilder.compile(tmpdir=cache_dir, verbose=True, debug=cffi_debug)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/cffi/api.py", line 725, in compile
    return recompile(self, module_name, source, tmpdir=tmpdir,
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/cffi/recompiler.py", line 1564, in recompile
    outputfilename = ffiplatform.compile('.', ext,
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/cffi/ffiplatform.py", line 20, in compile
    outputfilename = _build(tmpdir, ext, compiler_verbose, debug)
  File "/opt/anaconda3/envs/clean777/lib/python3.10/site-packages/cffi/ffiplatform.py", line 54, in _build
    raise VerificationError('%s: %s' % (e.__class__.__name__, e))
cffi.VerificationError: CompileError: command '/opt/anaconda3/envs/clean777/bin/arm64-apple-darwin20.0.0-clang' failed with exit code 1

The problem occurs when I have an coupled integral term that presents a spatial derivative of the test function.
The line a_10 = inner(p, u.dx(0)) * ds(10) leads to an error but the lines :

a_10 = inner(p, u) * ds(10)
a_10 = inner(p,dx(0), u) * ds(10)

doesn’t.
I know more or less why I have this problem, but my question is : Do you have any idea how I could implement this problematic term ?

Thank you,
Pierre.

2

This has been solved on main, ref: Resolve jacobians from duplicate meshes by jorgensd · Pull Request #733 · FEniCS/ffcx · GitHub
Here is a slightly modified version that runs on main:

from mpi4py import MPI
import gmsh

from dolfinx.io import gmshio

from ufl import (
    TestFunction,
    TrialFunction,
    TestFunctions,
    TrialFunctions,
    Measure,
    FacetNormal,
    CellNormal,
    inner,
    grad,
    div,
    Dn,
    MixedFunctionSpace,
    extract_blocks,
)
from basix.ufl import element
from dolfinx.fem import functionspace, form, petsc

import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio

import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio
import dolfinx.mesh as dmesh
import numpy as np


def generate_cube_gmshio(side_length=1.0, filename="cube.msh"):
    """
    Generate a cube with side length `side_length`,
    export the mesh to `filename`, and return:
      - final_mesh (full mesh)
      - cell_tags (tags for volume)
      - facet_tags (tags for surfaces)
      - submesh (mesh of one face)
      - entity_map (relation between submesh entities and parent mesh)
    """
    # Initialize the Gmsh API
    gmsh.initialize()

    comm = MPI.COMM_WORLD
    model_rank = 0
    model = gmsh.model
    model.add("cube")

    # Create the cube
    box_tag = model.occ.addBox(0, 0, 0, side_length, side_length, side_length)
    model.occ.synchronize()

    # Retrieve surfaces (dim=2)
    surfaces = model.occ.getEntities(dim=2)

    # Identify the face at z=0
    face_z0 = None
    for s in surfaces:
        com = model.occ.getCenterOfMass(s[0], s[1])
        if abs(com[2]) < 1e-8:  # z=0 face
            face_z0 = s[1]
            break
    if face_z0 is None:
        raise RuntimeError("Face at z=0 not found!")

    # Add physical groups
    model.addPhysicalGroup(3, [box_tag], tag=1)  # volume
    model.addPhysicalGroup(2, [face_z0], tag=10)  # surface z=0

    # Mesh options
    gmsh.option.setNumber("Mesh.ElementOrder", 2)
    gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
    model.mesh.generate(3)

    # Save mesh
    gmsh.write(filename)

    # Convert to FEniCSx
    # final_mesh, cell_tags, facet_tags = gmshio.model_to_mesh(model, comm, model_rank)
    mesh_data = gmshio.model_to_mesh(model, comm, model_rank)
    final_mesh = mesh_data.mesh
    cell_tags = mesh_data.cell_tags
    facet_tags = mesh_data.facet_tags
    gmsh.finalize()

    # --- Create submesh of face z=0 ---
    # Extract entities of the facet with tag=10
    facets = facet_tags.find(10)

    # Create the submesh (dim=2 since we want a surface)
    submesh, entity_map = dmesh.create_submesh(final_mesh, dim=2, entities=facets)[0:2]

    tdim = final_mesh.topology.dim
    fdim = tdim - 1
    submesh_tdim = submesh.topology.dim
    submesh_fdim = submesh_tdim - 1

    # Create the entity maps and an integration measure
    mesh_facet_imap = final_mesh.topology.index_map(fdim)
    mesh_num_facets = mesh_facet_imap.size_local + mesh_facet_imap.num_ghosts
    entity_maps_mesh = [entity_map]
    # extract = np.array(
    #     [
    #         entity_map.tolist().index(entity) if entity in entity_map else -1
    #         for entity in range(mesh_num_facets)
    #     ]
    # )
    # entity_maps_mesh = {submesh: extract}

    return final_mesh, cell_tags, facet_tags, submesh, entity_maps_mesh


# Example usage:
if __name__ == "__main__":
    # Generate a cube mesh and retrieve mesh data
    final_mesh, cell_tags, facet_tags, submesh, entity_maps_mesh = generate_cube_gmshio(
        1.0, "cube.msh"
    )

    # Define volume (dx) and surface (ds) measures
    dx = Measure("dx", domain=final_mesh)
    ds = Measure("ds", domain=final_mesh, subdomain_data=facet_tags)
    dx1 = Measure("dx", domain=submesh)

    # Define a Lagrange element of degree 3
    P1 = element("Lagrange", final_mesh.basix_cell(), 2)
    P = functionspace(final_mesh, P1)
    Q1 = element("Lagrange", submesh.basix_cell(), 2)
    Q = functionspace(submesh, Q1)

    # Define trial and test functions
    p, v = TrialFunction(P), TestFunction(P)
    q, u = TrialFunction(Q), TestFunction(Q)

    # Define the weak form
    a_00 = inner(p, v) * dx
    a_01 = inner(q, v) * ds(10)
    a_10 = inner(p, u.dx(0)) * ds(10)
    a_11 = inner(q, u) * dx1

    # Assemble the total form
    A_00 = form(a_00)
    A_01 = form(a_01, entity_maps=entity_maps_mesh)
    A_10 = form(a_10, entity_maps=entity_maps_mesh)
    A_11 = form(a_11)
    print("end")
3

Oh thanks ! I’ve already asked the question and the answer was no, but maybe it has changed since then : If I use FEniCSx in a conda environnement, as it is proposed on the FEniCSx website for MacOS, will I have access to the main branch of FEniCSx ? or will I have to download FEniCSx in pure local on my machine ? I’m not so good with those IT aspects…

4

Conda currently only supply stable releases.
There are several other ways of getting the main branch, such as

5

Many thanks ! As effective as always..