# Integral of inner surface wrongly evaluate to zero in 3D

**URL:** <https://fenicsproject.discourse.group/t/integral-of-inner-surface-wrongly-evaluate-to-zero-in-3d/8031>\
**Category:** Uncategorized\
**Created:** [April 12, 2022, 9:03am UTC](https://fenicsproject.discourse.group/t/integral-of-inner-surface-wrongly-evaluate-to-zero-in-3d/8031 "2022-04-12T09:03:54Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![RaphaelC](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/raphaelc/32/3394_2.png) [@RaphaelC](https://fenicsproject.discourse.group/u/RaphaelC)\
**Post date:** [April 12, 2022, 9:03am UTC](https://fenicsproject.discourse.group/t/integral-of-inner-surface-wrongly-evaluate-to-zero-in-3d/8031/1 "2022-04-12T09:03:54Z")

</div>

Hello,  
I have a 3D mesh containing an inner surface, and I want to add some force on it. By inner surface, I mean a surface between two blocks.  
I know how to do that with external forces.  
Let me give a MWE: first generate some example with gmsh, and get the xdmf files with meshio.

```python
import gmsh

gmsh.initialize()
tag_frame = 0
gmsh.model.occ.addBox(-1, -1, -1, 10, 10, 3, tag=tag_frame)
gmsh.model.occ.addBox(0, 0, -1, 8, 8, 3, tag=tag_frame+1)
tag_silicone = 2
tag_homogeneous = 3
gmsh.model.occ.cut([(3, tag_frame)], [(3, tag_frame+1)],tag=tag_silicone)
gmsh.model.occ.addBox(0, 0, -1, 8, 8, 3, tag=tag_homogeneous)
gmsh.model.occ.fragment([(3,tag_silicone)],[(3,tag_homogeneous)])
gmsh.model.occ.synchronize()
inner_marker = 2
xmin_marker = 3
xmin_interfaces=[]
inner_interfaces=[]
for surface in gmsh.model.getEntities(dim=2):
  com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
  if (abs(com[0]+1)<1.e-6):
    xmin_interfaces.append(surface[1])
  elif (abs(com[2]-0.5)<1.e-6):
    if (abs(com[0])<1.e-6): inner_interfaces.append(surface[1])
    elif (abs(com[0]-8)<1.e-6): inner_interfaces.append(surface[1])
    elif (abs(com[1])<1.e-6): inner_interfaces.append(surface[1])
    elif (abs(com[1]-8)<1.e-6): inner_interfaces.append(surface[1])
print(f"Nb of left surfaces: {len(xmin_interfaces)}")
print(f"Nb of inner surfaces: {len(inner_interfaces)}")
gmsh.model.addPhysicalGroup(2, inner_interfaces, inner_marker)
gmsh.model.setPhysicalName(2, inner_marker, "neutral interface")
gmsh.model.addPhysicalGroup(2, xmin_interfaces, xmin_marker)
gmsh.model.setPhysicalName(2, xmin_marker, "xmin interface")
volumes = gmsh.model.getEntities(dim=3)
silicon_marker = 11
homogeneous_marker = 12
gmsh.model.addPhysicalGroup(3, [tag_silicone], silicon_marker)
gmsh.model.setPhysicalName(3, silicon_marker, "Silicone")
gmsh.model.addPhysicalGroup(3, [tag_homogeneous], homogeneous_marker)
gmsh.model.setPhysicalName(3, homogeneous_marker, "HomogeneousMaterial")
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.15)
gmsh.model.mesh.generate(3)
gmsh.write("Box.msh")
gmsh.finalize()

import meshio
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh
msh = meshio.read("Box.msh")
tetra_mesh = create_mesh(msh, "tetra")
triangle_mesh = create_mesh(msh, "triangle")
meshio.write("Box.xdmf", tetra_mesh)
meshio.write("BoxT.xdmf", triangle_mesh)

```

Now if I open BoxT.xdmf with paraview, I can seem my interfaces, everything seems good.  
So let’s load this with fenicsX:

```python
import numpy as np
import ufl
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh
from dolfinx.io import XDMFFile
file_name = "Box.xdmf"
file_nameT = "BoxT.xdmf"

with XDMFFile(MPI.COMM_WORLD, file_name, "r") as xdmf:
  domain = xdmf.read_mesh(name="Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, file_nameT, "r") as xdmf:
  ft = xdmf.read_meshtags(domain, name = "Grid")
V = fem.VectorFunctionSpace(domain, ("CG", 2))

inner_marker = 2
xmin_marker = 3

attach_facets = ft.indices[ft.values==xmin_marker]
attach_dofs = fem.locate_dofs_topological(V, 3, attach_facets)
pressure_facets = ft.indices[ft.values==inner_marker]

marked_facets = np.hstack([attach_facets, pressure_facets])
marked_values = np.hstack([
    np.full(len(attach_facets), 1, dtype = np.int32), 
    np.full(len(pressure_facets), 2, dtype = np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, 2, marked_facets[sorted_facets], marked_values[sorted_facets])

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
SurfaceForm_inner = fem.form(fem.Constant(domain, PETSc.ScalarType(1.))*ds(2))
surface_inner = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner), op=MPI.SUM)
SurfaceForm_attach = fem.form(fem.Constant(domain, PETSc.ScalarType(1.))*ds(1))
surface_attach = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_attach), op=MPI.SUM)
print(f"Checking attach surface={surface_attach}")
print(f"Checking pressure surface={surface_inner}")

```

The output is:

```
Checking attach surface=30.00000000000009
Checking pressure surface=0.0

```

So something is wrong, the outer surface evaluates to correct value, but inner surface gets zero . I checked that pressure\_facets was not empty. Just taking one of the inner side instead also leads to zero value.  
One other strange thing is that changing mesh characteristic lengths

```auto
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.15)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.5)

```

leads to PETSC ERROR: Caught signal number 11 SEGV.

Would someone have a hint how to correct this, so that I can work with inner surfaces?

Thanks

---

<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 12, 2022, 10:14am UTC](https://fenicsproject.discourse.group/t/integral-of-inner-surface-wrongly-evaluate-to-zero-in-3d/8031/2 "2022-04-12T10:14:06Z")

</div>

> [@RaphaelC](#):
>
> `ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)`

The integration measure for facets that is shared between two cells is `"dS"` not `"ds"`, which is the integration measure for external facets, i.e. those only connected to a single cell.