diff --git a/felupe/__about__.py b/felupe/__about__.py index a73339bf..00ec2dcd 100644 --- a/felupe/__about__.py +++ b/felupe/__about__.py @@ -1 +1 @@ -__version__ = "0.0.8" +__version__ = "0.0.9" diff --git a/felupe/constitution/__init__.py b/felupe/constitution/__init__.py index 81af0ae8..c3d76b0e 100644 --- a/felupe/constitution/__init__.py +++ b/felupe/constitution/__init__.py @@ -13,3 +13,4 @@ from .df_da0 import MaterialFrom from .models import LinearElastic, NeoHooke from . import variation +from . import ad diff --git a/felupe/constitution/ad.py b/felupe/constitution/ad.py new file mode 100644 index 00000000..b486a03c --- /dev/null +++ b/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 . + +""" + +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] diff --git a/requirements.txt b/requirements.txt index 5d776556..970078c3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,4 +3,5 @@ scipy h5py meshio quadpy -numba \ No newline at end of file +numba +casadi \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 6c3ca9a4..39c41ba8 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = felupe -version = 0.0.8 +version = 0.0.9 author = Andreas Dutzler author_email = a.dutzler@gmail.com description = Finite Element Analysis @@ -35,4 +35,5 @@ install_requires = meshio quadpy numba + casadi python_requires = >=3.6 \ No newline at end of file diff --git a/tests/test_mixed_with_incsolve.py b/tests/test_mixed_with_incsolve.py index 265a0af5..7a8a7010 100644 --- a/tests/test_mixed_with_incsolve.py +++ b/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() diff --git a/tests/test_readme_with_nr.py b/tests/test_readme_with_nr.py index 02d4c879..ac803a31 100644 --- a/tests/test_readme_with_nr.py +++ b/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()