Merge pull request #224 from adtzlr/add-solidbody-tensor

Add SolidBodyTensor

README.md

@@ -50,9 +50,11 @@ All notable changes to this project will be documented in this file. The format
 
 ### Added
 - Add `SolidBody.evaluate.kirchhoff_stress()` method. Contrary to the Cauchy stress method, this gives correct results in incompressible plane stress.
+- Add `SolidBodyTensor` for tensor-based material definitions with state variables.
 
 ### Fixed
 - Fix `tovoigt()` helper for data with more or less than two trailing axes and 2D tensors.
+- Fix errors for `force()` and `moment()` helpers if the residuals are sparse.
 
 ## [3.1.0] - 2022-05-02
 

felupe/__init__.py

@@ -111,6 +111,7 @@
 from .mechanics import (
     SolidBody,
     SolidBodyPressure,
+    SolidBodyTensor,
 )
 
 __all__ = [

felupe/mechanics/__init__.py

@@ -1,2 +1,3 @@
 from ._solidbody import SolidBody
 from ._solidbody_pressure import SolidBodyPressure
+from ._solidbody_tensor import SolidBodyTensor

felupe/mechanics/_solidbody_tensor.py

@@ -0,0 +1,166 @@
+# -*- coding: utf-8 -*-
+"""
+ _______  _______  ___      __   __  _______  _______ 
+|       ||       ||   |    |  | |  ||       ||       |
+|    ___||    ___||   |    |  | |  ||    _  ||    ___|
+|   |___ |   |___ |   |    |  |_|  ||   |_| ||   |___ 
+|    ___||    ___||   |___ |       ||    ___||    ___|
+|   |    |   |___ |       ||       ||   |    |   |___ 
+|___|    |_______||_______||_______||___|    |_______|
+
+This file is part of felupe.
+
+Felupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+Felupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Felupe.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+from .._field import Field, FieldMixed, FieldsMixed, FieldAxisymmetric
+from .._assembly import IntegralForm, IntegralFormMixed, IntegralFormAxisymmetric
+from ..constitution import AreaChange
+from ..math import inv, dot, transpose, det
+from ._helpers import Assemble, Evaluate, Results
+
+
+class SolidBodyTensor:
+    """A SolidBody with methods for the assembly of sparse vectors/matrices
+    based on a ``MaterialTensor`` with state variables."""
+
+    def __init__(self, umat, field, statevars=None):
+
+        self.umat = umat
+        self.field = field
+
+        if isinstance(field, FieldMixed):
+            self._dV = self.field[0].region.dV
+        else:
+            self._dV = self.field.region.dV
+
+        self.results = Results(stress=True, elasticity=True)
+        self.results.kinematics = self._extract(self.field)
+        self.results.statevars = statevars
+
+        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._area_change = AreaChange()
+
+        self._form = {
+            Field: IntegralForm,
+            FieldMixed: IntegralFormMixed,
+            FieldsMixed: IntegralFormMixed,
+            FieldAxisymmetric: IntegralFormAxisymmetric,
+        }[type(self.field)]
+
+        self._kwargs = {
+            Field: dict(dV=self._dV, grad_v=True, grad_u=True),
+            FieldMixed: dict(dV=self._dV),
+            FieldsMixed: dict(dV=self._dV),
+            FieldAxisymmetric: dict(dV=self._dV, grad_v=True, grad_u=True),
+        }[type(self.field)]
+
+    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[slice(items)],
+            v=self.field,
+            **self._kwargs,
+        ).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)
+
+        self.results.stiffness = self._form(
+            fun=self.results.elasticity[slice(items)],
+            v=self.field,
+            u=self.field,
+            **self._kwargs,
+        ).assemble(parallel=parallel, jit=jit)
+
+        return self.results.stiffness
+
+    def _extract(self, field):
+
+        self.field = field
+
+        self.results.kinematics = self.field.extract()
+        if isinstance(self.field, Field):
+            self.results.kinematics = (self.results.kinematics,)
+
+        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)
+
+        function = self.umat.function(
+            *self.results.kinematics, self.results.statevars, *args, **kwargs
+        )
+
+        self.results.stress, self.results.statevars = function[:-1], function[-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.gradient(
+            *self.results.kinematics, self.results.statevars, *args, **kwargs
+        )
+
+        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

felupe/tools/_post.py

@@ -27,9 +27,12 @@
 
 import numpy as np
 from scipy.interpolate import interp1d
+from scipy.sparse import issparse
 
 
 def force(field, r, boundary, offsets=None):
+    if issparse(r):
+        r = r.toarray()
     if offsets is None:
         return ((r.reshape(-1, field.dim))[boundary.points]).sum(0)
     else:
@@ -39,7 +42,8 @@ def force(field, r, boundary, offsets=None):
 
 
 def moment(field, r, boundary, point=np.zeros(3), offsets=None):
-
+    if issparse(r):
+        r = r.toarray()
     point = point.reshape(1, 3)
 
     indices = np.array([(1, 2), (2, 0), (0, 1)])

tests/test_mechanics_tensor.py

@@ -0,0 +1,86 @@
+# -*- coding: utf-8 -*-
+"""
+ _______  _______  ___      __   __  _______  _______ 
+|       ||       ||   |    |  | |  ||       ||       |
+|    ___||    ___||   |    |  | |  ||    _  ||    ___|
+|   |___ |   |___ |   |    |  |_|  ||   |_| ||   |___ 
+|    ___||    ___||   |___ |       ||    ___||    ___|
+|   |    |   |___ |       ||       ||   |    |   |___ 
+|___|    |_______||_______||_______||___|    |_______|
+
+This file is part of felupe.
+
+Felupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+Felupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Felupe.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+import pytest
+import felupe as fe
+import numpy as np
+
+
+def pre_umat():
+
+    LE = fe.LinearElastic(E=210000, nu=0.3)
+
+    class LETensor:
+        def __init__(self, LE):
+            self.LE = LE
+
+        def function(self, F, statevars):
+            # return dummy state variables along with stress
+            return self.LE.stress(F), statevars
+
+        def gradient(self, F, statevars):
+            return self.LE.elasticity(F)
+
+    return LETensor(LE)
+
+
+def test_solidbody_tensor():
+
+    umat = pre_umat()
+
+    m = fe.Cube(n=3)
+    r = fe.RegionHexahedron(m)
+    u = fe.Field(r, dim=3)
+
+    sv = np.zeros((5, 16, r.quadrature.npoints, m.ncells))
+
+    b = fe.SolidBodyTensor(umat, u, sv)
+    r = b.assemble.vector()
+
+    K = b.assemble.matrix(u)
+    K = b.assemble.matrix()
+    r = b.assemble.vector(u)
+    r = b.assemble.vector()
+    F = b.results.kinematics[0]
+    P = b.results.stress[0]
+    s = b.evaluate.cauchy_stress()
+    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.shape == (3, 3, 8, 8)
+    assert P.shape == (3, 3, 8, 8)
+    assert s.shape == (3, 3, 8, 8)
+    assert t.shape == (3, 3, 8, 8)
+    assert C.shape == (3, 3, 3, 3, 8, 8)
+    assert z.shape == (5, 16, 8, 8)
+
+
+if __name__ == "__main__":
+    test_solidbody_tensor()

tests/test_tools.py

@@ -79,8 +79,9 @@ def test_solve_check():
         gradient=W.gradient(F),
     )
 
-    force = fe.tools.force(u, b, bounds["symx"])
-    moment = fe.tools.moment(u, b, bounds["symx"])
+    for b in [L.assemble(), L.assemble().toarray(), L.assemble().toarray()[:, 0]]:
+        force = fe.tools.force(u, b, bounds["symx"])
+        moment = fe.tools.moment(u, b, bounds["symx"])
 
     for a in [2, 3, 4, 5]:
         curve = fe.tools.curve(np.arange(a), np.ones(a) * force[0])