Enhance `SolidBodyNearlyIncompressible`

[adtzlr]

Currently, the secondary mixed fields $p$ and $\bar{J}$ are hardcoded as Q0-fields.

Try to generalize this to any kind of dual regions. The dual region should be obtained from FieldDual.

Preparation Tasks

• [x] Add a FieldDual(field) with the definitions from FieldsMixed. (see #629)
  • [x] ~Add it as method to Field.dual().~ Too risky because this does only work for template regions.
  • [x] Use FieldDual or Field.dual() inside FieldsMixed. (see #629)
• [x] Create a ThreeFieldVariationDistortional(material): A ThreeFieldVariation but simplified to $\psi(\boldsymbol{F}, p, \bar{J}) = \hat{\psi}(\boldsymbol{F}) + U(\bar{J}) + p\ (J - \bar{J})$. ~We could also use ThreeFieldVariation to start with but performance (runtimes) will be rather low.~ see #633
• [x] ~Create a helper StateDual with the dual fields.~ Probably not necessary.
• [ ] Enhance SolidBodyNearlyIncompressible or create a new class SolidBody...?. This should be very similar to ThreeFieldVariation but all fields are internally related to the displacement field. It more or less hides the ThreeFieldVariation behind the solid body class.

This is way more complicated... will take some time to be implemented ⏳

[adtzlr]

There must be some mistake which I do not see (yet). It is complicated to update the internal fields because they depend on the vector and the matrix. In the original implementation, the internal fields are updated on extract. This is not possible here!?

Probably some sign is wrong?

import felupe as fem
import numpy as np


class StateDual:
    def __init__(self, u=None, p=None, J=None):
        self.u = u
        self.p = p
        self.J = J


class SolidBodyX:
    def __init__(self, umat, field, bulk, statevars=None):
        self.umat = fem.NearlyIncompressible(umat, bulk=bulk)
        self.field = field

        self.results = fem.mechanics._helpers.Results(stress=True, elasticity=True)
        self.results.state = StateDual(
            u=field[0].copy(),
            p=fem.FieldDual(field.region),
            J=fem.FieldDual(field.region, values=1),
        )
        self.results.kinematics = self._extract(self.field)

        if statevars is not None:
            self.results.statevars = statevars
        else:
            statevars_shape = (0,)
            if hasattr(umat, "x"):
                statevars_shape = umat.x[-1].shape
            self.results.statevars = np.zeros(
                (
                    *statevars_shape,
                    field.region.quadrature.npoints,
                    field.region.mesh.ncells,
                )
            )

        self.assemble = fem.mechanics._helpers.Assemble(
            vector=self._vector, matrix=self._matrix
        )

        self.evaluate = fem.mechanics._helpers.Evaluate(
            gradient=self._gradient,
            hessian=self._hessian,
        )

        self._form = fem.IntegralForm

    def _extract(self, field):
        self.field = field
        self.results.kinematics = [
            *self.field.extract(),
            self.results.state.p.interpolate(),
            self.results.state.J.interpolate(),
        ]

        return self.results.kinematics

    def _gradient(self, field=None, args=(), kwargs={}):
        if field is not None:
            self.field = field
            self.results.kinematics = self._extract(self.field)

        gradient = self.umat.gradient(
            [*self.results.kinematics, self.results.statevars], *args, **kwargs
        )

        self.results.stress, self.results._statevars = gradient[:-1], gradient[-1]

        return self.results.stress

    def _hessian(self, field=None, args=(), kwargs={}):
        if field is not None:
            self.field = field
            self.results.kinematics = self._extract(self.field)

        self.results.elasticity = self.umat.hessian(
            [*self.results.kinematics, self.results.statevars], *args, **kwargs
        )

        return self.results.elasticity

    def _vector(self, field=None, parallel=False, items=None, args=(), kwargs={}):
        if field is not None:
            self.field = field

        self.results.stress = self._gradient(field, args=args, kwargs=kwargs)
        [fu, fp, fJ] = self._form(
            fun=self.results.stress[slice(items)],
            v=self.field[0] & self.results.state.p & self.results.state.J,
            dV=self.field.region.dV,
        ).assemble(parallel=parallel, block=False)

        from scipy.sparse.linalg import inv

        Kuu, Kup, KuJ, Kpp, KpJ, KJJ = self._matrix(self.field, reduce=False)
        self.results.force = (
            fu + Kup @ inv(KpJ) @ KJJ @ inv(KpJ).T @ fp - Kup @ inv(KpJ) @ fJ
        )

        du = self.field[0][:, None] - self.results.state.u[:, None]
        dJ = -inv(KpJ).T @ Kup.T @ du + inv(KpJ).T @ fp
        dp = -inv(KpJ) @ KJJ @ dJ + inv(KpJ) @ fJ

        self.results.state.u += du
        self.results.state.p += dp
        self.results.state.J += dJ

        return self.results.force

    def _matrix(
        self, field=None, parallel=False, items=None, reduce=True, args=(), kwargs={}
    ):
        if field is not None:
            self.field = field

        self.results.elasticity = self._hessian(field, args=args, kwargs=kwargs)

        form = self._form(
            fun=self.results.elasticity[slice(items)],
            v=self.field[0] & self.results.state.p & self.results.state.J,
            u=self.field[0] & self.results.state.p & self.results.state.J,
            dV=self.field.region.dV,
        )

        self.results.stiffness_values = form.integrate(
            parallel=parallel, out=self.results.stiffness_values
        )

        [Kuu, Kup, KuJ, Kpp, KpJ, KJJ] = form.assemble(
            values=self.results.stiffness_values, block=False
        )
        
        if not reduce:
            return Kuu, Kup, KuJ, Kpp, KpJ, KJJ

        from scipy.sparse.linalg import inv

        self.results.stiffness = Kuu + Kup @ inv(KpJ) @ KJJ @ inv(KpJ).T @ Kup.T

        return self.results.stiffness


field = fem.FieldsMixed(fem.RegionHexahedron(fem.Cube(n=6)), n=1)
boundaries, loadcase = fem.dof.uniaxial(field, clamped=False)
umat = fem.NeoHooke(mu=1)
solid = SolidBodyX(umat, field, bulk=5000)
step = fem.Step(items=[solid], boundaries=boundaries)
job = fem.Job(steps=[step]).evaluate(verbose=2)

[adtzlr]

Finally, this works for disconnected dual meshes in #640 (draft). For connected meshes, the volume-matrix can't be inverted (singular).

Considerations on merging #640

Is this generalization really useful if it doesn't work for connected dual meshes?