Merge pull request #260 from adtzlr/fix-planestrain-field-integration

Field Integration: Fix field-trimming for plane strain

felupe/_assembly/_base.py

@@ -163,8 +163,15 @@ def integrate(self, parallel=False, jit=False):
         dV = self.dV
         fun = self.fun
 
-        if len(fun) > v.dim:
-            fun = fun[tuple([slice(v.dim)] * (len(fun.shape) - 2))]
+        # plane strain
+        # trim 3d vector-valued functions to the dimension of the field
+        function_dimension = len(fun.shape) - 2
+        function_is_vector = function_dimension >= 1
+        function_is_3d = len(fun) == 3
+        field_is_2d = v.dim == 2
+
+        if function_is_vector and function_is_3d and field_is_2d:
+            fun = fun[tuple([slice(2)] * function_dimension)]
 
         if parallel:
             einsum = einsumt

felupe/mechanics/_solidbody_tensor.py

@@ -77,7 +77,6 @@ def _vector(
             self.field = field
 
         self.results.stress = self._gradient(field, args=args, kwargs=kwargs)
-
         self.results.force = self._form(
             fun=self.results.stress[slice(items)],
             v=self.field,
@@ -141,7 +140,7 @@ def _kirchhoff_stress(self, field=None):
 
         self._gradient(field)
 
-        P = self.results.stress[0][0]
+        P = self.results.stress[0]
         F = self.results.kinematics[0]
 
         return dot(P, transpose(F))
@@ -150,7 +149,7 @@ def _cauchy_stress(self, field=None):
 
         self._gradient(field)
 
-        P = self.results.stress[0][0]
+        P = self.results.stress[0]
         F = self.results.kinematics[0]
         J = det(F)
 

tests/test_mechanics_tensor.py

@@ -32,23 +32,43 @@
 
 def pre_umat():
 
-    LE = fe.LinearElastic(E=210000, nu=0.3)
+    NH = fe.NeoHooke(mu=1, bulk=5000)
 
-    class LETensor:
-        def __init__(self, LE):
-            self.LE = LE
-            self.x = [np.zeros((3, 3)), np.zeros((5, 16))]
+    class NHTensor:
+        def __init__(self, NH):
+            self.NH = NH
+            self.x = [np.zeros((3, 3)), np.zeros((1, 1))]
 
         def function(self, x):
-            F, statevars = x
+            F, statevars = x[:-1], x[-1]
             # return dummy state variables along with stress
-            return self.LE.stress(F), statevars
+            return [*self.NH.gradient(F), statevars]
 
         def gradient(self, x):
-            F, statevars = x
-            return self.LE.elasticity(F)
+            F, statevars = x[:-1], x[-1]
+            return self.NH.hessian(F)
 
-    return LETensor(LE)
+    return NHTensor(NH)
+
+def pre_umat_mixed():
+
+    NH = fe.ThreeFieldVariation(fe.NeoHooke(mu=1, bulk=5000))
+
+    class NHTensor:
+        def __init__(self, NH):
+            self.NH = NH
+            self.x = [np.zeros((3, 3)), np.zeros((1, 1)), np.zeros((1, 1)), np.zeros((1, 1))]
+
+        def function(self, x):
+            F, statevars = x[:-1], x[-1]
+            # return dummy state variables along with stress
+            return [*self.NH.gradient(F), statevars]
+
+        def gradient(self, x):
+            F, statevars = x[:-1], x[-1]
+            return self.NH.hessian(F)
+
+    return NHTensor(NH)
 
 
 def test_solidbody_tensor():
@@ -57,15 +77,15 @@ def test_solidbody_tensor():
 
     m = fe.Cube(n=3)
     r = fe.RegionHexahedron(m)
-    u = fe.Field(r, dim=3)
-    v = fe.FieldContainer([u])
-
-    sv = np.zeros((5, 16, r.quadrature.npoints, m.ncells))
+    v = fe.FieldsMixed(r, n=1)
+
+    sv = np.zeros((1, 1, r.quadrature.npoints, m.ncells))
 
     for statevars in [sv, None]:
+
         b = fe.SolidBodyTensor(umat, v, statevars)
         r = b.assemble.vector()
-
+    
         K = b.assemble.matrix(v)
         K = b.assemble.matrix()
         r = b.assemble.vector(v)
@@ -76,16 +96,53 @@ def test_solidbody_tensor():
         t = b.evaluate.kirchhoff_stress()
         C = b.results.elasticity
         z = b.results.statevars
-
+    
         assert K.shape == (81, 81)
         assert r.shape == (81, 1)
         assert F[0].shape == (3, 3, 8, 8)
-        assert P[0][0].shape == (3, 3, 8, 8)
+        assert P[0].shape == (3, 3, 8, 8)
+        assert s.shape == (3, 3, 8, 8)
+        assert t.shape == (3, 3, 8, 8)
+        assert C[0].shape == (3, 3, 3, 3, 8, 8)
+        assert z.shape == (1, 1, 8, 8)
+
+
+def test_solidbody_tensor_mixed():
+
+    umat = pre_umat_mixed()
+
+    m = fe.Cube(n=3)
+    r = fe.RegionHexahedron(m)
+    v = fe.FieldsMixed(r, n=3)
+
+    sv = np.zeros((1, 1, r.quadrature.npoints, m.ncells))
+
+    for statevars in [sv, None]:
+
+        b = fe.SolidBodyTensor(umat, v, statevars)
+        r = b.assemble.vector()
+
+        K = b.assemble.matrix(v)
+        K = b.assemble.matrix()
+        r = b.assemble.vector(v)
+        r = b.assemble.vector()
+        F = b.results.kinematics
+        P = b.results.stress
+        s = b.evaluate.cauchy_stress()
+        t = b.evaluate.kirchhoff_stress()
+        C = b.results.elasticity
+        z = b.results.statevars
+
+        assert K.shape == (97, 97)
+        assert r.shape == (97, 1)
+        assert F[0].shape == (3, 3, 8, 8)
+        assert P[0].shape == (3, 3, 8, 8)
         assert s.shape == (3, 3, 8, 8)
         assert t.shape == (3, 3, 8, 8)
         assert C[0].shape == (3, 3, 3, 3, 8, 8)
-        assert z.shape == (5, 16, 8, 8)
+        assert z.shape == (1, 1, 8, 8)
 
 
 if __name__ == "__main__":
     test_solidbody_tensor()
+    test_solidbody_tensor_mixed()