extract stress component at a point

Question

hi,

i am new to fenics. I want to plot stress component at a point.
Ex: I want to plot sigma_xx at point (1,1) within the domain.

sigma_w = D_mat*eps(u)
y = sigma_w[0] (pointA)

I will have now (3, ) vector. how to extract sigma_xx at a particular point.
getting error:

Thanks in advance

Answer 1

Try something like this:

T = TensorFunctionSpace(mesh, "DG", 0)
sigma = project(sigma_w, T)
[s11, s12, s21, s22] = sigma.split(True)
print s11(PointA)

or:

T = TensorFunctionSpace(mesh, "DG", 0)
sigma = project(sigma_w, T)
print (sigma.sub(0))(pointA)

Comments

i have got this error:

Traceback (most recent call last):

File "/home/hirshikesh/Research/Fenics/Thermomechanical/PlateHole_TwoWayNR_force.py", line 228, in
sigma = project(sigma_w,x)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/projection.py", line 107, in project
L = ufl.inner(w, v)dx
File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 130, in inner
return Inner(a, b)
File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 156, in new
ufl_assert(ash == bsh, "Shape mismatch.")
File "/usr/lib/python2.7/dist-packages/ufl/assertions.py", line 37, in ufl_assert
if not condition: error(message)
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 151, in error
raise self._exception_type(self._format_raw(*message))
UFLException: Shape mismatch.

Actually sigma_w is the vector of size (3, )

---

If sigma is a vector of tractions yo have to do this:

T = VectorFunctionSpace(mesh, "DG", 0) # you can use "CG" if it's better
sigma = project(sigma_w, T)
[s0, s1, s2] = sigma.split(True)
print s0(PointA)

or:

T = VectorFunctionSpace(mesh, "DG", 0)
sigma = project(sigma_w, T)
print (sigma.sub(0))(pointA)

---

I have done that also I am getting shape mismatch again.

---

Is a good idea show us a minimal working (or non-working) example of your code to reproduce the error message.

Regards.

---

Here I am attaching my code:

from dolfin import *
import numpy as np
import mshr
import matplotlib.pyplot as plt

import numpy as np

def D_Matrix(E, nu, stressState):

if stressState=='PlaneStress':
    a = E/(1-nu**2)
    D = np.array([[a, a*nu, 0],
                  [a*nu, a, 0],
                  [0, 0, a*(1-nu)/2],
                  ])
else:
    a = E/((1+nu)*(1-2*nu))
    D = np.array([[a*(1-nu), a*nu, 0],
                  [a*nu, a*(1-nu), 0],
                  [0, 0, a*(0.5-nu)],
                  ])
return D

domain = mshr.Rectangle(Point(-0.5,-0.50), Point(0.5,0.5))-mshr.Circle(Point(0.0,0.0), 0.1)
mesh = mshr.generate_mesh(domain,30)

plot(mesh,interactive = True)

stressState = 'PlaneStrain'

E = 19.25e09
nu = 0.3

V = VectorFunctionSpace(mesh,'CG',2)
u = TrialFunction(V)
v = TestFunction(V)

lmbda, mu = Enu/((1.0 + nu )(1.0-2.0nu)) , E/(2(1+nu))

class Left(SubDomain):
def inside(self,x,on_boundary):
tol = 1e-10
return (near(x[0],-0.5)) and on_boundary
class Right(SubDomain):
def inside(self,x,on_boundary):
tol = 1e-10
return (near(x[0],0.5)) and on_boundary

class Center(SubDomain):
def inside(self,x,on_boundary):
tol = 1e-10
return abs(x[0] + 0.5) < tol and abs(x[1] - 0.0) < tol

left = Left()
right = Right()
center = Center()

ccx = DirichletBC(V.sub(0),Constant(0.0),left)
ccy = DirichletBC(V.sub(1),Constant(0.0),left)

dirichlet_u = [ccx , ccy]

boundaries = FacetFunction("size_t",mesh)
boundaries.set_all(0)
center.mark(boundaries, 1)
left.mark(boundaries, 2)
right.mark(boundaries, 3)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

gR = Constant((1.0, 0.0))

D_mat = D_Matrix(E,nu,stressState)
D_mat = as_matrix(D_mat)

def eps(u):
return as_vector([u[i].dx(i) for i in range(2)] +[u[i].dx(j) + u[j].dx(i) for (i,j) in [(0,1)]])

a = (( inner(eps(v),D_mateps(u)))dx )
L = dot(gR,v)*ds(3)
u = Function(V)

solve(a==L,u,dirichlet_u)

plot(u,interactive = True)
pointA = (1,1)
sigma_w = D_mat*eps(u)
X = VectorFunctionSpace(mesh,'DG',0)
sigma = project(sigma_w,X)
[s0,s1,s2] = sigma.split(True)
y = s0(pointA)

---

Just try:

pointA = Point(0.2,0.2)
sigma_w = D_mat*eps(u)
X = VectorFunctionSpace(mesh,'DG',0,dim=3)
sigma = project(sigma_w,X)
[s0,s1,s2] = sigma.split(True)
y = s0(pointA)
plot(s0, interactive=True)

---

Thanks, It works.