Stress projection

[adtzlr]

The stress project(values, region) function is not well implemented. This must be enhanced. Currently, either the mean-cell values are averaged at the points or for some quadrature rules an extrapolation is done. Optionally, no averaging at points is done for disconnected meshes. The averaging part is fine, however, the projection from the quadrature to the mesh points should be changed.

Tasks

• rename project() to extrapolate(), see Add the L2-projection and replace it with the old tools.project() #683
• add a new project() function which works differently but supports the old arguments as well, see Add the L2-projection and replace it with the old tools.project() #683

[adtzlr]

Projection

A best-fit approximation is used to project values at quadrature-points to mesh-points.

$$ \int_v (y - x) \delta x \ dv = 0 $$

where $y$ is known at evaluated quadrature-points and the values $\hat{x}$, located at the mesh-points, are being looked for. The result is either per-cell (integrated and inverted per-cell - but not assembled) for discontinuous approximations across cells

$$ \hat{\boldsymbol{x}} = \left(\int_v \boldsymbol{h} \otimes \boldsymbol{h}\ dv \right)^{-1} \int_v \boldsymbol{h} \otimes y \ dv$$

or assembled as (system) sparse-matrix for continuous approximations.

$$ \left(\int_v \boldsymbol{h} \otimes \boldsymbol{h}\ dv \right) \hat{\boldsymbol{x}} = \int_v \boldsymbol{h} \otimes y \ dv$$

Note

The quadrature-scheme of the original region, used to evaluate $y$, must also be used for the evaluation of the shape (or ansatz / basis) functions $h$ of $x = \boldsymbol{h} \cdot \hat{\boldsymbol{x}}$.

import felupe as fem
import numpy as np

mesh = fem.Cube(n=11)
region = fem.RegionHexahedron(mesh)
field = fem.Field(region, dim=3).as_container()
y = field.extract()[0]

axes = y.shape[-2:]
dim = np.product(y.shape[:-2])
dV = region.dV
fun = np.eye(dim).reshape(dim, dim, 1, 1)

x = fem.Field(region, dim=dim).as_container()
A = fem.IntegralForm([fun], v=x, dV=dV, u=x, grad_v=[False], grad_u=[False]).assemble()
b = fem.IntegralForm([y.reshape(dim, *axes)], v=x, dV=dV, grad_v=[False]).assemble()

import pypardiso

xh = pypardiso.spsolve(A, b.toarray())

[adtzlr]

It seems to be way more efficient to do this component by component, as it drastically reduces the sparse-matrix size.

import felupe as fem
import numpy as np

mesh = fem.Cube(n=31)
region = fem.RegionHexahedron(mesh)
field = fem.Field(region, dim=3).as_container()

y = field.extract()[0][0, 0]

dV = region.dV
fun = np.ones((1, 1))

x = fem.Field(region).as_container()
A = fem.IntegralForm([fun], v=x, dV=dV, u=x, grad_v=[False], grad_u=[False]).assemble()
b = fem.IntegralForm([y], v=x, dV=dV, grad_v=[False]).assemble().toarray()

import pypardiso

xh = pypardiso.spsolve(A, b)

[adtzlr]

This one is the continuous projection. Hopefully efficient enough 🚀.

• works for quads, hexahedrons and all the lagrange bi/tri quadratic versions
• test serendipity elements
• works for triangles with quadrature order=2
• works for quadratic triangles with quadrature order=5 (4 is not implemented)
• test tetrahedrons

Some Thoughts

• The reference triangle is defined in (0, 1) whereas the rectangle is in (-1, 1).
• The sum of the volume matrix is the volume, correct result.
• For cell-wise constant values, switch to the existing method (mean=True).

import numpy as np
from scipy.sparse.linalg import spsolve
from felupe.assembly import IntegralFormCartesian
from felupe import Field


def project(values, region, solver=spsolve):
    "Return given array of values at quadrature-points projected to mesh-points."

    shape = values.shape[:-2]  # tensor-components
    dim = np.prod(shape).astype(int)  # tensor-size (enforce int for empty tuple)
    axes = values.shape[-2:]  # trailing axes

    # lhs (volume matrix) is constant for all items
    v = u = Field(region, dim=1)  # 1d-field for lhs
    A = IntegralFormCartesian(np.ones((1, 1)), v=v, dV=region.dV, u=u).assemble()

    # field of unknowns with projected values
    x = Field(region, dim=dim)

    # solve with individual right-hand-sides (per tensor-component of values)
    b = IntegralFormCartesian(values.reshape(dim, *axes), v=x, dV=region.dV).assemble()
    x.values[:] = spsolve(A, b.toarray().reshape(-1, dim)).reshape(x.values.shape)

    return x.values.reshape(-1, *shape)


import felupe as fem


mesh = fem.Rectangle(n=21)
region = fem.RegionQuad(mesh)
field = fem.FieldAxisymmetric(region, dim=2).as_container()
dxdX = field.extract()[0]

values_projected = project(dxdX, region)