Bilinear Material is crashing with both Fraser and LeapFrog Solver about weldform HOT 8 CLOSED

$ git diff ba5aff4
diff --git a/DEVLOG.md b/DEVLOG.md
index e70def6..a538fab 100644
--- a/DEVLOG.md
+++ b/DEVLOG.md
@@ -463,4 +463,7 @@
                                 - Fixed crash in Johnson Cook damage calc
                                 - Prevented sigma asterisk for damage be NAN at sigma_eq = 0
                                 - Added damage count
-                                - Added damage output
\ No newline at end of file
+                                - Added damage output
+                                - Corrected Johnson Cook damage expression
+                                - Calculate fracture strain and stresses per particle and not per pair
+                                - Corrected Acceleration expression
\ No newline at end of file
diff --git a/Source/Damage.cpp b/Source/Damage.cpp
index 571456d..ba85671 100644
--- a/Source/Damage.cpp
+++ b/Source/Damage.cpp
@@ -20,12 +20,23 @@ inline void Domain::CalcDamage(){
        //USED IN JOHNSON COOK DAMAGE
        #pragma omp parallel for schedule(static) num_threads(Nproc)
        #ifdef __GNUC__
-       for (size_t i=0; i< solid_part_count ; i++)     //Like in Domain::Move
+       for (size_t i=0; i< solid_part_count ; i++)     { //Like in Domain::Move
        #else
-       for (int i=0; i<Particles.Size(); i++)//Like in Domain::Move
+       for (int i=0; i<solid_part_count; i++) {//Like in Domain::Move
        #endif
-               Particles[i]->CalculateEquivalentStress();
-
+               if (Particles[i]->delta_pl_strain>0.0) {
+                       Particles[i]->CalculateEquivalentStress(); //Set Particles[i]->Sigma_eq
+                       if (Particles[i]->Sigma_eq>1.0e-3) {
+                               sig_as = 1.0/3.0* (Particles[i]->Sigma(0,0)+Particles[i]->Sigma(1,1)+Particles[i]->Sigma(2,2))/Particles[i]->Sigma_eq; //Stress triaxiality sig_m / sig_eff
+                               // Particles[i]->eps_f = Particles[i]->mat->damage->CalcFractureStrain(Particles[i]->eff_strain_rate, sig_as, PP[i]->T);
+                               Particles[i]->dam_D += Particles[i]->delta_pl_strain/Particles[i]->mat->damage->CalcFractureStrain(Particles[i]->eff_strain_rate, sig_as, Particles[i]->T);
+                               // if (Particles[i]->dam_D > 0.0 && Particles[i]->pl_strain ==0.0)
+                                       // cout <<"dam " <<Particles[i]->dam_D<<", Pl Strain: "<<Particles[i]->pl_strain<<endl;
+                               if (Particles[i]->dam_D > 1.0) Particles[i]->dam_D = 1.0;
+                       }
+               }//if incremental pl_strain
+       }
+
        #pragma omp parallel for schedule (static) private (PP, sig_as, i) num_threads(Nproc)
        #ifdef __GNUC__
        for (size_t k=0; k<Nproc;k++)
@@ -97,21 +108,21 @@ inline void Domain::CalcDamage(){

       //FIRST WE IMPLEMENT JOHNSON COOK FAILURE AS ISLAM 2017
       if (dam_D[k][p] <1.0){ //OR dam_df0[][] == 0.0
-        for (i=0;i<2 ;i++){
-                                       //IN CASE OF NO DAMAGE INCLUDED; IT IS NOT CALCULATED
-                                       // cout << "sig_00"<<PP[i]->Sigma(0,0)<<endl;
-                                       if (PP[i]->Sigma_eq>1.0e-3)
-                                               sig_as = 1.0/3.0* (PP[i]->Sigma(0,0)+PP[i]->Sigma(1,1)+PP[i]->Sigma(2,2))/PP[i]->Sigma_eq; //Stress triaxiality sig_m / sig_eff
-                                       else
-                                               sig_as = 0.0;
-                                       // cout << "Fr strain "<<PP[i]->mat->damage->CalcFractureStrain(PP[i]->eff_strain_rate, sig_as, PP[i]->T)<<endl;
-                                       // cout << "str_rate, sig_as, ppi "<< PP[i]->eff_strain_rate<<", " <<sig_as<< ", "<<PP[i]->T<<endl;
-                                       omp_set_lock(&PP[i]->my_lock);
-                                       PP[i]->dam_D += PP[i]->delta_pl_strain/PP[i]->mat->damage->CalcFractureStrain(PP[i]->eff_strain_rate, sig_as, PP[i]->T);
-                                       omp_unset_lock(&PP[i]->my_lock);
-                                       //cout << "PP[i]->dam_D: "<<PP[i]->dam_D<<endl;
-          //dam_D[k][p] += ;
-        }
+        // for (i=0;i<2 ;i++){
+                                       // //IN CASE OF NO DAMAGE INCLUDED; IT IS NOT CALCULATED
+                                       // // cout << "sig_00"<<PP[i]->Sigma(0,0)<<endl;
+                                       // if (PP[i]->Sigma_eq>1.0e-3) {
+                                               // sig_as = 1.0/3.0* (PP[i]->Sigma(0,0)+PP[i]->Sigma(1,1)+PP[i]->Sigma(2,2))/PP[i]->Sigma_eq; //Stress triaxiality sig_m / sig_eff
+                                               // // cout << "Fr strain "<<PP[i]->mat->damage->CalcFractureStrain(PP[i]->eff_strain_rate, sig_as, PP[i]->T)<<endl;
+                                               // // cout << "str_rate, sig_as, ppi "<< PP[i]->eff_strain_rate<<", " <<sig_as<< ", "<<PP[i]->T<<endl;
+                                               // // omp_set_lock(&PP[i]->my_lock);
+                                               // // PP[i]->dam_D += PP[i]->delta_pl_strain/PP[i]->eps_f;
+                                               // // if (PP[i]->dam_D > 1.0) PP[i]->dam_D = 1.0;
+                                               // // omp_unset_lock(&PP[i]->my_lock);
+                                       // }
+                                       // //cout << "PP[i]->dam_D: "<<PP[i]->dam_D<<endl;
+          // //dam_D[k][p] += ;
+        // }
                                dam_D[k][p] = (PP[0]->dam_D + PP[1]->dam_D) /2.0;
       } else {
                                dam_D[k][p] = 1.0;
@@ -122,9 +133,9 @@ inline void Domain::CalcDamage(){
   }//for pair p
   } //proc k

-       for (size_t i=0; i< solid_part_count ; i++)     //Like in Domain::Move
-               if (Particles[i]->dam_D==1.0)
-                       cout << "DAMAGE 1!"<<endl;
+       // for (size_t i=0; i< solid_part_count ; i++)  //Like in Domain::Move
+               // if (Particles[i]->dam_D==1.0)
+                       // cout << "DAMAGE 1!"<<endl;

 }

diff --git a/Source/InteractionAlt.cpp b/Source/InteractionAlt.cpp
index 2b49a33..13e4493 100644
--- a/Source/InteractionAlt.cpp
+++ b/Source/InteractionAlt.cpp
@@ -13,7 +13,7 @@ inline void Domain::CalcAccel() {
   Particle *P1, *P2;
   double dam_f = 1.0; //if not damage

-  #pragma omp parallel for schedule (static) private (P1,P2) num_threads(Nproc)
+  #pragma omp parallel for schedule (static) private (P1,P2,dam_f) num_threads(Nproc)^M
        #ifdef __GNUC__
        for (size_t k=0; k<Nproc;k++)
        #else
@@ -190,10 +190,10 @@ inline void Domain::CalcAccel() {
                // Locking the particle 1 for updating the properties
                omp_set_lock(&P1->my_lock);
                        if (!gradKernelCorr){
-                               P1->a                                   += mj * temp;
+                               P1->a                                   += dam_f*mj * temp;^M
                                //P1->dDensity  += mj * (di/dj) * temp1;
                        } else{
-                               P1->a                                   += mj * temp_c[0];
+                               P1->a                                   += dam_f *mj * temp_c[0];^M
                                //P1->dDensity  += mj * (di/dj) * temp1_c[0];
                        }

@@ -202,9 +202,9 @@ inline void Domain::CalcAccel() {
                // Locking the particle 2 for updating the properties
                omp_set_lock(&P2->my_lock);
                        if (!gradKernelCorr){
-                               P2->a                                   -= mi * temp;
+                               P2->a                                   -= dam_f * mi * temp;                           ^M
                        }else {
-                               P2->a                                   -= mi * temp_c[1];
+                               P2->a                                   -= dam_f * mi * temp_c[1];^M
                        }
                omp_unset_lock(&P2->my_lock);
     //#endif
@@ -462,7 +462,7 @@ inline void Domain::RateTensorsReduction(){
 inline void Domain::CalcDensInc() {
        double dam_f = 1.0; //if not damage
   Particle *P1, *P2;
-       #pragma omp parallel for schedule (static) private (P1,P2) num_threads(Nproc)
+       #pragma omp parallel for schedule (static) private (P1,P2,dam_f) num_threads(Nproc)^M
        #ifdef __GNUC__
        for (size_t k=0; k<Nproc;k++)
        #else
diff --git a/Source/Material.h b/Source/Material.h
index 9676a1f..7ff8ba9 100644
--- a/Source/Material.h
+++ b/Source/Material.h
@@ -56,7 +56,7 @@ public DamageModel {
   }
        std::string getDamageCriterion() {return string("JohnsonCook");}
        double CalcFractureStrain(const double &str_rate, const double &sig_as,const double &T) { //sig_as is stress triaxiality
-    return ( (m_D1 + pow(m_D2, m_D3*sig_as)) * (pow (1.0 + (str_rate/m_eps0), m_D4)) ); //Islam 2017 eq. 35,
+    return ( (m_D1 + m_D2 *exp(m_D3*sig_as)) * pow (1.0 + (str_rate/m_eps0), m_D4) * (1.0 + m_D5 * T)); //Islam 2017 eq. 35, ^M
   }
        virtual ~JohnsonCookDamage(){}
 };
diff --git a/Source/Output.cpp b/Source/Output.cpp
index 0cce6a1..0f3921a 100644
--- a/Source/Output.cpp
+++ b/Source/Output.cpp
@@ -275,6 +275,12 @@ inline void Domain::WriteXDMF (char const * FileKey)

     dsname.Printf("ps_en");
     H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,ps_en);
+^M
+               if (model_damage){^M
+                       dsname.Printf("Damage");^M
+                       H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Damage);                    ^M
+                       ^M
+               }^M

     dims[0] = 3*Particles.Size();
        dsname.Printf("Displacement");
@@ -282,11 +288,7 @@ inline void Domain::WriteXDMF (char const * FileKey)
        dsname.Printf("Contact Force");
     H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,ContForce);

-               if (model_damage){
-    // dsname.Printf("Damage");
-    // H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Damage);
-
-               }
+^M

        // dsname.Printf("Tg Dir");
     // H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,TgDir);
@@ -478,11 +480,11 @@ inline void Domain::WriteXDMF (char const * FileKey)
     oss << "       </DataItem>\n";
     oss << "     </Attribute>\n";
                if (model_damage){
-    // oss << "     <Attribute Name=\"damage\" AttributeType=\"Scalar\" Center=\"Node\">\n";
-    // oss << "       <DataItem Dimensions=\"" << Particles.Size() << "\" NumberType=\"Float\" Precision=\"10\"  Format=\"HDF\">\n";
-    // oss << "        " << fn.CStr() <<":/damage \n";
-    // oss << "       </DataItem>\n";
-    // oss << "     </Attribute>\n";
+                       oss << "     <Attribute Name=\"Damage\" AttributeType=\"Scalar\" Center=\"Node\">\n";^M
+                       oss << "       <DataItem Dimensions=\"" << Particles.Size() << "\" NumberType=\"Float\" Precision=\"10\"  Format=\"HDF\">\n";^M
+                       oss << "        " << fn.CStr() <<":/Damage \n";^M
+                       oss << "       </DataItem>\n";^M
+                       oss << "     </Attribute>\n";                   ^M
                }
   // oss << "     <Attribute Name=\"Tg Dir\" AttributeType=\"Vector\" Center=\"Node\">\n";
     // oss << "       <DataItem Dimensions=\"" << Particles.Size() << " 3\" NumberType=\"Float\" Precision=\"10\" Format=\"HDF\">\n";
diff --git a/Source/SolverLeapfrog.cpp b/Source/SolverLeapfrog.cpp
index 0884f1b..eb939cf 100644
--- a/Source/SolverLeapfrog.cpp
+++ b/Source/SolverLeapfrog.cpp
@@ -395,7 +395,7 @@ inline void Domain::SolveDiffUpdateLeapFrog (double tf, double dt, double dtOut,
                                for (int p=0;p<Particles.Size();p++){
                                        if (Particles[p]->dam_D>max_dam)
                                                max_dam = Particles[p]->dam_D;
-                                       if (Particles[p]->dam_D>1.0)
+                                       if (Particles[p]->dam_D==1.0) ^M
                                                dam_count++;
                                }