I would like to take consecutive dot product of a tensor with two constant vectors and integrate the resulting scalar along a line. For example, I want to integrate the scalar that results from the operation: (dot(dot(sigma(u), nx),ex) where ex and nx are two different constant vectors and sigma(u) is the stress tensor.

Here’s my attempt at it which gives as error in fenics.

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(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]

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

ex = Constant((1.0, 0.0))

ey = Constant((0.0, 1.0))

nx = Constant((0.7071, 0.0))

ny = Constant((0.0, 0.7071))

# Define variational problem

u = TrialFunction(V)

d = u.geometric_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)

#Want to calculate Traction

Tx = assemble(dot(dot(sigma(u), nx),ex)*dS(3))

print Tx

Ty = assemble(dot(dot(sigma(u), ny),ey)*dS(3))

print Ty

One suggestion would be to use facet normal but that does not give the right answer in unstructured mesh. It gives right answer only for structured mesh. So please fix the above code for me. Thank you very much!!

The multiplication symbols are not visible in the definition of epsilon and sigma. Please put those symbols before running the code. Thank you!