From dc9a7ccd507ee6b9c6a0abbd957e305ae62c24c1 Mon Sep 17 00:00:00 2001 From: Tian Sizhe Date: Wed, 9 Apr 2025 16:16:56 +0200 Subject: [PATCH 1/3] update draw function to draw shell as surface, and visualize the stress measured on the triangle as heatmap --- .../forcefield/TriangularShellForceField.h | 13 +- .../forcefield/TriangularShellForceField.inl | 173 +++++++++++++++--- 2 files changed, 157 insertions(+), 29 deletions(-) diff --git a/src/SofaShells/forcefield/TriangularShellForceField.h b/src/SofaShells/forcefield/TriangularShellForceField.h index 3caafd0..df77e6f 100644 --- a/src/SofaShells/forcefield/TriangularShellForceField.h +++ b/src/SofaShells/forcefield/TriangularShellForceField.h @@ -38,6 +38,7 @@ #include #include +#include // Uncomment the following to use quaternions instead of matrices for // rotations. Quaternions are slightly faster but numericaly quite unstable @@ -99,6 +100,12 @@ class TriangularShellForceField : public core::behavior::ForceField class TriangleInformation; + // Flags for the visualization + Data d_showArrow; + Data d_arrowRadius; + Data d_showTriangle; + Data d_showMeasuredValue; + protected: typedef Vec<9, Real> Displacement; // the displacement vector @@ -216,10 +223,14 @@ class TriangularShellForceField : public core::behavior::ForceField Data > d_measuredValues; Data d_isShellveryThin; Data d_use_rest_position; - Data d_arrow_radius; TRQSTriangleHandler* triangleHandler; +private: + sofa::helper::ColorMap* p_drawColorMap; ///< colormap to display the gradiant of stress if @sa showStressValue is set to true + Real m_minStress = 0; ///< min stress computed for @sa showStressValue + Real m_maxStress = 0; ///< max stress computed for @sa showStressValue + protected : // Selected elements diff --git a/src/SofaShells/forcefield/TriangularShellForceField.inl b/src/SofaShells/forcefield/TriangularShellForceField.inl index ff019be..15c1371 100644 --- a/src/SofaShells/forcefield/TriangularShellForceField.inl +++ b/src/SofaShells/forcefield/TriangularShellForceField.inl @@ -27,6 +27,10 @@ #include #include +#include +#include +#include + #include #include #include @@ -38,11 +42,12 @@ #include #include #include -#include #include #include #include +#include + #ifdef _WIN32 #include #endif @@ -85,8 +90,10 @@ TriangularShellForceField::TriangularShellForceField() "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")) , triangleInfo(initData(&triangleInfo, "triangleInfo", "Internal triangle data")) - , d_arrow_radius(initData(&d_arrow_radius, (Real)0.1, "arrow_radius", "the arrow radius")) - + , d_showArrow(initData(&d_showArrow, false, "showArrow", "Show the orientation of the triangle as arrow")) + , d_arrowRadius(initData(&d_arrowRadius, (Real)0.1, "arrow_radius", "the arrow radius")) + , d_showTriangle(initData(&d_showTriangle, false, "showTriangle", "Show the triangle")) + , d_showMeasuredValue(initData(&d_showMeasuredValue, false, "showMeasuredValue", "Show the measured value")) { d_membraneElement.beginEdit()->setNames(7, "None", // No membrane element @@ -117,6 +124,8 @@ TriangularShellForceField::TriangularShellForceField() d_measure.endEdit(); triangleHandler = new TRQSTriangleHandler(this, &triangleInfo); + + p_drawColorMap = new helper::ColorMap(256, "Blue to Red"); } @@ -701,6 +710,8 @@ void TriangularShellForceField::accumulateForce(VecDeriv &f, const Ve } + m_minStress = std::numeric_limits::max(); + m_maxStress = std::numeric_limits::lowest(); // Compute the measure (stress/strain) if (bMeasureStrain) { type::vector &values = *d_measuredValues.beginEdit(); @@ -710,6 +721,10 @@ void TriangularShellForceField::accumulateForce(VecDeriv &f, const Ve // NOTE: Shear strain is not included values[ tinfo->measure[i].id ] = helper::rsqrt( strain[0] * strain[0] + strain[1] * strain[1]); + if (values[ tinfo->measure[i].id ] < m_minStress) + m_minStress = values[ tinfo->measure[i].id ]; + if (values[ tinfo->measure[i].id ] > m_maxStress) + m_maxStress = values[ tinfo->measure[i].id ]; } d_measuredValues.endEdit(); } else if (bMeasureStress) { @@ -721,6 +736,10 @@ void TriangularShellForceField::accumulateForce(VecDeriv &f, const Ve values[ tinfo->measure[i].id ] = helper::rsqrt( stress[0] * stress[0] - stress[0] * stress[1] + stress[1] * stress[1] + 3 * stress[2] * stress[2]); + if (values[ tinfo->measure[i].id ] < m_minStress) + m_minStress = values[ tinfo->measure[i].id ]; + if (values[ tinfo->measure[i].id ] > m_maxStress) + m_maxStress = values[ tinfo->measure[i].id ]; } d_measuredValues.endEdit(); } @@ -1278,33 +1297,131 @@ void TriangularShellForceField::dktSD(StrainDisplacement &B, const Tr template void TriangularShellForceField::draw(const core::visual::VisualParams* vparams){ + helper::ScopedAdvancedTimer advancedTimer( "draw_TriangularShellForceField" ); + + if (!vparams->displayFlags().getShowForceFields()) { + return; + } + + // a flag indicating if there are infomation to draw + bool p_computeDrawInfo; + p_computeDrawInfo = d_showArrow.getValue() || + d_showTriangle.getValue() || + d_showMeasuredValue.getValue(); + + if (!p_computeDrawInfo) { + return; + } + + + const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); + + if (vparams->displayFlags().getShowWireFrame()) + vparams->drawTool()->setPolygonMode(0, true); + + vparams->drawTool()->disableLighting(); + + const VecCoord& x = this->mstate->read(core::ConstVecCoordId::position())->getValue(); + const type::vector& triangleInf = triangleInfo.getValue(); + const auto& triangles = _topology->getTriangles(); + const Size nbTriangles = triangles.size(); + // Draw arrows from using shell triangle info (tinfo.R) in the center of the triangle + if (d_showArrow.getValue()) + { + const VecCoord& positions =this->mstate->read(sofa::core::ConstVecCoordId::position())->getValue(); + const Real radius = d_arrowRadius.getValue(); + 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; + 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)); + } + triangleInfo.endEdit(); + } + + // Draw triangles using shell triangle info + if (d_showTriangle.getValue()) + { + std::vector vertices; + for (Size i = 0; i < nbTriangles; ++i) + { + const Triangle& tri = triangles[i]; + Index a = tri[0]; + Index b = tri[1]; + Index c = tri[2]; + + auto vertex_vector = x[a]; + auto vertex_position = sofa::type::Vec3(vertex_vector[0], vertex_vector[1], vertex_vector[2]); + vertices.push_back(sofa::type::Vec3(vertex_position)); + vertex_vector = x[b]; + vertex_position = sofa::type::Vec3(vertex_vector[0], vertex_vector[1], vertex_vector[2]); + vertices.push_back(sofa::type::Vec3(vertex_position)); + vertex_vector = x[c]; + vertex_position = sofa::type::Vec3(vertex_vector[0], vertex_vector[1], vertex_vector[2]); + vertices.push_back(sofa::type::Vec3(vertex_position)); + } + vparams->drawTool()->drawTriangles(vertices, sofa::type::RGBAColor(1, 0, 0, 0.9)); + vertices.clear(); + } + + // Visualize the stress on the shell + if (d_showMeasuredValue.getValue()) + { + std::vector vertices; + std::vector colorVector; + + auto evalColor = p_drawColorMap->getEvaluator(m_minStress, m_maxStress); + const auto& measuredValues = sofa::helper::getReadAccessor(d_measuredValues); + + if(measuredValues.size()){ + for (Size i = 0; i < nbTriangles; ++i) + { + const Triangle& tri = triangles[i]; + Index a = tri[0]; + Index b = tri[1]; + Index c = tri[2]; + + auto vertex_vector = x[a]; + auto vertex_position = sofa::type::Vec3(vertex_vector[0], vertex_vector[1], vertex_vector[2]); + auto measured_value = measuredValues[ triangleInf[i].measure[0].id ]; + colorVector.push_back(evalColor(measured_value)); + vertices.push_back(vertex_position); + + vertex_vector = x[b]; + vertex_position = sofa::type::Vec3(vertex_vector[0], vertex_vector[1], vertex_vector[2]); + measured_value = measuredValues[ triangleInf[i].measure[1].id ]; + colorVector.push_back(evalColor(measured_value)); + vertices.push_back(vertex_position); + + vertex_vector = x[c]; + vertex_position = sofa::type::Vec3(vertex_vector[0], vertex_vector[1], vertex_vector[2]); + measured_value = measuredValues[ triangleInf[i].measure[2].id ]; + colorVector.push_back(evalColor(measured_value)); + vertices.push_back(vertex_position); + + } + vparams->drawTool()->drawTriangles(vertices, colorVector); + vertices.clear(); + colorVector.clear(); + } + + } + - int nbTriangles=_topology->getNbTriangles(); - auto triangles = _topology->getTriangles(); - const VecCoord& positions =this->mstate->read(sofa::core::ConstVecCoordId::position())->getValue(); - type::vector& triangleInf = *(triangleInfo.beginEdit()); - const Real radius = d_arrow_radius.getValue(); - 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; - 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)); - } - triangleInfo.endEdit(); } From 1a97dfea9a132bc0229b14122d880510b2174287 Mon Sep 17 00:00:00 2001 From: Tian Sizhe Date: Thu, 10 Apr 2025 14:20:05 +0200 Subject: [PATCH 2/3] use a data member as color map instead of a pointer --- src/SofaShells/forcefield/TriangularShellForceField.h | 2 +- src/SofaShells/forcefield/TriangularShellForceField.inl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/SofaShells/forcefield/TriangularShellForceField.h b/src/SofaShells/forcefield/TriangularShellForceField.h index df77e6f..148c4bd 100644 --- a/src/SofaShells/forcefield/TriangularShellForceField.h +++ b/src/SofaShells/forcefield/TriangularShellForceField.h @@ -227,7 +227,7 @@ class TriangularShellForceField : public core::behavior::ForceField TRQSTriangleHandler* triangleHandler; private: - sofa::helper::ColorMap* p_drawColorMap; ///< colormap to display the gradiant of stress if @sa showStressValue is set to true + sofa::helper::ColorMap d_drawColorMap; ///< colormap to display the gradiant of stress if @sa showStressValue is set to true Real m_minStress = 0; ///< min stress computed for @sa showStressValue Real m_maxStress = 0; ///< max stress computed for @sa showStressValue diff --git a/src/SofaShells/forcefield/TriangularShellForceField.inl b/src/SofaShells/forcefield/TriangularShellForceField.inl index 15c1371..8d844ff 100644 --- a/src/SofaShells/forcefield/TriangularShellForceField.inl +++ b/src/SofaShells/forcefield/TriangularShellForceField.inl @@ -125,7 +125,7 @@ TriangularShellForceField::TriangularShellForceField() triangleHandler = new TRQSTriangleHandler(this, &triangleInfo); - p_drawColorMap = new helper::ColorMap(256, "Blue to Red"); + d_drawColorMap = helper::ColorMap(256, "Blue to Red"); } @@ -1384,7 +1384,7 @@ void TriangularShellForceField::draw(const core::visual::VisualParams std::vector vertices; std::vector colorVector; - auto evalColor = p_drawColorMap->getEvaluator(m_minStress, m_maxStress); + auto evalColor = d_drawColorMap.getEvaluator(m_minStress, m_maxStress); const auto& measuredValues = sofa::helper::getReadAccessor(d_measuredValues); if(measuredValues.size()){ From 7b6a73fc8d5b02507dcbfe6a2e2d9abfbd6c2702 Mon Sep 17 00:00:00 2001 From: Tian Sizhe Date: Thu, 10 Apr 2025 14:20:57 +0200 Subject: [PATCH 3/3] Update src/SofaShells/forcefield/TriangularShellForceField.inl use SCOPED_TIMER marco for Tracy support Co-authored-by: Alex Bilger --- src/SofaShells/forcefield/TriangularShellForceField.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SofaShells/forcefield/TriangularShellForceField.inl b/src/SofaShells/forcefield/TriangularShellForceField.inl index 15c1371..b2d86bf 100644 --- a/src/SofaShells/forcefield/TriangularShellForceField.inl +++ b/src/SofaShells/forcefield/TriangularShellForceField.inl @@ -1297,7 +1297,7 @@ void TriangularShellForceField::dktSD(StrainDisplacement &B, const Tr template void TriangularShellForceField::draw(const core::visual::VisualParams* vparams){ - helper::ScopedAdvancedTimer advancedTimer( "draw_TriangularShellForceField" ); + SCOPED_TIMER("draw_TriangularShellForceField"); if (!vparams->displayFlags().getShowForceFields()) { return;