[SOLVED] Extracting stiffness matrix from TriangleFEMForceField

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