[Hugo]

Dear Catherine,

Please have a look at :

• the introduction paper
• the documentation (Main principles are described)
• and the tutorial scenes that you can start with (in SOFA, examples/Tutorials/ folder.

Unfortunately SOFA is not yet super newcomer-friendly but this is our 2018 mission. Let me know if you would be interest in a specific training session for SOFA.

Best regards,

[Hugo]

Dear Catherine,

Welcome on the SOFA forum!
Yes it is possible. Physics simulation is what SOFA is made for.
Do you have already a SOFA scene ?

Hugo

[Hugo]

What’s your collision method and model exactly ? Or simpler : could you share with us the structure of your scene please ?

externalForce is a vector summing up forces coming from external interactions (usually fully explicit forces). Typically, forces due to the mouse interactions, JointSpringForceField or collision. With collision using an UncoupledConstraintCorrection for instance, the externalForce is filled when applying the contact force (as you can see UncoupledConstraintCorrection.inl line 484). But at the end of each time step, the externalForce vector is reset by the MechanicalObject (see MechanicalObject.inl line 1465)

Hope this help,

[Hugo]

Yes but I was talking about the “force” vector of the fixed mechanical object.
What about this one?

Hugo

[Hugo]

Dear @r0m1

Welcome on the SOFA forum!
Since the second object is fixed, all forces applied on this object should result from the collision only. Therefore, one way would be to access the “force” vector of this mechanical object in your python script and trigger any action regarding your threshold.

Let me know if this helps.

Hugo

[Hugo]

Dear @sdauto

the JointSpringForceField inherits from the PairInteractionForceField. This specific class of forcefield links two mechanical objects and apply forces on them. At the creation of the JointSpringForceField, the two mechanical objects are searched in the current node (where the JointSpringForceField is set) but the path to these mechanical objects can be manually set using “object1” and “object2“.

Then, to create springs between these objects all you need to specific is JointSpring (as described in JointSpringForceField.h line 83). For instance in XML format, you can write :

spring="BEGIN_SPRING  0 1  FREE_AXIS 0 0 0 0 1 0  KS_T 0 30000  KS_R 0 200000  KD 1  R_LIM_X -0.8 0.8  R_LIM_Y -1.57 1.57  R_LIM_Z 0 0  END_SPRING"

You can also have a look at the example in examples/Components/forcefield/JointSpringForceField.scn
Let me know if this helps.

Hugo

[Hugo]

Great @ant0nin
Thank you for your feedback good luck for the future developments and do not hesitate to ask any further question or help new users 😉

Hugo

[Hugo]

Now that you activated the plugin, you can try some example scenes.
Have you tried to start one of the scenes?

Make sure you have a RequiredPlugin in the scene:

<RequiredPlugin name='Xitact' />

And load the right driver in your scene.

Best,

Hugo

[Hugo]

Dear @sarath ,

It seems you are facing the same issue that this post, but I am not sure of the consequences of the active option: NO_BOOST_CMAKE.
Are you on CentOS ? it seems linked to this distribution.

Your error in SOFA comes from FileInfo.h but I guess it has a link with the above mentioned Cmake error:

CMake Error at /usr/lib64/boost/Boost.cmake:536 (message):
  The imported target "boost_date_time-static" references the file

     "/usr/lib64/lib64/libboost_date_time.a"

  but this file does not exist.

Let us know about your investigation.

Hugo

[Hugo]

Ah damed it, I just noticed as weel that the plugin is deprecated in the latest version of SOFA. The error seems to come from the code : VRPNImager.cpp: line 34

I will try to find some time to investigate, if anyone else got an input it would be kind.

Cheers

Hugo

[Hugo]

hey @ant0nin ,

Have you looked at the SlidingConstraint ? There is a scene example in examples/Components/constraint/SlidingConstraint.scn.

Cheers,

Hugo

[Hugo]

Dear @sdauto

I understand that you build your scene in C++.
But can you explain me why you need the addForce() function ?

This function is usually called by the ODE solver (integration scheme) to build the vector b of the Ax = b system.

Hugo

[Hugo]

Dear @sarath ,

First of all, I would advise to use a more recent version of SOFA (check out the latest 17.06 release). I did not use the VRPN plugin myself though, not sure to be really helpful.

Maybe anyone else in the community could help ?
I will activate and test the plugin.

Keep us updated. Cheers,

Hugo

[Hugo]

Dear @secretdevil,

Please find below an example with Python (Python plugin must be compiled) to dynamically create grids.

import sys
import Sofa

class RegularGridTopology (Sofa.PythonScriptController):

    def __init__(self, node, commandLineArguments) : 
        self.commandLineArguments = commandLineArguments
        print "Command line arguments for python : "+str(commandLineArguments)
        self.createGraph(node)
        return None;

    def createGraph(self,rootNode):

        # rootNode
        self.rootNode = rootNode
        rootNode.createObject('VisualStyle', displayFlags='showBehaviorModels showForceFields showVisual')
        rootNode.createObject('CollisionPipeline', verbose='0')
        rootNode.createObject('BruteForceDetection', name='N2')
        rootNode.createObject('CollisionResponse', response='default')
        rootNode.createObject('DiscreteIntersection')

        # rootNode/LiverFFD-lowres
        LiverFFD_lowres = rootNode.createChild('LiverFFD-lowres')
        self.LiverFFD_lowres = LiverFFD_lowres
        LiverFFD_lowres.createObject('EulerImplicit', rayleighStiffness='0.1', rayleighMass='0.1')
        LiverFFD_lowres.createObject('CGLinearSolver', threshold='1e-7', tolerance='1e-7', iterations='100')
        LiverFFD_lowres.createObject('MechanicalObject')
        LiverFFD_lowres.createObject('UniformMass', totalmass='100.0')
        LiverFFD_lowres.createObject('RegularGrid', nx='4', ny='3', nz='3', xmin='-10.25', ymin='0.25', zmin='-2', xmax='-3.25', ymax='5.25', zmax='3')
        LiverFFD_lowres.createObject('BoxROI', name='myboxROI', box='-9 0 0 -4 0.9 5', drawBoxes='1') # xmin,ymin,zmin, xmax,ymax,zmax
        LiverFFD_lowres.createObject('FixedConstraint', indices='@myboxROI.indices')
        LiverFFD_lowres.createObject('RegularGridSpringForceField', name='Springs', stiffness='4000', damping='4')
        # rootNode/LiverFFD-lowres/Visu
        Visu = LiverFFD_lowres.createChild('Visu')
        self.Visu = Visu
        Visu.createObject('OglModel', color='red', translation='-5 0 0', name='Visual', fileMesh='mesh/liver-smooth.obj')
        Visu.createObject('BarycentricMapping', input='@..', output='@Visual')

        #Compute the iterations
        self.iteration_number=0

        return 0;

    def onMouseButtonLeft(self, mouseX,mouseY,isPressed):
        ## usage e.g.
        #if isPressed : 
        #    print "Control+Left mouse button pressed at position "+str(mouseX)+", "+str(mouseY)
        return 0;

    def onKeyReleased(self, c):
        ## usage e.g.
        #if c=="A" :
        #    print "You released a"
        return 0;

    def initGraph(self, node):
        ## Please feel free to add an example for a simple usage in /data/Softwares/sofa/src/master/applications/plugins/SofaPython/scn2python.py
        return 0;

    def onKeyPressed(self, c):
        ## usage e.g.
        if c=="A" :
            self.iteration_number=self.iteration_number+1
            print "iteration "+str(self.iteration_number)

            yPositionMin = 0.25+self.iteration_number*7
            yPositionMax = yPositionMin+5.25

            yPositionROIMin = self.iteration_number*7
            yPositionROIMax = yPositionROIMin+0.9

            myGrid = self.rootNode.createChild('newGrid-'+str(self.iteration_number))

            newODE = myGrid.createObject('EulerImplicit', rayleighStiffness='0.1', rayleighMass='0.1')
            newLSolver = myGrid.createObject('CGLinearSolver', threshold='1e-7', tolerance='1e-7', iterations='100')
            newMO = myGrid.createObject('MechanicalObject', showObject='1')
            newMass = myGrid.createObject('UniformMass', totalmass='100.0')
            myRegularGrid = myGrid.createObject('RegularGrid', name='myRegularGrid', nx='4', ny='3', nz='3', xmax='-3.25', xmin='-10.25', ymin=yPositionMin, ymax=yPositionMax, zmin='-2', zmax='3')

            newBox = myGrid.createObject('BoxROI', name='myboxROI', box='-9 '+str(yPositionROIMin)+' 0 -4 '+str(yPositionROIMax)+' 5', drawBoxes='1')
            newConstraint = myGrid.createObject('FixedConstraint', indices='@myboxROI.indices')
            newFF = myGrid.createObject('RegularGridSpringForceField', name='Springs', stiffness='4000', damping='4', drawMode="1")

            newODE.init()
            newLSolver.init()
            newMO.init()
            newMass.init()
            myRegularGrid.init()

            newBox.init()
            newConstraint.init()
            newFF.init()

            # rootNode/LiverFFD-lowres/Visu
            Visu = myGrid.createChild('Visu')
            self.Visu = Visu
            newOGL = Visu.createObject('OglModel', color='red', translation='-5 '+str(yPositionROIMin)+' 0', name='Visual', fileMesh='mesh/liver-smooth.obj')
            newMapping = Visu.createObject('BarycentricMapping', input='@..', output='@Visual')
            newOGL.init()
            newMapping.init()

        return 0;

def createScene(rootNode):
    rootNode.findData('dt').value = '0.01'

    rootNode.createObject('RequiredPlugin', name='SofaPython')

    try : 
        sys.argv[0]
    except :
        commandLineArguments = []
    else :
        commandLineArguments = sys.argv
    myRegularGridTopology = RegularGridTopology(rootNode,commandLineArguments)
    return 0;

You can run this script using runSofa :

runSofa your_script.py

Let me know if it helps

Hugo

[Hugo]

Dear @sdauto,

There might be a misunderstanding of the physics resolution in SOFA. The mass-spring model defines a mass and a stiffness (through springs on each edge of the mesh) for the system. This is the physics. Then there is the resolution part where you need solvers to solve your linear system Ax=b

I invite you to take a look at the documentation for system resolution.

1-addForce() is a function called to fill the vector b of your linear system Ax=b
2-the addForce() is a function already implemented in the JointSpringForceField class. But again, this is not the way to run a simulation, you need an ODE solver and a linear solver. Check the doc on ODE solver and linear solver.

Let us know if it helps.

Hugo

[Hugo]

Dear @secretdevil,

I did not forget you, I am just gathering info to answer you and give you an example scene. Cheers

Hugo

[Hugo]

Dear @meekeech,

Unfortunately, the PLUGIN_XITACT is not much used currently at my best knowledge. It might need a bit of upgrade. I would be glad to help you with it.

The CMake error is normal since the file XiTrocarInterface.h is not in the plugin. You can try to remove the line 9 of the CMakeList.

Cheers,

Hugo

[Hugo]

Hey Xiaojuan,

I checked the file “mesh/liver2.msh” it looks pretty strange (look ids are not starting from 1). I wouldn’t use it.

In your last post (#10128), the first scene is working fine for me.

Best,

Hugo

[Hugo]

Hi Xiaojuan,

No worries, do not hesitate to set this topic as solved as soon as it works for you!
Cheers,

Hugo

[Hugo]

Hi Wendy,

Sorry for the delay 😉 but I did not forgot you, no worries!

1) did you solve your initial problem, the crash above mentioned

2) about your other question, no functions in classes cannot be called in XML. Using Python, functions in classes can be bound (slight difficult) but it could be feasible (but I am no expert).

Cheers,

Hugo

[Hugo]

Indeed, using a fine grid for Hexas and boundary conditions defined on it would work.
Let us know about the outcome of your work !

Cheers,

H

[Hugo]

Hey Wendy,

Well you can send it (the scene but mainly the associated files/meshes) through wetransfer (as mentioned above) or send it directly to the consortium email (available here, down the page).

Cheers,

Hugo

[Hugo]

Dear Xiaojuan,

First note that in your scene, tolerance and maxIt defined in EulerImplicitSolver are useless.

Moreover, your problem comes from the initialization of the MechanicalObject “independentParticles_dof”. Look at the position field, the 4th rigid particle is wrong defined (the quaternion is missing).

Here is a corrected version of your scene:

<Node    name="root"  gravity="0 -10 0"  dt="0.01"  time="0"  animate="0"  >
    <VisualStyle displayFlags="hideVisual showBehavior hideCollisionModels showMapping" />
    <DefaultAnimationLoop name="defaultAnimationLoop" />
    <DefaultVisualManagerLoop name="defaultVisualManagerLoop" />    
    
    <Node    name="simulatedScene"  gravity="0 -10 0"  dt="0.01"  time="0"  animate="0"  >
        <EulerImplicitSolver />
        <CGLinearSolver name="linear solver" iterations="25" tolerance="1e-09" threshold="1e-09" />
        
        <Node    name="rigidNode"  gravity="0 -10 0"  dt="0.01"  time="0"  animate="0"  >
            <MechanicalObject template="Rigid3d" name="rigidNode_dof"  position="0 0 0 0 0 0 1" /> 
            <FixedConstraint template="Rigid3d" name="rigidNode_fixedConstraint" indices="0"  />
            <Node    name="mappedParticles"  gravity="0 -10 0"  dt="0.01">
                <MechanicalObject template="Rigid3d" name="mappedParticles_dof"  position="0 0 0 0 0 0 1 " />
                <FixedConstraint template="Rigid3d" name="rigidNode_fixedConstraint" indices="0"  />
                <RigidRigidMapping template="Rigid3d,Rigid3d" name="mappedParticles_mapping"   input="@/simulatedScene/rigidNode/rigidNode_dof"  output="@mappedParticles_dof" />
            </Node>
        </Node>
        
        <Node    name="independentParticles"  gravity="0 -10 0"  dt="0.01"  time="0"  animate="0"  >
            <MechanicalObject template="Rigid3d" name="independentParticles_dof"  position="1 0 0 0 0 0 1 2 0 0 0 0 0 1 3 0 0 0 0 0 1 4 0 0 0 0 0 1 5 0 0 0 0 0 1 "   restScale="1" />
            <Node    name="deformableGrid"  gravity="0 -10 0"  dt="0.01"  time="0"  animate="0"  >
                <MechanicalObject template="Rigid3d" name="deformableGrid_dof"  position="0 0 0 0 0 0 1 1 0 0 0 0 0 1 2 0 0 0 0 0 1 3 0 0 0 0 0 1 4 0 0 0 0 0 1 5 0 0 0 0 0 1" />
                <Mesh name="lines" lines="0 1 1 2 2 3 3 4 4 5" />
                <SubsetMultiMapping template="Rigid3d,Rigid3d" name="deformableGrid_mapping"  indexPairs="1 0 0 0 0 1 0 2 0 3 0 4"  input="@../independentParticles_dof @../../rigidNode/mappedParticles/mappedParticles_dof"  output="@deformableGrid_dof" />
                
                <BeamFEMForceField name="deformableBeam"    poissonRatio="0.7"  youngModulus="10000000000"  radius="0.07" />
            </Node>
        </Node>
    </Node>
</Node>

Let me know if this helps (it should 😉 )

Hugo

[Hugo]

Hi @mmabrouk

I am sorry, I am not rock-solid about Rigid simulations and about the Compliant plugin (which has a specific API).

1) Yes they are doing the same, just the Compliant case is assembling (building) the matrices (even for mapping)

2) Looking from the code it seems that the format is : BEGIN_SPRING index_frame1 mass_frame1 index_frame2 mass_frame2 ENDSPRING (and additional info can be added about stiffness etc ex: KS_R 2000000 2000000 KD 1 RATIO 2.166667)

What you mention in your last post makes me think about a potential hint. The CompliantImplicitSolver implements a Newmark-beta implicit solver for constrained dynamics using assembly. Which means the system matrix A (from Ax=b system) will be built. With such direct approach, the functions that are called in the ForceField API are:

• addForce: for explicit forces stored in b vector
• addKToMatrix: for implicit forces stored in A

The CGLinearSolver is an iterative solver, in such cases, the function addDForce is called instead of addKToMatrix. This is why it works in one case and not the other one, I would guess. A solution would be to implement the same addKToMatrix in GearSpringForceField. I really hope this helps. Cheers,

Hugo

[Hugo]

Exactly @mmabrouk is 100% right, consider Python.
You’ll find examples in the SofaPython plugin, documentation. In the forum, I recently gave an example about Python here.

HTH,

Hugo.

[Author]

Posts