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

97

views

0

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

Thank you in advance for your invaluable help.

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:- 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?
- 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

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?`assemble_local(kappa('+')*dS(1),cell)`

- 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:

How one can formulate the following bilinear form over the internal boundary`#===== 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)`

and be sure that the values are taken from proper subdomains (sides of the boundary)?`n = FacetNormal(mesh) a = ( dp('+')*inner(u('+'),n('+')) + p('+')*inner(du('+'),n('+')) )*dS(1)`

Thank you in advance for your invaluable help.

```
# 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

Please login to add an answer/comment or follow this question.