# Bug either in legacy dolfin using dS or in pygmsh/gmsh mesh

**URL:** <https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030>\
**Category:** Legacy DOLFIN\
**Created:** [March 26, 2024, 2:01pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030 "2024-03-26T14:01:39Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![jacobkrich](https://avatars.discourse-cdn.com/v4/letter/j/f19dbf/32.png) [@jacobkrich](https://fenicsproject.discourse.group/u/jacobkrich)\
**Post date:** [March 26, 2024, 2:01pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030/1 "2024-03-26T14:01:39Z")

</div>

I do 2D Poisson/drift-diffusion semiconductor device modeling using legacy dolfin, as implemented in [Simudo](https://github.com/simudo/simudo). I have had a devil of a bug that I’ve spent months tracking down, and I have narrowed it down to the following simple MWE. The real problem I care about is a nonlinear problem on a mixed space; it occurs consistently when I use meshes that I generate with pygmsh/gmsh and does not occur when I use meshes that we make manually (by specifying explicit points/triangles).

With that said, here’s the issue. I make a very simple mesh, divided into three regions. In this case, I have only 4 triangles per region, and the figure shows the numbering of the cells.

 ![mesh](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/8e66ef9dbb40403e1dde0b6f2838d76e62626f26.jpeg)

I make a `dx` measure restricted to only some of the regions (I use 1,3 in the below example). If I add a term including `dS` into the weak form, one or more cells from one region behave as though they were in a different region.

My simple mesh is available here:  
[mesh.xdmf](https://www.dropbox.com/scl/fi/wsaaziy4xxd0dmy4w2zs7/mesh.xdmf?rlkey=rudqqpxpyx3j3f70xsq1gbjtl&dl=1)[mesh.h5](https://www.dropbox.com/scl/fi/wenh5xmiap4d1tvidvhkv/mesh.h5?rlkey=nm7rs7mpcq7y55qptwafw0r78&dl=1)  
And here is an MWE:

```python
import dolfin
import numpy as np

mesh = dolfin.Mesh()
with dolfin.XDMFFile("mesh.xdmf") as xdmf:
    xdmf.read(mesh)
    cell_markers = dolfin.MeshFunction('size_t', mesh, mesh.topology().dim())
    xdmf.read(cell_markers, "cv")

# Define integration measures on subdomains
dx = dolfin.dx(domain=mesh, subdomain_data=cell_markers)
dx_middle = dx(2)
dx_else = dx((1,3))
dS = dolfin.dS(domain=mesh) #internal facet measure

# Set up the DG0 space
DG0_space = dolfin.FunctionSpace(mesh, "DG", 0)
u = dolfin.Function(DG0_space)
v = dolfin.TestFunction(DG0_space)

# Here is our weak form. It's singular in this case, 
# so we're not going to try to solve it but just look at the assembled matrix.
F = u * v * dx_else

# The next term should have no effect at all. Not only is it
# multiplied by zero, but it contains no trial functions, so even
# if the constant were nonzero, it would still only contribute to the rhs vector b.
# But if it is included, the assembled matrix changes.
include_dS = True
if include_dS:
    F += dolfin.Constant(0.0) * v("+") * dS 

#This is a linear problem, but let's treat it as a nonlinear one:
J = dolfin.derivative(F,u)
bcs=[]
assembler = dolfin.fem.assembling.SystemAssembler(J,F,bcs)
A = dolfin.PETScMatrix()
assembler.assemble(A)
A_array = A.array()

#Find which dofs are in the middle region:
middle_region = [0.1,0.2] 
def region_dofs(space):
    dofmap = space.dofmap()
    region_dofs = []
    for cell in dolfin.cells(mesh):
        cell_dofs = dofmap.cell_dofs(cell.index())        
        cell_coords = cell.get_vertex_coordinates()
        cell_coords = np.reshape(cell_coords, (-1,2))        
        if np.all(middle_region[0] <= cell_coords[:,0]) and np.all(cell_coords[:,0]<=middle_region[1]):
            region_dofs.extend(cell_dofs)
    return list(set(region_dofs))

dofs = DG0_space.dofmap().dofs() # len 12
middle_dofs = region_dofs(DG0_space) # len 4

# Rows of the A matrix correspond to test functions v. 
# The test functions corresponding to the middle region (marked with a *) should 
# have no nonzeros. When include_dS is False, that happens correctly.
# When include_dS is True, however, two cells seem to swap places, and one cell in the middle
# region (#3) gets a nonzero while one cell in the outer region (#7) gets a zero
# Note that for other setups of this problem, the swapping has occurred between cells 2 and 7.
print(f"Include_dS: {include_dS}")
print(f"* marks cell in the middle region.")
for ind, i in enumerate(dofs):
    nonzero_row = np.nonzero(A_array[i,:])[0]
    extra = "*" if i in middle_dofs else ""
    print(f"Row corresponding to dof {i}{extra} has {len(nonzero_row)} nonzeros. Should be {0 if i in middle_dofs else 1}.")

```

When `include_dS` is `False`, the output is:

```python
Include_dS: False
* marks cell in the middle region.
Row corresponding to dof 0* has 0 nonzeros. Should be 0.
Row corresponding to dof 1* has 0 nonzeros. Should be 0.
Row corresponding to dof 2* has 0 nonzeros. Should be 0.
Row corresponding to dof 3* has 0 nonzeros. Should be 0.
Row corresponding to dof 4 has 1 nonzeros. Should be 1.
Row corresponding to dof 5 has 1 nonzeros. Should be 1.
Row corresponding to dof 6 has 1 nonzeros. Should be 1.
Row corresponding to dof 7 has 1 nonzeros. Should be 1.
Row corresponding to dof 8 has 1 nonzeros. Should be 1.
Row corresponding to dof 9 has 1 nonzeros. Should be 1.
Row corresponding to dof 10 has 1 nonzeros. Should be 1.
Row corresponding to dof 11 has 1 nonzeros. Should be 1.

```

As it should be.

When `include_dS` is `True`, however, the output should be identical but in fact is:

```python
Include_dS: True
* marks cell in the middle region.
Row corresponding to dof 0* has 0 nonzeros. Should be 0.
Row corresponding to dof 1* has 0 nonzeros. Should be 0.
Row corresponding to dof 2* has 1 nonzeros. Should be 0.
Row corresponding to dof 3* has 0 nonzeros. Should be 0.
Row corresponding to dof 4 has 1 nonzeros. Should be 1.
Row corresponding to dof 5 has 1 nonzeros. Should be 1.
Row corresponding to dof 6 has 1 nonzeros. Should be 1.
Row corresponding to dof 7 has 0 nonzeros. Should be 1.
Row corresponding to dof 8 has 1 nonzeros. Should be 1.
Row corresponding to dof 9 has 1 nonzeros. Should be 1.
Row corresponding to dof 10 has 1 nonzeros. Should be 1.
Row corresponding to dof 11 has 1 nonzeros. Should be 1.

```

Where you can see that row 2 and row 7 have exchanged. In this case, the `A` matrix is only 12x12, so one can inspect the whole thing and see that only the diagonal values are nonzeros, they are just in the wrong locations.

I am mystified. Is there a problem with this mesh that you can see? Or something that goes wrong when using both `dS` and domain-restricted `dx` in a weak form, but only with meshes from pygmsh/gmsh? I know that `dolfinx` is preferred at this point, but we have a large amount of code relying on legacy dolfin, so the switch would not be easy.

This problem also occurs without using `dolfin.derivative`. So, for example, letting

```python
u0 = dolfin.TrialFunction(DG0_space)
F0 = u0 * v * dx_else
F0 += dolfin.Constant(0.0) * v("+") * dS 
assembler = dolfin.fem.assembling.SystemAssembler(dolfin.lhs(F0),dolfin.rhs(F0),bcs)
A0 = dolfin.PETScMatrix()
assembler.assemble(A0)
A0_array = A0.array()

```

gives the same problem.

I am running dolfin 2019.2.0.dev0 installed with `apt` in an arm64 Ubuntu jammy docker container on Apple silicon. I have also reproduced this issue in a similar docker container on mantic running dolfin 2019.2.0.13.dev0.

Any help greatly appreciated!

---

<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:** [March 26, 2024, 2:57pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030/2 "2024-03-26T14:57:09Z")

</div>

> [@jacobkrich](#):
>
> ```auto
> dx = dolfin.dx(domain=mesh, subdomain_data=cell_markers)
> ...
> dS = dolfin.dS(domain=mesh) #internal facet measure
> 
> ```

Legacy dolfin used to rely on subdomain numbering to assign + and - sides in dS. What happens if you use

```auto
dS = dolfin.dS(domain=mesh, subdomain_data=cell_markers) #internal facet measure

```

instead?

---

<div class="post-metadata">

**Author:** ![jacobkrich](https://avatars.discourse-cdn.com/v4/letter/j/f19dbf/32.png) [@jacobkrich](https://fenicsproject.discourse.group/u/jacobkrich)\
**Post date:** [March 26, 2024, 3:11pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030/3 "2024-03-26T15:11:38Z")

</div>

Thank you for the suggestion. I have tried that amendment and the result is unchanged, though.

---

<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:** [March 26, 2024, 3:18pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030/4 "2024-03-26T15:18:08Z")

</div>

Could be related to:

> **[Assembly of integrals over interal boundaries messes up DirichletBCs](https://bitbucket.org/fenics-project/dolfin/issues/1075/assembly-of-integrals-over-interal)**
>
> In forms that contain integral over \*\*internal boundaries\*\* where also DirichletBCs are applied,   
> the solution does not agree with the given DirichletBC, when the assembly of matrix, rhs-vector and BCs is   
> done via \`assemble\_system\` or some...

I guess you are in a tricky situation, as there are no-one fixing bugs in legacy dolfin. Ive fixed a few over the last few years, relating to my research.

If you can track down what is the cause of the error (does it only happen on imported meshes?)  
and fix the issue i am happy to merge it into the master branch. However, it is unlikely that I have the time and bandwidth to find the underlying cause and add a fix myself.

---

<div class="post-metadata">

**Author:** ![jacobkrich](https://avatars.discourse-cdn.com/v4/letter/j/f19dbf/32.png) [@jacobkrich](https://fenicsproject.discourse.group/u/jacobkrich)\
**Post date:** [March 26, 2024, 4:35pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030/5 "2024-03-26T16:35:51Z")

</div>

Thank you for pointing me to that discussion, which does seem to be the same issue and points the finger at the SystemAssembler. I am not very good with C++, but I’ll see what I can do.

I assume the issue arises in `facet_wise_assembly` (invoked when `dS` is present) as opposed to `cell_wise_assembly`. It seems that some of the cell dofs get mixed up in there, somewhere, sometimes.

Does dolfinx have an equivalent of SystemAssembler, which I could possibly look at for reference?

---

<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:** [March 26, 2024, 7:53pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030/6 "2024-03-26T19:53:29Z")

</div>

The dolfinx code is quite different. The Closest you would get to systemassembler is the standard `assemble_matrix` and `assemble_vector` codes (maybe with `apply lifting`, as `assemble_system` in legacy dolfin does what is known as the lifting procedure, described in detail at  
[http://jsdokken.com/FEniCS23-tutorial/src/lifting.html#lifting](http://jsdokken.com/FEniCS23-tutorial/src/lifting.html#lifting)

---

<div class="post-metadata">

**Author:** ![jacobkrich](https://avatars.discourse-cdn.com/v4/letter/j/f19dbf/32.png) [@jacobkrich](https://fenicsproject.discourse.group/u/jacobkrich)\
**Post date:** [April 4, 2024, 9:23pm UTC](https://fenicsproject.discourse.group/t/bug-either-in-legacy-dolfin-using-ds-or-in-pygmsh-gmsh-mesh/14030/7 "2024-04-04T21:23:39Z")

</div>

I found the error in SystemAssembler and fixed it. I made a [pull request](https://bitbucket.org/fenics-project/dolfin/pull-requests/576), which describes the issue.