Home › Forum › SOFA › Using SOFA › SOFA scene with imposed displacements and Von Mises stress computation in 2D
- This topic has 13 replies, 5 voices, and was last updated 2 years, 8 months ago by
Hugo.
-
AuthorPosts
-
26 May 2020 at 18:31 #16386
pvillard
BlockedI am trying to write a simple SOFA scene including imposed displacements to a given geometry and export the results with the Von Mises stress for each node.
My scene is currently the following:
<?xml version='1.0'?> <Node name="root" dt="0.001" gravity="0 0 0" > <VisualStyle displayFlags="hideVisualModels hideBehaviorModels hideCollisionModels hideMappings showForceFields showInteractionForceFields showWireframe" /> <BackgroundSetting color="1 1 1 1" /> <Node name="object"> <EulerImplicit name="cg_odesolver" printLog="false" rayleighStiffness="0.01" rayleighMass="0.01" /> <CGLinearSolver /> <MeshObjLoader name="meshLoader" filename="object.obj" /> <TriangleSetTopologyContainer name="Container" src="@meshLoader" /> <TriangleSetGeometryAlgorithms name="GeomAlgo"/> <MechanicalObject name="MO" showIndicesScale="0.01" showIndices="false" /> <DiagonalMass massDensity="1050" /> <!-- 1.05 g/cm^3 according to https://bionumbers.hms.harvard.edu/files/Density%20and%20mass%20of%20each%20organ-tissue.pdf --> <TriangularFEMForceField name="FEM" youngModulus="100000" poissonRatio="0.3" method="large" /> <!-- data from the literature --> <MeshObjLoader name="displ" filename="displ.obj"/> <MeshObjLoader name="indices" filename="indices.obj"/> <AffineMovementConstraint indices="@indices.position" meshIndices = "@indices.position" translation="@displ.position" drawConstrainedPoints="1"/> </Node> </Node>
1) To impose the displacements, I have a list of indices given by the file “indices.obj”, I have the 2D displacement vectors given by the file “displ.obj”. I thought the command to use would be
AffineMovementConstraint
but it doesn’t seem to work as the translation should only be a 3D vector and not a list of 2D vectors as I have. Anyway, by using that, I have continuous displacements that are never-ending. What is the proper way to impose displacement on a list of nodes within a SOFA scene?2) I didn’t find how to export the stresses from the
WriteState
export command. Is there a way to do it without modifying SOFA code and recompiling it?Thank you in advance for any answer.
Best,
Pierre
31 May 2020 at 23:45 #16482Hugo
KeymasterHey @pvillard
1) The field indices of the
AffineMovementConstraint
is a list of integers corresponding to the indices of the constraint DOFs. I therefore don’t think that providing the field position of the loader as the input of indices is appropriate. If you want to constraint all points, you can create a box bounding the object and use the output indices.
As you noticed, the field translation is a Vector3 (you can make this translation in plane, i.e. 2D) and the rotation is a RotationMatrix. This one transformation will be applied on all selected indices.Do you want to constraint a different displacement for each point?
2)
WriteState
indeed exports the MechanicalState. To export any data, I would recommand to use the VTKExporter.Best,
Hugo
3 June 2020 at 18:04 #16530pvillard
BlockedHi Hugo,
Thank you very much for your (fast) answer.
1) Instead of trying to assign my 974 imposed displacements with 1 line (passing the indices with an obj file threw the attribute position), I wrote 974 lines as follow:
<AffineMovementConstraint indices="i" meshIndices = "i" translation="-x -y 0" drawConstrainedPoints="1"/>
where the vertex of indexi
has the displacement(x,y)
. I still have no warning and no error and the displayed vertices with the AffineMovementConstraint are the good ones. I have the same problem in the simulation as before.2) I managed to export the nodal forces with:
<VTKExporter name="vTKExporter" listening="true" pointsDataFields="MO.force" ...
. How can I get the stress instead?
Otherwise, is it possible to directly visualize Von-Mises stresses within runSofa?Thank you,
Pierre
12 June 2020 at 01:32 #16607Noura
BlockedHello Pierre,
I think that you should be able to visualize the stresses using the following parameters from “TriangularFEMForceField.h”
Data<bool> showStressValue; Data<bool> showStressVector;
Try this:
<TriangularFEMForceField name="FEM" youngModulus="100000" poissonRatio="0.3" method="large" showStressValue="true" showStressVector="false" />
You can refer to the function “computePrincipalStress” from “TriangularFEMForceField.h” to verify if this is what you need.
Otherwise, for Von Mises stresses, I’m not sure if they are implemented for the triangular FF. In the case of a tetrahedron FF you can visualize them for a FF assigned on a given volume mesh by:
<TetrahedronFEMForceField template="Vec3d" name="FEM" printLog="1" method="svd" poissonRatio="0.48" youngModulus="16.34" computeVonMisesStress="1" isToPrint="false" showVonMisesStressPerNode="false" showStressColorMap="blue 1 0 0 1 1 0.5 0.5 1" updateStiffness="true" />
Hope this helps!
Noura
12 June 2020 at 14:17 #16615pvillard
BlockedHi Noura,
Thank you very much for your help. I have missed these boolean data in
TriangularFEMForceField
. Now, I have a stress visualization within runSOFA.Do you know by any chance how I can export these data? I know I have to go through a
VTKExporter
and I guess from a field ofFEM
(the name of myTriangularFEMForceField
) but I didn’t figure out one that worksThank you in advance,
Pierre
13 June 2020 at 01:09 #16618Noura
BlockedHi Pierre,
I have no clue how to export the stress values directly via the “VTKExporter” component. Maybe it is possible !
Also, I’m not sure if it is possible to get the stresses directly in the scene since that there is no data field called “stress” in the “TriangularFEMForceField” (FF).Nevertheless, the stress values can be accessed and exported to an output stream by slightly modifying Sofa code.
I’ll explain in detail how to do it, but it is really simple.
1- There is a parameter called
f_computePrincipalStress
in “TriangularFEMForceField.h”2- The previous parameter is not activated by default in the constructor of the FF in “TriangularFEMForceField.inl”
,f_computePrincipalStress(initData(&f_computePrincipalStress,false,"computePrincipalStress","Compute principal stress for each triangle"))
3- In your scene, assign “true” to the attribute “computePrincipalStress” in the FF
<TriangularFEMForceField name="FEM" youngModulus="10000" poissonRatio="0.3" method="large" showStressValue="true" computePrincipalStress = "true" />
4- Note that a handy class called
TriangleInformation
is defined in “TriangularFEMForceField.h”. The latter contains many information related to the FF (including the stresses).
Further, this class has an output stream operatoroperator<<
which enables writing/exporting the values to an output file for instance.5- Then, you have to seek where to export them in Sofa code. An option would be to do that in the “AddForce” function where the boolean is examined, there you have
if (f_computePrincipalStress.getValue()) ... for(unsigned int i=0; i<nbTriangles; ++i) computePrincipalStress(i, triangleInf[i].stress); ...
You may replace it by
... std::ofstream file_principal_stresses ("stress.csv",std::ios::app ); if (f_computePrincipalStress.getValue()) for(unsigned int i=0; i<nbTriangles; ++i){ computePrincipalStress(i, triangleInf[i].stress); file_principal_stresses<<triangleInf[i].stress<<","; } ...
Finally, that was a quick way to get them, BUT I would not do it that way to avoid unnecessary writing to a stream each time “AddForce” is called!
Instead, you can define an additional parameter in “TriangularFEMForceField.h” to be activated only when writing to a file is desiredData<bool> f_exportStressValuesToFile;
Then, you initialize it in “TriangularFEMForceField.inl” as all other attributes:
TriangularFEMForceField<DataTypes>::TriangularFEMForceField() ... ,f_exportStressValuesToFile(initData(&f_exportStressValuesToFile,false,"exportPrincipalStress","Export principal stress for each triangle")) ...
You include exportPrincipalStress = “true” in the FF component in the scene, and you check if it is true somewhere in the code to write the stresses or not.
An additional parameter to hold the name of the file can be added in the same manner.Hope that this solves your issue in case you couldn’t figure out how to do it in a simpler way, otherwise I would be glad to learn how 🙂
Noura
15 June 2020 at 11:40 #16628Hugo
Keymaster15 June 2020 at 23:49 #16641Noura
BlockedHi @hugo,
I haven’t implemented it. I only had a look at the code to propose a quick solution. But I’ll try to dedicate some time to implement it (because I have to upgrade to the latest SOFA version first).
Noura
19 June 2020 at 16:10 #16676Rayan
BlockedHi Pierre,
I don’t know if it could help you, but I’ve wrote a component, linked to a mechanical object and which returns the current strain applied on this mechanical state. I am using deformation gradient. If you want to get the stress, maybe you can easily link it to the strain, depending on the model used.
Hope it will help you,
Rayan
19 June 2020 at 17:32 #16677pvillard
BlockedHi Rayan,
Thank you for your answer. I am not sure if I will use it but can you share it just in case?
Best,
Pierre
23 June 2020 at 22:50 #1670224 June 2020 at 15:57 #16715Rayan
BlockedHi,
Sorry for the wait 🙂
Here is the component:
#ifndef SOFA_COMPONENT_FORCEFIELD_CURRENTSTRAIN_INL #define SOFA_COMPONENT_FORCEFIELD_CURRENTSTRAIN_INL #include "CurrentStrain.h" #include <sofa/core/behavior/ForceField.inl> #include <SofaBaseTopology/TopologyData.inl> #include <sofa/defaulttype/VecTypes.h> #include <SofaBaseMechanics/MechanicalObject.h> #include <sofa/core/ObjectFactory.h> #include <fstream> // for reading the file #include <iostream> //for debugging #include <Eigen/Eigenvalues> #include <unsupported/Eigen/MatrixFunctions> namespace sofa { namespace component { namespace forcefield { using namespace sofa::defaulttype; using namespace sofa::component::topology; using namespace core::topology; template<class DataTypes> CurrentStrain<DataTypes>::CurrentStrain(sofa::core::behavior::MechanicalState<DataTypes> *ff) : core::objectmodel::BaseObject() , d_currentStrain(initData(&d_currentStrain,"CurrentStrain","Current strain per Element")) , l_mstate(initLink("mstate", "input MechanicalState"), ff) , m_topology(nullptr) , l_topology(initLink("topology", "link to the topology container")) { } template<class DataTypes> void CurrentStrain<DataTypes>::init() { std::cout<<"I'm in init()" <<std::endl; if (!l_mstate) { msg_error(this) << "currentStrain not linked to any mechanical state" ; return; } if (l_topology.empty()) { msg_info() << "link to Topology container should be set to ensure right behavior. First Topology found in current context will be used."; l_topology.set(this->getContext()->getMeshTopologyLink()); } m_topology = l_topology.get(); msg_info() << "Topology path used: '" << l_topology.getLinkedPath() << "'"; if (m_topology == nullptr) { msg_error() << "No topology component found at path: " << l_topology.getLinkedPath() << ", nor in current context: " << this->getContext()->name; sofa::core::objectmodel::BaseObject::d_componentstate.setValue(sofa::core::objectmodel::ComponentState::Invalid); return; } unsigned int i,j; unsigned int nbTetrahedra=m_topology->getNbTetrahedra(); typename DataTypes::Real volume; typename DataTypes::Coord point[4]; const typename DataTypes::VecCoord restPosition = l_mstate->read(sofa::core::ConstVecCoordId::restPosition())->getValue(); for(i=0; i<nbTetrahedra; i++ ) { m_shapeVector.resize(nbTetrahedra); const vector< Tetrahedron > &tetrahedronArray= m_topology->getTetrahedra(); ///describe the indices of the 4 tetrahedron vertices const Tetrahedron &t= tetrahedronArray[i]; for(j=0;j<4;++j) { point[j]=(restPosition)[t[j]]; } /// compute 6 times the rest volume volume=dot(cross(point[2]-point[0],point[3]-point[0]),point[1]-point[0]); for(j=0;j<4;++j) { if (!(j%2)) m_shapeVector[i].push_back(-cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/ volume); else m_shapeVector[i].push_back(cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/ volume); } } BaseData* pos_init = l_mstate->read(sofa::core::ConstVecCoordId::position())->getData(); c_callback.addInput(pos_init); c_callback.addCallback(std::bind(&CurrentStrain::strainplaceholder,this)); } template<class DataTypes> void CurrentStrain<DataTypes>::reinit() { this->init(); } template <class DataTypes> void CurrentStrain<DataTypes>::update() { } template <class DataTypes> void CurrentStrain<DataTypes>::strain(const core::MechanicalParams*, const DataVecCoord& d_x) { //std::cout<<"I'm in strain" <<std::endl; const VecCoord& x = d_x.getValue(); unsigned int i=0,j=0,k=0,l=0; unsigned int nbTetrahedra=m_topology->getNbTetrahedra(); assert(this->mstate); // ask why Coord dp[3],x0,sv; helper::vector<double>* currentStrain = d_currentStrain.beginEdit(); currentStrain->resize(nbTetrahedra); for(i=0; i<nbTetrahedra; i++ ) { Matrix3 deformationGradient; const Tetrahedron &ta= m_topology->getTetrahedron(i); x0=x[ta[0]]; dp[0]=x[ta[1]]-x0; sv=m_shapeVector[i][1]; //std::cout<<"m_shapeVector vaut : "<< m_shapeVector <<std::endl; for (k=0;k<3;++k) { for (l=0;l<3;++l) { deformationGradient[k][l]=dp[0][k]*sv[l]; } } for (j=1;j<3;++j) { dp[j]=x[ta[j+1]]-x0; sv=m_shapeVector[i][j+1]; for (k=0;k<3;++k) { for (l=0;l<3;++l) { deformationGradient[k][l]+=dp[j][k]*sv[l]; } } } Matrix3 Id,C; for (k=0;k<3;++k) { for (l=k;l<3;++l) { C(k,l) = (deformationGradient(0,k)*deformationGradient(0,l)+ deformationGradient(1,k)*deformationGradient(1,l)+ deformationGradient(2,k)*deformationGradient(2,l)); C(l,k) = C(k,l); } } //Create identity Id = Id*0; Id(0,0) = 1; Id(1,1) = 1; Id(2,2) = 1; Matrix3 E = 1/2.0*(C-Id); currentStrain->operator[](i) = sqrt(pow(E(0,0),2) + pow(E(1,1),2) + pow(E(2,2),2) + pow(E(1,2),2) + pow(E(0,1),2) + pow(E(0,2),2)); // std::cout<<"la force courante est: " << d_currentStrain<<std::endl; } d_currentStrain.endEdit(); } template<class DataTypes> void CurrentStrain<DataTypes>::strainplaceholder() { DataVecCoord d_pos; VecCoord &pos = *d_pos.beginEdit(); pos = this->l_mstate->read(core::ConstVecCoordId::position())->getValue(); d_pos.setValue(pos); strain(core::MechanicalParams::defaultInstance(), d_pos); d_pos.endEdit(); } } //namespace forcefield } //namespace component } //namespace sofa #endif
It’s quiet raw but don’t hesitate if you have questions 🙂
Regards,
Rayan
3 January 2022 at 07:04 #21210bani
BlockedHi Rayan,
Could you please also tell how to implement it using sofapython ?Thanks
Vani10 January 2022 at 17:18 #21350Hugo
Keymasterhey @vanisri
Rayan is not the area anymore, he was an intern.
You can get inspired from the C++ and translate it into python.All you need is to have an access to the topology and the mechanical object (mstate).
Then compute at init the shapeVectors.When you want it, compute the strain as it is done by the function
strainplaceholder
(itself calling thestrain()
) function in Rayan’s code.Best wishes,
Hugo
-
AuthorPosts
- You must be logged in to reply to this topic.