# Integrating a tensor along a line

#1

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

def sigma(u):
return lambda_nabla_div(u)Identity(d) + 2muepsilon(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!

#2

Strange that the experts do not know how to solve this fairly basic problem!