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

I have edited. But it is still giving en error. Please help!

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!!