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