# Marking internal facets at sharp angles

I have an issue integrating over an internal boundary that contains a 90° angle. I found a fix but wonder if there is not a more direct way of doing it.

``````import ufl
from mpi4py import MPI
import dolfinx
import numpy as np
import pyvista

mesh=dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,[(-5,0),(5,20)],(10,20))
ylower=10
yupper=15
x1=1
x2=4
def left(x):
lhoriz=np.isclose(x[1],ylower) & (x[0]>=-x2) & (x[0]<=-x1)
lvert=np.isclose(x[0],-x1) & (x[1]>=ylower) & (x[1]<=yupper)
return(lhoriz | lvert)
def right(x):
rhoriz=np.isclose(x[1],ylower) & (x[0]>=x1) & (x[0]<=x2)
rvert=np.isclose(x[0],x1) & (x[1]>=ylower) & (x[1]<=yupper)
return(rhoriz | rvert)

V = dolfinx.fem.functionspace(mesh, ("CG", 1))
left_facets=dolfinx.mesh.locate_entities(mesh,1,left)
right_facets=dolfinx.mesh.locate_entities(mesh,1,right)
fdim=1
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
# Fixing erroneously marked facets
cc=mesh.topology.connectivity(1,0)
xequal=mesh.geometry.x[cc.array.reshape(-1,2)[marked_facets]][:,0,0]==mesh.geometry.x[cc.array.reshape(-1,2)[marked_facets]][:,1,0]
yequal=mesh.geometry.x[cc.array.reshape(-1,2)[marked_facets]][:,0,1]==mesh.geometry.x[cc.array.reshape(-1,2)[marked_facets]][:,1,1]
marked_facets=marked_facets[xequal|yequal]
marked_values=marked_values[xequal|yequal]
# Fix end
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
dS=ufl.Measure("dS",mesh,subdomain_data=facet_tag)
lleft=dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dS(1)))
lright=dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dS(2)))
print(lleft,lright)
``````

Without the fix my left integral is sqrt(2) too large due to a diagonal facet being marked in the corner.

Is there a more Fenicx-onian way of doing this?

Thanks

You may want to consider using gmsh for generating the mesh, alongside its markers. There are several examples posted in the forum, as well as one demo.

1 Like

Thanks!
That is what I usually do in more complicated cases (or I create the mesh in Salome), but for this simple test case I thought it would be easier to do it with the tools dolfinx provide.

Seems to me that the artifact you see is due to using `|`. If I split your conditions in 4 separate functions I get the right markers without any workaround.

``````import ufl
from mpi4py import MPI
import dolfinx
import numpy as np
import pyvista

mesh=dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,[(-5,0),(5,20)],(10,20))
ylower=10
yupper=15
x1=1
x2=4

V = dolfinx.fem.functionspace(mesh, ("CG", 1))
left_facets_1=dolfinx.mesh.locate_entities(mesh,1,lambda x: np.isclose(x[1],ylower) & (x[0]>=-x2) & (x[0]<=-x1))
left_facets_2=dolfinx.mesh.locate_entities(mesh,1,lambda x: np.isclose(x[0],-x1) & (x[1]>=ylower) & (x[1]<=yupper))
right_facets_1=dolfinx.mesh.locate_entities(mesh,1,lambda x: np.isclose(x[1],ylower) & (x[0]>=x1) & (x[0]<=x2))
right_facets_2=dolfinx.mesh.locate_entities(mesh,1,lambda x: np.isclose(x[0],x1) & (x[1]>=ylower) & (x[1]<=yupper))
fdim=1
marked_facets = np.hstack([left_facets_1, left_facets_2, right_facets_1, right_facets_2])
marked_values = np.hstack([np.full_like(left_facets_1, 1), np.full_like(left_facets_2, 1), np.full_like(right_facets_1, 2), np.full_like(right_facets_2, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
``````