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