Skip to content

Commit

Permalink
Merge pull request #29 from SofaDefrost/add_new_models
Browse files Browse the repository at this point in the history 
new-equations + stability
  • Loading branch information
@pasqualeferr94
pasqualeferr94 authored Feb 19, 2025
2 parents 36c17b0 + d15108c commit 3a57ad1
Show file tree
Hide file tree
Showing 29 changed files with 488 additions and 746 deletions.
8 changes: 0 additions & 8 deletions Bindings/src/Binding_TetrahedronViscoElasticityFEMForceField.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,13 +95,6 @@ py::object getF(Myclass& self, int i)
return pybind11::cast(tmp);
}

//Vec 6
py::object getSPKStress(Myclass& self, int i)
{
auto tmp = self.m_tetrahedronInfo.getValue()[i].m_SPKStress;

return pybind11::cast(tmp);
}

//Vec 6
py::object getCauchyStress(Myclass& self, int i)
Expand Down Expand Up @@ -134,7 +127,6 @@ void moduleAddTetrahedronViscoelasticityFEMForceField(py::module &m)
p.def("getRestVolume", getRestVolume, "get the rest volume of a tetrahedron");
p.def("getVolScale", getVolScale, "get the volume Scale of a tetrahedron");
p.def("getF", getF, "get the deformation gradient of a tetrahedron");
p.def("getSPKStress", getSPKStress, "get the Second Piola Kirchhoff Stresses of a tetrahedron");
p.def("getCauchyStress", getCauchyStress, "get the Cauchy Stresses of a tetrahedron");
p.def("getVonMisesStress", getVonMisesStress, "get the Von Mises Stresses of a tetrahedron");

Expand Down
43 changes: 19 additions & 24 deletions examples/creep-test-viscohyperelastic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,10 +35,7 @@ def __init__(self, *args, **kwargs):
print(self.posmax1)

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()




Expand All @@ -60,34 +57,32 @@ def onAnimateBeginEvent(self,event):


## 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
self.node.cylinder.CFF.totalForce.value = [0,0,((3.141592653589793e2/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*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*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*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*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*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*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
#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()




Expand Down Expand Up @@ -133,9 +128,9 @@ def createScene(rootNode):
## SLS-MAXWELL FIRST ORDER Cylinder: Material F5000- R0.5 (BLUE)
cylinder = rootNode.addChild('cylinder')

cylinder.addObject('EulerImplicitSolver', name="Solver",rayleighMass = 0.1, rayleighStiffness = 0.1)
cylinder.addObject('EulerImplicitSolver', name="Solver",rayleighMass = 0.0, rayleighStiffness = 0.0)

cylinder.addObject('CGLinearSolver', name="ItSolver", iterations="250", tolerance="1e-30", threshold = '1e-12')
cylinder.addObject('CGLinearSolver', name="ItSolver", iterations="250", tolerance="1e-12", 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')
Expand Down
20 changes: 9 additions & 11 deletions examples/creep-test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,9 +36,6 @@ 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()



Expand Down Expand Up @@ -125,15 +122,16 @@ def createScene(rootNode):

cylinder.addObject('BoxROI', name="boxToPull", box=[-0.011, -0.011, 0.05, 0.011, 0.011, 0.051], drawBoxes=False)

E1 = 70e6
E2 = 20e6
E3 = 10e6
tau1 = 1e6/E1
tau2 = 1e6/E2
tau3 = 1e6/E3
G1 = 70e6
G2 = 20e6
G3 = 10e6
tau1 = 1e6/G1
tau2 = 1e6/G2
tau3 = 1e6/G3
nu = 0.44 ## Poisson's ratio
#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)
lamb = 2*G1*nu/(1-2*nu) ## relationship from wikipedia: https://en.wikipedia.org/wiki/Lam%C3%A9_parameters

cylinder.addObject('TetrahedronViscoelasticityFEMForceField', template='Vec3d', name='FEM', src ='@topo',materialName="SLSMaxwellFirstOrder", ParameterSet= str(G1)+' '+str(lamb)+' '+str(G2)+' '+str(tau2))

cylinder.addObject('ConstantForceField', name = "CFF", listening = True, totalForce =[0,0,0],template="Vec3d", src= "@topo", indices = cylinder.boxToPull.indices.linkpath)

Expand Down
20 changes: 10 additions & 10 deletions examples/stress-map.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,6 @@ def onAnimateBeginEvent(self,event):


stress = (self.node.cylinder.FEM.CauchyStress.value)/1e6 ## MPa
stressPerNode = np.array([[0,0,0,0],[0,0,0,0]])



Expand All @@ -48,7 +47,7 @@ def onAnimateBeginEvent(self,event):
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


print(self.stress[self.posmax1][2]/1e6)


## IN THIS CODE WE WILL DO A STRESS RELAXATION TEST, SO WE WILL APPLY A STEP AS INPUT, USING THE POSITIONCONSTRAINT.
Expand Down Expand Up @@ -84,7 +83,7 @@ def createScene(rootNode):
rootNode.addObject('FreeMotionAnimationLoop')
rootNode.addObject('GenericConstraintSolver', maxIterations=1e4, tolerance=1e-50)
rootNode.gravity = [0,0,-9.81]
rootNode.dt = (1e9/(20e9*1000))
rootNode.dt = (1e6/(20e6*100))

rootNode.addObject('VisualStyle', displayFlags='hideForceFields')
rootNode.addObject('OglSceneFrame', style='Arrows', alignment='TopRight')
Expand All @@ -101,15 +100,16 @@ def createScene(rootNode):
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
G1 = 70e6
tau1 = 1e6/G1
G2 = 20e6
tau2 = 1e6/G2
G3 = 10e6
tau3 = 1e9/G3
nu = 0.44
lamb = 30e6

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('TetrahedronViscoelasticityFEMForceField', template='Vec3d', name='FEM', src ='@topo',materialName="MaxwellFirstOrder", ParameterSet= str(G1)+' '+str(tau1)+' '+str(lamb))

cylinder.addObject('BoxROI', name='boxROI',box="-0.011 -0.011 -0.001 0.011 0.011 0.001", drawBoxes=True)
cylinder.addObject('FixedProjectiveConstraint', indices = '@boxROI.indices')
Expand Down
31 changes: 12 additions & 19 deletions examples/stress-relaxation-test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,6 @@ def __init__(self, *args, **kwargs):
print(self.pos3[self.posmax1][2], self.posmax1)

self.lin = self.node.cylinder.tetras.position.value[self.posmax1][2]

file1 = open(path + "SLS_Maxwell_relaxation.txt","w")
file1.write(str(0.0)+' '+str(0.0)+' '+str(0.0) +'\n')
file1.close()



def onAnimateEndEvent(self,event):
Expand All @@ -48,9 +43,6 @@ def onAnimateEndEvent(self,event):


## IN THIS CODE WE WILL DO A STRESS RELAXATION TEST, SO WE WILL APPLY A STEP AS INPUT, USING THE POSITIONCONSTRAINT.
file1 = open(path + "SLS_Maxwell_relaxation.txt","a")
file1.write(str(self.time)+' '+str(epsilon*100)+' '+str(self.stress[2]/1e6) +'\n' )
file1.close()


def createScene(rootNode):
Expand Down Expand Up @@ -80,9 +72,9 @@ def createScene(rootNode):


rootNode.addObject('FreeMotionAnimationLoop')
rootNode.addObject('GenericConstraintSolver', maxIterations=1e4, tolerance=1e-50)
rootNode.gravity = [0,-9810,0]
rootNode.dt = (1e9/(20e9*100))
rootNode.addObject('GenericConstraintSolver', maxIterations=1e4, tolerance=1e-12)
rootNode.gravity = [0,0, -9.81]
rootNode.dt = (1e6/(20e6*100))

rootNode.addObject('VisualStyle', displayFlags='hideForceFields')
rootNode.addObject('OglSceneFrame', style='Arrows', alignment='TopRight')
Expand All @@ -99,15 +91,16 @@ def createScene(rootNode):
cylinder.addObject('TetrahedronSetGeometryAlgorithms', template="Vec3d" ,name="GeomAlgo")

cylinder.addObject('UniformMass', totalMass="0.0126", src = '@topo')
E1 = 70e6
tau1 = 1e9/E1
E2 = 20e6
tau2 = 1e9/E2
E3 = 10e6
tau3 = 1e9/E3
G1 = 70e6
tau1 = 1e6/G1
G2 = 20e6
tau2 = 1e6/G2
G3 = 10e6
tau3 = 1e6/G3
nu = 0.44

cylinder.addObject('TetrahedronViscoelasticityFEMForceField', template='Vec3d', name='FEM', src ='@topo',materialName="MaxwellFirstOrder", ParameterSet= str(E1)+' '+str(tau1)+' '+str(nu))
#lamb = 2*G1*nu/(1-2*nu) ## relationship from wikipedia: https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
lamb = 0
cylinder.addObject('TetrahedronViscoelasticityFEMForceField', template='Vec3d', name='FEM', src='@topo', materialName="Burgers", ParameterSet=str(G1) + ' '+ str(tau1) + ' ' + str(G2)+' '+str(tau2)+ ' '+ str(lamb))

cylinder.addObject('BoxROI', name='boxROI',box="-0.011 -0.011 -0.001 0.011 0.011 0.001", drawBoxes=True)
cylinder.addObject('FixedProjectiveConstraint', indices = '@boxROI.indices')
Expand Down
3 changes: 0 additions & 3 deletions src/SofaViscoElastic/TetrahedronViscoHyperelasticityFEMForceField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,9 +114,6 @@ class TetrahedronViscoHyperelasticityFEMForceField : public core::behavior::Forc
/// Second Piola Kirchhoff stress tensor
MatrixSym m_SPKTensorGeneral;

/// Cauchy Green stress Tensor
MatrixSym m_CauchyStressTensor;

/// Von Mises Stress
Real m_VonMisesStress;

Expand Down
32 changes: 7 additions & 25 deletions src/SofaViscoElastic/TetrahedronViscoHyperelasticityFEMForceField.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,6 @@ template <class DataTypes> TetrahedronViscoHyperelasticityFEMForceField<DataType
, m_tetrahedronInfo(initData(&m_tetrahedronInfo, "tetrahedronInfo", "Internal tetrahedron data"))
, m_edgeInfo(initData(&m_edgeInfo, "edgeInfo", "Internal edge data"))
, d_stressSPK(initData(&d_stressSPK, "stressSPK","The stress of the Second Piola Kirchhoff Stress Tensor per Element"))
, d_Cauchystress(initData(&d_Cauchystress, "CauchyStress","The stress of the Cauchy Stress Tensor per Element"))
, d_stressVonMisesElement(initData(&d_stressVonMisesElement, "stressVonMisesElement","The stress of the Von Mises Stress per Element"))
, d_stressVonMisesNode(initData(&d_stressVonMisesNode, "stressVonMisesNode","The stress of the Von Mises Stress per Node"))

Expand Down Expand Up @@ -137,11 +136,7 @@ void TetrahedronViscoHyperelasticityFEMForceField<DataTypes>::instantiateMateria
{
m_myMaterial = std::make_unique<SLSOgdenSecondOrder<DataTypes>>();
}
//else if (material == "SLSKelvinVoigtSecondOrder")
//{

// m_myMaterial = std::make_unique<SLSKelvinVoigtSecondOrder<DataTypes>>();
//}

else
{
msg_error() << "material name " << material <<
Expand Down Expand Up @@ -340,7 +335,6 @@ void TetrahedronViscoHyperelasticityFEMForceField<DataTypes>::addForce(const cor
Coord dp[3], x0, sv;

vector<Vec6d> vecStressSPK;
vector<Vec6d> vecStressCauchy;

vector<Real> vecStressVonMisesElement;

Expand Down Expand Up @@ -414,38 +408,29 @@ void TetrahedronViscoHyperelasticityFEMForceField<DataTypes>::addForce(const cor
tetInfo->C(2, 2));



tetInfo->m_SPKTensorGeneral.clear();
tetInfo->m_CauchyStressTensor.clear();

tetInfo->m_SPKStress.clear();
tetInfo->m_CauchyStress.clear();

m_myMaterial->deriveSPKTensor(tetInfo, globalParameters, tetInfo->m_SPKTensorGeneral,tetInfo->m_CauchyStressTensor,dt);
m_myMaterial->deriveSPKTensor(tetInfo, globalParameters, tetInfo->m_SPKTensorGeneral,dt);

/// Stress map in simulation (Voigt notation of the Stress tensor)
vecStressSPK.push_back(tetInfo->m_SPKTensorGeneral.getVoigt());
vecStressCauchy.push_back(tetInfo->m_CauchyStressTensor.getVoigt());

/// Python Bindings
tetInfo->m_SPKStress = tetInfo->m_SPKTensorGeneral.getVoigt();
tetInfo->m_CauchyStress = tetInfo->m_CauchyStressTensor.getVoigt();

/// Calculation of the General Von mises Stress per Element
tetInfo->m_VonMisesStress = sqrt(0.5*((tetInfo->m_CauchyStress(0)-tetInfo->m_CauchyStress(1))*(tetInfo->m_CauchyStress(0)-tetInfo->m_CauchyStress(1)) +
(tetInfo->m_CauchyStress(1)-tetInfo->m_CauchyStress(2))*(tetInfo->m_CauchyStress(1)-tetInfo->m_CauchyStress(2)) +
(tetInfo->m_CauchyStress(2)-tetInfo->m_CauchyStress(0))*(tetInfo->m_CauchyStress(2)-tetInfo->m_CauchyStress(0)) +
3*(tetInfo->m_CauchyStress(3)*tetInfo->m_CauchyStress(3)+tetInfo->m_CauchyStress(4)*tetInfo->m_CauchyStress(4) +
tetInfo->m_CauchyStress(5)*tetInfo->m_CauchyStress(5))
tetInfo->m_VonMisesStress = sqrt(0.5*((tetInfo->m_SPKStress(0)-tetInfo->m_SPKStress(1))*(tetInfo->m_SPKStress(0)-tetInfo->m_SPKStress(1)) +
(tetInfo->m_SPKStress(1)-tetInfo->m_SPKStress(2))*(tetInfo->m_SPKStress(1)-tetInfo->m_SPKStress(2)) +
(tetInfo->m_SPKStress(2)-tetInfo->m_SPKStress(0))*(tetInfo->m_SPKStress(2)-tetInfo->m_SPKStress(0)) +
3*(tetInfo->m_SPKStress(3)*tetInfo->m_SPKStress(3)+tetInfo->m_SPKStress(4)*tetInfo->m_SPKStress(4) +
tetInfo->m_SPKStress(5)*tetInfo->m_SPKStress(5))
));

vecStressVonMisesElement.push_back(tetInfo->m_VonMisesStress);






for (l = 0; l < 4; ++l)
{
f[ta[l]] -= tetInfo->m_deformationGradient* (
Expand All @@ -460,9 +445,6 @@ void TetrahedronViscoHyperelasticityFEMForceField<DataTypes>::addForce(const cor
d_stressSPK.setValue(vecStressSPK);
d_stressSPK.endEdit();

d_Cauchystress.beginEdit();
d_Cauchystress.setValue(vecStressCauchy);
d_Cauchystress.endEdit();

d_stressVonMisesElement.beginEdit();
d_stressVonMisesElement.setValue(vecStressVonMisesElement);
Expand Down
4 changes: 0 additions & 4 deletions src/SofaViscoElastic/TetrahedronViscoelasticityFEMForceField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,8 +111,6 @@ class TetrahedronViscoelasticityFEMForceField : public core::behavior::ForceFiel
/// current tetrahedron volume
Real m_volScale{};
Real m_volume{};
/// Second Piola Kirchhoff stress tensor
MatrixSym m_SPKTensorGeneral;

/// Cauchy Green stress Tensor
MatrixSym m_CauchyStressTensor;
Expand All @@ -121,7 +119,6 @@ class TetrahedronViscoelasticityFEMForceField : public core::behavior::ForceFiel
Real m_VonMisesStress;

/// to export the stresses with Python Binding
Vec6 m_SPKStress;
Vec6 m_CauchyStress;

/// deformation gradient = gradPhi
Expand Down Expand Up @@ -164,7 +161,6 @@ class TetrahedronViscoelasticityFEMForceField : public core::behavior::ForceFiel

TetrahedronData<sofa::type::vector<TetrahedronRestInformation> > m_tetrahedronInfo; ///< Internal tetrahedron data
EdgeData<sofa::type::vector<EdgeInformation> > m_edgeInfo; ///< Internal edge data
Data<vector<Vec6d> > d_stressSPK;
Data<vector<Vec6d> > d_Cauchystress;
Data<vector<Real> > d_stressVonMisesElement;
Data<vector<Real> > d_stressVonMisesNode;
Expand Down
Loading

0 comments on commit 3a57ad1

Please sign in to comment.