Stress Projection for `TetraMINI` is complicated

adtzlr commented 7 months ago

This is not a classical code-bug, but it could be a bug from a user perspective.

Summary (for v7.19.1)

Adopt the README code-block to the MINI-element formulation ✅
Force-displacement curve is realistic ✅
Stresses at quadrature points are straightforward ✅
Stress-projection is error-prone for MINI-elements ⚠️

Possible Solutions

~Add a new argument to RegionTetraMINI(mesh, project=True) which raises the order of the quadrature scheme from 2 to 5!? This is a manual solution.~
~Raise the quadrature schemes of the template regions to always produce correct results for projection? This affects the performance (more computationally expensive assembly)!~
[x] Fall-back to tools.extrapolate(mean=True) for the projection of cell-means? --> This seems to be very promising compared to the results obtained with hexahedrons, see https://github.com/adtzlr/felupe/issues/690#issuecomment-1987001042! mean=True should enforce extrapolate() independent of the shape of the input array. The project method should enforce regions with quadrature schemes which are good for projection and if the shape of the input array is incompatible then an error is raised. This is implemented in PR #691

Detailed Description

The first part is easy...

import felupe as fem
from pypardiso import spsolve

mesh = fem.Cube(n=(16, 6, 6)).triangulate(3).add_midpoints_volumes(cell_type="tetra")
region = fem.RegionTetraMINI(mesh, fem.TetrahedronQuadrature(order=2))
field = fem.FieldsMixed(region)

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)

nh = fem.NeoHooke(mu=1)
solid = fem.SolidBody(fem.NearlyIncompressible(nh, bulk=500), field)

move = fem.math.linsteps([0, 2], num=3)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(solver=spsolve, parallel=True)
fig, ax = job.plot(
   xlabel="Displacement $u$ in mm $\longrightarrow$",
   ylabel="Normal Force $F$ in N $\longrightarrow$",
)

...this is more or less the README code-block ✅. Obtaining the Cauchy-stress at quadrature points is also straightforward ✅. But the projection produces weird stresses ⚠️ at the midvolume-points.

import numpy as np

# stress at quadrature points
cauchy = solid.evaluate.cauchy_stress(field)

# project the stress to mesh-points
stress = fem.project(cauchy, region)
stress[-mesh.ncells:] = np.nan

# view cell-means of stress xx
view = solid.view(cell_data=dict(Stress=np.mean(cauchy, axis=-2).T))
view.plot("Stress", show_undeformed=False).show()

# view projected point-data of stress xx
view = solid.view(point_data=dict(Stress=stress))
view.plot("Stress", show_undeformed=False).show()

So these points have to be ignored 🗑️. But the result is still very different to the cell-data (cell-means of stresses).


# ignore projected stresses at midvolume-points
stress[-mesh.ncells:] = np.nan

# set displacements of mid-volume points to zero
field[0].values[-mesh.ncells:] = 0

# view projected point-data of stress xx
view = solid.view(point_data=dict(Stress=stress))
view.plot("Stress", show_undeformed=False).show()

For RegionTetraMINI a quadrature scheme with order=5 has to be used to project values to mesh-points (order=2 is too low and the next implemented order is fifth-order). Repeat the above code-blocks with:

# ...
region = fem.RegionTetraMINI(mesh, fem.TetrahedronQuadrature(order=5))
# ...

~Another option would be to project the cell-means to the mesh-points (also use a fifth-order quadrature scheme!).~ This is probably not a good idea.

# project the stress to mesh-points
stress = fem.project(cauchy, region, mean=True)
stress[-mesh.ncells:] = np.nan

# view projected point-data of stress xx
view = solid.view(point_data=dict(Stress=stress))
view.plot("Stress", show_undeformed=False).show()

A similar task could also be done directly within the PyVista plotter. This seems to be the easiest operation. Same results may be obtained by extrapolate(mean=True).

view = solid.view(cell_data=dict(Stress=np.mean(cauchy, axis=-2).T))
view.mesh = view.mesh.cell_data_to_point_data()
view.plot("Stress", show_undeformed=False).show()

# alternative with equal result
stress = fem.tools.extrapolate(cauchy, region, mean=True)  # this is the old project-method
view = solid.view(point_data=dict(Stress=stress))
view.plot("Stress", show_undeformed=False).show()

Results obtained by a finer mesh:

(extrapolate cell-mean to point-data)

I have absolutely no idea how to simplify or make the stress projection for MINI-elements less error-prone. In the long term, we could separate global DOF from private DOF and then only the global DOF are considered for projection / plotting. Triangle and Tetra - Regions should be initiated with a higher-order quadrature scheme if it is used for interpolation or projection (and not only gradient evaluations). Something like Region(project=True) instead of manually setting the quadrature scheme.

adtzlr commented 7 months ago

Here are the results obtained with hexahedrons:

and with a finer mesh

and with an even finer mesh

and with an even even (...😄) finer mesh

adtzlr commented 7 months ago

With #691, an error is raised if the quadrature scheme's order is too low for projection.

import felupe as fem
from pypardiso import spsolve

mesh = fem.Cube(n=(11, 4, 4)).triangulate(3).add_midpoints_volumes(cell_type="tetra")
region = fem.RegionTetraMINI(mesh, fem.TetrahedronQuadrature(order=2))
field = fem.FieldsMixed(region)

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)

nh = fem.NeoHooke(mu=1)
solid = fem.SolidBody(fem.NearlyIncompressible(nh, bulk=500), field)

move = fem.math.linsteps([0, 2], num=3)
step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"])
job.evaluate(solver=spsolve, parallel=True)
fig, ax = job.plot(
   xlabel="Displacement $u$ in mm $\longrightarrow$",
   ylabel="Normal Force $F$ in N $\longrightarrow$",
)

import numpy as np

# stress at quadrature points
cauchy = solid.evaluate.cauchy_stress(field)

# project the stress to mesh-points
stress = fem.project(cauchy, region)

Hopefully this is useful and prevents some mistakes.

ValueError: Quadrature not supported (order is too low). Take the cell-means with `project(mean=True)` or use a higher-order quadrature scheme.

adtzlr / felupe