-
Notifications
You must be signed in to change notification settings - Fork 6
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
I added new material models for viscoelasticity
- Loading branch information
1 parent
0c88329
commit 85489e0
Showing
19 changed files
with
1,076 additions
and
28 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,131 @@ | ||
| import SofaRuntime | ||
| SofaRuntime.importPlugin("SofaComponentAll") | ||
|
|
||
| # to add elements like Node or objects | ||
| import Sofa.Core | ||
| root = Sofa.Core.Node() | ||
| import SofaViscoElastic | ||
| import math | ||
| import numpy as np | ||
|
|
||
| import os | ||
| path = os.path.dirname(os.path.abspath(__file__))+'/plot/' | ||
|
|
||
|
|
||
| class CylinderController(Sofa.Core.Controller): | ||
|
|
||
| def __init__(self, *args, **kwargs): | ||
| Sofa.Core.Controller.__init__(self,*args, **kwargs) #needed | ||
| self.time = 0.0 | ||
| self.node = kwargs['node'] | ||
| self.pos3 = kwargs['pos'] | ||
|
|
||
| self.max1 = 0 | ||
| self.posmax1 = 0 | ||
| self.max2 = 0 | ||
| self.posmax2 = 0 | ||
| for j in range(0, len(self.pos3)): | ||
| if self.pos3[j][2] >= self.max1 : | ||
| self.max1 = self.pos3[j][2] | ||
| self.posmax1 = j | ||
|
|
||
| print(self.pos3[self.posmax1][2], self.posmax1) | ||
|
|
||
| self.lin = self.node.cylinder.tetras.position.value[self.posmax1][2] | ||
|
|
||
|
|
||
|
|
||
| def onAnimateBeginEvent(self,event): | ||
|
|
||
|
|
||
| stress = (self.node.cylinder.FEM.CauchyStress.value)/1e6 ## MPa | ||
| stressPerNode = np.array([[0,0,0,0],[0,0,0,0]]) | ||
|
|
||
|
|
||
|
|
||
| self.time = self.node.time.value | ||
| epsilon = (self.node.cylinder.tetras.position.value[self.posmax1][2]-self.lin)/self.lin | ||
| self.node.cylinder.visu.display.pointData.value = (stress[0:len(self.pos3),2]) | ||
| self.node.cylinder.visu.Map.min.value = self.node.cylinder.visu.display.currentMin.value | ||
| self.node.cylinder.visu.Map.max.value = self.node.cylinder.visu.display.currentMax.value | ||
|
|
||
|
|
||
|
|
||
|
|
||
| ## IN THIS CODE WE WILL DO A STRESS RELAXATION TEST, SO WE WILL APPLY A STEP AS INPUT, USING THE POSITIONCONSTRAINT. | ||
|
|
||
|
|
||
|
|
||
| def createScene(rootNode): | ||
| rootNode.addObject('RequiredPlugin', name='SofaPython3') | ||
| rootNode.addObject('RequiredPlugin', name="Sofa.Component.Engine.Select") | ||
| rootNode.addObject('RequiredPlugin', name="Sofa.Component.IO.Mesh") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.LinearSolver.Iterative") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.Mapping.Linear") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.Mass") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.ODESolver.Forward") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.Setting") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.SolidMechanics.FEM.HyperElastic") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.SolidMechanics.FEM.Elastic") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.SolidMechanics.Spring") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.StateContainer") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.Topology.Container.Dynamic") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.Visual") | ||
| rootNode.addObject("RequiredPlugin", name= "Sofa.GL.Component.Rendering3D") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.AnimationLoop") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.Constraint.Lagrangian.Solver") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.MechanicalLoad") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.LinearSolver.Direct") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.Constraint.Lagrangian.Correction") | ||
| rootNode.addObject("RequiredPlugin", name = "Sofa.Component.Constraint.Projective") | ||
| rootNode.addObject("RequiredPlugin", name="Sofa.Component.ODESolver.Backward") | ||
|
|
||
|
|
||
| rootNode.addObject('FreeMotionAnimationLoop') | ||
| rootNode.addObject('GenericConstraintSolver', maxIterations=1e4, tolerance=1e-50) | ||
| rootNode.gravity = [0,0,-9.81] | ||
| rootNode.dt = (1e9/(20e9*1000)) | ||
|
|
||
| rootNode.addObject('VisualStyle', displayFlags='hideForceFields') | ||
| rootNode.addObject('OglSceneFrame', style='Arrows', alignment='TopRight') | ||
|
|
||
| cylinder = rootNode.addChild('cylinder') | ||
|
|
||
| cylinder.addObject('EulerImplicitSolver', name="Solver",rayleighMass = 0.0, rayleighStiffness = 0.0, firstOrder = True, trapezoidalScheme = False) | ||
| cylinder.addObject('SparseLDLSolver', name="directsolver") | ||
|
|
||
| cylinder.addObject('MeshVTKLoader', name='loader', filename='mesh/cylinder9178.vtk',) | ||
| cylinder.addObject('MechanicalObject', name='tetras', template='Vec3d', src = '@loader', rotation = [0, 0 , 0]) | ||
| cylinder.addObject('TetrahedronSetTopologyContainer', name="topo", src ='@loader') | ||
| cylinder.addObject('TetrahedronSetTopologyModifier' , name="Modifier") | ||
| cylinder.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d" ,name="GeomAlgo") | ||
|
|
||
| cylinder.addObject('UniformMass', totalMass="0.0126", src = '@topo') | ||
| E1 = 70e9 | ||
| tau1 = 1e9/E1 | ||
| E2 = 20e9 | ||
| tau2 = 1e9/E2 | ||
| E3 = 10e9 | ||
| tau3 = 1e9/E3 | ||
| nu = 0.44 | ||
|
|
||
| cylinder.addObject('TetrahedronViscoelasticityFEMForceField', template='Vec3d', name='FEM', src ='@topo',materialName="SLSKelvinVoigtSecondOrder", ParameterSet= str(E1)+' '+str(E2)+' '+str(tau2)+' '+str(E3)+' '+str(tau3)+' '+str(nu)) | ||
|
|
||
| cylinder.addObject('BoxROI', name='boxROI',box="-0.011 -0.011 -0.001 0.011 0.011 0.001", drawBoxes=True) | ||
| cylinder.addObject('FixedConstraint', indices = '@boxROI.indices') | ||
| cylinder.addObject('BoxROI', name="boxToPull", box=[-0.011, -0.011, 0.1, 0.011, 0.011, 0.101], drawBoxes=True) | ||
| cylinder.addObject('PartialFixedConstraint', indices=cylinder.boxToPull.indices.linkpath, fixedDirections=[1, 1, 0]) | ||
|
|
||
| ## STEP SIGNAL | ||
| cylinder.addObject('PositionConstraint', name = 'displacement', indices=cylinder.boxToPull.indices.linkpath,valueType="displacement", value = 1e-4 , useDirections=[0, 0, 1]) | ||
| cylinder.addObject('LinearSolverConstraintCorrection') | ||
|
|
||
| cylinder.addObject(CylinderController(node=rootNode, pos = rootNode.cylinder.tetras.position)) | ||
|
|
||
|
|
||
|
|
||
| modelVisu1 = cylinder.addChild('visu') | ||
| modelVisu1.addObject('DataDisplay', name = 'display') | ||
| modelVisu1.addObject('OglColorMap', name = 'Map', colorScheme = "Blue to Red", showLegend = True ) | ||
| modelVisu1.addObject('IdentityMapping',input="@..", output="@.") | ||
| return rootNode |
Oops, something went wrong.