Stress projection about felupe HOT 3 CLOSED

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())

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)

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)