Implementation of the Mean-Dilatation technique for nearly-incompressible materials

[adtzlr]

as described by Bonet & Wood, section 8.6.5 (p.231). Try to create a solid body SolidBodyMeanDilatation SolidBodyNearlyIncompressible, which only needs a user material for the distortional part.

[adtzlr]

This is super-fast but eventually super-wrong:

from felupe.math import (
    det,
    dot,
    transpose,
    trace,
    inv,
    ddot,
    identity,
    cdya_ik,
    cdya_il,
    dya,
)

from pypardiso import spsolve as solver

# undocumented, untested workaround if multiple blas libaries are installed
import os
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"


class NeoHooke:
    def __init__(self, mu):
        self.mu = mu

    def function(self, extract):
        F = extract[0]
        mu = self.mu
        J = det(F)
        C = dot(transpose(F), F)
        W = mu / 2 * (J ** (-2 / 3) * trace(C) - 3)
        return [W]

    def gradient(self, extract):
        F = extract[0]
        mu = self.mu
        J = det(F)
        iFT = transpose(inv(F, J))
        P = mu * (F - ddot(F, F) / 3 * iFT) * J ** (-2 / 3)
        return [P]

    def hessian(self, extract):
        F = extract[0]
        mu = self.mu
        J = det(F)
        iFT = transpose(inv(F, J))
        eye = identity(F)
        A4 = (
            mu
            * (
                cdya_ik(eye, eye)
                - 2 / 3 * dya(F, iFT)
                - 2 / 3 * dya(iFT, F)
                + 2 / 9 * ddot(F, F) * dya(iFT, iFT)
                + 1 / 3 * ddot(F, F) * cdya_il(iFT, iFT)
            )
            * J ** (-2 / 3)
        )

        return [A4]


class Volumetric:
    def function(self, extract):
        F = extract[0]
        W = det(F)
        return [W]

    def gradient(self, extract):
        F = extract[0]
        detF = det(F)
        iFT = transpose(inv(F, detF))
        return [detF * iFT]

    def hessian(self, extract):
        F = extract[0]
        detF = det(F)
        iFT = transpose(inv(F, detF))
        A4 = detF * (dya(iFT, iFT) - cdya_il(iFT, iFT))
        return [A4]


import felupe as fem
import numpy as np

# create a hexahedron-region on a cube
mesh = fem.Cube(n=21)
region = fem.RegionHexahedron(mesh)
region0 = fem.Region(
    fem.mesh.convert(mesh, order=0), 
    fem.ConstantHexahedron(), 
    fem.GaussLegendre(order=0, dim=3)
)
V = region.dV.sum(0)

# add a mixed field container (with displacements)
displacement = fem.FieldContainer([fem.Field(region, dim=3)])
pressure = fem.Field(region0, dim=1)

# apply a uniaxial elongation on the cube
boundaries, lc = fem.dof.uniaxial(displacement, clamped=True, move=0.4)

nh = NeoHooke(mu=1)
vol = Volumetric()

bulk = 5000
F = displacement.extract()

for iteration in range(8):
    F = displacement.extract()
    p = pressure.interpolate()

    L = fem.IntegralForm(nh.gradient(F), displacement, region.dV)
    a = fem.IntegralForm(nh.hessian(F), displacement, region.dV, displacement)

    Li = L.integrate()
    ai = a.integrate(parallel=True)

    Lv = fem.IntegralForm(vol.gradient(F), displacement, region.dV)
    av = fem.IntegralForm(vol.hessian(F), displacement, region.dV, displacement)

    M = Lv.integrate()[0]
    G = av.integrate(parallel=True)[0]

    r = L.assemble([Li[0] + p * M])
    K = a.assemble([ai[0] + p * G + bulk / V * dya(M, M)], parallel=True)

    system = fem.solve.partition(displacement, K, lc["dof1"], lc["dof0"], r)
    ddisplacement = fem.solve.solve(*system, lc["ext0"], solver=solver)

    displacement += ddisplacement
    
    # extract incremental deformation gradient from incremental displacement field
    dF = fem.Field(region, dim=3, values=ddisplacement.reshape(-1, 3)).grad()
    
    # incremental pressure and update of pressure field
    dp = bulk * np.einsum("ijqc,ijqc,qc->c", vol.gradient(F)[0], dF, region.dV) / V
    pressure += dp
    
    norm = fem.math.norm(ddisplacement)
    print(iteration, norm)

    if norm < 1e-4:
        break

fem.save(region, displacement)

[adtzlr]

Now, the additional vectors/matrices from the volumetric constraint have to be included in a new solid body definition. The pressure field has to be generated and updated inside.

import numpy as np

from .._field import Field
from .._assembly import IntegralFormMixed
from ..constitution import AreaChange, VolumeChange
from ..math import dot, transpose, det, dya
from ._helpers import Assemble, Evaluate, Results


class SolidBodyMeanDilatation:
    "A SolidBody with methods for the assembly of sparse vectors/matrices."

    def __init__(self, umat, field, bulk, regularize=False):

        self.umat = umat
        self.field = field
        self.bulk = bulk
        self.regularize = regularize

        self._volume_change = VolumeChange()
        self.volume = self.field.region.dV.sum(0)

        self.results = Results(stress=True, elasticity=True)
        self.results.kinematics, self.results.p = self._extract(self.field)

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

        self.evaluate = Evaluate(
            gradient=self._gradient,
            hessian=self._hessian,
            cauchy_stress=self._cauchy_stress,
            kirchhoff_stress=self._kirchhoff_stress,
        )

        self._form = IntegralFormMixed

    def _vector(
        self, field=None, parallel=False, jit=False, items=None, args=(), kwargs={}
    ):

        if field is not None:
            self.field = field

        self.results.stress = self._gradient(field, args=args, kwargs=kwargs)

        self.results.force = self._form(
            fun=self.results.stress,
            v=self.field,
            dV=self.field.region.dV,
        ).assemble(parallel=parallel, jit=jit)

        return self.results.force

    def _matrix(
        self, field=None, parallel=False, jit=False, items=None, args=(), kwargs={}
    ):

        if field is not None:
            self.field = field

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

        dJdF = self._volume_change.gradient(self.results.kinematics)
        H = self._form(fun=dJdF, v=self.field, dV=self.field.region.dV,).integrate(
            parallel=parallel, jit=jit
        )[0]

        bilinearform = self._form(
            fun=self.results.elasticity,
            v=self.field,
            u=self.field,
            dV=self.field.region.dV,
        )

        values = [
            bilinearform.integrate(parallel=parallel, jit=jit)[0]
            + self.bulk / self.volume * dya(H, H)
        ]

        self.results.stiffness = bilinearform.assemble(
            values=values, parallel=parallel, jit=jit
        )

        return self.results.stiffness

    def _extract(self, field):

        self.field = field
        self.results.kinematics = self.field.extract()
        mesh = self.field.region.mesh.copy()
        mesh.points += self.field[0].values
        v = type(self.field.region)(mesh).dV.sum(0)
        J = v / self.volume
        self.results.p = self.bulk * (J - 1)

        return self.results.kinematics, self.results.p

    def _gradient(self, field=None, args=(), kwargs={}):

        if field is not None:
            self.field = field
            self.results.kinematics, self.results.p = self._extract(self.field)

        dJdF = self._volume_change.gradient(self.results.kinematics)

        self.results.stress = [
            self.umat.gradient(self.results.kinematics, *args, **kwargs)[0]
            + self.results.p * dJdF
        ]

        return self.results.stress

    def _hessian(self, field=None, args=(), kwargs={}):

        if field is not None:
            self.field = field
            self.results.kinematics, self.results.p = self._extract(self.field)

        d2JdFdF = self._volume_change.hessian(self.results.kinematics)[0]
        q = 1 + self.regularize * (self.bulk - 1)

        self.results.elasticity = [
            self.umat.hessian(self.results.kinematics, *args, **kwargs)[0]
            + self.results.p / q * d2JdFdF
        ]

        return self.results.elasticity

    def _kirchhoff_stress(self, field=None):

        self._gradient(field)

        P = self.results.stress[0]
        F = self.results.kinematics[0]

        return dot(P, transpose(F))

    def _cauchy_stress(self, field=None):

        self._gradient(field)

        P = self.results.stress[0]
        F = self.results.kinematics[0]
        J = det(F)

        return dot(P, transpose(F)) / J

[adtzlr]

However, there seems to be some problems with convergence and wrong results for plane strain.

[adtzlr]

Let's save resources and stop the work on this unimplemented experiment.

[adtzlr]

Finally, I found a solution. The basic idea evolved out of #296. Usage will be basically

# create a hexahedron-region on a cube
mesh = fem.Cube(n=6)
region = fem.RegionHexahedron(mesh)

# add a field container (with a displacement field)
field = fem.FieldsMixed(region, n=1)

# define the constitutive material behaviour (distortional part only!)
# and create a solid body (define the bulk modulus here)
umat = fem.NeoHooke(mu=1)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)

# apply a uniaxial elongation on the cube
boundaries, lc = fem.dof.uniaxial(field, move=-0.1, clamped=True)
res = fem.newtonrhapson(items=[solid], **lc)