Feature/rotation

conda.recipe/meta.yaml

@@ -21,7 +21,7 @@ requirements:
     - numpy>=1.26
     - scipy
     - pyvista>=0.42
-    - simcoon>=1.9.6
+    - simcoon>=1.11.0
     - {{ solver_backend }}
 
 test:

environment.yml

@@ -12,7 +12,7 @@ dependencies:
   - scikit-umfpack # [(not x86) and (not x86_64)]
   - vtk>=9.2
   - pyvista>=0.39
-  - simcoon>=1.9.6
+  - simcoon>=1.11.0
   - pytest
   - pytest-xdist
   - sphinx_rtd_theme

environment_doc.yml

@@ -13,7 +13,7 @@ dependencies:
   - vtk>=9.2
   - pypandoc
   - pyvista>=0.39
-  - simcoon>=1.9.6
+  - simcoon>=1.11.0
   - pytest
   - sphinx
   - sphinx_rtd_theme

fedoo/constitutivelaw/composite_ud.py

@@ -3,6 +3,7 @@
 
 from fedoo.core.mechanical3d import Mechanical3D
 from fedoo.constitutivelaw.elastic_anisotropic import ElasticAnisotropic
+from simcoon import Rotation
 
 import numpy as np
 
@@ -194,34 +195,18 @@ def get_tangent_matrix(self, assembly, dimension=None):
         if not (
             np.isscalar(self.__parameters["angle"]) and self.__parameters["angle"] == 0
         ):
-            # angle in degree
-            angle_pli = self.__parameters["angle"] / 180.0 * np.pi
-            s = np.sin(angle_pli)
-            c = np.cos(angle_pli)
-            zero = 0 * s
-            one = zero + 1
-
-            R_epsilon = np.array(
-                [
-                    [c**2, s**2, zero, s * c, zero, zero],
-                    [s**2, c**2, zero, -s * c, zero, zero],
-                    [zero, zero, one, zero, zero, zero],
-                    [-2 * s * c, 2 * s * c, zero, c**2 - s**2, zero, zero],
-                    [zero, zero, zero, zero, c, s],
-                    [zero, zero, zero, zero, -s, c],
-                ]
+            rot = Rotation.from_euler(
+                "Z", np.atleast_1d(self.__parameters["angle"]), degrees=True
             )
-
-            if len(R_epsilon.shape) == 3:
-                R_epsilon = np.transpose(R_epsilon, [2, 0, 1])
-                R_sigma_inv = np.transpose(R_epsilon, [0, 2, 1])
-            else:
-                R_sigma_inv = R_epsilon.T
-            if len(H.shape) == 3:
+            QS = rot.as_voigt_stress_rotation()  # (N, 6, 6)
+            if H.ndim == 3:
                 H = np.rollaxis(H, 2, 0)
-            H = np.matmul(R_sigma_inv, np.matmul(H, R_epsilon))
-            if len(H.shape) == 3:
-                H = np.rollaxis(H, 0, 3)
+            H = np.matmul(QS, np.matmul(H, QS.transpose(0, 2, 1)))
+            if H.ndim == 3:
+                if H.shape[0] == 1:
+                    H = H[0]
+                else:
+                    H = np.rollaxis(H, 0, 3)
 
         H = self.local2global_H(H)
         if dimension == "2Dstress":

fedoo/core/mechanical3d.py

@@ -2,6 +2,7 @@
 
 # baseclass
 import numpy as np
+from simcoon import Rotation
 from fedoo.core.base import ConstitutiveLaw
 
 
@@ -68,42 +69,27 @@ def get_H_plane_stress(self, H):
             for i in range(6)
         ]
 
-    def local2global_H(self, H_global):
-        # Change of basis capability for laws on the form : StressTensor = H * StrainTensor
-        # StressTensor and StrainTensor are column vectors based on the voigt notation
+    def local2global_H(self, H):
+        """Rotate stiffness matrix from local material frame to global frame.
+
+        Uses simcoon.Rotation to build the 6x6 Voigt stress rotation matrix
+        QS from the local frame, then computes H_global = QS @ H @ QS^T.
+        """
         if self.local_frame is not None:
-            # building the matrix to change the basis of the stress and the strain
-            #            theta = np.pi/8
-            #            np.array([[np.cos(theta),np.sin(theta),0], [-np.sin(theta),np.cos(theta),0], [0,0,1]])
-            if len(self.local_frame.shape) == 2:
-                self.local_frame = self.local_frame.reshape(1, 3, 3)
-            R_epsilon = np.empty((len(self.local_frame), 6, 6))
-            R_epsilon[:, :3, :3] = self.local_frame**2
-            R_epsilon[:, :3, 3:6] = (
-                self.local_frame[:, :, [0, 2, 1]] * self.local_frame[:, :, [1, 0, 2]]
-            )
-            R_epsilon[:, 3:6, :3] = (
-                2 * self.local_frame[:, [0, 2, 1]] * self.local_frame[:, [1, 0, 2]]
-            )
-            R_epsilon[:, 3:6, 3:6] = (
-                self.local_frame[:, [[0], [2], [1]], [0, 2, 1]]
-                * self.local_frame[:, [[1], [0], [2]], [1, 0, 2]]
-                + self.local_frame[:, [[1], [0], [2]], [0, 2, 1]]
-                * self.local_frame[:, [[0], [2], [1]], [1, 0, 2]]
-            )
-            R_sigma_inv = R_epsilon.transpose(
-                0, 2, 1
-            )  # np.transpose(R_epsilon,[0,2,1])
-
-            if len(H_global.shape) == 3:
-                H_global = np.rollaxis(H_global, 2, 0)
-            H_local = np.matmul(R_sigma_inv, np.matmul(H_global, R_epsilon))
-            if len(H_local.shape) == 3:
-                if H_local.shape[0] == 1:
-                    return H_local[0, :, :]
-                H_local = np.rollaxis(H_local, 0, 3)
-
-            return H_local
-
-        else:
-            return H_global
+            local_frame = np.asarray(self.local_frame)
+            if local_frame.ndim == 2:
+                local_frame = local_frame[np.newaxis]
+            rot = Rotation.from_matrix(local_frame)
+            QS = rot.as_voigt_stress_rotation()  # (N, 6, 6)
+
+            H = np.asarray(H, dtype=float)
+            if H.ndim == 3:
+                H = np.rollaxis(H, 2, 0)  # (6,6,M) -> (M,6,6)
+            H = np.matmul(QS, np.matmul(H, QS.transpose(0, 2, 1)))
+            if H.ndim == 3:
+                if H.shape[0] == 1:
+                    return H[0]
+                return np.rollaxis(H, 0, 3)  # (N,6,6) -> (6,6,N)
+            return H
+
+        return H

fedoo/util/localframe.py

@@ -1,84 +1,52 @@
-# from fedoo.pgd.SeparatedArray import SeparatedArray
 from fedoo.core.mesh import Mesh
+from simcoon import Rotation
 
-# from fedoo.pgd.MeshPGD import MeshPGD
 import numpy as np
 
 
 class LocalFrame(np.ndarray):
-    """
-    The class LocalFrame is derived from the np.ndarray class including some additional function to treat local frame.
-    A LocalFrame object should be a (N,3,3) shaped array in 3D or (N,2,2) shaped array in 2D,
-    where N is the number of points (nodes, gauss point, elements, ...) in which local frames are defined.
-    If LF is a LocalFrame object:
-        - LF[i] gives the local frame associated to the ith point.
-        - LF[i][0], LF[i][1] and LF[i][2]  gives respectively the 3 vectors defining the local frame.
-          These 3 vectors are assumed unit and orthogonal.
+    """Array of local coordinate frames.
+
+    A LocalFrame object is an (N, 3, 3) shaped array in 3D or (N, 2, 2) in 2D,
+    where N is the number of points (nodes, gauss points, elements, ...).
+
+    If ``LF`` is a LocalFrame object:
+        - ``LF[i]`` gives the local frame at the i-th point.
+        - ``LF[i][0]``, ``LF[i][1]`` and ``LF[i][2]`` give the 3 unit
+          orthogonal vectors defining the local frame.
     """
 
-    def __new__(self, localFrame):
-        return np.asarray(localFrame).view(self)
+    def __new__(cls, localFrame):
+        return np.asarray(localFrame).view(cls)
 
     def Rotate(self, angle, axis="Z"):
-        # angle in degree
-        if angle is not 0:
-            angle = angle / 180.0 * np.pi
-
-            if axis.upper() == "X":
-                axis = (1, 0, 0)
-            elif axis.upper() == "Y":
-                axis = (0, 1, 0)
-            elif axis.upper() == "Z":
-                axis = (0, 0, 1)
-
-            s = np.sin(angle)
-            c = np.cos(angle)
-            x = axis[0]
-            y = axis[1]
-            z = axis[2]
-            axis = np.array([axis])
-
-            RotMatrix = np.array(
-                [
-                    [
-                        (1 - c) * x**2 + c,
-                        (1 - c) * x * y - z * s,
-                        (1 - c) * x * z + y * s,
-                    ],
-                    [
-                        (1 - c) * x * y + z * s,
-                        (1 - c) * y**2 + c,
-                        (1 - c) * y * z - x * s,
-                    ],
-                    [
-                        (1 - c) * x * z - y * s,
-                        (1 - c) * y * z + x * s,
-                        (1 - c) * z**2 + c,
-                    ],
-                ]
-            )
-
-            if len(RotMatrix.shape) == 3:
-                RotMatrix = np.transpose(RotMatrix, [2, 1, 0])
-            elif len(RotMatrix.shape) == 2:
-                RotMatrix = RotMatrix.T
-
-            self[:] = np.matmul(self, RotMatrix)
+        """Rotate all local frames by a given angle around a global axis.
+
+        Parameters
+        ----------
+        angle : float
+            Rotation angle in degrees.
+        axis : str, optional
+            Axis of rotation: 'X', 'Y' or 'Z' (default 'Z').
+
+        Returns
+        -------
+        self
+        """
+        if angle != 0:
+            rot = Rotation.from_euler(axis.upper(), angle, degrees=True)
+            self[:] = np.matmul(self, rot.as_matrix().T)
             return self
 
-    #            np.dot(axis.T,axis)*(1-c) \
-    #                 + np.array([[ c        ,-axis[2]*s, axis[1]*s],
-    #                             [ axis[2]*s, c        ,-axis[0]*s],
-    #                         [-axis[2]*s, axis[0]*s, c        ]])
+    def as_rotation(self):
+        """Convert local frames to a batch simcoon.Rotation object.
 
-    #            if len(RotMatrix.shape) == 3:
-    #                RotMatrix = np.transpose(R_epsilon,[2,0,1])
-    #            else: R_sigma_inv = R_epsilon.T
-    #            if len(H.shape) == 3: H = np.rollaxis(H,2,0)
-    #            H = np.matmul(R_sigma_inv, np.matmul(H,R_epsilon))
-    #            if len(H.shape) == 3: H = np.rollaxis(H,0,3)
-    #
-    #        return H
+        Returns
+        -------
+        simcoon.Rotation
+            Batch rotation from local to global frame.
+        """
+        return Rotation.from_matrix(np.asarray(self))
 
     def __getitem__(self, index):
         new = super(LocalFrame, self).__getitem__(index)
@@ -89,10 +57,28 @@ def __getitem__(self, index):
 
 
 def global_local_frame(n_points):
-    return LocalFrame([np.eye(3) for i in range(n_points)])
+    """Return identity local frames (global = local) for n_points."""
+    return LocalFrame(np.tile(np.eye(3), (n_points, 1, 1)))
 
 
 def GenerateCylindricalLocalFrame(crd, axis=2, origin=[0, 0, 0], dim=3):
+    """Generate cylindrical local frames (er, etheta, ez) at each node.
+
+    Parameters
+    ----------
+    crd : array_like or Mesh
+        Node coordinates, shape (N, 3) or (N, 2).
+    axis : int, optional
+        Cylinder axis: 0=X, 1=Y, 2=Z (default 2).
+    origin : array_like, optional
+        Origin of the cylindrical coordinate system.
+    dim : int, optional
+        Spatial dimension: 2 or 3 (default 3).
+
+    Returns
+    -------
+    LocalFrame
+    """
     if isinstance(crd, Mesh):
         crd = Mesh.nodes
 
@@ -106,7 +92,7 @@ def GenerateCylindricalLocalFrame(crd, axis=2, origin=[0, 0, 0], dim=3):
         crd = crd[:, 0:2]
         origin = np.array(origin)[0:2]
 
-    crd = crd - np.array(origin).reshape(1, -1)  # changement of origin
+    crd = crd - np.array(origin).reshape(1, -1)
 
     localFrame[:, 0, plane] = crd[:, plane] / np.sqrt(
         crd[:, plane[0]] ** 2 + crd[:, plane[1]] ** 2
@@ -115,68 +101,8 @@ def GenerateCylindricalLocalFrame(crd, axis=2, origin=[0, 0, 0], dim=3):
     if dim == 3:
         localFrame[:, 1] = np.cross(
             localFrame[:, 2], localFrame[:, 0]
-        )  # etheta is the cross product
+        )  # etheta = ez x er
     else:
         localFrame[:, 1, 0] = -localFrame[:, 0, 1]
         localFrame[:, 1, 1] = localFrame[:, 0, 0]
     return localFrame.view(LocalFrame)
-
-
-# MOVE separated_local_frame to the pgd folder
-
-# def separated_local_frame(localFrame, mesh, dimensions = ('X','Y','Z')):
-#     """
-#     Permit to automatically assign the localFrame to the appropriate submesh of the mesh object
-#     Generate a local frame under the form of the (3,3) shaped array dedicated of SeparatedArray objects
-#     This functions work only if the local frame is restricted to one subspace.
-#     """
-
-#     dim = localFrame.shape[-1]
-#     idmesh = mesh.FindCoordinatename(dimensions[0])
-#     if idmesh != mesh.FindCoordinatename(dimensions[1]): raise NameError("'{}' and '{}' coordinates should be associated to the same subMesh".format(dimensions[0], dimensions[1]))
-#     if dim == 3:
-#         if idmesh != mesh.FindCoordinatename(dimensions[3]):  raise NameError("'{}' and '{}' coordinates should be associated to the same subMesh. Consider using a 2D local frame.".format(dimensions[0], dimensions[2]))
-
-#     id_crd = []
-#     for label in dimensions:
-#         if label == 'X': id_crd.append(0)
-#         elif label == 'Y': id_crd.append(1)
-#         elif label == 'Z': id_crd.append(2)
-#         else: raise NameError("Coordinates for local frame should be 'X', 'Y' or 'Z'. '{}' unknown.".format(label))
-
-#     newLocalFrame = np.zeros((3, 3), dtype =object) #the resulting local frame is always in dim = 3
-
-#     for j in range(dim):
-#         for i in range(dim):
-#             newLocalFrame[i,id_crd[j]] = SeparatedArray([np.c_[localFrame[:,i,j]] if k==idmesh else np.array([[1.]]) for k in range(mesh.get_dimension())])
-#     return newLocalFrame
-
-
-#   TODO: Not finished
-#   if isinstance(crd, MeshPGD):
-#        mesh = crd
-#        crd_all = [] #list containing the values of coordinates for each nodes of the separated mesh using 3 SeparatedArray objects
-#        for ii, namecrd in enumerate(['X','Y','Z']):
-#            idmesh = mesh.FindCoordinatename(namecrd)
-#            subMesh = mesh.GetListMesh()[idmesh]
-#            crd = subMesh.nodes[:, subMesh.crd_name.index(namecrd)]
-#            crd_all.append(SeparatedArray([np.c_[crd] if i == idmesh else np.array([[1.]]) for i in range(mesh.get_dimension())]))
-#
-#        localFrame = np.zeros((dim, dim), dtype =object)
-#
-#        plane = [0,1,2] ; plane.pop(axis)
-#        localFrame[2, axis] =  SeparatedArray([np.array([[1.]]) for i in range(mesh.get_dimension())])#ez
-#
-#        crd_all[0] = crd_all[0] - origin[0] #changement of origin
-#        crd_all[1] = crd_all[1] - origin[1] #changement of origin
-#        crd_all[2] = crd_all[2] - origin[2] #changement of origin
-#
-#
-#
-#        localFrame[:, 0, plane] = crd[:,plane]/sqrt(crd_all[plane[0]]**2+crd_all[plane[1]]**2) #er
-#
-#        if dim == 3:
-#            localFrame[:, 1] = np.cross(localFrame[:,2], localFrame[:,0]) #etheta is the cross product
-#        else:
-#            localFrame[:, 1, 0] = -localFrame[:,0, 1]
-#            localFrame[:, 1, 1] = localFrame[:,0, 0]