diff --git a/CMakeLists.txt b/CMakeLists.txt index 235bf84..867e8ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,7 @@ sofa_find_package(Sofa.Component.Controller REQUIRED) sofa_find_package(Sofa.Component.Topology.Container.Dynamic REQUIRED) sofa_find_package(Sofa.Component.StateContainer REQUIRED) sofa_find_package(Sofa.Component.Mapping.Linear REQUIRED) +sofa_find_package(Sofa.Component.SolidMechanics.FEM.Elastic REQUIRED) sofa_find_package(Sofa.GL REQUIRED) set(README_FILE README.md) @@ -77,6 +78,7 @@ target_link_libraries(${PROJECT_NAME} Sofa.Component.Topology.Container.Dynamic Sofa.Component.StateContainer Sofa.Component.Mapping.Linear + Sofa.Component.SolidMechanics.FEM.Elastic Sofa.GL ) diff --git a/src/Shell/forcefield/TriangularShellForceField.h b/src/Shell/forcefield/TriangularShellForceField.h index 3caafd0..5027c27 100644 --- a/src/Shell/forcefield/TriangularShellForceField.h +++ b/src/Shell/forcefield/TriangularShellForceField.h @@ -30,7 +30,7 @@ #endif -#include +#include #include #include #include @@ -57,7 +57,7 @@ namespace forcefield using namespace sofa::type; using sofa::type::vector; using namespace sofa::core::topology; -using namespace sofa::core::behavior; +using namespace sofa::component::solidmechanics::fem::elastic; /// This class can be overridden if needed for additional storage within template specializations. template @@ -68,14 +68,14 @@ class TriangularShellForceFieldInternalData template -class TriangularShellForceField : public core::behavior::ForceField +class TriangularShellForceField : public BaseLinearElasticityFEMForceField { public: - SOFA_CLASS(SOFA_TEMPLATE(TriangularShellForceField,DataTypes), SOFA_TEMPLATE(core::behavior::ForceField,DataTypes)); + SOFA_CLASS(SOFA_TEMPLATE(TriangularShellForceField,DataTypes), SOFA_TEMPLATE(BaseLinearElasticityFEMForceField,DataTypes)); - typedef core::behavior::ForceField Inherited; - typedef typename DataTypes::VecCoord VecCoord; - typedef typename DataTypes::VecDeriv VecDeriv; + typedef BaseLinearElasticityFEMForceField Inherited; + typedef typename DataTypes::VecCoord VecCoord; + typedef typename DataTypes::VecDeriv VecDeriv; //typedef typename DataTypes::VecReal VecReal; typedef typename DataTypes::Coord Coord; @@ -206,8 +206,9 @@ class TriangularShellForceField : public core::behavior::ForceField sofa::core::topology::BaseMeshTopology* getTopology() {return _topology;} - Data d_poisson; - Data d_young; + using BaseLinearElasticityFEMForceField::getYoungModulusInElement; + using BaseLinearElasticityFEMForceField::getPoissonRatioInElement; + Data d_thickness; Data d_membraneElement; Data d_bendingElement; diff --git a/src/Shell/forcefield/TriangularShellForceField.inl b/src/Shell/forcefield/TriangularShellForceField.inl index fe379f7..cef57ba 100644 --- a/src/Shell/forcefield/TriangularShellForceField.inl +++ b/src/Shell/forcefield/TriangularShellForceField.inl @@ -73,9 +73,7 @@ void TriangularShellForceField::TRQSTriangleHandler::applyCreateFunct // -------------------------------------------------------------------------------------- template TriangularShellForceField::TriangularShellForceField() - : d_poisson(initData(&d_poisson,(Real)0.45,"poissonRatio","Poisson ratio in Hooke's law")) - , d_young(initData(&d_young,(Real)3000.,"youngModulus","Young modulus in Hooke's law")) - , d_thickness(initData(&d_thickness,(Real)0.1,"thickness","Thickness of the plates")) + : d_thickness(initData(&d_thickness,(Real)0.1,"thickness","Thickness of the plates")) , d_membraneElement(initData(&d_membraneElement, "membraneElement", "The membrane element to use")) , d_bendingElement(initData(&d_bendingElement, "bendingElement", "The bending plate element to use")) , d_corotated(initData(&d_corotated, true, "corotated", "Compute forces in corotational frame")) @@ -543,8 +541,8 @@ void TriangularShellForceField::computeRotation(Transformation& R, co template void TriangularShellForceField::computeMaterialStiffness() { - Real E = d_young.getValue(), - nu = d_poisson.getValue(), + Real E = getYoungModulusInElement(0), + nu = getPoissonRatioInElement(0), t = d_thickness.getValue(); materialMatrix[0][0] = 1.0; @@ -1131,25 +1129,25 @@ void TriangularShellForceField::andesTemplate(StiffnessMatrix &K, con template void TriangularShellForceField::computeStiffnessMatrixAll3I(StiffnessMatrix &K, TriangleInformation &tinfo) { - return andesTemplate(K, tinfo, 1.0, AndesBeta(4.0/9.0, 1.0/12.0, 5.0/12.0, 1.0/2.0, 0.0, 1.0/3.0, -1.0/3.0, -1.0/12.0, -1.0/2.0, -5.0/12.0)); + return andesTemplate(K, tinfo, 1.0, AndesBeta{4.0/9.0, 1.0/12.0, 5.0/12.0, 1.0/2.0, 0.0, 1.0/3.0, -1.0/3.0, -1.0/12.0, -1.0/2.0, -5.0/12.0}); } template void TriangularShellForceField::computeStiffnessMatrixAll3M(StiffnessMatrix &K, TriangleInformation &tinfo) { - return andesTemplate(K, tinfo, 1.0, AndesBeta(4.0/9.0, 1.0/4.0, 5.0/4.0, 3.0/2.0, 0.0, 1.0, -1.0, -1.0/4.0, -3.0/2.0, -5.0/4.0)); + return andesTemplate(K, tinfo, 1.0, AndesBeta{4.0/9.0, 1.0/4.0, 5.0/4.0, 3.0/2.0, 0.0, 1.0, -1.0, -1.0/4.0, -3.0/2.0, -5.0/4.0}); } template void TriangularShellForceField::computeStiffnessMatrixAllLS(StiffnessMatrix &K, TriangleInformation &tinfo) { - return andesTemplate(K, tinfo, 1.0, AndesBeta(4.0/9.0, 3.0/20.0, 3.0/4.0, 9.0/10.0, 0.0, 3.0/5.0, -3.0/5.0, -3.0/20.0, -9.0/10.0, -3.0/4.0)); + return andesTemplate(K, tinfo, 1.0, AndesBeta{4.0/9.0, 3.0/20.0, 3.0/4.0, 9.0/10.0, 0.0, 3.0/5.0, -3.0/5.0, -3.0/20.0, -9.0/10.0, -3.0/4.0}); } template void TriangularShellForceField::computeStiffnessMatrixLSTRet(StiffnessMatrix &K, TriangleInformation &tinfo) { - return andesTemplate(K, tinfo, 4.0/3.0, AndesBeta(1.0/2.0, 2.0/3.0, -2.0/3.0, 0.0, 0.0, -4.0/3.0, 4.0/3.0, -2.0/3.0, 0.0, 2.0/3.0)); + return andesTemplate(K, tinfo, 4.0/3.0, AndesBeta{1.0/2.0, 2.0/3.0, -2.0/3.0, 0.0, 0.0, -4.0/3.0, 4.0/3.0, -2.0/3.0, 0.0, 2.0/3.0}); } // Optimal ANDES membrane element @@ -1157,8 +1155,9 @@ void TriangularShellForceField::computeStiffnessMatrixLSTRet(Stiffnes template void TriangularShellForceField::computeStiffnessMatrixAndesOpt(StiffnessMatrix &K, TriangleInformation &tinfo) { - Real beta0 = helper::rmax(0.5 - 2.0*d_poisson.getValue()*d_poisson.getValue(), 0.01); - return andesTemplate(K, tinfo, 3.0/2.0, AndesBeta(beta0, 1.0, 2.0, 1.0, 0.0, 1.0, -1.0, -1.0, -1.0, -2.0)); + Real poissonRatio = getPoissonRatioInElement(0); + Real beta0 = helper::rmax(0.5 - 2.0 * poissonRatio * poissonRatio, 0.01); + return andesTemplate(K, tinfo, 3.0/2.0, AndesBeta{beta0, 1.0, 2.0, 1.0, 0.0, 1.0, -1.0, -1.0, -1.0, -2.0}); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -