Merge pull request #28 from SofaDefrost/add_new_models

Add new models

CMakeLists.txt

@@ -24,10 +24,18 @@ set(HEADER_FILES
     ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscoelastic/SLSKelvinVoigtFirstOrder.h
     ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscoelastic/SLSMaxwellFirstOrder.h
     ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscoelastic/Burgers.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscoelastic/LinearElastic.h
 
-
+    ## VISCOHYPERELASTIC MATERIALS 
     ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/BaseViscoHyperelasticMaterial.h
     ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSNeoHookeanFirstOrder.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSNeoHookeanSecondOrder.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSStableNeoHookeanFirstOrder.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSStableNeoHookeanSecondOrder.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSMooneyRivlinFirstOrder.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSMooneyRivlinSecondOrder.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSOgdenFirstOrder.h
+    ${SOFAVISCOELASTIC_SOURCE_DIR}/material/viscohyperelastic/SLSOgdenSecondOrder.h
 
 
 )

examples/Pendulum1.py

@@ -0,0 +1,182 @@
+# to be able to add sofa objects you need to first load the plugins that implement them.
+# For simplicity you can load the plugin "SofaComponentAll" that will load all most
+# common sofa objects.
+import SofaRuntime
+# SofaRuntime.importPlugin("SofaComponentAll")
+
+# to add elements like Node or objects
+import Sofa.Core
+root = Sofa.Core.Node()
+
+import math 
+import numpy as np
+import matplotlib.pyplot as plt
+
+import os
+path = os.path.dirname(os.path.abspath(__file__))+'/plot/'
+
+
+class PendulumController(Sofa.Core.Controller):
+
+	def __init__(self, *args, **kwargs):
+		Sofa.Core.Controller.__init__(self,*args, **kwargs) #needed
+		self.node = kwargs['node']		
+		# self.Positions =  kwargs['position']
+		# file3 = open("/data/Softwares/sofa/src/master/forum/Pasquale/plot/EulerExplicit.txt","w")
+		# file3.write(str(0.0)+' '+str(self.Positions[1][0])+' '+str(self.Positions[1][1])+' '+ str(self.Positions[1][2])+ ' '+'\n' )
+		# file3.close()
+		self.times = []
+		self.resultEE = []
+		self.resultRK = []
+		self.resultHHT = []
+		self.resultNewmark = []
+		self.resultEI = []
+		self.resultTRAP = []
+
+
+		file1 = open("/home/pasquale/EulerExplicit.txt","w")
+		file1.write(str(0.0)+' '+str(self.node.EulerExplicit.Particles.position.value[1][2])+'\n')
+		file1.close()
+
+		file2 = open("/home/pasquale/RungeKutta.txt","w")
+		file2.write(str(0.0)+' '+str(self.node.RungeKutta.Particles.position.value[1][2])+'\n')
+		file2.close()
+
+
+		file3 = open("/home/pasquale/HHT.txt","w")
+		file3.write(str(0.0)+' '+str(self.node.HHT.Particles.position.value[1][2])+'\n')
+		file3.close()
+
+		file4 = open("/home/pasquale/Newmark.txt","w")
+		file4.write(str(0.0)+' '+str(self.node.Newmark.Particles.position.value[1][2])+'\n')
+		file4.close()
+
+		file5 = open("/home/pasquale/EulerImplicit.txt","w")
+		file5.write(str(0.0)+' '+str(self.node.EulerImplicit.Particles.position.value[1][2])+'\n')
+		file5.close()
+
+		file5 = open("/home/pasquale/Trapezoidal.txt","w")
+		file5.write(str(0.0)+' '+str(self.node.Trapezoidal.Particles.position.value[1][2])+'\n')
+		file5.close()
+
+
+	def onAnimateBeginEvent(self, event): # called at each begin of animation step
+
+		self.time = self.node.time.value
+
+		file1 = open("/home/pasquale/EulerExplicit.txt","a")
+		file1.write(str(self.time)+' '+str(self.node.EulerExplicit.Particles.position.value[1][2])+'\n')
+		file1.close()
+
+		file2 = open("/home/pasquale/RungeKutta.txt","a")
+		file2.write(str(self.time)+' '+str(self.node.RungeKutta.Particles.position.value[1][2])+'\n')
+		file2.close()
+
+
+		file3 = open("/home/pasquale/HHT.txt","a")
+		file3.write(str(self.time)+' '+str(self.node.HHT.Particles.position.value[1][2])+'\n')
+		file3.close()
+
+		file4 = open("/home/pasquale/Newmark.txt","a")
+		file4.write(str(self.time)+' '+str(self.node.Newmark.Particles.position.value[1][2])+'\n')
+		file4.close()
+
+		file5 = open("/home/pasquale/EulerImplicit.txt","a")
+		file5.write(str(self.time)+' '+str(self.node.EulerImplicit.Particles.position.value[1][2])+'\n')
+		file5.close()
+
+		file5 = open("/home/pasquale/Trapezoidal.txt","a")
+		file5.write(str(self.time)+' '+str(self.node.Trapezoidal.Particles.position.value[1][2])+'\n')
+		file5.close()
+
+def createScene(rootNode):
+
+
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.Visual")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.LinearSolver.Iterative")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.ODESolver.Backward")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.Setting")
+	rootNode.addObject('RequiredPlugin', name="Sofa.GL.Component.Rendering3D")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.Constraint.Projective")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.Mass")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.SolidMechanics.Spring")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.StateContainer")
+	rootNode.addObject('RequiredPlugin', name="MyAwesomeComponents")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.Collision.Geometry")
+	rootNode.addObject('RequiredPlugin', name="Sofa.Component.ODESolver.Forward")
+	rootNode.addObject('VisualStyle', displayFlags='hideVisualModels showBehaviorModels showCollisionModels hideBoundingCollisionModels hideForceFields hideInteractionForceFields hideWireframe')
+
+	rootNode.gravity=[ 0, 0, 9.81]
+	rootNode.dt = 0.01
+	rootNode.name = 'rootNode'
+	rootNode.addObject('DefaultAnimationLoop', computeBoundingBox="0")
+
+	rootNode.addObject('BackgroundSetting', color=[0, 0.168627, 0.211765, 1])
+	rootNode.addObject('OglSceneFrame', style='Arrows', alignment='TopRight')
+
+	EE = rootNode.addChild("EulerExplicit")
+	EE.addObject('EulerExplicitSolver', name="Solver")
+	EE.addObject('CGLinearSolver', name="CG Solver", iterations="25", tolerance="1e-10", threshold="1e-10")
+	EE.addObject('MechanicalObject', name="Particles", template="Vec3d",position="0 0 0 1 0 0", velocity="0 0 0 0 0 0", showObject = True, showObjectScale = 10)
+	EE.addObject('UniformMass', name="Mass", totalMass="0.1")
+	EE.addObject('FixedConstraint', indices="0")
+	EE.addObject('StiffSpringForceField', name="Springs", spring="0 1 100 0.0 1")
+	# EE.addObject('SphereCollisionModel', radius="0.1")
+
+
+
+	RK = rootNode.addChild("RungeKutta")
+	RK.addObject('RungeKutta4Solver', name="Solver")
+	RK.addObject('CGLinearSolver', name="CG Solver", iterations="25", tolerance="1e-10", threshold="1e-10")
+	RK.addObject('MechanicalObject', name="Particles", template="Vec3d",position="0 0 0 1 0 0", velocity="0 0 0 0 0 0", showObject = True, showObjectScale = 10)
+	RK.addObject('UniformMass', name="Mass", totalMass="0.1")
+	RK.addObject('FixedConstraint', indices="0")
+	RK.addObject('StiffSpringForceField', name="Springs", spring="0 1 100 0.0 1")
+	# RK.addObject('SphereCollisionModel', radius="0.1")
+
+
+	RK = rootNode.addChild("HHT")
+	RK.addObject('HHTSolver', name="Solver")
+	RK.addObject('CGLinearSolver', name="CG Solver", iterations="25", tolerance="1e-10", threshold="1e-10")
+	RK.addObject('MechanicalObject', name="Particles", template="Vec3d",position="0 0 0 1 0 0", velocity="0 0 0 0 0 0", showObject = True, showObjectScale = 10)
+	RK.addObject('UniformMass', name="Mass", totalMass="0.1")
+	RK.addObject('FixedConstraint', indices="0")
+	RK.addObject('StiffSpringForceField', name="Springs", spring="0 1 100 0.0 1")
+	# RK.addObject('SphereCollisionModel', radius="0.1")
+
+
+
+	BDF = rootNode.addChild("Newmark")
+	BDF.addObject('NewmarkImplicitSolver', name="Solver")
+	BDF.addObject('CGLinearSolver', name="CG Solver", iterations="25", tolerance="1e-10", threshold="1e-10")
+	BDF.addObject('MechanicalObject', name="Particles", template="Vec3d",position="0 0 0 1 0 0", velocity="0 0 0 0 0 0", showObject = True, showObjectScale = 10)
+	BDF.addObject('UniformMass', name="Mass", totalMass="0.1")
+	BDF.addObject('FixedConstraint', indices="0")
+	BDF.addObject('StiffSpringForceField', name="Springs", spring="0 1 100 0.0 1")
+	# BDF.addObject('SphereCollisionModel', radius="0.1")
+
+
+
+	EI = rootNode.addChild("EulerImplicit")
+	EI.addObject('EulerImplicitSolver', name="Solver")
+	EI.addObject('CGLinearSolver', name="CG Solver", iterations="25", tolerance="1e-10", threshold="1e-10")
+	EI.addObject('MechanicalObject', name="Particles", template="Vec3d",position="0 0 0 1 0 0", velocity="0 0 0 0 0 0", showObject = True, showObjectScale = 10)
+	EI.addObject('UniformMass', name="Mass", totalMass="0.1")
+	EI.addObject('FixedConstraint', indices="0")
+	EI.addObject('StiffSpringForceField', name="Springs", spring="0 1 100 0.0 1")
+	# EI.addObject('SphereCollisionModel', radius="0.1")
+
+
+	TRAP = rootNode.addChild("Trapezoidal")
+	TRAP.addObject('EulerImplicitSolver', name="Solver", trapezoidalScheme=True)
+	TRAP.addObject('CGLinearSolver', name="CG Solver", iterations="25", tolerance="1e-10", threshold="1e-10")
+	TRAP.addObject('MechanicalObject', name="Particles", template="Vec3d",position="0 0 0 1 0 0", velocity="0 0 0 0 0 0", showObject = True, showObjectScale = 10)
+	TRAP.addObject('UniformMass', name="Mass", totalMass="0.1")
+	TRAP.addObject('FixedConstraint', indices="0")
+	TRAP.addObject('StiffSpringForceField', name="Springs", spring="0 1 100 0.0 1")
+	# TRAP.addObject('SphereCollisionModel', radius="0.1")
+
+
+	rootNode.addObject(PendulumController(node=rootNode))
+
+	return rootNode

examples/creep-test-viscohyperelastic.py

@@ -12,6 +12,7 @@
 from scipy import signal
 
 import os
+path = os.path.dirname(os.path.abspath(__file__))+'/plot/'
 
 class CylinderController(Sofa.Core.Controller):
 
@@ -35,6 +36,9 @@ def __init__(self, *args, **kwargs):
 
 		self.lin = self.node.cylinder.tetras.position.value[self.posmax1][2]
 
+		file1 = open(path + "SLS_Maxwell_cyclic.txt","w")
+		file1.write(str(0.0)+' '+str(0.0)+' '+str(0.0) +'\n')
+		file1.close()
 
 
 
@@ -47,15 +51,43 @@ def onAnimateBeginEvent(self,event):
 
 
 ##	STEP SIGNAL 
-		self.node.cylinder.CFF.totalForce.value = [0, 0, ((3.141592653589793e3))*7.8*np.heaviside(self.time, 0.0)] ## amplitude in force calculated with Matlab --> Step function 
+		#self.node.cylinder.CFF.totalForce.value = [0, 0, ((3.141592653589793e3))*7.8*np.heaviside(self.time, 0.0)] ## amplitude in force calculated with Matlab --> Step function 
 
-		if(self.time >= 1):
-			self.node.cylinder.CFF.totalForce.value = [0, 0, 0] ## amplitude in force calculated with Matlab --> Step function 
+		#if(self.time >= 1):
+		#	self.node.cylinder.CFF.totalForce.value = [0, 0, 0] ## amplitude in force calculated with Matlab --> Step function 
 
 		#self.node.cylinder.CFF.totalForce.value = [0, 0, 0] ## amplitude in force calculated with Matlab --> Step function 
 
+
+##	AMPLITUDE SWEEP 
+		self.node.cylinder.CFF.totalForce.value = [0,0,((3.141592653589793e3/2))*1*np.heaviside(self.time, 0.0)] ## amplitude in force calculated with Matlab --> Step function 
 
+		if(self.time >= self.tau*100 and self.time<= self.tau*150):
+			self.node.cylinder.CFF.totalForce.value = [0, 0, 0] ## amplitude in force calculated with Matlab --> Step function
+
 
+		if(self.time >= self.tau*150 and self.time <= self.tau*200):		
+			self.node.cylinder.CFF.totalForce.value = [0,0,((3.141592653589793e3/2))*3*np.heaviside(self.time, 0.0)] ## amplitude in force calculated with Matlab --> Step function 
+
+		if(self.time >= self.tau*200 and self.time <= self.tau*250):		
+			self.node.cylinder.CFF.totalForce.value = [0,0,0] ## amplitude in force calculated with Matlab --> Step function 
+
+		if(self.time >= self.tau*250 and self.time <= self.tau*300):		
+			self.node.cylinder.CFF.totalForce.value = [0,0,((3.141592653589793e3/2))*5*np.heaviside(self.time, 0.0)] ## amplitude in force calculated with Matlab --> Step function 
+
+		if(self.time >= self.tau*300 and self.time <= self.tau*350):		
+			self.node.cylinder.CFF.totalForce.value = [0,0,0] ## amplitude in force calculated with Matlab --> Step function 		
+
+		if(self.time >= self.tau*350 and self.time <= self.tau*400):		
+			self.node.cylinder.CFF.totalForce.value = [0,0,((3.141592653589793e3/2))*7*np.heaviside(self.time, 0.0)] ## amplitude in force calculated with Matlab --> Step function 
+
+		if(self.time >= self.tau*400):		
+			self.node.cylinder.CFF.totalForce.value = [0,0,0] ## amplitude in force calculated with Matlab --> Step function 
+
+		print(epsilon*100)
+		file1 = open(path + "SLS_Maxwell_cyclic.txt","a")
+		file1.write(str(self.time)+' '+str(self.node.cylinder.FEM.stressVonMisesElement.value[4])+' '+str(epsilon*100)+ '\n' )
+		file1.close()
 
 
 
@@ -101,7 +133,7 @@ def createScene(rootNode):
 ## SLS-MAXWELL FIRST ORDER Cylinder: Material F5000- R0.5  (BLUE)
 	cylinder = rootNode.addChild('cylinder')
 
-	cylinder.addObject('EulerImplicitSolver', name="Solver",rayleighMass = 0.0, rayleighStiffness = 0.0)
+	cylinder.addObject('EulerImplicitSolver', name="Solver",rayleighMass = 0.1, rayleighStiffness = 0.1)
 
 	cylinder.addObject('CGLinearSolver', name="ItSolver", iterations="250", tolerance="1e-30", threshold = '1e-12')
 	cylinder.addObject('MeshVTKLoader', name='loader', filename='mesh/cylinder5296.vtk', translation = [0, 0.0, 0])

examples/creep-test.py

@@ -12,6 +12,7 @@
 from scipy import signal
 
 import os
+path = os.path.dirname(os.path.abspath(__file__))+'/plot/'
 
 class CylinderController(Sofa.Core.Controller):
 
@@ -35,28 +36,33 @@ def __init__(self, *args, **kwargs):
 
 		self.lin = self.node.cylinder.tetras.position.value[self.posmax1][2]
 
-
+		file1 = open(path + "SLS_Maxwell_cyclic1.txt","w")
+		file1.write(str(0.0)+' '+str(0.0)+' '+str(0.0) +'\n')
+		file1.close()
 
 
 
 	def onAnimateBeginEvent(self,event):
 		self.time = self.node.time.value
-		self.tau = self.node.cylinder.FEM.ParameterSet.value[2] 
+		#self.tau = self.node.cylinder.FEM.ParameterSet.value[2] 
 		epsilon = (self.node.cylinder.tetras.position.value[4][2]-self.lin)/self.lin
 
 ## IN THIS CODE WE WILL DO A CREEP TEST, SO WE WILL APPLY A STEP AS INPUT, USING THE CONSTANTFORCEFIELD LOAD.
 
-
 ##	STEP SIGNAL 
 		self.node.cylinder.CFF.totalForce.value = [0,0,((3.141592653589793e2/2))*np.heaviside(self.time, 0.0)] ## amplitude in force calculated with Matlab --> Step function 
+## Ramp
 
-		if(self.time >= 1):
-			self.node.cylinder.CFF.totalForce.value = [0, 0, 0] ## amplitude in force calculated with Matlab --> Step function 
+		#self.node.cylinder.CFF.totalForce.value =  [0,0,((3.141592653589793e2/2))*(1-abs(signal.sawtooth(2 * np.pi*7*self.time)))]## amplitude in force calculated with Matlab 
+		#if(self.time>=0.1):
+		#	self.node.cylinder.CFF.totalForce.value =  [0,0,((3.141592653589793e2/2))]## amplitude in force calculated with Matlab 
 
 
 
 		print(epsilon*100)
-
+		#file1 = open(path + "SLS_Maxwell_cyclic1.txt","a")
+		#file1.write(str(self.time)+' '+str(self.node.cylinder.FEM.stressVonMisesElement.value[4])+' '+str(epsilon*100)+ '\n' )
+		#file1.close()
 
 
 
@@ -92,7 +98,7 @@ def createScene(rootNode):
 	rootNode.dt = (1e6/(20e6*100))
 	rootNode.name = 'rootNode'
 	rootNode.addObject('DefaultAnimationLoop', computeBoundingBox="0")
-	rootNode.addObject('GenericConstraintSolver', tolerance=1e-24, maxIterations=100000000)
+	rootNode.addObject('GenericConstraintSolver', tolerance=1e-24, maxIterations=1000)
 	rootNode.addObject('OglSceneFrame', style='Arrows', alignment='TopRight')
 
 
@@ -104,7 +110,7 @@ def createScene(rootNode):
 
 	cylinder.addObject('EulerImplicitSolver', name="Solver",rayleighMass = 0.0, rayleighStiffness = 0.0, firstOrder = True, trapezoidalScheme = False)
 
-	cylinder.addObject('CGLinearSolver', name="ItSolver", iterations="25000000", tolerance="1e-15", threshold = '1e-30')
+	cylinder.addObject('CGLinearSolver', name="ItSolver", iterations="2500", tolerance="1e-30", threshold = '1e-12')
 	cylinder.addObject('MeshVTKLoader', name='loader', filename='mesh/cylinder5296.vtk', translation = [0, 0.0, 0])
 	cylinder.addObject('MechanicalObject', name='tetras', template='Vec3d', src = '@loader')
 	cylinder.addObject('TetrahedronSetTopologyContainer', name="topo", src ='@loader')
@@ -126,7 +132,8 @@ def createScene(rootNode):
 	tau2 = 1e6/E2
 	tau3 = 1e6/E3
 	nu = 0.44 ## Poisson's ratio
-	cylinder.addObject('TetrahedronViscoelasticityFEMForceField', template='Vec3d', name='FEM', src ='@topo',materialName="SLSKelvinVoigtFirstOrder", ParameterSet= str(E1)+' '+str(E2)+' '+str(tau2)+' '+str(nu))
+	#cylinder.addObject('TetrahedronViscoelasticityFEMForceField', template='Vec3d', name='FEM', src ='@topo',materialName="SLSMaxwellFirstOrder", ParameterSet= str(E1)+' '+str(E2)+' '+str(tau2)+' '+str(nu))
+	cylinder.addObject('TetrahedronFEMForceField', template='Vec3d', name='FEM', src ='@topo',youngModulus = E1, poissonRatio = nu)
 
 	cylinder.addObject('ConstantForceField', name = "CFF", listening = True, totalForce =[0,0,0],template="Vec3d", src= "@topo", indices = cylinder.boxToPull.indices.linkpath)