Ideas for `SolidBodyContact`

[adtzlr]

This issue serves as a collection of ideas on how to implement a contact between two solid bodies (once more).

Requirements

• should be treated as an item, like SolidBody with methods to assemble the force vector- and the stiffness matrix

Roadmap

• Start with an elastic support - like trusses, but in a weak formulation and with a stiffness value per unit (which? deformed?) area and a fixed direction.

[adtzlr]

Elastic support

A part of the boundary $\Gamma$ of a solid body is connected to an (uniaxial) elastic support. The stiffness value is defined as force per deformed area. The direction of the springs are fixed, i.e. they do not vary w.r.t. the deformation of the solid body, no matter how distorted the solid body will be.

$$ \delta \Pi_\text{ext} = \int_{\partial v} \delta \boldsymbol{u} \cdot K \ (\boldsymbol{u} \cdot \boldsymbol{N}) \ \boldsymbol{n} \ da $$

Re-formulated by integration over the undeformed area:

$$ \delta \Pi_\text{ext} = \int_{\partial V} \delta \boldsymbol{u} \cdot K \ (\boldsymbol{u} \cdot \boldsymbol{N}) \ J \boldsymbol{F}^{-T} \boldsymbol{N} \ dA $$

This weak form requires both the deformation gradient and the displacement values. What is the most efficient way to implement this (without felupe.Form)? A mixed field is not directly useable, because the sub-block matrices / vectors must be summed instead of stacked. Just use one field, where a) the gradient and b) the values are extracted?

• A mixed-field with (u, F). No, just use one field with two extracts.
• Automatic differentiation? Yes.

[adtzlr]

Here's a first example.

code

import felupe as fem
import numpy as np
import pypardiso

mesh = fem.Cube(n=11)
region = fem.RegionHexahedron(mesh)
edge = fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 0)

field = fem.Field(region, dim=3).as_container()
displacement = fem.Field(edge, dim=3).as_container()


class SolidBodyElasticSupport:
    def __init__(self, field, stiffness=1.0, direction=[1.0, 0.0, 0.0]):
        self.field = field
        self.stiffness = stiffness
        u = self.field.extract(grad=False)[0]
        direction = np.array(direction) / np.linalg.norm(direction, axis=0)
        self.direction = np.broadcast_to(
            direction[..., np.newaxis, np.newaxis], u.shape
        )
        self.assemble = fem.mechanics.Assemble(
            self.vector, self.matrix, multiplier=-1.0
        )
        self.results = fem.mechanics.Results()

    def fun(self, F, u, K, M, N):
        import tensortrax as tr
        import tensortrax.math as tm

        x = F if isinstance(F, tr.Tensor) else u
        M = tm.array(M, shape=(3,), like=x)
        N = tm.array(N, shape=(3,), like=x)

        # n da = J F^-T N dA
        t = tm.linalg.det(F) * tm.transpose(tm.linalg.inv(F)) @ N
        sign = tm.sign(M @ t)
        return -K * u @ M * sign * t

    def extract(self, field=None):
        if field is not None:
            self.field = field
        F = self.field.extract()[0]
        u = self.field.extract(grad=False)[0]
        return F, u

    def vector(self, field=None, **kwargs):
        import tensortrax as tr

        F, u = self.extract(field=field)
        K, N, M = self.stiffness, self.field.region.normals, self.direction
        vector = fem.IntegralForm(
            fun=[tr.function(self.fun, ntrax=2)(F, u, K, M, N)],
            v=self.field,
            dV=self.field.region.dV,
            grad_v=[False],
        ).assemble()
        return vector

    def matrix(self, field=None, **kwargs):
        import tensortrax as tr

        F, u = self.extract(field=field)
        K, N, M = self.stiffness, self.field.region.normals, self.direction
        matrix = (
            fem.IntegralForm(
                fun=[tr.jacobian(self.fun, ntrax=2, wrt=1)(F, u, K, M, N)],
                v=self.field,
                dV=self.field.region.dV,
                u=self.field,
                grad_v=[False],
                grad_u=[False],
            ).assemble()
            + fem.IntegralForm(
                fun=[tr.jacobian(self.fun, ntrax=2, wrt=0)(F, u, K, M, N)],
                v=self.field,
                dV=self.field.region.dV,
                u=self.field,
                grad_v=[False],
                grad_u=[True],
            ).assemble()
        )
        return matrix


boundaries = fem.BoundaryDict(move=fem.Boundary(field[0], fx=1))
umat = fem.NeoHooke(mu=1)

solid = fem.SolidBodyNearlyIncompressible(umat=umat, field=field, bulk=5000)
support = SolidBodyElasticSupport(field=displacement, stiffness=10.0)

move = fem.math.linsteps([0, -0.4], num=5, axis=0, axes=3)
ramp = {boundaries["move"]: move}
step = fem.Step(items=[solid, support], ramp=ramp, boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate(solver=pypardiso.spsolve)

ax = solid.imshow("Principal Values of Cauchy Stress", project=fem.topoints)

[adtzlr]

Next improvement. Quite a lot of overhead because each and every Gauss-point on the active edges are tested against each other. Also, the contact stiffness seems to depend on the number of cells (Gauss-points).

code

import felupe as fem
import numpy as np
import pypardiso

mesh = fem.MeshContainer(
    [fem.Cube(n=3), fem.Cube(n=4).translate(1.05, 0)], merge=True
).stack()
region = fem.RegionHexahedron(mesh)
edges = [
    fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1),
    fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1.05),
]
field = fem.Field(region, dim=3).as_container()


class SolidBodyContact:
    def __init__(self, region, other_region, stiffness=1.0):
        cells = np.tile(region.mesh.cells, (other_region.mesh.ncells, 1))
        other_cells = np.tile(other_region.mesh.cells, (1, region.mesh.ncells))

        self.region = region.copy(mesh=region.mesh.copy(cells=cells))
        self.other_region = region.copy(
            mesh=region.mesh.copy(cells=other_cells.reshape(cells.shape))
        )

        (
            self.region.dA,
            self.region.dV,
            self.region.normals,
            self.region.tangents,
        ) = self.region._init_faces()
        (
            self.other_region.dA,
            self.other_region.dV,
            self.other_region.normals,
            self.other_region.tangents,
        ) = self.other_region._init_faces()

        self.field = fem.Field(self.region, dim=3).as_container()
        self.other_field = fem.Field(self.other_region, dim=3).as_container()

        self.stiffness = stiffness
        self.assemble = fem.mechanics.Assemble(
            self.vector, self.matrix, multiplier=-1.0
        )
        self.results = fem.mechanics.Results()

    def fun(self, F, u, u_ext, K, N):
        import tensortrax as tr
        import tensortrax.math as tm

        x = F if isinstance(F, tr.Tensor) else u if isinstance(u, tr.Tensor) else u_ext
        N = tm.array(N, shape=(3,), like=x)
        t = tm.linalg.det(F) * tm.transpose(tm.linalg.inv(F)) @ N
        n = t / tm.sqrt(tm.sum(t**2, axis=0))
        un = -(u_ext - u) @ n
        gn = tm.if_else(un.x > 0, un, 0 * un)
        return -K * gn * t

    def extract(self, field=None):
        if field is not None:
            self.field.link(field)
            self.other_field.link(field)

        F = self.field.extract()[0]
        u = self.field.extract(grad=False)[0]
        u_ext = self.other_field.extract(grad=False)[0]
        return F, u, u_ext

    def vector(self, field=None, **kwargs):
        import tensortrax as tr

        F, u, u_ext = self.extract(field=field)
        K, N = self.stiffness, self.field.region.normals
        vector = fem.IntegralForm(
            fun=[tr.function(self.fun, ntrax=2)(F, u, u_ext, K, N)],
            v=self.field,
            dV=self.field.region.dV,
            grad_v=[False],
        ).assemble()
        return vector

    def matrix(self, field=None, **kwargs):
        import tensortrax as tr

        F, u, u_ext = self.extract(field=field)
        K, N = self.stiffness, self.field.region.normals
        matrix = (
            fem.IntegralForm(
                fun=[tr.jacobian(self.fun, ntrax=2, wrt=1)(F, u, u_ext, K, N)],
                v=self.field,
                dV=self.field.region.dV,
                u=self.field,
                grad_v=[False],
                grad_u=[False],
            ).assemble()
            + fem.IntegralForm(
                fun=[tr.jacobian(self.fun, ntrax=2, wrt=2)(F, u, u_ext, K, N)],
                v=self.field,
                dV=self.field.region.dV,
                u=self.other_field,
                grad_v=[False],
                grad_u=[False],
            ).assemble()
            + fem.IntegralForm(
                fun=[tr.jacobian(self.fun, ntrax=2, wrt=0)(F, u, u_ext, K, N)],
                v=self.field,
                dV=self.field.region.dV,
                u=self.field,
                grad_v=[False],
                grad_u=[True],
            ).assemble()
        )
        return matrix


boundaries = fem.BoundaryDict(
    move=fem.Boundary(field[0], fx=2.05), fixed=fem.Boundary(field[0], fx=0)
)
umat = fem.NeoHooke(mu=1)

solid = fem.SolidBodyNearlyIncompressible(umat=umat, field=field, bulk=5000)
support = SolidBodyContact(*edges[::], stiffness=1.5)
support2 = SolidBodyContact(*edges[::-1], stiffness=1.5)

move = fem.math.linsteps([0, -0.3], num=10, axis=0, axes=3)
ramp = {boundaries["move"]: move}
step = fem.Step(items=[solid, support, support2], ramp=ramp, boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate(solver=pypardiso.spsolve)

solid.plot("Displacement", opacity=0.9).show()

[adtzlr]

see also adtzlr/matadi#143

import felupe as fem
import numpy as np

mesh = fem.Cube(n=3)
region = fem.RegionHexahedron(mesh)
field = fem.FieldsMixed(region, n=1)
u = fem.Field(fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1), dim=3)
fieldb = u & u
fieldb.link(field)

import matadi as mat
from matadi.math import det, inv, gradient, trace, zeros, eye, sqrt, sum1, dot, DM

# def fun(x):
#     F, u, statevars = x[0], x[1], x[-1]
#     return [det(F) * u[0], statevars]

# x = [mat.Variable("F", 3, 3), mat.Variable("x", 3), mat.Variable("z", 1)]
# umat = mat.Material(x, fun, statevars=1)

# statevars = field.extract()[-1][0]
# dWdF, statevars_new = umat.gradient([*field.extract(), statevars])
# dWdF, dWdu, statevars = umat.hessian([*field.extract(), statevars])


# def fun(x):
#     F, p, J, statevars = x[0], x[1], x[2], x[-1]
#     W = (
#         mat.models.neo_hooke(F * (J / det(F))**(1 / 3), C10=0.5)
#         + p * (det(F) - J)
#         + mat.models.volumetric(J, bulk=5000)
#     )
#     return gradient(W, F), gradient(W, p), gradient(W, J), statevars

# x = [mat.Variable("F", 3, 3), mat.Variable("p"), mat.Variable("J"), mat.Variable("z")]

def fun(x):
    F, u, N = x[0], x[1], x[-1]
    p = -1.0
    return zeros(3, 3), -p * det(F) * inv(F).T @ N, N


x = [mat.Variable("F", 3, 3), mat.Variable("u", 3), mat.Variable("N", 3)]

umatb = mat.MaterialTensor(x, fun, statevars=1, triu=True)

boundaries = dict(fixed=fem.Boundary(field[0], fx=0))
umat = fem.NeoHooke(mu=1, bulk=2)
solid = fem.SolidBody(umat=umat, field=field)


def apply(x):
    if len(x) == 2:
        return sum(x)
    else:
        return x[0] + x[1].T + x[2]


solidb = fem.SolidBody(
    umat=umatb,
    field=fieldb,
    block=False,
    apply=apply,
    statevars=u.region.normals,
    multiplier=-1.0,
)
v = u.as_container()
v.link(field)
solidc = fem.SolidBodyPressure(v, pressure=-1.0)

step = fem.Step(items=[solid, solidb], boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate(verbose=2, tol=1e3)

ax = solid.imshow("Principal Values of Cauchy Stress", factor=1)

[adtzlr]

Even with tensortrax, there is something strange here.

import felupe as fem
import numpy as np

mesh = fem.Cube(n=3)
region = fem.RegionHexahedron(mesh)
field = fem.FieldsMixed(region, n=1)
u = fem.Field(fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1), dim=3)
fieldb = u & u
fieldb.link(field)

import tensortrax.math as tm
import tensortrax as tr

N = u.region.normals


def stress(x, p, N):
    F, u = x[:2]

    def fun(F, u, p, N):
        J = tm.linalg.det(F)
        invF = tm.linalg.inv(F)
        
        return -p * J * invF.T @ tm.array(N, like=F[0])

    dWdu = tr.function(fun, wrt=0, ntrax=2)(F, u, p=p, N=N)
    return 0 * F, dWdu, x[-1]


def elasticity(x, p, N):
    F, u = x[:2]

    def fun(F, u, p, N):
        J = tm.linalg.det(F)
        invF = tm.linalg.inv(F)
        return -p * J * invF.T @ tm.array(N, like=F[0])

    d2WdudF = tr.jacobian(fun, wrt=0, ntrax=2)(F, u, p=p, N=N)
    return np.zeros((3, 3, 3, 3, 4, 4)), d2WdudF, np.zeros((3, 3, 4, 4))


umatb = fem.Material(stress, elasticity, p=-1.0, N=u.region.normals)

boundaries = dict(fixed=fem.Boundary(field[0], fx=0))
umat = fem.NeoHooke(mu=1, bulk=2)
solid = fem.SolidBody(umat=umat, field=field)


solidb = fem.SolidBody(
    umat=umatb,
    field=fieldb,
    block=False,
    apply=sum,
    multiplier=-1.0,
)
v = u.as_container()
v.link(field)
solidc = fem.SolidBodyPressure(v, pressure=-1.0)

step = fem.Step(items=[solid, solidb], boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate(verbose=2, tol=1e3)

ax = solid.imshow("Principal Values of Cauchy Stress", factor=1)

[adtzlr]

import felupe as fem
import numpy as np
import tensortrax.math as tm
import tensortrax as tr
import pypardiso

mesh = fem.MeshContainer(
    [fem.Cube(n=3), fem.Cube(n=4).translate(1.05, 0)], merge=True
).stack()
region = fem.RegionHexahedron(mesh)
edges = [
    fem.Field(fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1), dim=3),
    fem.Field(fem.RegionHexahedronBoundary(mesh, mask=mesh.x == 1.05), dim=3),
]
field = fem.Field(region, dim=3).as_container()


def fun(F, u, N):
    J = tm.linalg.det(F)
    invF = tm.linalg.inv(F)
    p = 5000 * (J - 1)
    return -p * J * invF.T @ tm.array(N, like=F[0])


def stress(x, N):
    F, u = x[:2]
    dWdu = tr.function(fun, wrt=0, ntrax=2)(F, u, N=N)
    return None, dWdu, x[-1]


def elasticity(x, N):
    F, u = x[:2]
    d2WdudF = tr.jacobian(fun, wrt=0, ntrax=2)(F, u, N=N)
    return None, None, d2WdudF, None

boundaries = fem.BoundaryDict(
    move=fem.Boundary(field[0], fx=2.05), fixed=fem.Boundary(field[0], fx=0)
)
umat = fem.NeoHooke(mu=1)

solid = fem.SolidBodyNearlyIncompressible(umat=umat, field=field, bulk=5000)
contact_1 = fem.SolidBody(
    umat=fem.Material(stress, elasticity, N=edges[0].region.normals),
    field=edges[0] & edges[0],
    block=False,
    apply=sum,
    multiplier=-1.0,
)
contact_2 = fem.SolidBody(
    umat=fem.Material(stress, elasticity, N=edges[1].region.normals),
    field=edges[1] & edges[1],
    block=False,
    apply=sum,
    multiplier=-1.0,
)

move = fem.math.linsteps([0, -0.3], num=10, axis=0, axes=3)
ramp = {boundaries["move"]: move}
step = fem.Step(items=[solid, contact_1, contact_2], ramp=ramp, boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate(solver=pypardiso.spsolve)

solid.plot("Displacement", opacity=0.9).show()

[adtzlr]

A contact between two solid bodies should be defined in the most high-level form as

contact = SolidBodyContact(solid, other_solid)

The displacement field of the solid bodies should be moved to the (free) boundaries to obtain and assemble the required weak-forms.

$$ \Pi = \int_V \psi\ dV + \int_{\partial V} \lambda\ g\ dA + \frac{1}{2} \int_{\partial V} \epsilon\ g^2\ dA $$

$g$ is the amount of the gap, i.e. the norm of the gap vector.

$$ \Pi = \int_V \psi\ dV + \lambda \left( \int_{\partial V_1} g_1\ dA_1 + \int_{\partial V_2} g_2\ dA_2 \right)$$

and

$$ g_1 = n_1 \cdot (x_2 - x_1)$$

For this kind of formulation, equal quadrature schemes and element formulations are required for the contact pairs.