Merge pull request #63 from adtzlr/add-uniaxial-loadcase

Add uniaxial loadcase

felupe/doftools.py

@@ -255,3 +255,21 @@ def __init__(
 
         self.dof = dof[self.mask]
         self.points = np.arange(mesh.npoints)[mask]
+
+
+def uniaxial(field, right=1, move=0.2, clamped=True):
+    "Define boundaries for uniaxial loading."
+
+    f1 = lambda x: np.isclose(x, right)
+
+    bounds = symmetry(field)
+
+    if clamped:
+        bounds["right"] = Boundary(field, fx=f1, skip=(1, 0, 0))
+
+    bounds["move"] = Boundary(field, fx=f1, skip=(0, 1, 1), value=move)
+
+    dof0, dof1 = partition(field, bounds)
+    ext0 = apply(field, bounds, dof0)
+
+    return bounds, dof0, dof1, ext0

felupe/tools.py

@@ -211,7 +211,8 @@ def newtonrhapson(
     """
 
     # copy x0
-    x = deepcopy(x0)
+    x = x0
+    # x = deepcopy(x0)
 
     # pre-evaluate function at given unknowns "x"
     f = fun(pre(x), *args, **kwargs)

tests/test_linearelastic.py

@@ -0,0 +1,77 @@
+# -*- 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 felupe as fe
+
+
+def test_linearelastic():
+
+    mesh = fe.mesh.Cube(n=3)
+    element = fe.element.Hexahedron()
+    quadrature = fe.quadrature.GaussLegendre(order=1, dim=3)
+    region = fe.Region(mesh, element, quadrature)
+    displacement = fe.Field(region, dim=3)
+    umat = fe.constitution.LinearElastic(E=1, nu=0.3)
+
+    boundaries, dof0, dof1, u0ext = fe.doftools.uniaxial(displacement, clamped=False)
+
+    def fun(strain):
+        stress = umat.stress(strain)
+        linearform = fe.IntegralForm(stress, displacement, region.dV, grad_v=True)
+        return linearform.assemble().toarray()[:, 0]
+
+    def jac(strain):
+        stress = umat.stress(strain)
+        elasti = umat.elasticity(strain, stress)
+        bilinearform = fe.IntegralForm(
+            elasti, displacement, region.dV, displacement, grad_v=True, grad_u=True
+        )
+        return bilinearform.assemble()
+
+    res = fe.tools.newtonrhapson(
+        fun,
+        displacement,
+        jac,
+        pre=fe.tools.strain,
+        solve=fe.tools.solve,
+        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,
+        },
+    )
+
+    fe.utils.save(region, res.x, filename="result.vtk")
+
+
+if __name__ == "__main__":
+    test_linearelastic()