Home › Forum › SOFA › Using SOFA › [SOLVED] Extracting stiffness matrix from TriangleFEMForceField
Tagged: 64_bits, Linux_ubuntu, SOFA_2006, stiffness matrix
- This topic has 2 replies, 2 voices, and was last updated 3 years, 11 months ago by nhnhan.
-
AuthorPosts
-
13 December 2020 at 03:42 #18024nhnhanBlocked
Dear,
I have been trying to extract the stiffness matrices of the triangle element by going into the TriangleFEMForceField.inl. As far as I know, the exact value of stiffness matrices for each time step can be calculated by using equations: J*K*Jt (without considering rotation motion of element) or R*J*K*Jt*Rt (taking rotation motion into account). I was successful to display the value of J*K*Jt on the terminal, however, when I tried to print these matrices (for example, my model has 4 triangle element) into a text file (using code below), only the last one appeared on the text file. It seems that the other three matrices were replaced right after they were written in the file. Since I’m very poor on C++ programming skills so could anyone have a better way to extract such information? Thank you very much for any suggestion or comment.
Code:template <class DataTypes> void TriangleFEMForceField<DataTypes>::computeForce( Displacement &F, const Displacement &Depl, const MaterialStiffness &K, const StrainDisplacement &J ) { defaulttype::Mat<3,6,Real> Jt; Jt.transpose( J ); // JtD = Jt * Depl; JtD[0] = Jt[0][0] * Depl[0] + /* Jt[0][1] * Depl[1] + */ Jt[0][2] * Depl[2] /* + Jt[0][3] * Depl[3] + Jt[0][4] * Depl[4] + Jt[0][5] * Depl[5] */ ; JtD[1] = /* Jt[1][0] * Depl[0] + */ Jt[1][1] * Depl[1] + /* Jt[1][2] * Depl[2] + */ Jt[1][3] * Depl[3] + /* Jt[1][4] * Depl[4] + */ Jt[1][5] * Depl[5]; JtD[2] = Jt[2][0] * Depl[0] + Jt[2][1] * Depl[1] + Jt[2][2] * Depl[2] + Jt[2][3] * Depl[3] + Jt[2][4] * Depl[4] /* + Jt[2][5] * Depl[5] */ ; defaulttype::Vec<3,Real> KJtD; // KJtD = K * JtD; KJtD[0] = K[0][0] * JtD[0] + K[0][1] * JtD[1] /* + K[0][2] * JtD[2] */; KJtD[1] = K[1][0] * JtD[0] + K[1][1] * JtD[1] /* + K[1][2] * JtD[2] */; KJtD[2] = /* K[2][0] * JtD[0] + K[2][1] * JtD[1] */ + K[2][2] * JtD[2]; // F = J * KJtD; F[0] = J[0][0] * KJtD[0] + /* J[0][1] * KJtD[1] + */ J[0][2] * KJtD[2]; F[1] = /* J[1][0] * KJtD[0] + */ J[1][1] * KJtD[1] + J[1][2] * KJtD[2]; F[2] = J[2][0] * KJtD[0] + /* J[2][1] * KJtD[1] + */ J[2][2] * KJtD[2]; F[3] = /* J[3][0] * KJtD[0] + */ J[3][1] * KJtD[1] + J[3][2] * KJtD[2]; F[4] = /* J[4][0] * KJtD[0] + J[4][1] * KJtD[1] + */ J[4][2] * KJtD[2]; F[5] = /* J[5][0] * KJtD[0] + */ J[5][1] * KJtD[1] /* + J[5][2] * KJtD[2] */ ; // Calculation of Stiffness matrice defaulttype::MatNoInit<6, 6, Real> JKJt; JKJt = J*K*Jt; // in-plane stiffness matrix, 6x6 std::cout <<"Value = " <<JKJt <<std::endl; std::ofstream myfile; myfile.open("/home/nhnhanbk/Desktop/Sofa/sofa/nhnhan/IoTouch/example.txt"); using namespace std::chrono; myfile.flush(); for (int i = 0; i < 6;i++) { auto now = std::chrono::system_clock::now(); auto ms = std::chrono::duration_cast<milliseconds>(now.time_since_epoch()) % 1000; auto timer = system_clock::to_time_t(now); std::tm bt = *std::localtime(&timer); std::ostringstream oss; oss << std::put_time(&bt, "%H:%M:%S"); // HH:MM:SS oss << '.' << std::setfill('0') << std::setw(3) << ms.count(); std::cout << "Time" << " "+oss.str() << std::endl; myfile <<JKJt[i] <<","<< " "+oss.str()<< "\n"; } myfile.close();
19 December 2020 at 00:35 #18064HugoKeymasterHey @nhnhan
would this not be due to the “open” options ?
Have you tried:outfile.open("example.txt", std::ios_base::app); // append instead of overwrite
Best
Hugo
19 December 2020 at 14:58 #18070 -
AuthorPosts
- You must be logged in to reply to this topic.