How to get value from sym(grad(uh))

General
1

Hello everyone.
I would like to ask that how can I access the value of the Strain, Stress as the code below?
I can check their shape using ufl_shape which are (2,2). But i dont know how to access the value of each element in the whole domain.
If you know the answer for this, please help me.
Thank you so much.

import numpy as np
from mpi4py import MPI

from dolfinx import default_scalar_type, plot
from dolfinx.fem import (Constant, Function, functionspace, dirichletbc,
                         Expression,locate_dofs_topological)
from dolfinx.mesh import create_rectangle, meshtags, CellType, locate_entities_boundary
from dolfinx.fem.petsc import LinearProblem

from ufl import (Measure, TestFunction, TrialFunction, Identity,
                 sqrt, dx, tr, grad, inner, sym)

L, H, Nx, Ny = 2., 1., 10, 5
domain = create_rectangle(MPI.COMM_WORLD, [[0,0], [L,H]], [Nx, Ny], CellType.triangle)
V = functionspace(domain, ("Lagrange", 2, (domain.geometry.dim,)))

E, nu = default_scalar_type(1.0e3), default_scalar_type(0.3)
mu = Constant(domain, E / (2 * (1 + nu)))
lam = Constant(domain, (E * nu) / ((1 + nu) * (1 - 2 * nu)))

def eps(u) :
    return sym(grad(u))

def sigma(u):
    return lam*tr(eps(u))*Identity(len(u)) + 2.0 * mu*eps(u)

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], L)

fdim = domain.topology.dim - 1
left_facets = locate_entities_boundary(domain, fdim, left)
right_facets = locate_entities_boundary(domain, fdim, right)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [dirichletbc(u_bc, left_dofs, V)]

T = Constant(domain, default_scalar_type((0, -1)))

u, v = TrialFunction(V), TestFunction(V)

metadata = {"quadrature_degree": 2}
ds = Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)

a = inner(sigma(u), eps(v)) * dx 
L1 = inner(v, T) * ds(2)

problem = LinearProblem(a, L1, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

Strain = eps(uh)
Stress = sigma(uh)
2

Thank you so much for your reading. I have solved the problem myself as the code below:

Stress_expr = Expression(Stress[0,0], V.element.interpolation_points())
Stresses = Function(V)
Stresses.interpolate(Stress_expr)

If you know any better ways, please let me know.
Thank you so much!

3

See for instance:
https://jsdokken.com/fenics22-tutorial/heat_eq.html#post-processing-without-projections
You can also use dolfinx.fem.Expression directly to evaluate Strain or Stress in a given point in any element dolfinx/python/test/unit/fem/test_expression.py at 9f3a1b5c4bb01b0c547340fa9881082a804e3aa5 · FEniCS/dolfinx · GitHub

4

Thank you so much :smiley:

5 
Stress_expr = Expression(Stress[0,0], V.element.interpolation_points())
Stresses = Function(V)
Stresses.interpolate(Stress_expr)

Why is [0,0] used here?

6

As he just wants to get the 0th component of the stress.
If you want to get the full tensor, you need to create a space with a blocked element, for instance
basix.ufl.element("Lagrange", mesh.topology.cell_name(), 0, shape=Stress.ufl_shape(), discontinuous=True)
or with a TensorElement("DG", mesh.topology.cell_name(), 0) if you are using an older version of DOLFINx.

2 Likes