Merge pull request #65 from adtzlr/add-automatic-differentation

Add automatic differentation capability in Constitution

felupe/__about__.py

@@ -1 +1 @@
-__version__ = "0.0.8"
+__version__ = "0.0.9"

felupe/constitution/__init__.py

@@ -13,3 +13,4 @@
 from .df_da0 import MaterialFrom
 from .models import LinearElastic, NeoHooke
 from . import variation
+from . import ad

felupe/constitution/ad.py

@@ -0,0 +1,163 @@
+# -*- 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 functools import partial
+import numpy as np
+import casadi as ca
+
+
+def apply(x, fun, fun_shape, trailing_axes=2):
+    "Helper function for the calculation of fun(x)."
+
+    # get shape of trailing axes
+    ax = x[0].shape[-trailing_axes:]
+
+    def rshape(z):
+        if len(z.shape) == trailing_axes:
+            return z.reshape(1, -1, order="F")
+        else:
+            return z.reshape(z.shape[0], -1, order="F")
+
+    # reshape input
+    y = [rshape(z) for z in x]
+
+    # map function to reshaped input
+    out = np.array(fun.map(np.product(ax))(*y))
+
+    return out.reshape(*fun_shape, *ax, order="F")
+
+
+class Material:
+    def __init__(self, W, *args, **kwargs):
+        "Material class (u) with Automatic Differentiation."
+
+        # init deformation gradient
+        F = ca.SX.sym("F", 3, 3)
+
+        # gradient and hessian of strain ernergy function W
+        d2WdFdF, dWdF = ca.hessian(W(F, *args, **kwargs), F)
+        # dWdF = ca.jacobian(W(F, *args, **kwargs), F)
+        # d2WdFdF = ca.jacobian(dWdF, F)
+
+        # generate casadi function objects
+        self._f_P = ca.Function("P", [F], [dWdF])
+        self._f_A = ca.Function("A", [F], [d2WdFdF])
+
+    # functions for stress P and elasticity A
+    def f(self, F):
+        fF = apply([F], fun=self._f_P, fun_shape=(3, 3))
+        return fF
+
+    def A(self, F):
+        AFF = apply([F], fun=self._f_A, fun_shape=(3, 3, 3, 3))
+        return AFF
+
+
+class Materialup:
+    def __init__(self, W, *args, **kwargs):
+        "Material class (u,p) with Automatic Differentiation."
+
+        # init deformation gradient
+        F = ca.SX.sym("F", 3, 3)
+        p = ca.SX.sym("p", 1)
+
+        # gradient and hessian of strain ernergy function W
+        w = W(F, p, *args, **kwargs)
+
+        d2WdF2, dWdF = ca.hessian(w, F)
+        d2Wdp2, dWdp = ca.hessian(w, p)
+
+        d2WdFdp = ca.jacobian(dWdF, p)
+
+        # generate casadi function objects
+        self._f_PF = ca.Function("PF", [F, p], [dWdF])
+        self._f_Pp = ca.Function("Pp", [F, p], [dWdp])
+
+        self._f_AFF = ca.Function("AFF", [F, p], [d2WdF2])
+        self._f_App = ca.Function("App", [F, p], [d2Wdp2])
+        self._f_AFp = ca.Function("AFp", [F, p], [d2WdFdp])
+
+    # functions for stress P and elasticity A
+    def f(self, F, p):
+        fF = apply([F, p], fun=self._f_PF, fun_shape=(3, 3))
+        fp = apply([F, p], fun=self._f_Pp, fun_shape=())
+        return [fF, fp]
+
+    def A(self, F, p):
+        AFF = apply([F, p], fun=self._f_AFF, fun_shape=(3, 3, 3, 3))
+        App = apply([F, p], fun=self._f_App, fun_shape=(1,))
+        AFp = apply([F, p], fun=self._f_AFp, fun_shape=(3, 3))
+        return [AFF, AFp, App]
+
+
+class MaterialupJ:
+    def __init__(self, W, *args, **kwargs):
+        "Material class (u,p,J) with Automatic Differentiation."
+
+        # init deformation gradient
+        F = ca.SX.sym("F", 3, 3)
+        p = ca.SX.sym("p", 1)
+        J = ca.SX.sym("J", 1)
+
+        # gradient and hessian of strain ernergy function W
+        w = W(F, p, J, *args, **kwargs)
+
+        d2WdF2, dWdF = ca.hessian(w, F)
+        d2Wdp2, dWdp = ca.hessian(w, p)
+        d2WdJ2, dWdJ = ca.hessian(w, J)
+
+        d2WdFdp = ca.jacobian(dWdF, p)
+        d2WdFdJ = ca.jacobian(dWdF, J)
+        d2WdpdJ = ca.jacobian(dWdp, J)
+
+        # generate casadi function objects
+        self._f_PF = ca.Function("PF", [F, p, J], [dWdF])
+        self._f_Pp = ca.Function("Pp", [F, p, J], [dWdp])
+        self._f_PJ = ca.Function("PJ", [F, p, J], [dWdJ])
+
+        self._f_AFF = ca.Function("AFF", [F, p, J], [d2WdF2])
+        self._f_App = ca.Function("App", [F, p, J], [d2Wdp2])
+        self._f_AJJ = ca.Function("AJJ", [F, p, J], [d2WdJ2])
+        self._f_AFp = ca.Function("AFp", [F, p, J], [d2WdFdp])
+        self._f_AFJ = ca.Function("AFp", [F, p, J], [d2WdFdJ])
+        self._f_ApJ = ca.Function("ApJ", [F, p, J], [d2WdpdJ])
+
+    # functions for stress P and elasticity A
+    def f(self, F, p, J):
+        fF = apply([F, p, J], fun=self._f_PF, fun_shape=(3, 3))
+        fp = apply([F, p, J], fun=self._f_Pp, fun_shape=(1,))
+        fJ = apply([F, p, J], fun=self._f_PJ, fun_shape=(1,))
+        return [fF, fp, fJ]
+
+    def A(self, F, p, J):
+        AFF = apply([F, p, J], fun=self._f_AFF, fun_shape=(3, 3, 3, 3))
+        App = apply([F, p, J], fun=self._f_App, fun_shape=(1,))
+        AJJ = apply([F, p, J], fun=self._f_AJJ, fun_shape=(1,))
+        AFp = apply([F, p, J], fun=self._f_AFp, fun_shape=(3, 3))
+        AFJ = apply([F, p, J], fun=self._f_AFJ, fun_shape=(3, 3))
+        ApJ = apply([F, p, J], fun=self._f_ApJ, fun_shape=(1,))
+        return [AFF, AFp, AFJ, App, ApJ, AJJ]

requirements.txt

@@ -3,4 +3,5 @@ scipy
 h5py
 meshio
 quadpy
-numba
+numba
+casadi

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = felupe
-version = 0.0.8
+version = 0.0.9
 author = Andreas Dutzler
 author_email = [email protected]
 description = Finite Element Analysis
@@ -35,4 +35,5 @@ install_requires =
     meshio
     quadpy
     numba
+    casadi
 python_requires = >=3.6

tests/test_mixed_with_incsolve.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 import felupe as fe
+import casadi as ca
 
 
 def test_mixed_with_incsolve():
@@ -42,5 +43,53 @@ def test_mixed_with_incsolve():
     fe.utils.save(region, fields, unstack=unstack, filename="result.vtk")
 
 
+def test_mixed_with_incsolve_ad_upJ():
+
+    mesh = fe.mesh.Cube(n=3)
+    mesh0 = fe.mesh.convert(mesh)
+
+    element0 = fe.element.ConstantHexahedron()
+    element1 = fe.element.Hexahedron()
+    quadrature = fe.quadrature.GaussLegendre(order=1, dim=3)
+
+    region0 = fe.Region(mesh0, element0, quadrature)
+    region = fe.Region(mesh, element1, quadrature)
+
+    displacement = fe.Field(region, dim=3)
+    pressure = fe.Field(region0, dim=1)
+    volumeratio = fe.Field(region0, dim=1, values=1)
+    fields = (displacement, pressure, volumeratio)
+
+    f1 = lambda x: np.isclose(x, 1)
+
+    boundaries = fe.doftools.symmetry(displacement)
+    boundaries["right"] = fe.Boundary(displacement, fx=f1, skip=(1, 0, 0))
+    boundaries["move"] = fe.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=0.1)
+
+    dof0, dof1, unstack = fe.doftools.partition(fields, boundaries)
+
+    def W(F, p, J, mu, bulk):
+        "Strain energy density function for Neo-Hookean material formulation."
+
+        detF = ca.det(F)
+
+        Fm = (J / detF) ** (1 / 3) * F
+
+        Cm = Fm.T @ Fm
+        I_Cm = ca.trace(Cm)
+
+        W_iso = mu / 2 * (I_Cm * J ** (-2 / 3) - 3)
+        W_vol = bulk / 2 * (J - 1) ** 2
+
+        return W_iso + W_vol + p * (detF - J)
+
+    umat = fe.constitution.ad.MaterialupJ(W, mu=1.0, bulk=500.0)
+
+    fe.tools.incsolve(fields, region, umat.f, umat.A, boundaries, [0.1, 0.2], verbose=0)
+
+    fe.utils.save(region, fields, unstack=unstack, filename="result.vtk")
+
+
 if __name__ == "__main__":
     test_mixed_with_incsolve()
+    test_mixed_with_incsolve_ad_upJ()

tests/test_readme_with_nr.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 import felupe as fe
+import casadi as ca
 
 
 def test_hex8_nh_nr():
@@ -66,5 +67,71 @@ def jac(F):
     fe.utils.save(region, displacement, filename=None)
 
 
+def test_hex8_nh_nr_ad():
+
+    mesh = fe.mesh.Cube(n=3)
+    element = fe.element.Hexahedron()
+    quadrature = fe.quadrature.GaussLegendre(order=1, dim=3)
+
+    region = fe.Region(mesh, element, quadrature)
+    dV = region.dV
+
+    displacement = fe.Field(region, dim=3)
+
+    def W(F, mu, bulk):
+        "Strain energy density function for Neo-Hookean material formulation."
+
+        J = ca.det(F)
+        C = F.T @ F
+        I_C = ca.trace(C)
+        return mu / 2 * (I_C * J ** (-2 / 3) - 3) + bulk / 2 * (J - 1) ** 2
+
+    umat = fe.constitution.ad.Material(W, mu=1.0, bulk=2.0)
+
+    f0 = lambda x: np.isclose(x, 0)
+    f1 = lambda x: np.isclose(x, 1)
+
+    boundaries = {}
+    boundaries["left"] = fe.Boundary(displacement, fx=f0)
+    boundaries["right"] = fe.Boundary(displacement, fx=f1, skip=(1, 0, 0))
+    boundaries["move"] = fe.Boundary(displacement, fx=f1, skip=(0, 1, 1), value=0.5)
+
+    dof0, dof1 = fe.doftools.partition(displacement, boundaries)
+    u0ext = fe.doftools.apply(displacement, boundaries, dof0)
+
+    def fun(F):
+        linearform = fe.IntegralForm(umat.f(F), displacement, dV, grad_v=True)
+        return linearform.assemble().toarray()[:, 0]
+
+    def jac(F):
+        bilinearform = fe.IntegralForm(
+            umat.A(F), displacement, dV, displacement, grad_v=True, grad_u=True
+        )
+        return bilinearform.assemble()
+
+    res = fe.tools.newtonrhapson(
+        fun,
+        displacement,
+        jac,
+        solve=fe.tools.solve,
+        pre=fe.tools.defgrad,
+        update=fe.tools.update,
+        check=fe.tools.check,
+        kwargs_solve={"field": displacement, "ext": u0ext, "dof0": dof0, "dof1": dof1},
+        kwargs_check={
+            "tol_f": 1e-3,
+            "tol_x": 1e-3,
+            "dof0": dof0,
+            "dof1": dof1,
+            "verbose": 0,
+        },
+    )
+
+    displacement = res.x
+
+    fe.utils.save(region, displacement, filename=None)
+
+
 if __name__ == "__main__":
     test_hex8_nh_nr()
+    test_hex8_nh_nr_ad()