Skip to content

[BeamInterpolation] Generalize parameters per edge #169

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Mar 26, 2025
Merged

[BeamInterpolation] Generalize parameters per edge #169

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
03d60b0
[BeamInterpolation] generalizes parameters per egde
EulalieCoevoet Feb 17, 2025
4889361
adds tests and fixes failing AdaptiveBeamForceFieldAndMassTest.checkV…
EulalieCoevoet Mar 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 71 additions & 1 deletion BeamAdapter_test/BeamInterpolation_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,8 @@ using sofa::component::statecontainer::MechanicalObject ;
using std::string;

#include <BeamAdapter/component/forcefield/AdaptiveBeamForceFieldAndMass.h>
#include <BeamAdapter/component/BeamInterpolation.h>
using sofa::component::fem::_beaminterpolation_::BeamInterpolation;

namespace sofa
{
Expand All @@ -71,6 +73,68 @@ struct BeamInterpolationTest : public sofa::testing::BaseSimulationTest,
sofa::simpleapi::importPlugin("BeamAdapter");
}

Node::SPtr createSingleBeam()
{
string scene =
"<?xml version='1.0'?>"
"<Node name='Root' gravity='0 -9.81 0' dt='0.01'>"
" <RequiredPlugin name='Sofa.Component.ODESolver.Backward' />"
" <RequiredPlugin name='Sofa.Component.LinearSolver.Direct' />"
" <RequiredPlugin name='Sofa.Component.Constraint.Projective' />"
" <RequiredPlugin name='Sofa.Component.StateContainer' />"
" <RequiredPlugin name='Sofa.Component.Topology.Container.Constant' />"
" <RequiredPlugin name='Sofa.Component.Topology.Container.Grid' />"
" <RequiredPlugin name='BeamAdapter' />"
" <DefaultAnimationLoop />"
" <DefaultVisualManagerLoop />"
" <Node name='BeamModel'/>"
" <EulerImplicitSolver rayleighStiffness='0.0' rayleighMass='0.0' />"
" <BTDLinearSolver />"
" <RegularGridTopology name='MeshLines' nx='200' ny='1' nz='1' xmax='100' xmin='0' ymin='0' ymax='0' zmax='0' zmin='0'/>"
" <MechanicalObject template='Rigid3d' name='DOFs' />"
" <FixedConstraint indices='0' />"
" <BeamInterpolation name='Interpol' radius='0.1'/>"
" </Node> "
"</Node> ";

Node::SPtr root = SceneLoaderXML::loadFromMemory("singleBeam", scene.c_str());
sofa::simulation::node::initRoot(root.get());

return root;
}

void checkDataInitialization(const Data<type::vector<SReal>>& data,
const SReal& defaultValue,
const sofa::Size& nbBeam,
const SReal& value)
{
const auto& vector = helper::getReadAccessor(data);
ASSERT_EQ(vector.size(), nbBeam);
for (const auto& v: vector)
ASSERT_FLOAT_EQ(v, value);
ASSERT_FLOAT_EQ(defaultValue, value);
}

void checkCreation()
{
Node::SPtr root = createSingleBeam();

// Search for Beam FF
BeamInterpolation<Rigid3dTypes>* beamInterpolation = nullptr;
root->getTreeObject(beamInterpolation);
ASSERT_NE(beamInterpolation, nullptr);

sofa::Size nbBeam = 199;

// Check component state and Data default values
checkDataInitialization(beamInterpolation->d_radius, beamInterpolation->m_defaultRadius, nbBeam, 0.1);
checkDataInitialization(beamInterpolation->d_innerRadius, beamInterpolation->m_defaultInnerRadius, nbBeam, 0.0);
checkDataInitialization(beamInterpolation->d_lengthY, beamInterpolation->m_defaultLengthY, nbBeam, 1.0);
checkDataInitialization(beamInterpolation->d_lengthZ, beamInterpolation->m_defaultLengthZ, nbBeam, 1.0);
checkDataInitialization(beamInterpolation->d_defaultYoungModulus, beamInterpolation->m_defaultYoungModulus, nbBeam, 1e5);
checkDataInitialization(beamInterpolation->d_poissonRatio, beamInterpolation->m_defaultPoissonRatio, nbBeam, 0.4);
}

void simpleScene(const std::vector<std::string>& lines)
{
assert(lines.size()==3);
Expand All @@ -96,7 +160,8 @@ struct BeamInterpolationTest : public sofa::testing::BaseSimulationTest,

root->init(core::ExecParams::defaultInstance());
root->reinit(core::ExecParams::defaultInstance()) ;
}else if(lines[2]=="W")
}
else if(lines[2]=="W")
{
EXPECT_MSG_EMIT(Error) ;
EXPECT_MSG_NOEMIT(Warning) ;
Expand All @@ -110,6 +175,7 @@ struct BeamInterpolationTest : public sofa::testing::BaseSimulationTest,
}
};


static std::vector<std::vector<std::string>> teststrings ={
{
"<MeshTopology name='meshSuture' edges='0 1' />"
Expand Down Expand Up @@ -143,6 +209,10 @@ TEST_P(BeamInterpolationTest, checkMinimalScene) {
ASSERT_NO_THROW(this->simpleScene(GetParam())) ;
}

TEST_F(BeamInterpolationTest, checkCreation) {
ASSERT_NO_THROW(this->checkCreation());
}

INSTANTIATE_TEST_SUITE_P(checkMinimalScene,
BeamInterpolationTest, ::testing::ValuesIn(teststrings) ) ;

Expand Down
7 changes: 6 additions & 1 deletion BeamAdapter_test/component/forcefield/AdaptiveBeamForceFieldAndMassTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,12 @@ struct AdaptiveBeamForceFieldAndMassTest : public sofa::testing::BaseSimulationT
// Check component state and Data default values
ASSERT_EQ(beamForceFieldMass->d_componentState.getValue(), sofa::core::objectmodel::ComponentState::Valid);
ASSERT_EQ(beamForceFieldMass->d_computeMass.getValue(), true);
ASSERT_FLOAT_EQ(beamForceFieldMass->d_massDensity.getValue(), 10.0);

const auto& massDensity = helper::getReadAccessor(beamForceFieldMass->d_massDensity);
ASSERT_EQ(massDensity.size(), 199);
for (const auto& m: massDensity)
ASSERT_FLOAT_EQ(m, 10.0);
ASSERT_FLOAT_EQ(beamForceFieldMass->m_defaultMassDensity, 10.0);
ASSERT_FLOAT_EQ(beamForceFieldMass->rayleighMass.getValue(), 0.0);
ASSERT_FLOAT_EQ(beamForceFieldMass->rayleighStiffness.getValue(), 0.0);
}
Expand Down
37 changes: 23 additions & 14 deletions src/BeamAdapter/component/BeamInterpolation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,29 +141,28 @@ class BeamInterpolation : public BaseBeamInterpolation<DataTypes>

const BeamSection& getBeamSection(sofa::Index beamId) override {
SOFA_UNUSED(beamId);
return this->m_constantSection;
if (m_constantBeam)
return this->m_constantSection;
return this->m_section[beamId];
}
void getInterpolationParameters(sofa::Index beamId, Real &_L, Real &_A, Real &_Iy , Real &_Iz,
Real &_Asy, Real &_Asz, Real &J) override;
void getMechanicalParameters(sofa::Index beamId, Real& youngModulus, Real& cPoisson, Real& massDensity) override;

void getTangentUsingSplinePoints(unsigned int edgeInList, const Real& baryCoord, const sofa::core::ConstVecCoordId &vecXId, Vec3& t );


/// computeActualLength => given the 4 control points of the spline, it provides an estimate of the length (using gauss points integration)




virtual void getCurvAbsAtBeam(const unsigned int& edgeInList_input, const Real& baryCoord_input, Real& x_output) {}
virtual void getBeamAtCurvAbs(const Real& x_input, unsigned int& edgeInList_output, Real& baryCoord_output, unsigned int start = 0) {}


Data<helper::OptionsGroup> crossSectionShape;

/// Circular Cross Section
Data<Real> d_radius;
Data<Real> d_innerRadius;
Real m_defaultRadius;
Data<type::vector<Real>> d_radius;
Real m_defaultInnerRadius;
Data<type::vector<Real>> d_innerRadius;

/// Square Cross Section
Data<Real> d_sideLength;
Expand All @@ -173,14 +172,17 @@ class BeamInterpolation : public BaseBeamInterpolation<DataTypes>
Data<Real> d_largeRadius;

/// Rectangular Cross Section
Data<Real> d_lengthY;
Data<Real> d_lengthZ;
Real m_defaultLengthY;
Data<type::vector<Real>> d_lengthY;
Real m_defaultLengthZ;
Data<type::vector<Real>> d_lengthZ;
Data<bool> d_dofsAndBeamsAligned;

Real m_defaultYoungModulus;
Real m_defaultPoissonRatio;
Data<type::vector<Real>> d_defaultYoungModulus;
Data<type::vector<Real>> d_poissonRatio;
Real m_defaultYoungModulus;
Data<type::vector<Real>> d_defaultYoungModulus;

Real m_defaultPoissonRatio;
Data<type::vector<Real>> d_poissonRatio;

Data<bool> d_straight;

Expand Down Expand Up @@ -214,7 +216,14 @@ protected :
Data< VecDeriv > d_InterpolatedVel;

/// GEOMETRICAL COMPUTATION (for now we suppose that the radius of the beam do not vary in space / in time)
bool m_constantBeam{true};
BeamSection m_constantSection;
type::vector<BeamSection> m_section;

void checkDataSize(Real& defaultValue, Data<type::vector<Real>>& dataList, const size_t& nbEdges);

void computeRectangularCrossSectionInertiaMatrix(const Real &Ly, const Real &Lz, BeamSection &section);
void computeCircularCrossSectionInertiaMatrix(const Real &r, const Real &rInner, BeamSection &section);
};

#if !defined(SOFA_PLUGIN_BEAMADAPTER_BEAMINTERPOLATION_CPP)
Expand Down
Loading
Loading