Obtaining the maximum strain energy over the history


I have the following weak form of the pde:

R_z = y*2*z*(psi_max)*dx+ k*(y*(z-1)/epsilon + epsilon*inner(grad(z),grad(y)))*dx

Here z is the trial function, y is the test function and k,epsilon are material parameters. psi_max is the maximum strain energy over the entire simulation history (I am applying load in steps). Does anyone have any ideas on how I can keep track of the maximum strain energy (at quadrature points essentially)?

I suppose I need to update the values of psi_max after each load step, but I am not quite sure how to do that.

Thank you!

Quadrature elements are currently a bit buggy, but can be used by using some deprecated functionality as follows:

from dolfin import *
from ufl import Max

# The following shuts off a deprecation warning for quadrature representation:
import warnings
from ffc.quadrature.deprecation \
    import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)

# A consistent quadrature degree must be used throughout, to avoid mismatch
# between the quadrature space and quadature of the variational form:  
dx = dx(metadata={"quadrature_degree":QUAD_DEG})

mesh = UnitSquareMesh(10,10)

# Function space with DoFs at quadrature points:
Qelement = FiniteElement("Quadrature", mesh.ufl_cell(),
Q = FunctionSpace(mesh, Qelement)

# Function to store previous max strain energy:
psi_max_old = Function(Q)

# Placeholder, to be replaced with whatever the expression for strain
# energy is in the problem:
psi = SpatialCoordinate(mesh)[0]

# Use this in the variational problem where psi_max is needed:
psi_max = Max(psi,psi_max_old)

# After solving in each load step, update psi_max_old.  When the unknowns of
# a variational problem are in a Quadrature-type space, you need to use the
# deprecated quadrature representation.  

# Plot (projection onto linears) to verify:
from matplotlib import pyplot as plt

Also check out Jeremy Bleyer’s plasticity tutorial, which uses quadrature spaces:


1 Like

Thank you David for your reply. I tried what you suggested. The updating of psi_max_old was really slow, so I tried to do the local projection as mentioned in the Jeremy Bleyer’s tutorial:

def local_project(v, V, u=None):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dxm
    b_proj = inner(v, v_)*dxm
    solver = LocalSolver(a_proj, b_proj)
    if u is None:
        u = Function(V)
        return u

#Call local_project function to update psi_max_old
local_project(psi_max, Q, psi_max_old)

but this gives me a compilation error. Do you have any ideas on what could be going wrong?

Thanks again!

If you add

parameters["form_compiler"]["representation"] = 'quadrature'

after importing dolfin it should work. The disadvantage of this is that quadrature representation will be used for all forms, and is less robust than the default ("uflacs") for very complicated forms. I usually try to use quadrature representation only where necessary.

Another trick I’ve found useful is to use an iterative method with a Jacobi preconditioner in a global L^2 projection to quadrature spaces, such that the preconditioner is the inverse of the (diagonal) matrix.

1 Like

Hi kamensky!

would you please tell me how you “usually try to use quadrature representation only
where necessary”?
For instance I want to use

import dolfin as df
df.parameters["form_compiler"]["representation"] = "uflacs"

for my global problem. But for any projection (involving quadrature spaces) I would like to

df.parameters["form_compiler"]["representation"] = "quadrature"

in the df.LocalSolver of the local_project function defined as above.

Do you know how to do that?


Maybe try something like this:

def local_project(v, V, u):
    parameters["form_compiler"]["representation"] = 'quadrature'
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    parameters["form_compiler"]["representation"] = 'uflacs'

local_project(psi_max, Q, psi_max_old)

Ah okay. I thought there might be an option to specify form compiler parameters to
the LocalSolver directly, but I guess this should also work.


Hi, your discussion is very helpful to me! But I’m not sure how to achieve it in the new version of dolfinx. For example, how to update the following line


Thank you!

you can now directly interpolate UFL expressions, see here:

Thanks for your reply. I tried (after solving after each load step)


But it generated an error:

RecursionError: maximum recursion depth exceeded in __instancecheck__

Do you have any suggestions? Thanks again!

Please provide a code that can reproduce this error, as it is impossible to guide you without the context of the error message.

Hi dokken, here is a code reproducing this recursion error.

import numpy as np
import dolfinx.fem as fem
import dolfinx.mesh as mesh
import ufl
from mpi4py import MPI

domain = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                             points=[[0,0], [1,1]],
                             n=[16,16], cell_type=mesh.CellType.quadrilateral)
ndim = domain.geometry.dim
fdim = ndim-1

element_u = ufl.VectorElement('Lagrange', domain.ufl_cell(), degree=1, dim=2)
V_u = fem.FunctionSpace(domain, element_u)
element_alpha = ufl.FiniteElement('Lagrange', domain.ufl_cell(), degree=1)
V_alpha = fem.FunctionSpace(domain, element_alpha)

# Measures
metadata = {"quadrature_degree": QUAD_DEG}
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
ds = ufl.Measure("ds", domain=domain, metadata=metadata)

u = fem.Function(V_u)
psi = 0.5 * ufl.inner(ufl.sym(ufl.grad(u)), ufl.sym(ufl.grad(u)))
# Function to store previous max value
psi_max_old = fem.Function(V_alpha)
# Use this in the variational problem where psi_max is needed
psi_max = ufl.max_value(psi, psi_max_old)

for i_t, t in enumerate(np.linspace(0, 0.01, 1)):
    # ... solve the problem ...

I guess that I made a mistake with the function space of psi_max_old. But I have no idea how to get through.


As seen in the example pointed out by @bleyerj, you need to wrap the psi_max as a dolfinx.fem.Expression before interpolating

Thanks for your help! I wrapped it up as follows and it worked.

V_psi_max = fem.FunctionSpace(domain, ("Lagrange", 1))
# After each load step
psi_max_expr = fem.Expression(psi_max, V_psi_max.element.interpolation_points())