Integrating function over internal marked boundary


Hello everyone,
I am solving simple linear elasticity problem, I want to calculate traction over the internal marked boundary, I am always getting a value zero, but I am successful in calculating over the external boundary. Also, I am unable to calculate the length of any internal marked line.

from fenics import *
from dolfin import *

Scaled variables

L = 1; W = 1
mu = 7.6923e10
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 11.53846e10
lambda_ = beta
g = gamma

Create mesh and define function space

mesh = RectangleMesh(Point(0,0), Point(1, 1), 50, 50, “left”)#UnitSquareMesh(25, 25)
V = VectorFunctionSpace(mesh, ‘P’, 1)

Define boundary condition

class bottom(SubDomain):
def inside(self,x, on_boundary):
tol = 1e-14
return (near(x[1],0.00)) and on_boundary

class top(SubDomain):
def inside(self,x, on_boundary):
tol = 1e-10
return near(x[1],1.00) and on_boundary

class Marked(SubDomain):
def inside(self,x, on_boundary):
tol = 1e-10
return near(x[0]+x[1],0.5,0.002)

boundaries = MeshFunction(“size_t”,mesh,mesh.topology().dim()-1)
bottom= bottom()
top = top()


bc_bottom = DirichletBC(V.sub(1),Constant((0)),bottom)
bc = [bc_bottom]
ds = Measure(‘ds’, domain=mesh, subdomain_data=boundaries)

Define strain and stress

def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))

def sigma(u):
return lambda_nabla_div(u)Identity(d) + 2muepsilon(u)

ex = Constant((1.0, 0.0))
ey = Constant((0.0, 1.0))

n = FacetNormal(mesh)

Define variational problem

u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0,0))
T = Constant((0, 100000))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds(2)

Compute solution

u = Function(V)
solve(a == L, u, bc)

L1 = assemble(1*ds(2)) #Length of Top Marked Boundary
print L1 #Calculating correct

L2 = assemble(1*ds(3)) ##Length of Interior Marked Boundary
print L2 #Alwayz giving Zero for any interior length

#Want to calculate Traction
Tx = assemble(dot(dot(sigma(u), n),ex)*ds(3))
print Tx #It always gives zero

Ty = assemble(dot(dot(sigma(u), n),ey)*ds(3))
print Ty #It always gives zero

Thanks in advance for your time


The facets you mark with 3 are interior facets of the mesh hence the measure to be used for surface itengration over them is dS, cf. ds for exterior facet integrals. The following post could be useful.


Your suggestion worked. Thank you MiroK!!