# Why is this simple traction problem giving an error

#1

I cannot understand why such a simple problem is giving an error. My code is given below:

``````from fenics import *
from dolfin import *
import ufl
import dolfin.cpp as cpp
import numpy
# 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(100, 100), 100, 100, "left")
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],100.00) and on_boundary

class Marked(SubDomain):
def inside(self,x, on_boundary):
tol = 1e-10
return near(x[0]+x[1],100,0.001) and x[0]>25 and x[1]>25

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]

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

dummy1 = as_tensor([\
[1.0, 0.0],\
[0.0, 0.0]\
])
dummy2 = as_tensor([\
[0.0, 1.0],\
[0.0, 0.0]\
])
dummy3 = as_tensor([\
[0.0, 0.0],\
[1.0, 0.0]\
])
dummy1 = as_tensor([\
[0.0, 0.0],\
[0.0, 1.0]\
])

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
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) + 2*mu*epsilon(u)

def sigel1(u):
return inner(sigma(u),dummy1)
def sigel2(u):
return inner(sigma(u),dummy2)
def sigel3(u):
return inner(sigma(u),dummy3)
def sigel4(u):
return inner(sigma(u),dummy4)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0,0))
T = Constant((0, 1000))
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)

Tx = assemble(0.7071*(sigel1(u)+sigel2(u))*dS(3))
print Tx
Ty = assemble(0.7071*(sigel3(u)+sigel4(u))*dS(3))
print Ty
``````

Please, please look into it. Thank you!`Preformatted text`

#2

Please format your code such that we can directly copy, paste and run it. In markdown, * and # etc. get translated into formatting of your text which is why running your code verbatim won’t work.

#3

Done! Could you please solve my problem.

#4

No. The function space `V` is not defined in your code.

#5

#6

Hello benny,

As far as I saw, you have two different mistake in your code.

The first one is, you have two seperate `dummy1` in your code. You have to change the latest one to `'dummy4'`

The second one is, you have to specify on which side of the facet you want to integrate in your dS as following:

``````Tx = assemble(0.7071*(sigel1(u('-'))+sigel2(u('-')))*dS(3))
print ('Tx',Tx)
Ty = assemble(0.7071*(sigel3(u('-'))+sigel4(u('-')))*dS(3))
print ('Ty',Ty)
``````

So there should be (’+’) or (’-’) side of the facet. Otherwise the result will be 0 probably.

#7

It worked and I am getting the correct result. Thank you so much CMA!!