Stress calculation for user defined materials

variational formulation
1

Hi Fenics community,

I am working on the implementation of a hyperelastic material model that needs operations like, for example, integrations over the unit sphere surface in order to get the stress in response to the deformation (gradient). Just something (like a UMAT or used defined material) of the form:

def getstress(F)
    ...
    return P

For the evaluation inside of this function, it would be convenient to use packages like quadpy for quadratures etc. However, the input F (being calculated from the displacement – the primary unknown u) is a UFL expression and cannot be treated like a numpy array or so.
Is there a “standard treatment” or suggestion how to handle this, i.e. does the entire workflow from F to P have to be via UFL expressions and corresponding UFL functions or is it possible to convert it for this concern into other formats to make it available to e.g. numpy functions or packages like the mentioned quadpy (or would this–if possible–make the computations super inefficient)?

Any help or suggestions would be much appreciated. Thanks in advance!

2

I’m not really sure what you want to do here.

In DOLFINx, you can use dolfinx.fem.Expression to evaluate a ufl expression at any point, see:
https://docs.fenicsproject.org/dolfinx/v0.4.1/python/demos/demo_elasticity.html#post-processing
or

The input to the expression operator can be any set of points, including points from quadpy (if they are formatted as a 2D array of shape (num_points, tdim).

You can also create quadrature elements (however, currently only with the quadrature schemes supported in basix (see basix — Basix 0.4.2 documentation for the rules).
However, this should hopefully be extended soon:
Support custom Quadrature-rules in QuadratureElement · Issue #510 · FEniCS/ffcx · GitHub

1 Like
3

Hi Jorgen,

Thanks a lot for your help and apologies for not making my point clear enough.

Essentially, I just want to code a constitutive material law that might contain some more computations at the material point level. Meaning, to replace, for example, the stress definition in your hyperelasticity tutorial:

#P = ufl.diff(psi, F)
P = getstress(F)

and calculating P via a separate function in which I might use functions of the package quadpy. The function might contain something like (very simplified and not physically meaningful):

def getstress(F):
   lambda = ufl.variable(ufl.tr(C))
   def a(theta):
      return sin(theta)*lambda
   return quadpy.quad(a,0,pi)*F

The function dolfinx.fem.Expression you suggested sounds suitable to “extract” the value of the UFL expression and to do further computations with, for example, numpy. I just wanted to ask if there is something like a best practice for the above mentioned task.

4

As you want to use P in a ufl formulation, you should use a quadrature element matching the quadrature rule of quadpy, as the degrees of freedom will simply be point evaluations of P at these points.

5

Hi @dokken ,

Thanks for your help. I have a follow-up question regarding this topic. Your answer emphasised that one should stay in the “ufl-space” for the way from deformation gradient F to stress P (so that one gets out P as ufl expression). Leaving aside my specific question with quadpy, is there a general suggestion or best-practice how to handle implementations where P relies on some expressions that cannot be written with ufl-compatible functions? In the extreme scenario some “black-box” behaviour with a function that takes F and gives P (maybe from an external library) or for example a neural net that defines the constitutive relation.

6

You could use the quadrature space. It is a space where data can only be evaluated at quadrature points, and thus is suitable for data that does not fall into a Finite element space.

7

There are some demonstrations by @michalhabera in dolfiny which could be a good starting point.

1 Like