- This topic has 11 replies, 5 voices, and was last updated 4 years, 6 months ago by
.
-
Posts
-
Hello 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>Hello 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,
ChristianeHello @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?
hello,@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!Hello 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à!
ChristianeHi,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>Hello 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,
ChristianeDear 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?
I’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,
ChristianeHi,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)?
Hello 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 packagesnumpyandpyquaternionthat 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.
- You must be logged in to reply to this topic.
