Strain/stress at gauss points of (in-)elastic problems in weak formulation

[loop]

Hey everybody,

I am new to Fenics but I think it is the most appropriate way to get the results I need instead of coding complete new FEM routines in python. My goal is to have 3 seperate Fenics routines at the end of the week, for linear elastic, non-linear elastic and linear in-elastic (history dependent) problems in the weak for. All of the routines should calculate strain and stress at the gauss points.

I started with the linear elastic demo ft06_elasticity. Of course feel free to reduce the problem to 2D if it is easier. A normal unitsqaure is also okay.

Since now I modified the code such that I can calculate the strain/stress at each mesh coordinate. But I could not modify it for Gauss Points. Any idea how to do that?

from __future__ import print_function
from fenics import *
import numpy as np

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

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

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

def sigma(u):
    return lambda_*div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u, title='Displacement', mode='displacement')


Vsig = TensorFunctionSpace(mesh,'DG',0)
Veps = TensorFunctionSpace(mesh, 'DG', 0)

PKTensor = project(sigma(u), Vsig)
GLFSTensor = project(epsilon(u), Vsig)


Strain = []
Stress = []
temp_vec = np.zeros(shape=[6, 1])
temp_vec2 = np.zeros(shape=[6, 1])

for i in range(mesh.coordinates().shape[0]):
    x = mesh.coordinates()[i,0]
    y = mesh.coordinates()[i,1]
    z = mesh.coordinates()[i,2]
    temp_vec[0] = GLFSTensor(x, y, z)[0]
    temp_vec[1] = GLFSTensor(x, y, z)[4]
    temp_vec[2] = GLFSTensor(x, y, z)[8]
    temp_vec[3] = 2*GLFSTensor(x, y, z)[5]
    temp_vec[4] = 2*GLFSTensor(x, y, z)[2]
    temp_vec[5] = 2*GLFSTensor(x, y, z)[1]

    temp_vec2[0] = PKTensor(x, y, z)[0]
    temp_vec2[1] = PKTensor(x, y, z)[4]
    temp_vec2[2] = PKTensor(x, y, z)[8]
    temp_vec2[3] = 2 * PKTensor(x, y, z)[5]
    temp_vec2[4] = 2 * PKTensor(x, y, z)[2]
    temp_vec2[5] = 2 * PKTensor(x, y, z)[1]

    Strain.append(temp_vec)
    Stress.append(temp_vec2)

print(Strain)
print(Stress)

[klunkean]

The inclusion of plasticity (or generally any constitutive model that needs to be solved at integration point level) is a lot more involved than one might expect from inspecting the linear elasticity demo.
What makes linear elasticity so easy to implement is that you can explicitly define the stresses as a function of the displacement gradient which is almost never possible in other cases such as plasticity.

For a first look into workarounds I recommend reading this tutorial on von-Mises plasticity in FEniCS.

To be honest: Since you are new to FEniCS I’d move away from the idea to be done with what you have planned by the end of the week.

[loop]

Thank you for your answer. Okay I admit that I overshoot the target with one week. Actually my idea is to have 3 different models i.e. linear elastic, non-linear elastic and inelastic/plastic for the same mesh. And I also know that it will be much harder to calculate it for the last one. Therefore I wanted to start with (non-)linear elastic problems. Later on I want to try to implement the plastic problem for the same mesh calculating strain and stresses at the same gaussian points as for the elastic problem.

Linear elastic (see code above): In a normal FEM code I would calculate the strain and stresses at the gaussian points by calculating the strain displacement matrices in the gaussian points and multiply it by it by the displacement. Stress results from the multiplication of the Youngs Modulus and the resulting strain.
Is there a way to calculate the strain displacement matrices in Fenics? Or is there another way? In my modified code I calculated the strain/stress at the nodes using the mesh-coordinates. This is only possible because it is already in the physical space. But Gaussian Points are defined in the parametric space.

[klunkean]

The way you proposed does not fit perfectly into the FEniCS design. The basic quantities are the primary degrees of freedom. In the case of a purely mechanical problem this would classically be the displacement vector at each node when you use standard linear Lagrange shape functions.

In the linear elastic case the stresses are a function of the displacement gradient and thus are piecewise constant for piecewise linear displacements. The simplest way to find your stresses in FEniCS is to solve your system for the displacement and then project the stress expression onto a piecewise constant function space.

The calculation of quantities at quadrature points can be done using quadrature elements as shown in the tutorial I linked in my answer above.

[loop]

Hey klunkean,

I took a deep look at the tutorial you told me. And I am trying to modify and play with it a little bit. I got a problem when I wanted to calculate the stress “sig” at a certain point e.g. (0,0). After using the local_project I am always getting the message:

evaluate_basis: Function not supported/implemented for QuadratureElement.

You have any idea how to resolve the problem?

[bleyerj]

Hi,
Quadrature elements cannot be evaluated at a given point since the fields are only defined at the quadrature points. For plotting or evaluation, you just need to project the stress on a suitable FunctionSpace like DG1 for instance, you can use local_project for this purpose.

[loop]

Hey Bleyerj,

If I cannot evaluate the quadrature elements at a given point, is it possible to get the evaluation on the quadrature points? I tried this by saving the quadrature points in a vector and then evaluate the quadrature elements on the entries of the vector, but I get the same error.

So how do I get the value of the stress/strain at each quadrature point?

[klunkean]

If you have

1. A list of all the quadrature point coordinates of your mesh, say for example something like

quadrature_points= [ [0.,1.],
                     [0.25,0.5],
                     ...        ]

and

2. Your stresses as a dolfin function created by projection onto a suitable function (e.g. DG1) space as @bleyerj suggested, say stress_DG1

You can evaluate the stresses at the quadrature points in a loop

for quad_point in quadrature_points:
      dfpt = Point(quad_point)
      stress_at_point = stress_DG1(dfpt)
      # then do whatever you like with that value

For more specific help maybe post some code snippet we can work with.

[loop]

Hi,

So I copy the loop of the tutorial and show you my modification which is not working:
First I evaluated the quadrature points by

gdim = mesh.geometry().dim()
xq = W.tabulate_dof_coordinates().reshape((-1, gdim))
# Remove the duplicates
xq0 = xq[W.sub(0).dofmap().dofs()]

After this I modified the end of the global Newton Raphson loop:

for (i, t) in enumerate(load_steps):
    loading.t = t
    A, Res = assemble_system(a_Newton, res, bc)
    nRes0 = Res.norm("l2")
    nRes = nRes0
    Du.interpolate(Constant((0, 0)))
    print("Increment:", str(i+1))
    niter = 0
    while nRes/nRes0 > tol and niter < Nitermax:
        solve(A, du.vector(), Res, "mumps")
        Du.assign(Du+du)
        deps = eps(Du)
        sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
        local_project(sig_, W, sig)
        local_project(n_elas_, W, n_elas)
        local_project(beta_, W0, beta)
        A, Res = assemble_system(a_Newton, res, bc)
        nRes = Res.norm("l2")
        print("    Residual:", nRes)
        niter += 1

Vstr = FunctionSpace(mesh, 'DG', 1)
strDG = local_project(sig, Vstr)
#Testing on random quadrature points
x = xq0[i, 0]
y = xq0[i, 1]

print(strDG(x,y))

But this is not working at all.

[klunkean]

What is the error message you get?

If I read your code correctly I see two problems:
You are trying to project a tensor valued object onto a scalar function space.
And you are not passing a dolfin.Point to the resulting function but two float arguments.

[loop]

Okay I changed the last part too:

Vstr = TensorFunctionSpace(mesh, 'DG', 1)
strDG = local_project(sig, Vstr)

for quad_point in xq0:
    dfpt = Point(quad_point)
    print(strDG(dfpt))

and I get the error: Shapes do not match: and .

[Sabarish_Gopi]

Hi,

Were you able to resolve the issue ? I’m facing the same problem.

I want to plot the stress values calculated using quadrature elements.

[bleyerj]

In the tutorial mentioned at the beginning of this post, the stress variable sig is represented as a vector of dimension 4 (\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{xy}), so you should project on:

Vstr = VectorFunctionSpace(mesh, "DG", 1, dim=4)

[Sabarish_Gopi]

I have been able to plot the 4 stress components separately. I want to use the stress values at an internal boundary in the domain for post processing.

How do I obtain the the value of a particular stress component at some x-y co-ordinate ?

[loop]

Is it possible to post your modified code, so the answer to the question can be found to others?

[bleyerj]

You can just do strDG(x, y)

[Sabarish_Gopi]

This is the snippet of code I had added to the one posted above

The change was from dim = 1 to dim = 4

Vstr = VectorFunctionSpace(mesh, “DG”, 1, dim=4)
strDG = local_project(sig, Vstr)

#For Plotting
import matplotlib.pyplot as plt
p = plot(strDG[1], mode=‘color’)
plt.colorbar§
plt.title(r"\sigma_{yy}",fontsize=26)
plt.show()

[bleyerj]

Hi,
an update on solving inelastic problems in FEniCS.
We have now included a FEniCS interface for the constitutive law code generator MFront in the MFrontGenericInterfaceSupport project. Documentation and commented demos are available here: Overview of the mgis.fenics module
This enables to solve problems involving complex non-linear constitutive equations inside FEniCS including plasticity of course.

[MingY]

Hi @bleyerj, I am a big fan of your Numerical tours of Computational Mechanics using FEniCS. Vary glad for the update of MFront project. Do you have a plan to create a damage model example in FEniCS with MFront? It’s would be great for me to have a example for ductile (or brittle) damage model of heterogenous in FEniCS. We have several projects on that model, but It’s challenge for me to create it myself.
Thanks,
Ming

[bleyerj]

Thanks for the enthusiasm,
you can find a phase-field example of brittle fracture in a homogeneous material. You can easily adapt that by specifying a Function instead of a Constant for the material properties.
You can also find in the MFront gallery some damage laws. In order to use them as non-reguarized local damage models you can simply adapt the small strain plasticity example for instance.