Skip to content

[forcefield] TriangularShellForceField: draw element showing the measured values (color map) #19

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

Open
EulalieCoevoet wants to merge 2 commits into
base: main
Choose a base branch
from
+98 −37
Open

[forcefield] TriangularShellForceField: draw element showing the measured values (color map) #19

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
f8d85dc
[forcefield] TriangularShellForceField: draw element showing the meas…
EulalieCoevoet May 21, 2026
6e58483
add option to draw the arrows
EulalieCoevoet May 21, 2026
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
2 changes: 1 addition & 1 deletion src/Shell/engine/FindClosePoints.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ void FindClosePoints<DataTypes>::doUpdate()
for (Index i=0; i<points.size()-1; i++) {
for (Index j=i+1; j<points.size(); j++) {
if ((points[i] - points[j]).norm() <= threshold) {
list.push_back(type::fixed_array<Index,2>(i,j));
list.push_back(type::fixed_array<Index,2>{i,j});
}
}
}
Expand Down
7 changes: 7 additions & 0 deletions src/Shell/forcefield/TriangularShellForceField.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@

#include <sofa/core/topology/BaseMeshTopology.h>
#include <sofa/core/topology/TopologyData.h>
#include <sofa/helper/ColorMap.h>


// Uncomment the following to use quaternions instead of matrices for
Expand Down Expand Up @@ -214,12 +215,18 @@ class TriangularShellForceField : public core::behavior::ForceField<DataTypes>
Data<bool> d_corotated;
Data<sofa::helper::OptionsGroup> d_measure;
Data<type::vector<Real> > d_measuredValues;
Data<bool> d_showMeasuredValuesColorMap;
Data<bool> d_showArrows;
Data<bool> d_isShellveryThin;
Data<bool> d_use_rest_position;
Data<Real> d_arrow_radius;

TRQSTriangleHandler* triangleHandler;

helper::ColorMap m_colorMap = helper::ColorMap(256, "Blue to Red");
Real m_minMeasuredValue;
Real m_maxMeasuredValue;

protected :

// Selected elements
Expand Down
126 changes: 90 additions & 36 deletions src/Shell/forcefield/TriangularShellForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,8 @@ TriangularShellForceField<DataTypes>::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_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"))
Expand Down Expand Up @@ -217,6 +219,8 @@ template <class DataTypes> void TriangularShellForceField<DataTypes>::reinit()
{
d_measuredValues.beginEdit()->resize(_topology->getNbPoints());
d_measuredValues.endEdit();
m_minMeasuredValue = 0;
m_maxMeasuredValue = 0;
}

/// Prepare to store info in the triangle array
Expand Down Expand Up @@ -702,26 +706,35 @@ void TriangularShellForceField<DataTypes>::accumulateForce(VecDeriv &f, const Ve
}

// Compute the measure (stress/strain)
if (bMeasureStrain) {
if (bMeasureStrain || bMeasureStress)
{
m_minMeasuredValue = std::numeric_limits<float>::max();
m_maxMeasuredValue = std::numeric_limits<float>::min();
type::vector<Real> &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<Real> &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();
}

Expand Down Expand Up @@ -1131,25 +1144,25 @@ void TriangularShellForceField<DataTypes>::andesTemplate(StiffnessMatrix &K, con
template <class DataTypes>
void TriangularShellForceField<DataTypes>::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 <class DataTypes>
void TriangularShellForceField<DataTypes>::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 <class DataTypes>
void TriangularShellForceField<DataTypes>::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 <class DataTypes>
void TriangularShellForceField<DataTypes>::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
Expand All @@ -1158,7 +1171,7 @@ template <class DataTypes>
void TriangularShellForceField<DataTypes>::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});
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Expand Down Expand Up @@ -1277,35 +1290,76 @@ void TriangularShellForceField<DataTypes>::dktSD(StrainDisplacement &B, const Tr
}

template <class DataTypes>
void TriangularShellForceField<DataTypes>::draw(const core::visual::VisualParams* vparams){
// Draw arrows from using shell triangle info (tinfo.R) in the center of the triangle
void TriangularShellForceField<DataTypes>::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<TriangleInformation>& triangleInf = *(triangleInfo.beginEdit());
const Real radius = d_arrow_radius.getValue();

std::vector<sofa::type::RGBAColor> colorVector;
std::vector<sofa::type::Vec3> vertices;
colorVector.reserve(nbTriangles * 3);
vertices.reserve(nbTriangles * 3);

helper::ReadAccessor<Data<type::vector<Real> > > values = d_measuredValues;

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));
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
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())
{
helper::ColorMap::evaluator<Real> 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));
}
}

// Draw triangles
vparams->drawTool()->drawTriangles(vertices, colorVector);

triangleInfo.endEdit();
}
}


} // namespace forcefield
Expand Down
Loading