### How to correctly integrate (assemble) over internal facets?

97
views
0
3 months ago by
Dear FEniCS community,

I am trying to compute (assemble) a simple integral over internal facets in order to deepen the topic available here (the meaning of '+' and '-' on interior facets), which is still not clear to me. I generated a simple 2D mesh with a straight-line vertical internal boundary marked by MeshFunction. It consists of four facets of equal size=1. As expected, the height of the mesh obtained from
assemble(1.*dS(1))​
is 4. However, if I define an expression called kappa which takes different values on both sides of the internal boundary, namely 1 to the left of it and .2 on the right side, and want to do similar calculations, the result is wrong regardless of the choice of '+' or '-' in
assemble(kappa('+')*dS(1))​​
What I did in order to circumvent this problem is a manual selection of proper kappa('+') and kappa('-') values depending on the side of the internal boundary they are taken from. For the simple MWE attached below the output is
assemble without loop:  1.6
assemble with loop:  4.0
Now I have three questions related to the topic:
1. Is it possible to compute such integrals correctly without looping over all facets in the mesh? In other words, how to make kappa('+') and kappa('-') purposeful and predictable? In the presented application kappa('+') took .2 three times and 1. once. Is it possible to somehow redefine the mesh to make it consistent?
2. Why kappa('+') and kappa('-') values are defined only for a cell on just one side of the internal boundary, not both? In particular, why do I get very big numbers when assembling a local integral
assemble_local(kappa('+')*dS(1),cell)​
using a cell on one side of the internal boundary, and reasonable, proper values when doing so with respect to a cell on the other side of the boundary?
3. How can I mirror such reasoning to the integration of test and trial functions over internal facets? Suppose one has the test and trial functions declared as follows on the considered mesh:
#===== FUNCTION SPACES =====
Uh = VectorElement('P', mesh.ufl_cell(), 2)
Ph = FiniteElement('P', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([Uh, Ph]))
(u, p) = TrialFunctions(W)
(du, dp) = TestFunctions(W)​​
How one can formulate the following bilinear form over the internal boundary
n = FacetNormal(mesh)
a = ( dp('+')*inner(u('+'),n('+')) + p('+')*inner(du('+'),n('+')) )*dS(1)​
and be sure that the values are taken from proper subdomains (sides of the boundary)?

# run using dolfin v2017.2.0
from dolfin import *
from mshr import *

def build_mesh(l1, l2, h, nl1, nl2, nh):
domain_1 = Rectangle(Point(.0, .0), Point(l1, h)) # https://fenicsproject.org/qa/11220/how-to-merge-concatenate-add-meshed-in-one-mesh/
domain_2 = Rectangle(Point(l1, h), Point(l1+l2, .0))
domain_tot = domain_1 + domain_2

domain_tot.set_subdomain(1, domain_1) # https://fenicsproject.org/qa/9512/consisting-geometries-boundary-between-geometries-contained/
domain_tot.set_subdomain(2, domain_2)
m = generate_mesh(domain_tot, 25) # FIXME: should use nl1, nl2, nh to control the mesh resolution

return m

#===== BORDERS =====
# Geometric parameters
H = 4.   # height
L2 = 24. # second thickness
L1 = L2  # first thickness

#===== MESH =====
# Create mesh
mesh = build_mesh(L1, L2, H, 10, 10, 8)

#===== SUBDOMAINS =====
# Define subdomain markers
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), mesh.domains())

#===== BOUNDARIES =====
# Define boundaries
tol = 1.E-14
def Interface(x):
return near(x[0], L1, tol)

# Create mesh function over cell facets
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_subdomains.set_all(0) # 0 for interior cell boundaries
interface_boundary = AutoSubDomain(Interface)
interface_boundary.mark(boundary_subdomains, 1) # 1 for interface boundary

# Define integration measure for boundary integrals
dS = Measure('dS', domain=mesh, subdomain_data=boundary_subdomains) # measure for internal boundaries

#===== MESH DATA =====
kappa = Expression('x[0] <= L1 + tol ? 1. : .2', degree=0, L1=L1, tol=tol)
dim = mesh.topology().dim()
mesh.init(dim-1,dim) # Build connectivity between facets and cells
print "assemble without loop: ",assemble(kappa('+')*dS(1))
assemble_value = 0.
for face in facets(mesh): # https://fenicsproject.org/qa/8830/compute-element-wise-integral-on-subdomain/
if boundary_subdomains[face]==1: # if a facet belongs to the interface_boundary subdomain
#		print face.index()
if subdomains[int(face.entities(dim)[0].item())]==1: # if a cell sharing the facet belongs to the first subdomain, take this cell
cell = Cell(mesh, face.entities(dim)[0])
assemble_value_plus = assemble_local(kappa('+')*dS(1),cell).item()
assemble_value_minus = assemble_local(kappa('-')*dS(1),cell).item()
if assemble_value_plus < 1.E100 and assemble_value_minus < 1.E100:
pass # do nothing (the values are reasonable)
elif assemble_value_plus > 1.E100 and assemble_value_minus > 1.E100:
cell = Cell(mesh, face.entities(dim)[1])
assemble_value_plus = assemble_local(kappa('-')*dS(1),cell).item() # now the cell lies on the opposite side of the interface and therefore has '+' and '-' of opposite meaning ('+' sign means the value at the cell on the current side of the interface)
assemble_value_minus = assemble_local(kappa('+')*dS(1),cell).item()
else:
print "WARNING: Alternating values..."
assemble_value += assemble_value_plus
else: # otherwise, take another cell
cell = Cell(mesh, face.entities(dim)[1])
assemble_value_plus = assemble_local(kappa('+')*dS(1),cell).item()
assemble_value_minus = assemble_local(kappa('-')*dS(1),cell).item()
if assemble_value_plus < 1.E100 and assemble_value_minus < 1.E100:
pass # do nothing (the values are reasonable)
elif assemble_value_plus > 1.E100 and assemble_value_minus > 1.E100:
cell = Cell(mesh, face.entities(dim)[0])
assemble_value_plus = assemble_local(kappa('-')*dS(1),cell).item() # now the cell lies on the opposite side of the interface and therefore has '+' and '-' of opposite meaning ('+' sign means the value at the cell on the current side of the interface)
assemble_value_minus = assemble_local(kappa('+')*dS(1),cell).item()
else:
print "WARNING: Alternating values..."
assemble_value += assemble_value_plus
print "assemble with loop: ",assemble_value
Community: FEniCS Project