[NEW PATCH] Minor Bug in the stress function

dG3D/src/dG3DMaterialLaw.cpp

@@ -5602,11 +5602,17 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   STensor33 dPdgradT; STensorOperation::zero(dPdgradT);
   // The tangent C_hom is combined with thermal tangents.
-  static fullMatrix<double> C_hom(_col_Di,_col_Di); if(C_hom.size1()!=_col_Di || C_hom.size2()!=_col_Di){C_hom.resize(_col_Di,_col_Di);}; C_hom.setAll(0.); 
-  static fullVector<double> PQ_hom(_col_Di); if(PQ_hom.size()!=_col_Di){PQ_hom.resize(_col_Di);}; PQ_hom.setAll(0.);
-  static fullVector<double> dPQdT_hom(_col_Di); if(dPQdT_hom.size()!=_col_Di){dPQdT_hom.resize(_col_Di);}; dPQdT_hom.setAll(0.);
+  fullMatrix<double> C_hom(_col_Di,_col_Di); 
+  // if(C_hom.size1()!=_col_Di || C_hom.size2()!=_col_Di){C_hom.resize(_col_Di,_col_Di);}; 
+  C_hom.setAll(0.); 
+  fullVector<double> PQ_hom(_col_Di); 
+  // if(PQ_hom.size()!=_col_Di){PQ_hom.resize(_col_Di);}; 
+  PQ_hom.setAll(0.);
+  fullVector<double> dPQdT_hom(_col_Di); 
+  // if(dPQdT_hom.size()!=_col_Di){dPQdT_hom.resize(_col_Di);}; 
+  dPQdT_hom.setAll(0.);
   double Cp_hom(0.), thermSrc_hom(0.), mechSource_hom(0.), dwdt_hom(0.), dmechSourcedt_hom(0.);
-  static fullVector<double> dwdf_hom(9), dmechSourcedf_hom(9); dwdf_hom.setAll(0.); dmechSourcedf_hom.setAll(0.);
+  fullVector<double> dwdf_hom(9), dmechSourcedf_hom(9); dwdf_hom.setAll(0.); dmechSourcedf_hom.setAll(0.);
   fullVector<double>& FHvec = ipvcur->getRefToCombinedGradientVect();
   fullVector<double>& PQvec = ipvcur->getRefToCombinedStressVect(); PQvec.setAll(0.0);
@@ -5618,13 +5624,27 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   dwdtvec.resize(_NLeaf); dwdtvec.setAll(0.); dmechSrcdTvec.resize(_NLeaf); dmechSrcdTvec.setAll(0.);
   dwdFvec.resize(9*_NLeaf); dwdFvec.setAll(0.); dmechSrcdFvec.resize(9*_NLeaf); dmechSrcdFvec.setAll(0.);
-  static fullVector<double> Res_node(_row_Di*_NInterface); if(Res_node.size()!=_row_Di*_NInterface){Res_node.resize(_row_Di*_NInterface);}; Res_node.setAll(0.);  
-  static fullVector<double> ab_step(_row_Di*_NInterface); if(ab_step.size()!=_row_Di*_NInterface){ab_step.resize(_row_Di*_NInterface);}; ab_step.setAll(0.0);
-  static fullMatrix<double> NTVC(_row_Di*_NInterface, _col_Di*_NLeaf); if(NTVC.size1()!=_row_Di*_NInterface || NTVC.size2()!=_col_Di*_NLeaf){NTVC.resize(_row_Di*_NInterface, _col_Di*_NLeaf);}; NTVC.setAll(0.0); 
-  static fullMatrix<double> Jacobian_inv(_row_Di*_NInterface, _row_Di*_NInterface); if(Jacobian_inv.size1()!=_row_Di*_NInterface || Jacobian_inv.size2()!=_row_Di*_NInterface){Jacobian_inv.resize(_row_Di*_NInterface, _row_Di*_NInterface);}; Jacobian_inv.setAll(0.0); 
-  static fullMatrix<double> tmp(_col_Di*_NLeaf, _col_Di); if(tmp.size1()!=_col_Di*_NLeaf || tmp.size2()!=_col_Di){tmp.resize(_col_Di*_NLeaf, _col_Di);}; tmp.setAll(0.0); 
-  static fullMatrix<double> dRdF(_row_Di*_NInterface, _col_Di); if(dRdF.size1()!=_row_Di*_NInterface || dRdF.size2()!=_col_Di){dRdF.resize(_row_Di*_NInterface, _col_Di);}; dRdF.setAll(0.0);
-  static fullMatrix<double> dadF(_row_Di*_NInterface, _col_Di); if(dadF.size1()!=_row_Di*_NInterface || dadF.size2()!=_col_Di){dadF.resize(_row_Di*_NInterface, _col_Di);}; dadF.setAll(0.0);
+  fullVector<double> Res_node(_row_Di*_NInterface); 
+  // if(Res_node.size()!=_row_Di*_NInterface){Res_node.resize(_row_Di*_NInterface);}; 
+  Res_node.setAll(0.);  
+  fullVector<double> ab_step(_row_Di*_NInterface); 
+  // if(ab_step.size()!=_row_Di*_NInterface){ab_step.resize(_row_Di*_NInterface);}; 
+  ab_step.setAll(0.0);
+  fullMatrix<double> NTVC(_row_Di*_NInterface, _col_Di*_NLeaf); 
+  // if(NTVC.size1()!=_row_Di*_NInterface || NTVC.size2()!=_col_Di*_NLeaf){NTVC.resize(_row_Di*_NInterface, _col_Di*_NLeaf);}; 
+  NTVC.setAll(0.0); 
+  fullMatrix<double> Jacobian_inv(_row_Di*_NInterface, _row_Di*_NInterface); 
+  // if(Jacobian_inv.size1()!=_row_Di*_NInterface || Jacobian_inv.size2()!=_row_Di*_NInterface){Jacobian_inv.resize(_row_Di*_NInterface, _row_Di*_NInterface);}; 
+  Jacobian_inv.setAll(0.0); 
+  fullMatrix<double> tmp(_col_Di*_NLeaf, _col_Di); 
+  // if(tmp.size1()!=_col_Di*_NLeaf || tmp.size2()!=_col_Di){tmp.resize(_col_Di*_NLeaf, _col_Di);}; 
+  tmp.setAll(0.0); 
+  fullMatrix<double> dRdF(_row_Di*_NInterface, _col_Di); 
+  // if(dRdF.size1()!=_row_Di*_NInterface || dRdF.size2()!=_col_Di){dRdF.resize(_row_Di*_NInterface, _col_Di);}; 
+  dRdF.setAll(0.0);
+  fullMatrix<double> dadF(_row_Di*_NInterface, _col_Di); 
+  // if(dadF.size1()!=_row_Di*_NInterface || dadF.size2()!=_col_Di){dadF.resize(_row_Di*_NInterface, _col_Di);}; 
+  dadF.setAll(0.0);
     
   _C.setAll(0.0);
   int ite = 0;
@@ -5639,7 +5659,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
     _Nv.mult(ab_cur, FHvec);
     for(int i=0; i<_NLeaf; i++){
-      Msg::Error("i is %d", i);
+      // Msg::Error("i is %d", i);
       
       if (_VA[pos_start+i] > 1.e-6){
         ThermoMechanicsDG3DIPVariableBase* ipv_i = dynamic_cast<ThermoMechanicsDG3DIPVariableBase*>(ipvcur->getIPv(i));
@@ -5799,7 +5819,6 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
   Jacobian.scale(-1.0);
   int pos0 = _Child[0][0];
   int pos1 = _Child[0][1];  
-
   if(r>_tol){
     Msg::Error("TMDMN residual (= %lf) didn't converge to tolenerce after %d iteration",r,ite);
     return;
@@ -5874,7 +5893,7 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
         for(int m=0;m<9;m++)
           for(int n=0;n<9;n++){
             dwdf_hom(n) += temp*dwdFvec(m+9*i)*tmp(9*i,n);
-            dmechSourcedf_hom(n) += temp*dmechSrcdTvec(m+9*i)*tmp(9*i,n);
+            dmechSourcedf_hom(n) += temp*dmechSrcdFvec(m+9*i)*tmp(9*i,n);
           }
       }
@@ -5919,7 +5938,10 @@ void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev,
       }
     }
   }
+  Msg::Error("Sizes of PQVec are %d", PQvec.size());
   ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+  Msg::Error("Sizes");
+  Msg::Error("Sizes of PQVec are %d", PQvec.size());
 } 
 void StochTMDMNDG3DMaterialLaw::setTime(const double t,const double dtime){