# Solve poisson equation with three submeshes

**URL:** <https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841>\
**Category:** dolfinx\
**Created:** [May 30, 2024, 9:49am UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841 "2024-05-30T09:49:23Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [May 30, 2024, 9:49am UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/1 "2024-05-30T09:49:23Z")

</div>

Hi all,

I’m trying to leverage the entity\_maps introduced in dolfinx 0.8.0 to solve a poisson equation on a unit square but with three submeshes codim 0.

Ultimately, I would like to model interface discontinuities using a DG method but have CG elements within the domains to obtain this:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/fdf30bcad2629360158712748713f3d561d76b99.png)

For now, I’m simply trying to understand how to solve a poisson equation (without these discontinuities) on three submeshes.

This is the code that I have:

```python
from dolfinx import mesh, fem
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from mpi4py import MPI
import ufl
from ufl import grad, dot, jump, avg
import numpy as np
from dolfinx.nls.petsc import NewtonSolver

msh = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)
V = fem.functionspace(msh, ("DG", 1))

tol = 1e-16

def left_domain(x):
    return x[0] < 0.5 + tol

def top_right_domain(x):
    return np.logical_and(x[0] > 0.5 - tol, x[1] > 0.5 - tol)

def bottom_right_domain(x):
    return np.logical_and(x[0] > 0.5 - tol, x[1] < 0.5 + tol)

tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
c_to_f = msh.topology.connectivity(tdim, tdim - 1)
f_to_c = msh.topology.connectivity(tdim - 1, tdim)
facet_map = msh.topology.index_map(tdim - 1)
cell_imap = msh.topology.index_map(tdim)
left_cells = mesh.locate_entities(msh, tdim, left_domain)
top_right_cells = mesh.locate_entities(msh, tdim, top_right_domain)
bottom_right_cells = mesh.locate_entities(msh, tdim, bottom_right_domain)

left_submesh, left_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
    msh, msh.topology.dim, left_cells
)
top_right_submesh, top_right_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
    msh, msh.topology.dim, top_right_cells
)
bottom_right_submesh, bottom_right_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
    msh, msh.topology.dim, bottom_right_cells
)

# Create three function spaces on each submesh
V_left = fem.functionspace(left_submesh, ("CG", 1))
V_top_right = fem.functionspace(top_right_submesh, ("CG", 1))
V_bottom_right = fem.functionspace(bottom_right_submesh, ("CG", 1))

# create functions and test functions
u_left = fem.Function(V_left)
v_left = ufl.TestFunction(V_left)
u_top_right = fem.Function(V_top_right)
v_top_right = ufl.TestFunction(V_top_right)
u_bottom_right = fem.Function(V_bottom_right)
v_bottom_right = ufl.TestFunction(V_bottom_right)

# create dx measures
cell_map = msh.topology.index_map(msh.topology.dim)
num_cells = cell_map.size_local + cell_map.num_ghosts
marker = np.arange(num_cells, dtype=np.int32)
mt = dolfinx.mesh.meshtags(msh, msh.topology.dim, marker, marker)

dx = ufl.Measure("dx", domain=msh, subdomain_data=mt)
dx_left = ufl.Measure("dx", domain=left_submesh)
dx_top_right = ufl.Measure("dx", domain=top_right_submesh)
dx_bottom_right = ufl.Measure("dx", domain=bottom_right_submesh)

# Define variational problem
F = 0
F += dot(grad(v_left), grad(u_left)) * dx_left
F += dot(grad(v_top_right), grad(u_top_right)) * dx_top_right
F += dot(grad(v_bottom_right), grad(u_bottom_right)) * dx_bottom_right

# source
f = 1
F += -v_left * f * dx_left
F += -v_top_right * f * dx_top_right
F += -v_bottom_right * f * dx_bottom_right

# boundary conditions
boundary_facets = mesh.exterior_facet_indices(msh.topology)
boundary_dofs = fem.locate_dofs_topological(V, tdim - 1, boundary_facets)
uD = fem.Function(V)
uD.interpolate(lambda x: np.full(x[0].shape, 0.0))
bc = fem.dirichletbc(uD, boundary_dofs)

# TODO add dS coupling terms for interfaces

# make entity maps
parent_to_sub_left = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_left[left_sub_to_parent] = np.arange(
    len(left_sub_to_parent), dtype=np.int32
)
entity_map = {left_submesh._cpp_object: parent_to_sub_left}

parent_to_sub_right = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_right[top_right_sub_to_parent] = np.arange(
    len(top_right_sub_to_parent), dtype=np.int32
)
entity_map.update({top_right_submesh._cpp_object: parent_to_sub_right})

parent_to_sub_bottom = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_bottom[bottom_right_sub_to_parent] = np.arange(
    len(bottom_right_sub_to_parent), dtype=np.int32
)
entity_map.update({bottom_right_submesh._cpp_object: parent_to_sub_bottom})

# NonLinearProblem
formulation = dolfinx.fem.form(F == 0, entity_maps=entity_map)

u = fem.Function(V)
problem = fem.petsc.NonlinearProblem(
    formulation,
    u,
    bcs=bc,
)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.solve(u)

```

At the moment the code doesn’t run and produces the following error:

```auto
Traceback (most recent call last):
  File "/home/remidm/FESTIM/interface_submeshes.py", line 118, in <module>
    problem = fem.petsc.NonlinearProblem(
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 905, in __init__
    J = ufl.derivative(F, u, du)
        ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/ufl/formoperators.py", line 404, in derivative
    raise ValueError(f"Invalid argument type {type(form)}.")
ValueError: Invalid argument type <class 'ufl.equation.Equation'>.

```

I couldn’t find a demo for this application case so I was wondering if the approach was at least the right one?

Thanks for any pointers!

Cheers  
Rem

---

<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:** [May 30, 2024, 11:19am UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/2 "2024-05-30T11:19:30Z")

</div>

Using non-linear problem for mixed dimensional methods is not stream-lined. Here is a NewtonSolver that I made recently for mixed dim problems:

```python
# Newton solver for mixed dim problems
# Author: Jørgen S. Dokken
# SPDX License identifier: MIT
class NewtonSolver():
    max_iterations: int
    bcs: list[dolfinx.fem.DirichletBC]
    A: PETSc.Mat
    b: PETSc.Vec
    J: dolfinx.fem.Form
    b: dolfinx.fem.Form
    dx: PETSc.Vec
    def __init__ (self, F:list[dolfinx.fem.form], J:list[list[dolfinx.fem.form]], w: list[dolfinx.fem.Function], 
                 bcs: list[dolfinx.fem.DirichletBC]|None=None, max_iterations:int=5,
                 petsc_options: dict[str, str|float|int|None]=None,
                 problem_prefix = "newton"):
        self.max_iterations = max_iterations
        self.bcs = [] if bcs is None else bcs
        self.b = dolfinx.fem.petsc.create_vector_block(F)
        self.F = F
        self.J = J
        self.A = dolfinx.fem.petsc.create_matrix_block(J)
        self.dx = self.A.createVecLeft()
        self.w = w
        self.x = dolfinx.fem.petsc.create_vector_block(F)

        # Set PETSc options
        opts = PETSc.Options()
        opts.prefixPush(problem_prefix)
        if petsc_options is not None:
            for k, v in petsc_options.items():
                opts[k] = v
        opts.prefixPop()

        # Define KSP solver    
        self._solver = PETSc.KSP().create(self.b.getComm().tompi4py())
        self._solver.setOperators(self.A)
        self._solver.setFromOptions()

        # Set matrix and vector PETSc options
        self.A.setOptionsPrefix(problem_prefix)
        self.A.setFromOptions()
        self.b.setOptionsPrefix(problem_prefix)
        self.b.setFromOptions()

   
    def solve(self, tol=1e-6, beta=1.0):
        i = 0

        while i < self.max_iterations:
            dolfinx.cpp.la.petsc.scatter_local_vectors(
                self.x,
                [si.x.petsc_vec.array_r for si in self.w],
                [
                    (si.function_space.dofmap.index_map, si.function_space.dofmap.index_map_bs)
                    for si in self.w
                ])
            self.x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

            # Assemble F(u_{i-1}) - J(u_D - u_{i-1}) and set du|_bc= u_D - u_{i-1}
            self.b.zeroEntries()
            dolfinx.fem.petsc.assemble_vector_block(self.b, self.F,self.J, bcs=self.bcs,x0=self.x, scale=-1.0)
            self.b.ghostUpdate(PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.FORWARD)
            # Assemble Jacobian
            self.A.zeroEntries()
            dolfinx.fem.petsc.assemble_matrix_block(self.A, self.J, bcs=self.bcs)
            self.A.assemble()
            
            
            self._solver.solve(self.b, self.dx)
            assert self._solver.getConvergedReason() > 0, "Linear solver did not converge"
            offset_start = 0
            for s in self.w:
                num_sub_dofs = s.function_space.dofmap.index_map.size_local * s.function_space.dofmap.index_map_bs
                s.x.array[:num_sub_dofs] -= beta*self.dx.array_r[offset_start:offset_start+num_sub_dofs]
                s.x.scatter_forward()
                offset_start += num_sub_dofs
            # Compute norm of update

            correction_norm = self.dx.norm(0)
            print(f"Iteration {i}: Correction norm {correction_norm}")
            if correction_norm < tol:
                break
            i += 1

```

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [May 30, 2024, 11:50am UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/3 "2024-05-30T11:50:14Z")

</div>

Thanks for this!

I’ve modified the code to the following:

```python
formulation = dolfinx.fem.form(F == 0, entity_maps=entity_map)

u = fem.Function(V)

solver = NewtonSolver(
    F=formulation, w=[u_left, u_top_right, u_bottom_right], bcs=bc, J=[] # not sure how to compute J here
)

solver.solve(u)

```

But now have this error:

```auto
Traceback (most recent call last):
  File "/home/remidm/FESTIM/interface_submeshes.py", line 125, in <module>
    solver = NewtonSolver(
             ^^^^^^^^^^^^^
  File "/home/remidm/FESTIM/newton_solver_mixed_dim.py", line 29, in __init__
    self.b = dolfinx.fem.petsc.create_vector_block(F)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 124, in create_vector_block
    maps = [
           ^
TypeError: 'Equation' object is not iterable

```

I thought maybe the `F==0` in `form` was causing the issue but when only setting `dolfinx.fem.form(F, entity_maps=entity_map)` I get:

```auto
Traceback (most recent call last):
  File "/home/remidm/FESTIM/interface_submeshes.py", line 117, in <module>
    formulation = dolfinx.fem.form(F, entity_maps=entity_map)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/dolfinx_0.8/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 178, in _form
    (domain,) = list(sd.keys()) # Assuming single domain
    ^^^^^^^^^
ValueError: too many values to unpack (expected 1)

```

---

<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:** [May 30, 2024, 12:22pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/4 "2024-05-30T12:22:39Z")

</div>

You would not make an equation in that sense, you would do something like:

```python
F0 = F00 + F01 
F1 = F10 + F11
residual_0 = dolfinx.fem.form(F0, entity_maps=entity_maps)
residual_1 = dolfinx.fem.form(F1, entity_maps=entity_maps)
jac00 = ufl.derivative(F0, u)
jac01 = ufl.derivative(F0, psi)
jac10 = ufl.derivative(F1, u)
jac11 = ufl.derivative(F1, psi)
J00 = dolfinx.fem.form(jac00, entity_maps=entity_maps)
J01 = dolfinx.fem.form(jac01, entity_maps=entity_maps)
J10 = dolfinx.fem.form(jac10, entity_maps=entity_maps)
J11 = dolfinx.fem.form(jac11, entity_maps=entity_maps)

J = [[J00, J01], [J10, J11]]
F = [residual_0, residual_1]

```

For a two-domain problem

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [May 30, 2024, 12:42pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/5 "2024-05-30T12:42:40Z")

</div>

It’s unclear to me what `F00` is with respect to `F01`

Moreover, is `psi` the test function?

is `u` a function defined over the whole domain not only a submesh?

---

<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:** [May 30, 2024, 2:10pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/6 "2024-05-30T14:10:04Z")

</div>

Hi Remi,

Sorry for the lack of clarity (my example is a bit experimental, using an exotic variational formulation Im not ready to share yet).  
However, I can give you the context:  
In my example `u` is a function on the whole mesh, while `psi` is a function on the subset of the boundary (using Joe’s meshview branches).  
So `u` and `psi` are the two unknowns

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [May 30, 2024, 2:14pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/7 "2024-05-30T14:14:21Z")

</div>

After a bit of trial and error this is where I am.  
I’ve simplified the problem to only a two domain problem (left and right).

- I created a new Function `psi` since the test function wasn’t accepted. But maybe I did this wrong…
- It’s unclear to me if I need to create functions on the submeshes and test functions on the whole mesh, or the opposite.

This is my code (apologies for the length…):

```python
from dolfinx import mesh, fem
import dolfinx.fem.petsc
import dolfinx.nls.petsc
from mpi4py import MPI
import ufl
from ufl import grad, dot, jump, avg
import numpy as np

from newton_solver_mixed_dim import NewtonSolver

msh = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100)
V = fem.functionspace(msh, ("DG", 1))

tol = 1e-16

def left_domain(x):
    return x[0] < 0.5 + tol

def right_domain(x):
    return x[0] > 0.5 - tol

def interface(x):
    return np.isclose(x[0], 0.5)

tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
c_to_f = msh.topology.connectivity(tdim, tdim - 1)
f_to_c = msh.topology.connectivity(tdim - 1, tdim)
facet_map = msh.topology.index_map(tdim - 1)
cell_imap = msh.topology.index_map(tdim)
left_cells = mesh.locate_entities(msh, tdim, left_domain)
right_cells = mesh.locate_entities(msh, tdim, right_domain)
boundary_facets = mesh.exterior_facet_indices(msh.topology)
interface_facets = mesh.locate_entities(msh, fdim, interface)

num_facets = facet_map.size_local + facet_map.num_ghosts
indices = np.arange(0, num_facets)
values = np.zeros(indices.shape, dtype=np.intc) # all facets are tagged with zero

interface_id = 2
values[boundary_facets] = 1
values[interface_facets] = interface_id

mesh_tags_facets = mesh.meshtags(msh, fdim, indices, values)

left_submesh, left_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
    msh, msh.topology.dim, left_cells
)
right_submesh, right_sub_to_parent, _, _ = dolfinx.mesh.create_submesh(
    msh, msh.topology.dim, right_cells
)

# Create three function spaces on each submesh
V_left = fem.functionspace(left_submesh, ("CG", 1))
V_right = fem.functionspace(right_submesh, ("CG", 1))

# create functions and test functions
u_left = fem.Function(V_left)
v_left = ufl.TestFunction(V_left)
u_right = fem.Function(V_right)
v_right = ufl.TestFunction(V_right)
u = fem.Function(V)
v = ufl.TestFunction(V)

# create measures
cell_map = msh.topology.index_map(msh.topology.dim)
num_cells = cell_map.size_local + cell_map.num_ghosts
indices = np.arange(0, num_cells)
values_cells = np.zeros(indices.shape, dtype=np.intc)
index_left = 1
index_right = 2
values_cells[left_cells] = index_left
values_cells[right_cells] = index_right

mesh_tags_cells = mesh.meshtags(msh, tdim, indices, values_cells)
mt = dolfinx.mesh.meshtags(msh, msh.topology.dim, indices, values_cells)
dx = ufl.Measure("dx", domain=msh, subdomain_data=mt)
dx_left = ufl.Measure("dx", domain=left_submesh)
dx_right = ufl.Measure("dx", domain=right_submesh)
dS = ufl.Measure("dS", domain=msh, subdomain_data=mesh_tags_facets)

# Define variational problem
F0 = dot(grad(v), grad(u_left)) * dx(index_left)
F1 = dot(grad(v), grad(u_right)) * dx(index_right)

# source
f = 1
F0 += -v * f * dx(index_left)
F1 += -v * f * dx(index_right)

# boundary conditions
boundary_facets = mesh.exterior_facet_indices(msh.topology)
boundary_dofs = fem.locate_dofs_topological(V, tdim - 1, boundary_facets)
uD = fem.Function(V)
uD.interpolate(lambda x: np.full(x[0].shape, 0.0))
bc = fem.dirichletbc(uD, boundary_dofs)

# TODO add dS coupling terms for interfaces
h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)
alpha = 1000
F0 += (
    -dot(avg(grad(v)), jump(u, n)) * dS(interface_id)
    - dot(jump(v, n), avg(grad(u))) * dS(interface_id)
    + alpha / avg(h) * dot(jump(v, n), jump(u, n)) * dS(interface_id)
)

# make entity maps
parent_to_sub_left = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_left[left_sub_to_parent] = np.arange(
    len(left_sub_to_parent), dtype=np.int32
)

parent_to_sub_right = np.full(num_cells, -1, dtype=np.int32)
parent_to_sub_right[right_sub_to_parent] = np.arange(
    len(right_sub_to_parent), dtype=np.int32
)

entity_maps = {
    left_submesh._cpp_object: parent_to_sub_left,
    right_submesh._cpp_object: parent_to_sub_right,
}

# NonLinearProblem
residual_0 = dolfinx.fem.form(F0, entity_maps=entity_maps)
residual_1 = dolfinx.fem.form(F1, entity_maps=entity_maps)

psi = dolfinx.fem.Function(V)

jac00 = ufl.derivative(F0, u)
jac01 = ufl.derivative(F0, psi)
jac10 = ufl.derivative(F1, u)
jac11 = ufl.derivative(F1, psi)
J00 = dolfinx.fem.form(jac00, entity_maps=entity_maps)
J01 = dolfinx.fem.form(jac01, entity_maps=entity_maps)
J10 = dolfinx.fem.form(jac10, entity_maps=entity_maps)
J11 = dolfinx.fem.form(jac11, entity_maps=entity_maps)

J = [[J00, J01], [J10, J11]]

solver = NewtonSolver(F=[residual_0, residual_1], w=[u_left, u_right], bcs=[bc], J=J)

solver.solve()

```

Giving:

```auto
Traceback (most recent call last):
  File "/home/remidm/FESTIM/interface_submeshes.py", line 148, in <module>
    solver.solve()
  File "/home/remidm/FESTIM/newton_solver_mixed_dim.py", line 88, in solve
    self._solver.solve(self.b, self.dx)
  File "petsc4py/PETSc/KSP.pyx", line 1774, in petsc4py.PETSc.KSP.solve
petsc4py.PETSc.Error: error code 73
[0] KSPSolve() at /home/conda/feedstock_root/build_artifacts/petsc_1717048443064/work/src/ksp/ksp/interface/itfunc.c:1078
[0] KSPSolve_Private() at /home/conda/feedstock_root/build_artifacts/petsc_1717048443064/work/src/ksp/ksp/interface/itfunc.c:831
[0] KSPSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1717048443064/work/src/ksp/ksp/interface/itfunc.c:415
[0] PCSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1717048443064/work/src/ksp/pc/interface/precon.c:1079
[0] PCSetUp_ILU() at /home/conda/feedstock_root/build_artifacts/petsc_1717048443064/work/src/ksp/pc/impls/factor/ilu/ilu.c:135
[0] MatILUFactorSymbolic() at /home/conda/feedstock_root/build_artifacts/petsc_1717048443064/work/src/mat/interface/matrix.c:7036
[0] MatILUFactorSymbolic_SeqAIJ() at /home/conda/feedstock_root/build_artifacts/petsc_1717048443064/work/src/mat/impls/aij/seq/aijfact.c:1745
[0] Object is in wrong state
[0] Matrix is missing diagonal entry 0

```

So obviously something is wrong ☹

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [May 30, 2024, 3:50pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/8 "2024-05-30T15:50:49Z")

</div>

Should be pretty much straightforward to with multiphenicsx, considering that the third tutorial (interface case) basically has the same structure.

However, the downside is that you would be installing an additional library, rather than using Joe’s functionalities which are built-in in dolfinx.

---

<div class="post-metadata">

**Author:** ![Caglar](https://avatars.discourse-cdn.com/v4/letter/c/77aa72/32.png) [@Caglar](https://fenicsproject.discourse.group/u/Caglar)\
**Post date:** [April 24, 2025, 2:53pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/9 "2025-04-24T14:53:34Z")

</div>

Hey, have you found a solution to this “missing diagonal entry”?

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [April 24, 2025, 3:09pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/10 "2025-04-24T15:09:16Z")

</div>

I don’t remember anything about this exact problem but with successfully implemented interface discontinuities using the mixed domain approach and entity maps. Ping @dokken

---

<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:** [May 6, 2025, 6:59pm UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/11 "2025-05-06T18:59:11Z")

</div>

There are several variations of the Poisson problem with coupling.  
One of them is within FESTIM itself ([GitHub - RemDelaporteMathurin/FESTIM: My personal fork of FESTIM](https://github.com/RemDelaporteMathurin/festim)) and might be hard to parse.  
We have several other versions floating around on discourse, including:  
[Discontinuity at interface using mixed domains - #2 by dokken](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/2)  
and the subsequent typo fix by @RemDelaporteMathurin  
as well as:  
[Poisson submesh DG coupling · GitHub](https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py)  
(which is also cross posted multiple times on discourse: [Search results for 'https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce' - FEniCS Project](https://fenicsproject.discourse.group/search?q=https%3A%2F%2Fgist.github.com%2Fjorgensd%2F8a5c32f491195e838f5863ca88b27bce))

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [May 8, 2025, 6:53am UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/12 "2025-05-08T06:53:04Z")

</div>

Hello dokken, your comment worries me a little, as I use the standard Newton solver without any visible issue solving a mixed elements problem.

Do you mean that the solution I get might not be the correct one? Or do you mean that it is hard to reach convergence or something else?

---

<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:** [May 8, 2025, 7:49am UTC](https://fenicsproject.discourse.group/t/solve-poisson-equation-with-three-submeshes/14841/13 "2025-05-08T07:49:15Z")

</div>

`MixedElement` is not the same as `MixedFunctionSpace`.

You shouldn’t be worried if the code runs and gives a solution for all variables.