diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml index 5f5e26a2..6096ea6e 100644 --- a/conda.recipe/meta.yaml +++ b/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: diff --git a/environment.yml b/environment.yml index 1eb321c8..b72cab6c 100644 --- a/environment.yml +++ b/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 diff --git a/environment_doc.yml b/environment_doc.yml index af7528f1..814b3824 100644 --- a/environment_doc.yml +++ b/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 diff --git a/fedoo/constitutivelaw/composite_ud.py b/fedoo/constitutivelaw/composite_ud.py index f6b710a9..6283ece8 100644 --- a/fedoo/constitutivelaw/composite_ud.py +++ b/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": diff --git a/fedoo/core/mechanical3d.py b/fedoo/core/mechanical3d.py index c703b0aa..92f4f6a7 100644 --- a/fedoo/core/mechanical3d.py +++ b/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 diff --git a/fedoo/util/localframe.py b/fedoo/util/localframe.py index b269c524..98c64a5b 100644 --- a/fedoo/util/localframe.py +++ b/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] diff --git a/fedoo/weakform/stress_equilibrium.py b/fedoo/weakform/stress_equilibrium.py index 0f729fca..9ebebde0 100644 --- a/fedoo/weakform/stress_equilibrium.py +++ b/fedoo/weakform/stress_equilibrium.py @@ -6,6 +6,7 @@ try: from simcoon import simmit as sim + from simcoon import Rotation as SimRotation USE_SIMCOON = True except ModuleNotFoundError: @@ -282,27 +283,23 @@ def set_start(self, assembly, pb): """Start a new time increment.""" if assembly._nlgeom: if "DStrain" in assembly.sv: - # rotate strain and stress -> need to be checked + # rotate strain and stress + rot = SimRotation.from_matrix(assembly.sv["DR"].transpose(2, 0, 1)) assembly.sv["Strain"] = StrainTensorList( - sim.rotate_strain_R( + rot.apply_strain( assembly.sv_start["Strain"].asarray(), - assembly.sv["DR"], ) + assembly.sv["DStrain"] ) assembly.sv["DStrain"] = StrainTensorList( np.zeros((6, assembly.n_gauss_points), order="F") ) - # or assembly.sv['DStrain'] = 0 perhaps more efficient to avoid a - # nul sum - - # update cauchy stress - if not (np.array_equal(assembly.sv["DispGradient"], 0)): - # True when the problem have been updated once - stress = assembly.sv["Stress"].asarray() - assembly.sv["Stress"] = StressTensorList( - sim.rotate_stress_R(stress, assembly.sv["DR"]) - ) + + # update cauchy stress + if not (np.array_equal(assembly.sv["DispGradient"], 0)): + # True when the problem have been updated once + stress = assembly.sv["Stress"].asarray() + assembly.sv["Stress"] = StressTensorList(rot.apply_stress(stress)) if assembly._nlgeom == "TL": assembly.sv["PK2"] = assembly.sv["Stress"].cauchy_to_pk2( assembly.sv["F"] diff --git a/pyproject.toml b/pyproject.toml index 3fcb1e4b..4dd671b2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,10 +27,11 @@ classifiers = [ dependencies = [ 'numpy>=2.00', 'scipy', + 'simcoon>=1.11.0', ] [project.optional-dependencies] -all = ['fedoo[solver,plot,simcoon,test,ipc]'] +all = ['fedoo[solver,plot,test,ipc]'] ipc = ['ipctk'] solver = [ 'pypardiso ; platform_machine!="arm64" and platform_machine!="aarch64"', @@ -40,9 +41,6 @@ plot = [ 'matplotlib', 'pyvista[io]', ] -simcoon = [ - 'simcoon' -] test = [ 'pytest', 'pytest-cov'