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