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.

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: