Merge branch 'master' into pr-add-field-to-surface-mesh

damienmarchal · web-flow · commit b769e33c333b · 2025-06-30T16:10:42.000+02:00
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
@@ -0,0 +1,7 @@
+version: 2
+updates:
+  - package-ecosystem: github-actions
+    directory: /
+    target-branch: devel
+    schedule:
+      interval: monthly
diff --git a/CMakePresets.json b/CMakePresets.json
@@ -155,6 +155,10 @@
           "type": "BOOL",
           "value": "ON"
         },
+        "PLUGIN_BEAMADAPTER_CUDA": {
+          "type": "BOOL",
+          "value": "ON"
+        },
         "PLUGIN_SOFAMATRIX": {
           "type": "BOOL",
           "value": "ON"
@@ -171,6 +175,10 @@
           "type": "BOOL",
           "value": "ON"
         },
+        "PLUGIN_SOFADISTANCEGRID_CUDA": {
+          "type": "BOOL",
+          "value": "ON"
+        },
         "PLUGIN_GEOMAGIC": {
           "type": "BOOL",
           "value": "ON"
@@ -272,6 +280,10 @@
           "type": "BOOL",
           "value": "ON"
         },
+        "PLUGIN_SOFASPHFLUID_CUDA": {
+          "type": "BOOL",
+          "value": "ON"
+        },
         "PLUGIN_SOFASPHFLUID": {
           "type": "BOOL",
           "value": "ON"
diff --git a/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.inl b/Sofa/Component/Constraint/Lagrangian/Model/src/sofa/component/constraint/lagrangian/model/BilateralLagrangianConstraint.inl
@@ -277,6 +277,8 @@ void BilateralLagrangianConstraint<DataTypes>::getConstraintResolution(const Con
                                                                         std::vector<ConstraintResolution*>& resTab,
                                                                         unsigned int& offset)
 {
+    if (!d_activate.getValue()) return;
+    
     SOFA_UNUSED(cParams);
     const unsigned minp=std::min(d_m1.getValue().size(), d_m2.getValue().size());
 
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FastTetrahedralCorotationalForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/FastTetrahedralCorotationalForceField.inl
@@ -51,7 +51,7 @@ void FastTetrahedralCorotationalForceField<DataTypes>::createTetrahedronRestInfo
 
     auto [lambda, mu] = Inherited::toLameParameters(_3DMat, youngModulusElement, poissonRatioElement);
 
-    typename DataTypes::Real volume,val;
+    typename DataTypes::Real val;
     typename DataTypes::Coord point[4]; //shapeVector[4];
     const typename DataTypes::VecCoord restPosition=this->mstate->read(core::vec_id::read_access::restPosition)->getValue();
 
@@ -63,20 +63,20 @@ void FastTetrahedralCorotationalForceField<DataTypes>::createTetrahedronRestInfo
     // store the point position
     for(j=0; j<4; ++j)
         point[j]=(restPosition)[t[j]];
-    /// compute 6 times the rest volume
-    volume=dot(cross(point[1]-point[0],point[2]-point[0]),point[0]-point[3]);
+
+    const auto tetrahedronVolume = -geometry::Tetrahedron::signedVolume(point[0],point[1],point[2],point[3]);
     /// store the rest volume
-    my_tinfo.restVolume=volume/6;
-    mu*=fabs(volume)/6;
-    lambda*=fabs(volume)/6;
+    my_tinfo.restVolume = tetrahedronVolume;
+    mu *= fabs(tetrahedronVolume);
+    lambda *= fabs(tetrahedronVolume);
 
     // store shape vectors at the rest configuration
     for(j=0; j<4; ++j)
     {
         if ((j%2)==0)
-            my_tinfo.shapeVector[j]=cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/volume;
+            my_tinfo.shapeVector[j] = cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/(tetrahedronVolume * 6);
         else
-            my_tinfo.shapeVector[j]= -cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/volume;
+            my_tinfo.shapeVector[j] = -cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/(tetrahedronVolume * 6);
     }
 
     /// compute the vertex stiffness of the linear elastic material, needed for addKToMatrix
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TetrahedralCorotationalFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TetrahedralCorotationalFEMForceField.inl
@@ -390,21 +390,10 @@ void TetrahedralCorotationalFEMForceField<DataTypes>::computeMaterialStiffness(i
 
     // divide by 36 times volumes of the element
     const VecCoord& X0=this->mstate->read(core::vec_id::read_access::restPosition)->getValue();
+    const auto tetrahedronVolume = geometry::Tetrahedron::volume(
+            X0[a], X0[b], X0[c], X0[d]);
 
-    Coord A = (X0)[b] - (X0)[a];
-    Coord B = (X0)[c] - (X0)[a];
-    Coord C = (X0)[d] - (X0)[a];
-    Coord AB = cross(A, B);
-    Real volumes6 = fabs( dot( AB, C ) );
-    if (volumes6<0)
-    {
-        msg_error() << "Negative volume for tetra " << a << ',' << b << ',' << c << ',' << d << "> = " << volumes6 / 6;
-    }
-    //	materialMatrix  /= (volumes6);//*6 christian
-    // @TODO: in TetrahedronFEMForceField, the stiffness matrix is divided by 6 compared to the code in TetrahedralCorotationalFEMForceField. Check which is the correct one...
-    // FF:  there is normally  a factor 1/6v in the strain-displacement matrix. Times transpose makes 1/36v². Integrating across the volume multiplies by v, so the factor is 1/36v
-    materialMatrix  /= (volumes6*6);
-
+   materialMatrix /= tetrahedronVolume * 36;
 }
 
 template<class DataTypes>
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TetrahedronFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TetrahedronFEMForceField.inl
@@ -283,18 +283,11 @@ void TetrahedronFEMForceField<DataTypes>::computeMaterialStiffness(Index i, Inde
 
     // divide by 36 times volumes of the element
     const VecCoord &initialPoints=d_initialPoints.getValue();
-    Coord A = initialPoints[b] - initialPoints[a];
-    Coord B = initialPoints[c] - initialPoints[a];
-    Coord C = initialPoints[d] - initialPoints[a];
-    Coord AB = cross(A, B);
-    Real volumes6 = fabs( dot( AB, C ) );
+    const auto tetrahedronVolume = geometry::Tetrahedron::volume(
+        initialPoints[a], initialPoints[b], initialPoints[c], initialPoints[d]);
 
-    m_restVolume += volumes6/6;
-    if (volumes6<0)
-    {
-        msg_error() << "Negative volume for tetra "<<i<<" <"<<a<<','<<b<<','<<c<<','<<d<<"> = "<<volumes6/6 ;
-    }
-    materialsStiffnesses[i] /= (volumes6*6); // 36*Volume in the formula
+    m_restVolume += tetrahedronVolume;
+    materialsStiffnesses[i] /= tetrahedronVolume * 36;
 }
 
 template<class DataTypes>
diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TriangleFEMForceField.inl b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/TriangleFEMForceField.inl
@@ -204,7 +204,8 @@ void TriangleFEMForceField<DataTypes>::computeMaterialStiffnesses()
         const Real Estrain = this->getYoungModulusInElement(i) / ((1 + _p) * (1 - 2 * _p));
         const Real Estress = this->getYoungModulusInElement(i) / (1 - _p * _p);
 
-        const Real triangleVolume = (Real)0.5 * d_thickness.getValue() * cross(p[b] - p[a], p[c] - p[a]).norm();
+        const Real triangleArea = sofa::geometry::Triangle::area(p[a], p[b], p[c]);
+        const Real triangleVolume = d_thickness.getValue() * triangleArea;
 
         if (d_planeStrain.getValue() == true)
         {
diff --git a/Sofa/Component/Visual/CMakeLists.txt b/Sofa/Component/Visual/CMakeLists.txt
@@ -15,6 +15,7 @@ set(HEADER_FILES
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/TrailRenderer.h
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/TrailRenderer.inl
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/Visual3DText.h
+    ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualBoundingBox.h
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualGrid.h
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualModelImpl.h
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualStyle.h
@@ -31,6 +32,7 @@ set(SOURCE_FILES
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/RecordedCamera.cpp
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/TrailRenderer.cpp
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/Visual3DText.cpp
+    ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualBoundingBox.cpp
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualGrid.cpp
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualModelImpl.cpp
     ${SOFACOMPONENTVISUAL_SOURCE_DIR}/VisualStyle.cpp
diff --git a/Sofa/Component/Visual/src/sofa/component/visual/VisualBoundingBox.cpp b/Sofa/Component/Visual/src/sofa/component/visual/VisualBoundingBox.cpp
@@ -0,0 +1,53 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+
+#include <sofa/component/visual/VisualBoundingBox.h>
+#include <sofa/core/visual/VisualParams.h>
+#include <sofa/core/ObjectFactory.h>
+#include <sofa/defaulttype/VecTypes.h>
+
+namespace sofa::component::visual
+{
+
+void registerVisualBoundingBox(sofa::core::ObjectFactory* factory)
+{
+    factory->registerObjects(core::ObjectRegistrationData("Display an Axis Aligned Bounding Box (AABB).")
+        .add< VisualBoundingBox >());
+}
+
+VisualBoundingBox::VisualBoundingBox()
+    : d_color(initData(&d_color, sofa::type::RGBAColor::yellow(),  "color", "Color of the lines of the box."))
+    , d_thickness(initData(&d_thickness, 1.0f,  "thickness", "Thickness of the lines of the box."))
+{
+
+}
+
+void VisualBoundingBox::doDrawVisual(const core::visual::VisualParams* vparams)
+{
+    
+    const auto& bbox = f_bbox.getValue();
+    vparams->drawTool()->disableLighting();
+    vparams->drawTool()->setMaterial(d_color.getValue());
+    vparams->drawTool()->drawBoundingBox(bbox.minBBox(), bbox.maxBBox(), d_thickness.getValue());
+}
+
+} // namespace sofa::component::visual
diff --git a/Sofa/Component/Visual/src/sofa/component/visual/VisualBoundingBox.h b/Sofa/Component/Visual/src/sofa/component/visual/VisualBoundingBox.h
@@ -0,0 +1,49 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2006 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+#include <sofa/component/visual/config.h>
+
+#include <sofa/core/visual/VisualModel.h>
+#include <sofa/type/RGBAColor.h>
+
+namespace sofa::component::visual
+{
+/*
+ * Display an Axis Orientated Bounding Box.
+ */
+class SOFA_COMPONENT_VISUAL_API VisualBoundingBox : public core::visual::VisualModel
+{
+public:
+    SOFA_CLASS(VisualBoundingBox, core::visual::VisualModel);
+
+    VisualBoundingBox();
+    ~VisualBoundingBox() override = default;
+
+    Data<sofa::type::RGBAColor> d_color; ///< Color of the lines in the grid. default=yellow
+    Data<float> d_thickness; ///< Thickness of the lines in the grid
+
+private:
+    void doDrawVisual(const core::visual::VisualParams*) override;
+    
+};
+
+} // namespace sofa::component::visual
diff --git a/Sofa/Component/Visual/src/sofa/component/visual/VisualStyle.cpp b/Sofa/Component/Visual/src/sofa/component/visual/VisualStyle.cpp
@@ -60,6 +60,16 @@ VisualStyle::VisualStyle()
     d_displayFlags.setWidget("widget_displayFlags");
 }
 
+void VisualStyle::init()
+{
+    updateVisualFlags(sofa::core::visual::VisualParams::defaultInstance());
+}
+
+void VisualStyle::bwdInit()
+{
+    applyBackupFlags(sofa::core::visual::VisualParams::defaultInstance());
+}
+
 void VisualStyle::updateVisualFlags(VisualParams* vparams)
 {
     backupFlags = vparams->displayFlags();
diff --git a/Sofa/Component/Visual/src/sofa/component/visual/VisualStyle.h b/Sofa/Component/Visual/src/sofa/component/visual/VisualStyle.h
@@ -66,6 +66,8 @@ class SOFA_COMPONENT_VISUAL_API VisualStyle : public sofa::core::visual::BaseVis
 protected:
     VisualStyle();
 public:
+    void init() override;
+    void bwdInit() override;
     void updateVisualFlags(VisualParams* ) override;
     void applyBackupFlags(VisualParams* ) override;
 
diff --git a/Sofa/Component/Visual/src/sofa/component/visual/init.cpp b/Sofa/Component/Visual/src/sofa/component/visual/init.cpp
@@ -33,6 +33,7 @@ extern void registerLineAxis(sofa::core::ObjectFactory* factory);
 extern void registerRecordedCamera(sofa::core::ObjectFactory* factory);
 extern void registerTrailRenderer(sofa::core::ObjectFactory* factory);
 extern void registerVisual3DText(sofa::core::ObjectFactory* factory);
+extern void registerVisualBoundingBox(sofa::core::ObjectFactory* factory);
 extern void registerVisualGrid(sofa::core::ObjectFactory* factory);
 extern void registerVisualModelImpl(sofa::core::ObjectFactory* factory);
 extern void registerVisualStyle(sofa::core::ObjectFactory* factory);
@@ -69,6 +70,7 @@ void registerObjects(sofa::core::ObjectFactory* factory)
     registerRecordedCamera(factory);
     registerTrailRenderer(factory);
     registerVisual3DText(factory);
+    registerVisualBoundingBox(factory);
     registerVisualGrid(factory);
     registerVisualModelImpl(factory);
     registerVisualStyle(factory);
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Tetrahedron.h b/Sofa/framework/Geometry/src/sofa/geometry/Tetrahedron.h
@@ -44,15 +44,33 @@ struct Tetrahedron
     * @remark	This function is not generic
     * @tparam   Node a container of the type sofa::type::Vec3 (needed for cross(), dot(), operator-)
     * @tparam   T scalar
-    * @param	n0,n1,n2,n3,n4 nodes of the tetrahedron
-    * @return	Volume of the hexahedron (a T scalar)
+    * @param	n0,n1,n2,n3 nodes of the tetrahedron
+    * @return	Volume of the tetrahedron (a T scalar)
     */
     template<typename Node,
              typename T = std::decay_t<decltype(*std::begin(std::declval<Node>()))>,
              typename = std::enable_if_t<std::is_scalar_v<T>>
     >
     [[nodiscard]]
     static constexpr auto volume(const Node& n0, const Node& n1, const Node& n2, const Node& n3)
+    {
+        return std::abs(signedVolume(n0, n1, n2, n3));
+    }
+
+    /**
+    * @brief	Compute the signed volume of a tetrahedron
+    * @remark	This function is not generic
+    * @tparam   Node a container of the type sofa::type::Vec3 (needed for cross(), dot(), operator-)
+    * @tparam   T scalar
+    * @param	n0,n1,n2,n3 nodes of the tetrahedron
+    * @return	Signed volume of the tetrahedron (a T scalar)
+    */
+    template<typename Node,
+             typename T = std::decay_t<decltype(*std::begin(std::declval<Node>()))>,
+             typename = std::enable_if_t<std::is_scalar_v<T>>
+    >
+    [[nodiscard]]
+    static constexpr auto signedVolume(const Node& n0, const Node& n1, const Node& n2, const Node& n3)
     {
         constexpr Node n{};
         static_assert(std::distance(std::begin(n), std::end(n)) == 3, "volume can only be computed in 3 dimensions.");
@@ -61,8 +79,7 @@ struct Tetrahedron
         const auto b = n2 - n0;
         const auto c = n3 - n0;
 
-        return std::abs(sofa::type::dot(sofa::type::cross(a, b), c) / static_cast<T>(6));
-
+        return sofa::type::dot(sofa::type::cross(a, b), c) / static_cast<T>(6);
     }
 };
 
diff --git a/applications/plugins/SofaCUDA/src/sofa/gpu/cuda/CudaTypes.h b/applications/plugins/SofaCUDA/src/sofa/gpu/cuda/CudaTypes.h
@@ -215,13 +215,13 @@ class Vec3r1 : public sofa::type::Vec<3,Real>
     typedef sofa::type::Vec<3,Real> Inherit;
     typedef Real real;
     enum { N=3 };
-    Vec3r1() : dummy((Real) 0.0) {}
+    constexpr Vec3r1() : dummy((Real) 0.0) {}
     template<class real2>
-    Vec3r1(const Vec<N,real2>& v): Inherit(v), dummy((Real) 0.0) {}
-    Vec3r1(real x, real y, real z) : Inherit(x,y,z), dummy((Real) 0.0) {}
+    constexpr Vec3r1(const Vec<N,real2>& v): Inherit(v), dummy((Real) 0.0) {}
+    constexpr Vec3r1(real x, real y, real z) : Inherit(x,y,z), dummy((Real) 0.0) {}
 
     /// Fast constructor: no initialization
-    explicit Vec3r1(NoInit n) : Inherit(n), dummy((Real) 0.0)
+    explicit constexpr Vec3r1(NoInit n) : Inherit(n), dummy((Real) 0.0)
     {
     }
 
diff --git a/examples/Component/Visual/VisualBoundingBox.scn b/examples/Component/Visual/VisualBoundingBox.scn
@@ -0,0 +1,10 @@
+<Node name="root" dt="0.01">
+    <RequiredPlugin name="Sofa.Component.IO.Mesh"/> <!-- Needed to use components [MeshOBJLoader] -->
+    <RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualBoundingBox] -->
+    <RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> <!-- Needed to use components [OglModel] -->
+    <DefaultAnimationLoop/>
+    
+    <MeshOBJLoader name="meshLoader_0" filename="mesh/liver-smooth.obj" handleSeams="1" />
+    <OglModel name="VisualModel" src="@meshLoader_0" color="red" />
+    <VisualBoundingBox color="green" thickness="20" bbox="@VisualModel.bbox"/>
+</Node>

Original file line number	Diff line number	Diff line change
`@@ -277,6 +277,8 @@ void BilateralLagrangianConstraint<DataTypes>::getConstraintResolution(const Con`
`277`	`277`	`std::vector<ConstraintResolution*>& resTab,`
`278`	`278`	`unsigned int& offset)`
`279`	`279`	`{`
	`280`	`+ if (!d_activate.getValue()) return;`
	`281`	`+`
`280`	`282`	`SOFA_UNUSED(cParams);`
`281`	`283`	`const unsigned minp=std::min(d_m1.getValue().size(), d_m2.getValue().size());`
`282`	`284`
Original file line number	Diff line number	Diff line change
`@@ -204,7 +204,8 @@ void TriangleFEMForceField<DataTypes>::computeMaterialStiffnesses()`
`204`	`204`	`const Real Estrain = this->getYoungModulusInElement(i) / ((1 + _p) * (1 - 2 * _p));`
`205`	`205`	`const Real Estress = this->getYoungModulusInElement(i) / (1 - _p * _p);`
`206`	`206`
`207`		`- const Real triangleVolume = (Real)0.5 * d_thickness.getValue() * cross(p[b] - p[a], p[c] - p[a]).norm();`
	`207`	`+ const Real triangleArea = sofa::geometry::Triangle::area(p[a], p[b], p[c]);`
	`208`	`+ const Real triangleVolume = d_thickness.getValue() * triangleArea;`
`208`	`209`
`209`	`210`	`if (d_planeStrain.getValue() == true)`
`210`	`211`	`{`