Home › Forum › SOFA › Using SOFA › [SOLVED] How can I construct a stable beam force field while the beams are not straight?
- This topic has 11 replies, 5 voices, and was last updated 4 years, 4 months ago by
MarcoMag.
-
AuthorPosts
-
5 October 2018 at 05:44 #12098
Wong
BlockedHello everyone,
I construct a scene using the component BeamFEMForceField where the beams are not straight. But the scene is not stable when it runs.
So how to solve the problem?
<Node name="root" dt="0.005"> <VisualStyle displayFlags="showBehaviorModels showForceFields showCollisionModels" /> <CollisionPipeline depth="6" verbose="0" draw="0" /> <BruteForceDetection name="N2" /> <CollisionResponse name="response" response="FrictionContact" /> <LocalMinDistance name="proximity" alarmDistance="0.15" contactDistance="0.05" angleCone="0.0" /> <FreeMotionAnimationLoop/> <LCPConstraintSolver tolerance="0.001" maxIt="1000"/> <!--<CollisionGroup name="Group" />--> <Node name="beam"> <EulerImplicit rayleighStiffness="0" printLog="false" rayleighMass="0.1" /> <CGLinearSolver iterations="100" tolerance="1e-5" threshold="1e-5"/> <MechanicalObject template="Rigid" name="DOFs" position="0 0 0 0 0 0 1 3.5 0 0 0 0 0 1 3.5 0 3.5 0 -0.707 0 0.707" /> <Mesh name="lines" lines="0 1 1 2" /> <FixedConstraint name="FixedConstraint" indices="0" /> <UniformMass totalMass = "1"/> <BeamFEMForceField name="FEM" radius="0.1" youngModulus="20000000" poissonRatio="0.49"/> <UncoupledConstraintCorrection/> </Node> </Node>
11 October 2018 at 19:46 #12168Christiane
BlockedHello Wong,
I had the same issue; I managed to stabilize it by setting the orientations of the DOFs in the mechanical object. Basically, you need to orient each node so that the x-axis is the slope of the curved beam.
I did this by creating a spline of my curved beams and then calculating the immediate slope at each node, aligning the x-axis with the spline. I didn’t bother with the y and z-axis, just oriented them randomly, because I was using beams with round sections anyway. I compared the resulting SOFA beams with an equivalent Ansys beam, and I had good agreement between the two.
Regards,
Christiane26 October 2018 at 18:30 #1227929 October 2018 at 02:48 #12281Wong
BlockedHello @hugo,
I think it should be a solution because I have tried it as well before. But now I do not use the beams in my scene any more.
By the way, can you answer my another question at this link?
13 November 2018 at 11:02 #12364shang
Blockedhello,@Christiane
I am wondering how you created the spline corresponding to the FEM beams, is the spline created by SOFA components? I am coming to the same issue with yours, I am appreciate that you can share your idea. Thank you!14 November 2018 at 08:52 #12432Christiane
BlockedHello Xiaojuan,
I used Matlab to define my spline and beam orientations. You can use the spline function with your beam points; it will help you define a continuous point series that represents your curved beam. Then you just use these points to determine the instantaneous slope at each point of your beam: that is the orientation you give to SOFA for your beam points.
Voilà!
Christiane14 November 2018 at 10:58 #12436shang
BlockedHi,Christiane
Does “beam points” the mesh for the beam? Does the “instantaneous slope” means the value of x-axis value of the position for the MechanicalObject in “beam” node? and What is your input of the spline function in matlab? is the x-axis value of the position of beam MechanicalObject? I am a little bit confused, thank you for helping me with this issue.Here is my definition of my beam:
<Node name="beam" gravity="0 9.8 0"> <EulerImplicit rayleighStiffness="0.1" printLog="false" /> <CGLinearSolver template="GraphScattered" iterations="20" threshold="1e-008" tolerance="1e-5" /> <MechanicalObject template="Rigid" name="catheterDOFs" rotation="0 -90 0" position=" 0 0 0 0 0 0 1 0.5 0 0 0 0 0 1 1 0 0 0 0 0 1 1.5 0 0 0 0 0 1 2 0 0 0 0 0 1 2.5 0 0 0 0 0 1 3 0 0 0 0 0 1 3.5 0 0 0 0 0 1 4 0 0 0 0 0 1 4.5 0 0 0 0 0 1 5 0 0 0 0 0 1 5.5 0 0 0 0 0 1 6 0 0 0 0 0 1 6.5 0 0 0 0 0 1 7 0 0 0 0 0 1 7.5 0 0 0 0 0 1 8 0 0 0 0 0 1 8.5 0 0 0 0 0 1 9 0 0 0 0 0 1 9.5 0 0 0 0 0 1 10 0 0 0 0 0 1 10.5 0 0 0 0 0 1 11 0 0 0 0 0 1 11.5 0 0 0 0 0 1 12 0 0 0 0 0 1 12.5 0 0 0 0 0 1 13 0 0 0 0 0 1 13.5 0 0 0 0 0 1 14 0 0 0 0 0 1 14.5 0 0 0 0 0 1 15 0 0 0 0 0 1 15.5 0 0 0 0 0 1 16 0 0 0 0 0 1 16.5 0 0 0 0 0 1 17 0 0 0 0 0 1 17.5 0 0 0 0 0 1 18 0 0 0 0 0 1 18.5 0 0 0 0 0 1 19 0 0 0 0 0 1 19.5 0 0 0 0 0 1 20 0 0 0 0 0 1 20.5 0 0 0 0 0 1 21 0 0 0 0 0 1 21.5 0 0 0 0 0 1 22 0 0 0 0 0 1 22.5 0 0 0 0 0 1 23 0 0 0 0 0 1 23.5 0 0 0 0 0 1 24 0 0 0 0 0 1 24.5 0 0 0 0 0 1 25 0 0 0 0 0 1 25.5 0 0 0 0 0 1 26 0 0 0 0 0 1 26.5 0 0 0 0 0 1 27 0 0 0 0 0 1 27.5 0 0 0 0 0 1 28 0 0 0 0 0 1 28.5 0 0 0 0 0 1 29 0 0 0 0 0 1 29.5 0 0 0 0 0 1 " /> <Mesh name="lines" lines=" 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31 32 32 33 33 34 34 35 35 36 36 37 37 38 38 39 39 40 40 41 41 42 42 43 43 44 44 45 45 46 46 47 47 48 48 49 49 50 50 51 51 52 52 53 53 54 54 55 55 56 56 57 57 58 58 59" /> <Monitor name="catheterInfo" template="Rigid" listening="1" indices="0 1 2 3 4" showPositions="0" PositionsColor="1 1 0 1" showMinThreshold="0.01" TrajectoriesPrecision="0.1" TrajectoriesColor="1 1 0 1" sizeFactor="1" /> <FixedConstraint name="FixedConstraint" indices="0" /> <UniformMass totalmass="6" /> <BeamFEMForceField name="FEM" radius="10" youngModulus="2000000" poissonRatio="0.9"/> <UncoupledConstraintCorrection /> <ConstantForceField /> <Node name="Collision"> <CylinderGridTopology name="coli" nx="6" ny="6" nz="59" length="59" radius="0.3" axis="1 0 0" /> <MechanicalObject /> <BeamLinearMapping isMechanical="true" /> <Triangle /> <Line /> <Point /> </Node> <Node name="visu"> <CylinderGridTopology name="coli" nx="6" ny="6" nz="59" length="59" radius="0.3" axis="1 0 0" /> <OglModel color="gray" /> <BeamLinearMapping isMechanical="true" /> </Node> <PythonScriptController name="catheter" filename="ScriptEvents.py" classname="catheter"/> </Node>
14 November 2018 at 11:22 #12437Christiane
BlockedHello again,
The beam points: I had 20 points that described my curved beams. Those are the points I used to mesh the beam, and for which I needed an orientation.
The Matlab spline function: I had the x, y and z coordinates of each point, so what I had was really 3 parametric splines, one for x, for y and for z, with t varying from 1 to 20. I simply input each spline into the spline function with t = 0:0.01:21. That gave me xx, yy and zz with a lot more points than my 20 initial points, and with extrapolation at both ends so I could determine the slope of those points as well.
With those points series, it was then very easy to calculate the 3D slope at my beam points. I then just transformed that into quaternions, and that is was I gave to SOFA as orientations. And yes, those orientations represent the x-axis of your beam.
Regards,
Christiane14 November 2018 at 12:37 #12438shang
BlockedDear Christiane
I have tried spline function with input:x-axis y-axis z-axis value of the initial position of my beams. And I calculate them in matlab. Because the x-axis value is linear as you can see in the scene file, and the y and z-axis value are all zeros. The spline interpolation was only done on x-values, and the spline result xx is same as x. The matlab code is list below:x = load('xx.txt')'; len=length(x); t=1:len; p=spline(t,x) tt = 0:0.1:len; xx=ppval(p,tt); plot(t,x) hold on plot(tt,xx)
and I also want to know whether the definition of “3D slope” that you mentioned is calculated with equation: sqrt(x^2 + y^2 + z^2).The result of this equation is a scalar.How to transform this value to a quaternion? And where to set the orientations you mentioned,which represent the x-axis of beams? Could the quaternions be set in “position” tag of the MechanicalObject of our beam mesh?
14 November 2018 at 13:00 #12441Christiane
BlockedI’m might be using totally the wrong English term, the slope is the first derivative of the spline. So if I want the local value of the slope at one point, I just calculate the vector between the two adjacent points. Because I used the Matlab spline function to create a series of closely set points, I can have a fairly localized estimation of the orientation. That worked well for me.
On the other hand, your current beam is set along the x-axis, so the orientation that you gave is correct.
Regards,
Christiane14 November 2018 at 14:46 #12443shang
BlockedHi,Christiane
You mean that my definition of mesh and MechanicalObject of the beams is totally right? But the catheter that I simulated using this node above is not stable at all, I am confused of this.PS:I am wondering whether your definition for the orientation of each point is the first 4 value of the position(which is type of Rigid)?
4 May 2020 at 17:28 #16012MarcoMag
BlockedHello everyone,
I know this is an old post but I came across similar problems recently.
Below is an example of how I orient my nodes.
Note that I make use of the python packagesnumpy
andpyquaternion
that can be installed viapip
.import sys import Sofa import numpy as np from pyquaternion import Quaternion def getCoordPointsAlongHellix(r, h, npt, tmax): # position vectors along the curve x = np.zeros((npt,3), dtype = float) # tangent vectors xp = np.zeros((npt,3), dtype = float) t = np.linspace(0,tmax, npt) x[:,0] = r * np.cos(t) x[:,1] = r * np.sin(t) x[:,2] = h * t xp[:,0] = -r * np.sin(t) xp[:,1] = r * np.cos(t) xp[:,2] = h # quaternions q = np.zeros((npt,4), dtype = float) a = np.array([1,0,0]) # trick: find the rotation matrix rotating global basis vector a into tangent vector b # extract the associated quaternion for i in range(npt): #tangent to the curve at the node b = xp[i]/np.linalg.norm(xp[i]) assert np.allclose(np.linalg.norm(b), 1.) v = np.cross(a,b) s = np.linalg.norm(v) if not np.allclose(s, 0): c = np.dot(a,b) vskew = np.array([[0, -v[2], v[1]], [v[2], 0 , -v[0]], [-v[1], v[0], 0] ]) # rotation matrix rotating a into b R = np.eye(3) + vskew + np.dot(vskew, vskew) * ((1-c)/(s**2)) else: R = np.eye(3) assert np.allclose(R.dot(a), b) # extract the quaternion qi = Quaternion(matrix=R).elements # components provided by pyquaternion are not in the same order # as the one expected by the constructor of Quater in SOFA # we must reorder them q[i]= [qi[1], qi[2], qi[3], qi[0]] return x, q class BeamFEMForceField (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 rootNode.createObject('RequiredPlugin', name='SofaOpenglVisual') rootNode.createObject('RequiredPlugin', name='SofaPython') showCollisionModels = 0 if showCollisionModels: rootNode.createObject('VisualStyle', displayFlags='showBehaviorModels showForceFields showCollisionModels') else: rootNode.createObject('VisualStyle', displayFlags='showBehaviorModels showForceFields') rootNode.createObject('DefaultPipeline', draw='0', depth='6', verbose='0') rootNode.createObject('BruteForceDetection', name='N2') rootNode.createObject('MinProximityIntersection', contactDistance='0.02', alarmDistance='0.03', name='Proximity') rootNode.createObject('DefaultContactManager', name='Response', response='default') # radius of the beam radius = 0.1 radiusSphereCollisionModel = str(1. * radius) # heigth of the strand h = 2 # rootNode/beamI --> straight beam nn = 5 nbeam = nn - 1 Coord = np.zeros((nn, 7),dtype=float) Coord[:,0] = 0 Coord[:,1] = 0 Coord[:,2] = np.linspace(0,h,nn) lines = np.zeros((nbeam, 2), dtype=int) lines[:,0] = np.arange(0, nn -1) lines[:,1] = np.arange(1, nn) # topology lines = str(lines.flatten()).replace('[', '').replace(']','') x, q = getCoordPointsAlongHellix(0, h, nn, tmax = 0.5 * np.pi) Coord[:, 3:] = q # import pdb; pdb.set_trace() strCoord = str(Coord.flatten()).replace('\n', '').replace('[', '').replace(']','') beamI = rootNode.createChild('beamI') self.beamI = beamI beamI.createObject('EulerImplicitSolver', rayleighStiffness='0', printLog='false', rayleighMass='0.1') beamI.createObject('BTDLinearSolver', printLog='false', template='BTDMatrix6d', verbose='false') beamI.createObject('MechanicalObject', position=strCoord, name='DOFs', template='Rigid3d') beamI.createObject('MeshTopology', lines=lines, name='lines') beamI.createObject('FixedConstraint', indices='0', name='FixedConstraint') beamI.createObject('UniformMass', vertexMass='1 1 0.01 0 0 0 0.1 0 0 0 0.1', printLog='false') beamI.createObject('BeamFEMForceField', radius=str(radius), name='FEM', poissonRatio='0.49', youngModulus='20000000') Collision = beamI.createChild('Collision') self.Collision = Collision beamI.createObject('SphereModel', radius=radiusSphereCollisionModel, name='SphereCollision') # rootNode/beam2 --> helicoidal beam nn = 10 nbeam = nn - 1 lines = np.zeros((nbeam, 2), dtype=int) lines[:,0] = np.arange(0, nn -1) lines[:,1] = np.arange(1, nn) # topology lines = str(lines.flatten()).replace('[', '').replace(']','') Coord = np.zeros((nn, 7),dtype=float) x, q = getCoordPointsAlongHellix(2 * radius, h, nn, tmax = 0.5 * np.pi) Coord[:,:3] = x Coord[:, 3:] = q strCoord = str(Coord.flatten()).replace('\n', '').replace('[', '').replace(']','') beam2 = rootNode.createChild('beam2') self.beam2 = beam2 beam2.createObject('EulerImplicitSolver', rayleighStiffness='0', printLog='false', rayleighMass='0.1') beam2.createObject('BTDLinearSolver', printLog='false', template='BTDMatrix6d', verbose='false') # beam2.createObject('MechanicalObject', position='0 -2 0.25 0 0 0 1 0 -1 0.25 0 0 0 1 0 0 0.25 0 0 0 1 0 1 0.25 0 0 0 1 0 2 0.25 0 0 0 1', name='DOFs', template='Rigid') beam2.createObject('MechanicalObject', position=strCoord, name='DOFs', template='Rigid3d') beam2.createObject('MeshTopology', lines=lines, name='lines') beam2.createObject('FixedConstraint', indices='0', name='FixedConstraint') beam2.createObject('UniformMass', vertexMass='1 1 0.01 0 0 0 0.1 0 0 0 0.1', printLog='false') beam2.createObject('BeamFEMForceField', radius=str(radius), name='FEM', poissonRatio='0.49', youngModulus='20000000') beam2.createObject('SphereModel', radius=radiusSphereCollisionModel, name='SphereCollision') 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 /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def onKeyPressed(self, c): ## usage e.g. #if c=="A" : # print "You pressed control+a" return 0; def onMouseWheel(self, mouseX,mouseY,wheelDelta): ## usage e.g. #if isPressed : # print "Control button pressed+mouse wheel turned at position "+str(mouseX)+", "+str(mouseY)+", wheel delta"+str(wheelDelta) return 0; def storeResetState(self): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def cleanup(self): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def onGUIEvent(self, strControlID,valueName,strValue): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def onEndAnimationStep(self, deltaTime): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def onLoaded(self, node): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def reset(self): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def onMouseButtonMiddle(self, mouseX,mouseY,isPressed): ## usage e.g. #if isPressed : # print "Control+Middle mouse button pressed at position "+str(mouseX)+", "+str(mouseY) return 0; def bwdInitGraph(self, node): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def onScriptEvent(self, senderNode, eventName,data): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def onMouseButtonRight(self, mouseX,mouseY,isPressed): ## usage e.g. #if isPressed : # print "Control+Right mouse button pressed at position "+str(mouseX)+", "+str(mouseY) return 0; def onBeginAnimationStep(self, deltaTime): ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py return 0; def createScene(rootNode): rootNode.findData('dt').value = '0.01' rootNode.findData('gravity').value = '0 0 -9.81' try : sys.argv[0] except : commandLineArguments = [] else : commandLineArguments = sys.argv myBeamFEMForceField = BeamFEMForceField(rootNode,commandLineArguments) return 0;
Hope this helps.
-
AuthorPosts
- You must be logged in to reply to this topic.