diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 7b13f941..e7debad0 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -72,3 +72,10 @@ jobs: - name: Run Python tests run: uv run pytest + + - name: Install fedoo and run integration tests (Unix) + if: runner.os != 'Windows' + run: | + # uv pip install "fedoo>=0.8.0" # uncomment when fedoo 0.8.0 is released + uv pip install "git+https://github.com/3MAH/fedoo.git@feature/rotation" + uv run pytest tests/test_fedoo_integration.py -v diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1e2aaee8..f7886d5c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -89,3 +89,12 @@ jobs: name: test-logs-${{ matrix.os }} path: build/Testing/Temporary/LastTest.log retention-days: 7 + + - name: Install simcoon + fedoo and run integration tests (Unix) + if: runner.os != 'Windows' + run: | + pip install scikit-build-core pybind11 numpy + pip install . --no-build-isolation + # pip install "fedoo>=0.8.0" pytest # uncomment when fedoo 0.8.0 is released + pip install "git+https://github.com/3MAH/fedoo.git@feature/rotation" pytest + pytest tests/test_fedoo_integration.py -v diff --git a/pyproject.toml b/pyproject.toml index 28087a7b..098696af 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,7 +8,7 @@ build-backend = "scikit_build_core.build" [project] name = "simcoon" -version = "1.10.2" +version = "1.11.0" description = "Simulation in Mechanics and Materials: Interactive Tools - A library for the simulation of multiphysics systems and heterogeneous materials" readme = "README.md" license = {text = "GPL-3.0-or-later"} @@ -187,7 +187,7 @@ before-all = [ archs = ["arm64"] # delocate bundles shared libraries -repair-wheel-command = "delocate-wheel --require-archs {delocate_archs} -w {dest_dir} -v {wheel}" +repair-wheel-command = "delocate-wheel --require-archs {delocate_archs} -w {dest_dir} -v --exclude libomp.dylib {wheel}" # Add macOS-specific env without replacing global environment (ARMADILLO_VERSION, etc.) [[tool.cibuildwheel.overrides]] diff --git a/python-setup/simcoon/rotation.py b/python-setup/simcoon/rotation.py index 83572195..cdc6cc7e 100644 --- a/python-setup/simcoon/rotation.py +++ b/python-setup/simcoon/rotation.py @@ -142,135 +142,212 @@ def from_scipy(cls, scipy_rot): return cls.from_quat(scipy_rot.as_quat()) # ------------------------------------------------------------------ - # Internal helper + # Internal helpers # ------------------------------------------------------------------ + @property + def _is_batch(self): + """True if this object stores more than one rotation.""" + return self.as_quat().ndim == 2 + def _to_cpp(self): - """Convert to a _CppRotation for C++ method dispatch.""" - return _CppRotation.from_quat(self.as_quat()) + """Convert to a _CppRotation for C++ method dispatch (single only).""" + q = self.as_quat() + if q.ndim == 2: + raise ValueError( + "Cannot convert a batch Rotation to a single _CppRotation. " + "Batch operations are handled automatically by the Python methods." + ) + return _CppRotation.from_quat(q) + + def _voigt_stress_matrices(self, active=True): + """Return QS matrices: (6,6) for single, (N,6,6) for batch.""" + q = self.as_quat() + if q.ndim == 1: + return _CppRotation.from_quat(q).as_voigt_stress_rotation(active) + return np.array([ + _CppRotation.from_quat(q[i]).as_voigt_stress_rotation(active) + for i in range(len(q)) + ]) + + def _voigt_strain_matrices(self, active=True): + """Return QE matrices: (6,6) for single, (N,6,6) for batch.""" + q = self.as_quat() + if q.ndim == 1: + return _CppRotation.from_quat(q).as_voigt_strain_rotation(active) + return np.array([ + _CppRotation.from_quat(q[i]).as_voigt_strain_rotation(active) + for i in range(len(q)) + ]) # ------------------------------------------------------------------ - # Mechanics methods (delegate to _CppRotation) + # Mechanics methods — support single and batch (Gauss-point) operations + # + # Single rotation: + # sigma (6,) → (6,) + # L (6,6) → (6,6) + # + # Batch of N rotations: + # sigma (6, N) → (6, N) one stress per rotation + # L (6, 6, N) → (6, 6, N) one stiffness per rotation # ------------------------------------------------------------------ def apply_stress(self, sigma, active=True): - """Apply rotation to a stress vector in Voigt notation. + """Apply rotation to stress vector(s) in Voigt notation. Parameters ---------- sigma : array_like - 6-component stress vector [s11, s22, s33, s12, s13, s23]. + Single (6,) stress vector, or (6, N) array for batch + (one column per Gauss point). active : bool, optional If True (default), active rotation. Returns ------- numpy.ndarray - Rotated stress vector. + Rotated stress: (6,) or (6, N). """ - return self._to_cpp().apply_stress(np.asarray(sigma, dtype=float), active).ravel() + sigma = np.asarray(sigma, dtype=float) + if not self._is_batch: + return self._to_cpp().apply_stress(sigma.ravel(), active).ravel() + QS = self._voigt_stress_matrices(active) # (N, 6, 6) + return np.einsum("nij,jn->in", QS, sigma) def apply_strain(self, epsilon, active=True): - """Apply rotation to a strain vector in Voigt notation. + """Apply rotation to strain vector(s) in Voigt notation. Parameters ---------- epsilon : array_like - 6-component strain vector [e11, e22, e33, 2*e12, 2*e13, 2*e23]. + Single (6,) strain vector, or (6, N) for batch. active : bool, optional If True (default), active rotation. Returns ------- numpy.ndarray - Rotated strain vector. + Rotated strain: (6,) or (6, N). """ - return self._to_cpp().apply_strain(np.asarray(epsilon, dtype=float), active).ravel() + epsilon = np.asarray(epsilon, dtype=float) + if not self._is_batch: + return self._to_cpp().apply_strain(epsilon.ravel(), active).ravel() + QE = self._voigt_strain_matrices(active) # (N, 6, 6) + return np.einsum("nij,jn->in", QE, epsilon) def apply_stiffness(self, L, active=True): - """Apply rotation to a 6x6 stiffness matrix. + """Apply rotation to 6x6 stiffness matrix/matrices. Parameters ---------- L : array_like - 6x6 stiffness matrix in Voigt notation. + Single (6, 6) stiffness matrix, or (6, 6, N) for batch. active : bool, optional If True (default), active rotation. Returns ------- numpy.ndarray - Rotated 6x6 stiffness matrix. + Rotated stiffness: (6, 6) or (6, 6, N). """ - return self._to_cpp().apply_stiffness(np.asarray(L, dtype=float), active) + L = np.asarray(L, dtype=float) + if not self._is_batch: + return self._to_cpp().apply_stiffness(L, active) + QS = self._voigt_stress_matrices(active) # (N, 6, 6) + # L_rot = QS @ L @ QS^T for each n + return np.einsum("nij,jkn,nlk->iln", QS, L, QS) def apply_compliance(self, M, active=True): - """Apply rotation to a 6x6 compliance matrix. + """Apply rotation to 6x6 compliance matrix/matrices. Parameters ---------- M : array_like - 6x6 compliance matrix in Voigt notation. + Single (6, 6) compliance matrix, or (6, 6, N) for batch. active : bool, optional If True (default), active rotation. Returns ------- numpy.ndarray - Rotated 6x6 compliance matrix. + Rotated compliance: (6, 6) or (6, 6, N). """ - return self._to_cpp().apply_compliance(np.asarray(M, dtype=float), active) + M = np.asarray(M, dtype=float) + if not self._is_batch: + return self._to_cpp().apply_compliance(M, active) + QE = self._voigt_strain_matrices(active) # (N, 6, 6) + # M_rot = QE @ M @ QE^T for each n + return np.einsum("nij,jkn,nlk->iln", QE, M, QE) def apply_strain_concentration(self, A, active=True): - """Apply rotation to a 6x6 strain concentration tensor. + """Apply rotation to 6x6 strain concentration tensor(s). Parameters ---------- A : array_like - 6x6 strain concentration tensor in Voigt notation. + Single (6, 6) tensor, or (6, 6, N) for batch. active : bool, optional If True (default), active rotation. Returns ------- numpy.ndarray - Rotated strain concentration tensor: QE * A * QS^T. + Rotated tensor: QE * A * QS^T. Shape (6, 6) or (6, 6, N). """ - return self._to_cpp().apply_strain_concentration(np.asarray(A, dtype=float), active) + A = np.asarray(A, dtype=float) + if not self._is_batch: + return self._to_cpp().apply_strain_concentration(A, active) + QE = self._voigt_strain_matrices(active) # (N, 6, 6) + QS = self._voigt_stress_matrices(active) # (N, 6, 6) + return np.einsum("nij,jkn,nlk->iln", QE, A, QS) def apply_stress_concentration(self, B, active=True): - """Apply rotation to a 6x6 stress concentration tensor. + """Apply rotation to 6x6 stress concentration tensor(s). Parameters ---------- B : array_like - 6x6 stress concentration tensor in Voigt notation. + Single (6, 6) tensor, or (6, 6, N) for batch. active : bool, optional If True (default), active rotation. Returns ------- numpy.ndarray - Rotated stress concentration tensor: QS * B * QE^T. + Rotated tensor: QS * B * QE^T. Shape (6, 6) or (6, 6, N). """ - return self._to_cpp().apply_stress_concentration(np.asarray(B, dtype=float), active) + B = np.asarray(B, dtype=float) + if not self._is_batch: + return self._to_cpp().apply_stress_concentration(B, active) + QS = self._voigt_stress_matrices(active) # (N, 6, 6) + QE = self._voigt_strain_matrices(active) # (N, 6, 6) + return np.einsum("nij,jkn,nlk->iln", QS, B, QE) def apply_tensor(self, m, inverse=False): - """Apply rotation to a 3x3 tensor (matrix). + """Apply rotation to 3x3 tensor(s). Parameters ---------- m : array_like - 3x3 tensor to rotate. + Single (3, 3) tensor, or (3, 3, N) for batch. inverse : bool, optional If True, apply inverse rotation. Default is False. Returns ------- numpy.ndarray - Rotated tensor: R * m * R^T (or R^T * m * R for inverse). + Rotated tensor: (3, 3) or (3, 3, N). """ - return self._to_cpp().apply_tensor(np.asarray(m, dtype=float), inverse) + m = np.asarray(m, dtype=float) + if not self._is_batch: + return self._to_cpp().apply_tensor(m, inverse) + # R matrices: (N, 3, 3) + R = self.as_matrix() + if inverse: + # R^T @ m @ R for each n + return np.einsum("nji,jkn,nkl->iln", R, m, R) + # R @ m @ R^T for each n + return np.einsum("nij,jkn,nlk->iln", R, m, R) def as_voigt_stress_rotation(self, active=True): """Get 6x6 rotation matrix for stress tensors in Voigt notation. @@ -283,9 +360,9 @@ def as_voigt_stress_rotation(self, active=True): Returns ------- numpy.ndarray - 6x6 stress rotation matrix (QS). + Single (6, 6) or batch (N, 6, 6) stress rotation matrix (QS). """ - return self._to_cpp().as_voigt_stress_rotation(active) + return self._voigt_stress_matrices(active) def as_voigt_strain_rotation(self, active=True): """Get 6x6 rotation matrix for strain tensors in Voigt notation. @@ -298,9 +375,9 @@ def as_voigt_strain_rotation(self, active=True): Returns ------- numpy.ndarray - 6x6 strain rotation matrix (QE). + Single (6, 6) or batch (N, 6, 6) strain rotation matrix (QE). """ - return self._to_cpp().as_voigt_strain_rotation(active) + return self._voigt_strain_matrices(active) # ------------------------------------------------------------------ # Compatibility helpers diff --git a/simcoon-python-builder/CMakeLists.txt b/simcoon-python-builder/CMakeLists.txt index a75cf0f9..9d211e64 100755 --- a/simcoon-python-builder/CMakeLists.txt +++ b/simcoon-python-builder/CMakeLists.txt @@ -74,8 +74,12 @@ if(CONDA_BUILD) else() # Wheel: libsimcoon is in same directory as _core if(APPLE) + # @loader_path → find libsimcoon.dylib next to _core.so + # @loader_path/../../.. → $CONDA_PREFIX/lib (conda's libomp.dylib) + # /usr/local/lib → CI / manual installs + # /opt/homebrew/opt/libomp/lib → Homebrew on Apple Silicon set_target_properties(_core PROPERTIES - INSTALL_RPATH "@loader_path" + INSTALL_RPATH "@loader_path;@loader_path/../../..;/usr/local/lib;/opt/homebrew/opt/libomp/lib" BUILD_WITH_INSTALL_RPATH ON ) elseif(UNIX) diff --git a/simcoon-python-builder/test/test_core/test_rotation.py b/simcoon-python-builder/test/test_core/test_rotation.py index a9f43a50..2771e9d9 100644 --- a/simcoon-python-builder/test/test_core/test_rotation.py +++ b/simcoon-python-builder/test/test_core/test_rotation.py @@ -291,3 +291,100 @@ def test_from_axis_angle_degrees(self): def test_from_axis_angle_invalid_axis(self): with pytest.raises(ValueError, match="axis must be 1, 2, or 3"): sim.Rotation.from_axis_angle(0.5, 4) + + +# =================================================================== +# 6. Batch (Gauss-point) operations +# =================================================================== + +class TestBatch: + """Batch rotation operations on (6, N) / (6, 6, N) arrays, + matching the Gauss-point data layout used by simcoon and fedoo.""" + + @pytest.fixture + def n_gauss(self): + return 64 + + @pytest.fixture + def batch_rot(self, n_gauss): + return sim.Rotation.random(n_gauss) + + def test_batch_apply_stress(self, batch_rot, n_gauss): + sigma = np.random.rand(6, n_gauss) * 100 + result = batch_rot.apply_stress(sigma) + assert result.shape == (6, n_gauss) + for i in range(n_gauss): + expected = batch_rot[i].apply_stress(sigma[:, i]) + np.testing.assert_allclose(result[:, i], expected, atol=1e-10) + + def test_batch_apply_strain(self, batch_rot, n_gauss): + eps = np.random.rand(6, n_gauss) * 0.01 + result = batch_rot.apply_strain(eps) + assert result.shape == (6, n_gauss) + for i in range(n_gauss): + expected = batch_rot[i].apply_strain(eps[:, i]) + np.testing.assert_allclose(result[:, i], expected, atol=1e-10) + + def test_batch_apply_stiffness(self, batch_rot, n_gauss): + L_single = sim.L_iso([70000, 0.3], "Enu") + L = np.repeat(L_single[:, :, np.newaxis], n_gauss, axis=2) + result = batch_rot.apply_stiffness(L) + assert result.shape == (6, 6, n_gauss) + for i in range(n_gauss): + expected = batch_rot[i].apply_stiffness(L_single) + np.testing.assert_allclose(result[:, :, i], expected, atol=1e-10) + + def test_batch_apply_compliance(self, batch_rot, n_gauss): + M_single = sim.M_iso([70000, 0.3], "Enu") + M = np.repeat(M_single[:, :, np.newaxis], n_gauss, axis=2) + result = batch_rot.apply_compliance(M) + assert result.shape == (6, 6, n_gauss) + for i in range(n_gauss): + expected = batch_rot[i].apply_compliance(M_single) + np.testing.assert_allclose(result[:, :, i], expected, atol=1e-10) + + def test_batch_apply_tensor(self, batch_rot, n_gauss): + T_single = np.array([[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]]) + T = np.repeat(T_single[:, :, np.newaxis], n_gauss, axis=2) + result = batch_rot.apply_tensor(T) + assert result.shape == (3, 3, n_gauss) + for i in range(n_gauss): + expected = batch_rot[i].apply_tensor(T_single) + np.testing.assert_allclose(result[:, :, i], expected, atol=1e-10) + + def test_batch_voigt_stress_rotation(self, batch_rot, n_gauss): + QS = batch_rot.as_voigt_stress_rotation() + assert QS.shape == (n_gauss, 6, 6) + for i in range(n_gauss): + expected = batch_rot[i].as_voigt_stress_rotation() + np.testing.assert_allclose(QS[i], expected, atol=1e-10) + + def test_batch_isotropic_invariance(self, batch_rot, n_gauss): + """Isotropic stiffness must be invariant under any rotation.""" + L_single = sim.L_iso([200e3, 0.3], "Enu") + L = np.repeat(L_single[:, :, np.newaxis], n_gauss, axis=2) + L_rot = batch_rot.apply_stiffness(L) + for i in range(n_gauss): + np.testing.assert_allclose(L_rot[:, :, i], L_single, atol=1e-6) + + def test_batch_roundtrip_stress(self, batch_rot, n_gauss): + sigma = np.random.rand(6, n_gauss) * 100 + sigma_rot = batch_rot.apply_stress(sigma) + sigma_back = batch_rot.apply_stress(sigma_rot, active=False) + np.testing.assert_allclose(sigma_back, sigma, atol=1e-10) + + def test_batch_from_matrix_dr_cube(self, n_gauss): + """Create batch Rotation from (3, 3, N) DR cube (fedoo layout).""" + rotations = sim.Rotation.random(n_gauss) + DR = rotations.as_matrix().transpose(1, 2, 0) # (3, 3, N) + assert DR.shape == (3, 3, n_gauss) + + rotations2 = sim.Rotation.from_matrix(DR.transpose(2, 0, 1)) + assert len(rotations2) == n_gauss + + sigma = np.random.rand(6, n_gauss) * 100 + np.testing.assert_allclose( + rotations.apply_stress(sigma), + rotations2.apply_stress(sigma), + atol=1e-10, + ) diff --git a/tests/test_fedoo_integration.py b/tests/test_fedoo_integration.py new file mode 100644 index 00000000..2c086250 --- /dev/null +++ b/tests/test_fedoo_integration.py @@ -0,0 +1,130 @@ +"""Integration tests for fedoo + simcoon. + +Tests both non-linear umat (OpenMP) and non-linear geometry +(tangent modulus transfer via deformation gradient). +""" + +import numpy as np +import pytest + +import fedoo +import simcoon +from simcoon import simmit as sim + + +@pytest.fixture(autouse=True) +def _reset_fedoo(): + """Reset fedoo global state between tests.""" + fedoo.ModelingSpace("3D") + yield + + +def _make_cube(name, nx=2): + """Create a unit cube mesh with standard Dirichlet BCs on min faces.""" + return fedoo.mesh.box_mesh( + nx=nx, ny=nx, nz=nx, + x_min=0, x_max=1, y_min=0, y_max=1, z_min=0, z_max=1, + elm_type="hex8", name=name, + ) + + +def _apply_compression_bc(pb, mesh, disp_z): + """Fix min faces, prescribe DispZ on top face.""" + pb.bc.add("Dirichlet", mesh.find_nodes("X", mesh.bounding_box.xmin), "DispX", 0) + pb.bc.add("Dirichlet", mesh.find_nodes("Y", mesh.bounding_box.ymin), "DispY", 0) + pb.bc.add("Dirichlet", mesh.find_nodes("Z", mesh.bounding_box.zmin), "DispZ", 0) + pb.bc.add("Dirichlet", mesh.find_nodes("Z", mesh.bounding_box.zmax), "DispZ", disp_z) + + +def test_epicp_small_strain(): + """Elastoplastic EPICP umat, small strain (exercises OpenMP).""" + mesh = _make_cube("epicp_ss") + props = np.array([200e3, 0.3, 0.0, 300.0, 1000.0, 0.5]) + mat = fedoo.constitutivelaw.Simcoon("EPICP", props, name="epicp_ss_mat") + + wf = fedoo.weakform.StressEquilibrium(mat, name="wf_epicp_ss") + asm = fedoo.Assembly.create(wf, mesh, elm_type="hex8", name="asm_epicp_ss") + + pb = fedoo.problem.NonLinear(asm, name="pb_epicp_ss") + _apply_compression_bc(pb, mesh, disp_z=0.01) + + pb.nlsolve(dt=0.5, tmax=1.0, update_dt=True, print_info=0) + + disp = pb.get_disp("DispZ") + assert np.max(np.abs(disp)) == pytest.approx(0.01, abs=1e-8) + p = asm.sv["Statev"][1] + assert np.max(p) > 0, "Expected plastic deformation" + + +def test_neohookean_nlgeom(): + """Neo-Hookean hyperelastic, large strain (tangent modulus transfer).""" + mesh = _make_cube("neohc_nl") + props = np.array([80e3, 400e3]) # mu, kappa + mat = fedoo.constitutivelaw.Simcoon("NEOHC", props, name="neohc_nl_mat") + + wf = fedoo.weakform.StressEquilibrium(mat, name="wf_neohc_nl", nlgeom=True) + asm = fedoo.Assembly.create(wf, mesh, elm_type="hex8", name="asm_neohc_nl") + + pb = fedoo.problem.NonLinear(asm, nlgeom=True, name="pb_neohc_nl") + _apply_compression_bc(pb, mesh, disp_z=0.1) + + pb.nlsolve(dt=0.5, tmax=1.0, update_dt=True, print_info=0) + + disp = pb.get_disp("DispZ") + assert np.max(np.abs(disp)) == pytest.approx(0.1, abs=1e-8) + F = asm.sv["F"] + assert not np.allclose(F, np.eye(3).reshape(3, 3, 1)) + + +def test_epicp_nlgeom(): + """EPICP + large strain (OpenMP umat + tangent modulus + F transfer).""" + mesh = _make_cube("epicp_nl") + props = np.array([200e3, 0.3, 0.0, 300.0, 1000.0, 0.5]) + mat = fedoo.constitutivelaw.Simcoon("EPICP", props, name="epicp_nl_mat") + + wf = fedoo.weakform.StressEquilibrium(mat, name="wf_epicp_nl", nlgeom=True) + asm = fedoo.Assembly.create(wf, mesh, elm_type="hex8", name="asm_epicp_nl") + + pb = fedoo.problem.NonLinear(asm, nlgeom=True, name="pb_epicp_nl") + _apply_compression_bc(pb, mesh, disp_z=0.05) + + pb.nlsolve(dt=0.25, tmax=1.0, update_dt=True, print_info=0) + + disp = pb.get_disp("DispZ") + assert np.max(np.abs(disp)) == pytest.approx(0.05, abs=1e-8) + p = asm.sv["Statev"][1] + assert np.max(p) > 0, "Expected plastic deformation" + F = asm.sv["F"] + assert not np.allclose(F, np.eye(3).reshape(3, 3, 1)) + Lt = asm.sv["TangentMatrix"] + sym_err = np.abs(Lt - Lt.transpose((1, 0, 2))).max() + assert sym_err < 1e-3, f"Tangent matrix not symmetric (err={sym_err})" + + +@pytest.mark.parametrize("n_threads", [1, 2, 4]) +def test_umat_openmp_threads(n_threads): + """Direct sim.umat with explicit n_threads to verify OpenMP.""" + n_gauss = 64 + props = np.asfortranarray( + np.tile([200e3, 0.3, 0.0, 300.0, 1000.0, 0.5], (n_gauss, 1)).T + ) + etot = np.zeros((6, n_gauss), order="F") + Detot = np.zeros((6, n_gauss), order="F") + Detot[2, :] = 5e-3 + + sigma = np.zeros((6, n_gauss), order="F") + statev = np.zeros((8, n_gauss), order="F") + Wm = np.zeros((4, n_gauss), order="F") + DR = np.empty((3, 3, n_gauss), order="F") + DR[...] = np.eye(3).reshape(3, 3, 1) + + stress, sv, wm, Lt = sim.umat( + "EPICP", etot, Detot, + np.array([]), np.array([]), + sigma, DR, props, statev, + 0.0, 1.0, Wm, n_threads=n_threads, + ) + + assert np.allclose(stress[:, 0:1], stress, atol=1e-6), \ + f"Stress not uniform across Gauss points with n_threads={n_threads}" + assert sv[1, 0] > 0, f"Expected plastic strain with n_threads={n_threads}"