[Prefabs] Add modifiers

bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_BaseComponent.h

@@ -1,4 +1,4 @@
-/******************************************************************************
+  /******************************************************************************
 *                 SOFA, Simulation Open-Framework Architecture                *
 *                    (c) 2021 INRIA, USTL, UJF, CNRS, MGH                     *
 *                                                                             *

bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_Node.cpp

@@ -619,6 +619,17 @@ py::object getForceField(Node *self, unsigned int index)
 }
 
 
+
+py::object getLinearSolver(Node *self, unsigned int index)
+{
+    auto* ls = self->linearSolver.get(index);
+    if (ls) {
+        return PythonFactory::toPython(sofa::core::castToBase(ls));
+    }
+    return py::none();
+}
+
+
 py::object getMechanicalState(Node *self)
 {
     sofa::core::behavior::BaseMechanicalState* state = self->mechanicalState.get();
@@ -721,6 +732,7 @@ void moduleAddNode(py::module &m) {
     p.def("getMass", &getMass, sofapython3::doc::sofa::core::Node::getMass);
     p.def("hasODESolver", &hasODESolver, sofapython3::doc::sofa::core::Node::hasODESolver);
     p.def("getForceField", &getForceField, sofapython3::doc::sofa::core::Node::getForceField);
+    p.def("getLinearSolver", &getLinearSolver, sofapython3::doc::sofa::core::Node::getLinearSolver);
     p.def("getMechanicalState", &getMechanicalState, sofapython3::doc::sofa::core::Node::getMechanicalState);
     p.def("getMechanicalMapping", &getMechanicalMapping, sofapython3::doc::sofa::core::Node::getMechanicalMapping);
     p.def("sendEvent", &sendEvent, sofapython3::doc::sofa::core::Node::sendEvent);

bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_Node_doc.h

@@ -373,6 +373,13 @@ static auto getForceField =
         :type index: unsigned int.
         )";
 
+static auto getLinearSolver =
+        R"(
+        Get the linear solver of a node, given an index.
+        :param index: index of the linear solver
+        :type index: unsigned int.
+        )";
+
 static auto getMechanicalState =
         R"(
         Get the mechanical state of the node.

examples/stlib/SofaScene.py

@@ -1,5 +1,3 @@
-from fontTools.afmLib import preferredAttributeOrder
-
 from stlib.geometries.plane import PlaneParameters
 from stlib.geometries.file import FileParameters
 from stlib.geometries.extract import ExtractParameters
@@ -99,7 +97,7 @@ def createScene(root):
     SParams.material.parameters = [200, 0.45]
 
     def SAddMaterial(node):
-        DeformableBehaviorParameters.addDeformableMaterial(node)
+        DeformableBehaviorParameters.addMaterial(node)
         #TODO deal with that is a more smooth way in the material directly
         node.addObject("LinearSolverConstraintCorrection", name="ConstraintCorrection", linearSolver=SNode.LinearSolver.linkpath, ODESolver=SNode.ODESolver.linkpath)
 

examples/stlib/node_modifier.py

@@ -0,0 +1,100 @@
+from stlib.geometries.plane import PlaneParameters
+from stlib.geometries.file import FileParameters
+from stlib.geometries.extract import ExtractParameters
+from stlib.materials.deformable import DeformableBehaviorParameters
+from stlib.collision import Collision, CollisionParameters
+from stlib.entities import Entity, EntityParameters
+from stlib.visual import Visual, VisualParameters
+from splib.core.enum_types import CollisionPrimitive, ElementType, ConstitutiveLaw
+from splib.simulation.headers import setupLagrangianCollision, setupDefaultHeader
+from splib.simulation.ode_solvers import addImplicitODE
+from splib.simulation.linear_solvers import addLinearSolver
+import dataclasses
+import numpy as np
+
+
+from stlib.node_modifiers import NodeModifier
+from stlib.node_modifiers.footers import SimulationSolversParameters, SimulationSettingsParameters
+from stlib.node_modifiers.attachments import FixConstraintParameters, AttachmentConstraintParameters
+
+def createScene(root):
+    root.gravity=[0,0,9.81]
+    root.dt=0.01
+
+    ##Environement
+    planes_lengthNormal = np.array([0, 1, 0])
+    planes_lengthNbEdge = 1
+    planes_widthNbEdge = 2
+    planes_lengthSize = 30
+    planes_widthSize = 70
+
+    plane1_collisionParams = CollisionParameters()
+    plane1_collisionParams.name = "UP"
+    plane1_collisionParams.primitives = [CollisionPrimitive.TRIANGLES]
+    plane1_collisionParams.kwargs = {"TriangleCollision" : {"moving" : False, "simulated" : False}}
+    plane1_collisionParams.geometry = PlaneParameters(np.array([15,0,5]), np.array([0,0,-1]),
+                                                      planes_lengthNormal, planes_lengthNbEdge, planes_widthNbEdge, planes_lengthSize, planes_widthSize)
+    plane1 = root.add(Collision, parameters = plane1_collisionParams)
+    # TODO being able to reuse already loaded geometry of current prefab to add any new sub prefab
+    # We need to enable to directly pass a link to an already existing prefab in place of a prefab parameter object
+    plane1_visu = plane1.addChild("Visu")
+    plane1_visu.addObject("OglModel", name="VisualModel", src="@../Geometry/container")
+
+
+    ### Logo
+    ModelsNode = root.addChild("ModelsNode")
+
+    LogoParams = EntityParameters(name = "Logo1",
+                                  geometry = FileParameters(filename="mesh/SofaScene/Logo.vtk"),
+                                  material = DeformableBehaviorParameters(),
+                                  collision = CollisionParameters(geometry = FileParameters(filename="mesh/SofaScene/LogoColli.sph")),
+                                  visual = VisualParameters(geometry = FileParameters(filename="mesh/SofaScene/LogoVisu.obj")))
+
+    LogoParams.geometry.elementType = ElementType.TETRAHEDRA
+    LogoParams.material.constitutiveLawType = ConstitutiveLaw.ELASTIC
+    LogoParams.material.parameters = [200, 0.4]
+    LogoParams.material.massDensity = 0.003261
+    LogoParams.collision.primitives = [CollisionPrimitive.SPHERES]
+    #TODO make this flawless with spheres. Here collisions elements are not in the topology and a link is to be made between the loader and the collision object
+    LogoParams.collision.kwargs = {"SphereCollision" : {"radius" : "@Geometry/loader.listRadius"}}
+    LogoParams.visual.color = [0.7, .35, 0, 0.8]
+
+    Logo = ModelsNode.add(Entity, parameters = LogoParams)
+
+
+
+    ### S
+    SParams = EntityParameters()
+    SParams.name = "S"
+    SParams.geometry = FileParameters(filename="mesh/SofaScene/S.vtk")
+    SParams.geometry.elementType = ElementType.TETRAHEDRA
+    SParams.material = DeformableBehaviorParameters()
+    SParams.material.constitutiveLawType = ConstitutiveLaw.ELASTIC
+    SParams.material.parameters = [200, 0.45]
+    SParams.material.massDensity = 0.011021
+    SParams.collision = CollisionParameters()
+    SParams.collision.primitives = [CollisionPrimitive.TRIANGLES]
+    # # #TODO: to fix link issues for extracted geometry, it might be better to give source geometry relative link + parameters
+    ## TODO: not working with static container because the init order is always wrong so that the triangle vector is always empty when initializing the container
+    SParams.collision.geometry = ExtractParameters(destinationType=ElementType.TRIANGLES, sourceParameters=SParams.geometry)
+    SParams.visual = VisualParameters()
+    SParams.visual.geometry = FileParameters(filename="mesh/SofaScene/SVisu.obj")
+    SParams.visual.color = [0.7, .7, 0.7, 0.8]
+
+    S = ModelsNode.add(Entity, parameters = SParams)
+
+    #TODO make the name automatically match the modifier type if none is given
+    root.add(NodeModifier, on = [ModelsNode], parameters = SimulationSolversParameters(constantSparsity=False))
+
+    root.add(NodeModifier, on = [root], parameters = SimulationSettingsParameters(displayFlags = ["showVisualModels", "showInteractionForceFields"],
+                                                                                  enableCollisionDetection = True,
+                                                                                  useLagrangian = True,
+                                                                                  parallelComputing = False,
+                                                                                  alarmDistance=0.3, contactDistance=0.02,
+                                                                                  frictionCoef=0.5, tolerance=1.0e-4, maxIterations=20))
+
+    Logo.add(NodeModifier, on = [Logo], parameters = FixConstraintParameters( boxROIs=[[-1, -2, -13, 3, 2, -7]]))
+    Logo.add(NodeModifier, on = [Logo], parameters = FixConstraintParameters( boxROIs=[[-100, -2, -13, -300, 2, -7]]))
+    ModelsNode.add(NodeModifier, on = [S, Logo], parameters = AttachmentConstraintParameters(name = "AttachmentConstraintParameters", indices1=[26,20,119,121], indices2=[722,732,574,573], stiffness=0.5, damping=0.0,
+                                                                                             length=[((9.43-9.35)**2 + (-.44-0.48)**2 + (-6.01+6.56)**2)**(0.5) ]))
+

splib/Testing.py

@@ -4,7 +4,7 @@
 from splib.simulation.linear_solvers import *
 from splib.mechanics.linear_elasticity import *
 from splib.mechanics.mass import *
-from splib.mechanics.fix_points import *
+from splib.mechanics.attachment import *
 from splib.topology.loader import *
 from splib.core.node_wrapper import *
 

splib/mechanics/__init__.py

@@ -1 +1 @@
-__all__ = ["linear_elasticity","hyperelasticity","fix_points","collision_model","mass"]
+__all__ = ["linear_elasticity","hyperelasticity","attachment","collision_model","mass"]

splib/mechanics/fix_points.py → splib/mechanics/attachment.py

@@ -28,3 +28,21 @@ def addFixation(node,type:ConstraintType,boxROIs=DEFAULT_VALUE, sphereROIs=DEFAU
         case _:
             print('Contraint type is either ConstraintType.PROJECTIVE, ConstraintType.WEAK or ConstraintType.LAGRANGIAN')
             return
+
+
+
+@ReusableMethod
+def attachObjects(node, type:ConstraintType, object1, object2, indices1, indices2, stiffness=DEFAULT_VALUE, damping=DEFAULT_VALUE, lengths=DEFAULT_VALUE, **kwargs):
+    match type:
+        case ConstraintType.WEAK:
+            node.addObject("SpringForceField", name="attachment", object1=object1, object2=object2, indices1=indices1, indices2=indices2, stiffness=stiffness, damping=damping, lengths=lengths, **kwargs)
+            return
+        case ConstraintType.PROJECTIVE:
+            node.addObject("AttachProjectiveConstraint", name="attachment", object1=object1, object2=object2, indices1=indices1, indices2=indices2, **kwargs)
+            return
+        case ConstraintType.LAGRANGIAN:
+            node.addObject("BilateralLagrangianConstraint", name="attachment", object1=object1, object2=object2, first_point=indices1, second_point=indices2, **kwargs)
+            return
+        case _:
+            print('Contraint type is either ConstraintType.PROJECTIVE, ConstraintType.WEAK or ConstraintType.LAGRANGIAN')
+            return

splib/simulation/headers.py

@@ -67,16 +67,16 @@ def setupPenalityCollisionHeader(node,  displayFlags = "showVisualModels",backgr
     node.addObject(parallelPrefix+'BVHNarrowPhase',  name="narrowPhase", **kwargs)
 
     if(stick):
-        node.addObject('CollisionResponse',name="ContactManager", response="BarycentricStickContact",**kwargs)
+        node.addObject('CollisionResponse',name="ContactManager", response="StickContactForceField",**kwargs)
     else:
-        node.addObject('CollisionResponse',name="ContactManager", response="BarycentricPenalityContact",**kwargs)
+        node.addObject('CollisionResponse',name="ContactManager", response="PenalityContactForceField",**kwargs)
     node.addObject('LocalMinDistance' ,name="Distance", alarmDistance=alarmDistance, contactDistance=contactDistance, **kwargs)
 
     return node
 
 
 @ReusableMethod
-def setupLagrangianCollision(node,  displayFlags = "showVisualModels",backgroundColor=[1,1,1,1], parallelComputing=False, stick=False, alarmDistance=DEFAULT_VALUE, contactDistance=DEFAULT_VALUE, frictionCoef=0.0, tolerance=0.0, maxIterations=100, **kwargs):
+def setupLagrangianCollision(node,  enableCollision = True, displayFlags = "showVisualModels",backgroundColor=[1,1,1,1], parallelComputing=False, stick=False, alarmDistance=DEFAULT_VALUE, contactDistance=DEFAULT_VALUE, frictionCoef=0.0, tolerance=0.0, maxIterations=100, **kwargs):
     node.addObject('VisualStyle', displayFlags=displayFlags)
     node.addObject('BackgroundSetting', color=backgroundColor)
 
@@ -107,21 +107,24 @@ def setupLagrangianCollision(node,  displayFlags = "showVisualModels",background
     if(parallelComputing):
         parallelPrefix="Parallel"
 
-    node.addObject('CollisionPipeline', name="collisionPipeline",
-                   **kwargs)
+    if enableCollision:
+        node.addObject('CollisionPipeline', name="collisionPipeline",
+                       **kwargs)
 
-    node.addObject(parallelPrefix+'BruteForceBroadPhase', name="broadPhase",
-                   **kwargs)
+        node.addObject(parallelPrefix+'BruteForceBroadPhase', name="broadPhase",
+                       **kwargs)
 
-    node.addObject(parallelPrefix+'BVHNarrowPhase',  name="narrowPhase",
-                   **kwargs)
+        node.addObject(parallelPrefix+'BVHNarrowPhase',  name="narrowPhase",
+                       **kwargs)
+
+        if(stick):
+            node.addObject('CollisionResponse',name="ContactManager", response="StickContactConstraint", responseParams="tol="+str(tolerance),**kwargs)
+        else:
+            node.addObject('CollisionResponse',name="ContactManager", response="FrictionContactConstraint", responseParams="mu="+str(frictionCoef),**kwargs)
+
+        node.addObject('NewProximityIntersection' ,name="Distance", alarmDistance=alarmDistance, contactDistance=contactDistance, **kwargs)
 
-    if(stick):
-        node.addObject('CollisionResponse',name="ContactManager", response="StickContactConstraint", responseParams="tol="+str(tolerance),**kwargs)
-    else:
-        node.addObject('CollisionResponse',name="ContactManager", response="FrictionContactConstraint", responseParams="mu="+str(frictionCoef),**kwargs)
 
-    node.addObject('NewProximityIntersection' ,name="Distance", alarmDistance=alarmDistance, contactDistance=contactDistance, **kwargs)
     node.addObject('BlockGaussSeidelConstraintSolver',name="ConstraintSolver", tolerance=tolerance, maxIterations=maxIterations, multithreading=parallelComputing,**kwargs)
     node.addObject("ConstraintAttachButtonSetting")
 

stlib/__init__.py

@@ -1,7 +1,8 @@
-__all__ = ["core","entities","geometries","materials","collision","visual"]
+__all__ = ["core","entities","geometries","materials","collision","visual", "node_modifiers"]
 
 import Sofa.Core
 from stlib.core.basePrefab import BasePrefab
+from stlib.node_modifiers import NodeModifier
 
 def __genericAdd(self : Sofa.Core.Node, typeName, **kwargs):
     def findName(cname, names):
@@ -25,9 +26,13 @@ def checkName(context : Sofa.Core.Node, name):
 
     # Check if a name is provided, if not, use the one of the class
     params = kwargs.copy()
-    if isinstance(typeName, type) and issubclass(typeName, BasePrefab): #Only for prefabs
-        if len(params.keys()) > 1 or (len(params.keys()) == 1 and "parameters" not in params):
-            raise RuntimeError("Invalid argument, a prefab takes only the \"parameters\" kwargs as input")
+    if isinstance(typeName, type) :
+        if issubclass(typeName, BasePrefab): #Only for prefabs
+            if len(params.keys()) > 1 or (len(params.keys()) == 1 and "parameters" not in params):
+                raise RuntimeError("Invalid argument, a prefab takes only the \"parameters\" kwargs as input")
+        elif issubclass(typeName, NodeModifier):
+            if len(params.keys()) > 2 or (len(params.keys()) == 2 and (("parameters" not in params) or ("on" not in params ))):
+                raise RuntimeError("Invalid argument, a node modifier takes two parameters: the node modifier parameters structure as \"parameters\" kwarg and the \"on\" kwarg as a list of nodes on which to apply the modifier")
 
     elif "name" not in params : #This doesn't apply to prefab
         if isinstance(typeName, str):
@@ -43,7 +48,7 @@ def checkName(context : Sofa.Core.Node, name):
         else:
             raise RuntimeError("Invalid argument ", typeName)
 
-    if isinstance(typeName, type) and issubclass(typeName, BasePrefab) and len(params.keys()) == 1:
+    if isinstance(typeName, type) and ((issubclass(typeName, BasePrefab) and len(params.keys()) == 1) or (issubclass(typeName, NodeModifier) and  len(params.keys()) == 2)):
         params["parameters"].name = checkName(self, params["parameters"].name)
     else:
         params["name"] = checkName(self, params["name"])
@@ -53,6 +58,12 @@ def checkName(context : Sofa.Core.Node, name):
         pref = self.addChild(typeName(**params))
         pref.init()
         pref.postInit()
+    elif isinstance(typeName, type) and issubclass(typeName, NodeModifier):
+        if "Modifiers" not in self.children:
+            self.addChild("Modifiers")
+        pref = self.Modifiers.addObject(typeName(params["parameters"]))
+        touchedNodes = pref.apply(self, params["on"])
+        pref.register(self, touchedNodes)
     elif isinstance(typeName, Sofa.Core.Node):
         pref = self.addChild(typeName(**params))
     elif isinstance(typeName, type) and issubclass(typeName, Sofa.Core.Object):

stlib/geometries/__geometry__.py

@@ -1,5 +1,5 @@
 from stlib.core.basePrefab import BasePrefab
-from stlib.core.baseParameters import BaseParameters, Optional, dataclasses, Any, Callable
+from stlib.core.baseParameters import BaseParameters, Optional, dataclasses, Any
 from splib.topology.dynamic import addDynamicTopology
 from splib.topology.static import addStaticTopology
 from splib.core.enum_types import ElementType

stlib/geometries/extract.py

@@ -48,6 +48,5 @@ def __init__(self,
         self.postInit = partial(extractGeometry, sourceParameters.elementType)
 
         if(not (sourceParameters.elementType == ElementType.TETRAHEDRA and destinationType == ElementType.TRIANGLES)
-            and not (sourceParameters.elementType == ElementType.HEXAHEDRA and destinationType == ElementType.QUADS) ):
-            raise ValueError("Only configuration possible are 'Tetrahedra to Triangles' and 'Hexahedra to Quads'")
-
+                and not (sourceParameters.elementType == ElementType.HEXAHEDRA and destinationType == ElementType.QUADS) ):
+            raise ValueError("Only configuration possible are 'Tetrahedra to Triangles' and 'Hexahedra to Quads'")

stlib/node_modifiers/__init__.py

@@ -0,0 +1 @@
+from .__node_modifier__ import *