Specifying quantities on interior surfaces with dolfinx

This question has been clearly answered for old dolfin, but I cannot find anything yet on specifying discontinuous quantities on interior surfaces in dolfinx.

I define an interior surface dS and then wish to integrate over this surface using discontinuous quantities (such as FacetNormal, which has opposite sign on either side of the facet). How do I specify which side of the facet I am looking at?

Here is a working example using the ‘old dolfin’ method, where I define my subdomains and the surface between them, then integrate a vector dotted with the normal, whose resulting sign depends on the normal direction:

``````from mpi4py import MPI
from dolfinx import mesh, fem, io
import ufl
from petsc4py import PETSc
import numpy as np

domain = mesh.create_unit_square(MPI.COMM_WORLD, 4,4)

def Omega_0(x):
return x[0] <= 0.5

def Omega_1(x):
return x[0] >= 0.5

#mark left as 1, right as 2
cells_0 = mesh.locate_entities(domain, domain.topology.dim, Omega_0)
cells_0_mark = np.zeros_like(cells_0) + 1
cells_1 = mesh.locate_entities(domain, domain.topology.dim, Omega_1)
cells_1_mark = np.zeros_like(cells_1) + 2
cells= np.concatenate((cells_0,cells_1))
cells_mark = np.concatenate((cells_0_mark,cells_1_mark))
ct = mesh.meshtags(domain,domain.topology.dim, np.array(cells).astype(np.int32),np.array(cells_mark).astype(np.int32))

def d_Omega(x):
return abs(x[0] - 0.5) < 1e-7

d_facets = mesh.locate_entities(domain, domain.topology.dim -1 , d_Omega)
d_facets_mark = np.zeros_like(d_facets) + 1
ft = mesh.meshtags(domain,domain.topology.dim-1, np.array(d_facets).astype(np.int32),np.array(d_facets_mark).astype(np.int32))
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)

dS = ufl.Measure('dS',domain=domain, subdomain_data=ft)

v_cg2 = ufl.VectorElement('CG', domain.ufl_cell(),2)
V2 = fem.FunctionSpace(domain,v_cg2)

uD2 = fem.Function(V2)
def u_exact(x):
return np.array([x[0],x[0]-x[0]], dtype=PETSc.ScalarType)
uD2.interpolate(u_exact)

uD0 = fem.Function(V2)
def u_zero(x):
return np.array([x[0]-x[0],x[0]-x[0]], dtype=PETSc.ScalarType)
uD0.interpolate(u_zero)
n = ufl.FacetNormal(domain)

# here the larger mesh tags are on the right, so the normals should point left, and the
# value should be -0.5, which we get
domain_constant = (ufl.inner(uD0,uD0))*ufl.dx(domain=domain, subdomain_data=ct)
test= fem.form( ufl.inner(uD2, n('+'))*dS(1) + domain_constant)
test_val = fem.assemble_scalar(test)
print(test_val) #-0.5

#mark left as 2, right as 1
cells_0_mark_reverse = np.zeros_like(cells_0) + 2
cells_1_mark_reverse = np.zeros_like(cells_1) + 1
cells_mark_reverse = np.concatenate((cells_0_mark_reverse,cells_1_mark_reverse))
ct_reverse = mesh.meshtags(domain,domain.topology.dim, np.array(cells).astype(np.int32),np.array(cells_mark_reverse).astype(np.int32))

# here the larger mesh tags are on the left, so the normals should point right, and the
# value should be +0.5, however we still get -0.5
reverse_domain_constant = (ufl.inner(uD0,uD0))*ufl.dx(domain=domain, subdomain_data=ct_reverse)
test = fem.form( ufl.inner(uD2, n('+'))*dS(1) + reverse_domain_constant)
test_val = fem.assemble_scalar(test)
print(test_val) #-0.5
``````

Thanks,
Jennifer

There is a patch for this coming up (by @jpdean in Let integration entities be specified manually by jpdean · Pull Request #2269 · FEniCS/dolfinx · GitHub).

For now I would suggest the following:

``````from mpi4py import MPI
from dolfinx import mesh, fem
import ufl
from petsc4py import PETSc
import numpy as np

domain = mesh.create_unit_square(MPI.COMM_WORLD, 4, 4)

def Omega_0(x):
return x[0] <= 0.5

def Omega_1(x):
return x[0] >= 0.5

def d_Omega(x):
return abs(x[0] - 0.5) < 1e-7

d_facets = mesh.locate_entities(domain, domain.topology.dim - 1, d_Omega)
d_facets_mark = np.zeros_like(d_facets) + 1
ft = mesh.meshtags(domain, domain.topology.dim-1, np.array(
d_facets).astype(np.int32), np.array(d_facets_mark).astype(np.int32))
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)

dS = ufl.Measure('dS', domain=domain, subdomain_data=ft)

v_cg2 = ufl.VectorElement('CG', domain.ufl_cell(), 2)
V2 = fem.FunctionSpace(domain, v_cg2)

uD2 = fem.Function(V2)

def u_exact(x):
return np.array([x[0], x[0]-x[0]], dtype=PETSc.ScalarType)

uD2.interpolate(u_exact)

n = ufl.FacetNormal(domain)

# mark left as 1, right as 2
cells_0 = mesh.locate_entities(domain, domain.topology.dim, Omega_0)
cells_1 = mesh.locate_entities(domain, domain.topology.dim, Omega_1)

V_DG = fem.FunctionSpace(domain, ("DG", 0))
weight = fem.Function(V_DG)
weight.x.array[cells_0] = 1
weight.x.array[cells_1] = 0
weight.x.scatter_forward()

# Integrate from right side, n = (1, 0)
test = fem.form(2*ufl.avg(weight*ufl.inner(uD2, n))*dS(1))
test_val = fem.assemble_scalar(test)
assert np.isclose(domain.comm.allreduce(test_val, op=MPI.SUM), 0.5)

# Integrate from left side, n = (-1, 0)
weight.x.array[cells_0] = 0
weight.x.array[cells_1] = 1
test_val = fem.assemble_scalar(test)
assert np.isclose(domain.comm.allreduce(test_val, op=MPI.SUM), -0.5)
``````
1 Like

Thank you, this works for my applications! I look forward to the new patch and the functionality it will add.