Variable material coefficients as quadrature functions (in plasticity)

Hi there!

I am currently implementing classical elastoplasticity involving two different materials (in contact), with different elastic constants, yield stresses, etc. The subdomains corresponding to each material are given by a meshfunction (inherited from gmsh).

It is well-known that variables like plastic deformation, etc. should be represented as quadrature functions and that the return mapping algorithm (RMA) should be applied at the quadrature points. (In fact, I implemented the RMA at the level of the quadrature functions’ numpy arrays.)

Since the material constants are part of the RMA it would be very convenient and efficient to implement the variable coefficients (which are, in fact, constant in each subdomain) as quadrature functions as well. Hence, I am interested in an efficient way of transforming the meshfunction and two constants (corresponding to each material) to a quadrature function.

you can use an expression from mesh domains (see for instance The expression can then be interpolated to a Quadrature function.

Thank you for the link which was very helpful. Previously, I had found only an outdated version of this tutorial.

With the help of the hints given there I found an efficient solution to my problem. Perhaps it is of some use to other people, too.

def meshfunc2quadfunc(subdomains,k0,k1):


    return local_project(k,W0)

Here, subdomains is the meshfunction giving 0 for subdomain 0 and 1 for subdomain 1 and k0, k1 are the respective material constants. The spaces P0, W0 and the function local_project are taken from the standard fenics plasticity example: