adtzlr · adtzlr · Oct 28, 2024 · Oct 28, 2024 · Oct 28, 2024 · Oct 28, 2024
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -8,6 +8,8 @@ All notable changes to this project will be documented in this file. The format
 - Add the `order`-argument to `FieldContainer.extract(order="C")` as well as for `Field`, `FieldAxisymmetric`, `FieldPlaneStrain` to return C-contiguous arrays by default.
 - Add an optional multiplier to `Laplace(multiplier=1.0)`.
 - Add optional keyword-arguments to `math.transpose(**kwargs)` to support optional `out` and `order`-keywords.
+- Add the attribute `RegionBoundary.tangents`, which contains a list of tangent unit vectors. For `quad` cell-types the length of this list is one and for `hexahedron` cell-types it is of length two.
+- Add `math.inplane(A, vectors)` to return the in-plane components of a symmetric tensor `A`, where the plane is defined by its standard unit vectors.
 
 ### Changed
 - Change default `np.einsum(..., order="K")` to `np.einsum(..., order="C")` in the methods of `Field`, `FieldAxisymmetric`, `FieldPlaneStrain` and `FieldContainer`.

diff --git a/docs/felupe/math.rst b/docs/felupe/math.rst
@@ -34,6 +34,7 @@ This module contains math functions.
    math.identity
    math.sym
    math.dya
+   math.inplane
    math.inv
    math.det
    math.dev
@@ -65,4 +66,4 @@ This module contains math functions.
 **Detailed API Reference**
 
 .. automodule:: felupe.math
-   :members: linsteps, rotation_matrix, deformation_gradient, strain, strain_stretch_1d, extract, values, norm, interpolate, grad, identity, sym, dya, inv, det, dev,cof, eig, eigh, eigvals, eigvalsh, equivalent_von_mises, transpose, majortranspose, trace, cdya_ik, cdya_il, cdya, cross, dot, ddot, tovoigt, reshape, ravel, solve_2d
+   :members: linsteps, rotation_matrix, deformation_gradient, strain, strain_stretch_1d, extract, values, norm, interpolate, grad, identity, sym, dya, inplane, inv, det, dev,cof, eig, eigh, eigvals, eigvalsh, equivalent_von_mises, transpose, majortranspose, trace, cdya_ik, cdya_il, cdya, cross, dot, ddot, tovoigt, reshape, ravel, solve_2d
diff --git a/src/felupe/math/__init__.py b/src/felupe/math/__init__.py
@@ -33,6 +33,7 @@
     eigvalsh,
     equivalent_von_mises,
     identity,
+    inplane,
     inv,
     majortranspose,
     ravel,
@@ -74,6 +75,7 @@
     "eigvalsh",
     "equivalent_von_mises",
     "identity",
+    "inplane",
     "inv",
     "majortranspose",
     "ravel",

diff --git a/src/felupe/math/_tensor.py b/src/felupe/math/_tensor.py
@@ -1481,3 +1481,55 @@ def equivalent_von_mises(A):
 
     devA = dev(A)
     return np.sqrt(3 / 2 * ddot(devA, devA))
+
+
+def inplane(A, vectors, **kwargs):
+    r"""Return the in-plane components of symmetric 2x2 or 3x3 tensors, where the planes
+    are defined by their standard unit vectors.
+
+    Parameters
+    ----------
+    A : ndarray of shape (N, N, ...)
+        Symmetric second-order tensors.
+    vectors : list of ndarray of shape (N, ...)
+        List of standard unit vectors of the planes.
+    **kwargs : dict, optional
+        Optional keyword-arguments for :func:`numpy.einsum`.
+
+    Notes
+    -----
+    The first two axes of the tensor and the first axis of each vector are the
+    tensor/vector dimensions and all remaining trailing axes are treated as batch
+    dimensions.
+
+    ..  math::
+
+        \boldsymbol{A}_{N-1} = \sum_{\alpha, \beta}
+            \left( \boldsymbol{E}_\alpha^T \boldsymbol{A} \boldsymbol{E}_\beta\right)
+            \boldsymbol{E}_\alpha \otimes \boldsymbol{E}_\beta
+
+    Returns
+    -------
+    ndarray of shape (N - 1, N - 1, ...)
+        The in-plane components of symmetric (N, N) tensors.
+
+    Examples
+    --------
+    >>> import felupe as fem
+    >>>
+    >>> A = np.array([[0, 3, 5], [3, 1, 4], [5, 4, 2]])
+    >>> vectors = [[1, 0, 0], [0, 1, 0]]
+    >>>
+    >>> A_inplane = fem.math.inplane(A, vectors)
+    >>> A_inplane
+    array([[0, 3],
+           [3, 1]])
+    """
+
+    vectors = np.array(vectors)
+
+    i, j = np.triu_indices(len(vectors))
+    ij = np.zeros((len(vectors), len(vectors)), dtype=int)
+    ij[i, j] = ij[j, i] = np.arange(len(i), dtype=int)
+
+    return np.einsum("ij...,ai...,aj...->a...", A, vectors[i], vectors[j], **kwargs)[ij]
diff --git a/src/felupe/region/_boundary.py b/src/felupe/region/_boundary.py
@@ -435,11 +435,13 @@ def __init__(
         super().__init__(mesh_boundary_cells, element, quadrature, grad=grad)
 
         if grad:
-            self.dA, self.dV, self.normals = self._init_faces()
+            self.dA, self.dV, self.normals, self.tangents = self._init_faces()
 
     def _init_faces(self):
         "Initialize (norm of) face normals of cells."
 
+        tangents = []
+
         if (
             self.mesh.cell_type == "quad"
             or self.mesh.cell_type == "quad8"
@@ -448,13 +450,18 @@ def _init_faces(self):
             dA_1 = self.dXdr[:, 0][::-1]
             dA_1[0] = -dA_1[0]
 
+            tangents.append(self.dXdr[:, 0] / np.linalg.norm(self.dXdr[:, 0], axis=0))
+
         elif (
             self.mesh.cell_type == "hexahedron"
             or self.mesh.cell_type == "hexahedron20"
             or self.mesh.cell_type == "hexahedron27"
         ):
             dA_1 = cross(self.dXdr[:, 0], self.dXdr[:, 1])
 
+            tangents.append(self.dXdr[:, 0] / np.linalg.norm(self.dXdr[:, 0], axis=0))
+            tangents.append(self.dXdr[:, 1] / np.linalg.norm(self.dXdr[:, 1], axis=0))
+
         dA = -dA_1 * self.quadrature.weights.reshape(-1, 1)
 
         # norm and unit normal vector
@@ -466,7 +473,7 @@ def _init_faces(self):
                 dA = np.pad(dA, ((0, 1), (0, 0), (0, 0)))
                 normals = np.pad(normals, ((0, 1), (0, 0), (0, 0)))
 
-        return dA, dV, normals
+        return dA, dV, normals, tangents
 
     def mesh_faces(self):
         "Return a Mesh with face-cells on the selected boundary."

diff --git a/tests/test_math.py b/tests/test_math.py
@@ -185,7 +185,18 @@ def test_math_linsteps():
     assert np.allclose(steps[-1], (0, 5))
 
 
+def test_inplane():
+    import felupe as fem
+
+    A = np.array([[0, 3, 5], [3, 1, 4], [5, 4, 2]], dtype=float)
+    vectors = [[1, 0, 0], [0, 1, 0]]
+
+    A_inplane = fem.math.inplane(A, vectors)
+    assert np.allclose(A_inplane.ravel(), [0, 3, 3, 1])
+
+
 if __name__ == "__main__":
     test_math()
     test_math_field()
     test_math_linsteps()
+    test_inplane()