Compilation problem with mixed-dimensional PDE problem

dolfinx
1

Hello everyone,

I was trying to implement a mixed-dimensional problem. For a better understanding, I wanted to build a prototype example that solves a Poisson equation simultaneously on a 2d mesh and a 1d mesh consisting of one boundary edge. So far I am stuck at compiling the form of the mixed problem. Here’s the MWE:

import dolfinx as dfx
from mpi4py.MPI import COMM_WORLD as comm
import numpy as np
import ufl

# Define the base mesh
mesh = dfx.mesh.create_unit_square(comm, 24, 24)

# Define the submesh as bottom boundary of the mesh
bottom_tag = 42

def bottom_marker(x, tol=1e-6):
    return np.isclose(x[1], 0.0, atol=tol)

gdim = mesh.geometry.dim
tdim = mesh.topology.dim
fdim = tdim - 1

mesh.topology.create_connectivity(fdim, gdim)
facet_map = mesh.topology.index_map(fdim)
num_facets_local = facet_map.size_local + facet_map.num_ghosts

facet_values = np.zeros(num_facets_local, dtype=np.int32)

facet_values[dfx.mesh.locate_entities_boundary(mesh, fdim, bottom_marker)] = bottom_tag

all_facets_idx = np.arange(num_facets_local, dtype=np.int32)

facet_tags = dfx.mesh.meshtags(mesh, fdim, all_facets_idx, values=facet_values)

mesh0, mesh0_to_mesh = dfx.mesh.create_submesh(mesh, fdim, facet_tags.find(bottom_tag))[:2]

# Entity maps to be used for compiling the problem
facet_imap = mesh.topology.index_map(facet_tags.dim)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
mesh_to_mesh0 = np.full(num_facets, -1)
mesh_to_mesh0[mesh0_to_mesh] = np.arange(len(mesh0_to_mesh))
entity_maps = {mesh0: mesh_to_mesh0}

# Function spaces, solution and test functions
V = dfx.fem.functionspace(mesh, ("Lagrange", 1))
V0 = dfx.fem.functionspace(mesh0, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(V, V0)

u, u0 = dfx.fem.Function(V), dfx.fem.Function(V0)
v, v0 = ufl.TestFunctions(W)  # type: ignore

dx = ufl.Measure("dx", domain=mesh)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags,
                 subdomain_id=bottom_tag)  # type: ignore

# Define forms for each sub-problem
form1 = ufl.dot(ufl.grad(u), ufl.grad(v)) * dx + 4. * v * dx  # type: ignore

form2 = ufl.dot(ufl.grad(u0), ufl.grad(v0)) * ds + 4. * v0 * ds  # type: ignore

form = form1 + form2  # # type: ignore

# Compile the forms: Here's the problem
dfx.fem.form([form1, form2], entity_maps=entity_maps)

du, du0 = ufl.TrialFunctions(W)  # type: ignore

jac = ufl.derivative(form, u, du) + ufl.derivative(form, u0, du0)

J = dfx.fem.form(ufl.extract_blocks(jac), entity_maps=entity_maps)  #type: ignore

The code is mostly adopted from the two tutorials https://jsdokken.com/FEniCS-workshop/src/multiphysics/coupling.html and https://jsdokken.com/FEniCS-workshop/src/multiphysics/submeshes.html, but I wanted to adapt it to the case where the boundary is governed by PDE that gets some information of the bulk problem.

When running the code throws an error at the compilation stage:

libffcx_forms_55d25df265ff6b1ba7ec9f9d1142822678c2fb5d.c:619:8: error: redefinition of 'J_c0'
  619 | double J_c0 = 0.0;
      |        ^
libffcx_forms_55d25df265ff6b1ba7ec9f9d1142822678c2fb5d.c:617:8: note: previous definition is here
  617 | double J_c0 = 0.0;
      |        ^
libffcx_forms_55d25df265ff6b1ba7ec9f9d1142822678c2fb5d.c:620:8: error: redefinition of 'J_c1'
  620 | double J_c1 = 0.0;
      |        ^
libffcx_forms_55d25df265ff6b1ba7ec9f9d1142822678c2fb5d.c:618:8: note: previous definition is here
  618 | double J_c1 = 0.0;
      |        ^
2 errors generated.

Although there are a few related topics in this forum, I couldn’t make sense of them which is why I am asking. Can anybody give me some tipps what is wrong with the code?

Best,
Tom