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.