# Integrating function over internal marked boundary

#1

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)
boundaries.set_all(0)
bottom= bottom()
top = top()
bottom.mark(boundaries,1)
top.mark(boundaries,2)
Marked=Marked()
Marked.mark(boundaries,3)

meshfile=File(‘mesh.pvd’)
meshfile<<boundaries

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

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

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.