How to compute the strain & stress matrices?

1

Hello everyone,

I am new to FEniCS, and I’m currently struggling on the following problem. I simply want to print (if possible) the strain & stress matrices for the elasticity problem at each point. I worked with the linear elasticity example, and for now I added the following lines of code.


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u, lambd, mu):
    # mu and lambd are lame's parameters
    d = u.geometric_dimension()  # space dimension
    return lambd*tr(epsilon(u))*Identity(d) + 2*mu*epsilon(u)

# Compute strain & stress
T = TensorFunctionSpace(mesh, 'P', 2)
stress = sigma(u, lambd, mu)
strain = epsilon(u)
s1 = project(stress, T)
s2 = project(strain, T)

Now I know I can obtain stress & strain using for example np.array(s1.vector()), but it returns a lot of values (35 * 9) not ordered at all. I’d like to directly obtain the matrices at each point, or to understand how they are ordered to recreate these matrices.
I don’t understand much to FEniCS yet so it is highly possible that what I’m asking is really simple or in the contrary not possible, I don’t know.

And second question, how can I print the coordinates of all points ? (meaning the points in the mesh, but also the Gauss points)

2

Hi, you can just do s1(x,y,z)at the point coordinate of your choice, assuming that you are in 3D. It should return 9 values ordered as follows:

[[0, 1, 2],
 [3, 4, 5],
 [6, 7, 8]]

If you print s1.vector().get_local() you get a vector of 9 values at the first dof vertex, then 9 values at the second, etc… Beware that there is a reordering of the vertex numbering of your mesh and the dof numbering. You can access the dof coordinates using

T.tabulate_dof_coordinates()

The mesh coordinates can be accessed from mesh.coordinates().

3

Oh so I was almost there already … well thanks :slight_smile:

And thanks for the mesh.coordinates(), but how can I also display the Gauss points ?

4

This is somehow hidden in FEniCS. It will depend on how you defined your integration measures, in the bilinear form for instance. Gauss points are not attached to a mesh or an interpolation space but to an integration measure such as dx. Without specifying anything, FEniCS will automatically deduce the required quadrature degree for the form you have defined.
If you have a “standard” bilinear form such as a = dot(grad(u), grad(v))*dx and u, v are of degree 2 let’s say, then it will be enough to compute the bilinear form with a degree 2 quadrature (i.e. 3 Gauss points on a triangle).
You can pass yourself the required degree d by using dx(metadata={"quadrature_degree": d}).
Why do you want to display the Gauss points exactly ?

If you want to see where are the Gauss points of degree d=2 for instance on a 2D mesh you can use a Quadrature function space with the corresponding degree such as:

from dolfin import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(1,1, "crossed")

Ve = FiniteElement("Quadrature", mesh.ufl_cell(), degree=2, quad_scheme='default')
V = FunctionSpace(mesh, Ve)

plot(mesh)
dofs = V.tabulate_dof_coordinates()
for (i, dof) in enumerate(dofs):
    plt.plot(*dof, "xk")
    plt.annotate(str(i), xy=dof, xytext=dof+0.02)
5

Thank you for your answer

I am comparing the results (coordinates) of a simulation software with Fenics, to check the accuracy. And I want to have the Gauss points to see if they are the same, and if not I’d like to see how it might influence my results.

6

Hi,

is there a way of doing this in Fenicsx? As far as I know the ‘project’ function is not there anymore in Fenicsx.

7

A projection is simply solving
\int_\Omega u \cdot v ~\mathrm{d}x = \int_\Omega f \cdot v ~\mathrm{d} x
which is easy to implement in DOLFINx:

A more scalable approach than the route taken in DOLFINy is to create a projector class, that caches the mass matrix, so that you do not have to re-assemble (and refactorize in the case of LU) if you use it multiple times).

For the quadrature points (in physical space), consider the first bit of: Numerical values from ufl.SpatialCoordinate - #2 by dokken

8

Hi dokken,

I tried to use your code to extract the strains for all points of my mesh. I don’t know if I have done it correctly.If I want to extract the strains matrices at each point, what kind of modification should I do to the code?

I would appreciate anyhelp.

# Linear Elasticty Problem

from petsc4py.PETSc import ScalarType
import ufl
import numpy as np
from dolfinx.fem import *
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.fem import Function, FunctionSpace, Constant
from ufl import TestFunction, dx, inner

L = 10.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[L,1,1]],[10,5,5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))

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 = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

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 = mesh.meshtags(domain, fdim, marked_facets[sorted_facets],marked_values[sorted_facets])

u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V,facet_tag.dim, facet_tag.find(1))
bc_l = fem.dirichletbc(u_bc, left_dofs, V)


right_dofs = fem.locate_dofs_topological(V.sub(0),facet_tag.dim, facet_tag.find(2))
bc_r = fem.dirichletbc(ScalarType(1.0),right_dofs,V.sub(0))
bcs =  [bc_l, bc_r]

L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

def eps(u):
    return ufl.sym(ufl.grad(u))
def sig(u):
    return lambda_* ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*eps(u)




u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0,0,0)))
a = ufl.inner(sig(u), eps(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx #+ ufl.dot(T,v) * ds

Problem = fem.petsc.LinearProblem(a, L, bcs =bcs, petsc_options={"ksp_type":"preonly", "pc_type":"lu"})
uh = Problem.solve()

# Projection of Strain Tensor

import dolfinx
import dolfinx.io

from mpi4py import MPI
from petsc4py import PETSc

epsilon = eps(uh)
V_eps = fem.TensorFunctionSpace(domain,("DG",0))
strain_expr = fem.Expression(epsilon,V_eps.element.interpolation_points())
strains = fem.Function(V_eps)
strains.interpolate(strain_expr)

# Project DG function to CG space
W= dolfinx.fem.TensorFunctionSpace(domain, ("CG",1))
u, w = ufl.TrialFunction(W), ufl.TestFunction(W)
a= dolfinx.fem.form(ufl.inner(u, w) * ufl.dx)
l = dolfinx.fem.form(ufl.inner(strains, w) * ufl.dx)

eps_h = dolfinx.fem.Function(W)

bcs = []
A = dolfinx.fem.petsc.assemble_matrix(a,bcs)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector(l)

solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.setOperators(A)
solver.solve(b, eps_h.vector)
eps_h.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

print(eps_h.vector.array)

#