From f8d85dc74f3b953fb872a9be81360d2bd2cf4d1c Mon Sep 17 00:00:00 2001 From: EulalieCoevoet Date: Thu, 21 May 2026 13:09:03 +0200 Subject: [PATCH 1/2] [forcefield] TriangularShellForceField: draw element showing the measured values (color map) --- src/Shell/engine/FindClosePoints.inl | 2 +- .../forcefield/TriangularShellForceField.h | 6 ++ .../forcefield/TriangularShellForceField.inl | 102 +++++++++++++----- 3 files changed, 83 insertions(+), 27 deletions(-) diff --git a/src/Shell/engine/FindClosePoints.inl b/src/Shell/engine/FindClosePoints.inl index 88b989c..5c870e1 100644 --- a/src/Shell/engine/FindClosePoints.inl +++ b/src/Shell/engine/FindClosePoints.inl @@ -59,7 +59,7 @@ void FindClosePoints::doUpdate() for (Index i=0; i(i,j)); + list.push_back(type::fixed_array{i,j}); } } } diff --git a/src/Shell/forcefield/TriangularShellForceField.h b/src/Shell/forcefield/TriangularShellForceField.h index 3caafd0..c1e23bc 100644 --- a/src/Shell/forcefield/TriangularShellForceField.h +++ b/src/Shell/forcefield/TriangularShellForceField.h @@ -37,6 +37,7 @@ #include #include +#include // Uncomment the following to use quaternions instead of matrices for @@ -214,12 +215,17 @@ class TriangularShellForceField : public core::behavior::ForceField Data d_corotated; Data d_measure; Data > d_measuredValues; + Data d_showMeasuredValuesColorMap; Data d_isShellveryThin; Data d_use_rest_position; Data d_arrow_radius; TRQSTriangleHandler* triangleHandler; + helper::ColorMap m_colorMap = helper::ColorMap(256, "Blue to Red"); + Real m_minMeasuredValue; + Real m_maxMeasuredValue; + protected : // Selected elements diff --git a/src/Shell/forcefield/TriangularShellForceField.inl b/src/Shell/forcefield/TriangularShellForceField.inl index fe379f7..4808d49 100644 --- a/src/Shell/forcefield/TriangularShellForceField.inl +++ b/src/Shell/forcefield/TriangularShellForceField.inl @@ -81,6 +81,7 @@ TriangularShellForceField::TriangularShellForceField() , d_corotated(initData(&d_corotated, true, "corotated", "Compute forces in corotational frame")) , d_measure(initData(&d_measure, "measure", "Compute the strain or stress")) , d_measuredValues(initData(&d_measuredValues, "measuredValues", "Measured values for stress or strain")) + , d_showMeasuredValuesColorMap(initData(&d_showMeasuredValuesColorMap, false, "showMeasuredValuesColorMap", "draw elements showing the measured values (color map from blue to red)")) , d_isShellveryThin(initData(&d_isShellveryThin, false, "isShellveryThin", "This bool is to adapt " "computation in case we are using verry tiny(thickness) shell element")) , d_use_rest_position(initData(&d_use_rest_position, true, "use_rest_position", "Use the rest position inteat of using postion to update the restposition")) @@ -217,6 +218,8 @@ template void TriangularShellForceField::reinit() { d_measuredValues.beginEdit()->resize(_topology->getNbPoints()); d_measuredValues.endEdit(); + m_minMeasuredValue = 0; + m_maxMeasuredValue = 0; } /// Prepare to store info in the triangle array @@ -702,26 +705,35 @@ void TriangularShellForceField::accumulateForce(VecDeriv &f, const Ve } // Compute the measure (stress/strain) - if (bMeasureStrain) { + if (bMeasureStrain || bMeasureStress) + { + m_minMeasuredValue = std::numeric_limits::max(); + m_maxMeasuredValue = std::numeric_limits::min(); type::vector &values = *d_measuredValues.beginEdit(); - for (unsigned int i=0; i< tinfo->measure.size(); i++) { - Vec3 strain = tinfo->measure[i].B * Dm + tinfo->measure[i].Bb * Db; - // Norm from strain in x and y - // NOTE: Shear strain is not included - values[ tinfo->measure[i].id ] = helper::rsqrt( - strain[0] * strain[0] + strain[1] * strain[1]); + + if (bMeasureStrain) { + for (unsigned int i=0; i< tinfo->measure.size(); i++) { + Vec3 strain = tinfo->measure[i].B * Dm + tinfo->measure[i].Bb * Db; + // Norm from strain in x and y + // NOTE: Shear strain is not included + values[ tinfo->measure[i].id ] = helper::rsqrt(strain[0] * strain[0] + strain[1] * strain[1]); + } + } else if (bMeasureStress) { + for (unsigned int i=0; i< tinfo->measure.size(); i++) { + Vec3 stress = materialMatrix * tinfo->measure[i].B * Dm + + materialMatrix * tinfo->measure[i].Bb * Db; + // Von Mises stress criterion (plane stress) + values[ tinfo->measure[i].id ] = helper::rsqrt(stress[0] * stress[0] - stress[0] * stress[1] + + stress[1] * stress[1] + 3 * stress[2] * stress[2]); + } } - d_measuredValues.endEdit(); - } else if (bMeasureStress) { - type::vector &values = *d_measuredValues.beginEdit(); - for (unsigned int i=0; i< tinfo->measure.size(); i++) { - Vec3 stress = materialMatrix * tinfo->measure[i].B * Dm - + materialMatrix * tinfo->measure[i].Bb * Db; - // Von Mises stress criterion (plane stress) - values[ tinfo->measure[i].id ] = helper::rsqrt( - stress[0] * stress[0] - stress[0] * stress[1] - + stress[1] * stress[1] + 3 * stress[2] * stress[2]); + + for (size_t i = 0; i < values.size(); i++) + { + m_minMeasuredValue = std::min(values[i], m_minMeasuredValue); + m_maxMeasuredValue = std::max(values[i], m_maxMeasuredValue); } + d_measuredValues.endEdit(); } @@ -1131,25 +1143,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 @@ -1158,7 +1170,7 @@ 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)); + 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}); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -1277,14 +1289,30 @@ void TriangularShellForceField::dktSD(StrainDisplacement &B, const Tr } template -void TriangularShellForceField::draw(const core::visual::VisualParams* vparams){ - // Draw arrows from using shell triangle info (tinfo.R) in the center of the triangle +void TriangularShellForceField::draw(const core::visual::VisualParams* vparams) +{ + if (!vparams->displayFlags().getShowForceFields()) + return; + + const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); + vparams->drawTool()->disableLighting(); + + if (vparams->displayFlags().getShowWireFrame()) + vparams->drawTool()->setPolygonMode(0, true); - int nbTriangles=_topology->getNbTriangles(); + unsigned int nbTriangles=_topology->getNbTriangles(); auto triangles = _topology->getTriangles(); const VecCoord& positions =this->mstate->read(sofa::core::vec_id::read_access::position)->getValue(); type::vector& triangleInf = *(triangleInfo.beginEdit()); const Real radius = d_arrow_radius.getValue(); + + std::vector colorVector; + std::vector vertices; + colorVector.reserve(nbTriangles * 3); + vertices.reserve(nbTriangles * 3); + + helper::ReadAccessor > > values = d_measuredValues; + for (unsigned int i=0; i< nbTriangles; i++) { const TriangleInformation &tinfo = triangleInf[i]; @@ -1300,12 +1328,34 @@ void TriangularShellForceField::draw(const core::visual::VisualParams // compute the center of the triangle Vec3 center = (positions[a].getCenter() + positions[b].getCenter() + positions[c].getCenter())/3; + // Draw arrows from using shell triangle info (tinfo.R) in the center of the triangle vparams->drawTool()->drawArrow(center, center + vx, radius, type::RGBAColor(1.0, 0.0, 0.0, 1.0)); vparams->drawTool()->drawArrow(center, center + vy, radius, type::RGBAColor(0.0, 1.0, 0.0, 1.0)); vparams->drawTool()->drawArrow(center, center + vz, radius, type::RGBAColor(0.0, 0.0, 1.0, 1.0)); + + if(d_showMeasuredValuesColorMap.getValue() & !values.empty()) + { + helper::ColorMap::evaluator evalColor = m_colorMap.getEvaluator(m_minMeasuredValue, m_maxMeasuredValue); + colorVector.push_back(evalColor(values[triangle[0]])); + colorVector.push_back(evalColor(values[triangle[1]])); + colorVector.push_back(evalColor(values[triangle[2]])); + } + else + { + colorVector.push_back(sofa::type::RGBAColor::green()); + colorVector.push_back(sofa::type::RGBAColor(0, 0.5, 0.5, 1)); + colorVector.push_back(sofa::type::RGBAColor(0, 0, 1, 1)); + } + + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[a]))); + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[b]))); + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[c]))); } + // Draw triangles + vparams->drawTool()->drawTriangles(vertices, colorVector); + triangleInfo.endEdit(); - } +} } // namespace forcefield From 6e584835ef3fe1c1e0e385131d5d1142d41f22c5 Mon Sep 17 00:00:00 2001 From: EulalieCoevoet Date: Thu, 21 May 2026 18:47:06 +0200 Subject: [PATCH 2/2] add option to draw the arrows --- .../forcefield/TriangularShellForceField.h | 1 + .../forcefield/TriangularShellForceField.inl | 32 +++++++++++-------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/Shell/forcefield/TriangularShellForceField.h b/src/Shell/forcefield/TriangularShellForceField.h index c1e23bc..d23534f 100644 --- a/src/Shell/forcefield/TriangularShellForceField.h +++ b/src/Shell/forcefield/TriangularShellForceField.h @@ -216,6 +216,7 @@ class TriangularShellForceField : public core::behavior::ForceField Data d_measure; Data > d_measuredValues; Data d_showMeasuredValuesColorMap; + Data d_showArrows; Data d_isShellveryThin; Data d_use_rest_position; Data d_arrow_radius; diff --git a/src/Shell/forcefield/TriangularShellForceField.inl b/src/Shell/forcefield/TriangularShellForceField.inl index 4808d49..22dac48 100644 --- a/src/Shell/forcefield/TriangularShellForceField.inl +++ b/src/Shell/forcefield/TriangularShellForceField.inl @@ -82,6 +82,7 @@ TriangularShellForceField::TriangularShellForceField() , d_measure(initData(&d_measure, "measure", "Compute the strain or stress")) , d_measuredValues(initData(&d_measuredValues, "measuredValues", "Measured values for stress or strain")) , d_showMeasuredValuesColorMap(initData(&d_showMeasuredValuesColorMap, false, "showMeasuredValuesColorMap", "draw elements showing the measured values (color map from blue to red)")) + , d_showArrows(initData(&d_showArrows, false, "showArrows", "Show the orientation of the triangle as arrow.")) , d_isShellveryThin(initData(&d_isShellveryThin, false, "isShellveryThin", "This bool is to adapt " "computation in case we are using verry tiny(thickness) shell element")) , d_use_rest_position(initData(&d_use_rest_position, true, "use_rest_position", "Use the rest position inteat of using postion to update the restposition")) @@ -1315,23 +1316,29 @@ void TriangularShellForceField::draw(const core::visual::VisualParams for (unsigned int i=0; i< nbTriangles; i++) { - const TriangleInformation &tinfo = triangleInf[i]; - Vec3 vx = tinfo.R * Vec3(1, 0, 0); - Vec3 vy = tinfo.R * Vec3(0, 1, 0); - Vec3 vz = tinfo.R * Vec3(0, 0, 1); - auto triangle = _topology->getTriangle(i); // get triangle i indices auto a = triangle[0]; auto b = triangle[1]; auto c = triangle[2]; - // compute the center of the triangle - Vec3 center = (positions[a].getCenter() + positions[b].getCenter() + positions[c].getCenter())/3; + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[a]))); + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[b]))); + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[c]))); + // Draw arrows from using shell triangle info (tinfo.R) in the center of the triangle - vparams->drawTool()->drawArrow(center, center + vx, radius, type::RGBAColor(1.0, 0.0, 0.0, 1.0)); - vparams->drawTool()->drawArrow(center, center + vy, radius, type::RGBAColor(0.0, 1.0, 0.0, 1.0)); - vparams->drawTool()->drawArrow(center, center + vz, radius, type::RGBAColor(0.0, 0.0, 1.0, 1.0)); + if (d_showArrows.getValue()) + { + // compute the center of the triangle + Vec3 center = (positions[a].getCenter() + positions[b].getCenter() + positions[c].getCenter())/3; + const TriangleInformation &tinfo = triangleInf[i]; + Vec3 vx = tinfo.R * Vec3(1, 0, 0); + Vec3 vy = tinfo.R * Vec3(0, 1, 0); + Vec3 vz = tinfo.R * Vec3(0, 0, 1); + vparams->drawTool()->drawArrow(center, center + vx, radius, type::RGBAColor(1.0, 0.0, 0.0, 1.0)); + vparams->drawTool()->drawArrow(center, center + vy, radius, type::RGBAColor(0.0, 1.0, 0.0, 1.0)); + vparams->drawTool()->drawArrow(center, center + vz, radius, type::RGBAColor(0.0, 0.0, 1.0, 1.0)); + } if(d_showMeasuredValuesColorMap.getValue() & !values.empty()) { @@ -1346,11 +1353,8 @@ void TriangularShellForceField::draw(const core::visual::VisualParams colorVector.push_back(sofa::type::RGBAColor(0, 0.5, 0.5, 1)); colorVector.push_back(sofa::type::RGBAColor(0, 0, 1, 1)); } - - vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[a]))); - vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[b]))); - vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(positions[c]))); } + // Draw triangles vparams->drawTool()->drawTriangles(vertices, colorVector);