diff --git a/Sofa/Component/Engine/Generate/CMakeLists.txt b/Sofa/Component/Engine/Generate/CMakeLists.txt index 62df2348889..e7fecbb054b 100644 --- a/Sofa/Component/Engine/Generate/CMakeLists.txt +++ b/Sofa/Component/Engine/Generate/CMakeLists.txt @@ -24,6 +24,8 @@ set(HEADER_FILES ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/GroupFilterYoungModulus.inl ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/JoinPoints.h ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/JoinPoints.inl + ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/LinearToHigherOrderElements.h + ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/LinearToHigherOrderElements.inl ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/MergeMeshes.h ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/MergeMeshes.inl ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/MergePoints.h @@ -62,6 +64,7 @@ set(SOURCE_FILES ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/GenerateSphere.cpp ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/GroupFilterYoungModulus.cpp ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/JoinPoints.cpp + ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/LinearToHigherOrderElements.cpp ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/MergeMeshes.cpp ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/MergePoints.cpp ${SOFACOMPONENTENGINEGENERATE_SOURCE_DIR}/MergeSets.cpp diff --git a/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.cpp b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.cpp new file mode 100644 index 00000000000..eca84271cab --- /dev/null +++ b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.cpp @@ -0,0 +1,35 @@ +/****************************************************************************** +* 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 + +namespace sofa::component::engine::generate +{ + +void registerLinearToHigherOrderElements(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(core::ObjectRegistrationData("Merge several meshes.") + .add< LinearToHigherOrderElements >(true) + ); +} + +} diff --git a/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.h b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.h new file mode 100644 index 00000000000..a5056c2087a --- /dev/null +++ b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.h @@ -0,0 +1,71 @@ +/****************************************************************************** +* 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 +#include + +namespace sofa::component::engine::generate +{ + +template +class LinearToHigherOrderElements : + public sofa::core::DataEngine, + public sofa::core::behavior::TopologyAccessor, + public sofa::core::behavior::SingleStateAccessor +{ +public: + SOFA_CLASS3(LinearToHigherOrderElements, + sofa::core::DataEngine, + sofa::core::behavior::TopologyAccessor, + sofa::core::behavior::SingleStateAccessor); + + template //e.g. sofa::geometry::Edge + using TopologyElement = sofa::topology::Element; + + template //e.g. sofa::geometry::Edge + using SeqElement = sofa::type::vector>; + + void init() override; + + sofa::DataVecCoord_t d_position; + + Data> d_quadraticEdges; + Data> d_quadraticTriangles; + Data> d_quadraticQuads; + Data> d_quadraticTetrahedra; + Data> d_quadraticHexahedra; + + Data d_computeFromEdges; + Data d_computeFromTriangles; + Data d_computeFromQuads; + Data d_computeFromTetrahedra; + Data d_computeFromHexahedra; + +protected: + void doUpdate() override; + + LinearToHigherOrderElements(); +}; + +} // namespace sofa::component::engine::generate diff --git a/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.inl b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.inl new file mode 100644 index 00000000000..bc7596535dc --- /dev/null +++ b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/LinearToHigherOrderElements.inl @@ -0,0 +1,226 @@ +/****************************************************************************** +* 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 + +namespace sofa::component::engine::generate +{ + +template +LinearToHigherOrderElements::LinearToHigherOrderElements() + : d_position(initData(&d_position, sofa::VecCoord_t{}, "position", "Output position")) + , d_quadraticEdges(initData(&d_quadraticEdges, SeqElement{}, + "quadratic_edges", "List of quadratic edges")) + , d_quadraticTriangles(initData(&d_quadraticTriangles, SeqElement{}, + "quadratic_triangles", "List of quadratic triangles")) + , d_quadraticQuads(initData(&d_quadraticQuads, SeqElement{}, + "quadratic_quads", "List of quadratic quads")) + , d_quadraticTetrahedra(initData(&d_quadraticTetrahedra, SeqElement{}, + "quadratic_tetrahedra", "List of quadratic tetrahedra")) + , d_quadraticHexahedra(initData(&d_quadraticHexahedra, SeqElement{}, + "quadratic_hexahedra", "List of quadratic hexahedra")) + , d_computeFromEdges(initData(&d_computeFromEdges, true, "computeFromEdges", + "Use edges from the input to generate higher-order edges")) + , d_computeFromTriangles(initData(&d_computeFromTriangles, true, "computeFromTriangles", + "Use triangles from the input to generate higher-order triangles")) + , d_computeFromQuads(initData(&d_computeFromQuads, true, "computeFromQuads", + "Use quads from the input to generate higher-order quads")) + , d_computeFromTetrahedra(initData(&d_computeFromTetrahedra, true, "computeFromTetrahedra", + "Use tetrahedra from the input to generate higher-order tetrahedra")) + , d_computeFromHexahedra(initData(&d_computeFromHexahedra, true, "computeFromHexahedra", + "Use hexahedra from the input to generate higher-order hexahedra")) +{ + this->addOutput(&d_position); + this->addOutput(&d_quadraticEdges); + this->addOutput(&d_quadraticTriangles); + this->addOutput(&d_quadraticQuads); + this->addOutput(&d_quadraticTetrahedra); + this->addOutput(&d_quadraticHexahedra); +} + +template +void LinearToHigherOrderElements::init() +{ + DataEngine::init(); + + if (!this->isComponentStateInvalid()) + { + this->validateTopology(); + } + + if (!this->isComponentStateInvalid()) + { + sofa::core::behavior::SingleStateAccessor::init(); + } +} + + +struct TetrahedronEdge +{ + TetrahedronEdge(sofa::Index p, sofa::Index q) : a(std::min(p, q)), b(std::max(p, q)) {} + sofa::Index a, b; + bool operator==(const TetrahedronEdge& e) const + { + return a == e.a && b == e.b; + } +}; + +struct Face { + std::array indices; + Face(sofa::Index a, sofa::Index b, sofa::Index c, sofa::Index d) { + indices = {a, b, c, d}; + std::sort(indices.begin(), indices.end()); + } + bool operator==(const Face& other) const { return indices == other.indices; } +}; + + +template +void LinearToHigherOrderElements::doUpdate() +{ + if (!l_topology) + return; + + //input elements + const auto& edges = l_topology->getElements(); + const auto& triangles = l_topology->getElements(); + const auto& quads = l_topology->getElements(); + const auto& tetrahedra = l_topology->getElements(); + const auto& hexahedra = l_topology->getElements(); + + //input position + const auto inPosition = this->mstate->readPositions(); + + //output elements + auto quadraticEdges = sofa::helper::getWriteOnlyAccessor(d_quadraticEdges); + auto quadraticTriangles = sofa::helper::getWriteOnlyAccessor(d_quadraticTriangles); + auto quadraticQuads = sofa::helper::getWriteOnlyAccessor(d_quadraticQuads); + auto quadraticTetrahedra = sofa::helper::getWriteOnlyAccessor(d_quadraticTetrahedra); + auto quadraticHexahedra = sofa::helper::getWriteOnlyAccessor(d_quadraticHexahedra); + + //output position + auto outPosition = sofa::helper::getWriteOnlyAccessor(d_position); + outPosition.wref() = inPosition.ref(); + + auto edgeHash = [&inPosition](const TetrahedronEdge& edge){ return edge.a + inPosition.size() * edge.b; }; + std::unordered_map midEdgePointsMap(10, edgeHash); + + auto getOrAddMidEdgePoint = [&](sofa::Index a, sofa::Index b) + { + TetrahedronEdge edge{a, b}; + const auto [it, success] = midEdgePointsMap.insert(std::make_pair(edge, outPosition.size())); + if (success) + { + outPosition.push_back(0.5 * (inPosition[a] + inPosition[b])); + } + return it->second; + }; + + auto faceHash = [](const Face& f) { + size_t h = 0; + for (auto i : f.indices) h ^= std::hash{}(i) + 0x9e3779b9 + (h << 6) + (h >> 2); + return h; + }; + std::unordered_map midFacePointsMap(10, faceHash); + + auto getOrAddMidFacePoint = [&](sofa::Index a, sofa::Index b, sofa::Index c, sofa::Index d) + { + Face face{a, b, c, d}; + const auto [it, success] = midFacePointsMap.insert(std::make_pair(face, outPosition.size())); + if (success) + { + outPosition.push_back(0.25 * (inPosition[a] + inPosition[b] + inPosition[c] + inPosition[d])); + } + return it->second; + }; + + if (d_computeFromEdges.getValue()) + for (const auto& element : edges) + { + quadraticEdges.emplace_back(element[0], element[1], getOrAddMidEdgePoint(element[0], element[1])); + } + + if (d_computeFromTriangles.getValue()) + for (const auto& element : triangles) + { + quadraticTriangles.emplace_back(element[0], element[1], element[2], + getOrAddMidEdgePoint(element[0], element[1]), + getOrAddMidEdgePoint(element[1], element[2]), + getOrAddMidEdgePoint(element[2], element[0])); + } + + if (d_computeFromQuads.getValue()) + for (const auto& element : quads) + { + quadraticQuads.emplace_back(element[0], element[1], element[2], element[3], + getOrAddMidEdgePoint(element[0], element[1]), + getOrAddMidEdgePoint(element[1], element[2]), + getOrAddMidEdgePoint(element[2], element[3]), + getOrAddMidEdgePoint(element[3], element[0]), + getOrAddMidFacePoint(element[0], element[1], element[2], element[3])); + } + + if (d_computeFromTetrahedra.getValue()) + for (const auto& element : tetrahedra) + { + std::array newIndices; + for (std::size_t i = 0; i < 6; ++i) + { + const auto [a, b] = core::topology::edgesInTetrahedronArray[i]; + newIndices[i] = getOrAddMidEdgePoint(element[a], element[b]); + } + + quadraticTetrahedra.emplace_back(element[0], element[1], element[2], element[3], + newIndices[0], newIndices[1], newIndices[2], newIndices[3], + newIndices[4], newIndices[5]); + } + + if (d_computeFromHexahedra.getValue()) + for (const auto& element : hexahedra) + { + std::array midEdges; + for (std::size_t i = 0; i < 12; ++i) + { + midEdges[i] = + getOrAddMidEdgePoint(element[core::topology::edgesInHexahedronArray[i][0]], + element[core::topology::edgesInHexahedronArray[i][1]]); + } + + std::array midFaces; + for (std::size_t i = 0; i < 6; ++i) + { + const auto quad = core::topology::quadsOrientationInHexahedronArray[i]; + midFaces[i] = getOrAddMidFacePoint(element[quad[0]], element[quad[1]], element[quad[2]], element[quad[3]]); + } + + sofa::Index midVolume = outPosition.size(); + outPosition.push_back(0.125 * (inPosition[element[0]] + inPosition[element[1]] + inPosition[element[2]] + inPosition[element[3]] + + inPosition[element[4]] + inPosition[element[5]] + inPosition[element[6]] + inPosition[element[7]])); + + quadraticHexahedra.emplace_back(element[0], element[1], element[2], element[3], element[4], element[5], element[6], element[7], + midEdges[0], midEdges[1], midEdges[2], midEdges[3], midEdges[4], midEdges[5], midEdges[6], midEdges[7], midEdges[8], midEdges[9], midEdges[10], midEdges[11], + midFaces[0], midFaces[1], midFaces[2], midFaces[3], midFaces[4], midFaces[5], midVolume); + } +} + +} // namespace sofa::component::engine::generate diff --git a/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/init.cpp b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/init.cpp index ebdca5b8b3c..be36cc6d2c1 100644 --- a/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/init.cpp +++ b/Sofa/Component/Engine/Generate/src/sofa/component/engine/generate/init.cpp @@ -35,6 +35,7 @@ extern void registerGenerateRigidMass(sofa::core::ObjectFactory* factory); extern void registerGenerateSphere(sofa::core::ObjectFactory* factory); extern void registerGroupFilterYoungModulus(sofa::core::ObjectFactory* factory); extern void registerJoinPoints(sofa::core::ObjectFactory* factory); +extern void registerLinearToHigherOrderElements(sofa::core::ObjectFactory* factory); extern void registerMergeMeshes(sofa::core::ObjectFactory* factory); extern void registerMergePoints(sofa::core::ObjectFactory* factory); extern void registerMergeSets(sofa::core::ObjectFactory* factory); @@ -42,12 +43,12 @@ extern void registerMergeVectors(sofa::core::ObjectFactory* factory); extern void registerMeshBarycentricMapperEngine(sofa::core::ObjectFactory* factory); extern void registerMeshClosingEngine(sofa::core::ObjectFactory* factory); extern void registerMeshTetraStuffing(sofa::core::ObjectFactory* factory); -extern void registerNormalsFromPoints(sofa::core::ObjectFactory* factory); extern void registerNormEngine(sofa::core::ObjectFactory* factory); +extern void registerNormalsFromPoints(sofa::core::ObjectFactory* factory); extern void registerRandomPointDistributionInSurface(sofa::core::ObjectFactory* factory); extern void registerSpiral(sofa::core::ObjectFactory* factory); -extern void registerVolumeFromTriangles(sofa::core::ObjectFactory* factory); extern void registerVolumeFromTetrahedrons(sofa::core::ObjectFactory* factory); +extern void registerVolumeFromTriangles(sofa::core::ObjectFactory* factory); extern "C" { SOFA_EXPORT_DYNAMIC_LIBRARY void initExternalModule(); @@ -82,6 +83,7 @@ void registerObjects(sofa::core::ObjectFactory* factory) registerGenerateSphere(factory); registerGroupFilterYoungModulus(factory); registerJoinPoints(factory); + registerLinearToHigherOrderElements(factory); registerMergeMeshes(factory); registerMergePoints(factory); registerMergeSets(factory); diff --git a/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.cpp b/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.cpp index fc194d76fde..9e52da4fa7f 100644 --- a/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.cpp +++ b/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.cpp @@ -40,6 +40,18 @@ template class SOFA_COMPONENT_MASS_API FEMMass; template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; + void registerFEMMass(sofa::core::ObjectFactory* factory) { factory->registerObjects(sofa::core::ObjectRegistrationData("Finite-element mass (inertia and body force)") @@ -54,6 +66,18 @@ void registerFEMMass(sofa::core::ObjectFactory* factory) .add< FEMMass >() .add< FEMMass >() .add< FEMMass >() + + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() + .add< FEMMass >() ); } diff --git a/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.h b/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.h index 46eab79274a..d132c121b8d 100644 --- a/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.h +++ b/Sofa/Component/Mass/src/sofa/component/mass/FEMMass.h @@ -261,6 +261,18 @@ template class SOFA_COMPONENT_MASS_API FEMMass; template class SOFA_COMPONENT_MASS_API FEMMass; template class SOFA_COMPONENT_MASS_API FEMMass; + +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; +template class SOFA_COMPONENT_MASS_API FEMMass; #endif } // namespace sofa::component::mass diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.cpp index da1a5d7dcc3..4a1e84d42d8 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/CorotationalFEMForceField.cpp @@ -43,6 +43,17 @@ void registerCorotationalFEMForceField(sofa::core::ObjectFactory* factory) .add< CorotationalFEMForceField >() .add< CorotationalFEMForceField >() .add< CorotationalFEMForceField >() + + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() + .add< CorotationalFEMForceField >() ); } @@ -58,4 +69,15 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForc template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; + } 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 a74f6d0b002..32a553f2ac8 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 @@ -109,6 +109,21 @@ struct RotationMethods : RotationMethodsC } }; +//partial specialization for quadratic tetrahedron +template +struct RotationMethods : RotationMethodsContainer, IdentityRotation +> +{ + using Inherit = RotationMethodsContainer, IdentityRotation>; + + explicit RotationMethods(sofa::core::objectmodel::BaseObject* parent) : Inherit(parent) + { + this->d_rotationMethod.setValue(PolarDecomposition::getItem().key); + } +}; + template class CorotationalFEMForceField : @@ -186,6 +201,17 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API Corotational extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API CorotationalFEMForceField; #endif } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.cpp index e25229b29eb..abe930ca8c6 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.cpp @@ -43,6 +43,18 @@ void registerLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* factory) .add< LinearSmallStrainFEMForceField >() .add< LinearSmallStrainFEMForceField >() .add< LinearSmallStrainFEMForceField >() + + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() + .add< LinearSmallStrainFEMForceField >() ); } @@ -58,4 +70,15 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFE template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h index 899edf4c084..f4b3d37b306 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/LinearSmallStrainFEMForceField.h @@ -93,6 +93,18 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallS extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; + +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API LinearSmallStrainFEMForceField; #endif } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/TriangleFEMForceField_stepTest.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/TriangleFEMForceField_stepTest.cpp index c29bf739f08..aa0ae7e23b1 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/tests/TriangleFEMForceField_stepTest.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/tests/TriangleFEMForceField_stepTest.cpp @@ -45,7 +45,7 @@ struct TriangleFEMForceField_stepTest : public ForceField_test<_TriangleFEMForce { auto topology = modeling::addNew(this->node); topology->setName("topology"); - topology->d_seqTriangles.setValue({{0,1,2}}); + topology->d_seqTriangles.toData().setValue({{0,1,2}}); topology->d_seqPoints.setParent(&this->dof->x); //Position diff --git a/Sofa/Component/Topology/Container/Constant/CMakeLists.txt b/Sofa/Component/Topology/Container/Constant/CMakeLists.txt index 958c2faaa08..ebf509cff7d 100644 --- a/Sofa/Component/Topology/Container/Constant/CMakeLists.txt +++ b/Sofa/Component/Topology/Container/Constant/CMakeLists.txt @@ -7,6 +7,7 @@ set(HEADER_FILES ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/config.h.in ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/init.h ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/MeshTopology.h + ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/MeshTopologyContainer.h ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/CubeTopology.h ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/SphereQuadTopology.h ) @@ -14,6 +15,7 @@ set(HEADER_FILES set(SOURCE_FILES ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/init.cpp ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/MeshTopology.cpp + ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/MeshTopologyContainer.cpp ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/CubeTopology.cpp ${SOFACOMPONENTTOPOLOGYCONTAINERCONSTANT_SOURCE_DIR}/SphereQuadTopology.cpp ) diff --git a/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.cpp b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.cpp index b0e14ac3a76..a4beb496ea7 100644 --- a/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.cpp +++ b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.cpp @@ -39,16 +39,16 @@ MeshTopology::EdgeUpdate::EdgeUpdate(MeshTopology* t) { if( topology->hasVolume() ) { - addInput(&topology->d_seqHexahedra); - addInput(&topology->d_seqTetrahedra); - addOutput(&topology->d_seqEdges); + addInput(&topology->d_seqHexahedra.toData()); + addInput(&topology->d_seqTetrahedra.toData()); + addOutput(&topology->d_seqEdges.toData()); setDirtyValue(); } else if( topology->hasSurface() ) { - addInput(&topology->d_seqTriangles); - addInput(&topology->d_seqQuads); - addOutput(&topology->d_seqEdges); + addInput(&topology->d_seqTriangles.toData()); + addInput(&topology->d_seqQuads.toData()); + addOutput(&topology->d_seqEdges.toData()); setDirtyValue(); } @@ -221,8 +221,8 @@ void MeshTopology::EdgeUpdate::updateFromSurface() MeshTopology::TriangleUpdate::TriangleUpdate(MeshTopology *t) :PrimitiveUpdate(t) { - addInput(&topology->d_seqTetrahedra); - addOutput(&topology->d_seqTriangles); + addInput(&topology->d_seqTetrahedra.toData()); + addOutput(&topology->d_seqTriangles.toData()); setDirtyValue(); } @@ -275,8 +275,8 @@ void MeshTopology::TriangleUpdate::doUpdate() MeshTopology::QuadUpdate::QuadUpdate(MeshTopology *t) :PrimitiveUpdate(t) { - addInput(&topology->d_seqHexahedra); - addOutput(&topology->d_seqQuads); + addInput(&topology->d_seqHexahedra.toData()); + addOutput(&topology->d_seqQuads.toData()); setDirtyValue(); } @@ -486,13 +486,6 @@ void registerMeshTopology(sofa::core::ObjectFactory* factory) MeshTopology::MeshTopology() : d_seqPoints(initData(&d_seqPoints, "position", "List of point positions")) - , d_seqEdges(initData(&d_seqEdges, "edges", "List of edge indices")) - , d_seqTriangles(initData(&d_seqTriangles, "triangles", "List of triangle indices")) - , d_seqQuads(initData(&d_seqQuads, "quads", "List of quad indices")) - , d_seqTetrahedra(initData(&d_seqTetrahedra, "tetrahedra", "List of tetrahedron indices")) - , d_seqHexahedra(initData(&d_seqHexahedra, "hexahedra", "List of hexahedron indices")) - , d_seqPrisms(initData(&d_seqPrisms, "prisms", "List of prisms indices")) - , d_seqPyramids(initData(&d_seqPyramids, "pyramids", "List of pyramids indices")) , d_seqUVs(initData(&d_seqUVs, "uv", "List of uv coordinates")) , d_computeAllBuffers(initData(&d_computeAllBuffers, false, "computeAllBuffers", "Option to compute all crossed topology buffers at init. False by default")) , nbPoints(0) @@ -506,9 +499,9 @@ MeshTopology::MeshTopology() { m_upperElementType = sofa::geometry::ElementType::EDGE; addAlias(&d_seqPoints, "points"); - addAlias(&d_seqEdges, "lines"); - addAlias(&d_seqTetrahedra, "tetras"); - addAlias(&d_seqHexahedra, "hexas"); + addAlias(&d_seqEdges.toData(), "lines"); + addAlias(&d_seqTetrahedra.toData(), "tetras"); + addAlias(&d_seqHexahedra.toData(), "hexas"); addAlias(&d_seqUVs, "texcoords"); } @@ -747,51 +740,6 @@ void MeshTopology::addUV(SReal u, SReal v) nbPoints = (int)d_seqUVs.getValue().size(); } -const MeshTopology::SeqEdges& MeshTopology::getEdges() -{ - return d_seqEdges.getValue(); -} - -const MeshTopology::SeqTriangles& MeshTopology::getTriangles() -{ - return d_seqTriangles.getValue(); -} - -const MeshTopology::SeqQuads& MeshTopology::getQuads() -{ - return d_seqQuads.getValue(); -} - -const MeshTopology::SeqTetrahedra& MeshTopology::getTetrahedra() -{ - if (!validTetrahedra) - { - updateTetrahedra(); - validTetrahedra = true; - } - return d_seqTetrahedra.getValue(); -} - -const MeshTopology::SeqHexahedra& MeshTopology::getHexahedra() -{ - if (!validHexahedra) - { - updateHexahedra(); - validHexahedra = true; - } - return d_seqHexahedra.getValue(); -} - -const BaseMeshTopology::SeqPrisms& MeshTopology::getPrisms() -{ - return d_seqPrisms.getValue(); -} - -const BaseMeshTopology::SeqPyramids& MeshTopology::getPyramids() -{ - return d_seqPyramids.getValue(); -} - const MeshTopology::SeqUV& MeshTopology::getUVs() { return d_seqUVs.getValue(); diff --git a/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.h b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.h index 281f5ceeca1..ed2dd459cee 100644 --- a/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.h +++ b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopology.h @@ -23,17 +23,17 @@ #include #include -#include +#include #include #include namespace sofa::component::topology::container::constant { -class SOFA_COMPONENT_TOPOLOGY_CONTAINER_CONSTANT_API MeshTopology : public core::topology::BaseMeshTopology +class SOFA_COMPONENT_TOPOLOGY_CONTAINER_CONSTANT_API MeshTopology : public MeshTopologyContainer { public: - SOFA_CLASS(MeshTopology,core::topology::BaseMeshTopology); + SOFA_CLASS(MeshTopology, MeshTopologyContainer); protected: class PrimitiveUpdate : public sofa::core::DataEngine @@ -91,16 +91,6 @@ class SOFA_COMPONENT_TOPOLOGY_CONTAINER_CONSTANT_API MeshTopology : public core: void setNbPoints(Size n) override; - // Complete sequence accessors - - const SeqEdges& getEdges() override; - const SeqTriangles& getTriangles() override; - const SeqQuads& getQuads() override; - const SeqTetrahedra& getTetrahedra() override; - const SeqHexahedra& getHexahedra() override; - const SeqPrisms& getPrisms() override; - const SeqPyramids& getPyramids() override; - // If using STEP loader, include also uv coordinates typedef Index UVID; typedef type::Vec2 UV; @@ -297,17 +287,19 @@ class SOFA_COMPONENT_TOPOLOGY_CONTAINER_CONSTANT_API MeshTopology : public core: typedef type::vector SeqPoints; Data< SeqPoints > d_seqPoints; ///< List of point positions - Data d_seqEdges; ///< List of edge indices - Data d_seqTriangles; ///< List of triangle indices - Data d_seqQuads; ///< List of quad indices - Data d_seqTetrahedra; ///< List of tetrahedron indices - Data d_seqHexahedra; ///< List of hexahedron indices - Data d_seqPrisms; - Data d_seqPyramids; + SeqElementProxy d_seqEdges { this }; + SeqElementProxy d_seqTriangles { this }; + SeqElementProxy d_seqQuads { this }; + SeqElementProxy d_seqTetrahedra { this }; + SeqElementProxy d_seqHexahedra { this }; + SeqElementProxy d_seqPrisms { this }; + SeqElementProxy d_seqPyramids { this }; + Data d_seqUVs; ///< List of uv coordinates Data d_computeAllBuffers; ///< Option to call method computeCrossElementBuffers. False by default protected: + Size nbPoints; template diff --git a/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopologyContainer.cpp b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopologyContainer.cpp new file mode 100644 index 00000000000..fc62f37f5bf --- /dev/null +++ b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopologyContainer.cpp @@ -0,0 +1,163 @@ +/****************************************************************************** +* 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 + +namespace sofa::component::topology::container::constant +{ + +using sofa::geometry::ElementType; + +constexpr const char* elementTypeToDataName(ElementType type) +{ + switch (type) + { + case ElementType::POINT: { return "position"; } + case ElementType::EDGE: { return "edges"; } + case ElementType::TRIANGLE: { return "triangles"; } + case ElementType::QUAD: { return "quads"; } + case ElementType::TETRAHEDRON: { return "tetrahedra"; } + case ElementType::HEXAHEDRON: { return "hexahedra"; } + case ElementType::PRISM: { return "prisms"; } + case ElementType::PYRAMID: { return "pyramids"; } + + case ElementType::QUADRATIC_EDGE: { return "quadratic_edges"; } + case ElementType::QUADRATIC_TRIANGLE: { return "quadratic_triangles"; } + case ElementType::QUADRATIC_QUAD: { return "quadratic_quads"; } + case ElementType::QUADRATIC_TETRAHEDRON: { return "quadratic_tetrahedra"; } + case ElementType::QUADRATIC_HEXAHEDRON: { return "quadratic_hexahedra"; } + case ElementType::QUADRATIC_PRISM: { return "quadratic_prisms"; } + case ElementType::QUADRATIC_PYRAMID: { return "quadratic_pyramids"; } + + default: + return "Unknown"; + } +} +template +void registerSingleData(core::objectmodel::BaseComponent* meshTopology, auto& dataContainer) +{ + using AllElements = std::remove_cvref_t; + using DataSeqElementPtr = std::tuple_element_t; + using DataSeqElement = typename DataSeqElementPtr::element_type; + using SeqElement = typename DataSeqElement::value_type; + using TopologyElement = typename SeqElement::value_type; + static constexpr auto elementType = TopologyElement::Element_type; + + const char* elementName = sofa::geometry::elementTypeToString(elementType); + auto elementNameStr = sofa::helper::downcaseString(std::string(elementName)); + + auto dataName = elementTypeToDataName(elementType); + std::string dataNameStr { dataName }; + sofa::helper::replaceAll(dataNameStr, " ", "_"); + + auto& dataPtr = std::get(dataContainer); + dataPtr = std::make_unique("List of " + elementNameStr); + + sofa::core::objectmodel::BaseData& data = *dataPtr; + data.setName(dataName); + + meshTopology->addData(dataPtr.get(), dataName); +} + +template +void registerData(core::objectmodel::BaseComponent* meshTopology, auto& dataContainer, std::index_sequence) +{ + (registerSingleData(meshTopology, dataContainer), ...); +} + +const void* MeshTopologyContainer::getElementsRaw( + const sofa::geometry::ElementType& elementType) const noexcept +{ + switch (elementType) + { + case ElementType::UNKNOWN: + case ElementType::POINT: + break; + case ElementType::EDGE: + return &get()->getValue(); + case ElementType::TRIANGLE: + return &get()->getValue(); + case ElementType::QUAD: + return &get()->getValue(); + case ElementType::TETRAHEDRON: + return &get()->getValue(); + case ElementType::HEXAHEDRON: + return &get()->getValue(); + case ElementType::PRISM: + return &get()->getValue(); + case ElementType::PYRAMID: + return &get()->getValue(); + case ElementType::QUADRATIC_EDGE: + return &get()->getValue(); + case ElementType::QUADRATIC_TRIANGLE: + return &get()->getValue(); + case ElementType::QUADRATIC_QUAD: + return &get()->getValue(); + case ElementType::QUADRATIC_TETRAHEDRON: + return &get()->getValue(); + case ElementType::QUADRATIC_HEXAHEDRON: + return &get()->getValue(); + case ElementType::QUADRATIC_PRISM: + return &get()->getValue(); + case ElementType::QUADRATIC_PYRAMID: + return &get()->getValue(); + case ElementType::SIZE: + default: + return nullptr; + } + return nullptr; +} + +MeshTopologyContainer::MeshTopologyContainer() : Inherit1() +{ + registerData(this, m_container, std::make_index_sequence>{}); +} + +const core::topology::BaseMeshTopology::SeqEdges& MeshTopologyContainer::getEdges() +{ + return get()->getValue(); +} +const core::topology::BaseMeshTopology::SeqTriangles& MeshTopologyContainer::getTriangles() +{ + return get()->getValue(); +} +const core::topology::BaseMeshTopology::SeqQuads& MeshTopologyContainer::getQuads() +{ + return get()->getValue(); +} +const core::topology::BaseMeshTopology::SeqTetrahedra& MeshTopologyContainer::getTetrahedra() +{ + return get()->getValue(); +} +const core::topology::BaseMeshTopology::SeqHexahedra& MeshTopologyContainer::getHexahedra() +{ + return get()->getValue(); +} +const core::topology::BaseMeshTopology::SeqPrisms& MeshTopologyContainer::getPrisms() +{ + return get()->getValue(); +} +const core::topology::BaseMeshTopology::SeqPyramids& MeshTopologyContainer::getPyramids() +{ + return get()->getValue(); +} + +} // namespace sofa::component::topology::container::constant diff --git a/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopologyContainer.h b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopologyContainer.h new file mode 100644 index 00000000000..a5e7c8f1295 --- /dev/null +++ b/Sofa/Component/Topology/Container/Constant/src/sofa/component/topology/container/constant/MeshTopologyContainer.h @@ -0,0 +1,184 @@ +/****************************************************************************** +* 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 + +namespace sofa::component::topology::container::constant +{ + +class SOFA_COMPONENT_TOPOLOGY_CONTAINER_CONSTANT_API MeshTopologyContainer : + public sofa::core::topology::BaseMeshTopology +{ +public: + SOFA_CLASS(MeshTopologyContainer, core::topology::BaseMeshTopology); + + template //e.g. sofa::geometry::Edge + using TopologyElement = sofa::topology::Element; + + template //e.g. sofa::geometry::Edge + using SeqElement = sofa::type::vector>; + + template //e.g. sofa::geometry::Edge + using DataSeqElement = sofa::Data>; + + template //e.g. sofa::geometry::Edge + using DataSeqElementPtr = std::unique_ptr>; + + using AllElements = std::tuple< + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr, + DataSeqElementPtr + >; + + AllElements m_container; + + const SeqEdges& getEdges() override; + const SeqTriangles& getTriangles() override; + const SeqQuads& getQuads() override; + const SeqTetrahedra& getTetrahedra() override; + const SeqHexahedra& getHexahedra() override; + const SeqPrisms& getPrisms() override; + const SeqPyramids& getPyramids() override; + + template + const DataSeqElementPtr& get() const + { + return std::get>(m_container); + } + + template + DataSeqElementPtr& get() + { + return std::get>(m_container); + } + +protected: + + const void* getElementsRaw(const sofa::geometry::ElementType& elementType) const noexcept override; + + MeshTopologyContainer(); + +public: + // A proxy for legacy data members + template + struct SeqElementProxy + { + MeshTopologyContainer* m_meshTopologyContainer { nullptr }; + SeqElementProxy(MeshTopologyContainer* meshTopologyContainer) : m_meshTopologyContainer(meshTopologyContainer) {} + + DataSeqElement& toData() + { + return *m_meshTopologyContainer->get(); + } + + const DataSeqElement& toData() const + { + return *m_meshTopologyContainer->get(); + } + + operator const DataSeqElement&() const + { + return toData(); + } + + operator DataSeqElement&() + { + return toData(); + } + + SeqElement* beginEdit() + { + return toData().beginEdit(); + } + + void endEdit() + { + toData().endEdit(); + } + + const SeqElement& getValue() const + { + return toData().getValue(); + } + + void setValue(const DataSeqElement& value) + { + toData().setValue(value); + } + + core::BaseData* getParent() const + { + return toData().getParent(); + } + + bool setParent(core::BaseData* parent, const std::string& path = std::string()) + { + return toData().setParent(parent, path); + } + + void delInput(core::objectmodel::DDGNode* n) + { + toData().delInput(n); + } + + SeqElement* beginWriteOnly() + { + return toData().beginWriteOnly(); + } + }; +}; + +} + +namespace sofa::helper +{ +template +auto getWriteOnlyAccessor(sofa::component::topology::container::constant::MeshTopologyContainer::SeqElementProxy& m_container) +{ + return getWriteOnlyAccessor(m_container.toData()); +} + +template +auto getWriteAccessor(sofa::component::topology::container::constant::MeshTopologyContainer::SeqElementProxy& m_container) +{ + return getWriteAccessor(m_container.toData()); +} + +template +auto getReadAccessor(const sofa::component::topology::container::constant::MeshTopologyContainer::SeqElementProxy& m_container) +{ + return getReadAccessor(m_container.toData()); +} + +} diff --git a/Sofa/Component/Topology/Container/Constant/tests/MeshTopology_test.cpp b/Sofa/Component/Topology/Container/Constant/tests/MeshTopology_test.cpp index e580a134430..b538206d6e0 100644 --- a/Sofa/Component/Topology/Container/Constant/tests/MeshTopology_test.cpp +++ b/Sofa/Component/Topology/Container/Constant/tests/MeshTopology_test.cpp @@ -129,6 +129,7 @@ bool MeshTopology_test::testHexahedronTopology() // Check tetrahedra container buffers size EXPECT_EQ(topo->getNbHexahedra(), nbrHexahedron); EXPECT_EQ(topo->getHexahedra().size(), nbrHexahedron); + EXPECT_EQ(topo->getElements().size(), nbrHexahedron); //// check tetrahedron buffer const sofa::type::vector& hexahedra1 = topoCon->getHexahedronArray(); @@ -243,6 +244,7 @@ bool MeshTopology_test::testTetrahedronTopology() // Check tetrahedra container buffers size EXPECT_EQ(topo->getNbTetrahedra(), nbrTetrahedron); EXPECT_EQ(topo->getTetrahedra().size(), nbrTetrahedron); + EXPECT_EQ(topo->getElements().size(), nbrTetrahedron); //// check tetrahedron buffer const sofa::type::vector& tetrahedra1 = topoCon->getTetrahedronArray(); @@ -357,6 +359,7 @@ bool MeshTopology_test::testQuadTopology() // Check quads container buffers size EXPECT_EQ(topo->getNbQuads(), nbrQuad); EXPECT_EQ(topo->getQuads().size(), nbrQuad); + EXPECT_EQ(topo->getElements().size(), nbrQuad); //// check quad buffer const sofa::type::vector& quads1 = topoCon->getQuadArray(); @@ -453,6 +456,7 @@ bool MeshTopology_test::testTriangleTopology() // Check triangles container buffers size EXPECT_EQ(topo->getNbTriangles(), nbrTriangle); EXPECT_EQ(topo->getTriangles().size(), nbrTriangle); + EXPECT_EQ(topo->getElements().size(), nbrTriangle); //// check triangle buffer const sofa::type::vector& triangles1 = topoCon->getTriangleArray(); @@ -549,6 +553,7 @@ bool MeshTopology_test::testEdgeTopology() // Check edge container buffers size EXPECT_EQ(topo->getNbEdges(), nbrEdge); EXPECT_EQ(topo->getEdges().size(), nbrEdge); + EXPECT_EQ(topo->getElements().size(), nbrEdge); //// check edge buffer const sofa::type::vector& edges1 = topoCon->getEdgeArray(); diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.cpp b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.cpp index 420b894798d..39243139685 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.cpp +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.cpp @@ -599,4 +599,14 @@ bool EdgeSetTopologyContainer::unlinkTopologyHandlerToData(core::topology::Topol } } +const void* EdgeSetTopologyContainer::getElementsRaw( + const sofa::geometry::ElementType& elementType) const noexcept +{ + if (elementType == sofa::geometry::ElementType::EDGE) + { + return &d_edge.getValue(); + } + return nullptr; +} + } //namespace sofa::component::topology::container::dynamic diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.h b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.h index 782642c3ffb..fb40e3a71dd 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.h +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/EdgeSetTopologyContainer.h @@ -180,6 +180,8 @@ class SOFA_COMPONENT_TOPOLOGY_CONTAINER_DYNAMIC_API EdgeSetTopologyContainer : p protected: + const void* getElementsRaw(const sofa::geometry::ElementType& elementType) const noexcept override; + /** \brief Creates the EdgeSet array. * * This function must be implemented by derived classes to create a list of edges from a set of triangles or tetrahedra diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.cpp b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.cpp index 54b02aa80fe..4acd009af55 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.cpp +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.cpp @@ -1246,4 +1246,14 @@ bool HexahedronSetTopologyContainer::unlinkTopologyHandlerToData(core::topology: } } +const void* HexahedronSetTopologyContainer::getElementsRaw( + const sofa::geometry::ElementType& elementType) const noexcept +{ + if (elementType == sofa::geometry::ElementType::HEXAHEDRON) + { + return &d_hexahedron.getValue(); + } + return QuadSetTopologyContainer::getElementsRaw(elementType); +} + } //namespace sofa::component::topology::container::dynamic diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.h b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.h index 9bb8f012910..0ac09e7bd65 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.h +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/HexahedronSetTopologyContainer.h @@ -318,6 +318,8 @@ class SOFA_COMPONENT_TOPOLOGY_CONTAINER_DYNAMIC_API HexahedronSetTopologyContain protected: + const void* getElementsRaw(const sofa::geometry::ElementType& elementType) const noexcept override; + /** \brief Creates the EdgeSet array. * * Create the set of edges when needed. diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.cpp b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.cpp index 3eb8784df6c..3f5a240ee40 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.cpp +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.cpp @@ -795,5 +795,14 @@ bool QuadSetTopologyContainer::unlinkTopologyHandlerToData(core::topology::Topol return EdgeSetTopologyContainer::unlinkTopologyHandlerToData(topologyHandler, elementType); } } +const void* QuadSetTopologyContainer::getElementsRaw( + const sofa::geometry::ElementType& elementType) const noexcept +{ + if (elementType == sofa::geometry::ElementType::QUAD) + { + return &d_quad.getValue(); + } + return EdgeSetTopologyContainer::getElementsRaw(elementType); +} } //namespace sofa::component::topology::container::dynamic diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.h b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.h index 58efedcc105..b29f65225f9 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.h +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/QuadSetTopologyContainer.h @@ -226,6 +226,8 @@ class SOFA_COMPONENT_TOPOLOGY_CONTAINER_DYNAMIC_API QuadSetTopologyContainer : p protected: + const void* getElementsRaw(const sofa::geometry::ElementType& elementType) const noexcept override; + /** \brief Creates the QuadSet array. * * This function must be implemented by derived classes to create a list of quads from a set of hexahedra for instance. diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.cpp b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.cpp index a55c321fccb..43521067823 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.cpp +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.cpp @@ -102,6 +102,16 @@ void TetrahedronSetTopologyContainer::createTetrahedronSetArray() msg_error() << "createTetrahedronSetArray method must be implemented by a child topology."; } +const void* TetrahedronSetTopologyContainer::getElementsRaw( + const sofa::geometry::ElementType& elementType) const noexcept +{ + if (elementType == sofa::geometry::ElementType::TETRAHEDRON) + { + return &d_tetrahedron.getValue(); + } + return TriangleSetTopologyContainer::getElementsRaw(elementType); +} + void TetrahedronSetTopologyContainer::createEdgeSetArray() { if(!hasTetrahedra()) // this method should only be called when tetrahedra exist diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.h b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.h index 2376324fbb2..29df1a2c3fb 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.h +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TetrahedronSetTopologyContainer.h @@ -294,6 +294,9 @@ class SOFA_COMPONENT_TOPOLOGY_CONTAINER_DYNAMIC_API TetrahedronSetTopologyContai friend std::istream& operator>>(std::istream& in, TetrahedronSetTopologyContainer& t); protected: + + const void* getElementsRaw(const sofa::geometry::ElementType& elementType) const noexcept override; + /** \brief Creates the EdgeSet array. * * Create the set of edges when needed. diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.cpp b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.cpp index 9ff133653cc..32d92295584 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.cpp +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.cpp @@ -1078,5 +1078,14 @@ bool TriangleSetTopologyContainer::unlinkTopologyHandlerToData(core::topology::T return EdgeSetTopologyContainer::unlinkTopologyHandlerToData(topologyHandler, elementType); } } +const void* TriangleSetTopologyContainer::getElementsRaw( + const sofa::geometry::ElementType& elementType) const noexcept +{ + if (elementType == sofa::geometry::ElementType::TRIANGLE) + { + return &d_triangle.getValue(); + } + return EdgeSetTopologyContainer::getElementsRaw(elementType); +} } //namespace sofa::component::topology::container::dynamic diff --git a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.h b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.h index ba0bfa3c8f9..f5d2639e1c3 100644 --- a/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.h +++ b/Sofa/Component/Topology/Container/Dynamic/src/sofa/component/topology/container/dynamic/TriangleSetTopologyContainer.h @@ -263,6 +263,8 @@ class SOFA_COMPONENT_TOPOLOGY_CONTAINER_DYNAMIC_API TriangleSetTopologyContainer protected: + const void* getElementsRaw(const sofa::geometry::ElementType& elementType) const noexcept override; + /** \brief Creates the TriangleSet array. * * This function must be implemented by derived classes to create a list of triangles from a set of tetrahedra for instance diff --git a/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/GridTopology.cpp b/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/GridTopology.cpp index 86c7551c976..2a23965f9b7 100644 --- a/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/GridTopology.cpp +++ b/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/GridTopology.cpp @@ -37,9 +37,9 @@ GridTopology::GridUpdate::GridUpdate(GridTopology *t): m_topology(t) { addInput(&t->d_n); - addOutput(&t->d_seqEdges); - addOutput(&t->d_seqQuads); - addOutput(&t->d_seqHexahedra); + addOutput(&t->d_seqEdges.toData()); + addOutput(&t->d_seqQuads.toData()); + addOutput(&t->d_seqHexahedra.toData()); setDirtyValue(); } diff --git a/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridRamificationTopology.cpp b/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridRamificationTopology.cpp index cdd3e7bc6ac..73e59b30d3c 100644 --- a/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridRamificationTopology.cpp +++ b/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridRamificationTopology.cpp @@ -848,8 +848,8 @@ void SparseGridRamificationTopology::buildVirtualFinerLevels() { _virtualFinerLevels[0]->d_seqPoints.setParent(&this->d_seqPoints); _virtualFinerLevels[0]->d_facets.setParent(&this->d_facets); - _virtualFinerLevels[0]->d_seqTriangles.setParent(&this->d_seqTriangles); - _virtualFinerLevels[0]->d_seqQuads.setParent(&this->d_seqQuads); + _virtualFinerLevels[0]->d_seqTriangles.setParent(&this->d_seqTriangles.toData()); + _virtualFinerLevels[0]->d_seqQuads.setParent(&this->d_seqQuads.toData()); } else _virtualFinerLevels[0]->load(fileTopology.c_str()); diff --git a/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridTopology.cpp b/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridTopology.cpp index 3ad22bb854a..de497b8e115 100644 --- a/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridTopology.cpp +++ b/Sofa/Component/Topology/Container/Grid/src/sofa/component/topology/container/grid/SparseGridTopology.cpp @@ -1243,8 +1243,8 @@ void SparseGridTopology::buildVirtualFinerLevels() { _virtualFinerLevels[0]->d_seqPoints.setParent(&this->d_seqPoints); _virtualFinerLevels[0]->d_facets.setParent(&this->d_facets); - _virtualFinerLevels[0]->d_seqTriangles.setParent(&this->d_seqTriangles); - _virtualFinerLevels[0]->d_seqQuads.setParent(&this->d_seqQuads); + _virtualFinerLevels[0]->d_seqTriangles.setParent(&this->d_seqTriangles.toData()); + _virtualFinerLevels[0]->d_seqQuads.setParent(&this->d_seqQuads.toData()); } else _virtualFinerLevels[0]->load(fileTopology.c_str()); diff --git a/Sofa/Component/Topology/Container/Grid/tests/SparseGridTopology_test.cpp b/Sofa/Component/Topology/Container/Grid/tests/SparseGridTopology_test.cpp index 406356b5568..5d1a0d64414 100644 --- a/Sofa/Component/Topology/Container/Grid/tests/SparseGridTopology_test.cpp +++ b/Sofa/Component/Topology/Container/Grid/tests/SparseGridTopology_test.cpp @@ -93,7 +93,7 @@ bool SparseGridTopology_test::buildFromMeshParams() //Pyramid centered on 0 0 0 sparseGrid1->d_seqPoints.setValue({{0, 0, 1}, {-1, 0, -1}, {0, 1, -1}, {1, 0, -1}, {0, -1, -1} }); - sparseGrid1->d_seqTriangles.setValue({{0, 1, 2}, {0, 2, 3}, {0, 3, 4}, {0, 4, 1}, {1, 2, 3}, {3, 4, 1} }); + sparseGrid1->d_seqTriangles.toData().setValue({{0, 1, 2}, {0, 2, 3}, {0, 3, 4}, {0, 4, 1}, {1, 2, 3}, {3, 4, 1} }); sparseGrid1->setN({ 10,10,10 }); sparseGrid1->init(); diff --git a/Sofa/framework/Core/src/sofa/core/topology/BaseMeshTopology.h b/Sofa/framework/Core/src/sofa/core/topology/BaseMeshTopology.h index 7910c16640f..9ba88a6e1ef 100644 --- a/Sofa/framework/Core/src/sofa/core/topology/BaseMeshTopology.h +++ b/Sofa/framework/Core/src/sofa/core/topology/BaseMeshTopology.h @@ -108,6 +108,32 @@ class SOFA_CORE_API BaseMeshTopology : public core::topology::Topology virtual const SeqHexahedra& getHexahedra() = 0; virtual const SeqPrisms& getPrisms() = 0; virtual const SeqPyramids& getPyramids() = 0; + + /** + * Generic elements' accessor. + * + * Relies on the virtual method getElementsRaw. If this method is not implemented by the derived + * class, or if the element type is not supported by the derived class, the method returns an + * empty element list. + * + * @tparam ElementType the type of elements, e.g. sofa::geometry::Edge + * @return The element list contained in the topology + */ + template // e.g. sofa::geometry::Edge + const auto& getElements() const + { + using TopologyElement = sofa::topology::Element; + using SeqElement = sofa::type::vector; + + const void* rawSequence = getElementsRaw(ElementType::Element_type); + if (!rawSequence) + { + static SeqElement empty; + return empty; + } + return *static_cast(rawSequence); + } + /// @} /// Random accessors @@ -335,6 +361,13 @@ class SOFA_CORE_API BaseMeshTopology : public core::topology::Topology sofa::core::objectmodel::DataFileName fileTopology; + + virtual const void* getElementsRaw(const sofa::geometry::ElementType& elementType) const noexcept + { + SOFA_UNUSED(elementType); + return nullptr; + } + public: bool insertInNode( objectmodel::BaseNode* node ) override; diff --git a/Sofa/framework/Core/src/sofa/core/topology/Topology.cpp b/Sofa/framework/Core/src/sofa/core/topology/Topology.cpp index 05312726dc8..43a6db3c8dc 100644 --- a/Sofa/framework/Core/src/sofa/core/topology/Topology.cpp +++ b/Sofa/framework/Core/src/sofa/core/topology/Topology.cpp @@ -40,6 +40,18 @@ bool Topology::removeInNode(objectmodel::BaseNode* node) return true; } + +// Quadratic triangle subdivision +const unsigned int trianglesInQuadraticTriangles[4][3] = { {0,3,4}, {3,1,5}, {3,5,4}, {4,5,2}}; + +// Quadratic Quad subdivision into 4 linear quads +const unsigned int quadsInQuadraticQuads[4][4] = { + {0, 4, 8, 7}, // Sub-quad 0 + {4, 1, 5, 8}, // Sub-quad 1 + {8, 5, 2, 6}, // Sub-quad 2 + {7, 8, 6, 3} // Sub-quad 3 +}; + // Tetrahedron const unsigned int edgesInTetrahedronArray[6][2] = {{0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3}}; const unsigned int trianglesOrientationInTetrahedronArray[4][3] = {{1,2,3}, {0,3,2}, {1,3,0}, {0,2,1}}; @@ -49,4 +61,19 @@ const unsigned int edgesInHexahedronArray[12][2] = {{0,1},{0,3},{0,4},{1,2},{1,5 const unsigned int quadsOrientationInHexahedronArray[6][4] = {{0,3,2,1}, {4,5,6,7}, {0,1,5,4}, {1,2,6,5}, {2,3,7,6}, {3,0,4,7}}; const unsigned int verticesInHexahedronArray[2][2][2] = { {{0,4}, {3,7}}, {{1,5}, {2,6}} }; +// Quadratic tetrahedron +const unsigned int quadraticTrianglesInQuadraticTetrahedronArray[4][6] = { {0,1,3,4,6,8}, {2,0,3,5,9,6}, {1,2,3,7,8,9}, {2,1,0,7,5,4} }; + +// Quadratic hexahedron (6 quadratic quads of 9 nodes each) +const unsigned int quadraticQuadsInQuadraticHexahedronArray[6][9] = { + {0,4,7,3,10,17,15,9,25}, // Face 0: -X + {1,2,6,5,11,14,18,12,23}, // Face 1: +X + {0,1,5,4,8,12,16,10,22}, // Face 2: -Y + {2,3,7,6,13,15,19,14,24}, // Face 3: +Y + {1,0,3,2,8,9,13,11,20}, // Face 4: -Z + {6,7,4,5,19,17,16,18,21} // Face 5: +Z +}; + + + } // namespace sofa::core::topology diff --git a/Sofa/framework/Core/src/sofa/core/topology/Topology.h b/Sofa/framework/Core/src/sofa/core/topology/Topology.h index 6d27818cf18..3bd51b979ed 100644 --- a/Sofa/framework/Core/src/sofa/core/topology/Topology.h +++ b/Sofa/framework/Core/src/sofa/core/topology/Topology.h @@ -119,6 +119,12 @@ class SOFA_CORE_API Topology : public virtual sofa::core::objectmodel::BaseCompo virtual SReal getPZ(Index /*i*/) const { return 0.0; } }; +/// List of linear sub-triangles in quadratic triangle +SOFA_CORE_API extern const unsigned int trianglesInQuadraticTriangles[4][3]; + +/// List of linear sub-quads in a quadratic quad +SOFA_CORE_API extern const unsigned int quadsInQuadraticQuads[4][4]; + /// List of pair of vertex indices (edge) in a tetrahedron SOFA_CORE_API extern const unsigned int edgesInTetrahedronArray[6][2]; /// List of 3 vertex indices (triangle) in a tetrahedron @@ -132,6 +138,12 @@ SOFA_CORE_API extern const unsigned int quadsOrientationInHexahedronArray[6][4]; // List of vertex indices in a hexahedron SOFA_CORE_API extern const unsigned int verticesInHexahedronArray[2][2][2]; +/// List of quadratic triangles in a quadratic tetrahedron +SOFA_CORE_API extern const unsigned int quadraticTrianglesInQuadraticTetrahedronArray[4][6]; + +/// List of quadratic quads in a quadratic hexahedron +SOFA_CORE_API extern const unsigned int quadraticQuadsInQuadraticHexahedronArray[6][9]; + } // namespace sofa::core::topology namespace sofa::geometry diff --git a/Sofa/framework/Core/src/sofa/core/visual/DrawMesh.h b/Sofa/framework/Core/src/sofa/core/visual/DrawMesh.h index f68114a8e84..48f1d70674b 100644 --- a/Sofa/framework/Core/src/sofa/core/visual/DrawMesh.h +++ b/Sofa/framework/Core/src/sofa/core/visual/DrawMesh.h @@ -26,6 +26,7 @@ #include #include #include +#include #include @@ -66,7 +67,7 @@ struct BaseDrawMesh sofa::core::topology::BaseMeshTopology* topology, const ColorContainer& colors = Derived::defaultColors) { - const auto totalNbElements = topology->getNbElements(); + const auto totalNbElements = topology->getElements().size(); const auto elementsToDraw = sofa::helper::IotaView(static_cast(0), totalNbElements); drawSomeElements(drawTool, position, topology, elementsToDraw, colors); } @@ -595,6 +596,562 @@ struct SOFA_CORE_API DrawElementMesh } }; + +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 2> +{ + using ElementType = sofa::geometry::QuadraticEdge; + friend BaseDrawMesh; + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor::lime(), + sofa::type::RGBAColor::silver() + }; + +private: + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getElements(); + + // Draw as 2 line segments per quadratic edge + constexpr std::size_t nbSegments = 2; + const auto size = elementIndices.size() * nbSegments * 2; + for (auto& p : renderedPoints) + { + p.resize(size); + } + + std::array renderedPointId {}; + for (auto i : elementIndices) + { + const auto& element = elements[i]; + const auto center = this->elementCenter(position, sofa::topology::Edge{element[0], element[1]}); + + // Segment 0-4 (vertex 0 to mid-edge) + renderedPoints[i % NumberColors][renderedPointId[i % NumberColors]++] = + sofa::type::toVec3(this->applyElementSpace(position[element[0]], center)); + renderedPoints[i % NumberColors][renderedPointId[i % NumberColors]++] = + sofa::type::toVec3(this->applyElementSpace(position[element[2]], center)); + + // Segment 4-1 (mid-edge to vertex 1) + renderedPoints[i % NumberColors][renderedPointId[i % NumberColors]++] = + sofa::type::toVec3(this->applyElementSpace(position[element[2]], center)); + renderedPoints[i % NumberColors][renderedPointId[i % NumberColors]++] = + sofa::type::toVec3(this->applyElementSpace(position[element[1]], center)); + } + + for (std::size_t j = 0; j < NumberColors; ++j) + { + drawTool->drawLines(renderedPoints[j], 2.f, colors[j]); + } + } +}; + + +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 3> +{ + using ElementType = sofa::geometry::QuadraticTriangle; + friend BaseDrawMesh; + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor::green(), + sofa::type::RGBAColor::teal(), + sofa::type::RGBAColor::blue() + }; + +private: + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getElements(); + + // Each quadratic triangle is subdivided into 4 triangles + constexpr std::size_t nbSubTriangles = 4; + for (auto& p : renderedPoints) + { + p.resize(elementIndices.size() * nbSubTriangles * sofa::geometry::Triangle::NumberOfNodes); + } + + // Sub-triangles: corner triangles + center triangle + static constexpr std::array, 4> subTriangles {{ + {{0, 3, 5}}, // corner at vertex 0 + {{3, 1, 4}}, // corner at vertex 1 + {{5, 4, 2}}, // corner at vertex 2 + {{3, 4, 5}} // center triangle + }}; + + for (std::size_t i = 0; i < elementIndices.size(); ++i) + { + const auto& element = elements[elementIndices[i]]; + const auto center = this->elementCenter(position, sofa::topology::Triangle{element[0], element[1], element[2]}); + + for (std::size_t t = 0; t < nbSubTriangles; ++t) + { + for (std::size_t v = 0; v < 3; ++v) + { + const auto vertexId = element[subTriangles[t][v]]; + const auto p = this->applyElementSpace(position[vertexId], center); + renderedPoints[(i + t) % NumberColors][(i * nbSubTriangles + t) * 3 + v] = sofa::type::toVec3(p); + } + } + } + + for (std::size_t j = 0; j < NumberColors; ++j) + { + drawTool->drawTriangles(renderedPoints[j], colors[j]); + } + } +}; + + +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 2> +{ + using ElementType = sofa::geometry::QuadraticQuad; + friend BaseDrawMesh; + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor::green(), + sofa::type::RGBAColor::orange() + }; + +private: + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getElements(); + + // Each quadratic quad is subdivided into 4 quads + constexpr std::size_t nbSubQuads = 4; + const auto size = elementIndices.size() * nbSubQuads * sofa::geometry::Quad::NumberOfNodes; + for (auto& p : renderedPoints) + { + p.resize(size); + } + + // Sub-quads using corner vertices, mid-edges, and center + static constexpr std::array, 4> subQuads {{ + {{0, 4, 8, 7}}, // bottom-left + {{4, 1, 5, 8}}, // bottom-right + {{8, 5, 2, 6}}, // top-right + {{7, 8, 6, 3}} // top-left + }}; + + std::array renderedPointId {}; + for (auto i : elementIndices) + { + const auto& element = elements[i]; + const auto center = this->elementCenter(position, sofa::topology::Quad{element[0], element[1], element[2], element[3]}); + + for (std::size_t q = 0; q < nbSubQuads; ++q) + { + for (std::size_t v = 0; v < 4; ++v) + { + const auto vertexId = element[subQuads[q][v]]; + const auto p = this->applyElementSpace(position[vertexId], center); + renderedPoints[i % NumberColors][renderedPointId[i % NumberColors]++] = sofa::type::toVec3(p); + } + } + } + + for (std::size_t j = 0; j < NumberColors; ++j) + { + drawTool->drawQuads(renderedPoints[j], colors[j]); + } + } +}; + +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 4> +{ + using ElementType = sofa::geometry::QuadraticTetrahedron; + friend BaseDrawMesh; + + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor::blue(), + sofa::type::RGBAColor::black(), + sofa::type::RGBAColor::azure(), + sofa::type::RGBAColor::cyan() + }; + +private: + template + void doDraw(sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getElements(); + + // Each QuadraticTetrahedron has 4 quadratic faces. + // Each quadratic face is drawn as 4 triangles. + constexpr std::size_t nbTrianglesPerFace = 4; + constexpr std::size_t nbFaces = 4; + + for (auto& p : renderedPoints) + { + p.resize(elementIndices.size() * nbTrianglesPerFace * sofa::geometry::Triangle::NumberOfNodes); + } + + for (std::size_t i = 0; i < elementIndices.size(); ++i) + { + const auto& element = elements[elementIndices[i]]; + const auto center = this->elementCenter(position, sofa::topology::Tetrahedron{element[0], element[1], element[2], element[3]}); + + for (std::size_t f = 0; f < nbFaces; ++f) + { + const auto& face = topology::quadraticTrianglesInQuadraticTetrahedronArray[f]; + const sofa::topology::QuadraticTriangle FACE( + element[face[0]], element[face[1]], element[face[2]], + element[face[3]], element[face[4]], element[face[5]]); + + for (std::size_t t = 0; t < nbTrianglesPerFace; ++t) + { + const auto& triangle = topology::trianglesInQuadraticTriangles[t]; + const sofa::topology::Triangle TRIANGLE(FACE[triangle[0]], FACE[triangle[1]], FACE[triangle[2]]); + + for (std::size_t v = 0; v < 3; ++v) + { + const auto vertexId = TRIANGLE[v]; + const auto p = this->applyElementSpace(position[vertexId], center); + renderedPoints[f][(i * nbTrianglesPerFace + t) * 3 + v] = sofa::type::toVec3(p); + } + } + } + } + + for (std::size_t f = 0; f < nbFaces; ++f) + { + // Generate slightly different colors for the 4 sub-triangles of a face + std::array subColors; + for (int t = 0; t < 4; ++t) + { + subColors[t] = type::RGBAColor::lighten(colors[f], t * 0.15_sreal); + } + + // Draw each sub-triangle set with its specific color + for (std::size_t t = 0; t < nbTrianglesPerFace; ++t) + { + // We need a temporary view of the buffer for this specific sub-triangle across all tetrahedra + sofa::type::vector subBuffer; + subBuffer.resize(elementIndices.size() * 3); + for(std::size_t i = 0; i < elementIndices.size(); ++i) + { + subBuffer[i*3 + 0] = renderedPoints[f][(i * nbTrianglesPerFace + t) * 3 + 0]; + subBuffer[i*3 + 1] = renderedPoints[f][(i * nbTrianglesPerFace + t) * 3 + 1]; + subBuffer[i*3 + 2] = renderedPoints[f][(i * nbTrianglesPerFace + t) * 3 + 2]; + } + drawTool->drawTriangles(subBuffer, subColors[t]); + } + } + } +}; + + +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 6> +{ + using ElementType = sofa::geometry::QuadraticHexahedron; + friend BaseDrawMesh; + static constexpr std::size_t NumberQuadsInHexahedron = 6; + + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor(0.7f,0.7f,0.1f,1.f), + sofa::type::RGBAColor(0.7f,0.0f,0.0f,1.f), + sofa::type::RGBAColor(0.0f,0.7f,0.0f,1.f), + sofa::type::RGBAColor(0.0f,0.0f,0.7f,1.f), + sofa::type::RGBAColor(0.1f,0.7f,0.7f,1.f), + sofa::type::RGBAColor(0.7f,0.1f,0.7f,1.f) + }; + +private: + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getElements(); + + // Each face (6 total) is subdivided into 4 quads + constexpr std::size_t nbQuadsPerFace = 4; + for (auto& p : renderedPoints) + { + p.resize(elementIndices.size() * nbQuadsPerFace * sofa::geometry::Quad::NumberOfNodes); + } + + for (std::size_t i = 0; i < elementIndices.size(); ++i) + { + const auto& element = elements[elementIndices[i]]; + const auto center = this->elementCenter(position, + sofa::topology::Hexahedron{element[0], element[1], element[2], element[3], element[4], element[5], element[6], element[7]}); + + for (std::size_t f = 0; f < NumberQuadsInHexahedron; ++f) + { + const auto& faceIndices = topology::quadraticQuadsInQuadraticHexahedronArray[f]; + + for (std::size_t q = 0; q < nbQuadsPerFace; ++q) + { + const auto& subQuad = topology::quadsInQuadraticQuads[q]; + for (std::size_t v = 0; v < 4; ++v) + { + const auto vertexId = element[faceIndices[subQuad[v]]]; + const auto p = this->applyElementSpace(position[vertexId], center); + renderedPoints[f][(i * nbQuadsPerFace + q) * 4 + v] = sofa::type::toVec3(p); + } + } + } + } + + for (std::size_t f = 0; f < NumberQuadsInHexahedron; ++f) + { + drawTool->drawQuads(renderedPoints[f], colors[f]); + } + } +}; + + +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 5> +{ + using ElementType = sofa::geometry::QuadraticPrism; + friend BaseDrawMesh; + + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor::green(), + sofa::type::RGBAColor::teal(), + sofa::type::RGBAColor::navy(), + sofa::type::RGBAColor::gold(), + sofa::type::RGBAColor::purple() + }; + +private: + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getElements(); + + // 2 triangular faces (each subdivided into 4 triangles) + 3 quad faces (each subdivided into 4 quads) + constexpr std::size_t nbSubTriangles = 4; + constexpr std::size_t nbSubQuads = 4; + + renderedPoints[0].resize(elementIndices.size() * nbSubTriangles * 3); + renderedPoints[1].resize(elementIndices.size() * nbSubTriangles * 3); + renderedPoints[2].resize(elementIndices.size() * nbSubQuads * 4); + renderedPoints[3].resize(elementIndices.size() * nbSubQuads * 4); + renderedPoints[4].resize(elementIndices.size() * nbSubQuads * 4); + + // Bottom triangle (0,1,2) subdivided + static constexpr std::array, 4> bottomTriangles {{ + {{0, 6, 8}}, {{6, 1, 7}}, {{8, 7, 2}}, {{6, 7, 8}} + }}; + // Top triangle (3,4,5) subdivided + static constexpr std::array, 4> topTriangles {{ + {{3, 12, 14}}, {{12, 4, 13}}, {{14, 13, 5}}, {{12, 13, 14}} + }}; + // Quad faces subdivided (using corner vertices, mid-edges, and face centers) + static constexpr std::array, 4> quadSubdivision {{ + {{0, 6, 15, 9}}, {{6, 1, 10, 15}}, {{15, 10, 4, 12}}, {{9, 15, 12, 3}} + }}; + + std::array pointId{}; + + for (auto idx : elementIndices) + { + const auto& element = elements[idx]; + const auto center = this->elementCenter(position, + sofa::topology::Prism{element[0], element[1], element[2], element[3], element[4], element[5]}); + + // Bottom triangle + for (const auto& tri : bottomTriangles) + { + for (std::size_t v = 0; v < 3; ++v) + { + const auto p = this->applyElementSpace(position[element[tri[v]]], center); + renderedPoints[0][pointId[0]++] = sofa::type::toVec3(p); + } + } + + // Top triangle + for (const auto& tri : topTriangles) + { + for (std::size_t v = 0; v < 3; ++v) + { + const auto p = this->applyElementSpace(position[element[tri[v]]], center); + renderedPoints[1][pointId[1]++] = sofa::type::toVec3(p); + } + } + + // Three quad faces (0-1-4-3), (1-2-5-4), (2-0-3-5) + std::array, 3> quadFaceNodes {{ + {{0, 1, 4, 3, 6, 10, 12, 9, 15}}, // face 0-1-4-3 + {{1, 2, 5, 4, 7, 11, 13, 10, 16}}, // face 1-2-5-4 + {{2, 0, 3, 5, 8, 9, 14, 11, 17}} // face 2-0-3-5 + }}; + + for (std::size_t faceIdx = 0; faceIdx < 3; ++faceIdx) + { + const auto& faceNodes = quadFaceNodes[faceIdx]; + for (const auto& quad : quadSubdivision) + { + for (std::size_t v = 0; v < 4; ++v) + { + const auto localIdx = quad[v]; + const auto vertexId = element[faceNodes[localIdx]]; + const auto p = this->applyElementSpace(position[vertexId], center); + renderedPoints[2 + faceIdx][pointId[2 + faceIdx]++] = sofa::type::toVec3(p); + } + } + } + } + + drawTool->drawTriangles(renderedPoints[0], colors[0]); + drawTool->drawTriangles(renderedPoints[1], colors[1]); + drawTool->drawQuads(renderedPoints[2], colors[2]); + drawTool->drawQuads(renderedPoints[3], colors[3]); + drawTool->drawQuads(renderedPoints[4], colors[4]); + } +}; + + +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 5> +{ + using ElementType = sofa::geometry::QuadraticPyramid; + friend BaseDrawMesh; + + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor::green(), + sofa::type::RGBAColor::teal(), + sofa::type::RGBAColor::navy(), + sofa::type::RGBAColor::gold(), + sofa::type::RGBAColor::purple() + }; + +private: + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getElements(); + + // 1 quad base (subdivided into 4 quads) + 4 triangular faces (each subdivided into 4 triangles) + constexpr std::size_t nbSubQuads = 4; + constexpr std::size_t nbSubTriangles = 4; + + renderedPoints[0].resize(elementIndices.size() * nbSubQuads * 4); + renderedPoints[1].resize(elementIndices.size() * nbSubTriangles * 3); + renderedPoints[2].resize(elementIndices.size() * nbSubTriangles * 3); + renderedPoints[3].resize(elementIndices.size() * nbSubTriangles * 3); + renderedPoints[4].resize(elementIndices.size() * nbSubTriangles * 3); + + // Base quad (0,1,2,3) subdivided with center node 13 + static constexpr std::array, 4> baseQuadSubdivision {{ + {{0, 5, 13, 8}}, {{5, 1, 6, 13}}, {{13, 6, 2, 7}}, {{8, 13, 7, 3}} + }}; + + // Triangular faces subdivided + static constexpr std::array, 4>, 4> triangleFaces {{ + {{{0, 5, 9}, {5, 1, 10}, {9, 10, 4}, {5, 10, 9}}}, // face 0-1-4 + {{{1, 6, 10}, {6, 2, 11}, {10, 11, 4}, {6, 11, 10}}}, // face 1-2-4 + {{{2, 7, 11}, {7, 3, 12}, {11, 12, 4}, {7, 12, 11}}}, // face 2-3-4 + {{{3, 8, 12}, {8, 0, 9}, {12, 9, 4}, {8, 9, 12}}} // face 3-0-4 + }}; + + std::array pointId{}; + + for (auto idx : elementIndices) + { + const auto& element = elements[idx]; + const auto center = this->elementCenter(position, + sofa::topology::Pyramid{element[0], element[1], element[2], element[3], element[4]}); + + // Draw base quad + for (const auto& quad : baseQuadSubdivision) + { + for (std::size_t v = 0; v < 4; ++v) + { + const auto p = this->applyElementSpace(position[element[quad[v]]], center); + renderedPoints[0][pointId[0]++] = sofa::type::toVec3(p); + } + } + + // Draw 4 triangular faces + for (std::size_t faceIdx = 0; faceIdx < 4; ++faceIdx) + { + for (const auto& tri : triangleFaces[faceIdx]) + { + for (std::size_t v = 0; v < 3; ++v) + { + const auto p = this->applyElementSpace(position[element[tri[v]]], center); + renderedPoints[1 + faceIdx][pointId[1 + faceIdx]++] = sofa::type::toVec3(p); + } + } + } + } + + drawTool->drawQuads(renderedPoints[0], colors[0]); + drawTool->drawTriangles(renderedPoints[1], colors[1]); + drawTool->drawTriangles(renderedPoints[2], colors[2]); + drawTool->drawTriangles(renderedPoints[3], colors[3]); + drawTool->drawTriangles(renderedPoints[4], colors[4]); + } +}; + class SOFA_CORE_API DrawMesh { public: @@ -615,6 +1172,7 @@ class SOFA_CORE_API DrawMesh void drawLine(sofa::helper::visual::DrawTool* drawTool, const PositionContainer& position, sofa::core::topology::BaseMeshTopology* topology) { drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); } template @@ -622,15 +1180,21 @@ class SOFA_CORE_API DrawMesh { drawElements(drawTool, position, topology); drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); } template void drawVolume(sofa::helper::visual::DrawTool* drawTool, const PositionContainer& position, sofa::core::topology::BaseMeshTopology* topology) { drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); } template @@ -642,16 +1206,22 @@ class SOFA_CORE_API DrawMesh } const auto hasTriangles = !topology->getTriangles().empty(); + const auto hasQTriangle = !topology->getElements().empty(); const auto hasQuads = !topology->getQuads().empty(); + const auto hasQQuad = !topology->getElements().empty(); - const auto hasSurfaceElements = hasTriangles || hasQuads; + const auto hasSurfaceElements = hasTriangles || hasQTriangle || hasQuads || hasQQuad; const auto hasTetra = !topology->getTetrahedra().empty(); + const auto hasQTetra = !topology->getElements().empty(); const auto hasHexa = !topology->getHexahedra().empty(); + const auto hasQHexa = !topology->getElements().empty(); const auto hasPrism = !topology->getPrisms().empty(); + const auto hasQPrism = !topology->getElements().empty(); const auto hasPyramid = !topology->getPyramids().empty(); + const auto hasQPyramid = !topology->getElements().empty(); - const bool hasVolumeElements = hasTetra || hasHexa || hasPrism || hasPyramid; + const bool hasVolumeElements = hasTetra || hasQTetra || hasHexa || hasQHexa || hasPrism || hasQPrism || hasPyramid || hasQPyramid; if (!hasSurfaceElements && !hasVolumeElements) { @@ -673,12 +1243,19 @@ class SOFA_CORE_API DrawMesh private: std::tuple< DrawElementMesh, + DrawElementMesh, DrawElementMesh, + DrawElementMesh, DrawElementMesh, + DrawElementMesh, DrawElementMesh, + DrawElementMesh, DrawElementMesh, + DrawElementMesh, DrawElementMesh, - DrawElementMesh + DrawElementMesh, + DrawElementMesh, + DrawElementMesh > m_meshes; }; diff --git a/Sofa/framework/FEM/CMakeLists.txt b/Sofa/framework/FEM/CMakeLists.txt index 5e88eb85484..68318f15885 100644 --- a/Sofa/framework/FEM/CMakeLists.txt +++ b/Sofa/framework/FEM/CMakeLists.txt @@ -9,6 +9,7 @@ set(HEADER_FILES ${SOFAFEMSRC_ROOT}/FiniteElement.h ${SOFAFEMSRC_ROOT}/FiniteElement[all].h + ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].h ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].h ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].h @@ -16,6 +17,14 @@ set(HEADER_FILES ${SOFAFEMSRC_ROOT}/FiniteElement[Quad].h ${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].h ${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].h + + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticEdge].h + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticHexahedron].h + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticPrism].h + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticPyramid].h + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticQuad].h + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticTetrahedron].h + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticTriangle].h ) set(SOURCE_FILES @@ -26,8 +35,17 @@ set(SOURCE_FILES ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Pyramid].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Quad].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticTetrahedron].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].cpp + + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticEdge].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticHexahedron].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticPrism].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticPyramid].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticQuad].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticTetrahedron].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[QuadraticTriangle].cpp ) find_package(Sofa.Core REQUIRED) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].cpp new file mode 100644 index 00000000000..87e6cdf7cd5 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].cpp @@ -0,0 +1,33 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_QUADRATIC_EDGE_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; +template struct SOFA_FEM_API FiniteElement; +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h new file mode 100644 index 00000000000..15bae01c757 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h @@ -0,0 +1,93 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_EDGE_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1); + + constexpr static std::array referenceElementNodes {{ + ReferenceCoord{0}, // vertex 0 + ReferenceCoord{1}, // vertex 1 + ReferenceCoord{0.5} // mid-edge node + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getElements(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + const Real xi = q[0]; + return { + 2 * xi * (xi - 0.5), // vertex 0: (2*xi - 1) * xi + 2 * xi * (xi - 0.5) + 1 - 2 * xi, // vertex 1: (2*xi - 1) * (xi - 1) = 2*xi^2 - 3*xi + 1 + 4 * xi * (1 - xi) // mid-edge: 4*xi*(1-xi) + }; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const Real xi = q[0]; + return { + {4 * xi - 3}, // vertex 0 + {4 * xi - 1}, // vertex 1 + {4 - 8 * xi} // mid-edge + }; + } + + static constexpr std::array quadraturePoints() + { + // constexpr Real a = (1. - 1. / std::sqrt(3.)) / 2.; + constexpr Real a { 0.211324865405 }; + // constexpr Real b = (1. + 1. / std::sqrt(3.)) / 2.; + constexpr Real b { 0.788675134595 }; + constexpr Real w = 0.5; + + constexpr sofa::type::Vec q0(a); + constexpr sofa::type::Vec q1(b); + + constexpr std::array q { + std::make_pair(q0, w), + std::make_pair(q1, w) + }; + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_EDGE_CPP) +extern template struct SOFA_FEM_API FiniteElement; +extern template struct SOFA_FEM_API FiniteElement; +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].cpp new file mode 100644 index 00000000000..70e60f9cec5 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].cpp @@ -0,0 +1,31 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_QUADRATIC_HEXAHEDRON_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h new file mode 100644 index 00000000000..733b39707ba --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h @@ -0,0 +1,255 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_HEXAHEDRON_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::QuadraticHexahedron, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Quadratic Hexahedrons are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + + // 8 corner vertices + /*00*/ {-1, -1, -1}, // vertex 0 + /*01*/ { 1, -1, -1}, // vertex 1 + /*02*/ { 1, 1, -1}, // vertex 2 + /*03*/ {-1, 1, -1}, // vertex 3 + /*04*/ {-1, -1, 1}, // vertex 4 + /*05*/ { 1, -1, 1}, // vertex 5 + /*06*/ { 1, 1, 1}, // vertex 6 + /*07*/ {-1, 1, 1}, // vertex 7 + + // 12 mid-edge nodes + /*08*/ { 0, -1, -1}, // mid-edge 0-1 + /*09*/ {-1, 0, -1}, // mid-edge 0-3 + /*10*/ {-1, -1, 0}, // mid-edge 0-4 + /*11*/ { 1, 0, -1}, // mid-edge 1-2 + /*12*/ { 1, -1, 0}, // mid-edge 1-5 + /*13*/ { 0, 1, -1}, // mid-edge 2-3 + /*14*/ { 1, 1, 0}, // mid-edge 2-6 + /*15*/ {-1, 1, 0}, // mid-edge 3-7 + /*16*/ { 0, -1, 1}, // mid-edge 4-5 + /*17*/ {-1, 0, 1}, // mid-edge 7-4 + /*18*/ { 1, 0, 1}, // mid-edge 5-6 + /*19*/ { 0, 1, 1}, // mid-edge 6-7 + + // 6 face-center nodes + /*20*/ { 0, 0, -1}, // face center bottom (0-1-2-3) + /*21*/ { 0, 0, 1}, // face center top (4-5-6-7) + /*22*/ { 0, -1, 0}, // face center front (0-1-5-4) + /*23*/ { 1, 0, 0}, // face center right (1-2-6-5) + /*24*/ { 0, 1, 0}, // face center back (2-3-7-6) + /*25*/ {-1, 0, 0}, // face center left (3-0-4-7) + + // 1 volume center node + /*26*/ {0, 0, 0} // volume center + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getElements(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + const auto N0 = (q[0] - 1) * (q[1] - 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 8; + const auto N1 = (q[0] + 1) * (q[1] - 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 8; + const auto N2 = (q[0] + 1) * (q[1] + 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 8; + const auto N3 = (q[0] - 1) * (q[1] + 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 8; + const auto N4 = (q[0] - 1) * (q[1] - 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 8; + const auto N5 = (q[0] + 1) * (q[1] - 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 8; + const auto N6 = (q[0] + 1) * (q[1] + 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 8; + const auto N7 = (q[0] - 1) * (q[1] + 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 8; + const auto N8 = -(q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[2] - 1) * q[1] * q[2] / 4; + const auto N9 = -(q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * q[0] * q[2] / 4; + const auto N10 = -(q[0] - 1) * (q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1] / 4; + const auto N11 = -(q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * q[0] * q[2] / 4; + const auto N12 = -(q[0] + 1) * (q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1] / 4; + const auto N13 = -(q[0] - 1) * (q[0] + 1) * (q[1] + 1) * (q[2] - 1) * q[1] * q[2] / 4; + const auto N14 = -(q[0] + 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1] / 4; + const auto N15 = -(q[0] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1] / 4; + const auto N16 = -(q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[2] + 1) * q[1] * q[2] / 4; + const auto N17 = -(q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (q[2] + 1) * q[0] * q[2] / 4; + const auto N18 = -(q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] + 1) * q[0] * q[2] / 4; + const auto N19 = -(q[0] - 1) * (q[0] + 1) * (q[1] + 1) * (q[2] + 1) * q[1] * q[2] / 4; + const auto N20 = (q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * q[2] / 2; + const auto N21 = (q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] + 1) * q[2] / 2; + const auto N22 = (q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[1] / 2; + const auto N23 = (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] / 2; + const auto N24 = (q[0] - 1) * (q[0] + 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[1] / 2; + const auto N25 = (q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] / 2; + const auto N26 = + -(q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1); + + return {N0, N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, N11, N12, N13, N14, N15, N16, N17, N18, N19, N20, N21, N22, N23, N24, N25, N26}; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const type::VecNoInit N_0 = {(2 * q[0] - 1) * (q[1] - 1) * (q[2] - 1) * q[1] * q[2] / 8, + (q[0] - 1) * (2 * q[1] - 1) * (q[2] - 1) * q[0] * q[2] / 8, + (q[0] - 1) * (q[1] - 1) * (2 * q[2] - 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_1 = {(2 * q[0] + 1) * (q[1] - 1) * (q[2] - 1) * q[1] * q[2] / 8, + (q[0] + 1) * (2 * q[1] - 1) * (q[2] - 1) * q[0] * q[2] / 8, + (q[0] + 1) * (q[1] - 1) * (2 * q[2] - 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_2 = {(2 * q[0] + 1) * (q[1] + 1) * (q[2] - 1) * q[1] * q[2] / 8, + (q[0] + 1) * (2 * q[1] + 1) * (q[2] - 1) * q[0] * q[2] / 8, + (q[0] + 1) * (q[1] + 1) * (2 * q[2] - 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_3 = {(2 * q[0] - 1) * (q[1] + 1) * (q[2] - 1) * q[1] * q[2] / 8, + (q[0] - 1) * (2 * q[1] + 1) * (q[2] - 1) * q[0] * q[2] / 8, + (q[0] - 1) * (q[1] + 1) * (2 * q[2] - 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_4 = {(2 * q[0] - 1) * (q[1] - 1) * (q[2] + 1) * q[1] * q[2] / 8, + (q[0] - 1) * (2 * q[1] - 1) * (q[2] + 1) * q[0] * q[2] / 8, + (q[0] - 1) * (q[1] - 1) * (2 * q[2] + 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_5 = {(2 * q[0] + 1) * (q[1] - 1) * (q[2] + 1) * q[1] * q[2] / 8, + (q[0] + 1) * (2 * q[1] - 1) * (q[2] + 1) * q[0] * q[2] / 8, + (q[0] + 1) * (q[1] - 1) * (2 * q[2] + 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_6 = {(2 * q[0] + 1) * (q[1] + 1) * (q[2] + 1) * q[1] * q[2] / 8, + (q[0] + 1) * (2 * q[1] + 1) * (q[2] + 1) * q[0] * q[2] / 8, + (q[0] + 1) * (q[1] + 1) * (2 * q[2] + 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_7 = {(2 * q[0] - 1) * (q[1] + 1) * (q[2] + 1) * q[1] * q[2] / 8, + (q[0] - 1) * (2 * q[1] + 1) * (q[2] + 1) * q[0] * q[2] / 8, + (q[0] - 1) * (q[1] + 1) * (2 * q[2] + 1) * q[0] * q[1] / 8}; + const type::VecNoInit N_8 = { + -(q[1] - 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 2, + -(q[0] - 1) * (q[0] + 1) * (2 * q[1] - 1) * (q[2] - 1) * q[2] / 4, + -(q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (2 * q[2] - 1) * q[1] / 4}; + const type::VecNoInit N_9 = { + -(2 * q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * q[2] / 4, + -(q[0] - 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 2, + -(q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (2 * q[2] - 1) * q[0] / 4}; + const type::VecNoInit N_10 = { + -(2 * q[0] - 1) * (q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[1] / 4, + -(q[0] - 1) * (2 * q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[0] / 4, + -(q[0] - 1) * (q[1] - 1) * q[0] * q[1] * q[2] / 2}; + const type::VecNoInit N_11 = { + -(2 * q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * q[2] / 4, + -(q[0] + 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 2, + -(q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (2 * q[2] - 1) * q[0] / 4}; + const type::VecNoInit N_12 = { + -(2 * q[0] + 1) * (q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[1] / 4, + -(q[0] + 1) * (2 * q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[0] / 4, + -(q[0] + 1) * (q[1] - 1) * q[0] * q[1] * q[2] / 2}; + const type::VecNoInit N_13 = { + -(q[1] + 1) * (q[2] - 1) * q[0] * q[1] * q[2] / 2, + -(q[0] - 1) * (q[0] + 1) * (2 * q[1] + 1) * (q[2] - 1) * q[2] / 4, + -(q[0] - 1) * (q[0] + 1) * (q[1] + 1) * (2 * q[2] - 1) * q[1] / 4}; + const type::VecNoInit N_14 = { + -(2 * q[0] + 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[1] / 4, + -(q[0] + 1) * (2 * q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] / 4, + -(q[0] + 1) * (q[1] + 1) * q[0] * q[1] * q[2] / 2}; + const type::VecNoInit N_15 = { + -(2 * q[0] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[1] / 4, + -(q[0] - 1) * (2 * q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] / 4, + -(q[0] - 1) * (q[1] + 1) * q[0] * q[1] * q[2] / 2}; + const type::VecNoInit N_16 = { + -(q[1] - 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 2, + -(q[0] - 1) * (q[0] + 1) * (2 * q[1] - 1) * (q[2] + 1) * q[2] / 4, + -(q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (2 * q[2] + 1) * q[1] / 4}; + const type::VecNoInit N_17 = { + -(2 * q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (q[2] + 1) * q[2] / 4, + -(q[0] - 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 2, + -(q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (2 * q[2] + 1) * q[0] / 4}; + const type::VecNoInit N_18 = { + -(2 * q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] + 1) * q[2] / 4, + -(q[0] + 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 2, + -(q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (2 * q[2] + 1) * q[0] / 4}; + const type::VecNoInit N_19 = { + -(q[1] + 1) * (q[2] + 1) * q[0] * q[1] * q[2] / 2, + -(q[0] - 1) * (q[0] + 1) * (2 * q[1] + 1) * (q[2] + 1) * q[2] / 4, + -(q[0] - 1) * (q[0] + 1) * (q[1] + 1) * (2 * q[2] + 1) * q[1] / 4}; + const type::VecNoInit N_20 = { + (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * q[0] * q[2], + (q[0] - 1) * (q[0] + 1) * (q[2] - 1) * q[1] * q[2], + (q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (2 * q[2] - 1) / 2}; + const type::VecNoInit N_21 = { + (q[1] - 1) * (q[1] + 1) * (q[2] + 1) * q[0] * q[2], + (q[0] - 1) * (q[0] + 1) * (q[2] + 1) * q[1] * q[2], + (q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (2 * q[2] + 1) / 2}; + const type::VecNoInit N_22 = { + (q[1] - 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1], + (q[0] - 1) * (q[0] + 1) * (2 * q[1] - 1) * (q[2] - 1) * (q[2] + 1) / 2, + (q[0] - 1) * (q[0] + 1) * (q[1] - 1) * q[1] * q[2]}; + const type::VecNoInit N_23 = { + (2 * q[0] + 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) / 2, + (q[0] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1], + (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * q[0] * q[2]}; + const type::VecNoInit N_24 = { + (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1], + (q[0] - 1) * (q[0] + 1) * (2 * q[1] + 1) * (q[2] - 1) * (q[2] + 1) / 2, + (q[0] - 1) * (q[0] + 1) * (q[1] + 1) * q[1] * q[2]}; + const type::VecNoInit N_25 = { + (2 * q[0] - 1) * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) / 2, + (q[0] - 1) * (q[2] - 1) * (q[2] + 1) * q[0] * q[1], + (q[0] - 1) * (q[1] - 1) * (q[1] + 1) * q[0] * q[2]}; + const type::VecNoInit N_26 = { + -2 * (q[1] - 1) * (q[1] + 1) * (q[2] - 1) * (q[2] + 1) * q[0], + -2 * (q[0] - 1) * (q[0] + 1) * (q[2] - 1) * (q[2] + 1) * q[1], + -2 * (q[0] - 1) * (q[0] + 1) * (q[1] - 1) * (q[1] + 1) * q[2]}; + + return sofa::type::Mat( + N_0, N_1, N_2, N_3, N_4, N_5, N_6, N_7, N_8, N_9, N_10, N_11, N_12, N_13, N_14, N_15, + N_16, N_17, N_18, N_19, N_20, N_21, N_22, N_23, N_24, N_25, N_26); + } + + static constexpr std::array quadraturePoints() + { + constexpr Real a = 0.7745966692414834; // sqrt(3/5) + constexpr Real w_a = 5.0 / 9.0; + constexpr Real w_0 = 8.0 / 9.0; + + constexpr std::array pts = {-a, 0.0, a}; + constexpr std::array wts = {w_a, w_0, w_a}; + + std::array q {}; + int index = 0; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + for (int k = 0; k < 3; ++k) + { + q[index++] = std::make_pair( + sofa::type::Vec(pts[i], pts[j], pts[k]), + wts[i] * wts[j] * wts[k] + ); + } + } + } + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_HEXAHEDRON_CPP) +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPrism].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPrism].cpp new file mode 100644 index 00000000000..9561494cab0 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPrism].cpp @@ -0,0 +1,31 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_QUADRATIC_PRISM_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPrism].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPrism].h new file mode 100644 index 00000000000..1a0490e76b5 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPrism].h @@ -0,0 +1,219 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_PRISM_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::QuadraticPrism, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Quadratic Prisms are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + // 6 corner vertices + {0, 0, -1}, // vertex 0 + {1, 0, -1}, // vertex 1 + {0, 1, -1}, // vertex 2 + {0, 0, 1}, // vertex 3 + {1, 0, 1}, // vertex 4 + {0, 1, 1}, // vertex 5 + // 9 mid-edge nodes + {0.5, 0, -1}, // mid-edge 0-1 + {0.5, 0.5, -1}, // mid-edge 1-2 + {0, 0.5, -1}, // mid-edge 2-0 + {0, 0, 0}, // mid-edge 0-3 + {1, 0, 0}, // mid-edge 1-4 + {0, 1, 0}, // mid-edge 2-5 + {0.5, 0, 1}, // mid-edge 3-4 + {0.5, 0.5, 1}, // mid-edge 4-5 + {0, 0.5, 1}, // mid-edge 5-3 + // 3 face-center nodes (for the 3 quad faces) + {0.5, 0, 0}, // face center (0-1-4-3) + {0.5, 0.5, 0}, // face center (1-2-5-4) + {0, 0.5, 0} // face center (2-0-3-5) + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getElements(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + const Real l0 = static_cast(1) - q[0] - q[1]; + const Real l1 = q[0]; + const Real l2 = q[1]; + const Real zeta = q[2]; + + // Corner vertices (6) + const Real N0 = 0.125 * l0 * (2 * l0 - 1) * (zeta * zeta - zeta); + const Real N1 = 0.125 * l1 * (2 * l1 - 1) * (zeta * zeta - zeta); + const Real N2 = 0.125 * l2 * (2 * l2 - 1) * (zeta * zeta - zeta); + const Real N3 = 0.125 * l0 * (2 * l0 - 1) * (zeta * zeta + zeta); + const Real N4 = 0.125 * l1 * (2 * l1 - 1) * (zeta * zeta + zeta); + const Real N5 = 0.125 * l2 * (2 * l2 - 1) * (zeta * zeta + zeta); + + // Mid-edge nodes on triangular faces (6) + const Real N6 = 0.5 * l0 * l1 * (zeta * zeta - zeta); + const Real N7 = 0.5 * l1 * l2 * (zeta * zeta - zeta); + const Real N8 = 0.5 * l2 * l0 * (zeta * zeta - zeta); + const Real N12 = 0.5 * l0 * l1 * (zeta * zeta + zeta); + const Real N13 = 0.5 * l1 * l2 * (zeta * zeta + zeta); + const Real N14 = 0.5 * l2 * l0 * (zeta * zeta + zeta); + + // Mid-edge nodes connecting triangular faces (3) + const Real N9 = 0.25 * l0 * (2 * l0 - 1) * (1 - zeta * zeta); + const Real N10 = 0.25 * l1 * (2 * l1 - 1) * (1 - zeta * zeta); + const Real N11 = 0.25 * l2 * (2 * l2 - 1) * (1 - zeta * zeta); + + // Face-center nodes on quad faces (3) + const Real N15 = 0.5 * l0 * l1 * (1 - zeta * zeta); + const Real N16 = 0.5 * l1 * l2 * (1 - zeta * zeta); + const Real N17 = 0.5 * l2 * l0 * (1 - zeta * zeta); + + return {N0, N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, N11, N12, N13, N14, N15, N16, N17}; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const Real l0 = static_cast(1) - q[0] - q[1]; + const Real l1 = q[0]; + const Real l2 = q[1]; + const Real zeta = q[2]; + + return { + // Corner vertices (6) + {0.125 * (-4 * l0 + 1) * (zeta * zeta - zeta), + 0.125 * (-4 * l0 + 1) * (zeta * zeta - zeta), + 0.125 * l0 * (2 * l0 - 1) * (2 * zeta - 1)}, + + {0.125 * (4 * l1 - 1) * (zeta * zeta - zeta), + 0, + 0.125 * l1 * (2 * l1 - 1) * (2 * zeta - 1)}, + + {0, + 0.125 * (4 * l2 - 1) * (zeta * zeta - zeta), + 0.125 * l2 * (2 * l2 - 1) * (2 * zeta - 1)}, + + {0.125 * (-4 * l0 + 1) * (zeta * zeta + zeta), + 0.125 * (-4 * l0 + 1) * (zeta * zeta + zeta), + 0.125 * l0 * (2 * l0 - 1) * (2 * zeta + 1)}, + + {0.125 * (4 * l1 - 1) * (zeta * zeta + zeta), + 0, + 0.125 * l1 * (2 * l1 - 1) * (2 * zeta + 1)}, + + {0, + 0.125 * (4 * l2 - 1) * (zeta * zeta + zeta), + 0.125 * l2 * (2 * l2 - 1) * (2 * zeta + 1)}, + + // Mid-edge nodes on triangular faces (6) + {0.5 * (l1 - l0) * (zeta * zeta - zeta), + -0.5 * l1 * (zeta * zeta - zeta), + 0.5 * l0 * l1 * (2 * zeta - 1)}, + + {0.5 * l2 * (zeta * zeta - zeta), + 0.5 * l1 * (zeta * zeta - zeta), + 0.5 * l1 * l2 * (2 * zeta - 1)}, + + {-0.5 * l2 * (zeta * zeta - zeta), + 0.5 * (l0 - l2) * (zeta * zeta - zeta), + 0.5 * l2 * l0 * (2 * zeta - 1)}, + + {0.25 * (-4 * l0 + 1) * (1 - zeta * zeta), + 0.25 * (-4 * l0 + 1) * (1 - zeta * zeta), + -0.5 * l0 * (2 * l0 - 1) * zeta}, + + {0.25 * (4 * l1 - 1) * (1 - zeta * zeta), + 0, + -0.5 * l1 * (2 * l1 - 1) * zeta}, + + {0, + 0.25 * (4 * l2 - 1) * (1 - zeta * zeta), + -0.5 * l2 * (2 * l2 - 1) * zeta}, + + {0.5 * (l1 - l0) * (zeta * zeta + zeta), + -0.5 * l1 * (zeta * zeta + zeta), + 0.5 * l0 * l1 * (2 * zeta + 1)}, + + {0.5 * l2 * (zeta * zeta + zeta), + 0.5 * l1 * (zeta * zeta + zeta), + 0.5 * l1 * l2 * (2 * zeta + 1)}, + + {-0.5 * l2 * (zeta * zeta + zeta), + 0.5 * (l0 - l2) * (zeta * zeta + zeta), + 0.5 * l2 * l0 * (2 * zeta + 1)}, + + // Face-center nodes on quad faces (3) + {0.5 * (l1 - l0) * (1 - zeta * zeta), + -0.5 * l1 * (1 - zeta * zeta), + -l0 * l1 * zeta}, + + {0.5 * l2 * (1 - zeta * zeta), + 0.5 * l1 * (1 - zeta * zeta), + -l1 * l2 * zeta}, + + {-0.5 * l2 * (1 - zeta * zeta), + 0.5 * (l0 - l2) * (1 - zeta * zeta), + -l2 * l0 * zeta} + }; + } + + static constexpr std::array quadraturePoints() + { + constexpr Real a = 2. / 3.; + constexpr Real b = 1. / 6.; + // constexpr Real c = 1. / std::sqrt(3.); + constexpr Real c { 0.57735026919 }; + constexpr Real w = 1. / 6.; + + constexpr sofa::type::Vec q0(b, b, -c); + constexpr sofa::type::Vec q1(a, b, -c); + constexpr sofa::type::Vec q2(b, a, -c); + constexpr sofa::type::Vec q3(b, b, c); + constexpr sofa::type::Vec q4(a, b, c); + constexpr sofa::type::Vec q5(b, a, c); + + constexpr std::array q { + std::make_pair(q0, w), + std::make_pair(q1, w), + std::make_pair(q2, w), + std::make_pair(q3, w), + std::make_pair(q4, w), + std::make_pair(q5, w) + }; + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_PRISM_CPP) +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPyramid].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPyramid].cpp new file mode 100644 index 00000000000..bb6af354ec7 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPyramid].cpp @@ -0,0 +1,31 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_QUADRATIC_PYRAMID_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPyramid].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPyramid].h new file mode 100644 index 00000000000..39aefe7f772 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticPyramid].h @@ -0,0 +1,192 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_PYRAMID_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::QuadraticPyramid, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Quadratic Pyramids are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + // 5 corner vertices + {-1, -1, 0}, // vertex 0 + {1, -1, 0}, // vertex 1 + {1, 1, 0}, // vertex 2 + {-1, 1, 0}, // vertex 3 + {0, 0, 1}, // vertex 4 (apex) + // 8 mid-edge nodes + {0, -1, 0}, // mid-edge 0-1 + {1, 0, 0}, // mid-edge 1-2 + {0, 1, 0}, // mid-edge 2-3 + {-1, 0, 0}, // mid-edge 3-0 + {-0.5, -0.5, 0.5}, // mid-edge 0-4 + {0.5, -0.5, 0.5}, // mid-edge 1-4 + {0.5, 0.5, 0.5}, // mid-edge 2-4 + {-0.5, 0.5, 0.5}, // mid-edge 3-4 + // 1 face-center node (on the quad base) + {0, 0, 0} // face center of base quad + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getElements(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + const Real xi = q[0]; + const Real eta = q[1]; + const Real zeta = q[2]; + const Real r = (zeta < 0.999) ? (1 - zeta) : 0.001; // Avoid singularity at apex + + // Corner vertices (5) + const Real N0 = 0.125 * (xi * xi - xi) * (eta * eta - eta) / r; + const Real N1 = 0.125 * (xi * xi + xi) * (eta * eta - eta) / r; + const Real N2 = 0.125 * (xi * xi + xi) * (eta * eta + eta) / r; + const Real N3 = 0.125 * (xi * xi - xi) * (eta * eta + eta) / r; + const Real N4 = zeta * (2 * zeta - 1); + + // Mid-edge nodes on base (4) + const Real N5 = 0.25 * (1 - xi * xi) * (eta * eta - eta) / r; + const Real N6 = 0.25 * (xi * xi + xi) * (1 - eta * eta) / r; + const Real N7 = 0.25 * (1 - xi * xi) * (eta * eta + eta) / r; + const Real N8 = 0.25 * (xi * xi - xi) * (1 - eta * eta) / r; + + // Mid-edge nodes connecting to apex (4) + const Real N9 = 0.5 * (xi * xi - xi) * (eta * eta - eta) * zeta / r; + const Real N10 = 0.5 * (xi * xi + xi) * (eta * eta - eta) * zeta / r; + const Real N11 = 0.5 * (xi * xi + xi) * (eta * eta + eta) * zeta / r; + const Real N12 = 0.5 * (xi * xi - xi) * (eta * eta + eta) * zeta / r; + + // Face-center node on base (1) + const Real N13 = 0.5 * (1 - xi * xi) * (1 - eta * eta) / r; + + return {N0, N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, N11, N12, N13}; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const Real xi = q[0]; + const Real eta = q[1]; + const Real zeta = q[2]; + const Real r = (zeta < 0.999) ? (1 - zeta) : 0.001; + + return { + // Corner vertices (5) + {0.125 * (2 * xi - 1) * (eta * eta - eta) / r, + 0.125 * (xi * xi - xi) * (2 * eta - 1) / r, + 0.125 * (xi * xi - xi) * (eta * eta - eta) / (r * r)}, + + {0.125 * (2 * xi + 1) * (eta * eta - eta) / r, + 0.125 * (xi * xi + xi) * (2 * eta - 1) / r, + 0.125 * (xi * xi + xi) * (eta * eta - eta) / (r * r)}, + + {0.125 * (2 * xi + 1) * (eta * eta + eta) / r, + 0.125 * (xi * xi + xi) * (2 * eta + 1) / r, + 0.125 * (xi * xi + xi) * (eta * eta + eta) / (r * r)}, + + {0.125 * (2 * xi - 1) * (eta * eta + eta) / r, + 0.125 * (xi * xi - xi) * (2 * eta + 1) / r, + 0.125 * (xi * xi - xi) * (eta * eta + eta) / (r * r)}, + + {0, 0, 4 * zeta - 1}, + + // Mid-edge nodes on base (4) + {-0.5 * xi * (eta * eta - eta) / r, + 0.25 * (1 - xi * xi) * (2 * eta - 1) / r, + 0.25 * (1 - xi * xi) * (eta * eta - eta) / (r * r)}, + + {0.25 * (2 * xi + 1) * (1 - eta * eta) / r, + -0.5 * (xi * xi + xi) * eta / r, + 0.25 * (xi * xi + xi) * (1 - eta * eta) / (r * r)}, + + {-0.5 * xi * (eta * eta + eta) / r, + 0.25 * (1 - xi * xi) * (2 * eta + 1) / r, + 0.25 * (1 - xi * xi) * (eta * eta + eta) / (r * r)}, + + {0.25 * (2 * xi - 1) * (1 - eta * eta) / r, + -0.5 * (xi * xi - xi) * eta / r, + 0.25 * (xi * xi - xi) * (1 - eta * eta) / (r * r)}, + + // Mid-edge nodes connecting to apex (4) + {0.5 * (2 * xi - 1) * (eta * eta - eta) * zeta / r, + 0.5 * (xi * xi - xi) * (2 * eta - 1) * zeta / r, + 0.5 * (xi * xi - xi) * (eta * eta - eta) * (r - zeta) / (r * r)}, + + {0.5 * (2 * xi + 1) * (eta * eta - eta) * zeta / r, + 0.5 * (xi * xi + xi) * (2 * eta - 1) * zeta / r, + 0.5 * (xi * xi + xi) * (eta * eta - eta) * (r - zeta) / (r * r)}, + + {0.5 * (2 * xi + 1) * (eta * eta + eta) * zeta / r, + 0.5 * (xi * xi + xi) * (2 * eta + 1) * zeta / r, + 0.5 * (xi * xi + xi) * (eta * eta + eta) * (r - zeta) / (r * r)}, + + {0.5 * (2 * xi - 1) * (eta * eta + eta) * zeta / r, + 0.5 * (xi * xi - xi) * (2 * eta + 1) * zeta / r, + 0.5 * (xi * xi - xi) * (eta * eta + eta) * (r - zeta) / (r * r)}, + + // Face-center node on base (1) + {-xi * (1 - eta * eta) / r, + -eta * (1 - xi * xi) / r, + 0.5 * (1 - xi * xi) * (1 - eta * eta) / (r * r)} + }; + } + + static constexpr std::array quadraturePoints() + { + // 5-point quadrature for pyramid + constexpr Real a = 0.584237394672974; + constexpr Real b = 0.138196601125011; + constexpr Real w1 = 0.133333333333333; + constexpr Real w2 = 0.075000000000000; + + constexpr sofa::type::Vec q0(0, 0, b); + constexpr sofa::type::Vec q1(a, 0, 0.5); + constexpr sofa::type::Vec q2(0, a, 0.5); + constexpr sofa::type::Vec q3(-a, 0, 0.5); + constexpr sofa::type::Vec q4(0, -a, 0.5); + + constexpr std::array q { + std::make_pair(q0, w1), + std::make_pair(q1, w2), + std::make_pair(q2, w2), + std::make_pair(q3, w2), + std::make_pair(q4, w2) + }; + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_PYRAMID_CPP) +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].cpp new file mode 100644 index 00000000000..b6b68400801 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].cpp @@ -0,0 +1,32 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_QUADRATIC_QUAD_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h new file mode 100644 index 00000000000..e1f312c21f5 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h @@ -0,0 +1,125 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_QUAD_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2); + + constexpr static std::array referenceElementNodes {{ + {-1, -1}, // vertex 0 + {1, -1}, // vertex 1 + {1, 1}, // vertex 2 + {-1, 1}, // vertex 3 + {0, -1}, // mid-edge 0-1 + {1, 0}, // mid-edge 1-2 + {0, 1}, // mid-edge 2-3 + {-1, 0}, // mid-edge 3-0 + {0, 0} // face center + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getElements(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + const Real xi = q[0]; + const Real eta = q[1]; + + return { + 0.25 * (xi * xi - xi) * (eta * eta - eta), // vertex 0 + 0.25 * (xi * xi + xi) * (eta * eta - eta), // vertex 1 + 0.25 * (xi * xi + xi) * (eta * eta + eta), // vertex 2 + 0.25 * (xi * xi - xi) * (eta * eta + eta), // vertex 3 + 0.5 * (1 - xi * xi) * (eta * eta - eta), // mid-edge 0-1 + 0.5 * (xi * xi + xi) * (1 - eta * eta), // mid-edge 1-2 + 0.5 * (1 - xi * xi) * (eta * eta + eta), // mid-edge 2-3 + 0.5 * (xi * xi - xi) * (1 - eta * eta), // mid-edge 3-0 + (1 - xi * xi) * (1 - eta * eta) // face center + }; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const Real xi = q[0]; + const Real eta = q[1]; + + return { + // vertex 0 + {0.25 * (2 * xi - 1) * (eta * eta - eta), 0.25 * (xi * xi - xi) * (2 * eta - 1)}, + // vertex 1 + {0.25 * (2 * xi + 1) * (eta * eta - eta), 0.25 * (xi * xi + xi) * (2 * eta - 1)}, + // vertex 2 + {0.25 * (2 * xi + 1) * (eta * eta + eta), 0.25 * (xi * xi + xi) * (2 * eta + 1)}, + // vertex 3 + {0.25 * (2 * xi - 1) * (eta * eta + eta), 0.25 * (xi * xi - xi) * (2 * eta + 1)}, + // mid-edge 0-1 + {-xi * (eta * eta - eta), 0.5 * (1 - xi * xi) * (2 * eta - 1)}, + // mid-edge 1-2 + {0.5 * (2 * xi + 1) * (1 - eta * eta), -(xi * xi + xi) * eta}, + // mid-edge 2-3 + {-xi * (eta * eta + eta), 0.5 * (1 - xi * xi) * (2 * eta + 1)}, + // mid-edge 3-0 + {0.5 * (2 * xi - 1) * (1 - eta * eta), -(xi * xi - xi) * eta}, + // face center + {-2 * xi * (1 - eta * eta), -2 * eta * (1 - xi * xi)} + }; + } + + static constexpr std::array quadraturePoints() + { + // constexpr Real a = 1. / std::sqrt(3.); + constexpr Real a { 0.57735026919 }; + constexpr Real w = 1.0; + + constexpr sofa::type::Vec q0(-a, -a); + constexpr sofa::type::Vec q1(a, -a); + constexpr sofa::type::Vec q2(a, a); + constexpr sofa::type::Vec q3(-a, a); + + constexpr std::array q { + std::make_pair(q0, w), + std::make_pair(q1, w), + std::make_pair(q2, w), + std::make_pair(q3, w) + }; + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_QUAD_CPP) +extern template struct SOFA_FEM_API FiniteElement; +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].cpp new file mode 100644 index 00000000000..71b32966867 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].cpp @@ -0,0 +1,31 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TETAHEDRON_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h new file mode 100644 index 00000000000..2cd90ee677e --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h @@ -0,0 +1,136 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TETAHEDRON_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::QuadraticTetrahedron, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Quadratic Tetrahedrons are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + {0, 0, 0}, // vertex 0 + {1, 0, 0}, // vertex 1 + {0, 1, 0}, // vertex 2 + {0, 0, 1}, // vertex 3 + {0.5, 0, 0}, // mid-edge 0-1 + {0, 0.5, 0}, // mid-edge 0-2 + {0, 0, 0.5}, // mid-edge 0-3 + {0.5, 0.5, 0}, // mid-edge 1-2 + {0.5, 0, 0.5}, // mid-edge 1-3 + {0, 0.5, 0.5} // mid-edge 2-3 + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getElements(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + const Real l0 = static_cast(1) - q[0] - q[1] - q[2]; + const Real l1 = q[0]; + const Real l2 = q[1]; + const Real l3 = q[2]; + + return { + l0 * (2 * l0 - 1), // vertex 0 + l1 * (2 * l1 - 1), // vertex 1 + l2 * (2 * l2 - 1), // vertex 2 + l3 * (2 * l3 - 1), // vertex 3 + 4 * l0 * l1, // mid-edge 0-1 + 4 * l0 * l2, // mid-edge 0-2 + 4 * l0 * l3, // mid-edge 0-3 + 4 * l1 * l2, // mid-edge 1-2 + 4 * l1 * l3, // mid-edge 1-3 + 4 * l2 * l3 // mid-edge 2-3 + }; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const Real l0 = static_cast(1) - q[0] - q[1] - q[2]; + const Real l1 = q[0]; + const Real l2 = q[1]; + const Real l3 = q[2]; + + return { + // vertex 0: derivatives with respect to q[0], q[1], q[2] + {-3 + 4 * (q[0] + q[1] + q[2]), -3 + 4 * (q[0] + q[1] + q[2]), -3 + 4 * (q[0] + q[1] + q[2])}, + // vertex 1 + {4 * q[0] - 1, 0, 0}, + // vertex 2 + {0, 4 * q[1] - 1, 0}, + // vertex 3 + {0, 0, 4 * q[2] - 1}, + // mid-edge 0-1 + {4 * (l0 - l1), -4 * l1, -4 * l1}, + // mid-edge 0-2 + {-4 * l2, 4 * (l0 - l2), -4 * l2}, + // mid-edge 0-3 + {-4 * l3, -4 * l3, 4 * (l0 - l3)}, + // mid-edge 1-2 + {4 * l2, 4 * l1, 0}, + // mid-edge 1-3 + {4 * l3, 0, 4 * l1}, + // mid-edge 2-3 + {0, 4 * l3, 4 * l2} + }; + } + + static constexpr std::array quadraturePoints() + { + // 4-point quadrature rule for quadratic tetrahedron + // constexpr Real a = (5. + 3. * std::sqrt(5.)) / 20.; + constexpr Real a { 0.585410196625 }; + // constexpr Real b = (5. - std::sqrt(5.)) / 20.; + constexpr Real b { 0.138196601125 }; + constexpr Real w = 1. / 24.; + + constexpr sofa::type::Vec q0(a, b, b); + constexpr sofa::type::Vec q1(b, a, b); + constexpr sofa::type::Vec q2(b, b, a); + constexpr sofa::type::Vec q3(b, b, b); + + constexpr std::array q { + std::make_pair(q0, w), + std::make_pair(q1, w), + std::make_pair(q2, w), + std::make_pair(q3, w) + }; + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TETAHEDRON_CPP) +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].cpp new file mode 100644 index 00000000000..d6aa7153971 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].cpp @@ -0,0 +1,32 @@ +/****************************************************************************** +* 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 * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TRIANGLE_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h new file mode 100644 index 00000000000..c77ebd73cac --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h @@ -0,0 +1,113 @@ +/****************************************************************************** +* 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 + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TRIANGLE_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::QuadraticTriangle, DataTypes, 2); + + constexpr static std::array referenceElementNodes {{ + {0, 0}, // vertex 0 + {1, 0}, // vertex 1 + {0, 1}, // vertex 2 + {0.5, 0}, // mid-edge 0-1 + {0.5, 0.5}, // mid-edge 1-2 + {0, 0.5} // mid-edge 2-0 + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getElements(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + const Real l0 = static_cast(1) - q[0] - q[1]; + const Real l1 = q[0]; + const Real l2 = q[1]; + + return { + l0 * (2 * l0 - 1), // vertex 0 + l1 * (2 * l1 - 1), // vertex 1 + l2 * (2 * l2 - 1), // vertex 2 + 4 * l0 * l1, // mid-edge 0-1 + 4 * l1 * l2, // mid-edge 1-2 + 4 * l2 * l0 // mid-edge 2-0 + }; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + const Real l0 = static_cast(1) - q[0] - q[1]; + const Real l1 = q[0]; + const Real l2 = q[1]; + + return { + // vertex 0 + {-3 + 4 * (q[0] + q[1]), -3 + 4 * (q[0] + q[1])}, + // vertex 1 + {4 * q[0] - 1, 0}, + // vertex 2 + {0, 4 * q[1] - 1}, + // mid-edge 0-1 + {4 * (l0 - l1), -4 * l1}, + // mid-edge 1-2 + {4 * l2, 4 * l1}, + // mid-edge 2-0 + {-4 * l2, 4 * (l0 - l2)} + }; + } + + static constexpr std::array quadraturePoints() + { + constexpr Real a = 2. / 3.; + constexpr Real b = 1. / 6.; + constexpr Real w = 1. / 6.; + + constexpr sofa::type::Vec q0(b, b); + constexpr sofa::type::Vec q1(a, b); + constexpr sofa::type::Vec q2(b, a); + + constexpr std::array q { + std::make_pair(q0, w), + std::make_pair(q1, w), + std::make_pair(q2, w) + }; + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TRIANGLE_CPP) +extern template struct SOFA_FEM_API FiniteElement; +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h index 163df52198c..ffd3dace02b 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h @@ -28,3 +28,11 @@ #include #include #include + +#include +#include +#include +#include +#include +#include +#include diff --git a/Sofa/framework/FEM/test/FiniteElement_test.cpp b/Sofa/framework/FEM/test/FiniteElement_test.cpp index 68744ae5f85..0b2100347c3 100644 --- a/Sofa/framework/FEM/test/FiniteElement_test.cpp +++ b/Sofa/framework/FEM/test/FiniteElement_test.cpp @@ -93,6 +93,74 @@ TEST(FiniteElement, pyramid3dWeights) testSumWeights(8); } +TEST(FiniteElement, quadraticHexa3dWeights) +{ + testSumWeights(8); +} + +/** + * Checks that the shape function i at the node j is equal to delta_ij + */ +template +void testShapeFunctionsAtNodes() +{ + using FE = sofa::fem::FiniteElement; + + for (unsigned int j = 0; j < FE::NumberOfNodesInElement; ++j) + { + const auto& node_j = FE::referenceElementNodes[j]; + const auto N = FE::shapeFunctions(node_j); + + for (unsigned int i = 0; i < FE::NumberOfNodesInElement; ++i) + { + if (i == j) + EXPECT_NEAR(N[i], 1.0, 1e-12) << "Shape function " << i << " should be 1 at node " << j; + else + EXPECT_NEAR(N[i], 0.0, 1e-12) << "Shape function " << i << " should be 0 at node " << j; + } + } +} + +TEST(FiniteElement, edge1dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + +TEST(FiniteElement, triangle2dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + +TEST(FiniteElement, quad2dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + +TEST(FiniteElement, tetra3dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + +TEST(FiniteElement, hexa3dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + +TEST(FiniteElement, prism3dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + +TEST(FiniteElement, pyramid3dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + +TEST(FiniteElement, quadraticHexa3dShapeFunctionsAtNodes) +{ + testShapeFunctionsAtNodes(); +} + /** * Checks that the sum of the gradients of shape functions is zero at the evaluation point. */ @@ -180,4 +248,173 @@ TEST(FiniteElement, pyramid3dGradientShapeFunctions) testGradientShapeFunctions(sofa::type::Vec3(0.5, 0.5, 0.5)); } +TEST(FiniteElement, quadraticHexa3dGradientShapeFunctions) +{ + testGradientShapeFunctions(sofa::type::Vec3(0., 0., 0.)); + testGradientShapeFunctions(sofa::type::Vec3(1,1,1)); +} + +// Partition of unity at arbitrary (non-node) points, not just at nodes +template +void testPartitionOfUnity() +{ + using FE = sofa::fem::FiniteElement; + using Real = sofa::Real_t; + + for (const auto& [q, w] : FE::quadraturePoints()) + { + const auto N = FE::shapeFunctions(q); + Real sum = 0; + for (unsigned int i = 0; i < FE::NumberOfNodesInElement; ++i) + sum += N[i]; + EXPECT_NEAR(sum, Real(1), 1e-10) << "Shape functions must sum to 1 everywhere"; + } +} + +TEST(FiniteElement, edge1dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +TEST(FiniteElement, triangle2dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +TEST(FiniteElement, quad2dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +TEST(FiniteElement, tetra3dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +TEST(FiniteElement, hexa3dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +TEST(FiniteElement, prism3dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +TEST(FiniteElement, pyramid3dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +TEST(FiniteElement, quadraticHexa3dPartitionOfUnity) +{ + testPartitionOfUnity(); +} + +// Finite-difference cross-check: analytical gradient vs numerical gradient +// Catches algebra/sign errors without needing hand-derived expected values. +template +void testGradientMatchesFiniteDifference( + const sofa::type::Vec::TopologicalDimension, sofa::Real_t>& q) +{ + using FE = sofa::fem::FiniteElement; + using Real = sofa::Real_t; + constexpr Real h = 1e-6; + + const auto analyticalGrad = FE::gradientShapeFunctions(q); + + for (unsigned int d = 0; d < FE::TopologicalDimension; ++d) + { + auto qPlus = q, qMinus = q; + qPlus[d] += h; + qMinus[d] -= h; + + const auto Nplus = FE::shapeFunctions(qPlus); + const auto Nminus = FE::shapeFunctions(qMinus); + + for (unsigned int i = 0; i < FE::NumberOfNodesInElement; ++i) + { + const Real fd = (Nplus[i] - Nminus[i]) / (2 * h); + EXPECT_NEAR(analyticalGrad[i][d], fd, 1e-5) + << "Gradient mismatch for node " << i << ", direction " << d; + } + } +} + + +TEST(FiniteElement, edge1dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec1(0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec1(1.)); + testGradientMatchesFiniteDifference(sofa::type::Vec1(-1.)); +} + +TEST(FiniteElement, edge2dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec1(0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec1(1.)); + testGradientMatchesFiniteDifference(sofa::type::Vec1(-1.)); +} + +TEST(FiniteElement, edge3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec1(0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec1(1.)); + testGradientMatchesFiniteDifference(sofa::type::Vec1(-1.)); +} + +TEST(FiniteElement, triangle2dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec2(0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec2(1., 1.)); +} + +TEST(FiniteElement, triangle3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec2(0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec2(1., 1.)); +} + +TEST(FiniteElement, quad2dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec2(0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec2(1., 1.)); +} + +TEST(FiniteElement, quad3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec2(0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec2(1., 1.)); +} + +TEST(FiniteElement, tetra3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec3(0., 0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec3(1., 1., 1.)); +} + +TEST(FiniteElement, hexa3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec3(0., 0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec3(1., 1., 1.)); +} + +TEST(FiniteElement, prism3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec3(0., 0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec3(1., 1., 1.)); +} + +TEST(FiniteElement, pyramid3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec3(0., 0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec3(0.5, 0.5, 0.5)); +} + +TEST(FiniteElement, quadraticHexa3dGradientMatchesFiniteDifference) +{ + testGradientMatchesFiniteDifference(sofa::type::Vec3(0., 0., 0.)); + testGradientMatchesFiniteDifference(sofa::type::Vec3(1, 1, 1)); +} + } diff --git a/Sofa/framework/Geometry/CMakeLists.txt b/Sofa/framework/Geometry/CMakeLists.txt index 311e613ff22..e92ec8ac4f0 100644 --- a/Sofa/framework/Geometry/CMakeLists.txt +++ b/Sofa/framework/Geometry/CMakeLists.txt @@ -15,6 +15,7 @@ set(HEADER_FILES ${SOFAGEOMETRYSRC_ROOT}/Prism.h ${SOFAGEOMETRYSRC_ROOT}/Pyramid.h ${SOFAGEOMETRYSRC_ROOT}/Quad.h + ${SOFAGEOMETRYSRC_ROOT}/QuadraticElements.h ${SOFAGEOMETRYSRC_ROOT}/Tetrahedron.h ${SOFAGEOMETRYSRC_ROOT}/Triangle.h diff --git a/Sofa/framework/Geometry/src/sofa/geometry/ElementType.h b/Sofa/framework/Geometry/src/sofa/geometry/ElementType.h index 57ece622444..d850562e0df 100644 --- a/Sofa/framework/Geometry/src/sofa/geometry/ElementType.h +++ b/Sofa/framework/Geometry/src/sofa/geometry/ElementType.h @@ -30,7 +30,9 @@ namespace sofa::geometry enum class ElementType : sofa::Size { UNKNOWN, + POINT, + EDGE, TRIANGLE, QUAD, @@ -38,6 +40,15 @@ enum class ElementType : sofa::Size HEXAHEDRON, PRISM, PYRAMID, + + QUADRATIC_EDGE, + QUADRATIC_TRIANGLE, + QUADRATIC_QUAD, + QUADRATIC_TETRAHEDRON, + QUADRATIC_HEXAHEDRON, + QUADRATIC_PRISM, + QUADRATIC_PYRAMID, + SIZE }; @@ -55,6 +66,15 @@ constexpr const char* elementTypeToString(ElementType type) case ElementType::HEXAHEDRON: { return "Hexahedron"; } case ElementType::PRISM: { return "Prism"; } case ElementType::PYRAMID: { return "Pyramid"; } + + case ElementType::QUADRATIC_EDGE: { return "Quadratic Edge"; } + case ElementType::QUADRATIC_TRIANGLE: { return "Quadratic Triangle"; } + case ElementType::QUADRATIC_QUAD: { return "Quadratic Quad"; } + case ElementType::QUADRATIC_TETRAHEDRON: { return "Quadratic Tetrahedron"; } + case ElementType::QUADRATIC_HEXAHEDRON: { return "Quadratic Hexahedron"; } + case ElementType::QUADRATIC_PRISM: { return "Quadratic Prism"; } + case ElementType::QUADRATIC_PYRAMID: { return "Quadratic Pyramid"; } + default: return "Unknown"; } diff --git a/Sofa/framework/Geometry/src/sofa/geometry/QuadraticElements.h b/Sofa/framework/Geometry/src/sofa/geometry/QuadraticElements.h new file mode 100644 index 00000000000..03913448bb4 --- /dev/null +++ b/Sofa/framework/Geometry/src/sofa/geometry/QuadraticElements.h @@ -0,0 +1,79 @@ +/****************************************************************************** +* 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 + +namespace sofa::geometry +{ + +struct QuadraticEdge +{ + static constexpr sofa::Size NumberOfNodes = 3; + static constexpr sofa::Size PolynomialOrder = 2; + static constexpr ElementType Element_type = ElementType::QUADRATIC_EDGE; +}; + +struct QuadraticTriangle +{ + static constexpr sofa::Size NumberOfNodes = 6; + static constexpr sofa::Size PolynomialOrder = 2; + static constexpr ElementType Element_type = ElementType::QUADRATIC_TRIANGLE; +}; + +struct QuadraticQuad +{ + static constexpr sofa::Size NumberOfNodes = 9; + static constexpr sofa::Size PolynomialOrder = 2; + static constexpr ElementType Element_type = ElementType::QUADRATIC_QUAD; +}; + +struct QuadraticTetrahedron +{ + static constexpr sofa::Size NumberOfNodes = 10; + static constexpr sofa::Size PolynomialOrder = 2; + static constexpr ElementType Element_type = ElementType::QUADRATIC_TETRAHEDRON; +}; + +struct QuadraticHexahedron +{ + static constexpr sofa::Size NumberOfNodes = 27; + static constexpr sofa::Size PolynomialOrder = 2; + static constexpr ElementType Element_type = ElementType::QUADRATIC_HEXAHEDRON; +}; + +struct QuadraticPrism +{ + static constexpr sofa::Size NumberOfNodes = 18; + static constexpr sofa::Size PolynomialOrder = 2; + static constexpr ElementType Element_type = ElementType::QUADRATIC_PRISM; +}; + +struct QuadraticPyramid +{ + static constexpr sofa::Size NumberOfNodes = 14; + static constexpr sofa::Size PolynomialOrder = 2; + static constexpr ElementType Element_type = ElementType::QUADRATIC_PYRAMID; +}; + +} diff --git a/Sofa/framework/Topology/CMakeLists.txt b/Sofa/framework/Topology/CMakeLists.txt index 9c53e061317..50bbe31b776 100644 --- a/Sofa/framework/Topology/CMakeLists.txt +++ b/Sofa/framework/Topology/CMakeLists.txt @@ -12,6 +12,7 @@ set(HEADER_FILES ${SOFATOPOLOGYSRC_ROOT}/Edge.h ${SOFATOPOLOGYSRC_ROOT}/Triangle.h ${SOFATOPOLOGYSRC_ROOT}/Quad.h + ${SOFATOPOLOGYSRC_ROOT}/QuadraticElements.h ${SOFATOPOLOGYSRC_ROOT}/Tetrahedron.h ${SOFATOPOLOGYSRC_ROOT}/Pentahedron.h ${SOFATOPOLOGYSRC_ROOT}/Prism.h diff --git a/Sofa/framework/Topology/src/sofa/topology/QuadraticElements.h b/Sofa/framework/Topology/src/sofa/topology/QuadraticElements.h new file mode 100644 index 00000000000..4ce6ae75351 --- /dev/null +++ b/Sofa/framework/Topology/src/sofa/topology/QuadraticElements.h @@ -0,0 +1,39 @@ +/****************************************************************************** +* 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 + +namespace sofa::topology +{ + +using QuadraticEdge = sofa::topology::Element; +using QuadraticTriangle = sofa::topology::Element; +using QuadraticQuad = sofa::topology::Element; +using QuadraticTetrahedron = sofa::topology::Element; +using QuadraticHexahedron = sofa::topology::Element; +using QuadraticPrism = sofa::topology::Element; +using QuadraticPyramid = sofa::topology::Element; + +} diff --git a/Sofa/framework/Topology/src/sofa/topology/Topology.h b/Sofa/framework/Topology/src/sofa/topology/Topology.h index 16a54ffcb80..108828dd683 100644 --- a/Sofa/framework/Topology/src/sofa/topology/Topology.h +++ b/Sofa/framework/Topology/src/sofa/topology/Topology.h @@ -30,4 +30,5 @@ #include #include #include +#include #include diff --git a/examples/Component/Engine/Generate/LinearToHigherOrderElements.scn b/examples/Component/Engine/Generate/LinearToHigherOrderElements.scn new file mode 100644 index 00000000000..ddd76503b8b --- /dev/null +++ b/examples/Component/Engine/Generate/LinearToHigherOrderElements.scn @@ -0,0 +1,71 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/parallel/assembled/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/parallel/assembled/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..b1d35919472 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/parallel/assembled/CorotationalFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/parallel/matrixfree/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/parallel/matrixfree/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..2766817da67 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/parallel/matrixfree/CorotationalFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/sequential/assembled/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/sequential/assembled/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..1178c9d0d91 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/sequential/assembled/CorotationalFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/sequential/matrixfree/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/sequential/matrixfree/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..9429cbf762e --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/corotational/sequential/matrixfree/CorotationalFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/parallel/assembled/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/parallel/assembled/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..1b88e2a5b8d --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/parallel/assembled/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/parallel/matrixfree/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/parallel/matrixfree/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..719c77d6f73 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/parallel/matrixfree/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/sequential/assembled/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/sequential/assembled/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..ec1fb2dc36b --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/sequential/assembled/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/sequential/matrixfree/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/sequential/matrixfree/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..3d89220c35c --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_hexahedron/linear/sequential/matrixfree/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/QuadraticTopology.xml b/examples/Validation/cantilever_beam/quadratic_tetrahedron/QuadraticTopology.xml new file mode 100644 index 00000000000..12a5c6cb23e --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/QuadraticTopology.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/parallel/assembled/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/parallel/assembled/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..3991ca5a482 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/parallel/assembled/CorotationalFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/parallel/matrixfree/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/parallel/matrixfree/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..9fa5a384412 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/parallel/matrixfree/CorotationalFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/sequential/assembled/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/sequential/assembled/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..310274c85f3 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/sequential/assembled/CorotationalFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/sequential/matrixfree/CorotationalFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/sequential/matrixfree/CorotationalFEMForceField.scn new file mode 100644 index 00000000000..cba4d4f08ce --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/corotational/sequential/matrixfree/CorotationalFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/parallel/assembled/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/parallel/assembled/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..7b6666c5513 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/parallel/assembled/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/parallel/matrixfree/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/parallel/matrixfree/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..1a15b3d914c --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/parallel/matrixfree/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/sequential/assembled/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/sequential/assembled/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..7b6666c5513 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/sequential/assembled/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/sequential/matrixfree/LinearSmallStrainFEMForceField.scn b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/sequential/matrixfree/LinearSmallStrainFEMForceField.scn new file mode 100644 index 00000000000..baaf8a3fa26 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/linear/sequential/matrixfree/LinearSmallStrainFEMForceField.scn @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Validation/cantilever_beam/quadratic_tetrahedron/plugins.xml b/examples/Validation/cantilever_beam/quadratic_tetrahedron/plugins.xml new file mode 100644 index 00000000000..9de8a5fb9d3 --- /dev/null +++ b/examples/Validation/cantilever_beam/quadratic_tetrahedron/plugins.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + + + + +