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

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) + 2mu*epsilon(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