# Facets orientation in inner surface with fenicsX

**URL:** <https://fenicsproject.discourse.group/t/facets-orientation-in-inner-surface-with-fenicsx/8034>\
**Category:** Uncategorized\
**Created:** [April 12, 2022, 3:39pm UTC](https://fenicsproject.discourse.group/t/facets-orientation-in-inner-surface-with-fenicsx/8034 "2022-04-12T15:39:47Z")\
**Posts on this page:** 4\
**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, 3:39pm UTC](https://fenicsproject.discourse.group/t/facets-orientation-in-inner-surface-with-fenicsx/8034/1 "2022-04-12T15:39:47Z")

</div>

Hello,  
I’ve been trying to add some force on an inner surface. However, the facets seems to have incoherent orientation.  
According to [Integrating over an interior surface - #4 by MiroK](https://fenicsproject.discourse.group/t/integrating-over-an-interior-surface/247/4) using different mesh domains on both sides of the surface should be sufficient to get a coherent orientation. Let me introduce a MWE. The mesh is generated with gmsh and consist of two blocks of same height, one inserted in the other. I have one tag for each: tag\_silicone and tag\_homogeneous. I will also tag each interface surface separately. Here is the code:

```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_1 = 10
inner_marker_2 = 11
inner_marker_3 = 12
inner_marker_4 = 13
inner_interfaces_1=[]
inner_interfaces_2=[]
inner_interfaces_3=[]
inner_interfaces_4=[]
for surface in gmsh.model.getEntities(dim=2):
  com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
  if (abs(com[2]-0.5)<1.e-6):
    if (abs(com[0])<1.e-6): inner_interfaces_1.append(surface[1])
    elif (abs(com[0]-8)<1.e-6): inner_interfaces_2.append(surface[1])
    elif (abs(com[1])<1.e-6): inner_interfaces_3.append(surface[1])
    elif (abs(com[1]-8)<1.e-6): inner_interfaces_4.append(surface[1])
gmsh.model.addPhysicalGroup(2, inner_interfaces_1, inner_marker_1)
gmsh.model.setPhysicalName(2, inner_marker_1, "neutral interface")
gmsh.model.addPhysicalGroup(2, inner_interfaces_2, inner_marker_2)
gmsh.model.setPhysicalName(2, inner_marker_2, "neutral interface")
gmsh.model.addPhysicalGroup(2, inner_interfaces_3, inner_marker_3)
gmsh.model.setPhysicalName(2, inner_marker_3, "neutral interface")
gmsh.model.addPhysicalGroup(2, inner_interfaces_4, inner_marker_4)
gmsh.model.setPhysicalName(2, inner_marker_4, "neutral 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 I will compute the surface of each interface surface of 4 different manners.

- Integrate 1dS over the corresponding domain
- Integrate n\cdot n\_0 dS where n\_0 is the reference normal, and n the computed normal
- Integrate n\cdot n\_0 dS + 0 dx i.e. try the trick to add dx with subdomains to get coherent orientation of facets
- Integrate (n\cdot n\_0)^2 dS

Note: In practice, I don’t need the surface, but I have to use the normal.  
The resulting script is

```python
from weakref import ref
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")
  ct = xdmf.read_meshtags(domain, 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))

silicon_marker = 11
homogeneous_marker = 12
silicon_tets = ct.indices[ct.values==silicon_marker]
silicon_dofs = fem.locate_dofs_topological(V, 3, silicon_tets)
homogeneous_tets = ct.indices[ct.values==homogeneous_marker]
homogeneous_dofs = fem.locate_dofs_topological(V, 3, homogeneous_tets)

marked_tets = np.hstack([silicon_tets, homogeneous_tets])
marked_values = np.hstack([np.full(len(silicon_tets), 1, dtype = np.int32), np.full(len(homogeneous_tets), 2, dtype = np.int32)])
sorted_tets = np.argsort(marked_tets)
tet_tag = mesh.meshtags(domain, 3, marked_tets[sorted_tets], marked_values[sorted_tets])

inner_marker_1 = 10
inner_marker_2 = 11
inner_marker_3 = 12
inner_marker_4 = 13

pressure_facets_1 = ft.indices[ft.values==inner_marker_1]
pressure_facets_2 = ft.indices[ft.values==inner_marker_2]
pressure_facets_3 = ft.indices[ft.values==inner_marker_3]
pressure_facets_4 = ft.indices[ft.values==inner_marker_4]

marked_facets = np.hstack([pressure_facets_1, pressure_facets_2, pressure_facets_3, pressure_facets_4])
marked_values = np.hstack([
    np.full(len(pressure_facets_1), 0, dtype = np.int32),
    np.full(len(pressure_facets_2), 1, dtype = np.int32),
    np.full(len(pressure_facets_3), 2, dtype = np.int32),
    np.full(len(pressure_facets_4), 3, 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)
dx = ufl.Measure("dx", domain=domain, metadata=metadata, subdomain_data=tet_tag)

u = fem.Function(V)
I = ufl.variable(ufl.Identity(3))
F = ufl.variable(I+ufl.grad(u))
J = ufl.variable(ufl.det(F))
N = ufl.FacetNormal(domain)
Finv = ufl.inv(F)
n = J*Finv.T*N

vec_x = fem.Constant(domain,PETSc.ScalarType((1.0,0,0)))
vec_y = fem.Constant(domain,PETSc.ScalarType((0.0,1.0,0)))
vecs=[-vec_x,vec_x,-vec_y,vec_y]
reference = 8.0*3.0
for i in range(0,4):
    SurfaceForm_inner = fem.form(fem.Constant(domain, PETSc.ScalarType(1.))*dS(i))
    surface_inner = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner), op=MPI.SUM)
    SurfaceForm_inner_2 = fem.form(ufl.inner(n('-'),vecs[i])*dS(i))
    surface_inner_2 = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner_2), op=MPI.SUM)
    SurfaceForm_inner_3 = fem.form(ufl.inner(n('-'),vecs[i])*dS(i)+ fem.Constant(domain, PETSc.ScalarType(0.0))*dx)
    surface_inner_3 = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner_3), op=MPI.SUM)
    SurfaceForm_inner_4 = fem.form(ufl.inner(n('-'),vecs[i])**2*dS(i))
    surface_inner_4 = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner_4), op=MPI.SUM)

    print(f"Checking pressure on surface {i}: reference={reference}, {surface_inner} {surface_inner_2} {surface_inner_3} {surface_inner_4}")

```

Each output line should contains 4 output values, first and last should be coherent as they do not depends on the facets orientation. A correct last value should indicate the n is normal to the facet. If I understood the trick to get coherent orientation, third value should be correct.  
However my output is

```auto
Checking pressure on surface 0: 
   reference=24.0, 23.999999999999897 23.80972755772505 23.80972755772505 23.999999999999897
Checking pressure on surface 1: 
   reference=24.0, 23.99999999999983 -18.501788312810014 -18.501788312810014 23.99999999999983
Checking pressure on surface 2: 
   reference=24.0, 23.99999999999986 23.736066937512884 23.736066937512884 23.99999999999986
Checking pressure on surface 3: 
   reference=24.0, 24.0 -18.570081136289787 -18.570081136289787 24.0

```

Any idea what am I doing wrong ?  
Best

---

<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 13, 2022, 11:39am UTC](https://fenicsproject.discourse.group/t/facets-orientation-in-inner-surface-with-fenicsx/8034/2 "2022-04-13T11:39:16Z")

</div>

After further inspection of the source code, it seems like there is a bug that makes it impossible to order the interior facet integrals (for now).  
Im working on a fix.

---

<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 13, 2022, 12:15pm UTC](https://fenicsproject.discourse.group/t/facets-orientation-in-inner-surface-with-fenicsx/8034/3 "2022-04-13T12:15:14Z")

</div>

I found some way to fix my problem, but it does not generalize very well. As I know the real normal, I can then use the ufl conditional operator:

```python
vec_x = fem.Constant(domain,PETSc.ScalarType((1.0,0,0)))
vec_y = fem.Constant(domain,PETSc.ScalarType((0.0,1.0,0)))
vecs=[-vec_x,vec_x,-vec_y,vec_y]
for i in range(0,4):
    n_corr=ufl.conditional(ufl.le(ufl.inner(N('+'),vecs[i]),0.0),n('-'),n('+'))
    SurfaceForm_inner = fem.form(ufl.inner(n_corr,vecs[i])*dS(i))#+ fem.Constant(domain, PETSc.ScalarType(0.0))*dx)
    surface_inner = domain.comm.allreduce(fem.assemble_scalar(SurfaceForm_inner), op=MPI.SUM)

    print(f"Checking pressure on surface {i}: reference={reference}, {surface_inner} ")

```

However it works only as my normal vector is very specific (constant on each side).

---

<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 13, 2022, 2:27pm UTC](https://fenicsproject.discourse.group/t/facets-orientation-in-inner-surface-with-fenicsx/8034/4 "2022-04-13T14:27:59Z")

</div>

I have pushed a fix: [Fix consistent restriction of one-sided interior facet integrals by jorgensd · Pull Request #2127 · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/pull/2127)  
Hopefully it will get merged soon, and the example in the issue should illustrate how to use the restrictions.