Improve AreaChange

README.md

@@ -56,6 +56,7 @@ All notable changes to this project will be documented in this file. The format
 - Add parallel versions of math helpers (`dya`, `cdya`, `dot`, `ddot`) using [einsumt](https://pypi.org/project/einsumt/).
 - Add `parallel` keyword to constitutive models (`NeoHooke`, `LinearElasticTensorNotation` and `ThreeFieldVariation`).
 - Add `RegionBoundary` along with template regions for `Quad` and `Hexahedron` and `GaussLegendreBoundary`.
+- Add optional normal vector argument for function and gradient methods of `AreaChange`.
 
 ### Changed
 - Enforce consistent arguments for functions inside `mesh.tools` (`points, cells, cell_data` or `Mesh`).

felupe/constitution/_kinematics.py

@@ -25,8 +25,17 @@
 
 """
 
+import numpy as np
+
+try:
+    from einsumt import einsumt
+except:
+    print("ImportWarning: Module `einsumt` not found. Fall back to `np.einsum()`.")
+    from numpy import einsum as einsumt
+
 from ..math import (
     transpose,
+    dot,
     inv,
     dya,
     cdya_ik,
@@ -104,29 +113,39 @@ class AreaChange:
     def __init__(self, parallel=False):
         self.parallel = parallel
 
-    def function(self, F):
+    def function(self, F, N=None):
         """Area change.
 
         Arguments
         ---------
         F : ndarray
             Deformation gradient
+        N : ndarray or None, optional
+            Area normal vector (default is None)
 
         Returns
         -------
         ndarray
             Cofactor matrix of the deformation gradient
         """
         J = det(F)
-        return J * transpose(inv(F, J))
 
-    def gradient(self, F):
+        Fs = J * transpose(inv(F, J))
+
+        if N is None:
+            return Fs
+        else:
+            return dot(Fs, N, parallel=self.parallel)
+
+    def gradient(self, F, N=None):
         """Gradient of area change.
 
         Arguments
         ---------
         F : ndarray
             Deformation gradient
+        N : ndarray or None, optional
+            Area normal vector (default is None)
 
         Returns
         -------
@@ -135,12 +154,23 @@ def gradient(self, F):
         """
 
         J = det(F)
+
         dJdF = self.function(F)
-        return (
+        dFsdF = (
             dya(dJdF, dJdF, parallel=self.parallel)
             - cdya_il(dJdF, dJdF, parallel=self.parallel)
         ) / J
 
+        if self.parallel:
+            einsum = einsumt
+        else:
+            einsum = np.einsum
+
+        if N is None:
+            return dFsdF
+        else:
+            return einsum("ijkl...,j...->ikl...", dFsdF, N)
+
 
 class VolumeChange:
     r"""Volume Change.

tests/test_constitution.py

@@ -197,6 +197,8 @@ def test_linear_planestrain():
 def test_kinematics():
     r, F = pre(sym=False, add_identity=True)
 
+    N = F[:, 0]
+
     for parallel in [False, True]:
 
         lc = fe.constitution.LineChange(parallel=parallel)
@@ -206,6 +208,9 @@ def test_kinematics():
         xf = lc.function(F)
         xg = lc.gradient(F)
 
+        Yf = ac.function(F, N)
+        Yg = ac.gradient(F, N)
+
         yf = ac.function(F)
         yg = ac.gradient(F)
 
@@ -221,6 +226,9 @@ def test_kinematics():
         assert yf.shape == (3, 3, *F.shape[-2:])
         assert yg.shape == (3, 3, 3, 3, *F.shape[-2:])
 
+        assert Yf.shape == (3, *F.shape[-2:])
+        assert Yg.shape == (3, 3, 3, *F.shape[-2:])
+
         assert zf.shape == F.shape[-2:]
         assert zg.shape == (3, 3, *F.shape[-2:])
         assert zh.shape == (3, 3, 3, 3, *F.shape[-2:])

tests/test_region.py

@@ -37,7 +37,7 @@ def test_region():
     r = fe.RegionQuad(mesh)
     r = fe.RegionQuadBoundary(mesh)
     r = fe.RegionConstantQuad(mesh)
-    
+
     mesh.cell_type = "some_fancy_cell_type"
     with pytest.raises(NotImplementedError):
         r = fe.RegionBoundary(mesh, fe.Quad(), fe.GaussLegendreBoundary(order=1, dim=2))