diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt index 85447e27cdc..e62a5419031 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/CMakeLists.txt @@ -62,7 +62,6 @@ set(HEADER_FILES ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/IdentityRotation.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/PolarDecomposition.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/RotationMethodsContainer.h - ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/StablePolarDecomposition.h ${SOFACOMPONENTSOLIDMECHANICSFEMELASTIC_SOURCE_DIR}/impl/rotations/TriangleRotation.h ) diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h index beaaf5995ef..4f610a84d83 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.h @@ -27,7 +27,6 @@ #include #include #include -#include #include #include #include @@ -54,10 +53,13 @@ namespace sofa::component::solidmechanics::fem::elastic */ template struct RotationMethods : RotationMethodsContainer, PolarDecomposition, IdentityRotation + RigidFitStablePolarDecomposition, RigidFitPolarDecomposition, IdentityRotation, + AffineFitStablePolarDecomposition, AffineFitPolarDecomposition > { - using Inherit = RotationMethodsContainer, PolarDecomposition, IdentityRotation>; + using Inherit = RotationMethodsContainer, RigidFitPolarDecomposition, IdentityRotation, + AffineFitStablePolarDecomposition, AffineFitPolarDecomposition>; explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) {} }; @@ -68,11 +70,13 @@ using Vec3Real = sofa::defaulttype::StdVectorTypes, sof //partial specialization for linear triangle in 3D template struct RotationMethods, sofa::geometry::Triangle> : RotationMethodsContainer, sofa::geometry::Triangle, - StablePolarDecomposition>, PolarDecomposition>, IdentityRotation, TriangleRotation> + RigidFitStablePolarDecomposition, sofa::geometry::Triangle>, RigidFitPolarDecomposition, sofa::geometry::Triangle>, IdentityRotation, TriangleRotation>, + AffineFitStablePolarDecomposition, sofa::geometry::Triangle>, AffineFitPolarDecomposition, sofa::geometry::Triangle> > { using Inherit = RotationMethodsContainer, sofa::geometry::Triangle, - StablePolarDecomposition>, PolarDecomposition>, IdentityRotation, TriangleRotation> >; + RigidFitStablePolarDecomposition, sofa::geometry::Triangle>, RigidFitPolarDecomposition, sofa::geometry::Triangle>, IdentityRotation, TriangleRotation>, + AffineFitStablePolarDecomposition, sofa::geometry::Triangle>, AffineFitPolarDecomposition, sofa::geometry::Triangle>>; explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) { @@ -83,11 +87,13 @@ struct RotationMethods, sofa::geometry::Triangle> : RotationMetho //partial specialization for linear tetrahedron template struct RotationMethods : RotationMethodsContainer, PolarDecomposition, IdentityRotation, TriangleRotation + RigidFitStablePolarDecomposition, RigidFitPolarDecomposition, IdentityRotation, TriangleRotation, + AffineFitStablePolarDecomposition, AffineFitPolarDecomposition > { using Inherit = RotationMethodsContainer, PolarDecomposition, IdentityRotation, TriangleRotation >; + RigidFitStablePolarDecomposition, RigidFitPolarDecomposition, IdentityRotation, TriangleRotation, + AffineFitStablePolarDecomposition, AffineFitPolarDecomposition>; explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) { @@ -98,11 +104,13 @@ struct RotationMethods : RotationMethods //partial specialization for linear hexahedron template struct RotationMethods : RotationMethodsContainer, PolarDecomposition, IdentityRotation, HexahedronRotation + RigidFitStablePolarDecomposition, RigidFitPolarDecomposition, IdentityRotation, HexahedronRotation, + AffineFitStablePolarDecomposition, AffineFitPolarDecomposition > { using Inherit = RotationMethodsContainer, PolarDecomposition, IdentityRotation, HexahedronRotation >; + RigidFitStablePolarDecomposition, RigidFitPolarDecomposition, IdentityRotation, HexahedronRotation, + AffineFitStablePolarDecomposition, AffineFitPolarDecomposition>; explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) { @@ -113,15 +121,17 @@ struct RotationMethods : RotationMethodsC //partial specialization for quadratic tetrahedron template struct RotationMethods : RotationMethodsContainer, IdentityRotation + RigidFitPolarDecomposition, IdentityRotation, + AffineFitPolarDecomposition > { using Inherit = RotationMethodsContainer, IdentityRotation>; + RigidFitPolarDecomposition, IdentityRotation, + AffineFitPolarDecomposition>; explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) { - this->d_rotationMethod.setValue(PolarDecomposition::getItem().key); + this->d_rotationMethod.setValue(RigidFitPolarDecomposition::getItem().key); } }; diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h index fc96cd6db28..f9b413d911a 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/PolarDecomposition.h @@ -22,17 +22,31 @@ #pragma once #include +#include #include #include namespace sofa::component::solidmechanics::fem::elastic { -template +/// How the motion of the nodes is fitted before the polar decomposition. +enum class Fit { Rigid, Affine }; + +/// How the polar decomposition is computed: +/// Standard calls Decompose::polarDecomposition +/// Stable calls Decompose::polarDecomposition_stable +enum class Algorithm { Standard, Stable }; + +template struct PolarDecomposition { using RotationMatrix = sofa::type::Mat>; + /// The remainder of the dimensions of an embedded element within a higher-dimensional space. + // TODO: embedded elements (Codimension > 0) are not handled. + static constexpr sofa::Size Codimension = + DataTypes::spatial_dimensions - sofa::fem::FiniteElement::TopologicalDimension; + template sofa::Coord_t computeCentroid(const std::array, NumberOfNodesInElement>& nodes) { @@ -64,15 +78,73 @@ struct PolarDecomposition Q[j] = nodesRestPosition[j] - t0; } - const auto H = P.multTranspose(Q); + const auto F = [&P, &Q]() -> RotationMatrix + { + const auto H = P.multTranspose(Q); + + if constexpr (Codimension > 0) + { + // Not handled: H is decomposed as it is. + return H; + } + else if constexpr (FitType == Fit::Affine) + { + const auto M = Q.multTranspose(Q); + return H * (M / sofa::type::trace(M)).inverted(); + } + else + { + return H; + } + }(); - sofa::helper::Decompose>::polarDecomposition(H, rotationMatrix); + if constexpr (AlgorithmType == Algorithm::Standard) + { + sofa::helper::Decompose>::polarDecomposition(F, rotationMatrix); + } + else + { + sofa::helper::Decompose>::polarDecomposition_stable(F, rotationMatrix); + } } static constexpr sofa::helper::Item getItem() { - return {"polar", "Polar decomposition"}; + if constexpr (AlgorithmType == Algorithm::Standard) + { + if constexpr (FitType == Fit::Rigid) + { + return {"polar", "Polar decomposition of the rigid fit"}; + } + else + { + return {"polar_affine", "Polar decomposition of the affine fit"}; + } + } + else + { + if constexpr (FitType == Fit::Rigid) + { + return {"stable_polar", "Stable polar decomposition of the rigid fit"}; + } + else + { + return {"stable_polar_affine", "Stable polar decomposition of the affine fit"}; + } + } } }; +template +using RigidFitPolarDecomposition = PolarDecomposition; + +template +using AffineFitPolarDecomposition = PolarDecomposition; + +template +using RigidFitStablePolarDecomposition = PolarDecomposition; + +template +using AffineFitStablePolarDecomposition = PolarDecomposition; + } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/StablePolarDecomposition.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/StablePolarDecomposition.h deleted file mode 100644 index 96f84d4569f..00000000000 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/impl/rotations/StablePolarDecomposition.h +++ /dev/null @@ -1,78 +0,0 @@ -/****************************************************************************** -* SOFA, Simulation Open-Framework Architecture * -* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * -* * -* This program is free software; you can redistribute it and/or modify it * -* under the terms of the GNU Lesser General Public License as published by * -* the Free Software Foundation; either version 2.1 of the License, or (at * -* your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, but WITHOUT * -* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * -* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * -* for more details. * -* * -* You should have received a copy of the GNU Lesser General Public License * -* along with this program. If not, see . * -******************************************************************************* -* Authors: The SOFA Team and external contributors (see Authors.txt) * -* * -* Contact information: contact@sofa-framework.org * -******************************************************************************/ -#pragma once - -#include -#include -#include - -namespace sofa::component::solidmechanics::fem::elastic -{ - -template -struct StablePolarDecomposition -{ - using RotationMatrix = sofa::type::Mat>; - - template - sofa::Coord_t computeCentroid(const std::array, NumberOfNodesInElement>& nodes) - { - sofa::Coord_t centroid; - for (const auto& node : nodes) - { - centroid += node; - } - centroid /= static_cast>(NumberOfNodesInElement); - return centroid; - } - - template - void computeRotation(RotationMatrix& rotationMatrix, const RotationMatrix& initialRotationMatrix, - const std::array, NumberOfNodesInElement>& nodesPosition, - const std::array, NumberOfNodesInElement>& nodesRestPosition) - { - SOFA_UNUSED(initialRotationMatrix); - - const auto t = computeCentroid(nodesPosition); - const auto t0 = computeCentroid(nodesRestPosition); - - sofa::type::Mat> P(sofa::type::NOINIT); - sofa::type::Mat> Q(sofa::type::NOINIT); - - for (sofa::Size j = 0; j < NumberOfNodesInElement; ++j) - { - P[j] = nodesPosition[j] - t; - Q[j] = nodesRestPosition[j] - t0; - } - - const auto H = P.multTranspose(Q); - - sofa::helper::Decompose>::polarDecomposition_stable(H, rotationMatrix); - } - - static constexpr sofa::helper::Item getItem() - { - return {"stable_polar", "Stable polar decomposition"}; - } -}; - -} diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt index 756b079090b..3c2c042bae7 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/CMakeLists.txt @@ -9,6 +9,7 @@ set(SOURCE_FILES HexahedralFEMForceField_test.cpp HexahedronFEMForceField_test.cpp LinearSmallStrainFEMForceField_test.cpp + RotationMethods_test.cpp StrainDisplacement_test.cpp TetrahedralCorotationalFEMForceField_test.cpp TetrahedronFEMForceField_stepTest.cpp diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/RotationMethods_test.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/RotationMethods_test.cpp new file mode 100644 index 00000000000..59ba5091d2d --- /dev/null +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/RotationMethods_test.cpp @@ -0,0 +1,554 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace sofa +{ + +using namespace sofa::component::solidmechanics::fem::elastic; + +namespace +{ + +using sofa::defaulttype::Vec2Types; +using sofa::defaulttype::Vec3Types; + +template +using Matrix = sofa::type::Mat>; + +template +using Nodes = std::array, N>; + +constexpr SReal degreeToRadian = std::numbers::pi_v / 180; + +Matrix rotation(SReal angle) +{ + return {{std::cos(angle), -std::sin(angle)}, {std::sin(angle), std::cos(angle)}}; +} + +/// Rotation by an angle about an axis. +Matrix rotation(const sofa::type::Vec3& axis, SReal angle) +{ + const auto K = sofa::type::crossProductMatrix(axis.normalized()); + return Matrix::Identity() + K * std::sin(angle) + K * K * (1 - std::cos(angle)); +} + +// Shapes + +/// Maps the plane z = 0 to a plane that contains no coordinate axis. +const Matrix tilt = rotation(sofa::type::Vec3(1.0, -1.0, 2.0), 40 * degreeToRadian); + +/// The 2D element placed in the plane z = 0. +template +Nodes embedElement(const Nodes& nodes) +{ + Nodes embeddedNodes; + for (std::size_t j = 0; j < N; ++j) + { + embeddedNodes[j] = sofa::type::Vec3(nodes[j][0], nodes[j][1], 0.0); + } + return embeddedNodes; +} + +/// The 2D element placed in a tilted plane. +template +Nodes tiltElement(const Nodes& nodes) +{ + Nodes tiltedNodes = embedElement(nodes); + for (auto& node : tiltedNodes) + { + node = tilt * node; + } + return tiltedNodes; +} + +// Regular shapes: M = sum_j q_j q_j^T is proportional to the identity. + +const Nodes regularQuad2D { + sofa::type::Vec2(0.0, 0.0), sofa::type::Vec2(1.0, 0.0), + sofa::type::Vec2(1.0, 1.0), sofa::type::Vec2(0.0, 1.0)}; +const Nodes regularTriangle2D { + sofa::type::Vec2(0.0, 0.0), sofa::type::Vec2(1.0, 0.0), sofa::type::Vec2(0.5, std::sqrt(3.0) / 2)}; +const Nodes regularTriangleFlat = embedElement(regularTriangle2D); +const Nodes regularTriangleTilted = tiltElement(regularTriangle2D); +const Nodes regularTetrahedron { + sofa::type::Vec3(1.0, 1.0, 1.0), sofa::type::Vec3(1.0, -1.0, -1.0), + sofa::type::Vec3(-1.0, 1.0, -1.0), sofa::type::Vec3(-1.0, -1.0, 1.0)}; +const Nodes regularHexahedron { + sofa::type::Vec3(0.0, 0.0, 0.0), sofa::type::Vec3(1.0, 0.0, 0.0), + sofa::type::Vec3(1.0, 1.0, 0.0), sofa::type::Vec3(0.0, 1.0, 0.0), + sofa::type::Vec3(0.0, 0.0, 1.0), sofa::type::Vec3(1.0, 0.0, 1.0), + sofa::type::Vec3(1.0, 1.0, 1.0), sofa::type::Vec3(0.0, 1.0, 1.0)}; + +// Irregular shapes: M is not proportional to the identity. + +const Nodes irregularQuad2D { + sofa::type::Vec2(0.0, 0.0), sofa::type::Vec2(3.0, 0.0), + sofa::type::Vec2(3.0, 1.0), sofa::type::Vec2(0.0, 1.0)}; +const Nodes irregularTriangle2D { + sofa::type::Vec2(0.0, 0.0), sofa::type::Vec2(1.0, 0.0), sofa::type::Vec2(0.0, 1.0)}; +const Nodes irregularTriangleFlat = embedElement(irregularTriangle2D); +const Nodes irregularTriangleTilted = tiltElement(irregularTriangle2D); +const Nodes irregularTetrahedron { + sofa::type::Vec3(0.0, 0.0, 0.0), sofa::type::Vec3(1.0, 0.0, 0.0), + sofa::type::Vec3(0.0, 1.0, 0.0), sofa::type::Vec3(0.0, 0.0, 1.0)}; +const Nodes irregularHexahedron { + sofa::type::Vec3(0.0, 0.0, 0.0), sofa::type::Vec3(3.0, 0.0, 0.0), + sofa::type::Vec3(3.0, 2.0, 0.0), sofa::type::Vec3(0.0, 2.0, 0.0), + sofa::type::Vec3(0.0, 0.0, 1.0), sofa::type::Vec3(3.0, 0.0, 1.0), + sofa::type::Vec3(3.0, 2.0, 1.0), sofa::type::Vec3(0.0, 2.0, 1.0)}; + +// Rigid motion + +const Matrix rotation2D = rotation(30 * degreeToRadian); +const Matrix rotation3D = rotation(sofa::type::Vec3(1.0, 2.0, 3.0), 30 * degreeToRadian); + +// Stretch patch: stretches along a direction that follows no element edge + +const Matrix stretch2D {{1.5, 0.0}, {0.0, 1.0}}; +const Matrix stretch3D {{1.5, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; + +const Matrix offAxisStretch2D = rotation2D * stretch2D * rotation2D.transposed(); +const Matrix offAxisStretch3D = rotation3D * stretch3D * rotation3D.transposed(); + +/// offAxisStretch2D within the plane z = 0. +const Matrix flatStretch3D = [] +{ + Matrix stretch = Matrix::Identity(); + for (std::size_t i = 0; i < 2; ++i) + { + for (std::size_t j = 0; j < 2; ++j) + { + stretch(i, j) = offAxisStretch2D(i, j); + } + } + return stretch; +}(); + +/// offAxisStretch2D within the tilted plane. +const Matrix tiltedStretch3D = tilt * flatStretch3D * tilt.transposed(); + +// Renumbering: the stretches followed by a rotation, and renumberings that keep the element valid + +const Matrix rotatedStretch2D = rotation2D * offAxisStretch2D; +const Matrix rotatedStretch3D = rotation3D * offAxisStretch3D; +const Matrix rotatedFlatStretch3D = rotation3D * flatStretch3D; +const Matrix rotatedTiltedStretch3D = rotation3D * tiltedStretch3D; + +const std::array triangleNumbering {1, 2, 0}; +const std::array quadNumbering {1, 2, 3, 0}; +const std::array tetrahedronNumbering {1, 2, 0, 3}; +const std::array hexahedronNumbering {1, 2, 3, 0, 5, 6, 7, 4}; + +/// Whether R equals the expected rotation, entry by entry, within 1e-10. +template +::testing::AssertionResult sameRotation(const RotationMatrix& R, const RotationMatrix& expected) +{ + for (std::size_t i = 0; i < RotationMatrix::nbLines; ++i) + { + for (std::size_t j = 0; j < RotationMatrix::nbCols; ++j) + { + if (std::abs(R(i, j) - expected(i, j)) > 1e-10) + { + return ::testing::AssertionFailure() + << "R(" << i << ", " << j << ") = " << R(i, j) << ", expected " << expected(i, j); + } + } + } + return ::testing::AssertionSuccess() << "R equals the expected rotation"; +} + +/// Whether R is orthogonal with determinant +1, within 1e-10. +template +::testing::AssertionResult isProperRotation(const RotationMatrix& R) +{ + const auto RtR = R.multTranspose(R); + for (std::size_t i = 0; i < RotationMatrix::nbLines; ++i) + { + for (std::size_t j = 0; j < RotationMatrix::nbCols; ++j) + { + if (std::abs(RtR(i, j) - (i == j ? 1 : 0)) > 1e-10) + { + return ::testing::AssertionFailure() + << "R is not orthogonal: (R^T R)(" << i << ", " << j << ") = " << RtR(i, j); + } + } + } + if (std::abs(sofa::type::determinant(R) - 1) > 1e-10) + { + return ::testing::AssertionFailure() << "R is not a rotation: det R = " << sofa::type::determinant(R); + } + return ::testing::AssertionSuccess(); +} + +/// Rotation from the rest to the current positions, expected to be a proper rotation. Some methods +/// need the initial rotation, which is computed first from the rest positions. +template +typename Method::RotationMatrix extractRotation(const std::array& positions, + const std::array& restPositions) +{ + using RotationMatrix = typename Method::RotationMatrix; + + Method method; + RotationMatrix initialRotation; + method.computeRotation(initialRotation, RotationMatrix::Identity(), restPositions, restPositions); + initialRotation.transpose(); + + RotationMatrix R; + method.computeRotation(R, initialRotation, positions, restPositions); + EXPECT_TRUE(isProperRotation(R)); + return R; +} + +/// A rigid motion x = R0 X + c gives R = R0. +template +::testing::AssertionResult meetsRigidMotion(const std::array& restPositions, + const typename Method::RotationMatrix& R0) +{ + Coord c; + for (std::size_t i = 0; i < Coord::spatial_dimensions; ++i) + { + c[i] = 0.5 * (i + 1); + } + + std::array positions; + for (std::size_t j = 0; j < N; ++j) + { + positions[j] = R0 * restPositions[j] + c; + } + + return sameRotation(extractRotation(positions, restPositions), R0); +} + +/// A uniform stretch x = U X gives R = I. +template +::testing::AssertionResult meetsStretchPatch(const std::array& restPositions, + const typename Method::RotationMatrix& U) +{ + std::array positions; + for (std::size_t j = 0; j < N; ++j) + { + positions[j] = U * restPositions[j]; + } + + return sameRotation(extractRotation(positions, restPositions), Method::RotationMatrix::Identity()); +} + +/// Renumbering the nodes of a deformed element x = F X leaves R unchanged. +template +::testing::AssertionResult meetsRenumbering(const std::array& restPositions, + const typename Method::RotationMatrix& F, + const std::array& numbering) +{ + std::array positions; + for (std::size_t j = 0; j < N; ++j) + { + positions[j] = F * restPositions[j]; + } + + std::array renumberedPositions; + std::array renumberedRestPositions; + for (std::size_t j = 0; j < N; ++j) + { + renumberedPositions[j] = positions[numbering[j]]; + renumberedRestPositions[j] = restPositions[numbering[j]]; + } + + return sameRotation(extractRotation(renumberedPositions, renumberedRestPositions), + extractRotation(positions, restPositions)); +} + +} + +TEST(RotationMethods, regularQuad2DRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(regularQuad2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(regularQuad2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(regularQuad2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(regularQuad2D, rotation2D))); +} + +TEST(RotationMethods, regularTriangle2DRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(regularTriangle2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(regularTriangle2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(regularTriangle2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(regularTriangle2D, rotation2D))); +} + +TEST(RotationMethods, regularTriangleFlatRigidMotion) +{ + EXPECT_TRUE(meetsRigidMotion>(regularTriangleFlat, rotation3D)); +} + +TEST(RotationMethods, regularTriangleTiltedRigidMotion) +{ + EXPECT_TRUE(meetsRigidMotion>(regularTriangleTilted, rotation3D)); +} + +TEST(RotationMethods, regularTetrahedronRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(regularTetrahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(regularTetrahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(regularTetrahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(regularTetrahedron, rotation3D))); + EXPECT_TRUE(meetsRigidMotion>(regularTetrahedron, rotation3D)); +} + +TEST(RotationMethods, regularHexahedronRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(regularHexahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(regularHexahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(regularHexahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(regularHexahedron, rotation3D))); + EXPECT_TRUE(meetsRigidMotion>(regularHexahedron, rotation3D)); +} + +TEST(RotationMethods, irregularQuad2DRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(irregularQuad2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(irregularQuad2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(irregularQuad2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(irregularQuad2D, rotation2D))); +} + +TEST(RotationMethods, irregularTriangle2DRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(irregularTriangle2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(irregularTriangle2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(irregularTriangle2D, rotation2D))); + EXPECT_TRUE((meetsRigidMotion>(irregularTriangle2D, rotation2D))); +} + +TEST(RotationMethods, irregularTriangleFlatRigidMotion) +{ + EXPECT_TRUE(meetsRigidMotion>(irregularTriangleFlat, rotation3D)); +} + +TEST(RotationMethods, irregularTriangleTiltedRigidMotion) +{ + EXPECT_TRUE(meetsRigidMotion>(irregularTriangleTilted, rotation3D)); +} + +TEST(RotationMethods, irregularTetrahedronRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(irregularTetrahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(irregularTetrahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(irregularTetrahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(irregularTetrahedron, rotation3D))); + EXPECT_TRUE(meetsRigidMotion>(irregularTetrahedron, rotation3D)); +} + +TEST(RotationMethods, irregularHexahedronRigidMotion) +{ + EXPECT_TRUE((meetsRigidMotion>(irregularHexahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(irregularHexahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(irregularHexahedron, rotation3D))); + EXPECT_TRUE((meetsRigidMotion>(irregularHexahedron, rotation3D))); + EXPECT_TRUE(meetsRigidMotion>(irregularHexahedron, rotation3D)); +} + +TEST(RotationMethods, regularQuad2DStretchPatch) +{ + EXPECT_TRUE((meetsStretchPatch>(regularQuad2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(regularQuad2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(regularQuad2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(regularQuad2D, offAxisStretch2D))); +} + +TEST(RotationMethods, regularTriangle2DStretchPatch) +{ + EXPECT_TRUE((meetsStretchPatch>(regularTriangle2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(regularTriangle2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(regularTriangle2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(regularTriangle2D, offAxisStretch2D))); +} + +TEST(RotationMethods, regularTriangleFlatStretchPatch) +{ + EXPECT_FALSE(meetsStretchPatch>(regularTriangleFlat, flatStretch3D)); +} + +TEST(RotationMethods, regularTriangleTiltedStretchPatch) +{ + EXPECT_FALSE(meetsStretchPatch>(regularTriangleTilted, tiltedStretch3D)); +} + +TEST(RotationMethods, regularTetrahedronStretchPatch) +{ + EXPECT_TRUE((meetsStretchPatch>(regularTetrahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(regularTetrahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(regularTetrahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(regularTetrahedron, offAxisStretch3D))); + EXPECT_FALSE(meetsStretchPatch>(regularTetrahedron, offAxisStretch3D)); +} + +TEST(RotationMethods, regularHexahedronStretchPatch) +{ + EXPECT_TRUE((meetsStretchPatch>(regularHexahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(regularHexahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(regularHexahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(regularHexahedron, offAxisStretch3D))); + EXPECT_FALSE(meetsStretchPatch>(regularHexahedron, offAxisStretch3D)); +} + +TEST(RotationMethods, irregularQuad2DStretchPatch) +{ + EXPECT_FALSE((meetsStretchPatch>(irregularQuad2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(irregularQuad2D, offAxisStretch2D))); + EXPECT_FALSE((meetsStretchPatch>(irregularQuad2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(irregularQuad2D, offAxisStretch2D))); +} + +TEST(RotationMethods, irregularTriangle2DStretchPatch) +{ + EXPECT_FALSE((meetsStretchPatch>(irregularTriangle2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(irregularTriangle2D, offAxisStretch2D))); + EXPECT_FALSE((meetsStretchPatch>(irregularTriangle2D, offAxisStretch2D))); + EXPECT_TRUE((meetsStretchPatch>(irregularTriangle2D, offAxisStretch2D))); +} + +TEST(RotationMethods, irregularTriangleFlatStretchPatch) +{ + EXPECT_FALSE(meetsStretchPatch>(irregularTriangleFlat, flatStretch3D)); +} + +TEST(RotationMethods, irregularTriangleTiltedStretchPatch) +{ + EXPECT_FALSE(meetsStretchPatch>(irregularTriangleTilted, tiltedStretch3D)); +} + +TEST(RotationMethods, irregularTetrahedronStretchPatch) +{ + EXPECT_FALSE((meetsStretchPatch>(irregularTetrahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(irregularTetrahedron, offAxisStretch3D))); + EXPECT_FALSE((meetsStretchPatch>(irregularTetrahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(irregularTetrahedron, offAxisStretch3D))); + EXPECT_FALSE(meetsStretchPatch>(irregularTetrahedron, offAxisStretch3D)); +} + +TEST(RotationMethods, irregularHexahedronStretchPatch) +{ + EXPECT_FALSE((meetsStretchPatch>(irregularHexahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(irregularHexahedron, offAxisStretch3D))); + EXPECT_FALSE((meetsStretchPatch>(irregularHexahedron, offAxisStretch3D))); + EXPECT_TRUE((meetsStretchPatch>(irregularHexahedron, offAxisStretch3D))); + EXPECT_FALSE(meetsStretchPatch>(irregularHexahedron, offAxisStretch3D)); +} + +TEST(RotationMethods, regularQuad2DRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(regularQuad2D, rotatedStretch2D, quadNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularQuad2D, rotatedStretch2D, quadNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularQuad2D, rotatedStretch2D, quadNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularQuad2D, rotatedStretch2D, quadNumbering))); +} + +TEST(RotationMethods, regularTriangle2DRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(regularTriangle2D, rotatedStretch2D, triangleNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularTriangle2D, rotatedStretch2D, triangleNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularTriangle2D, rotatedStretch2D, triangleNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularTriangle2D, rotatedStretch2D, triangleNumbering))); +} + +TEST(RotationMethods, regularTriangleFlatRenumbering) +{ + EXPECT_FALSE(meetsRenumbering>(regularTriangleFlat, rotatedFlatStretch3D, triangleNumbering)); +} + +TEST(RotationMethods, regularTriangleTiltedRenumbering) +{ + EXPECT_FALSE(meetsRenumbering>(regularTriangleTilted, rotatedTiltedStretch3D, triangleNumbering)); +} + +TEST(RotationMethods, regularTetrahedronRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(regularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_FALSE(meetsRenumbering>(regularTetrahedron, rotatedStretch3D, tetrahedronNumbering)); +} + +TEST(RotationMethods, regularHexahedronRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(regularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(regularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_FALSE(meetsRenumbering>(regularHexahedron, rotatedStretch3D, hexahedronNumbering)); +} + +TEST(RotationMethods, irregularQuad2DRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(irregularQuad2D, rotatedStretch2D, quadNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularQuad2D, rotatedStretch2D, quadNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularQuad2D, rotatedStretch2D, quadNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularQuad2D, rotatedStretch2D, quadNumbering))); +} + +TEST(RotationMethods, irregularTriangle2DRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(irregularTriangle2D, rotatedStretch2D, triangleNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularTriangle2D, rotatedStretch2D, triangleNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularTriangle2D, rotatedStretch2D, triangleNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularTriangle2D, rotatedStretch2D, triangleNumbering))); +} + +TEST(RotationMethods, irregularTriangleFlatRenumbering) +{ + EXPECT_FALSE(meetsRenumbering>(irregularTriangleFlat, rotatedFlatStretch3D, triangleNumbering)); +} + +TEST(RotationMethods, irregularTriangleTiltedRenumbering) +{ + EXPECT_FALSE(meetsRenumbering>(irregularTriangleTilted, rotatedTiltedStretch3D, triangleNumbering)); +} + +TEST(RotationMethods, irregularTetrahedronRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(irregularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularTetrahedron, rotatedStretch3D, tetrahedronNumbering))); + EXPECT_FALSE(meetsRenumbering>(irregularTetrahedron, rotatedStretch3D, tetrahedronNumbering)); +} + +TEST(RotationMethods, irregularHexahedronRenumbering) +{ + EXPECT_TRUE((meetsRenumbering>(irregularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_TRUE((meetsRenumbering>(irregularHexahedron, rotatedStretch3D, hexahedronNumbering))); + EXPECT_FALSE(meetsRenumbering>(irregularHexahedron, rotatedStretch3D, hexahedronNumbering)); +} + +}