diff --git a/CMakeLists.txt b/CMakeLists.txt index c6569a363..2ecb75afd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,10 +58,19 @@ set(HEADER_FILES ${BEAMADAPTER_SRC}/component/mapping/MultiAdaptiveBeamMapping.h ${BEAMADAPTER_SRC}/component/mapping/MultiAdaptiveBeamMapping.inl + ${BEAMADAPTER_SRC}/component/model/BaseRodSectionMaterial.h + ${BEAMADAPTER_SRC}/component/model/BaseRodSectionMaterial.inl + ${BEAMADAPTER_SRC}/component/model/RodMeshSection.h + ${BEAMADAPTER_SRC}/component/model/RodMeshSection.inl + ${BEAMADAPTER_SRC}/component/model/RodSpireSection.h + ${BEAMADAPTER_SRC}/component/model/RodSpireSection.inl + ${BEAMADAPTER_SRC}/component/model/RodStraightSection.h + ${BEAMADAPTER_SRC}/component/model/RodStraightSection.inl + ${BEAMADAPTER_SRC}/utils/BeamSection.h ${BEAMADAPTER_SRC}/utils/BeamActions.h ${BEAMADAPTER_SRC}/utils/deprecatedcomponent.h - ) +) set(SOURCE_FILES ${BEAMADAPTER_SRC}/initBeamAdapter.cpp @@ -86,7 +95,11 @@ set(SOURCE_FILES ${BEAMADAPTER_SRC}/component/mapping/AdaptiveBeamMapping.cpp ${BEAMADAPTER_SRC}/component/mapping/BeamLengthMapping.cpp ${BEAMADAPTER_SRC}/component/mapping/MultiAdaptiveBeamMapping.cpp - ) + + ${BEAMADAPTER_SRC}/component/model/RodMeshSection.cpp + ${BEAMADAPTER_SRC}/component/model/RodSpireSection.cpp + ${BEAMADAPTER_SRC}/component/model/RodStraightSection.cpp +) if(SofaImplicitField_FOUND) set(HEADER_FILES ${HEADER_FILES} diff --git a/doc/AdaptiveBeamMapping_00000001.png b/doc/AdaptiveBeamMapping_00000001.png new file mode 100644 index 000000000..3137425a2 Binary files /dev/null and b/doc/AdaptiveBeamMapping_00000001.png differ diff --git a/doc/BeamAdapter.pdf b/doc/BeamAdapter.pdf index 3ebc90859..43ccf3cb9 100644 Binary files a/doc/BeamAdapter.pdf and b/doc/BeamAdapter.pdf differ diff --git a/examples/3instruments.scn b/examples/3instruments.scn index 2e1843149..8d876cdad 100644 --- a/examples/3instruments.scn +++ b/examples/3instruments.scn @@ -24,25 +24,34 @@ - + + + + + - + + + + + - - + + + + + + diff --git a/examples/3instruments_collis.scn b/examples/3instruments_collis.scn index 4a5511faa..7c30267b6 100644 --- a/examples/3instruments_collis.scn +++ b/examples/3instruments_collis.scn @@ -38,24 +38,33 @@ - + + + + + - + + + + + - + + + + + @@ -83,7 +92,7 @@ - diff --git a/examples/SingleBeam.scn b/examples/SingleBeam.scn index 8899cf388..4f9d73e1c 100644 --- a/examples/SingleBeam.scn +++ b/examples/SingleBeam.scn @@ -6,7 +6,9 @@ + + @@ -15,8 +17,8 @@ - - + + diff --git a/examples/SingleBeamDeployment.scn b/examples/SingleBeamDeployment.scn index ca9739617..6b6287acc 100644 --- a/examples/SingleBeamDeployment.scn +++ b/examples/SingleBeamDeployment.scn @@ -15,11 +15,10 @@ - + + + + diff --git a/examples/SingleBeamDeploymentCollision.scn b/examples/SingleBeamDeploymentCollision.scn index b455d0f44..4797b11e3 100644 --- a/examples/SingleBeamDeploymentCollision.scn +++ b/examples/SingleBeamDeploymentCollision.scn @@ -32,11 +32,11 @@ - + + + + + diff --git a/examples/Tool_from_loader.scn b/examples/Tool_from_loader.scn index 544bcba57..e76e8fd76 100644 --- a/examples/Tool_from_loader.scn +++ b/examples/Tool_from_loader.scn @@ -23,11 +23,12 @@ - - + + + + + + diff --git a/src/BeamAdapter/component/controller/InterventionalRadiologyController.inl b/src/BeamAdapter/component/controller/InterventionalRadiologyController.inl index 1a6026511..8ccbb5ab7 100644 --- a/src/BeamAdapter/component/controller/InterventionalRadiologyController.inl +++ b/src/BeamAdapter/component/controller/InterventionalRadiologyController.inl @@ -248,6 +248,9 @@ void InterventionalRadiologyController::onMouseEvent(MouseEvent * mev template void InterventionalRadiologyController::onKeyPressedEvent(KeypressedEvent *kev) { + + return; + /// Control keys for interventonal Radiology simulations: switch(kev->getKey()) { @@ -308,56 +311,56 @@ template void InterventionalRadiologyController::onBeginAnimationStep(const double dt) { SOFA_UNUSED(dt); - - BaseContext* context = getContext(); - auto xInstrTip = sofa::helper::getWriteOnlyAccessor(d_xTip); - if(m_FF || m_RW) - { - int id = d_controlledInstrument.getValue(); - if (id >= (int)xInstrTip.size()) - { - msg_warning()<<"Controlled Instument num "<getDt(); - else - { - unsigned int newSensorData = m_currentSensorData + 1; - - while( m_sensorMotionData[newSensorData][0] < context->getTime() ) - { - m_currentSensorData = newSensorData; - newSensorData++; - } - if(newSensorData >= m_sensorMotionData.size()) - { - xInstrTip[id] = 0; - } - else - { - xInstrTip[id] += m_sensorMotionData[m_currentSensorData][1]; - } - } - } - if (m_RW) - { - xInstrTip[id] -= d_speed.getValue()* context->getDt(); - // verif min x : - if ( xInstrTip[id] < 0.0) - { - xInstrTip[id] = 0.0; - m_RW = false; - } - } - } - - /// The tip of the instrument can not be further than its total length - for (unsigned int i=0; i m_instrumentsList[i]->getRestTotalLength() ) - xInstrTip[i] = m_instrumentsList[i]->getRestTotalLength(); + + //BaseContext* context = getContext(); + //auto xInstrTip = sofa::helper::getWriteOnlyAccessor(d_xTip); + //if(m_FF || m_RW) + //{ + // int id = d_controlledInstrument.getValue(); + // if (id >= (int)xInstrTip.size()) + // { + // msg_warning()<<"Controlled Instument num "<getDt(); + // else + // { + // unsigned int newSensorData = m_currentSensorData + 1; + + // while( m_sensorMotionData[newSensorData][0] < context->getTime() ) + // { + // m_currentSensorData = newSensorData; + // newSensorData++; + // } + // if(newSensorData >= m_sensorMotionData.size()) + // { + // xInstrTip[id] = 0; + // } + // else + // { + // xInstrTip[id] += m_sensorMotionData[m_currentSensorData][1]; + // } + // } + // } + // if (m_RW) + // { + // xInstrTip[id] -= d_speed.getValue()* context->getDt(); + // // verif min x : + // if ( xInstrTip[id] < 0.0) + // { + // xInstrTip[id] = 0.0; + // m_RW = false; + // } + // } + //} + + ///// The tip of the instrument can not be further than its total length + //for (unsigned int i=0; i m_instrumentsList[i]->getRestTotalLength() ) + // xInstrTip[i] = m_instrumentsList[i]->getRestTotalLength(); applyInterventionalRadiologyController(); } diff --git a/src/BeamAdapter/component/engine/SteerableCatheter.h b/src/BeamAdapter/component/engine/SteerableCatheter.h index 17ceebb98..27819a66a 100644 --- a/src/BeamAdapter/component/engine/SteerableCatheter.h +++ b/src/BeamAdapter/component/engine/SteerableCatheter.h @@ -104,9 +104,10 @@ class SteerableCatheter : public WireRestShape /// Bring inherited attributes and function in the current lookup context. /// otherwise any access to the base::attribute would require /// the "this->" approach. - using Inherit1::d_spireDiameter; - using Inherit1::d_length; - using Inherit1::d_straightLength; + Data d_length; + Data d_straightLength; + Data d_spireDiameter; + Data d_spireHeight; /////////////////////////////////////////////////////////////////////////// }; diff --git a/src/BeamAdapter/component/engine/WireRestShape.cpp b/src/BeamAdapter/component/engine/WireRestShape.cpp index dacba394e..34572378a 100644 --- a/src/BeamAdapter/component/engine/WireRestShape.cpp +++ b/src/BeamAdapter/component/engine/WireRestShape.cpp @@ -52,7 +52,7 @@ using namespace sofa::defaulttype; /// //////////////////////////////////////////////////////////////////////////////////////////////////// -static int WireRestShapeClass = core::RegisterObject("Describe the shape functions on multiple segments using curvilinear abscissa") +const int WireRestShapeClass = core::RegisterObject("Describe the shape functions on multiple segments using curvilinear abscissa") .add< WireRestShape >(true) ; diff --git a/src/BeamAdapter/component/engine/WireRestShape.h b/src/BeamAdapter/component/engine/WireRestShape.h index 8c77de028..2ab349ce9 100644 --- a/src/BeamAdapter/component/engine/WireRestShape.h +++ b/src/BeamAdapter/component/engine/WireRestShape.h @@ -33,11 +33,11 @@ #include -#include +#include + #include #include #include -#include #include namespace sofa::component::engine @@ -48,14 +48,17 @@ namespace _wirerestshape_ using sofa::core::topology::TopologyContainer; using sofa::component::topology::container::dynamic::EdgeSetTopologyModifier; -using sofa::component::topology::mapping::Edge2QuadTopologicalMapping; using sofa::core::loader::MeshLoader; +using namespace sofa::beamadapter; + /** * \class WireRestShape * \brief Describe the shape functions on multiple segments - * - * Describe the shape functions on multiple segments using curvilinear abscissa + * + * Describe the full shape of a Wire with a given set of @sa BaseRodSectionMaterial. The wire is discretized by a set of beams (given by the keyPoints and the relatives Beam density) + * @sa d_keyPoints and @d_density are computed by method @sa initLengths using the set of rod sections description. + * This component compute the beam discretization and the shape functions on multiple segments using curvilinear abscissa. */ template class WireRestShape : public core::objectmodel::BaseObject @@ -69,8 +72,6 @@ class WireRestShape : public core::objectmodel::BaseObject using Vec3 = sofa::type::Vec<3, Real>; using Quat = sofa::type::Quat; - using BeamSection = sofa::beamadapter::BeamSection; - /** * @brief Default Constructor. */ @@ -82,25 +83,19 @@ class WireRestShape : public core::objectmodel::BaseObject virtual ~WireRestShape() = default; /////////////////////////// Inherited from BaseObject ////////////////////////////////////////// - void parse(core::objectmodel::BaseObjectDescription* arg) override; void init() override ; - - void draw(const core::visual::VisualParams * vparams) override ; - - /////////////////////////// Methods of WireRestShape ////////////////////////////////////////// - /// For coils: a part of the coil instrument can be brokenIn2 (by default the point of release is the end of the straight length) - Real getReleaseCurvAbs() const {return d_straightLength.getValue();} + /////////////////////////// Methods of WireRestShape ////////////////////////////////////////// /// This function is called by the force field to evaluate the rest position of each beam void getRestTransformOnX(Transform &global_H_local, const Real &x); /// This function gives the Young modulus and Poisson's coefficient of the beam depending on the beam position - void getYoungModulusAtX(const Real& x_curv, Real& youngModulus, Real& cPoisson); + void getYoungModulusAtX(const Real& x_curv, Real& youngModulus, Real& cPoisson) const; /// This function gives the mass density and the BeamSection data depending on the beam position - void getInterpolationParam(const Real& x_curv, Real &_rho, Real &_A, Real &_Iy , Real &_Iz, Real &_Asy, Real &_Asz, Real &_J); + void getInterpolationParam(const Real& x_curv, Real &_rho, Real &_A, Real &_Iy , Real &_Iz, Real &_Asy, Real &_Asz, Real &_J) const; /** * This function provides a type::vector with the curviliar abscissa of the noticeable point(s) @@ -110,77 +105,48 @@ class WireRestShape : public core::objectmodel::BaseObject /// Functions enabling to load and use a geometry given from OBJ external file - void initRestConfig(); - void getRestPosNonProcedural(Real& abs, Coord &p); void computeOrientation(const Vec3& AB, const Quat& Q, Quat &result); - void initFromLoader(); - bool checkTopology(); - - [[nodiscard]] bool fillTopology(); + + Real getLength() ; - void getCollisionSampling(Real &dx, const Real &x_curv) ; + void getCollisionSampling(Real &dx, const Real &x_curv); void getNumberOfCollisionSegment(Real &dx, unsigned int &numLines) ; - //TODO(dmarchal 2017-05-17) Please specify who and when it will be done either a time after wich - //we can remove the todo. - // todo => topological change ! - void releaseWirePart(); - void rotateFrameForAlignX(const Quat &input, Vec3 &x, Quat &output); + + /////////////////////////// Deprecated Methods ////////////////////////////////////////// + + /// For coils: a part of the coil instrument can be brokenIn2 (by default the point of release is the end of the straight length) + Real getReleaseCurvAbs() const { + msg_warning() << "Releasing catheter or brokenIn2 mode is not anymore supported. Feature has been removed after release v23.06"; + return 0.0; + } + + void releaseWirePart() { + msg_warning() << "Releasing catheter or brokenIn2 mode is not anymore supported. Feature has been removed after release v23.06"; + } + +protected: + /// Internal method to init Lengths vector @sa d_keyPoints using the length of each materials @sa l_sectionMaterials. + void initLengths(); + /// Internal method to init Edge Topology @sa _topology using the list of materials @sa l_sectionMaterials. Returns false if init can't be performed. + bool initTopology(); public: - /// Analitical creation of wire shape... - Data d_isAProceduralShape; - Data d_nonProceduralScale; - Data d_length; - Data d_straightLength; - Data d_spireDiameter; - Data d_spireHeight; Data > d_density; Data > d_keyPoints; - Data< int > d_numEdges; - Data > d_numEdgesCollis; - - /// User Data about the Young modulus - Data d_poissonRatio; - Data d_youngModulus1; - Data d_youngModulus2; - - /// Radius - Data d_radius1; - Data d_radius2; - Data d_innerRadius1; - Data d_innerRadius2; - - Data d_massDensity1; - Data d_massDensity2; - - /// broken in 2 case - Data d_brokenIn2; - Data d_drawRestShape; - -private: - /// Data required for the File loading - type::vector m_localRestPositions; - type::vector m_localRestTransforms; - type::vector m_curvAbs ; - double m_absOfGeometry {0}; - BeamSection beamSection1; - BeamSection beamSection2; + /// Vector or links to the Wire section material. The order of the linked material will define the WireShape structure. + MultiLink, BaseRodSectionMaterial, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_sectionMaterials; +private: /// Link to be set to the topology container in the component graph. SingleLink, TopologyContainer, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology; /// Pointer to the topology container, should be set using @sa l_topology, otherwise will search for one in current Node. TopologyContainer* _topology{ nullptr }; /// Pointer to the topology modifier. Will be set at init by searching one in @sa _topology context. EdgeSetTopologyModifier* edgeMod{ nullptr }; - - /// Link to be set to the topology container in the component graph. - SingleLink, MeshLoader, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_loader; - /// Pointer to the MeshLoader, should be set using @sa l_loader, otherwise will search for one in current Node. - MeshLoader* loader{ nullptr }; }; diff --git a/src/BeamAdapter/component/engine/WireRestShape.inl b/src/BeamAdapter/component/engine/WireRestShape.inl index dda4ad4b5..1ae3484e6 100644 --- a/src/BeamAdapter/component/engine/WireRestShape.inl +++ b/src/BeamAdapter/component/engine/WireRestShape.inl @@ -35,12 +35,11 @@ #include #include -#include #include #include -#define EPSILON 0.0000000001 +#define EPSILON 0.0001 #define VERIF 1 namespace sofa::component::engine @@ -57,83 +56,15 @@ using sofa::core::objectmodel::BaseContext ; * @brief Default Constructor. */ template -WireRestShape::WireRestShape() : - //TODO(dmarchal 2017-05-17) not sure that procedural & nonProceduralScale are very understandable name...are they exclusives ? - //if so have look in my comment in the init section. - d_isAProceduralShape( initData(&d_isAProceduralShape,(bool)true,"isAProceduralShape","is the guidewire shape mathemetically defined ?") ) - , d_nonProceduralScale( initData ( &d_nonProceduralScale, (Real)1.0, "nonProceduralScale", "scale of the model defined by file" ) ) - , d_length(initData(&d_length, (Real)1.0, "length", "total length of the wire instrument")) - , d_straightLength(initData(&d_straightLength, (Real)0.0, "straightLength", "length of the initial straight shape")) - , d_spireDiameter(initData(&d_spireDiameter, (Real)0.1, "spireDiameter", "diameter of the spire")) - , d_spireHeight(initData(&d_spireHeight, (Real)0.01, "spireHeight", "height between each spire")) - , d_density(initData(&d_density, "densityOfBeams", "density of beams between key points")) - , d_keyPoints(initData(&d_keyPoints,"keyPoints","key points of the shape (curv absc)")) - , d_numEdges(initData(&d_numEdges, 10, "numEdges","number of Edges for the visual model")) - , d_numEdgesCollis(initData(&d_numEdgesCollis,"numEdgesCollis", "number of Edges for the collision model" )) - , d_poissonRatio(initData(&d_poissonRatio,(Real)0.49,"poissonRatio","Poisson Ratio")) - , d_youngModulus1(initData(&d_youngModulus1,(Real)5000,"youngModulus","Young Modulus")) - , d_youngModulus2(initData(&d_youngModulus2,(Real)3000,"youngModulusExtremity","youngModulus for beams at the extremity\nonly if not straight")) - , d_radius1(initData(&d_radius1,(Real)1.0f,"radius","radius")) - , d_radius2(initData(&d_radius2,(Real)1.0f,"radiusExtremity","radius for beams at the extremity\nonly if not straight")) - , d_innerRadius1(initData(&d_innerRadius1,(Real)0.0f,"innerRadius","inner radius if it applies")) - , d_innerRadius2(initData(&d_innerRadius2,(Real)0.0f,"innerRadiusExtremity","inner radius for beams at the extremity\nonly if not straight")) - , d_massDensity1(initData(&d_massDensity1,(Real)1.0,"massDensity", "Density of the mass (usually in kg/m^3)" )) - , d_massDensity2(initData(&d_massDensity2,(Real)1.0,"massDensityExtremity", "Density of the mass at the extremity\nonly if not straight" )) - , d_brokenIn2(initData(&d_brokenIn2, (bool)false, "brokenIn2", "")) - , d_drawRestShape(initData(&d_drawRestShape, (bool)false, "draw", "draw rest shape")) - , l_topology(initLink("topology", "link to the topology container")) - , l_loader(initLink("loader", "link to the MeshLoader")) +WireRestShape::WireRestShape() + : d_density(initData(&d_density, "densityOfBeams", "density of beams between key points")) + , d_keyPoints(initData(&d_keyPoints,"keyPoints","key points of the shape (curv absc)")) + , l_sectionMaterials(initLink("wireMaterials", "link to Wire Section Materials (to be ordered according to the instrument, from handle to tip)")) + , l_topology(initLink("topology", "link to the topology container")) { - d_spireDiameter.setGroup("Procedural"); - d_spireHeight.setGroup("Procedural"); -} - -template -void WireRestShape::rotateFrameForAlignX(const Quat &input, Vec3 &x, Quat &output) -{ - x.normalize(); - Vec3 x0=input.inverseRotate(x); - - Real cTheta=x0[0]; - Real theta; - if (cTheta>(1-EPSILON)) - { - output = input; - } - else - { - theta=acos(cTheta); - // axis of rotation - Vec3 dw(0,-x0[2],x0[1]); - dw.normalize(); - - // computation of the rotation - Quat inputRoutput; - inputRoutput.axisToQuat(dw, theta); - output=input*inputRoutput; - } } -template -void WireRestShape::parse(core::objectmodel::BaseObjectDescription* args) -{ - const char* arg = args->getAttribute("procedural") ; - if(arg) - { - msg_warning() << "The attribute 'procedural' has been renamed into 'isAProceduralShape'. " << msgendl - << "To remove this warning you need to update your scene and replace 'procedural' with 'isAProceduralShape'" ; - - /// As arg is owned by the "procedural" attribute it cannot be removed before - /// being copied in the "isAProceduralShape". So please keep the ordering of the - /// two following functions. - args->setAttribute("isAProceduralShape", arg) ; - args->removeAttribute("procedural") ; - - } - - Inherit1::parse(args) ; -} template void WireRestShape::init() @@ -161,146 +92,95 @@ void WireRestShape::init() return; } - if (!d_isAProceduralShape.getValue()) - { - // Get meshLoader, check first if loader has been set using link. Otherwise will search in current context. - loader = l_loader.get(); - - if (!loader) - this->getContext()->get(loader); - if (!loader) { - msg_error() << "Cannot find a mesh loader. Please insert a MeshObjLoader in the same node or use l_loader to specify the path in the scene graph."; - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; - } - else - { - msg_info() << "Found a mesh with " << loader->d_edges.getValue().size() << " edges"; - initFromLoader(); - } - } - else + if (l_sectionMaterials.empty()) { - if (!fillTopology()) - { - msg_error() << "Error while trying to fill the associated topology, setting the state to Invalid"; - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; - } + msg_error() << "No BaseRodSectionMaterial set. At least one material should be set and link using wireMaterials."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; } + + //////////////////////////////////////////////////////// + ////////// keyPoint list and Density Assignement /////// + //////////////////////////////////////////////////////// + + initLengths(); + // Get pointer to the topology Modifier (for topological changes) _topology->getContext()->get(edgeMod); - if (edgeMod == nullptr) { msg_warning() << "No EdgeSetTopologyModifier found in the same node as the topology container: " << _topology->getName() << ". This wire won't support topological changes."; } + initTopology(); + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + msg_info() << "WireRestShape end init"; +} - //////////////////////////////////////////////////////// - ////////// keyPoint list and Density Assignement /////// - //////////////////////////////////////////////////////// - auto keyPointList = sofa::helper::getWriteOnlyAccessor(d_keyPoints); - if(!keyPointList.size()) - { - keyPointList.push_back(0.0); - if(d_straightLength.getValue()>= 0.001*this->d_length.getValue() && d_straightLength.getValue() <= 0.999*d_length.getValue()) - keyPointList.push_back(d_straightLength.getValue()); - keyPointList.push_back(d_length.getValue()); - } - if( d_density.getValue().size() != keyPointList.size()-1) - { - auto densityList = sofa::helper::getWriteOnlyAccessor(d_density); +template +void WireRestShape::initLengths() +{ + auto keyPointList = sofa::helper::getWriteOnlyAccessor(d_keyPoints); + auto densityList = sofa::helper::getWriteOnlyAccessor(d_density); - if(densityList.size() > keyPointList.size()-1 ) - densityList.resize(keyPointList.size()-1); - else - { - densityList.clear(); - - if(d_straightLength.getValue()>= 0.001*this->d_length.getValue() ) - { - int numNodes = (int) floor(5.0*d_straightLength.getValue() / d_length.getValue() ); - densityList.push_back(numNodes); - } - if( d_straightLength.getValue() <= 0.999*d_length.getValue()) - { - int numNodes = (int) floor(20.0*(1.0 - d_straightLength.getValue() / d_length.getValue()) ); - densityList.push_back(numNodes); - } - } - } + keyPointList.resize(l_sectionMaterials.size() + 1); + keyPointList[0] = Real(0.0); - if(!d_numEdgesCollis.getValue().size()) + densityList.resize(l_sectionMaterials.size()); + + for (unsigned int i = 0; i < l_sectionMaterials.size(); ++i) { - auto densityCol = sofa::helper::getWriteOnlyAccessor(d_numEdgesCollis); - densityCol.resize(keyPointList.size()-1); - for (unsigned int i=0; igetLength(); + densityList[i] = rodSection->getNbCollisionEdges(); } - - msg_info() <<"WireRestShape end init" ; - - // Prepare beam sections - double r = this->d_radius1.getValue(); - double rInner = this->d_innerRadius1.getValue(); - this->beamSection1._r = r; - this->beamSection1._rInner = rInner; - this->beamSection1._Iz = M_PI*(r*r*r*r - rInner*rInner*rInner*rInner)/4.0; - this->beamSection1._Iy = this->beamSection1._Iz ; - this->beamSection1._J = this->beamSection1._Iz + this->beamSection1._Iy; - this->beamSection1._A = M_PI*(r*r - rInner*rInner); - this->beamSection1._Asy = 0.0; - this->beamSection1._Asz = 0.0; - - r = this->d_radius2.getValue(); - rInner = this->d_innerRadius2.getValue(); - this->beamSection2._r = r; - this->beamSection2._rInner = rInner; - this->beamSection2._Iz = M_PI*(r*r*r*r - rInner*rInner*rInner*rInner)/4.0; - this->beamSection2._Iy = this->beamSection2._Iz ; - this->beamSection2._J = this->beamSection2._Iz + this->beamSection2._Iy; - this->beamSection2._A = M_PI*(r*r - rInner*rInner); - this->beamSection2._Asy = 0.0; - this->beamSection2._Asz = 0.0; - - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); } template -void WireRestShape::releaseWirePart(){ - - d_brokenIn2.setValue(true); +bool WireRestShape::initTopology() +{ + /// fill topology : + _topology->clear(); + _topology->cleanup(); - if ( edgeMod == nullptr ) + const type::vector& keyPts = d_keyPoints.getValue(); + if (l_sectionMaterials.size() != keyPts.size() - 1) { - msg_error() << "no edgeSetModifier in the node -> cannot do the topological change"; - return; + msg_error() << "Wrong number of inputs. Component can't be init. Number of input materials: " << l_sectionMaterials.size() << ", should be equal to keyPointList.size()-1. keyPointList.size() is equal to: " << keyPts.size(); + return false; } - ///////// remove the edge that is cut ////// - for ( sofa::Size i=0; i<_topology->getNbPoints(); i++) - { - if( _topology->getPX(i) > this->getReleaseCurvAbs() + EPSILON ) - { - type::vector edge_remove; - edge_remove.push_back( i-1 ); - - msg_info() << "releaseWirePart() -> remove edge number "<< i ; - - edgeMod->removeEdges(edge_remove,false); // remove the single edge and do not remove any point... + + Real prev_length = 0.0; + int prev_edges = 0; + int startPtId = 0; + for (auto i = 0; i < l_sectionMaterials.size(); ++i) + { + // Add topology of the material + int nbrVisuEdges = l_sectionMaterials.get(i)->getNbVisualEdges(); + Real length = fabs(keyPts[i + 1] - keyPts[i]); + Real dx = length / nbrVisuEdges; + + // add points from the material + for (int i = startPtId; i < nbrVisuEdges + 1; i++) { + _topology->addPoint(prev_length + i * dx, 0, 0); + } - msg_info() << "WireRestShape _topology name="<<_topology->getName()<<" - numEdges ="<<_topology->getNbEdges() ; - - return; + // add segments from the material + for (int i = prev_edges; i < prev_edges + nbrVisuEdges; i++) { + _topology->addEdge(i, i + 1); } + + prev_length = length; + prev_edges = nbrVisuEdges; + startPtId = 1; // Assume the last point of mat[n] == first point of mat[n+1] } - dmsg_info() <<" Wire Part is brokenIn2... should implement a topo change !" ; + return true; } @@ -308,64 +188,55 @@ template void WireRestShape::getSamplingParameters(type::vector& xP_noticeable, type::vector& nbP_density) const { - - xP_noticeable.clear(); - nbP_density.clear(); - - if (d_brokenIn2.getValue()) - { - for (unsigned int i=0; i getReleaseCurvAbs() ) - break; - xP_noticeable.push_back(x); - nbP_density.push_back(d_density.getValue()[i]); - } - xP_noticeable.push_back( getReleaseCurvAbs()); - - dmsg_info() <<"getSamplingParameters brokenIn2 detected - return xP_noticeable ="<& keyPts = d_keyPoints.getValue(); + + // verify that size of number of materials == size of keyPoints-1 + if (l_sectionMaterials.size() != keyPts.size() - 1) + { + msg_error() << "Problem size of number of materials: " << l_sectionMaterials.size() + << " != size of keyPoints-1 " << keyPts.size()-1 + << ". Returning default values."; + numLines = 20; + dx = totalLength / numLines; return; } - - - for (unsigned int i=1; id_keyPoints.getValue().size(); i++) + + // Check in which section x_used belongs to and get access to this section material + for (auto i = 1; i< keyPts.size(); ++i) { - if( x_used < this->d_keyPoints.getValue()[i] ) + if (x_used <= keyPts[i]) { - numLines = (unsigned int)d_numEdgesCollis.getValue()[i-1]; - dx=(this->d_keyPoints.getValue()[i] - this->d_keyPoints.getValue()[i-1])/numLines; + numLines = l_sectionMaterials.get(i-1)->getNbCollisionEdges(); + + Real length = fabs(keyPts[i] - keyPts[i-1]); + dx = length / numLines; return; } } - dx=d_length.getValue()/20; - msg_error() << " problem is getCollisionSampling : x_curv "<::getRestTransformOnX(Transform &global_H_local, co { Real x_used = x - EPSILON; - if(x_used>d_length.getValue()) - x_used=d_length.getValue(); - - if(x_used<0.0) - x_used=0.0; - - if( x_used < d_straightLength.getValue()) - { - global_H_local.set(Vec3(x_used, 0.0, 0.0 ), Quat()); - return; - } - - if(d_isAProceduralShape.getValue()) - { - Real projetedLength = d_spireDiameter.getValue()*M_PI; - Real lengthSpire=sqrt(d_spireHeight.getValue()*d_spireHeight.getValue() + projetedLength*projetedLength ); - // angle in the z direction - Real phi= atan(d_spireHeight.getValue()/projetedLength); - - Quat Qphi; - Qphi.axisToQuat(Vec3(0,0,1),phi); - - // spire angle (if theta=2*PI, there is a complete spire between startx and x_used) - Real lengthCurve= x_used-d_straightLength.getValue(); - Real numSpire=lengthCurve/lengthSpire; - Real theta= 2*M_PI*numSpire; - - // computation of the Quat - Quat Qtheta; - Qtheta.axisToQuat(Vec3(0,1,0),theta); - Quat newSpireQuat = Qtheta*Qphi; - - - // computation of the position - Real radius=d_spireDiameter.getValue()/2.0; - Vec3 PosEndCurve(radius*sin(theta), numSpire*d_spireHeight.getValue(), radius*(cos(theta)-1) ); - Vec3 SpirePos=PosEndCurve + Vec3(d_straightLength.getValue(),0,0); - - global_H_local.set(SpirePos,newSpireQuat); - } - else - { - x_used = x_used - d_straightLength.getValue(); - x_used = x_used/(d_length.getValue()-d_straightLength.getValue()) * m_absOfGeometry; - - Coord p; - this->getRestPosNonProcedural(x_used,p); - Vec3 PosEndCurve = p.getCenter(); - Quat ExtremityQuat = p.getOrientation(); - Vec3 ExtremityPos = PosEndCurve + Vec3(d_straightLength.getValue(),0,0); - - global_H_local.set(ExtremityPos,ExtremityQuat); - } -} - - -template -void WireRestShape::getYoungModulusAtX(const Real& x_curv, Real& youngModulus, Real& cPoisson) -{ - //Initialization - Real _E1, _E2; - youngModulus = 0.0; - cPoisson = 0.0; - - //Get the two possible values of the Young modulus - _E1 = this->d_youngModulus1.getValue(); - _E2 = this->d_youngModulus2.getValue(); - - //Get User data - cPoisson = this->d_poissonRatio.getValue(); - - //Depending on the position of the beam, determine the Young modulus - if(x_curv <= this->d_straightLength.getValue()) - { - youngModulus = _E1; - } - else - { - if(_E2 == 0.0) - youngModulus = _E1; - else - youngModulus = _E2; + const Real totalLength = this->getLength(); + if (x_used > totalLength) { + x_used = totalLength; } - return; -} - - -template -void WireRestShape::getInterpolationParam(const Real& x_curv, Real &_rho, Real &_A, Real &_Iy , Real &_Iz, Real &_Asy, Real &_Asz, Real &_J) -{ - if(x_curv <= this->d_straightLength.getValue()) - { - if(d_massDensity1.isSet()) - _rho = d_massDensity1.getValue(); - - if(d_radius1.isSet()) - { - _A =beamSection1._A; - _Iy =beamSection1._Iy; - _Iz =beamSection1._Iz; - _Asy =beamSection1._Asy; - _Asz =beamSection1._Asz; - _J =beamSection1._J; - } + else if (x_used < 0.0) { + x_used = 0.0; } - else + + const type::vector& keyPts = d_keyPoints.getValue(); + for (auto i = 1; i < keyPts.size(); ++i) { - if(d_massDensity2.isSet()) - _rho = d_massDensity2.getValue(); - else if(d_massDensity1.isSet()) - _rho = d_massDensity1.getValue(); - - if(d_radius2.isSet()) - { - _A =beamSection2._A; - _Iy =beamSection2._Iy; - _Iz =beamSection2._Iz; - _Asy =beamSection2._Asy; - _Asz =beamSection2._Asz; - _J =beamSection2._J; - } - else if(d_radius1.isSet()) + if (x_used <= keyPts[i]) { - _A =beamSection1._A; - _Iy =beamSection1._Iy; - _Iz =beamSection1._Iz; - _Asy =beamSection1._Asy; - _Asz =beamSection1._Asz; - _J =beamSection1._J; + return l_sectionMaterials.get(i - 1)->getRestTransformOnX(global_H_local, x_used, keyPts[i - 1]); } } -} - -template -bool WireRestShape::checkTopology() -{ - if (!loader->d_edges.getValue().size()) - { - msg_error() << "There is no edges in the topology loaded by " << loader->getName() ; - return false; - } - if (loader->d_triangles.getValue().size()) - { - msg_error() << "There are triangles in the topology loaded by " << loader->getName() ; - return false; - } - - if (loader->d_quads.getValue().size()) - { - msg_error() << "There are quads in the topology loaded by " << loader->getName() ; - return false; - } - - if (loader->d_polygons.getValue().size()) - { - msg_error() << "There are polygons in the topology loaded by " << loader->getName() ; - return false; - } - - //TODO(dmarchal 2017-05-17) when writing a TODO please specify: - // who will do that - // when it will be done - /// \todo check if the topology is like a wire - - - return true; + msg_error() << " problem in getRestTransformOnX : x_curv " << x_used << " is not between keyPoints" << d_keyPoints.getValue(); } template -bool WireRestShape::fillTopology() +void WireRestShape::getYoungModulusAtX(const Real& x_curv, Real& youngModulus, Real& cPoisson) const { - if (!_topology) - { - msg_error() << "Topology is null"; - return false; - } - - const auto length = this->d_length.getValue(); - if (length <= Real(0.0)) - { - msg_error() << "Length is 0 (or negative), check if d_length has been given or computed."; - return false; - } + const Real x_used = x_curv - Real(EPSILON); + const type::vector& keyPts = d_keyPoints.getValue(); - int nbrEdges = d_numEdges.getValue(); - if (nbrEdges <= 0) + // Depending on the position of the beam, determine the corresponding section material and returning its Young modulus + for (auto i = 1; i < keyPts.size(); ++i) { - msg_warning() << "Number of edges has been set to an invalid value: " << nbrEdges << ". Value should be a positive integer. Setting to default value: 10"; - nbrEdges = 10; + if (x_used <= keyPts[i]) + { + return l_sectionMaterials.get(i - 1)->getYoungModulusAtX(youngModulus, cPoisson); + } } - /// fill topology : - _topology->clear(); - _topology->cleanup(); - - Real dx = this->d_length.getValue() / nbrEdges; - - /// add points - for (int i = 0; i < d_numEdges.getValue() + 1; i++) - _topology->addPoint(i * dx, 0, 0); - - /// add segments - for (int i = 0; i < d_numEdges.getValue(); i++) - _topology->addEdge(i, i + 1); - - return true; + msg_error() << " problem in getYoungModulusAtX : x_curv " << x_curv << " is not between keyPoints" << keyPts; } - template -void WireRestShape::initFromLoader() +void WireRestShape::getInterpolationParam(const Real& x_curv, Real &_rho, Real &_A, Real &_Iy , Real &_Iz, Real &_Asy, Real &_Asz, Real &_J) const { - if (!checkTopology()) - { - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; - } - - // this is... dirty but d_straightLength itself is not relevant for loader-created wire - // but the code uses it actively and expects it to not be zero. - if (d_straightLength.getValue() <= 0.0) - { - msg_warning() << "straightLength cannot be 0 (or negative...). Setting a minimum value."; - d_straightLength.setValue(0.0001); - } - - type::vector vertices; - sofa::core::topology::BaseMeshTopology::SeqEdges edges; - - //get the topology position - auto topoVertices = sofa::helper::getReadAccessor(loader->d_positions); - - //copy the topology edges in a local vector - auto topoEdges = sofa::helper::getReadAccessor(loader->d_edges); - edges = topoEdges.ref(); - - /** renumber the vertices **/ - type::vector verticesConnexion; //gives the number of edges connected to a vertex - for(unsigned int i =0; i < topoVertices.size(); i++) - verticesConnexion.push_back(2); - - for(const auto& ed : edges) - { - verticesConnexion[ed[0]]--; - verticesConnexion[ed[1]]--; - } + const Real x_used = x_curv - Real(EPSILON); + const type::vector& keyPts = d_keyPoints.getValue(); - msg_info() << "Successfully compute the vertex connexion" ; - - // check for the first corner of the edge - unsigned int firstIndex = 0; - bool found = false; - while((firstIndex < verticesConnexion.size()) && !found) - { - if(verticesConnexion[firstIndex] == 1) - found = true; - else - firstIndex++; - } - - if(firstIndex == verticesConnexion.size()) - { - msg_error() << "The first vertex of the beam structure is not found, probably because of a closed structure" ; - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; - } - - vertices.push_back(topoVertices[firstIndex]); - - while(edges.size() > 0) + // Check in which section x_used belongs to and get access to this section material + for (auto i = 1; i < keyPts.size(); ++i) { - auto it = edges.begin(); - auto end = edges.end(); - - bool notFound = true; - while (notFound && (it != end)) + if (x_used <= keyPts[i]) { - const auto& ed = (*it); - auto toDel = it; - it++; - if(ed[0] == firstIndex) - { - vertices.push_back(topoVertices[ed[1]]); - firstIndex = ed[1]; - edges.erase(toDel); - notFound = false; - - } - else if(ed[1] == firstIndex) - { - vertices.push_back(topoVertices[ed[0]]); - firstIndex = ed[0]; - edges.erase(toDel); - notFound = false; - } + return l_sectionMaterials.get(i - 1)->getInterpolationParam(_rho, _A, _Iy, _Iz, _Asy, _Asz, _J); } } - msg_info() << "Successfully computed the topology" ; - - m_localRestPositions = vertices; - - for(unsigned int i = 0; i < m_localRestPositions.size() - 1; i++) - m_localRestPositions[i] *= d_nonProceduralScale.getValue(); - - initRestConfig(); - // TODO epernod 2022-08-05: Init from loader seems quite buggy, need to check if this is still needed and working - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); -} - - -template -void WireRestShape::initRestConfig() -{ - m_curvAbs.clear(); - double tot = 0; - m_curvAbs.push_back(0); - Quat input, output; - input.identity(); - m_localRestTransforms.resize(m_localRestPositions.size()); - m_localRestTransforms[0].setOrigin(Vec3(0,0,0)); - m_localRestTransforms[0].setOrientation(input); - - for(unsigned int i = 0; i < m_localRestPositions.size() - 1; i++) - { - Vec3 vec = m_localRestPositions[i+1] - m_localRestPositions[i]; - double norm = vec.norm(); - tot += norm; - - this->rotateFrameForAlignX(input, vec, output); - - input = output; - - m_localRestTransforms[i+1].setOrientation(output); - - Vec3 localPos = m_localRestPositions[i+1] - m_localRestPositions[0]; - - m_localRestTransforms[i+1].setOrigin(localPos); - - m_curvAbs.push_back(tot); - } - m_absOfGeometry = tot; - - Real newLength = d_straightLength.getValue() + m_absOfGeometry; - d_length.setValue(newLength); - - msg_info() <<"Length of the loaded shape = "<< m_absOfGeometry << ", total length with straight length = " << newLength ; - - if (!fillTopology()) - { - msg_error() << "Error while trying to fill the associated topology, setting the state to Invalid"; - this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); - return; - } + msg_error() << " problem in getInterpolationParam : x_curv " << x_curv << " is not between keyPoints" << keyPts; } -template -void WireRestShape::getRestPosNonProcedural(Real& abs, Coord &p) -{ - /*** find the range which includes the "requested" abs ***/ - double startingAbs = 0; unsigned int index = 0; - - while ((startingAbs < abs) && (index < m_localRestPositions.size())) - { - index++; - startingAbs = m_curvAbs[index]; - } - - /*** OOB ***/ - if(abs > startingAbs) - { - msg_error() << "abs = "< typename WireRestShape::Real WireRestShape::getLength() { - if(d_brokenIn2.getValue()) - return d_straightLength.getValue(); - else - return d_length.getValue(); + return d_keyPoints.getValue().back(); } + template void WireRestShape::getNumberOfCollisionSegment(Real &dx, unsigned int &numLines) { numLines = 0; - for (unsigned i=0; igetNbCollisionEdges(); } - dx=d_length.getValue()/numLines; + dx = getLength() / numLines; } + template void WireRestShape::computeOrientation(const Vec3& AB, const Quat& Q, Quat &result) { @@ -824,26 +355,6 @@ void WireRestShape::computeOrientation(const Vec3& AB, const Quat& Q, } -template -void WireRestShape::draw(const core::visual::VisualParams* vparams) -{ - if (!d_drawRestShape.getValue()) - return; - - vparams->drawTool()->saveLastState(); - vparams->drawTool()->setLightingEnabled(false); - - std::vector< sofa::type::Vec3 > points; - points.reserve(m_localRestPositions.size()); - - for (unsigned int i = 0; i < m_localRestPositions.size(); i++) - { - points.emplace_back(m_localRestPositions[i][0], m_localRestPositions[i][1], m_localRestPositions[i][2]); - } - - vparams->drawTool()->drawPoints(points, 10, sofa::type::RGBAColor(1, 0.5, 0.5, 1)); - vparams->drawTool()->restoreLastState(); -} } // namespace _wirerestshape_ using _wirerestshape_::WireRestShape; diff --git a/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.h b/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.h index ac681e6b8..d2a0b910f 100644 --- a/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.h +++ b/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.h @@ -154,6 +154,8 @@ class AdaptiveBeamMapping : public Mapping Data> d_segmentsCurvAbs; /*!< (output) the abscissa of each created point on the collision model */ Data d_parallelMapping; /*!< flag to enable parallel internal computation of apply/applyJ */ + Data d_onlyVisual; + SingleLink, BInterpolation, BaseLink::FLAG_STOREPATH|BaseLink::FLAG_STRONGLINK> l_adaptativebeamInterpolation; diff --git a/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.inl b/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.inl index 8a9029c12..f4df8bdd9 100644 --- a/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.inl +++ b/src/BeamAdapter/component/mapping/AdaptiveBeamMapping.inl @@ -77,6 +77,7 @@ AdaptiveBeamMapping::AdaptiveBeamMapping(State< In >* from, State< Out , d_nbPointsPerBeam(initData(&d_nbPointsPerBeam, 0.0, "nbPointsPerBeam", "if non zero, we will adapt the points depending on the discretization, with this num of points per beam (compatible with useCurvAbs)")) , d_segmentsCurvAbs(initData(&d_segmentsCurvAbs, "segmentsCurvAbs", "the abscissa of each point on the collision model", true, true)) , d_parallelMapping(initData(&d_parallelMapping, false, "parallelMapping", "flag to enable parallel internal computation")) + , d_onlyVisual(initData(&d_onlyVisual, (bool)false, "onlyVisual", "Really not mechanical mapping")) , l_adaptativebeamInterpolation(initLink("interpolation", "Path to the Interpolation component on scene"), interpolation) , m_inputMapping(nullptr) , m_isSubMapping(isSubMapping) @@ -292,6 +293,9 @@ void AdaptiveBeamMapping< TIn, TOut>::apply(const MechanicalParams* mparams, Dat template void AdaptiveBeamMapping< TIn, TOut>::applyJ(const core::MechanicalParams* mparams, Data& dOut, const Data& dIn) { + if (d_onlyVisual.getValue()) + return; + SOFA_UNUSED(mparams); ScopedAdvancedTimer timer("AdaptiveBeamMapping_applyJ"); @@ -377,6 +381,9 @@ void AdaptiveBeamMapping< TIn, TOut>::applyJ(const core::MechanicalParams* mpara template void AdaptiveBeamMapping< TIn, TOut>::applyJT(const core::MechanicalParams* mparams, Data& dOut, const Data& dIn) { + if (d_onlyVisual.getValue()) + return; + SOFA_UNUSED(mparams); ScopedAdvancedTimer timer("AdaptiveBeamMapping_applyJT"); diff --git a/src/BeamAdapter/component/model/BaseRodSectionMaterial.h b/src/BeamAdapter/component/model/BaseRodSectionMaterial.h new file mode 100644 index 000000000..d023b6a69 --- /dev/null +++ b/src/BeamAdapter/component/model/BaseRodSectionMaterial.h @@ -0,0 +1,120 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace sofa::beamadapter +{ + +using sofa::core::loader::MeshLoader; + +/** + * \class BaseRodSectionMaterial + * \brief Base class describing a Rod section which will define a set of beam elements. + * + * This class provide an api to define a rod/wire section using physical and geometry parameters. + * The section will then be modellized by a set of beam elements. Inheriting class should provide the geometry structure: + * @sa RodMeshSection to define a rod using a mesh file, @sa RodSpireSection or @sa RodStraightSection to define procedural shapes. + * Method @sa initSection and @sa getRestTransformOnX should be overriden to provide the correct creation and interpolation. + * + * The rod section is described by: + * - Topology parameters: vertices and edges @sa d_nbEdgesVisu and @sa d_nbEdgesCollis + * - Geometry parameters: radius @sa d_radius, @sa d_innerRadius and length @sa d_length + * - Mechanical parameters: @sa d_poissonRatio and @sa d_youngModulus + */ +template +class BaseRodSectionMaterial : public core::objectmodel::BaseObject +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(BaseRodSectionMaterial, DataTypes), core::objectmodel::BaseObject); + + using Coord = typename DataTypes::Coord; + using Real = typename Coord::value_type; + using Transform = typename sofa::defaulttype::SolidTypes::Transform; + using Vec3 = sofa::type::Vec<3, Real>; + using Quat = sofa::type::Quat; + using Size = sofa::Size; + + /////////////////////////// Inherited from BaseObject ////////////////////////////////////////// + + /// Default Constructor + BaseRodSectionMaterial(); + + /// init method from BaseObject API. Will call internal @see initSection to be overriden by children + void init() override; + + + /////////////////////////// Geometry and physics Getter ////////////////////////////////////////// + + /// Returns the number of visual edges of this section. To be set or computed by child. + [[nodiscard]] int getNbVisualEdges() const { return d_nbEdgesVisu.getValue(); } + + /// Returns the number of collision edges of this section. To be set or computed by child. + [[nodiscard]] int getNbCollisionEdges() const { return d_nbEdgesCollis.getValue(); } + + /// Returns the total length of this section. To be set or computed by child. + [[nodiscard]] Real getLength() const { return d_length.getValue(); } + + /// Returns the Young modulus and Poisson's coefficient of this section + void getYoungModulusAtX(Real& youngModulus, Real& cPoisson) const; + + /// Returns the mass density and the BeamSection of this section + void getInterpolationParam(Real& _rho, Real& _A, Real& _Iy, Real& _Iz, Real& _Asy, Real& _Asz, Real& _J) const; + + + + /// This function is called to get the rest position of the beam depending on the current curved abscisse given in parameter + virtual void getRestTransformOnX(Transform& global_H_local, const Real& x_used, const Real& x_start) + { + SOFA_UNUSED(global_H_local); + SOFA_UNUSED(x_used); + SOFA_UNUSED(x_start); + } + +protected: + /// Internal method to init the section. to be overidden by child. + virtual bool initSection() { return false; } + +public: + Data d_poissonRatio; ///< Data defining the mehcanical Poisson ratio of this section + Data d_youngModulus; ///< Data defining the mehcanical Young Modulus of this section + Data d_massDensity; ///< Data defining the mehcanical mass density of this section + + Data d_radius; ///< Data defining the geometry radius of this section + Data d_innerRadius; ///< Data defining the geometry internal radius of this section is hollow + Data d_length; ///< Data defining the geometry length of this section + + Data d_nbEdgesVisu; ///< Data defining the number of visual edges composing this section + Data d_nbEdgesCollis; ///< Data defining the number of collision edges composing this section + +private: + /// Internal structure to store physical parameter of the a beam section + BeamSection beamSection; +}; + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/BaseRodSectionMaterial.inl b/src/BeamAdapter/component/model/BaseRodSectionMaterial.inl new file mode 100644 index 000000000..ffef90e69 --- /dev/null +++ b/src/BeamAdapter/component/model/BaseRodSectionMaterial.inl @@ -0,0 +1,96 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +namespace sofa::beamadapter +{ + +template +BaseRodSectionMaterial::BaseRodSectionMaterial() + : d_poissonRatio(initData(&d_poissonRatio, (Real)0.49, "poissonRatio", "Poisson Ratio of this section")) + , d_youngModulus(initData(&d_youngModulus, (Real)5000, "youngModulus", "Young Modulus of this section")) + , d_massDensity(initData(&d_massDensity, (Real)1.0, "massDensity", "Density of the mass (usually in kg/m^3)")) + , d_radius(initData(&d_radius, (Real)1.0, "radius", "Full radius of this section")) + , d_innerRadius(initData(&d_innerRadius, (Real)0.0, "innerRadius", "Inner radius of this section if hollow")) + , d_length(initData(&d_length, (Real)1.0, "length", "Total length of this section")) + , d_nbEdgesVisu(initData(&d_nbEdgesVisu, (Size)10, "nbEdgesVisu", "number of Edges for the visual model")) + , d_nbEdgesCollis(initData(&d_nbEdgesCollis, (Size)20, "nbEdgesCollis", "number of Edges for the collision model")) +{ + +} + + +template +void BaseRodSectionMaterial::init() +{ + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Loading); + + // Prepare beam sections + double r = this->d_radius.getValue(); + double rInner = this->d_innerRadius.getValue(); + this->beamSection._r = r; + this->beamSection._rInner = rInner; + this->beamSection._Iz = M_PI * (r * r * r * r - rInner * rInner * rInner * rInner) / 4.0; + this->beamSection._Iy = this->beamSection._Iz; + this->beamSection._J = this->beamSection._Iz + this->beamSection._Iy; + this->beamSection._A = M_PI * (r * r - rInner * rInner); + this->beamSection._Asy = 0.0; + this->beamSection._Asz = 0.0; + + // call delegate method to init the section + bool res = initSection(); + + if (res) + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + else + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); +} + + +template +void BaseRodSectionMaterial::getYoungModulusAtX(Real& youngModulus, Real& cPoisson) const +{ + youngModulus = this->d_youngModulus.getValue(); + cPoisson = this->d_poissonRatio.getValue(); +} + + +template +void BaseRodSectionMaterial::getInterpolationParam(Real& _rho, Real& _A, Real& _Iy, Real& _Iz, Real& _Asy, Real& _Asz, Real& _J) const +{ + if (d_massDensity.isSet()) + _rho = d_massDensity.getValue(); + + if (d_radius.isSet()) + { + _A = beamSection._A; + _Iy = beamSection._Iy; + _Iz = beamSection._Iz; + _Asy = beamSection._Asy; + _Asz = beamSection._Asz; + _J = beamSection._J; + } +} + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodMeshSection.cpp b/src/BeamAdapter/component/model/RodMeshSection.cpp new file mode 100644 index 000000000..58b506884 --- /dev/null +++ b/src/BeamAdapter/component/model/RodMeshSection.cpp @@ -0,0 +1,39 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_PLUGIN_BEAMADAPTER_RODMESHSECTION_CPP + +#include +#include +#include +#include + +namespace sofa::beamadapter +{ + +using namespace sofa::defaulttype; + +const int RodMeshSectionClass = core::RegisterObject("Class defining a Rod Section using a MeshLoader and material parameters.") + .add< RodMeshSection >(true); + +template class SOFA_BEAMADAPTER_API RodMeshSection; + +}// namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodMeshSection.h b/src/BeamAdapter/component/model/RodMeshSection.h new file mode 100644 index 000000000..0c75b9c8e --- /dev/null +++ b/src/BeamAdapter/component/model/RodMeshSection.h @@ -0,0 +1,86 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace sofa::beamadapter +{ + +using sofa::core::loader::MeshLoader; + +/** + * \class RodMeshSection + * \brief Specialization class of @sa BaseRodSectionMaterial describing a rod section created using a Mesh file + * + * This class will describe a rod section defined by a mesh file structure using the link @sa l_loader + * Method @sa initFromLoader and @sa initRestConfig will define the beam structure using the geometry of the given mesh + * as well as the Length. Mechanical parameters are set using the @sa BaseRodSectionMaterial Data + * Method @sa getRestTransformOnX will return the current position of the curviline abscisse along the mesh structure. + */ +template +class RodMeshSection : public sofa::beamadapter::BaseRodSectionMaterial +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(RodMeshSection, DataTypes), SOFA_TEMPLATE(BaseRodSectionMaterial, DataTypes)); + + /// Default Constructor + RodMeshSection(); + + /// Override method to get the rest position of the beam. In this implementation, it will interpolate along the loaded mesh geometry + void getRestTransformOnX(Transform& global_H_local, const Real& x_used, const Real& x_start) override; + +protected: + /// Internal method to init the section. Called by @sa BaseRodSectionMaterial::init() method + bool initSection() override; + + /// Internal method called by initSection to init from a linked MeshLoader @sa l_loader + bool initFromLoader(); + /// Internal method called by initFromLoader to compute @sa m_localRestPositions and @sa m_localRestTransforms given the mesh structure + void initRestConfig(); + /// Method to check if the given loader has a edge set structure + bool checkLoaderTopology(); + + /// Tool method to rotate the input frame @param input given an axis @param x. Result is set in @param output + void rotateFrameForAlignX(const Quat& input, Vec3& x, Quat& output); + +public: + /// Link to be set to the topology container in the component graph. + SingleLink, MeshLoader, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_loader; + +private: + /// Pointer to the MeshLoader, should be set using @sa l_loader, otherwise will search for one in current Node. + MeshLoader* p_loader{ nullptr }; + + type::vector m_localRestPositions; ///< rest position of the key points interpolated on the mesh geometry + type::vector m_localRestTransforms; ///< rest transform of the key points interpolated on the mesh geometry + type::vector m_curvAbs; ///< set of absciss curviline points + Real m_absOfGeometry{ 0 }; ///< max curv absciss of this mesh structure +}; + +#if !defined(SOFA_PLUGIN_BEAMADAPTER_RODMESHSECTION_CPP) +extern template class SOFA_BEAMADAPTER_API RodMeshSection; +#endif + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodMeshSection.inl b/src/BeamAdapter/component/model/RodMeshSection.inl new file mode 100644 index 000000000..1867c195b --- /dev/null +++ b/src/BeamAdapter/component/model/RodMeshSection.inl @@ -0,0 +1,298 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +#define EPSILON 0.0001 + +namespace sofa::beamadapter +{ + +template +RodMeshSection::RodMeshSection() + : BaseRodSectionMaterial() + , l_loader(initLink("loader", "link to the MeshLoader")) +{ + +} + + +template +bool RodMeshSection::initSection() +{ + // Get meshLoader, check first if loader has been set using link. Otherwise will search in current context. + p_loader = l_loader.get(); + + if (!p_loader) + this->getContext()->get(p_loader); + + if (!p_loader) { + msg_error() << "Cannot find a mesh loader. Please insert a MeshObjLoader in the same node or use l_loader to specify the path in the scene graph."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return false; + } + else + { + msg_info() << "Found a mesh with " << p_loader->d_edges.getValue().size() << " edges"; + return initFromLoader(); + } +} + + +template +void RodMeshSection::getRestTransformOnX(Transform& global_H_local, const Real& x_used, const Real& x_start) +{ + Real abs_curr = x_used - x_start; + abs_curr = abs_curr /(d_length.getValue()) * m_absOfGeometry; + + Coord p; + + /*** find the range which includes the "requested" abs ***/ + Real startingAbs = 0; + unsigned int index = 0; + while ((startingAbs < abs_curr) && (index < m_localRestPositions.size())) + { + index++; + startingAbs = m_curvAbs[index]; + } + + /*** OOB ***/ + if (abs_curr > startingAbs) + { + msg_error() << "Out of bound position requested= " << abs_curr << " with startingAbs = " << startingAbs; + return; + } + else /*** Expected case ***/ + { + const Real alpha = (abs_curr - m_curvAbs[index - 1]) / (m_curvAbs[index] - m_curvAbs[index - 1]); + const Real one_minus_alpha = 1 - alpha; + const Vec3 result = m_localRestTransforms[index - 1].getOrigin() * one_minus_alpha + m_localRestTransforms[index].getOrigin() * alpha; + Quat slerp; + slerp.slerp(m_localRestTransforms[index - 1].getOrientation(), m_localRestTransforms[index].getOrientation(), alpha, true); + slerp.normalize(); + + p.getCenter() = result; + p.getOrientation() = slerp; + } + + Vec3 PosEndCurve = p.getCenter(); + Quat ExtremityQuat = p.getOrientation(); + Vec3 ExtremityPos = PosEndCurve + Vec3(x_start, 0, 0); + + global_H_local.set(ExtremityPos, ExtremityQuat); +} + + +template +bool RodMeshSection::initFromLoader() +{ + if (!checkLoaderTopology()) + { + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return false; + } + + type::vector vertices; + sofa::core::topology::BaseMeshTopology::SeqEdges edges; + + //get the topology position + auto topoVertices = sofa::helper::getReadAccessor(p_loader->d_positions); + + //copy the topology edges in a local vector + auto topoEdges = sofa::helper::getReadAccessor(p_loader->d_edges); + edges = topoEdges.ref(); + + /** renumber the vertices **/ + type::vector verticesConnexion; //gives the number of edges connected to a vertex + for (unsigned int i = 0; i < topoVertices.size(); i++) + verticesConnexion.push_back(2); + + for (const auto& ed : edges) + { + verticesConnexion[ed[0]]--; + verticesConnexion[ed[1]]--; + } + + msg_info() << "Successfully compute the vertex connexion"; + + // check for the first corner of the edge + unsigned int firstIndex = 0; + bool found = false; + while ((firstIndex < verticesConnexion.size()) && !found) + { + if (verticesConnexion[firstIndex] == 1) + found = true; + else + firstIndex++; + } + + if (firstIndex == verticesConnexion.size()) + { + msg_error() << "The first vertex of the beam structure is not found, probably because of a closed structure"; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return false; + } + + vertices.push_back(topoVertices[firstIndex]); + + while (edges.size() > 0) + { + auto it = edges.begin(); + auto end = edges.end(); + + bool notFound = true; + while (notFound && (it != end)) + { + const auto& ed = (*it); + auto toDel = it; + it++; + if (ed[0] == firstIndex) + { + vertices.push_back(topoVertices[ed[1]]); + firstIndex = ed[1]; + edges.erase(toDel); + notFound = false; + + } + else if (ed[1] == firstIndex) + { + vertices.push_back(topoVertices[ed[0]]); + firstIndex = ed[0]; + edges.erase(toDel); + notFound = false; + } + } + } + + msg_info() << "Successfully computed the topology"; + + m_localRestPositions = vertices; + + initRestConfig(); + + return true; +} + + +template +void RodMeshSection::initRestConfig() +{ + m_curvAbs.clear(); + Real tot = 0; + m_curvAbs.push_back(0); + Quat input, output; + input.identity(); + m_localRestTransforms.resize(m_localRestPositions.size()); + m_localRestTransforms[0].setOrigin(Vec3(0, 0, 0)); + m_localRestTransforms[0].setOrientation(input); + + for (unsigned int i = 0; i < m_localRestPositions.size() - 1; i++) + { + Vec3 vec = m_localRestPositions[i + 1] - m_localRestPositions[i]; + Real norm = vec.norm(); + tot += norm; + + this->rotateFrameForAlignX(input, vec, output); + + input = output; + + m_localRestTransforms[i + 1].setOrientation(output); + + Vec3 localPos = m_localRestPositions[i + 1] - m_localRestPositions[0]; + + m_localRestTransforms[i + 1].setOrigin(localPos); + + m_curvAbs.push_back(tot); + } + m_absOfGeometry = tot; + + d_length.setValue(m_absOfGeometry); + + msg_info() << "Length of the loaded shape = " << m_absOfGeometry; +} + + +template +bool RodMeshSection::checkLoaderTopology() +{ + if (!p_loader->d_edges.getValue().size()) + { + msg_error() << "There is no edges in the topology loaded by " << p_loader->getName(); + return false; + } + + if (p_loader->d_triangles.getValue().size()) + { + msg_error() << "There are triangles in the topology loaded by " << p_loader->getName(); + return false; + } + + if (p_loader->d_quads.getValue().size()) + { + msg_error() << "There are quads in the topology loaded by " << p_loader->getName(); + return false; + } + + if (p_loader->d_polygons.getValue().size()) + { + msg_error() << "There are polygons in the topology loaded by " << p_loader->getName(); + return false; + } + + return true; +} + + +template +void RodMeshSection::rotateFrameForAlignX(const Quat& input, Vec3& x, Quat& output) +{ + x.normalize(); + Vec3 x0 = input.inverseRotate(x); + + Real cTheta = x0[0]; + Real theta; + if (cTheta > (1 - EPSILON)) + { + output = input; + } + else + { + theta = acos(cTheta); + // axis of rotation + Vec3 dw(0, -x0[2], x0[1]); + dw.normalize(); + + // computation of the rotation + Quat inputRoutput; + inputRoutput.axisToQuat(dw, theta); + + output = input * inputRoutput; + } +} + + + + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodSpireSection.cpp b/src/BeamAdapter/component/model/RodSpireSection.cpp new file mode 100644 index 000000000..04b64d83a --- /dev/null +++ b/src/BeamAdapter/component/model/RodSpireSection.cpp @@ -0,0 +1,39 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_PLUGIN_BEAMADAPTER_RODSPIRESECTION_CPP + +#include +#include +#include +#include + +namespace sofa::beamadapter +{ + +using namespace sofa::defaulttype; + +const int RodSpireSectionClass = core::RegisterObject("Class defining a rod spire section, defining material and geometry parameters.") + .add< RodSpireSection >(true); + +template class SOFA_BEAMADAPTER_API RodSpireSection; + +}// namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodSpireSection.h b/src/BeamAdapter/component/model/RodSpireSection.h new file mode 100644 index 000000000..91a793d88 --- /dev/null +++ b/src/BeamAdapter/component/model/RodSpireSection.h @@ -0,0 +1,65 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + +namespace sofa::beamadapter +{ + +using sofa::core::loader::MeshLoader; + +/** + * \class RodSpireSection + * \brief Specialization class of @sa BaseRodSectionMaterial describing a rod spire section. + * + * This class will describe a rod spire section using spire diameter and height between each spire. Length and mechanical + * parameters are the same as @sa BaseRodSectionMaterial Data + * Method @sa getRestTransformOnX will return the current position of the curviline abscisse along the spire. + */ +template +class RodSpireSection : public sofa::beamadapter::BaseRodSectionMaterial +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(RodSpireSection, DataTypes), SOFA_TEMPLATE(BaseRodSectionMaterial, DataTypes)); + + /// Default Constructor + RodSpireSection(); + + /// Override method to get the rest position of the beam. In this implementation, it will compute the current position given the spire parameters + void getRestTransformOnX(Transform& global_H_local, const Real& x_used, const Real& x_start) override; + +protected: + /// Internal method to init the section. Called by @sa BaseRodSectionMaterial::init() method + bool initSection() override; + +public: + Data d_spireDiameter; ///< Data defining the diameter of the spire + Data d_spireHeight; ///< Data defining the height between each spire +}; + +#if !defined(SOFA_PLUGIN_BEAMADAPTER_RODSPIRESECTION_CPP) +extern template class SOFA_BEAMADAPTER_API RodSpireSection; +#endif + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodSpireSection.inl b/src/BeamAdapter/component/model/RodSpireSection.inl new file mode 100644 index 000000000..cb9ff3230 --- /dev/null +++ b/src/BeamAdapter/component/model/RodSpireSection.inl @@ -0,0 +1,97 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace sofa::beamadapter +{ + +template +RodSpireSection::RodSpireSection() + : BaseRodSectionMaterial() + , d_spireDiameter(initData(&d_spireDiameter, (Real)0.1, "spireDiameter", "diameter of the spire")) + , d_spireHeight(initData(&d_spireHeight, (Real)0.01, "spireHeight", "height between each spire")) +{ + +} + + +template +bool RodSpireSection::initSection() +{ + const auto length = this->d_length.getValue(); + if (length <= Real(0.0)) + { + msg_error() << "Length is 0 (or negative), check if d_length has been given or well computed."; + return false; + } + + if (int nbrEdgesVisu = d_nbEdgesVisu.getValue() <= 0) + { + msg_warning() << "Number of visual edges has been set to an invalid value: " << nbrEdgesVisu << ". Value should be a positive integer. Setting to default value: 10"; + d_nbEdgesVisu.setValue(10); + } + + + if (int nbEdgesCollis = d_nbEdgesCollis.getValue() <= 0) + { + msg_warning() << "Number of collision edges has been set to an invalid value: " << nbEdgesCollis << ". Value should be a positive integer. Setting to default value: 20"; + d_nbEdgesCollis.setValue(10); + } + + return true; +} + + +template +void RodSpireSection::getRestTransformOnX(Transform& global_H_local, const Real& x_used, const Real& x_start) +{ + Real projetedLength = d_spireDiameter.getValue() * M_PI; + Real lengthSpire = sqrt(d_spireHeight.getValue() * d_spireHeight.getValue() + projetedLength * projetedLength); + // angle in the z direction + Real phi = atan(d_spireHeight.getValue() / projetedLength); + + Quat Qphi; + Qphi.axisToQuat(Vec3(0, 0, 1), phi); + + // spire angle (if theta=2*PI, there is a complete spire between startx and x) + Real lengthCurve = x_used - x_start; + Real numSpire = lengthCurve / lengthSpire; + Real theta = 2 * M_PI * numSpire; + + // computation of the Quat + Quat Qtheta; + Qtheta.axisToQuat(Vec3(0, 1, 0), theta); + Quat newSpireQuat = Qtheta * Qphi; + + // computation of the position + Real radius = d_spireDiameter.getValue() / 2.0; + Vec3 PosEndCurve(radius * sin(theta), numSpire * d_spireHeight.getValue(), radius * (cos(theta) - 1)); + Vec3 SpirePos = PosEndCurve + Vec3(x_start, 0, 0); + + global_H_local.set(SpirePos, newSpireQuat); +} + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodStraightSection.cpp b/src/BeamAdapter/component/model/RodStraightSection.cpp new file mode 100644 index 000000000..777466bf0 --- /dev/null +++ b/src/BeamAdapter/component/model/RodStraightSection.cpp @@ -0,0 +1,39 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_PLUGIN_BEAMADAPTER_RODSTRAIGHTSECTION_CPP + +#include +#include +#include +#include + +namespace sofa::beamadapter +{ + +using namespace sofa::defaulttype; + +const int RodStraightSectionClass = core::RegisterObject("Class defining a rod straight section Material, defining material and geometry parameters.") + .add< RodStraightSection >(true); + +template class SOFA_BEAMADAPTER_API RodStraightSection; + +}// namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodStraightSection.h b/src/BeamAdapter/component/model/RodStraightSection.h new file mode 100644 index 000000000..629c946fb --- /dev/null +++ b/src/BeamAdapter/component/model/RodStraightSection.h @@ -0,0 +1,58 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + +namespace sofa::beamadapter +{ + +/** + * \class RodStraightSection + * \brief Specialization class of @sa BaseRodSectionMaterial describing a rod straight section. + * + * This class will describe a rod straight section which will parametrized only using the @sa BaseRodSectionMaterial Data + * Method @sa getRestTransformOnX will return: Vec3(current_x, 0 0) + */ +template +class RodStraightSection : public sofa::beamadapter::BaseRodSectionMaterial +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(RodStraightSection, DataTypes), SOFA_TEMPLATE(BaseRodSectionMaterial, DataTypes)); + + /// Default Constructor + RodStraightSection(); + + /// Override method to get the rest position of the beam. In this implementation, it will basically returns Vec3(x_start + x_used, 0 0) + void getRestTransformOnX(Transform& global_H_local, const Real& x_used, const Real& x_start) override; + +protected: + /// Internal method to init the section. Called by @sa BaseRodSectionMaterial::init() method + bool initSection() override; +}; + +#if !defined(SOFA_PLUGIN_BEAMADAPTER_RODSTRAIGHTSECTION_CPP) +extern template class SOFA_BEAMADAPTER_API RodStraightSection; +#endif + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/component/model/RodStraightSection.inl b/src/BeamAdapter/component/model/RodStraightSection.inl new file mode 100644 index 000000000..cc499ca9b --- /dev/null +++ b/src/BeamAdapter/component/model/RodStraightSection.inl @@ -0,0 +1,73 @@ +/****************************************************************************** +* BeamAdapter plugin * +* (c) 2006 Inria, University of Lille, CNRS * +* * +* 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: see Authors.md * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +namespace sofa::beamadapter +{ + +template +RodStraightSection::RodStraightSection() + : BaseRodSectionMaterial() +{ + +} + + +template +bool RodStraightSection::initSection() +{ + const auto length = this->d_length.getValue(); + if (length <= Real(0.0)) + { + msg_error() << "Length is 0 (or negative), check if d_length has been given or well computed."; + return false; + } + + if (int nbrEdgesVisu = d_nbEdgesVisu.getValue() <= 0) + { + msg_warning() << "Number of visual edges has been set to an invalid value: " << nbrEdgesVisu << ". Value should be a positive integer. Setting to default value: 10"; + d_nbEdgesVisu.setValue(10); + } + + + if (int nbEdgesCollis = d_nbEdgesCollis.getValue() <= 0) + { + msg_warning() << "Number of collision edges has been set to an invalid value: " << nbEdgesCollis << ". Value should be a positive integer. Setting to default value: 20"; + d_nbEdgesCollis.setValue(10); + } + + return true; +} + + +template +void RodStraightSection::getRestTransformOnX(Transform& global_H_local, const Real& x_used, const Real& x_start) +{ + global_H_local.set(Vec3(x_start + x_used, 0.0, 0.0), Quat()); +} + + +} // namespace sofa::beamadapter diff --git a/src/BeamAdapter/gpu/cuda/CudaInstantiations.cpp b/src/BeamAdapter/gpu/cuda/CudaInstantiations.cpp index 5441f7f7d..94079918e 100644 --- a/src/BeamAdapter/gpu/cuda/CudaInstantiations.cpp +++ b/src/BeamAdapter/gpu/cuda/CudaInstantiations.cpp @@ -28,6 +28,9 @@ #include #include #include +#include +#include +#include #include #include @@ -86,6 +89,19 @@ namespace sofa::component::mapping #endif } // namespace sofa::component::mapping +namespace sofa::beamadapter +{ + template class SOFA_BEAMADAPTER_API RodMeshSection; + template class SOFA_BEAMADAPTER_API RodSpireSection; + template class SOFA_BEAMADAPTER_API RodStraightSection; + +#ifdef SOFA_GPU_CUDA_DOUBLE + template class SOFA_BEAMADAPTER_API RodMeshSection; + template class SOFA_BEAMADAPTER_API RodSpireSection; + template class SOFA_BEAMADAPTER_API RodStraightSection; +#endif +} // namespace sofa::beamadapter + namespace sofa::gpu::cuda { @@ -139,6 +155,29 @@ int CudaMultiAdaptiveBeamMappingClass = core::RegisterObject("Set the positions #endif ; +const int CudaRodMeshSectionClass = core::RegisterObject("Class defining a Rod Section using a MeshLoader and material parameters using CUDA.") + .add< sofa::beamadapter::RodMeshSection >() +#ifdef SOFA_GPU_CUDA_DOUBLE + .add< sofa::beamadapter::RodMeshSection >() +#endif + ; + +const int CudaRodSpireSectionClass = core::RegisterObject("Class defining a rod spire section, defining material and geometry parameters using CUDA.") + .add< sofa::beamadapter::RodSpireSection >() +#ifdef SOFA_GPU_CUDA_DOUBLE + .add< sofa::beamadapter::RodSpireSection >() +#endif +; + +const int CudaRodStraightSectionClass = core::RegisterObject("Class defining a rod straight section Material, defining material and geometry parameters using CUDA.") + .add< sofa::beamadapter::RodStraightSection >() +#ifdef SOFA_GPU_CUDA_DOUBLE + .add< sofa::beamadapter::RodStraightSection >() +#endif +; + + + } // namespace sofa::gpu::cuda