adtzlr · adtzlr · Feb 20, 2024 · Feb 20, 2024 · Feb 20, 2024
diff --git a/src/felupe/assembly/_integral.py b/src/felupe/assembly/_integral.py
@@ -21,6 +21,7 @@
 
 from ..field._axi import FieldAxisymmetric
 from ..field._base import Field
+from ..field._dual import FieldDual
 from ..field._planestrain import FieldPlaneStrain
 from ._axi import IntegralFormAxisymmetric
 from ._cartesian import IntegralFormCartesian
@@ -207,6 +208,7 @@ def __init__(self, fun, v, dV, u=None, grad_v=None, grad_u=None):
 
         IntForm = {
             Field: IntegralFormCartesian,
+            FieldDual: IntegralFormCartesian,
             FieldPlaneStrain: IntegralFormCartesian,
             FieldAxisymmetric: IntegralFormAxisymmetric,
         }[type(self.v[0])]

diff --git a/src/felupe/field/_dual.py b/src/felupe/field/_dual.py
@@ -56,6 +56,11 @@ class FieldDual(Field):
     mesh: Mesh or None, optional
         A mesh which is used for the dual region (default is None). If None, the mesh
         is taken from the region.
+    disconnect : bool, optional
+        A flag to disconnect the dual mesh (each cell has its own points). Default is
+        None. If None, a disconnected mesh is used except for regions with
+        quadratic-triangle and quadratic-tetra elements as well as for regions with the
+        MINI-elements.
     **kwargs : dict, optional
         Optional keyword arguments for the dual region.
 
@@ -92,6 +97,7 @@ def __init__(
         offset=0,
         npoints=None,
         mesh=None,
+        disconnect=None,
         **kwargs,
     ):
         # create dual regions
@@ -120,7 +126,9 @@ def __init__(
             RegionTetraMINI: {"disconnect": False},
             RegionTriangleMINI: {"disconnect": False},
             RegionLagrange: {},
-        }
+        }[type(region)]
+        if disconnect is not None:
+            mesh_kwargs["disconnect"] = disconnect
         points_per_cell = {
             RegionConstantHexahedron: 1,
             RegionConstantQuad: 1,
@@ -145,7 +153,7 @@ def __init__(
                 points_per_cell=points_per_cell[RegionDual],
                 offset=offset,
                 npoints=npoints,
-                **mesh_kwargs[type(region)],
+                **mesh_kwargs,
             )
 
         region_dual = RegionDual(mesh, **{**kwargs0, **kwargs})

diff --git a/src/felupe/math/_tensor.py b/src/felupe/math/_tensor.py
@@ -223,28 +223,24 @@ def dot(A, B, mode=(2, 2), parallel=False, **kwargs):
 
     if mode == (2, 2):
         return einsum("ik...,kj...->ij...", A, B, **kwargs)
-
     elif mode == (1, 1):
         return einsum("i...,i...->...", A, B, **kwargs)
-
     elif mode == (4, 4):
         return einsum("ijkp...,plmn...->ijklmn...", A, B, **kwargs)
-
     elif mode == (2, 1):
         return einsum("ij...,j...->i...", A, B, **kwargs)
-
     elif mode == (1, 2):
         return einsum("i...,ij...->j...", A, B, **kwargs)
-
+    elif mode == (2, 3):
+        return einsum("im...,mjk...->ijk...", A, B, **kwargs)
+    elif mode == (3, 2):
+        return einsum("ijm...,mk...->ijk...", A, B, **kwargs)
     elif mode == (4, 1):
         return einsum("ijkl...,l...->ijk...", A, B, **kwargs)
-
     elif mode == (1, 4):
         return einsum("i...,ijkl...->jkl...", A, B, **kwargs)
-
     elif mode == (2, 4):
         return einsum("im...,mjkl...->ijkl...", A, B, **kwargs)
-
     elif mode == (4, 2):
         return einsum("ijkm...,ml...->ijkl...", A, B, **kwargs)
 
@@ -262,6 +258,10 @@ def ddot(A, B, mode=(2, 2), parallel=False, **kwargs):
 
     if mode == (2, 2):
         return einsum("ij...,ij...->...", A, B, **kwargs)
+    if mode == (2, 3):
+        return einsum("ij...,ijk...->k...", A, B, **kwargs)
+    if mode == (3, 2):
+        return einsum("ijk...,ij...->k...", A, B, **kwargs)
     elif mode == (2, 4):
         return einsum("ij...,ijkl...->kl...", A, B, **kwargs)
     elif mode == (4, 2):

diff --git a/src/felupe/mechanics/__init__.py b/src/felupe/mechanics/__init__.py
@@ -1,12 +1,13 @@
 from ._curve import CharacteristicCurve
-from ._helpers import StateNearlyIncompressible
+from ._helpers import StateNearlyIncompressible, StateNearlyIncompressibleX
 from ._item import FormItem
 from ._job import Job
 from ._multipoint import MultiPointConstraint, MultiPointContact
 from ._pointload import PointLoad
 from ._solidbody import SolidBody
 from ._solidbody_gravity import SolidBodyGravity
 from ._solidbody_incompressible import SolidBodyNearlyIncompressible
+from ._solidbody_incompressible_x import SolidBodyNearlyIncompressibleX
 from ._solidbody_pressure import SolidBodyPressure
 from ._step import Step
 
@@ -19,6 +20,7 @@
     "SolidBody",
     "SolidBodyGravity",
     "SolidBodyNearlyIncompressible",
+    "SolidBodyNearlyIncompressibleX",
     "SolidBodyPressure",
     "Step",
     "MultiPointConstraint",

diff --git a/src/felupe/mechanics/_helpers.py b/src/felupe/mechanics/_helpers.py
@@ -17,9 +17,9 @@
 """
 import numpy as np
 
-from ..assembly import IntegralForm
+from ..assembly import IntegralForm, IntegralFormCartesian
 from ..constitution import AreaChange
-from ..field import FieldAxisymmetric
+from ..field import FieldAxisymmetric, FieldDual
 from ..math import det
 
 
@@ -113,3 +113,94 @@ def v(self):
             dA = self.field.region.dV
             dV = 2 * np.pi * R * dA
         return (det(self.F[0]) * dV).sum(0)
+
+
+class StateNearlyIncompressibleX:
+    "A State with internal fields for (nearly) incompressible solid bodies."
+
+    def __init__(self, field, pressure=None, volume_ratio=None, **kwargs):
+        self.field = field
+        self.u = self.field[0].values
+
+        self.pressure = pressure
+        self.volume_ratio = volume_ratio
+
+        if self.pressure is None:
+            self.pressure = FieldDual(field.region, **kwargs)
+        if self.volume_ratio is None:
+            self.volume_ratio = FieldDual(field.region, values=1, **kwargs)
+
+        self.dJdF = AreaChange().function
+
+        # deformation gradient
+        self.F = field.extract()
+
+        # inverse of volume matrix
+        self.inv_V = np.linalg.inv(self.volume().T).T
+
+    def integrate_shape_function_gradient(self, parallel=False):
+        r"""Return the Integrated shape function gradient matrix w.r.t. the deformed
+        coordinates.
+
+        Notes
+        -----
+        ..  math::
+
+            h = \int_V \delta \boldsymbol{F} : \frac{\partial J}{\partial\boldsymbol{F}}
+                \ \Delta p ~ dV
+
+        """
+        return IntegralForm(
+            fun=self.dJdF(self.F),
+            v=self.field,
+            u=self.pressure & None,
+            dV=self.field.region.dV,
+            grad_v=[True],
+            grad_u=[False],
+        ).integrate(parallel=parallel)[0]
+
+    def volume(self, parallel=False):
+        r"""Return integrated differential (undeformed) volumes matrix with dual-trial
+        and dual-test fields.
+
+        Notes
+        -----
+        ..  math::
+
+            V = \int_V \delta p \Delta p ~ dV
+
+        """
+        return IntegralForm(
+            fun=[np.ones((1, 1))],
+            v=self.pressure & None,
+            dV=self.field.region.dV,
+            u=self.pressure & None,
+            grad_v=[False],
+            grad_u=[False],
+        ).integrate(parallel=parallel)[0]
+
+    def fp(self, parallel=False):
+        dV = self.field.region.dV
+        J = self.volume_ratio.interpolate()
+        v = self.pressure & None
+        return IntegralForm(
+            fun=[det(self.F[0]) - J], v=v, dV=dV, grad_v=[False]
+        ).integrate(parallel=parallel)[0]
+
+    def fJ(self, bulk, parallel=False):
+        dV = self.field.region.dV
+        p = self.pressure.interpolate()
+        J = self.volume_ratio.interpolate()
+        v = self.pressure & None
+        return IntegralForm(
+            fun=[bulk * (J - 1) - p], v=v, dV=dV, grad_v=[False]
+        ).integrate(parallel=parallel)[0]
+
+    def constraint(self, bulk, parallel=False):
+        dV = self.field.region.dV
+        p = self.pressure.interpolate()
+        detF = det(self.F[0])
+        v = self.pressure & None
+        return IntegralForm(
+            fun=[bulk * (detF - 1) - p], v=v, dV=dV, grad_v=[False]
+        ).integrate(parallel=parallel)[0]
diff --git a/src/felupe/mechanics/_solidbody_incompressible.py b/src/felupe/mechanics/_solidbody_incompressible.py
@@ -281,7 +281,6 @@ def _vector(self, field=None, parallel=False, items=None, args=(), kwargs={}):
         self.results.stress = self._gradient(
             field, parallel=parallel, args=args, kwargs=kwargs
         )
-
         form = self._form(
             fun=self.results.stress,
             v=self.field,
@@ -294,7 +293,7 @@ def _vector(self, field=None, parallel=False, items=None, args=(), kwargs={}):
 
         values = form.integrate(parallel=parallel)
         np.add(values[0], h * (self.bulk * (v / self.V - 1) - p), out=values[0])
-
+        
         self.results.force = form.assemble(values=values)
 
         return self.results.force
@@ -333,7 +332,7 @@ def _extract(self, field, parallel=False):
         # change of state variables due to change of displacement field
         dJ = ddot(h, du, mode=(2, 2)) / self.V + (v / self.V - J)
         dp = self.bulk * (dJ + J - 1) - p
-
+        
         self.field = field
         self.results.kinematics = self.results.state.F = self.field.extract(
             out=self.results.kinematics
@@ -366,7 +365,7 @@ def _gradient(self, field=None, parallel=False, args=(), kwargs={}):
     def _hessian(self, field=None, parallel=False, args=(), kwargs={}):
         if field is not None:
             self.results.kinematics = self._extract(field, parallel=parallel)
-
+        
         d2JdF2 = self._area_change.gradient
         F = self.results.kinematics[0]
         statevars = self.results.statevars