Marking internal facets at sharp angles

dolfinx
1

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

2

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
3

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.

4

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])