# Solving two Poisson problems on subdomains of a parent mesh yielding non-deterministic success

**URL:** <https://fenicsproject.discourse.group/t/solving-two-poisson-problems-on-subdomains-of-a-parent-mesh-yielding-non-deterministic-success/19661>\
**Category:** General\
**Created:** [April 18, 2026, 7:20pm UTC](https://fenicsproject.discourse.group/t/solving-two-poisson-problems-on-subdomains-of-a-parent-mesh-yielding-non-deterministic-success/19661 "2026-04-18T19:20:39Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [April 18, 2026, 7:20pm UTC](https://fenicsproject.discourse.group/t/solving-two-poisson-problems-on-subdomains-of-a-parent-mesh-yielding-non-deterministic-success/19661/1 "2026-04-18T19:20:39Z")

</div>

I’m trying to solve a pretty straightforward problem:

Find u\_0 \in V\_0(\Omega\_0) and u\_1 \in V\_1(\Omega\_1) such that

\begin{align} (\nabla u\_0, \nabla v\_0)\_{\Omega\_0} &= (1, v\_0)\_{\Omega\_0} \\ (\nabla u\_1, \nabla v\_1)\_{\Omega\_1} &= (1, v\_1)\_{\Omega\_1} \end{align}

for all v\_0 \in V\_{0,E}(\Omega\_0) and v\_1 \in V\_{1,E}(\Omega\_1). Here V\_{\cdot, E}(\cdot) corresponds to the space modified to accommodate homogeneous Dirichlet boundary data.

The following code is my attempt to approximate the solution of the above problem.

```Python
from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import numpy as np
import ufl

# Parent mesh
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)

# Extract submeshes, mesh0 left, mesh1 right
cells = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)
sub_cells0 = cells[dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, cells)[:, 0] <= 0.5]
sub_cells1 = cells[dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, cells)[:, 0] >= 0.5]
sub_mesh0, mapping_entity0, mapping_vertex0, mapping_geom0 = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, sub_cells0)
sub_mesh1, mapping_entity1, mapping_vertex1, mapping_geom1 = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, sub_cells1)

# Create domain markers
sub_mt0 = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, sub_cells0, np.ones_like(sub_cells0))
sub_mt1 = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, sub_cells1, np.ones_like(sub_cells1))

# Create spaces on each subdomain and mixed element
V0 = dolfinx.fem.functionspace(sub_mesh0, ("CG", 1))
V1 = dolfinx.fem.functionspace(sub_mesh1, ("CG", 1))
W = ufl.MixedFunctionSpace(V0, V1)

# Integration measures on each subdomain
dx0 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt0)
dx1 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt1)

# FE formulation
u0, u1 = ufl.TrialFunctions(W)
v0, v1 = ufl.TestFunctions(W)

a = ufl.inner(ufl.grad(u0), ufl.grad(v0)) * dx0
a += ufl.inner(ufl.grad(u1), ufl.grad(v1)) * dx1

L = ufl.inner(dolfinx.fem.Constant(sub_mesh0, 1.0), v0) * dx0
L += ufl.inner(dolfinx.fem.Constant(sub_mesh1, 1.0), v1) * dx1

# Extract block system
a_blocked = ufl.extract_blocks(a)
L_blocked = ufl.extract_blocks(L)

# Create BCs
sub_mesh0.topology.create_connectivity(sub_mesh0.topology.dim, sub_mesh0.topology.dim - 1)
sub_mesh0.topology.create_connectivity(sub_mesh0.topology.dim - 1, sub_mesh0.topology.dim)
u0_bdry = dolfinx.fem.Function(V0)
u0_bc = dolfinx.fem.dirichletbc(
    u0_bdry, dolfinx.fem.locate_dofs_topological(
        V0, sub_mesh0.topology.dim - 1, dolfinx.mesh.exterior_facet_indices(sub_mesh0.topology)
))

sub_mesh1.topology.create_connectivity(sub_mesh1.topology.dim, sub_mesh1.topology.dim - 1)
sub_mesh1.topology.create_connectivity(sub_mesh1.topology.dim - 1, sub_mesh1.topology.dim)
u1_bdry = dolfinx.fem.Function(V1)
u1_bc = dolfinx.fem.dirichletbc(
    u1_bdry, dolfinx.fem.locate_dofs_topological(
        V1, sub_mesh1.topology.dim - 1, dolfinx.mesh.exterior_facet_indices(sub_mesh1.topology)
))

# Solve problem
bcs = [u0_bc, u1_bc]
u0h = dolfinx.fem.Function(V0)
u1h = dolfinx.fem.Function(V1)
problem = dolfinx.fem.petsc.LinearProblem(
    a_blocked,
    L_blocked,
    u=[u0h, u1h],
    bcs=bcs,
    entity_maps=[mapping_entity0, mapping_entity1],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
        "ksp_monitor": None,
    },
    petsc_options_prefix="monolithic_",
)
problem.solve()

```

I notice that I’m getting non-deterministic success. Sometimes this solves giving me the expected result, and sometimes it does not, yielding the error:

```auto
Traceback (most recent call last):
  File "$SCRIPT_PATH/poisson_subdomain_coupled.py", line 65, in <module>
    problem = dolfinx.fem.petsc.LinearProblem(
        a_blocked,
    ...<11 lines>...
        petsc_options_prefix="monolithic_",
    )
  File "$LIB_PATH/python3.14/site-packages/dolfinx/fem/petsc.py", line 804, in __init__
    self._A = create_matrix(self._a, kind=kind)
              ~~~~~~~~~~~~~ ^^^^^^^^^^^^^^^^^^^^
  File "$LIB_PATH/python3.14/site-packages/dolfinx/fem/petsc.py", line 193, in create_matrix
    return _cpp.fem.petsc.create_matrix_block(_a, kind)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ^^^^^^^^^^
RuntimeError: Cannot insert rows that do not exist in the IndexMap.
Exception ignored while calling deallocator <function LinearProblem. __del__ at 0x12aee5380>:
Traceback (most recent call last):
  File "$LIB_PATH/python3.14/site-packages/dolfinx/fem/petsc.py", line 882, in __del__
    lambda obj: obj is not None, (self._solver, self._A, self._b, self._x, self._P_mat)
AttributeError: 'LinearProblem' object has no attribute '_solver'

```

Am I doing something silly in my original script?

```auto
>>> import dolfinx
>>> dolfinx.git_commit_hash
'338f32f692a1500571c910a29e69a3ff1d5c6549'
>>> dolfinx. __version__'0.11.0.dev0'

```

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [April 18, 2026, 10:31pm UTC](https://fenicsproject.discourse.group/t/solving-two-poisson-problems-on-subdomains-of-a-parent-mesh-yielding-non-deterministic-success/19661/2 "2026-04-18T22:31:38Z")

</div>

I think I may have poorly defined the subdomain integration measures. Making the following change

```Python
# Integration measures on each subdomain
dx0 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt0)
dx1 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt1)

```

```Python
# Integration measures on each subdomain
dx0 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt0, subdomain_id=1)
dx1 = ufl.Measure("dx", domain=mesh, subdomain_data=sub_mt1, subdomain_id=1)

```

Seems to give me consistently good results

---

<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:** [April 21, 2026, 3:30pm UTC](https://fenicsproject.discourse.group/t/solving-two-poisson-problems-on-subdomains-of-a-parent-mesh-yielding-non-deterministic-success/19661/3 "2026-04-21T15:30:28Z")

</div>

The first definition should in theory just integrate over all cells of the submesh. (As it should ignore all cells not in the submap.)  
It is a bit suprising it doesn’t work consistently.

To me it is a bit unclear why you are not using dx0 = ufl.Measure(“dx”, domain=submesh0)  
And similarly for dx1 as there is no coupling between the two submeshes. Then you wouldn’t need the meshtags.