Mixed Dimensional Integration Measure "dS"

1

Hi all,

I am using the mixed-dimensional branch to simulate a subsurface fault/fracture in a continuum rock (@cdaversin). The “rock” is a 2D surface and the “fault” is a 1D line that is completely embedded in the surface. I am solving the equation of linear elasticity with contact mechanics imposed with Lagrange multipliers on the internal line. Here is the equation:

∫Ω (𝜎:∇𝑣) dΩ + ∫Γ (𝜆−𝑝[[𝑢]])⋅[[𝑣]] dΓ = 0 
∫Γ (−𝑝[[𝑢]])𝜇 dΓ = 0

Where lambda is the Lagrange multiplier, mu is the virtual Lagrange multiplier, u is the displacement, v is the virtual displacement, p is the penalty for contact, and [] represents the jump operator.

My question is this: Does the mixed dimensional branch allow for jump operations in the displacement field, and if so do I use the “dx” measure of the Lagrange multiplier or the “dS” measure of the displacement field? Here is a MWE:

from dolfin import *
import matplotlib.pyplot as plt

class Roller_X(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 0) \
        or near(x[0], 100)

class Roller_Y(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 0) \
        or near(x[1], 100)

class Fracture(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 50) and x[0] > 25-DOLFIN_EPS and x[0] < 75+DOLFIN_EPS

# Mesh generation
n_x = 30
n_y = 30
mesh = RectangleMesh(Point(0,0,0),Point(100,100,0),n_x,n_y, "left/right")

# Interface markers
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Fracture().mark(facets,1)

# Meshes
submesh = MeshView.create(facets, 1)

# Function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
LM = FunctionSpace(submesh, "CG", 1)
W = MixedFunctionSpace(V, LM)

# Elasticity parameters
E, nu = 5e9, 0.22
lmbda = Constant(E*nu/(1+nu)/(1-2*nu))
mu = Constant(E/2/(1+nu))
pN = Constant(1e11) # Contact penalty [Pa/m]
n = FacetNormal(mesh)

# Neumann condition
g = Constant((0.0, -1.0e6))

# Functions
u, ulm = TrialFunctions(W)
v, vlm = TestFunctions(W)

bc_x = DirichletBC(W.sub_space(0).sub(0), 0.0, Roller_X())
bc_y = DirichletBC(W.sub_space(0).sub(1), 0.0, Roller_Y())

#---------------------------------------------------------------
# QUESTION HERE ------------------------------------------------
#---------------------------------------------------------------
dx = Measure("dx", domain=W.sub_space(0).mesh())
dS  = Measure("dS", domain=W.sub_space(0).mesh(), subdomain_data=facets)
dxf = Measure("dx", domain=W.sub_space(1).mesh())
#---------------------------------------------------------------

# Stress tensor
def sigma(u):
    return 2.0*mu*epsilon(u) + lmbda*tr(sym(grad(u)))*Identity(2)

# Strain tensor
def epsilon(u):
    return sym(grad(u))

nn = Constant((0,1))

# Variational problem
a = inner(sigma(u), epsilon(v))*dx + (ulm('+')-(pN*jump(u,nn)))*jump(v,nn)*dS(1) -  + pN*jump(u,nn)*vlm*dxf
l = inner(g, v)*dx + Constant(0.0)*vlm*dxf

# Solve
u_sol = Function(W)
solve(a==l, u_sol, [bc_x, bc_y], solver_parameters={"linear_solver":"direct"})
(uu, lam) = u_sol.split()

stress = sigma(uu)
p = plot(uu.sub(1), mode='color')
plt.colorbar(p)
#plot(stress[1,1])
plt.title(r"$u_y [m]$",fontsize=16)
plt.show()