I am trying to get hierarchical FEMxDEM coupling to work, similar to Guo & Zhao (2014). I am using YADE for the DEM part.
The idea of this technique is to get the stiffness or tangent operator from the DEM simulations at every Gauss (quadrature) point as a replacement for having to specify an explicit constitutive law. Example code using Escrypt and YADE can be found here and here.
I have managed to get everything to work for FEniCSx with YADE except for one part. The DEM simulations provide the stiffness tensor and stresses at the Gauss points, but I need to project them back to the nodes before calling any type of ‘solve’ on the FEM level. (My previous mistake was directly assigning the solutions at the Gauss points to the nodes, which clearly results in strange results.)
The method I currently use (same as in the paper) is an implicit/iterative scheme using a Newton-Raphson method (similar to this one) to reach a stable displacement u. Ideally, I won’t have to project the stresses from the Gauss points to the nodes at each and every (implicit) iteration that is tuning du.
How do I project all necessary quantities (e.g. stresses and stiffness tensor) from the Gauss point back to the nodes? Is there a way to do this efficiently without doing many projections per time step?
Let me know if you require any further information. Unfortunately, even a MWE is rather involved.
Your help is much appreciated.
With kind regards,
Danny
Hi,
in the demo you mentionned, you can see that the stress is defined as a Quadrature function, meaning that its array of values are directly the values at the Gauss points and not the nodes. You can do the same for the tangent stiffness. You can thus directly update their values with what comes out of your DEM. No projection should be involved. You just have use the same quadrature rule and figure out the correct ordering in each element.
In msFEM2D.py it indeed seems the case that they directly use stress quadrature functions. However, in biaxialSmooth.py they do seem to use “proj = Projector(mydomain)” and “sig = proj(stress) # project Gauss point value to nodal value” for every loading step. I suppose that this is then only to print / visualise the results.
In FEniCSx, how do I make sure that the defined fem.Function variables (such as the tangent stiffness) are defined at the Gauss points and that the solver uses only these values for iteratively improving du?
In my current scheme, I defined the fem.Functions for e.g. the stress like
We = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(gdim,gdim), scheme="default", degree=deg_quad)
and directly assigned the Gauss point values to these fem.Function on each update. Unfortunately, it looks like the solves may be update u at the node values, since it looks as if the deformation is being underestimated. Could you hint me in the right direction?