diff --git a/NonLinearSolver/internalPoints/ipVEVPMFH.cpp b/NonLinearSolver/internalPoints/ipVEVPMFH.cpp
index a353b569e774e3d2e015767dababce1fa3b0872b..8887d1596c3e3f9d563e231cea8eac7362adab8f 100644
--- a/NonLinearSolver/internalPoints/ipVEVPMFH.cpp
+++ b/NonLinearSolver/internalPoints/ipVEVPMFH.cpp
@@ -13,64 +13,68 @@
 #include "fonctions_pratiques.h"
 
 
-double IPVEVPMFH::defoEnergy() const
-{
+
+double IPVEVPMFH::defoEnergy() const{
   return _elasticEne;
 }
-double IPVEVPMFH::plasticEnergy() const
-{
+
+
+
+double IPVEVPMFH::plasticEnergy() const{
   return _plasticEne;
 }
 
 
-double IPVEVPMFH::get(int comp) const
-{
-
 
-  if(comp == IPField::SIG_XX) // Cauchy stress field
-	  return _state->strs[0];
+double IPVEVPMFH::get(int comp) const{
+  
+  if(comp == IPField::SIG_XX)			// Cauchy stress field 
+	  return _state.strs[0];
   else if(comp == IPField::SIG_YY)
-	  return _state->strs[1];
+	  return _state.strs[1];
   else if(comp == IPField::SIG_ZZ)
-	  return _state->strs[2];	  
+	  return _state.strs[2];	  
   else if(comp==IPField::SIG_XY) 
-	  return (_state->strs[3] / sqrt(2));
+	  return (_state.strs[3] / sqrt(2));
   else if(comp==IPField::SIG_XZ) 
-	  return (_state->strs[5] / sqrt(2));
+	  return (_state.strs[5] / sqrt(2));
   else if(comp==IPField::SIG_YZ) 
-	  return (_state->strs[4] / sqrt(2));
-  else if(comp == IPField::SVM) // von Mises
-  {
-     double svm= VEVPMFH_FP::vonmises(_state->strs);
+	  return (_state.strs[4] / sqrt(2));
+  
+  else if(comp == IPField::SVM){		// von Mises
+     double svm= VEVPMFH_FP::vonmises(_state.strs);
      return svm;
   }
 
-  else if(comp == IPField::PLASTICSTRAIN) // Accumulated plasticity 
-     return _state->p;
-  else if(comp == IPField::STRAIN_XX) // True strain field
-	  return _state->strn[0];
+  else if(comp == IPField::PLASTICSTRAIN) 	// Accumulated plasticity 
+     return _state.p;
+  else if(comp == IPField::STRAIN_XX) 	// True strain field
+	  return _state.strn[0];
   else if(comp == IPField::STRAIN_YY) 
-	  return _state->strn[1];
+	  return _state.strn[1];
   else if(comp == IPField::STRAIN_ZZ) 
-	  return _state->strn[2];
+	  return _state.strn[2];
   else if(comp == IPField::STRAIN_XY)
-	  return _state->strn[3]/sqrt(2);
- else if(comp == IPField::STRAIN_XZ)
-	  return _state->strn[5]/sqrt(2);
- else if(comp == IPField::STRAIN_XZ)
-	  return _state->strn[4]/sqrt(2);
+	  return _state.strn[3]/sqrt(2);
+  else if(comp == IPField::STRAIN_XZ)
+	  return _state.strn[5]/sqrt(2);
+  else if(comp == IPField::STRAIN_XZ)
+	  return _state.strn[4]/sqrt(2);
   else
      return 0.;
 }
 
-void IPVEVPMFH::restart()
-{
+
+
+
+void IPVEVPMFH::restart(){
+  
   IPVariableMechanics::restart();
   
-  restartManager::restart(_state->p); 
-  restartManager::restart(_state->strs,6);
-  restartManager::restart(_state->strn,6);
-  restartManager::restart(_state->pstrn,6);
+  restartManager::restart(_state.p); 
+  restartManager::restart(_state.strs,6);
+  restartManager::restart(_state.strn,6);
+  restartManager::restart(_state.pstrn,6);
   
 
   restartManager::restart(_elasticEne);
@@ -79,7 +83,6 @@ void IPVEVPMFH::restart()
   restartManager::restart(_irreversibleEnergy);
   restartManager::restart(_DirreversibleEnergyDF);
   
-
   return;
 }
 
diff --git a/NonLinearSolver/internalPoints/ipVEVPMFH.h b/NonLinearSolver/internalPoints/ipVEVPMFH.h
index bf4d043d9adcab670262e14bfc8bb3ae466c1478..5c615a4642bacc31ca966ea67e9dc8e8c0860dfa 100644
--- a/NonLinearSolver/internalPoints/ipVEVPMFH.h
+++ b/NonLinearSolver/internalPoints/ipVEVPMFH.h
@@ -15,82 +15,76 @@
 #include "fullMatrix.h"
 #include "fonctions_pratiques.h"
 
-class IPVEVPMFH : public IPVariableMechanics
-{
+class IPVEVPMFH : public IPVariableMechanics{
 
- public: // All data public to avoid the creation of function to access
-
-
-
- // state variables
-  VEVPMFH_FP::STATE* _state;
-  double _elasticEne;
-  double _plasticEne;
+  public: // All data public to avoid the creation of function to access
+  // state variables
+    VEVPMFH_FP::STATE _state;
+    double _elasticEne;
+    double _plasticEne;
   
 
- protected:
-// parthFollowing
-  double _irreversibleEnergy;
-  STensor3 _DirreversibleEnergyDF;    // Derivatives in terms of deformation gradient
-
-  
+  protected:
+  // parthFollowing
+    double _irreversibleEnergy;
+    STensor3 _DirreversibleEnergyDF;    // Derivatives in terms of deformation gradient
 
-
- public:
-  IPVEVPMFH() : IPVariableMechanics(), _elasticEne(0.),_plasticEne(0.),
-                            _irreversibleEnergy(0),  _DirreversibleEnergyDF(0.)
-    { 
+  public:
+    IPVEVPMFH() : IPVariableMechanics(), _elasticEne(0.),_plasticEne(0.), _irreversibleEnergy(0),  _DirreversibleEnergyDF(0.){ 
       //create and initialize
-       VEVPMFH_FP::mallocstate(_state);
+       VEVPMFH_FP::mallocstate(&_state); 
     }
 
-  IPVEVPMFH(const IPVEVPMFH &source) : IPVariableMechanics(source)
-    {
-         //create and copy
-        VEVPMFH_FP::mallocstate(_state);
-        VEVPMFH_FP::copystate(source._state, _state);
+
+    IPVEVPMFH(const IPVEVPMFH &source) : IPVariableMechanics(source){
+        //create and copy
+        VEVPMFH_FP::mallocstate(&_state);
+        VEVPMFH_FP::copystate(&(source._state), &_state);
 	_elasticEne=source._elasticEne;
  	_plasticEne=source._plasticEne;
         _irreversibleEnergy = source._irreversibleEnergy;
-        _DirreversibleEnergyDF = source._DirreversibleEnergyDF;  
-
+        _DirreversibleEnergyDF = source._DirreversibleEnergyDF;
     }
 
-  IPVEVPMFH &operator = (const IPVariable &_source)
-  {
+
+    IPVEVPMFH &operator = (const IPVariable &_source){
       IPVariableMechanics::operator=(_source);
-      const IPVEVPMFH *source=static_cast<const IPVEVPMFH *> (&_source);
+      const IPVEVPMFH* source=static_cast<const IPVEVPMFH *> (&_source);
+      
+      //copy 
+      VEVPMFH_FP::copystate(&(source->_state), &_state);
       _elasticEne=source->_elasticEne;
       _plasticEne=source->_plasticEne;
-
-      //copy
-      VEVPMFH_FP::copystate(source->_state, _state);
       _irreversibleEnergy = source->_irreversibleEnergy;
-      _DirreversibleEnergyDF = source->_DirreversibleEnergyDF;  
-
-      return *this;
-  }
-  virtual ~IPVEVPMFH()
-  {
-      //delete le state
-      VEVPMFH_FP::freestate(*_state);
-  }
-  virtual double get(const int i) const;
-  virtual double defoEnergy() const;
-  virtual double plasticEnergy() const;
+      _DirreversibleEnergyDF = source->_DirreversibleEnergyDF; 
+      
+      return *this; 
+    }
+  
+  
+  
+    virtual ~IPVEVPMFH(){
+      //Destructor 
+      VEVPMFH_FP::freestate(_state);
+    }
+  
+  
+    virtual double get(const int i) const;
+    virtual double defoEnergy() const;
+    virtual double plasticEnergy() const;
 
 
-  virtual IPVariable* clone() const {return new IPVEVPMFH(*this);};
-  virtual void restart();
+    virtual IPVariable* clone() const {return new IPVEVPMFH(*this);};
+    virtual void restart();
 
-  // for path following
-  virtual double irreversibleEnergy() const {return _irreversibleEnergy;};
-  virtual double& getRefToIrreversibleEnergy() {return _irreversibleEnergy;};  
-  virtual const STensor3& getConstRefToDIrreversibleEnergyDF() const{return _DirreversibleEnergyDF;};
-  virtual STensor3& getRefToDIrreversibleEnergyDF() {return _DirreversibleEnergyDF;};
+    // for path following
+    virtual double irreversibleEnergy() const {return _irreversibleEnergy;};
+    virtual double& getRefToIrreversibleEnergy() {return _irreversibleEnergy;};  
+    virtual const STensor3& getConstRefToDIrreversibleEnergyDF() const{return _DirreversibleEnergyDF;};
+    virtual STensor3& getRefToDIrreversibleEnergyDF() {return _DirreversibleEnergyDF;};
 
-  virtual VEVPMFH_FP::STATE* getSTATE() {return _state;}
-  virtual const VEVPMFH_FP::STATE* getSTATE() const {return _state;}
+    virtual VEVPMFH_FP::STATE* getSTATE() {return &_state;}
+    virtual const VEVPMFH_FP::STATE* getSTATE() const {return &_state;}
   
   
 };
diff --git a/NonLinearSolver/materialLaw/box_secant.cpp b/NonLinearSolver/materialLaw/box_secant.cpp
index 08c10fd76869ae917cea1f980ebc0a939869a056..19ec142a1ce6fb1c776a5cfbb2add41ff7a1206f 100644
--- a/NonLinearSolver/materialLaw/box_secant.cpp
+++ b/NonLinearSolver/materialLaw/box_secant.cpp
@@ -6,34 +6,23 @@
 #include <stdio.h>
 #include <string.h>
 
-//#include "global.h"
-
-//#include "matprop.h"
-//#include "state.h"
-
 #include "box_secant.h"
-//#include "fonctions_secant.h"
-//#include "fonctions_pratiques.h"
-using namespace VEVPMFH_FP;
-
 
-/**********************************************************************************************************************/
-/* Resolution of a time step of an viscoelastic-viscoplastic ISOTROPIC phase given the residual strain increment (J2) */
-/**********************************************************************************************************************/
-/** November 2020 - M.HADDAD  **/
 
-/* INPUT: material : material properties with a MATPROP struct
-          dt: time step
-          STATE tn : State at time tn
-          STATE tn_res : Residual state at time tn 
-          dstrn_r: Residual strain increment on the time step */
-
-/* OUTPUT: strs: stress at time t=tau
-           p: accumulated plasticity at time t=tau 
-           pstrntn: plastic strain tensor */
+//Description: Resolution of a time step of an viscoelastic-viscoplastic ISOTROPIC phase given the residual strain increment (J2) s(see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
+//
+//	** INPUT: 	material : material properties with a MATPROP struct
+//         	 	dt: time step
+//			STATE tn : State at time tn with a STATE struct
+//		        dstrn: strain increment on the time step
+//
+//	** OUTPUT: 	STATE tn+1 : State at time tn with a STATE struct
+//           		Update strs, p and pstrn 
            
-int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(MATPROP *material, STATE *tn, double dt, double *dstrn, STATE *tau, double **Calgo)
-{
+int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(const VEVPMFH_FP::MATPROP *material, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn, STATE *tau, double **Calgo){
 
   /**	Declaration	**/
   /***********************/
@@ -49,7 +38,7 @@ int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(MATPROP *material, STATE
   double R, dRdp;
   double *S_pred;
   double *N_pred;
-  double strseq, p, G_tld;
+  double strseq_tau, ptau=0, G_tld;
   double *strstau;
   double *pstrntau; 
   double **tens4;
@@ -107,6 +96,7 @@ int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(MATPROP *material, STATE
   /** yield function  **/
   f_pred = strseq_pred - R - (material->yldstrs);
   
+  
   /** VE Predictor **/
   /******************/
   if (f_pred <= 0){
@@ -118,68 +108,89 @@ int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(MATPROP *material, STATE
   
   /** VP Corrector  **/
   /*******************/
-  else{ 
-    
-    dev(strs_pred, S_pred);
-    addtens2((3/(2*strseq_pred)), S_pred, 0, S_pred, N_pred);
-    
-    strseq = strseq_pred;
-    p = tn->p;
+  else{
+  
+    /* compute N_pred */
+     dev(strs_pred, S_pred);
+     addtens2((3/(2*strseq_pred)), S_pred, 0, S_pred, N_pred);
+
+
+    /*Initial values of p and equivelent stress*/
+    ptau=tn->p;
+    strseq_tau = strseq_pred;
+
+    /*compute G_tilde (dt)*/
     G_tld = compute_Gtilde(material,dt);
-    computevisc (strseq, p, material, &gv, &dgvdstrseq, &dgvdp);
-    
-    /** residuals **/
-    kp = p - (tn->p) - gv*dt;
-    kstrs = strseq + 3*G_tld*(p-(tn->p)) - strseq_pred;
-    
-    res = max(fabs(kp),fabs(kstrs/(material->yldstrs)));
-    
+
+    computevisc (strseq_tau, ptau, material, &gv, &dgvdstrseq, &dgvdp);
+  
+  
+    /*Initial values of kp and k_sigma*/
+    kp = ptau - (tn->p) - (gv * dt);
+    kstrs = strseq_tau + (3 * G_tld * (ptau-(tn->p))) - strseq_pred;
+
+    res = max(fabs(kp),fabs(kstrs/(material->E0)));
+
+
+
     for (m=1; res > tol; m++) { 
+    
+      /* Compute correction on p = cp */
+      /*cp =  - (kp + dgvdstrseq*dt*kstrs)/(1+(3*G_tld+dRdp)*dgvdstrseq*dt);*/
       cp =  - (kp + dgvdstrseq*dt*kstrs)/(1-(dgvdp*dt)+ (3*G_tld*dgvdstrseq*dt));
-      p += cp;
-      strseq += -kstrs-3*G_tld*cp;
-      
-      
-      compute_R (p, material, &R, &dRdp);
-      computevisc (strseq, p, material, &gv, &dgvdstrseq, &dgvdp);
-      
-      kp = p - (tn->p) - gv*dt;
-      kstrs = strseq + 3*G_tld*(p-(tn->p)) - strseq_pred;
-      res = max(fabs(kp),fabs(kstrs/(material->yldstrs)));
-      
+
+
+      /*update accumulated plastic strain p and eq_stress */
+      ptau += cp;
+      strseq_tau += -kstrs-(3 * G_tld * cp);
+
+      compute_R (ptau, material, &R, &dRdp);
+      computevisc (strseq_tau, ptau, material, &gv, &dgvdstrseq, &dgvdp);
+    
+
+      /* update kp and k_strs for the next iteration if res still > tol*/
+      kp = ptau - (tn->p) - (gv * dt);
+      kstrs = strseq_tau + (3 * G_tld * (ptau-(tn->p))) - strseq_pred;
+      res = max(fabs(kp),fabs(kstrs/(material->E0)));
+
+
       if (m>100) {
         /*printf("\nWARNING: Problems of convergence in constbox_VEVP !\n");
         printf("dstrn[0]=%g; ptn:%g; strstn[0]: %g; dt:%g; G_tld:%g; strseqtr:%g\n", dstrn[0], tn->p, (tn->strs)[0], dt, G_tld, strseq_tr); 
         printf("p:%g; dRdp:%g; gv:%g; dgvdstrseq:%g; dgvdp:%g\n", *p, dRdp, gv, dgvdstrseq, dgvdp);
-        printf("ftr: %g; strseq:%g; material->yldstrs: %g; R: %g; f: %g\n", ftr, strseq, material->yldstrs, R, strseq - material->yldstrs-R);*/ 
+        printf("ftr: %g; strseq:%g; material->yldstrs: %g; R: %g; f: %g\n", ftr, strseq, material->yldstrs, R, strseq - material->yldstrs-R); */
         printf("fabs(kp):%g;  fabs(kstrs/(material->yldstrs))=%g; res:%g; tol:%g\n", fabs(kp),fabs(kstrs/(material->yldstrs)), res, tol); 
         return (0);
-      } 
-    }
+      }    
+   }
   
-    addtens2(((2*strseq)/3), N_pred, 0, N_pred, Stau);
-    trace_strs_pred = trace(strs_pred);
-  
-    /** Update the Cauchy stress **/
+ 
+
+    addtens2((2.0*strseq_tau)/3.0, N_pred, 0, N_pred, Stau);
+    trace_strs_pred = trace(strs_pred);     
+   
+   
+    /* Update stress*/
     for (i=0;i<3;i++) {
-      strstau[i] = Stau[i] + (1/3)*trace_strs_pred;
-      strstau[i+3] = Stau[i+3];
+      (tau->strs)[i] = Stau[i] + (1.0/3.0)*trace_strs_pred;
+      (tau->strs)[i+3] = Stau[i+3];
     }
-    copyvect(strstau, tau->strs, 6);
-  
-    /** Update the plastic strain **/
-    tau->p = p;
-    dp = p - tn->p;
-    addtens2(dp, N_pred, 1, tn->pstrn, pstrntau);
-    copyvect(pstrntau, tau->pstrn, 6);
-  
-    /** Update Calgo  **/
+
+    /* Update p and plastic strain*/
+    tau->p = ptau;
+    addtens2((ptau-(tn->p)), N_pred, 1, tn->pstrn, tau->pstrn);
+    
+    
+    /* Update Calgo **/
+    computevisc (strseq_tau, ptau, material, &gv, &dgvdstrseq, &dgvdp);
+    dp = ptau - (tn->p);
+    
     hv = (1 / (dt * dgvdstrseq)) + (3*G_tld) - (dgvdp/dgvdstrseq);
     produit_tensoriel(N_pred, N_pred, NcrossN);
     compute_Idev(Idev);
-    addtens4((3/(2*strseq)), Idev, (-1/strseq), NcrossN, dNdstrs);  
+    addtens4((3/(2*strseq_tau)), Idev, (-1/strseq_tau), NcrossN, dNdstrs);  
     addtens4(1, E_tilde, -(pow((2*G_tld),2)/hv), NcrossN, tens4);
-    addtens4(1, tens4, -(pow((2*G_tld),2)*((strseq*dp)/(strseq + (3*G_tld*dp)))), dNdstrs, Calgo);
+    addtens4(1, tens4, -(pow((2*G_tld),2)*((strseq_tau*dp)/(strseq_tau + (3*G_tld*dp)))), dNdstrs, Calgo);
   }
   
 
diff --git a/NonLinearSolver/materialLaw/box_secant.h b/NonLinearSolver/materialLaw/box_secant.h
index 0137dda0ada94999a59af125d6a5605470b0bee9..c788822d1eda62e408a4509ebdbef6ad470a442a 100644
--- a/NonLinearSolver/materialLaw/box_secant.h
+++ b/NonLinearSolver/materialLaw/box_secant.h
@@ -1,12 +1,11 @@
 #ifndef BOX_SECANT_H
 #define BOX_SECANT_H 1
 
-//#include "matprop.h"
-//#include "state.h"
-
 #include "fonctions_pratiques.h"
+
 using namespace VEVPMFH_FP;
-/* Resolution of the constitutive box in viscoelasto-viscoplasticity with the secant formulation */
-int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(MATPROP *material, STATE *tn, double dt, double *dstrn, STATE *tau, double **Calgo);
+
+
+int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(const VEVPMFH_FP::MATPROP *material, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn, STATE *tau, double **Calgo);
 
 #endif
diff --git a/NonLinearSolver/materialLaw/fonctions_pratiques.cpp b/NonLinearSolver/materialLaw/fonctions_pratiques.cpp
index f925ce1d66c929140256c0f54dd5d52c7fac5567..d99a2d05d56fd98bf76418f32abeb125dd468b92 100644
--- a/NonLinearSolver/materialLaw/fonctions_pratiques.cpp
+++ b/NonLinearSolver/materialLaw/fonctions_pratiques.cpp
@@ -4,59 +4,76 @@
 #include <math.h>
 #include <stdlib.h>
 #include <stdio.h>
-#include <cstring>
+#include <string.h>
+
+
 #include "fonctions_pratiques.h"
 #include "box_secant.h"
 
-//using namespace VEVPMFH_FP;
 
-/********************************************************
-// ** CONVERTIT UN TENSEUR D'ORDRE 2 EN UN VECTEUR 6*1 **
-// ******************************************************/
 
-void VEVPMFH_FP::convertmatrvect (double **matr, double *vect)
-{
+
+
+
+//Description: Convert a matrix(3*3) to a vector(6*1) (Issam Doghri's convention)
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::convertmatrvect (double **matr, double *vect){
+  
+  /** Declaration **/
   int i;
 
-  for (i=0; i<3; i++) {
-    vect[i] = matr[i][i];
-  }
+  /**** Start ****/
+  for (i=0; i<3; i++) {vect[i] = matr[i][i];}
+  
   vect[3]=sqrt(2)*matr[0][1];
   vect[4]=sqrt(2)*matr[1][2];
   vect[5]=sqrt(2)*matr[0][2];
+  
 }
 
-/***************************************************************
-// ** CONVERTIT UN TENSEUR SOUS FORME 6*1 EN UNE MATRICE 3*3 **
-// *************************************************************/
 
-void VEVPMFH_FP::convertvectmatr (double *vect, double **matr)
-{
-  int i;
 
-  for (i=0; i<3; i++) {
-    matr[i][i] = vect[i];
-  }
+//Description: Convert a vector(6*1) to a matrix (3*3) (Issam Doghri's convention)
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::convertvectmatr (double *vect, double **matr){
+  
+  /** Declaration **/
+  int i;
+  
+  /**** Start ****/
+  for (i=0; i<3; i++) {matr[i][i] = vect[i];}
+  
   matr[0][1] = vect[3]/sqrt(2);
   matr[1][0] = matr[0][1];
   matr[1][2] = vect[4]/sqrt(2);
   matr[2][1] = matr[1][2];
   matr[0][2] = vect[5]/sqrt(2);
   matr[2][0] = matr[0][2];
+  
 }
 
-/* ********************************************
-// ** SUM OF TWO TENSORS OF THE FOURTH ORDER **
-// ** EACH ONE CAN BE PREVIOUSLY BY A SCALAR **
-// ********************************************/
 
-void VEVPMFH_FP::addtens4 (double a, double **A, double b, double **B, double **res)
-{
+
+
+//Description: Sum of two 4th order tensors (stored as matrices)
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::addtens4 (double a, double **A, double b, double **B, double **res){
+  
+  /** Declaration **/
   int i,j;
 
+  /**** Start ****/
   if (b!=0) {      /* b != 0 */
- 
-   for (i=0; i<6;i++) {
+    for (i=0; i<6;i++) {
       for (j=0; j<6; j++) {
 	res[i][j] = a*A[i][j] + b*B[i][j];
       }
@@ -64,7 +81,6 @@ void VEVPMFH_FP::addtens4 (double a, double **A, double b, double **B, double **
   }
 
   else {           /* b == 0 */
-
     for (i=0; i<6;i++) {
       for (j=0; j<6; j++) {
 	res[i][j] = a*A[i][j];
@@ -75,15 +91,19 @@ void VEVPMFH_FP::addtens4 (double a, double **A, double b, double **B, double **
 }
 
 
-/* ********************************************
-// ** SUM OF TWO TENSORS OF THE SECOND ORDER **
-// ** EACH ONE CAN BE PREVIOUSLY BY A SCALAR **
-// ********************************************/
 
-void VEVPMFH_FP::addtens2 (double a, double *A, double b, double *B, double *res)
-{
+
+//Description: Sum of two 2nd order tensors (stored as vectors)
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::addtens2 (double a, double *A, double b, double *B, double *res){
+  
+  /** Declaration **/
   int i;
 
+  /**** Start ****/
   if (b!=0) {                /* b!=0 */
     for (i=0; i<6; i++) {
       res[i] = a*A[i] + b*B[i];
@@ -99,37 +119,38 @@ void VEVPMFH_FP::addtens2 (double a, double *A, double b, double *B, double *res
 }
 
 
-/* ************************************************
-// ** CONTRACTION DE DEUX TENSEURS D'ORDRE 2 MIS **
-// ** SOUS FORME VECTORIELLE EN UN SCALAIRE      **
-// ***********************************************/
 
-double VEVPMFH_FP::contraction22 (const double *TA, const double *TB)
+//Description: Double contraction of two 2nd order tensors (stored as vectors)
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-{
+double VEVPMFH_FP::contraction22 (const double *TA, const double *TB){
 
+  /** Declaration **/
   int i;
   double res = 0;
   
-  for(i=0; i<6; i++) {
-    res = res + TA[i]*TB[i];
-  };
-  
+  /**** Start ****/
+  for(i=0; i<6; i++) {res = res + TA[i]*TB[i];}
+
   return res;
+
 }
     
-/* *************************************************************      
-// ** CONTRACTION D'UN TENSEUR D'ORDRE 4 (mis sous forme 6*6) **
-// **      ET D'UN TENSEUR D'ORDRE 2 (sous forme 6*1)         **
-// **  POUR DONNER UN TENSEUR D'ORDRE 2 (mis sous forme 6*1)  **
-// ************************************************************/
-
-void VEVPMFH_FP::contraction42 (double **T4, double *T2, double *Sol)
+    
+    
+//Description: Double contraction of 4th order tensor with a 2nd order tensors
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-{
+void VEVPMFH_FP::contraction42 (double **T4, double *T2, double *Sol){
 
+  /** Declaration **/
   int i,j;
   
+  /**** Start ****/
   for (i=0; i<6; i++) {
     Sol[i] = 0;
 
@@ -140,15 +161,19 @@ void VEVPMFH_FP::contraction42 (double **T4, double *T2, double *Sol)
   
 }
 
-/* ************************************************************** 
-// ** CALCUL DE LA CONTRACTION DE DEUX TENSEURS D'ORDRE 4(6*6) **
-// **  RENVOIE UN TENSEUR D'ORDRE 4 SOUS FORME 6*6             **
-// **************************************************************/
 
-void VEVPMFH_FP::contraction44 (double **TA, double **TB, double **Tres)
-{
+
+//Description: Double contraction of two 4th order tensors
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::contraction44 (double **TA, double **TB, double **Tres){
+  
+  /** Declaration **/
   int i,j,k;
 
+  /**** Start ****/
   for (i=0; i<6; i++) {
 
     for (j=0; j<6; j++) {
@@ -162,53 +187,66 @@ void VEVPMFH_FP::contraction44 (double **TA, double **TB, double **Tres)
 
 }
 
-/* ********************************************************************* 
-// ** CALCUL DE LA DOUBLE CONTRACTION DE DEUX TENSEURS D'ORDRE 4(6*6) **
-// **  RENVOIE UN SCALAIRE                                            **
-// *********************************************************************/
 
-double VEVPMFH_FP::dcontr44 (double **TA, double **TB)
-{
+
+//Description: Double contraction of two 4th order tensors
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+double VEVPMFH_FP::dcontr44 (double **TA, double **TB){
+  
+  /** Declaration **/
   int i,j;
   double sol;
 
+  /**** Start ****/
   sol = 0;
 
   for (i=0; i<6; i++) {
     for (j=0; j<6; j++) {
-      sol+=TA[i][j]*TB[i][j];
+      sol = sol + TA[i][j]*TB[i][j];
     }
   }
+  
   return(sol);
+
 }
 
-/* ***************************************************************************************************************** 
-// ** CALCUL D'UNE CONTRACTION DE DEUX TENSEURS D'ORDRE 4(6*6) SUIVIE D'UNE CONTRACTION AVEC UN TENSEUR D'ORDRE 2 **
-// **  RENVOIE UN TENSEUR D'ORDRE 2                                                                               **
-// *****************************************************************************************************************/
 
-void VEVPMFH_FP::contraction442 (double **A, double **B, double *C, double *Tres)
-{
+
+//Description: Double contraction of two 4th order tensors followed by a double contraction with 2nd order tensor 
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::contraction442 (double **A, double **B, double *C, double *Tres){
+  
+  /** Declaration **/
   double temp[6];
 
+  /**** Start ****/
   VEVPMFH_FP::contraction42(B,C,temp);
   VEVPMFH_FP::contraction42(A, temp, Tres);
+  
 }
 
-/* ************************************************************************************************************
-// ** DOUBLE CONTRACTION OF 3 FOURTH ORDER TENSORS                                                           **
-// ** 2 FIRST INDICES OF THE FIRST TENSOR ARE REPEATED - 2 LAST INDICES OF THE LAST TENSOR ARE REPEATED ALSO **
-// ** THIS OPERATION CAN THUS BE VIEWED AS CONTRACTIONS I2-TENS4-TENS4-TENS4-I2                              **
-// ** A IS SUPPOSED FULLY SYMMETRIC SO THAT I2-TENS4 = TENS4-I2                                              **
-// ************************************************************************************************************/
 
-double VEVPMFH_FP::contr444_reduced (double **A, double **B, double **C)
-{
+
+//Description: Double contraction of three 4th order tensors
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+double VEVPMFH_FP::contr444_reduced (double **A, double **B, double **C){
+
+  /** Declaration **/
   int i;
   double I2[6];
   double *temp1, *temp2, *temp3;
   double sol;
 
+  /** Allocation & Init. **/
   for (i=0; i<3; i++) {
     I2[i] = 1;
     I2[i+3] = 0;
@@ -217,81 +255,115 @@ double VEVPMFH_FP::contr444_reduced (double **A, double **B, double **C)
   temp1 = (double*)malloc(6*sizeof(double));
   temp2 = (double*)malloc(6*sizeof(double));
   temp3 = (double*)malloc(6*sizeof(double));
-
+  
+  /**** Start ****/
   VEVPMFH_FP::contraction42(A, I2, temp1);
   VEVPMFH_FP::contraction42(C, I2, temp2);
   VEVPMFH_FP::contraction42(B, temp2, temp3);
   sol=VEVPMFH_FP::contraction22(temp1, temp3);
-
+  
+  /** Freeing **/
   free(temp1);
   free(temp2);
   free(temp3);
 
   return(sol);
+
 }
 
-/* ************************************************************************************************
-// ** DOUBLE CONTRACTION OF THREE FOURTH ORDER TENSORD; LAST CONTRACTION BEING OVER FOUR INDICES **
-// ************************************************************************************************/
 
-double VEVPMFH_FP::contraction444 (double **A, double **B, double **C)
-{
+
+//Description: Double contraction of three 4th order tensors
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+double VEVPMFH_FP::contraction444 (double **A, double **B, double **C){
+  
+  /** Declaration **/
   double **temp;
   double sol;
 
+  /** Allocation & Init. **/
   VEVPMFH_FP::mallocmatrix(&temp, 6, 6);
 
+  /**** Start ****/
   VEVPMFH_FP::contraction44(A, B, temp);
   sol=VEVPMFH_FP::dcontr44(temp, C);
 
+  /** Freeing **/
   VEVPMFH_FP::freematrix(temp, 6);
 
   return(sol);
+  
 }
 
 
-/* ************************************************************ 
-// ** CALCUL DU PRODUIT TENSORIEL DE DEUX TENSEURS D'ORDRE 2 **
-// **  RENVOIE UN TENSEUR D'ORDRE 4 SOUS FORME MATRICIELLE   **
-// ***********************************************************/
 
-void VEVPMFH_FP::produit_tensoriel (double *TA, double *TB, double **Tres)
+//Description: Tensor product of two 2nd order tensors  
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-{
+void VEVPMFH_FP::produit_tensoriel (double *TA, double *TB, double **Tres){
+  
+  /** Declaration **/
   int i,j;
 
+  /**** Start ****/
   for (i=0; i<6; i++) {
     for (j=0; j<6; j++) {
       Tres[i][j] = TA[i]*TB[j];
-    };
-  };
+    }
+  }
   
 }
 
 
 
 
-/* *************************************************************       
-// ** DECOMPOSITION LU D'UNE MATRICE CAREE D'ORDRE QUELCONQUE **
-// *************************************************************/
-/*
- * "ludcmp": computes the LU decomposition of a square matrix.
+
+//Description: Max of two numbers 
+//Author: Mohamed Haddad
+//Created: april 2020
+//Copyright: See COPYING file that comes with this distribution
+
+double VEVPMFH_FP::max(double a, double b){
+
+  if (a>b) return(a);
+  else return (b);
+
+}
+
+
+
+
+
+
+
+
+
+//Description: DECOMPOSITION LU D'UNE MATRICE CAREE D'ORDRE QUELCONQUE
+//Author: Olivier Pierard
+//Created: april 2003
+//Copyright: See COPYING file that comes with this distribution
+
+/* "ludcmp": computes the LU decomposition of a square matrix.
  * Original: Olivier Pierard (april 2003)
  * Source : Numerical recipies in C pp44-50(available from internet)
- * Some modifications about indices were rapported
- */
+ * Some modifications about indices were rapported */
 
 #define TINY 1.0e-14
 
-void VEVPMFH_FP::ludcmp(double **a, int n, int *indx, double *d, int *error)
-     /* Given a matrix a(n*n), this routine replaces it by the LU decomposition of a rows permutation 
+void VEVPMFH_FP::ludcmp(double **a, int n, int *indx, double *d, int *error){
+/* Given a matrix a(n*n), this routine replaces it by the LU decomposition of a rows permutation 
 itself.  a and n are input.  a is output containing L and U; indx[1..n] si an output vector that records 
 the row permutation effected y the partial pivoting, d is outpus as +/-1 depending on whether the number
 of row interchanges was even or odd, respectively.  This routine is used in combination with lubksb to 
 solve linear equations or invert a matrix.  error=1 if matrix is singular(but LU decomposition possible).
 If one row is null, the function is automatically stopped*/
 
-{
+
   int i,imax,j,k;
   double big, dum, sum, temp;
   double *vv;
@@ -362,11 +434,15 @@ If one row is null, the function is automatically stopped*/
   free(vv);
 }
 
-/**************************************
-** FORWARD AND BACKWARD SUBSTITUTION **
-***************************************/
 
-void VEVPMFH_FP::lubksb (double **a, int n, int *indx, double b[])
+
+
+//Description: FORWARD AND BACKWARD SUBSTITUTION
+//Author: Olivier Pierard
+//Created: april 2003
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::lubksb (double **a, int n, int *indx, double b[]){
 
      /* Solves the set of n linear equations A.X=b.  Here a(n*n) is input, not as the matrix A but rather 
 as its LU decomposition, determined by the routine ludcmp.  indx[1..n] is input as the permutation vector 
@@ -375,7 +451,7 @@ vector X.  a,n, and indx are not modified by this routine and can be left in pla
 with different right-hand sides b.  This routine takes into account the possibility that b will begin 
 with many zero elements, do it is efficient for use in matrix inversion. */
 
-{
+
   int i, ii=-1, ip, j;
   double sum;
 
@@ -403,13 +479,15 @@ with many zero elements, do it is efficient for use in matrix inversion. */
 
 }
 
-/************************************************
-** CALCUL DE L'VEVPMFH_FP::inverse D'UNE MATRICE D'ORDRE n **
-*************************************************/
 
-void VEVPMFH_FP::inverse (double **a, double **y, int *error, int n)
 
-{
+//Description: Compute the inverse of matrix
+//Author: Olivier Pierard
+//Created: april 2003
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::inverse (double **a, double **y, int *error, int n){
+
   double d;
   double *col;
   int i,j,*indx;
@@ -452,18 +530,17 @@ void VEVPMFH_FP::inverse (double **a, double **y, int *error, int n)
 }
 
 
-/* ***********************************************************       
-// ** Extraction of isotropic part of a fourth order tensor **
-// ***********************************************************/
-/*
- * "extractiso": Extract the isotropic part of a fourth order tensor.
- * Original: Olivier Pierard (April 2003)
- * INPUT : A fourth order tensor (6*6)
- * OUTPUT : An isotropic fourth order tensor
- */
 
-void VEVPMFH_FP::extractiso (double **tensani, double **tensiso)
 
+//Description: Extract the isotropic part of a 4th order tensor 
+//Author: Olivier Pierard
+//Created: april 2003
+//Modified:  02/04/2020 by M.Haddad
+//Copyright: See COPYING file that comes with this distribution
+//** INPUT : A fourth order tensor (6*6)
+//** OUTPUT : An isotropic fourth order tensor
+
+void VEVPMFH_FP::extractiso (double **tensani, double **tensiso)
 {
   double deltaij[6], deltakl[6];
   double I[6][6]={0};
@@ -519,24 +596,20 @@ void VEVPMFH_FP::extractiso (double **tensani, double **tensiso)
 
 
 
-/* ***********************************************       
-// ** Memory allocation for a table of pointers **
-// ***********************************************/
-/*
- * "VEVPMFH_FP::mallocmatrix": Memory allocation for a table of pointers.
- * Original: April 2003
- * INPUT : m : the number of lines(number of boxes for pointers)
- *         n : the number of rows
- *         **A : A declared pointer to pointers not yet initialized
- * OUTPUT : 
- */
-
-void VEVPMFH_FP::mallocmatrix (double ***A, int m, int n)
-
-{
+//Description: Memory allocation for a table of pointers 
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+ 
+void VEVPMFH_FP::mallocmatrix (double ***A, int m, int n){
+  
+  /** Dec **/
   int i,j;
 
+  /** Alloc. & Init. **/
   *A = (double**)malloc(m*sizeof(double*));
+  
+  /**** Start ****/
   for(i=0;i<m;i++){
     (*A)[i] = (double*)malloc(n*sizeof(double));
   }
@@ -546,53 +619,49 @@ void VEVPMFH_FP::mallocmatrix (double ***A, int m, int n)
       (*A)[i][j]=0;
     }
   }
+  
 }
 
 
 
 
-/* *******************************************       
-// ** Freeing memory of a table of pointers **
-// *******************************************/
-/*
- * "VEVPMFH_FP::freematrix": Freeing memory of a table of pointers.
- * Original: April 2003
- * INPUT : Tensor which has to be freed, number of lines
- * OUTPUT : nothing; memory cleaned
- */
-
-void VEVPMFH_FP::freematrix (double **A, int m)
+//Description: Freeing memory of a table of pointers
+//Author: Mohamed Haddad
+//Created: 01/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-{
+void VEVPMFH_FP::freematrix (double **A, int m){
+  
+  /** Dec **/
   int i;
 
+  /**** Start ****/
   for (i=0;i<m;i++) {
     free (A[i]);
   }
+  
+  /** Freeing **/
   free (A);
 
 }
 
 
-/* ***************************************       
-// ** Memory allocation for a 3D matrix **
-// ***************************************/
-/*
- * "VEVPMFH_FP::mallocmatrix3D": Memory allocation for a 3D matrix.
- * Original: July 2003
- * INPUT : p : first and main indice (e.g.: number of colloc. points)
- *         m : the number of lines(number of boxes for pointers)
- *         n : the number of rows
- *         ***A : A declared pointer to pointers to pointers not yet initialized
- * OUTPUT : A initialized
- */
 
-void VEVPMFH_FP::mallocmatrix3D (double ****A, int p, int m, int n)
 
-{
+//Description: Memory allocation for a 3D matrix
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::mallocmatrix3D (double ****A, int p, int m, int n){
+  
+  /** Dec **/
   int i,j,k;
 
+  /** Alloc. & Init. **/
   *A = (double***)malloc(p*sizeof(double**));
+  
+  /**** Start ****/
   for(k=0;k<p;k++){
     (*A)[k] = (double**)malloc(m*sizeof(double*));
   }
@@ -605,50 +674,53 @@ void VEVPMFH_FP::mallocmatrix3D (double ****A, int p, int m, int n)
       }
     }
   }
+
 }
 
 
-/* ***********************************       
-// ** Freeing memory of a 3D matrix **
-// ***********************************/
-/*
- * "VEVPMFH_FP::freematrix3D": Freeing memory of a 3D matrix.
- * Original: July 2003
- * INPUT : 3D matrix which has to be freed, p and m : main and second indices
- * OUTPUT : nothing; memory cleaned
- */
 
-void VEVPMFH_FP::freematrix3D (double ***A, int p, int m)
 
-{
+
+//Description: Freeing memory of a 3D matrix
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::freematrix3D (double ***A, int p, int m){
+
+  /** Dec **/
   int k, i;
 
+  /**** Start ****/
   for (k=0; k<p; k++) {
     for (i=0; i<m; i++) {
       free (A[k][i]);
     }
     free (A[k]);
   }
+  
+  /** Freeing **/
   free(A);
+  
 }
 
-/*******************************************************************************
-// ** COMPUTE THE ISOTROPIC HOOKE'S ELASTIC OPERATOR FROM E AND NU **
-********************************************************************************/
 
-/* June 2003
-// INPUT : E : Elastic young's modulus in the longitudinal direction (same in other directions if isotropic)
-//         nu : Poisson's ratio in the longitudinal direction (same in other directions if isotropic)
-// OUTPUT : Isotropic stiffness tensor stored as in book of ID p.555
-*/
 
-void VEVPMFH_FP::compute_Hooke(double E, double nu, double **C)
 
-{
+
+//Description: COMPUTE THE ISOTROPIC HOOKE'S ELASTIC OPERATOR FROM E AND NU 
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::compute_Hooke(double E, double nu, double **C){
+  
+  /** Dec **/
   double lambda;
   double mu;
   int i,j;
 
+  /**** Start ****/
   if(nu==0.5) {
     printf("nu = 0.5, Impossible to compute lambda\n");
     exit(0);
@@ -673,46 +745,67 @@ void VEVPMFH_FP::compute_Hooke(double E, double nu, double **C)
     C[i][i] = C[i][i] + 2*mu;
   }
      
+}
 
 
-}
 
-/********************************************/
-/** Print on the screen values of a vector **/
-/********************************************/
 
-void VEVPMFH_FP::printvect(double *vect, int n)
-{
+
+
+//Description:  Print on the screen values of a vector
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::printvect(double *vect, int n){
+  
+  /** Dec **/
   int i;
 
+  /**** Start */
   for (i=0;i<n;i++) {
     printf("\t %g\t",vect[i]);
   }
   printf("\n");
+  
 }
 
-/**************************************************/
-/** Print values of a vector to a specified file **/
-/**************************************************/
 
-void VEVPMFH_FP::fprintvect(FILE *pfile, double *vect, int n)
-{
-  int i;
 
+
+//Description: Print values of a vector to a specified file 
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::fprintvect(FILE *pfile, double *vect, int n){
+
+  /** Dec **/
+  int i;
+ 
+  /** Start **/ 
   for (i=0; i<n; i++) {
     fprintf(pfile, "\t %g\t",vect[i]);
   }
   fprintf(pfile, "\n");
+  
 }
 
-/************************************************/
-/** Print on the screen values of a m*n matrix **/
-/************************************************/
 
-void VEVPMFH_FP::printmatr(double **matr4, int m, int n)
-{
+
+
+
+//Description: Print on the screen values of a m*n matrix  
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::printmatr(double **matr4, int m, int n){
+  
+  /** Dec **/
   int i, j;
 
+  /**** Start **/
   for (i=0; i<m; i++) {
     for (j=0; j<n; j++) {
       printf ("\t %g\t ", matr4[i][j]);
@@ -720,16 +813,23 @@ void VEVPMFH_FP::printmatr(double **matr4, int m, int n)
     printf("\n");
   }
   printf("\n");
+  
 }   
 
-/******************************************************/
-/** Print values of a m*n matrix to a specified file **/
-/******************************************************/
 
-void VEVPMFH_FP::fprintmatr(FILE *pfile, double **matr4, int m, int n)
-{
+
+
+//Description: Print values of a m*n matrix to a specified file   
+//Author: Mohamed Haddad
+//Created: 03/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::fprintmatr(FILE *pfile, double **matr4, int m, int n){
+
+  /** Dec **/
   int i, j;
 
+  /**** Start **/
   for (i=0; i<m; i++) {
     for (j=0; j<n; j++) {
       fprintf (pfile, "\t %g\t ", matr4[i][j]);
@@ -737,38 +837,47 @@ void VEVPMFH_FP::fprintmatr(FILE *pfile, double **matr4, int m, int n)
     fprintf(pfile, "\n");
   }
   fprintf(pfile, "\n");
+  
 }   
 
-/************************************************/
-/** Compute the trace of a second order tensor **/
-/************************************************/
 
-/* October 2003 */
-/* New prototype in May 2004 */
 
-double VEVPMFH_FP::trace(double *tens2)
-{
+
+//Description: Compute the trace of a second order tensor   
+//Author: Mohamed Haddad
+//Created: 05/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+double VEVPMFH_FP::trace(double *tens2){
+
+  /** Dec **/
   int i;
   double tr=0;
 
   for (i=0;i<3;i++) {
     tr = tr + tens2[i];
   }
+  
   return (tr);
+
 }
 
-/**********************************************************/
-/** Extract the deviatoric part of a second order tensor **/
-/**********************************************************/
 
 
-/* October 2003 */
 
-void VEVPMFH_FP::dev(double *tens2, double *dev_tens2)
-{
+
+//Description: Extract the deviatoric part of a second order tensor  
+//Author: Mohamed Haddad
+//Created: 05/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::dev(double *tens2, double *dev_tens2){
+
+  /** Dec **/
   int i;
   double tr_tens2=0;
 
+  /**** Start ****/
   for (i=0;i<3;i++) {
     tr_tens2 = tr_tens2 + tens2[i];
   }
@@ -784,88 +893,50 @@ void VEVPMFH_FP::dev(double *tens2, double *dev_tens2)
 }
 
 
-/*********************************************/
-/** Compute the von Mises equivalent stress **/
-/*********************************************/
-
-/* October 2003 */
-/* New prototype in May 2004 */
-
-double VEVPMFH_FP::vonmises(double *tens2)
-{
-  return(pow(0.5*(pow(tens2[0]-tens2[1],2) + pow(tens2[1]-tens2[2],2) + pow(tens2[2]-tens2[0],2)) + 3*pow(tens2[3]/sqrt(2),2) + 3*pow(tens2[4]/sqrt(2),2) + 3*pow(tens2[5]/sqrt(2),2),0.5));
-
-}
-
-/******************************************/
-/** Compute the equivalent strain scalar **/
-/******************************************/
-
-/* November 2005 */
-
-double VEVPMFH_FP::eqstrn(double *tens2)
-{
-  return(pow(2./9.*(pow(tens2[0]-tens2[1],2) + pow(tens2[1]-tens2[2],2) + pow(tens2[2]-tens2[0],2)) + 4./3.*pow(tens2[3]/sqrt(2),2) + 4./3.*pow(tens2[4]/sqrt(2),2) + 4./3.*pow(tens2[5]/sqrt(2),2),0.5));
 
-}
 
-/***********************/
-/* Compare two vectors */
-/***********************/
 
-/* Return 1 if the two vectors are identical, 0 at least one component of the two vectors is different */
 
-/* November 2003 */
-/* Condition of equality modified in April 2004 */
+//Description: Compute the von Mises equivalent stress  
+//Author: Mohamed Haddad
+//Created: 05/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-int VEVPMFH_FP::compare(double *A, double *B, int size)
-{
-  int i;
+double VEVPMFH_FP::vonmises(double *tens2){
 
-  for (i=0; i<size; i++) {
-    if (fabs(A[i]-B[i])>(1e-10*fabs(A[i]))) {
-      return(0);
-    }
-  }
-  return(1); 
+  return(pow(0.5*(pow(tens2[0]-tens2[1],2) + pow(tens2[1]-tens2[2],2) + pow(tens2[2]-tens2[0],2)) + 3*pow(tens2[3]/sqrt(2),2) + 3*pow(tens2[4]/sqrt(2),2) + 3*pow(tens2[5]/sqrt(2),2),0.5));
 
 }
 
 
-/************************/
-/* Compare two matrices */
-/************************/
 
-/* Return 1 if the two matrices are identical, 0 otherwise */
 
-/* April 2003 */
 
-int VEVPMFH_FP::comparematr(double **A, double **B, int m, int n)
-{
-  int i;
+//Description: Compute the equivalent strain   
+//Author: Mohamed Haddad
+//Created: 05/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-  for (i=0; i<m; i++) {
-    if (!compare(A[i], B[i], n))
-      return(0);
-  }
-  return(1); 
+double VEVPMFH_FP::eqstrn(double *tens2){
+  return(pow(2./9.*(pow(tens2[0]-tens2[1],2) + pow(tens2[1]-tens2[2],2) + pow(tens2[2]-tens2[0],2)) + 4./3.*pow(tens2[3]/sqrt(2),2) + 4./3.*pow(tens2[4]/sqrt(2),2) + 4./3.*pow(tens2[5]/sqrt(2),2),0.5));
 
 }
 
 
 
-/*****************************/
-/* Check if a vector is null */
-/*****************************/
 
-/* Return 1 il all the elements are nil, 0 otherwise */
 
-/* November 2003 */
+//Description: Check if a vector is null 
+//Author: Mohamed Haddad
+//Created: 05/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-int VEVPMFH_FP::isnull(double *A, int size)
-{
+int VEVPMFH_FP::isnull(double *A, int size){
+
+  /** Dec & Init. **/
   int i, nil=1;
 
+  /**** Start ****/
   for (i=0; i<size; i++) {
     if (A[i]!=0) {
       nil=0;
@@ -877,21 +948,28 @@ int VEVPMFH_FP::isnull(double *A, int size)
 
 }
 
-/***********************************************/
-/* Check if a fourth order tensor is isotropic */
-/***********************************************/
 
-/* April 2004 - O. Pierard */
-/* Input: A: a fourth order tensor
-   Output: res: 1 if A is isotropic; 0 otherwise */
 
-int VEVPMFH_FP::isisotropic(double**A)
-{
+
+
+
+//Description: Check if a fourth order tensor is isotropic
+//Author: O. Pierard
+//Created: April 2004
+//Copyright: See COPYING file that comes with this distribution
+//
+//**Input: A: a fourth order tensor
+//**Output: res: 1 if A is isotropic; 0 otherwise 
+
+int VEVPMFH_FP::isisotropic(double**A){
+
+  /** Dec & Init. **/ 
   double lambda, mu;
   int i, j;
   int res = 1;
   double tol=1e-2;
 
+  /**** Start ****/
   lambda = A[0][1];
   mu = A[3][3]/2;
 
@@ -918,138 +996,43 @@ int VEVPMFH_FP::isisotropic(double**A)
   }
 
   return res;
+  
 }
 
-/******************************/
-/* Copy a vector into another */
-/******************************/
-
-/* November 2003 */
-
-void VEVPMFH_FP::copyvect(double *A, double *B, int size)
-{
-  int i;
-
-  for (i=0; i<size; i++) {
-    B[i] = A[i];
-  }
-}
-
-/***********************************************/
-/* Find the maximum positive value of a vector */
-/***********************************************/
-
-/* April 2004 */
-
-double VEVPMFH_FP::maxvect(double *vect, int n)
-/* vect : vector n*1
-   n : length of the vector
-   return a double, the maximum value of the vector */
-{
-
-  int i;
-  double max = -1e20;
-
-  for (i=0; i<n; i++) {
-    if (vect[i]>max)
-      max = vect[i];
-  }
-  return max;
-
-}
-
-
-/********************************************************/
-/* Find the maximum value in absolute value of a vector */
-/********************************************************/
-
-/* April 2004 */
-
-double VEVPMFH_FP::maxvectabs(double *vect, int n)
-/* vect : vector n*1
-   n : length of the vector
-   return a double, the maximum value in absolute value of the vector */
-{
-
-  int i;
-  double max = 0;
 
-  for (i=0; i<n; i++) {
-    if (fabs(vect[i])>max)
-      max = fabs(vect[i]);
-  }
-  return max;
 
-}
 
-/****************************************************************/
-/* Find the minimum non nil value in absolute value of a vector */
-/****************************************************************/
+//Description: Copy a vector into another 
+//Author: Mohamed Haddad 
+//Created: 06/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
-/* May 2004 */
+void VEVPMFH_FP::copyvect(double *A, double *B, int size){
 
-/* matr : matrix m*n
-   m : size of the vector */
-
-double VEVPMFH_FP::minvectabs(double *vect, int m)
-
-{
+  /** Dec **/
   int i;
-  double min = 1e100;
-  for (i=0; i<6; i++) {
-    if ((fabs(vect[i])<min)&&(vect[i]!=0))
-      min = fabs(vect[i]);
-  }
-
-  return min;
-
-}
-
-/****************************************************************/
-/* Find the minimum non nil value in absolute value of a matrix */
-/****************************************************************/
-
-/* May 2004 */
 
-/* matr : matrix m*n
-   m, n : size of the matrix */
-
-double VEVPMFH_FP::minmatrabs(double **matr, int m, int n)
-
-{
-  int i, j;
-  double min = 1e100;
-  for (i=0; i<6; i++) {
-    for (j=0; j<6; j++) {
-      if ((fabs(matr[i][j])<min)&&(matr[i][j]!=0)) 
-	min = fabs(matr[i][j]);
-    }
+  for (i=0; i<size; i++) {
+    B[i] = A[i];
   }
-
-  return min;
-
+  
 }
 
-double VEVPMFH_FP::max(double a, double b){
-
-if (a>b) return(a);
-else return (b);
-
-}
 
 
-/*****************************************************************/
-/* 		Compute Identity tensor I			 */
-/*****************************************************************/
 
-/* November 2020 */
 
+//Description: Compute Identity tensor I
+//Author: Mohamed Haddad 
+//Created: November 2020 
+//Copyright: See COPYING file that comes with this distribution
 
 void VEVPMFH_FP::compute_I(double **I){
-
+  
+  /** Dec **/
   int i,j;
 
-
+  /**** Start ****/
   for(j=0;j<6;j++){
     for(i=0;i<6;i++)
       I[i][j] = 0;  
@@ -1060,17 +1043,26 @@ void VEVPMFH_FP::compute_I(double **I){
 
 }
 
-/** FUNCTION TO COMPUTE Ivol (6*6) TENSOR **/
-/*******************************************/
-void VEVPMFH_FP::compute_Ivol(double **I_vol){
 
+
+
+//Description: FUNCTION TO COMPUTE Ivol (6*6) TENSOR
+//Author: Mohamed Haddad 
+//Created: November 2020 
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::compute_Ivol(double **I_vol){
+  
+  /** Dec **/
   int i,j;
   double *deltaij;
   double *deltakl;
-
+  
+  /** Alloc. **/
   deltaij=(double*)malloc(6*sizeof(double));
   deltakl=(double*)malloc(6*sizeof(double));
 
+  /**** Start ****/
   for (i=0;i<3;i++) {
     deltaij[i] = 1;
     deltakl[i] = 1;
@@ -1087,64 +1079,129 @@ void VEVPMFH_FP::compute_Ivol(double **I_vol){
       }
     }
 
-free(deltaij);
-free(deltakl);
+  /** Freeing **/
+  free(deltaij);
+  free(deltakl);
 
 }
 
-/** FUNCTION TO COMPUTE Idev (6*6) TENSOR **/
-/*******************************************/
-void VEVPMFH_FP::compute_Idev(double **I_dev){
 
-double **I;
-double **Ivol;
 
 
-VEVPMFH_FP::mallocmatrix(&I,6,6);
-VEVPMFH_FP::mallocmatrix(&Ivol,6,6);
 
-VEVPMFH_FP::compute_I(I);
-VEVPMFH_FP::compute_Ivol(Ivol);
 
-VEVPMFH_FP::addtens4(1,I,-1,Ivol,I_dev);
+//Description: FUNCTION TO COMPUTE Idev (6*6) TENSOR
+//Author: Mohamed Haddad 
+//Created: 06/04/2020
+//Copyright: See COPYING file that comes with this distribution
 
+void VEVPMFH_FP::compute_Idev(double **I_dev){
+  
+  /** Declaration **/
+  double **I, **Ivol;
 
+  /** Allocation **/
+  VEVPMFH_FP::mallocmatrix(&I,6,6);
+  VEVPMFH_FP::mallocmatrix(&Ivol,6,6);
 
-VEVPMFH_FP::freematrix(I,6);
-VEVPMFH_FP::freematrix(Ivol,6);
+  /**** Start ****/
+  VEVPMFH_FP::compute_I(I);
+  VEVPMFH_FP::compute_Ivol(Ivol);
 
+  VEVPMFH_FP::addtens4(1,I,-1,Ivol,I_dev);
 
 
+  /** Freeing **/
+  VEVPMFH_FP::freematrix(I,6);
+  VEVPMFH_FP::freematrix(Ivol,6);
 
 }
 
 
-/** FUNCTION TO COMPUTE ZERO VECTOR (SIZE * 1) **/
-/************************************************/
 
-void VEVPMFH_FP::zero_vect(double size, double *vect){
-   
-   int i;
-   
-   for(i=0;i<size;i++)
-     vect[i] = 0;
 
+
+//Description: Read material properties & simulation parameters from .dat file
+//Author: Mohamed Haddad 
+//Created: 06/04/2020
+//Copyright: See COPYING file that comes with this distribution
+void VEVPMFH_FP::mallocMATPROP(VEVPMFH_FP::MATPROP* matr, VEVPMFH_FP::MATPROP* incl, const char* inputfilename){
+  FILE *inputfile;
+  char c;
+  int i;
+  
+  int matr_nbrGi, matr_nbrKj, incl_nbrGi, incl_nbrKj;
+  
+  /**** Start ****/
+  // Open input .dat file
+  inputfile = fopen(inputfilename,"r");
+
+  if (inputfile==NULL) {
+    printf("\n The file %s cannot be opened, check that the file currently exists. \n", inputfilename);
+    exit(1);
+  }
+  
+  // Read kine by line 
+  // Matrix properties
+  
+  for(i=0;i<18;i++) {
+    while((c=fgetc(inputfile))!='\n');
+  }
+  fscanf(inputfile,"%i",&matr_nbrGi);
+  
+  for(i=0;i<(2*matr_nbrGi+1);i++) {
+    while((c=fgetc(inputfile))!='\n');
+  }
+  fscanf(inputfile,"%i",&matr_nbrKj);
+  
+  for(i=0;i<26;i++) {
+    while((c=fgetc(inputfile))!='\n');
+  }
+  fscanf(inputfile,"%i",&incl_nbrGi);
+  
+  for(i=0;i<(2*incl_nbrGi+1);i++) {
+    while((c=fgetc(inputfile))!='\n');
+  }
+  fscanf(inputfile,"%i",&incl_nbrKj);
+  
+  matr->gi = (double*)calloc(matr_nbrGi,sizeof(double));
+  matr->gi = (double*)calloc(matr_nbrGi,sizeof(double));
+  
+  matr->Kj = (double*)calloc(matr_nbrKj,sizeof(double));
+  matr->kaj = (double*)calloc(matr_nbrKj,sizeof(double));
+  
+  
+  incl->gi = (double*)calloc(incl_nbrGi,sizeof(double));
+  incl->gi = (double*)calloc(incl_nbrGi,sizeof(double));
+  
+  incl->Kj = (double*)calloc(incl_nbrKj,sizeof(double));
+  incl->kaj = (double*)calloc(incl_nbrKj,sizeof(double));
+  
+  
+  fclose(inputfile);
 }
 
 
 
 
-/** !! Copyin fonctions_secant here **/
-/**************************************/
 
-void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl, const char* inputfilename)
-{
+
+//Description: Read material properties & simulation parameters from .dat file
+//Author: Mohamed Haddad 
+//Created: 06/04/2020
+//Copyright: See COPYING file that comes with this distribution
+
+void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl, const char* inputfilename){
+
+  /** Dec **/
   FILE *inputfile;
   char c;
   int i;
   int temp=0;
   double temp1=0;
   
+  /**** Start ****/
+  // Open input .dat file
   inputfile = fopen(inputfilename,"r");
 
   if (inputfile==NULL) {
@@ -1152,11 +1209,11 @@ void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl
     exit(1);
   }
 
+  // Read input file line by line 
+  // Matrix properties
   for(i=0;i<6;i++) {
     while((c=fgetc(inputfile))!='\n');
   }
-  /**	Matrix proprieties	**/
-  /*******************************/
   fscanf(inputfile,"%lf",&(matr->E0));
   while((c=fgetc(inputfile))!='\n');
   fscanf(inputfile,"%lf",&(matr->nu));
@@ -1244,8 +1301,8 @@ void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl
     exit(0);
   }
 
-  /**	Inclusions proprieties		**/
-  /***************************************/
+ 
+  // Inclusion properties 
   for(i=0;i<4;i++) {
     while((c=fgetc(inputfile))!='\n');
   }
@@ -1281,7 +1338,7 @@ void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl
   fscanf(inputfile,"%i",&(incl->nbrGi));
 
 
- incl->Gi = (double*) malloc((incl->nbrGi)*sizeof(double));
+  incl->Gi = (double*) malloc((incl->nbrGi)*sizeof(double));
 
   for(i=0;i<incl->nbrGi;i++) {
     while((c=fgetc(inputfile))!='\n');
@@ -1330,16 +1387,16 @@ void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl
     temp1+=incl->Kj[i];
 
   incl->K_inf = incl->K0 - temp1;
-
   
-  /**	Matrix volume fraction	**/
+  // Matrix volume fraction
   matr->v = 1 - incl->v;
-  matr->ar=0;              /* Aspect Ratio has no sense for the matrix */
+  matr->ar=0;         /* Aspect Ratio has no sense for the matrix */
   matr->orient[0]=0;  /* Could represent the orientation of the RVE in a global referential */
   matr->orient[1]=0;
   matr->orient[2]=0;
 
 
+  // Avoid case of zero relaxation times (if yes exit(0)) 
   for(i=0;i<incl->nbrGi;i++)
   if (incl->gi[i] == 0){
     printf("FATAL ERROR: gi must be not nil for the inclusions\n");
@@ -1353,9 +1410,8 @@ void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl
   }
 
  
-  /**	Inclusions orientation		**/
-  /***************************************/
-
+  
+  // Inclusion's orientation 
   for(i=0;i<4;i++) {
     while((c=fgetc(inputfile))!='\n');
   }
@@ -1366,84 +1422,29 @@ void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl
   while((c=fgetc(inputfile))!='\n');
   fscanf(inputfile,"%lf",&(incl->orient[2]));
   
-  /**		End read inputs			**/
-  /*****************************************************/
+  // End reading & close input file 
   fclose(inputfile);
+  
 }
 
  
  
-/* *************************************************
-// ** COMPUTATION OF STRESS IN LINEAR ELASTICTIY  **
-// *************************************************/
-
-/* May 2004 - O. Pierard
-// INPUT: Young modulus and Poisson's ratio 
-//        Strain tensor
-// OUTPUT: Stress tensor */
-
-void VEVPMFH_FP::computestrslinel(double young, double nu, double *strn, double *strs) {
-
-  double **C;
-
-  VEVPMFH_FP::mallocmatrix(&C, 6, 6);
-
-  VEVPMFH_FP::compute_Hooke(young, nu, C);
-  VEVPMFH_FP::contraction42(C, strn, strs);
-  
-  VEVPMFH_FP::freematrix(C, 6);
-}
-
-
-
-/* *************************************************
-// ** COMPUTATION OF STRAIN IN LINEAR ELASTICTIY  **
-// *************************************************/
-
-/* May 2004 - O. Pierard
-// INPUT: Young modulus and Poisson's ratio 
-//        Stress tensor
-// OUTPUT: Strain tensor */
-
-void computestrnlinel(double young, double nu, double *strs, double *strn) {
-
-  double **C, **S;
-  int error = 0;
-
-  VEVPMFH_FP::mallocmatrix(&C, 6, 6);
-  VEVPMFH_FP::mallocmatrix(&S, 6, 6);
-
-  VEVPMFH_FP::compute_Hooke(young, nu, C);
-  VEVPMFH_FP::inverse(C, S, &error, 6);
-  if (error != 0) {
-    printf("FATAL ERROR: Error occured in the inversion of C in computestrnlinel\n");
-    exit(0);
-  }
-
-  VEVPMFH_FP::contraction42(S, strs, strn);
-  
-  VEVPMFH_FP::freematrix(C, 6);
-  VEVPMFH_FP::freematrix(S, 6);
-}
-
-
+ 
 
+//Description: COMPUTE HARDENING
+//Author: Mohamed Haddad 
+//Created: 06/04/2020
+//Modified: December 2020 by M.Haddad (Addition of Mixed power and exponential law) 
+//Copyright: See COPYING file that comes with this distribution 
+// 										
+/**INPUT:	-ACUMULTED PLASTICITY
+//	  	-MATERIAL PARAMETERS
+//**OUTPUT: 	-R(P)	
+//	   	-dR/dp*/							
 
-/**************************************************************************************
-** 		COMPUTE HARDENING						   	**
-**************************************************************************************/
-/*2003 - O. Pierard 									*/	
-/*Modified in november 2005: order of arguments 					*/
-/*MODIFIED IN DECEMBER 2020 BY M.HADDAD: Addition of Mixed power and exponential law */
-/* 											*/
-/*   INPUT : 	ACUMULTED PLASTICITY							*/
-/*		MATERIAL PARAMETERS							*/
-/*   OUTPUT :	R(P)									*/
-/*		dR/dp									*/
-/*************************************************************************************/
+void VEVPMFH_FP::compute_R (double p, const VEVPMFH_FP::MATPROP *mat, double *Rp, double *dRdp){
 
-void VEVPMFH_FP::compute_R (double p, VEVPMFH_FP::MATPROP *mat, double *Rp, double *dRdp)
-{
+  // Invalid hardening properties 
   if ((p<0.)||((mat->hardexp)<=0.)) {
     printf("Invalid parameter(s) for the hardening law in function VEVPMFH_FP::compute_R:\np: %g, hardmod1: %g, hardmod2: %g, hardexp: %g\n", p, mat->hardmodulus1 , mat->hardmodulus2, mat->hardexp);
     exit(0);
@@ -1452,7 +1453,7 @@ void VEVPMFH_FP::compute_R (double p, VEVPMFH_FP::MATPROP *mat, double *Rp, doub
 
   switch (mat->hardmodel) {
   
-  /* Simple power Law */
+  // Law 1: Power law
   case 1:
          
     /*if (hardexp<=1) 
@@ -1461,34 +1462,35 @@ void VEVPMFH_FP::compute_R (double p, VEVPMFH_FP::MATPROP *mat, double *Rp, doub
     if ((mat->hardexp)>1) 
       printf("WARNING: Exponent of the power law hardening function should be inferior to 1 to be physical (Lemaitre/Chaboche p168/177)!\n");
 
-    /** When p=0 **/
-    /**************/
+    // p=0
     if (p==0) {
       *Rp=0;
-      if ((mat->hardexp)>1.0+1e-10) /* n>1 */
+      if ((mat->hardexp)>1.0+1e-10) 					/* n>1 */
 	*dRdp = 0;
-      else if (fabs((mat->hardexp)-1.0)<1e-10) /* n=1 */
+      else if (fabs((mat->hardexp)-1.0)<1e-10) 			/* n=1 */
 	*dRdp = (mat->hardmodulus1)*(mat->hardexp);
-      else if (((mat->hardexp)>0)&&((mat->hardexp)<1.0-1e-10)) /* n < 1*/
+      else if (((mat->hardexp)>0)&&((mat->hardexp)<1.0-1e-10)) 	/* n < 1*/
 	*dRdp = (mat->hardmodulus1)*(mat->hardexp)/(pow(1e-10, 1-(mat->hardexp)));
       else
 	printf("Bad choice of hardening exponent with the Power law hardening model\n");
     }
-    /** When p is not nil **/
-    /************************/
-    else {      /* p>0 */
+    
+    // p!= 0
+    else {      
       *Rp = (mat->hardmodulus1)*pow(p,(mat->hardexp));
       *dRdp = (mat->hardexp)*(*Rp)/p ; 
     }
     break;
 
 
-  /* Exponential law */
-  case 2:       
+  // Exponential law 
+  case 2:      
+    // p=0 
     if (p==0) {
       *Rp=0;
       *dRdp = (mat->hardexp)*(mat->hardmodulus1);
     }
+    // p!= 0
     else {
       *Rp = (mat->hardmodulus1)*(1-exp(-(mat->hardexp)*p));
       *dRdp = (mat->hardexp)*(mat->hardmodulus1)*exp(-(mat->hardexp)*p);
@@ -1497,12 +1499,14 @@ void VEVPMFH_FP::compute_R (double p, VEVPMFH_FP::MATPROP *mat, double *Rp, doub
     
     
     
-  /** Mixed power and exponential law **/  
+  // Mixed power and exponential law 
   case 3:
+    // p=0
     if (p==0.) {
       *Rp = 0.;
       *dRdp = (mat->hardmodulus1) + (mat->hardmodulus2)*(mat->hardexp);
     }
+    // p!= 0
     else {
       *Rp = (mat->hardmodulus1)*p + (mat->hardmodulus2)*(1-exp(-(mat->hardexp)*p));
       *dRdp = (mat->hardmodulus1) +  (mat->hardexp)*(mat->hardmodulus2)*exp(-(mat->hardexp)*p) ;
@@ -1514,7 +1518,9 @@ void VEVPMFH_FP::compute_R (double p, VEVPMFH_FP::MATPROP *mat, double *Rp, doub
   default :
     printf("parameter chhard invalid in function VEVPMFH_FP::compute_R\n");
     exit(0);
-  }             /* End of switch */
+  
+  // End switch
+  }             
 
 
 }
@@ -1522,11 +1528,13 @@ void VEVPMFH_FP::compute_R (double p, VEVPMFH_FP::MATPROP *mat, double *Rp, doub
 
 
 
-/* ************************************************************
-// ** EVALUATION OF THE VISCOPLASTIC FUNCTION (GENERAL FORM) ** 
-// ************************************************************/
 
-/*    April 2004 - O. Pierard
+
+//Description: EVALUATION OF THE VISCOPLASTIC FUNCTION (GENERAL FORM)
+//Author: Mohamed Haddad 
+//Created: 06/04/2020
+//Copyright: See COPYING file that comes with this distribution 
+//
 //    INPUT : strseq : Von Mises equivalent stress
 //            p : Accumulated plasticity
 //            mat : material properties and choice of the viscoplastic function
@@ -1534,22 +1542,25 @@ void VEVPMFH_FP::compute_R (double p, VEVPMFH_FP::MATPROP *mat, double *Rp, doub
 //             dgvdstrseq : Derivative of gv w.r.t. strseq
 //             dgvdp : Derivative of gv w.r.t. p */
 
-void VEVPMFH_FP::computevisc (double strseq, double p, VEVPMFH_FP::MATPROP *mat, double *gv, double *dgvdstrseq, double *dgvdp)
-{
+void VEVPMFH_FP::computevisc (double strseq, double p, const VEVPMFH_FP::MATPROP *mat, double *gv, double *dgvdstrseq, double *dgvdp){
+  
+  /** Dec **/
   double Rp, dRdp, f;
 
+  /**** Start ****/
   VEVPMFH_FP::compute_R(p, mat, &Rp, &dRdp);
   f = strseq - (mat->yldstrs) - Rp;
 
+
+  
   switch (mat->viscmodel) {
-    /**** Norton's power law ****/
+  // Norton's power law 
   case 1:     
     if ((mat->visccoef<0)||(mat->viscexp1<1)) {
       printf("Invalid value of the creep coefficient (must be > 0) or creep exponent1 (must be > 1 for continuous derivatives - and also to be physical (Lemaitre/Chaboche p257)) in the Norton's viscous function\n");
       /*exit(0);*/
     }
 
-
     if (f <= 0) {
       *gv = 0;
       *dgvdstrseq = 0;
@@ -1564,7 +1575,7 @@ void VEVPMFH_FP::computevisc (double strseq, double p, VEVPMFH_FP::MATPROP *mat,
     break;
 
 
-    /**** 3 parameters creep function (power_law3 in Digimat and defined in ABAQUS) ****/
+  // 3 parameters creep function (power_law3 in Digimat and defined in ABAQUS)
   case 2:      
     if ((mat->visccoef<=0)||(mat->viscexp1<=0)) {
       printf("Invalid value of the creep coefficient (must be > 0) or creep exponent1 (must be > 0) in the 3 parameters creep function\n");
@@ -1577,14 +1588,16 @@ void VEVPMFH_FP::computevisc (double strseq, double p, VEVPMFH_FP::MATPROP *mat,
     }
     else {
       if (p==0)
-	p=1e-12;   /* Why p==0 && f>0 ???? */
+       p=1e-12;   /* Why p==0 && f>0 ???? */
       *gv = pow((mat->visccoef)*pow(strseq, mat->viscexp1), 1/(mat->viscexp2+1)) * pow((mat->viscexp2+1)*p, (mat->viscexp2)/(mat->viscexp2+1));
       *dgvdstrseq = (mat->viscexp1 * mat->visccoef * pow(strseq, mat->viscexp1 - 1) / (mat->viscexp2+1)) * pow(mat->visccoef*pow(strseq, mat->viscexp1), -mat->viscexp2/(mat->viscexp2 + 1)) * pow((mat->viscexp2+1)*p, (mat->viscexp2)/(mat->viscexp2+1));
       *dgvdp = mat->viscexp2 * pow((mat->visccoef * pow(strseq, mat->viscexp1))/((mat->viscexp2 + 1)*p), 1/(mat->viscexp2 + 1));
     }
     break;
+    
+    
 
-    /**** 2 parameters creep function (power_law2 in Digimat and defined in J. Li and G.J. Weng(1998), CST, 58/11, 1803-1810 ****/
+  // 2 parameters creep function (power_law2 in Digimat and defined in J. Li and G.J. Weng(1998), CST, 58/11, 1803-1810 ****/
   case 3:     
     if ((mat->visccoef<=0)||(mat->viscexp1<=0)) {   
       printf("Invalid value of the creep coefficient (must be > 0) or creep exponent1 (must be > 0) in the 2 parameters creep function\n");
@@ -1604,7 +1617,7 @@ void VEVPMFH_FP::computevisc (double strseq, double p, VEVPMFH_FP::MATPROP *mat,
 
 
 
-    /**** From Abaqus *RATE DEPENDENT, TYPE=POWERLAW  ******/
+  // From Abaqus *RATE DEPENDENT, TYPE=POWERLAW  ******/
   case 4:
     if (f<=0) {
       *gv = 0;
@@ -1619,65 +1632,84 @@ void VEVPMFH_FP::computevisc (double strseq, double p, VEVPMFH_FP::MATPROP *mat,
     break;
 
 
+
   default :
     printf("parameter choice_viscoplastic invalid in function VEVPMFH_FP::computevisc in fonctions_affine.c\n");
     exit(0);
-  }           /* End of switch */
+    
+  //End of switch  
+  }          
 
 }
 
 
 
 
-/**************************************
- ** Constructor of a VEVPMFH_FP::STATE structure **
- **************************************/
 
-/* 2003 - O. Pierard */
+//Description: Constructor of a state structure
+//Author: Mohamed Haddad 
+//Created: 07/04/2020
+//Copyright: See COPYING file that comes with this distribution 
 
 void VEVPMFH_FP::mallocstate(VEVPMFH_FP::STATE *st) {
 
+  /** Dec **/
   int i, j;
+  
+  /** Alloc. & Init. **/ 
   st->strn = (double*)calloc(6,sizeof(double));
   st->strs = (double*)calloc(6,sizeof(double));
   st->pstrn = (double*)calloc(6,sizeof(double));
   st->p = 0;
+  
+  /**** Start ****/
   VEVPMFH_FP::mallocmatrix(&st->Si, VEVPMFH_FP::max_nbrGi, 6);
   VEVPMFH_FP::mallocmatrix(&st->Sigma_Hj, VEVPMFH_FP::max_nbrKj, 6);
+  
 }
   
 
 
-/**************************************
- ** Destructor of a VEVPMFH_FP::STATE structure **
- **************************************/
 
-/* 2003 - O. Pierard */
+
+
+
+//Description: Destructor of a state structure 
+//Author: Mohamed Haddad 
+//Created: 07/04/2020
+//Copyright: See COPYING file that comes with this distribution 
 
 void VEVPMFH_FP::freestate(VEVPMFH_FP::STATE st) {
 
+  /** Dec **/
   int i;
+  
+  /**** Start ****/
   free(st.strn);
   free(st.strs);
   free(st.pstrn);
   VEVPMFH_FP::freematrix(st.Si,VEVPMFH_FP::max_nbrGi);
   VEVPMFH_FP::freematrix(st.Sigma_Hj,VEVPMFH_FP::max_nbrKj);
+
 }
 
 
 
 
-/*******************************************
- ** Copy a VEVPMFH_FP::STATE structure A into a B one **
- *******************************************/
 
-/* 2003 - O. Pierard */
 
-void VEVPMFH_FP::copystate(VEVPMFH_FP::STATE *A, VEVPMFH_FP::STATE *B)  {
+//Description: Copy a state structure A into a B one
+//Author: Mohamed Haddad 
+//Created: 07/04/2020
+//Copyright: See COPYING file that comes with this distribution 
 
+void VEVPMFH_FP::copystate(const VEVPMFH_FP::STATE *A, VEVPMFH_FP::STATE *B) {
+
+  /** Dec **/
   int i,j;
   B->p = A->p;
 
+  /**** Start ****/
   for (i=0; i<6; i++) {
     B->strn[i] = A->strn[i];
     B->strs[i] = A->strs[i];
@@ -1698,24 +1730,28 @@ void VEVPMFH_FP::copystate(VEVPMFH_FP::STATE *A, VEVPMFH_FP::STATE *B)  {
      }
    }
 
-
 }
 
 
 
 
-/***********************************************************
- ** Print contents of a VEVPMFH_FP::STATE structure to the file pfile **
- ***********************************************************/
 
-/* 2003 - O. Pierard */
+
+//Description: Print contents of a state structure to the file pfile
+//Author: Mohamed Haddad 
+//Created: 07/04/2020
+//Copyright: See COPYING file that comes with this distribution 
 
 void VEVPMFH_FP::fprintstate(FILE *pfile, VEVPMFH_FP::STATE *A) {
 
+  /** Dec **/
   int i;
   double temp;
+  
+  /** Init. **/
   temp = A->p;
 
+  /**** Start ****/
   fprintf(pfile, "Strain:");
   VEVPMFH_FP::fprintvect(pfile,A->strn,6);
   fprintf(pfile, "Stress:");
@@ -1730,18 +1766,23 @@ void VEVPMFH_FP::fprintstate(FILE *pfile, VEVPMFH_FP::STATE *A) {
 
 
 
-/*******************************************************
- ** Print contents of a VEVPMFH_FP::STATE structure to the screen **
- *******************************************************/
 
-/* June 2004 - O. Pierard */
+
+//Description: Print contents of a state structure to the screen
+//Author: Mohamed Haddad 
+//Created: 07/04/2020
+//Copyright: See COPYING file that comes with this distribution 
 
 void VEVPMFH_FP::printstate(VEVPMFH_FP::STATE *A) {
 
+  /** Dec **/
   int i;
   double temp;
 
+  /** Init. **/
   temp = A->p;
+  
+  /**** Start ****/
   printf("\nStrain:");
   VEVPMFH_FP::printvect(A->strn,6);
   printf("Stress:");
@@ -1764,34 +1805,37 @@ void VEVPMFH_FP::printstate(VEVPMFH_FP::STATE *A) {
 
 
 
-/********************************************************
- ** Print results of a the matrix phase in output file **
- ********************************************************/
 
-/* September 2022 - M. Haddad */
+//Description: Print results of a the matrix phase in output file
+//Author: Mohamed Haddad 
+//Created: September 2022
+//Copyright: See COPYING file that comes with this distribution 
 
-void fileprintresults(FILE *pfile, double time, VEVPMFH_FP::STATE *st)
-{
-fprintf(pfile, "%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e \n", time, st->strn[0], st->strn[1], st->strn[2],\
-st->strs[0], st->strs[1], st->strs[2],sqrt(2)*st->strn[3],sqrt(2)*st->strn[4], sqrt(2)*st->strn[5],\
-st->strs[3]/sqrt(2), st->strs[4]/sqrt(2),st->strs[5]/sqrt(2), st->p, VEVPMFH_FP::vonmises(st->strs), VEVPMFH_FP::eqstrn(st->strn),\
-st->pstrn[0], st->pstrn[1], st->pstrn[2], sqrt(2)*st->pstrn[3], sqrt(2)*st->pstrn[4], sqrt(2)*st->pstrn[5]); 
+void fileprintresults(FILE *pfile, double time, VEVPMFH_FP::STATE *st){
+
+  fprintf(pfile, "%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e \n", time, st->strn[0], st->strn[1], st->strn[2],\
+  st->strs[0], st->strs[1], st->strs[2],sqrt(2)*st->strn[3],sqrt(2)*st->strn[4], sqrt(2)*st->strn[5],\
+  st->strs[3]/sqrt(2), st->strs[4]/sqrt(2),st->strs[5]/sqrt(2), st->p, VEVPMFH_FP::vonmises(st->strs), VEVPMFH_FP::eqstrn(st->strn),\
+  st->pstrn[0], st->pstrn[1], st->pstrn[2], sqrt(2)*st->pstrn[3], sqrt(2)*st->pstrn[4], sqrt(2)*st->pstrn[5]); 
+  
 }
 
 
 
 
 
-/*****************************************************
- ** Print contents of a VEVPMFH_FP::MATPROP structure in a file **
- *****************************************************/
 
-/* 2020 - M. Haddad */
+//Description: Print contents of a MATPROP structure in a file
+//Author: Mohamed Haddad 
+//Created: 08/04/2020
+//Copyright: See COPYING file that comes with this distribution 
 
 void VEVPMFH_FP::fprintmatprop(FILE *pfile, VEVPMFH_FP::MATPROP *A) {
 
+  /** Dec **/
   int i;
   
+  /**** Start ****/
   fprintf(pfile, "\n E(t=0): %e\n",A->E0);
   fprintf(pfile, "nu: %g\n",A->nu);
   fprintf(pfile, "yldstrs: %e\n",A->yldstrs);
@@ -1812,43 +1856,47 @@ void VEVPMFH_FP::fprintmatprop(FILE *pfile, VEVPMFH_FP::MATPROP *A) {
   fprintf(pfile, "number of shear relaxation times: %d\n",A->nbrGi);
   fprintf(pfile, "Gi:");
   for(i=0;i<A->nbrGi;i++){
-   fprintf(pfile, "%e\t",A->Gi[i]);
-   fprintf(pfile, "\n");
-   }
+    fprintf(pfile, "%e\t",A->Gi[i]);
+    fprintf(pfile, "\n");
+  }
 
   fprintf(pfile, "gi:");
   for(i=0;i<A->nbrGi;i++){
-   fprintf(pfile, "%e\t",A->gi[i]);
-   fprintf(pfile, "\n");
-   }
+    fprintf(pfile, "%e\t",A->gi[i]);
+    fprintf(pfile, "\n");
+  }
 
   fprintf(pfile, "K0: %e\n",A->K0);
   fprintf(pfile, "K_inf: %e\n",A->K_inf);
   fprintf(pfile, "number of bulk relaxation times: %d\n",A->nbrKj);
   fprintf(pfile, "Kj:");
   for(i=0;i<A->nbrKj;i++){
-   fprintf(pfile, "%e\t",A->Kj[i]);
-   fprintf(pfile, "\n");
-   }
+    fprintf(pfile, "%e\t",A->Kj[i]);
+    fprintf(pfile, "\n");
+  }
 
   fprintf(pfile, "kj:");
   for(i=0;i<A->nbrKj;i++){
-   fprintf(pfile, "%e\t",A->kaj[i]);
-   fprintf(pfile, "\n");
-   }
+    fprintf(pfile, "%e\t",A->kaj[i]);
+    fprintf(pfile, "\n");
+  }
+  
 }
 
 
 
 
-/*********************************************************
- ** Print contents of a VEVPMFH_FP::MATPROP structure on the screen **
- *********************************************************/
 
-/* 2020 - M. Haddad */
+
+
+//Description: Print contents of a MATPROP structure on the screen
+//Author: Mohamed Haddad 
+//Created: 08/04/2020
+//Copyright: See COPYING file that comes with this distribution 
 
 void VEVPMFH_FP::printmatprop(VEVPMFH_FP::MATPROP *A) {
 
+  /**** Start ****/
   printf("\n E(t=0): %g\n",A->E0);
   printf("nu: %g\n",A->nu);
   printf("yldstrs: %g\n",A->yldstrs);
@@ -1885,31 +1933,40 @@ void VEVPMFH_FP::printmatprop(VEVPMFH_FP::MATPROP *A) {
 
 
 
-/**	COMPUTE G_tilde(dt)	**/
-/*********************************/
-/** November 2020 - M.HADDAD	**/
 
 
-double VEVPMFH_FP::compute_Gtilde(VEVPMFH_FP::MATPROP *material_prop, double dt){
+//Description: COMPUTE G_tilde(dt) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
+
+double VEVPMFH_FP::compute_Gtilde(const VEVPMFH_FP::MATPROP *material_prop, double dt){
 
+  /** Dec. & Init. **/ 
   int i;
   double result=material_prop->G_inf;
 
+  /**** Start ****/
   for(i=0;i<material_prop->nbrGi;i++){
-    /*printf("part %d = %g \n", i, material_prop->Gi[i] * (1 - exp(-dt / (material_prop->gi[i])) ) * (dt / material_prop->gi[i]));*/
     result += material_prop->Gi[i] * (1 - exp(-dt / (material_prop->gi[i])) ) * (material_prop->gi[i] / dt);
   }
 
   return(result);
+  
 }
 
 
 
 
-/**	COMPUTE K_tilde(dt)	**/
-/*********************************/
-/** November 2020 - M.HADDAD	**/
-double VEVPMFH_FP::compute_Ktilde(VEVPMFH_FP::MATPROP *material_prop, double dt){
+
+
+
+//Description: COMPUTE K_tilde(dt) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
+
+double VEVPMFH_FP::compute_Ktilde(const VEVPMFH_FP::MATPROP *material_prop, double dt){
 
   int i;
   double result=material_prop->K_inf;
@@ -1924,26 +1981,33 @@ double VEVPMFH_FP::compute_Ktilde(VEVPMFH_FP::MATPROP *material_prop, double dt)
 
 
 
-/**	COMPUTE E_tilde(dt)	**/
-/*********************************/
-/** November 2020 - M.HADDAD	**/
 
-void VEVPMFH_FP::compute_Etilde(VEVPMFH_FP::MATPROP *material_prop, double dt, double **E_tilde){
+
+//Description: COMPUTE E_tilde(dt) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
+
+void VEVPMFH_FP::compute_Etilde(const VEVPMFH_FP::MATPROP *material_prop, double dt, double **E_tilde){
   
+  /** Dec. **/
   double Gtilde, Ktilde;
   double **Idev, **Ivol;
 
-  Gtilde = compute_Gtilde(material_prop,dt);
-  Ktilde = compute_Ktilde(material_prop,dt);
-
+  /** Alloc. & Init. **/
   VEVPMFH_FP::mallocmatrix(&Idev,6,6);
   VEVPMFH_FP::mallocmatrix(&Ivol,6,6);
+  
+  /**** Start ****/
+  Gtilde = compute_Gtilde(material_prop,dt);
+  Ktilde = compute_Ktilde(material_prop,dt);
 
   VEVPMFH_FP::compute_Idev(Idev);
   VEVPMFH_FP::compute_Ivol(Ivol);
 
   VEVPMFH_FP::addtens4(2*Gtilde,Idev,3*Ktilde,Ivol,E_tilde);
-
+  
+  /** Freeing **/
   VEVPMFH_FP::freematrix(Idev,6);
   VEVPMFH_FP::freematrix(Ivol,6); 
 
@@ -1952,11 +2016,14 @@ void VEVPMFH_FP::compute_Etilde(VEVPMFH_FP::MATPROP *material_prop, double dt, d
 
 
 
-/**	COMPUTE E_inf(dt)	**/
-/*********************************/
-/** November 2020 - M.HADDAD	**/
 
-void VEVPMFH_FP::compute_Einf(VEVPMFH_FP::MATPROP *material_prop, double **E_inf){
+
+//Description: COMPUTE E_inf(dt) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
+
+void VEVPMFH_FP::compute_Einf(const VEVPMFH_FP::MATPROP *material_prop, double **E_inf){
 
   double **Idev, **Ivol;
 
@@ -1977,55 +2044,67 @@ void VEVPMFH_FP::compute_Einf(VEVPMFH_FP::MATPROP *material_prop, double **E_inf
 
 
 
-/**	DECOMPOSITION OF VE STRAIN TO VOL AND DEV PARTS 		**/
-/*************************************************************************/
+
+//Description: DECOMPOSITION OF VE STRAIN TO VOL AND DEV PARTS (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
 
 void VEVPMFH_FP::decomposition_dstrn_ve(double *dstrn_ve, double *dstrn_ve_volpart, double *dstrn_ve_devpart){
 
+  /** Dec. **/
   double **Idev, **Ivol;
 
+  /** Alloc. and Init. **/ 
   VEVPMFH_FP::mallocmatrix(&Idev,6,6);
   VEVPMFH_FP::mallocmatrix(&Ivol,6,6);
 
+  /**** Start ****/
   VEVPMFH_FP::compute_Ivol(Ivol);
   VEVPMFH_FP::compute_Idev(Idev);
  
   VEVPMFH_FP::contraction42(Ivol, dstrn_ve, dstrn_ve_volpart);
   VEVPMFH_FP::contraction42(Idev, dstrn_ve, dstrn_ve_devpart);
 
+  /** Freeing **/
   VEVPMFH_FP::freematrix(Idev,6);
   VEVPMFH_FP::freematrix(Ivol,6); 
 
-
 }
 
 
 
 
 
-/**	COMPUTE Si(tn+1)		**/
-/*****************************************/
-/** November 2020 - M.HADDAD	**/
 
+//Description: COMPUTE Si(tn+1) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
 
-void VEVPMFH_FP::compute_Si(VEVPMFH_FP::MATPROP *material_prop, VEVPMFH_FP::STATE *tn, double dt, double *dstrn_ve, double **Si){
+void VEVPMFH_FP::compute_Si(const VEVPMFH_FP::MATPROP *material_prop, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn_ve, double **Si){
 
+   /** Dec. **/
    int i,j;
    double *dstrn_ve_vol, *dstrn_ve_dev;
 
+   /** Alloc. **/
    dstrn_ve_vol= (double*) calloc(6,sizeof(double));
    dstrn_ve_dev= (double*) calloc(6,sizeof(double));
 
+   /**** Start ****/
+   // comute dev and vol parts of the strain 
    VEVPMFH_FP::decomposition_dstrn_ve(dstrn_ve,dstrn_ve_vol,dstrn_ve_dev);
 
+   // Si_n+1 = exp(-dt/gi)*Si_n + (1 - exp(-dt/gi)) * (gi/dt) * (dstrn_ve)_dev 
    for(i=0;i<material_prop->nbrGi;i++){
-
       for(j=0;j<6;j++){
-      Si[i][j] = ( exp(-dt / material_prop->gi[i]) * tn->Si[i][j] ) + ( 2 * material_prop->Gi[i] * (1 - exp(-dt / (material_prop->gi[i])) ) * (material_prop->gi[i] / dt) * dstrn_ve_dev[j]);
+        Si[i][j] = ( exp(-dt / material_prop->gi[i]) * tn->Si[i][j] ) + ( 2 * material_prop->Gi[i] * (1 - exp(-dt / (material_prop->gi[i])) ) * (material_prop->gi[i] / dt) * dstrn_ve_dev[j]);
       }
 
-    }
+   }
 
+  /** Freeing **/
   free(dstrn_ve_vol);
   free(dstrn_ve_dev); 
    
@@ -2035,54 +2114,63 @@ void VEVPMFH_FP::compute_Si(VEVPMFH_FP::MATPROP *material_prop, VEVPMFH_FP::STAT
 
 
 
-/**	COMPUTE Sigma_Hj(tn+1)		**/
-/*****************************************/
-/** November 2020 - M.HADDAD	**/
 
-void VEVPMFH_FP::compute_SigmaHj(VEVPMFH_FP::MATPROP *material_prop, VEVPMFH_FP::STATE *tn, double dt, double *dstrn_ve, double **Sigma_Hj){
+//Description: COMPUTE Sigma_Hj(tn+1) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
+
+void VEVPMFH_FP::compute_SigmaHj(const VEVPMFH_FP::MATPROP *material_prop, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn_ve, double **Sigma_Hj){
 
+   /** Dec. **/
    int i,j;
    double *dstrn_ve_vol, *dstrn_ve_dev;
 
+   /** Alloc. **/
    dstrn_ve_vol= (double*) calloc(6,sizeof(double));
    dstrn_ve_dev= (double*) calloc(6,sizeof(double));
 
+   /**** Start ****/
    VEVPMFH_FP::decomposition_dstrn_ve(dstrn_ve,dstrn_ve_vol,dstrn_ve_dev);
 
    for(i=0;i<material_prop->nbrKj;i++){
-
       for(j=0;j<6;j++){
-       Sigma_Hj[i][j] = ( exp(-dt / (material_prop->kaj[i])) * tn->Sigma_Hj[i][j] ) + ( 3 * material_prop->Kj[i] * (1 - exp(-dt / (material_prop->kaj[i]))) * (material_prop->kaj[i] / dt) * dstrn_ve_vol[j]);
-       }
-
-    }
+         Sigma_Hj[i][j] = ( exp(-dt / (material_prop->kaj[i])) * tn->Sigma_Hj[i][j] ) + ( 3 * material_prop->Kj[i] * (1 - exp(-dt / (material_prop->kaj[i]))) * (material_prop->kaj[i] / dt) * dstrn_ve_vol[j]);
+      }
+   }
 
+  /** Freeing **/
   free(dstrn_ve_vol);
-  free(dstrn_ve_dev); 
+  free(dstrn_ve_dev);
+   
 }
 
 
 
-/**	COMPUTE sum of Si(tn)*exp(-dt/gi)		**/
-/*********************************************************/
-/** November 2020 - M.HADDAD	**/
 
 
-void VEVPMFH_FP::compute_sum_Siexp(VEVPMFH_FP::MATPROP *material_prop, double dt, double **Si, double *sum_Siexp){
 
-  int i,j;
 
-  for(j=0;j<6;j++) 
-   sum_Siexp[j] = 0;
+//Description: COMPUTE sum of Si(tn)*exp(-dt/gi) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
 
+void VEVPMFH_FP::compute_sum_Siexp(const VEVPMFH_FP::MATPROP *material_prop, double dt, double **Si, double *sum_Siexp){
   
-   for(j=0;j<6;j++){
+  /** Dec. **/
+  int i,j;
 
-     for(i=0;i<material_prop->nbrGi;i++){
-       sum_Siexp[j] += Si[i][j] * exp( -dt / (material_prop->gi[i])) ;
+  /** Init. **/
+  for(j=0;j<6;j++){sum_Siexp[j] = 0;}
 
-     }
+  
+  /**** Start ****/
+  for(j=0;j<6;j++){
+    for(i=0;i<material_prop->nbrGi;i++){
+       sum_Siexp[j] = sum_Siexp[j] + (Si[i][j] * exp( -dt / (material_prop->gi[i]))) ;
     }
+  }
 
 
 }
@@ -2090,35 +2178,43 @@ void VEVPMFH_FP::compute_sum_Siexp(VEVPMFH_FP::MATPROP *material_prop, double dt
 
 
 
-/**	COMPUTE sum of SigmaHj(tn)*exp(-dt/kj)		**/
-/*********************************************************/
-/** November 2020 - M.HADDAD	**/
-void VEVPMFH_FP::compute_sum_SigmaHjexp(VEVPMFH_FP::MATPROP *material_prop, double dt, double **Sigma_Hj, double *sum_SigmaHjexp){
 
- int i,j;
 
-  for(j=0;j<6;j++) 
-   sum_SigmaHjexp[j] = 0;
+//Description: COMPUTE sum of SigmaHj(tn)*exp(-dt/kj) (see B.Miled & I.Doghri, 2011, IJSS.)
+//Author: Mohamed Haddad 
+//Created: November 2020
+//Copyright: See COPYING file that comes with this distribution 
 
+void VEVPMFH_FP::compute_sum_SigmaHjexp(const VEVPMFH_FP::MATPROP *material_prop, double dt, double **Sigma_Hj, double *sum_SigmaHjexp){
 
-   for(j=0;j<6;j++){
-   
-     for(i=0;i<material_prop->nbrKj;i++){
+  /** Dec. **/
+  int i,j;
 
-       sum_SigmaHjexp[j] += Sigma_Hj[i][j] * exp( -dt / (material_prop->kaj[i])) ;
+  /** Init. **/
+  for(j=0;j<6;j++){sum_SigmaHjexp[j] = 0;}
 
+  /**** Starting ****/
+  for(j=0;j<6;j++){
+     for(i=0;i<material_prop->nbrKj;i++){
+       sum_SigmaHjexp[j] = sum_SigmaHjexp[j]  + (Sigma_Hj[i][j] * exp( -dt / (material_prop->kaj[i]))) ;
      }
-    }
+  }
 
 }
 
 
 
 
-/** COMPUTE Calgo USING PERTURBATION METHOD	      **/
-/******************************************************/
 
-void CalgoUsingPerturbation(VEVPMFH_FP::MATPROP *matrprop, VEVPMFH_FP::STATE *matrtn, VEVPMFH_FP::STATE *matrtau, double dt, double *dstrn, double **Calgo){
+
+//Description: COMPUTE Calgo USING PERTURBATION METHOD
+//Author: Mohamed Haddad 
+//Created: September 2022
+//Copyright: See COPYING file that comes with this distribution 
+
+void CalgoUsingPerturbation(const VEVPMFH_FP::MATPROP *matrprop, VEVPMFH_FP::STATE *matrtn, VEVPMFH_FP::STATE *matrtau, double dt, double *dstrn, double **Calgo){
+
+  /** Dec. **/
   VEVPMFH_FP::STATE matrtaupert1, matrtaupert2;
   
   int cvarstrn, i;
@@ -2126,6 +2222,7 @@ void CalgoUsingPerturbation(VEVPMFH_FP::MATPROP *matrprop, VEVPMFH_FP::STATE *ma
   double *vardstrn, *dstrnpert1, *dstrnpert2;
   double **Calgopert1, **Calgopert2;
   
+  /** Alloc. & Init. **/
   VEVPMFH_FP::mallocstate(&matrtaupert1);
   VEVPMFH_FP::mallocstate(&matrtaupert2);
   
@@ -2135,28 +2232,32 @@ void CalgoUsingPerturbation(VEVPMFH_FP::MATPROP *matrprop, VEVPMFH_FP::STATE *ma
   VEVPMFH_FP::mallocmatrix(&Calgopert1,6,6);
   VEVPMFH_FP::mallocmatrix(&Calgopert2,6,6);
   
-
-  
-  for (cvarstrn=0; cvarstrn<6; cvarstrn++) { /* Computation of the tangent operator via a perturbation method*/
+  /**** Starting ****/
+  for (cvarstrn=0; cvarstrn<6; cvarstrn++) { 
     for (i=0; i<6; i++) vardstrn[i]=0;
+    // Compute the perturation 
     vardstrn[cvarstrn] = dstrn[cvarstrn]*1.e-3 + dstrn[0]*1.e-5 + 1.e-10;  /* Last term to get a perturbation even if nil component*/ 
     
-    addtens2(1, dstrn, 1, vardstrn, dstrnpert1);                            /* Perturbated strain increment */
+    // Perurbation of dstrn increment 
+    addtens2(1, dstrn, 1, vardstrn, dstrnpert1);                            
+    // Compute the new mechanical state after perturbation 
     ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(matrprop, matrtn, dt, dstrnpert1, &matrtaupert1, Calgopert1);
     
-    addtens2(1, dstrn, -1, vardstrn, dstrnpert1);                            /* Perturbated strain increment */
+    addtens2(1, dstrn, -1, vardstrn, dstrnpert1);                            
     ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(matrprop, matrtn, dt, dstrnpert2, &matrtaupert2, Calgopert2);
 
+    // C = d dstrs/ d dstrn
     for (i=0;i<6;i++) {
-      /*Calgo[i][cvarstrn] = ((matrtaupert1.strs)[i]-(matrtaupert2.strs)[i])/(2*vardstrn[cvarstrn]);*/
-      Calgo[i][cvarstrn] = ((matrtaupert1.strs)[i]-(matrtau->strs)[i])/vardstrn[cvarstrn];
+      Calgo[i][cvarstrn] = ((matrtaupert1.strs)[i]-(matrtaupert2.strs)[i])/(2*vardstrn[cvarstrn]);
     }
 
-  }   /** Tangent operator corresponding to the macro strain increment found */
-
+  }   
+  
+  /** Freeing **/
+  VEVPMFH_FP::freestate(matrtaupert1); VEVPMFH_FP::freestate(matrtaupert2);
   free(vardstrn); free(dstrnpert1); free(dstrnpert2);
   VEVPMFH_FP::freematrix(Calgopert1,6);VEVPMFH_FP::freematrix(Calgopert2,6); 
-  VEVPMFH_FP::freestate(matrtaupert1); VEVPMFH_FP::freestate(matrtaupert2);
+  
 }
 
 
diff --git a/NonLinearSolver/materialLaw/fonctions_pratiques.h b/NonLinearSolver/materialLaw/fonctions_pratiques.h
index f3e9d46e65b198a37e51935e15334e7291790aa7..746840ce5a43e4110d13a27810568c42ed758345 100644
--- a/NonLinearSolver/materialLaw/fonctions_pratiques.h
+++ b/NonLinearSolver/materialLaw/fonctions_pratiques.h
@@ -1,251 +1,157 @@
 #ifndef FONCTIONS_PRATIQUES_H
 #define FONCTIONS_PRATIQUES_H 1
 
-namespace VEVPMFH_FP{
-
-const double max_nbrGi = 4;
-const double max_nbrKj = 4;
-
-/* Convertit un tenseur d'ordre 2 en un vecteur 6*1 */
-void convertmatrvect (double **matr, double *vect);
-
-/* Convertit un tenseur sous forme 6*1 en une matrice 3*3 */
-void convertvectmatr (double *vect, double **matr);
-
-/* Addition of two tensors of the fourth order, each one can be previously multiplied by a scalar */
-void addtens4 (double a, double **A, double b, double **B, double **res);
-
-/* Addition of two tensors of the second order, each one can be previously multiplied by a scalar */
-void addtens2 (double a, double *A, double b, double *B, double *res);
-
-/* Contraction de deux tenseurs d'ordre 2 mis sous forme 6*1 */
-double contraction22 (const double *, const double *);
-
-/* Contraction d'un tenseur d'ordre4(6*6) et d'ordre2(6*1) */
-void contraction42 (double **T4, double *T2, double *Sol);
-
-/* Contraction de deux tenseurs d'ordre 4(6*6) sur deux indices */
-void contraction44 (double **, double **, double **);
-
-/* Contraction de deux tenseurs d'ordre 4(6*6) sur quatre indices */
-double dcontr44 (double **TA, double **TB);
-
-/* Double contraction of three fourth order tensors; last contraction being over four indices */
-double contraction444 (double **A, double **B, double **C);
-
-/* Double contraction of a two fourth order tensors and one of the second order */
-void contraction442 (double **A, double **B, double *C, double *Tres);
-
-/* Contraction of three 4th order tensor - same two first indices of the first tensor - idem for the two last ones of the last tensor */
-double contr444_reduced (double **A, double **B, double **C);
-
-/* Produit vectoriel de deux tenseurs d'ordre2(6*1) */
-void produit_tensoriel (double *TA, double *TB, double **Tres);
-
-/* Compute the rotation matrix q which transform a Sd order tensor from the local to the global system of axes */
-void compute_q (double phi, double theta, double psi, double **q);
-
-/* Changement de repere d'un tenseur d'ordre 2 du repere global au repere local */
-void rotatetens2gloloc(double *, double, double, double, double *);
-
-/* Changement de repere d'un tenseur d'ordre 2 du repere local au repere global */
-void rotatetens2locglo (double *tenseur, double phi, double theta, double psi, double *tenseur_rot);
-
-/* Changement de repere d'un tenseur d'ordre 4 du repere global au repere local */
-void rotatetens4gloloc (double **, double, double, double, double **);
-
-/* Changement de repere d'un tenseur d'ordre 4 du repere local au repere global */
-void rotatetens4locglo (double **tenseur, double phi, double theta, double psi, double **tenseur_rot);
-
-/* Factorisatin LU d'une matrice d'ordre n */
-void ludcmp(double **a, int n, int *indx, double *d, int *error);
-
-/* Forward and backward substitution */
-void lubksb (double **a, int n, int *indx, double b[]);
-
-/* Computes the inverse of a matrix of size n */
-void inverse (double **a, double **y, int *error, int n);
-
-/* Extract isotropic part of a fourth order tensor */
-void extractiso (double **tensani, double **tensiso);
-
-/* Memory allocation of a table of pointers */
-void mallocmatrix (double ***A, int m, int n);
-
-/* Freeing memory of a table of pointers */
-void freematrix (double **A, int m);
-
-/* Memory allocation for a 3D matrix */
-void mallocmatrix3D (double ****A, int p, int m, int n);
-
-/* Freeing memory of a 3D matrix */
-void freematrix3D (double ***A, int p, int m);
-
-/* Compute the transverse isotropic Hooke's elastic operator from E and nu */
-void compute_Hooke(double E, double nu, double **C);
-
-/* Print on the screen values of a vector */
-void printvect(double *vect,int  n);
-
-/* Print values of a vector to a specified file */
-void fprintvect(FILE *pfile, double *vect, int n);
 
-/* Print on the screen values of a m*n matrix */
-void printmatr(double **matr4, int m, int n);
 
-/* Print values of a m*n matrix to a specified file */
-void fprintmatr(FILE *pfile, double **matr4, int m, int n);
 
-/* Compute the trace of a second order tensor */
-double trace(double *tens2);
-
-/* Extract the deviatoric part of a second order tensor */
-void dev(double *tens2, double *dev_tens2);
-
-/* Compute the von Mises equivalent stress */
-double vonmises(double *tens2);
-
-/* Compute the equivalent strain scalar */
-double eqstrn(double *tens2);
-
-/* Compare two vectors */
-int compare(double *A, double *B, int size);
-
-/* Compare two matrices */
-int comparematr(double **A, double **B, int m, int n);
-
-/* Check if a vector is null */
-int isnull(double *A, int size);
-
-/* Check if a fourth order tensor is isotropic */
-int isisotropic(double**A);
-
-/* Copy a vector into another */
-void copyvect(double *A, double *B, int size);
-
-/* Find the maximum value of a vector */
-double maxvect(double *vect, int n);
-
-/* Find the maximum value in absolute value of a vector */
-double maxvectabs(double *vect, int n);
-
-/* Find the minimum non nil value in absolute value of a vector */
-double minvectabs(double *vect, int m);
-
-/* Find the smallest non-nil value in absolute value of a matrix */
-double minmatrabs(double **matr, int m, int n);
-
-/* Find max of two values */
-double max(double a, double b);
-
-
-/* Compute I */
-void compute_I(double **I);
+namespace VEVPMFH_FP{
 
-/* Compute Ivol */
-void compute_Ivol(double **I_vol);
+  const int max_nbrGi = 4;  
+  const int max_nbrKj = 4;
+  
+  // Structure containing the mechanical variables
+  typedef struct State {
+    double *strn;  		/* Strain tensor */
+    double *strs;  		/* Stress tensor */
+    double p;      		/* Accumulated plasticity - scalar */
+    double *pstrn;  		/* Plastic strain */
+    double **Si;        	/* matrix nbrGi*6 : deviatoric part of viscous stress relative to Gi */
+    double **Sigma_Hj;  	/* matrix nbrGi*6 : volumetric part of viscous stress relative to kj */
+  } STATE;
+
+
+
+
+  // Structure containing the macro mechanical variables (Only for MFH)
+  typedef struct Statemacro {
+    double *strn;       	 /* Macroscopic strain tensor */
+    double *strs;        	 /* Macroscopic stress tensor */
+    double **Cmacro;		 /* Macroscopic  elastic stiffness (or equivalent) */
+    double *Beta_macro;    	 /* Macroscopic thermal stress (or equivalent) */
+  } STATEMACRO;
+
+
+
+
+  // Structure containing the VE-VP material properties
+  typedef struct Matprop{ 
+    double E0;
+    double nu;
+    double yldstrs;
+    int hardmodel;
+    double hardmodulus1;
+    double hardmodulus2;
+    double hardexp;
+    int viscmodel;
+    double visccoef;
+    double viscexp1;
+    double viscexp2;
+    double v;
+    double ar;
+    double orient[3];
+    
+ 
+    // VE properties 
+    double G0;
+    double G_inf;
+    int nbrGi;  		/* number of shear relaxation times */
+    double *Gi; 		/* shear relaxation weights */ 
+    double *gi; 		/* shear relaxation times */
+  
+    double K0;
+    double K_inf;
+    int nbrKj;   		/* number of bulk relaxation times */ 
+    double *Kj;  		/* bulk relaxation weights */
+    double *kaj; 		/* bulk relaxation times */ 
+  }MATPROP;
+
+
+  void convertmatrvect (double **matr, double *vect);
+  void convertvectmatr (double *vect, double **matr);
+
+  void addtens4 (double a, double **A, double b, double **B, double **res);
+  void addtens2 (double a, double *A, double b, double *B, double *res);
+  double contraction22 (const double *, const double *);
+  void contraction42 (double **T4, double *T2, double *Sol);
+  void contraction44 (double **, double **, double **);
+  double dcontr44 (double **TA, double **TB);
+  double contraction444 (double **A, double **B, double **C);
+  void contraction442 (double **A, double **B, double *C, double *Tres);
+  double contr444_reduced (double **A, double **B, double **C);
+  void produit_tensoriel (double *TA, double *TB, double **Tres);
+
+  void ludcmp(double **a, int n, int *indx, double *d, int *error);
+  void lubksb (double **a, int n, int *indx, double b[]);
+  void inverse (double **a, double **y, int *error, int n);
+  void extractiso (double **tensani, double **tensiso);
+
+  void mallocmatrix (double ***A, int m, int n);
+  void freematrix (double **A, int m);
+  void mallocmatrix3D (double ****A, int p, int m, int n);
+  void freematrix3D (double ***A, int p, int m);
+  void mallocMATPROP(VEVPMFH_FP::MATPROP* matr, VEVPMFH_FP::MATPROP* incl, const char* inputfilename);
+  
+  void compute_Hooke(double E, double nu, double **C);
 
-/* Compute Idev */
-void compute_Idev(double **I_dev);
+  void printvect(double *vect,int  n);
+  void fprintvect(FILE *pfile, double *vect, int n);
+  void printmatr(double **matr4, int m, int n);
+  void fprintmatr(FILE *pfile, double **matr4, int m, int n);
 
-/* Compute zero vector of size size **/
-void zero_vect(double size, double *vect);
+  double max(double a, double b);
+  int isnull(double *A, int size);
 
+  double trace(double *tens2);
+  void dev(double *tens2, double *dev_tens2);
+  double vonmises(double *tens2);
+  double eqstrn(double *tens2);
 
-/* Structure containing the mechanical variables */
-typedef struct State {
-  double *strn;  /* Strain tensor */
-  double *strs;  /* Stress tensor */
-  double p;      /* Accumulated plasticity - scalar */
-  double *pstrn;  /* Plastic strain */
-  double **Si;        /* matrix nbrGi*6 : deviatoric part of viscous stress relative to Gi */
-  double **Sigma_Hj;  /* matrix nbrGi*6 : volumetric part of viscous stress relative to kj */
-} STATE;
+  int isisotropic(double**A);
+  void copyvect(double *A, double *B, int size);
 
+  void compute_I(double **I);
+  void compute_Ivol(double **I_vol);
+  void compute_Idev(double **I_dev);
 
 
+  void readinput (MATPROP *matr, MATPROP *incl, const char* inputfilename);
+  void computestrslinel(double young, double nu, double *strn, double *strs);
+  void computestrnlinel(double young, double nu, double *strs, double *strn);
 
-/* Structure containing the macro mechanical variables (Only for MFH) */
-typedef struct Statemacro {
-  double *strn;        /* Macroscopic strain tensor */
-  double *strs;        /* Macroscopic stress tensor */
-  double **Cmacro;
-  double *Beta_macro;   /* Macroscopic eigenstrain tensor */
-} STATEMACRO;
+  void compute_R (double p, const MATPROP *mat, double *Rp, double *dRdp);
+  void computevisc (double strseq, double p, const MATPROP *mat, double *gv, double *dgvdstrseq, double *dgvdp);
 
+  void mallocstate(STATE *st);
+  void freestate(STATE st);
+  void copystate(const STATE *A, STATE *B);
+  void fprintstate(FILE *pfile, STATE *A);
+  void printstate(STATE *A);
 
+  void fileprintresults(FILE *pfile, double time, STATE *st);
 
+  void printmatprop(MATPROP *A);
+  void fprintmatprop(FILE *pfile, MATPROP *A);
 
-/* Structure containing the VE-VP material properties */
-typedef struct Matprop
-{ 
-  double E0;
-  double nu;
-  double yldstrs;
-  int hardmodel;
-  double hardmodulus1;
-  double hardmodulus2;
-  double hardexp;
-  int viscmodel;
-  double visccoef;
-  double viscexp1;
-  double viscexp2;
-  double v;
-  double ar;
-  double orient[3];
-  /***********************
-  **	VE parameters	**
-  ***********************/
-  double G0;
-  double G_inf;
-  int nbrGi;  /* number of shear relaxation times */
-  double *Gi; /* shear relaxation weights */ 
-  double *gi; /* shear relaxation times */
+  double compute_Gtilde(const VEVPMFH_FP::MATPROP *material_prop, double dt);
+  double compute_Ktilde(const VEVPMFH_FP::MATPROP *material_prop, double dt);
+  void compute_Etilde(const VEVPMFH_FP::MATPROP *material_prop, double dt, double **E_tilde);
+  void compute_Einf(const VEVPMFH_FP::MATPROP *material_prop, double **E_inf);
   
-  double K0;
-  double K_inf;
-  int nbrKj;   /* number of bulk relaxation times */ 
-  double *Kj;  /* bulk relaxation weights */
-  double *kaj; /* bulk relaxation times */ 
-}MATPROP;
-
-void readinput (MATPROP *matr, MATPROP *incl, const char* inputfilename);
-void computestrslinel(double young, double nu, double *strn, double *strs);
-void computestrnlinel(double young, double nu, double *strs, double *strn);
-
-void compute_R (double p, MATPROP *mat, double *Rp, double *dRdp);
-void computevisc (double strseq, double p, MATPROP *mat, double *gv, double *dgvdstrseq, double *dgvdp);
-
-void mallocstate(STATE *st);
-void freestate(STATE st);
-void copystate(STATE *A, STATE *B);
-void fprintstate(FILE *pfile, STATE *A);
-void printstate(STATE *A);
-
-
-void fileprintresults(FILE *pfile, double time, STATE *st);
-
-void printmatprop(MATPROP *A);
-void fprintmatprop(FILE *pfile, MATPROP *A);
-
-
-
+  void decomposition_dstrn_ve(double *dstrn_ve, double *dstrn_ve_volpart, double *dstrn_ve_devpart);
+  
+  void compute_Si(const VEVPMFH_FP::MATPROP *material_prop, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn_ve, double **Si);
+  void compute_SigmaHj(const VEVPMFH_FP::MATPROP *material_prop, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn_ve, double **Sigma_Hj);
 
-double compute_Gtilde(MATPROP *material_prop, double dt);
-double compute_Ktilde(MATPROP *material_prop, double dt);
-void compute_Etilde(MATPROP *material_prop, double dt, double **E_tilde);
-void compute_Einf(MATPROP *material_prop, double **E_inf);
-void decomposition_dstrn_ve(double *dstrn_ve, double *dstrn_ve_volpart, double *dstrn_ve_devpart);
-void compute_Si(MATPROP *material_prop, STATE *tn, double dt, double *dstrn_ve, double **Si);
-void compute_SigmaHj(MATPROP *material_prop, STATE *tn, double dt, double *dstrn_ve, double **Sigma_Hj);
+  void compute_sum_Siexp(const VEVPMFH_FP::MATPROP *material_prop, double dt, double **Si, double *sum_Siexp);
+  void compute_sum_SigmaHjexp(const VEVPMFH_FP::MATPROP *material_prop, double dt, double **Sigma_Hj, double *sum_SigmaHjexp);
+  
+  //void compute_sum_Si2exp(MATPROP *material_prop, double dt, double dt_unload, double **Si, double *sum_Si2exp);
+  //void compute_sum_SigmaHj2exp(MATPROP *material_prop, double dt, double dt_unload, double **Sigma_Hj, double *sum_SigmaHj2exp);
 
-void compute_sum_Siexp(MATPROP *material_prop, double dt, double **Si, double *sum_Siexp);
-void compute_sum_SigmaHjexp(MATPROP *material_prop, double dt, double **Sigma_Hj, double *sum_SigmaHjexp);
-void compute_sum_Si2exp(MATPROP *material_prop, double dt, double dt_unload, double **Si, double *sum_Si2exp);
-void compute_sum_SigmaHj2exp(MATPROP *material_prop, double dt, double dt_unload, double **Sigma_Hj, double *sum_SigmaHj2exp);
+  void CalgoUsingPerturbation(const MATPROP *matrprop, STATE *matrtn, STATE *matrtau, double dt, double *dstrn, double **Calgo);
 
-void CalgoUsingPerturbation(MATPROP *matrprop, STATE *matrtn, STATE *matrtau, double dt, double *dstrn, double **Calgo);
+/** End namespace **/
 }
 
- /** End namespace **/
+ 
 #endif
diff --git a/NonLinearSolver/materialLaw/mlawVEVPMFH.cpp b/NonLinearSolver/materialLaw/mlawVEVPMFH.cpp
index 0b629acfa81d72f96fbcf9956dda41c2589ab5f9..76ed95bba2d0f0ab946822d300d05ad83a919e74 100644
--- a/NonLinearSolver/materialLaw/mlawVEVPMFH.cpp
+++ b/NonLinearSolver/materialLaw/mlawVEVPMFH.cpp
@@ -16,129 +16,146 @@
 #include "STensorOperations.h"
 #include "fonctions_pratiques.h"
 #include "box_secant.h"
-//#include "nonLinearMechSolver.h"
+#include <cstring>
 
-mlawVEVPMFH::mlawVEVPMFH(const int num, const double rho,
-                   const char *propName) : materialLaw(num,true), _rho(rho)
-  {
+#include "nonLinearMechSolver.h"
+
+mlawVEVPMFH::mlawVEVPMFH(const int num, const double rho, const char *propName) : materialLaw(num,true), _rho(rho){
         
         std::ifstream in(propName);
-        if(!in) Msg::Error("Cannot open the file %s! Maybe is missing   ",propName);
-
-       //fill propmat
-        VEVPMFH_FP::readinput (matrprop, inclprop, propName);
-  }
+          if(!in) Msg::Error("Cannot open the file %s! Maybe is missing   ",propName);
+        
+        inputfilename = propName;
+        VEVPMFH_FP::readinput(&matrprop, &inclprop, propName);
+        /*VEVPMFH_FP::max_nbrGi = VEVPMFH_FP::max(matrprop.nbrGi, inclprop.nbrGi);
+        VEVPMFH_FP::max_nbrKj = VEVPMFH_FP::max(matrprop.nbrKj, inclprop.nbrKj);*/
+}
+  
   
   
-mlawVEVPMFH::mlawVEVPMFH(const mlawVEVPMFH &source) : materialLaw(source), _rho(source._rho)
-  {
+  
+mlawVEVPMFH::mlawVEVPMFH(const mlawVEVPMFH &source) : materialLaw(source), _rho(source._rho){
+
      // create and copy
      _rho = source._rho;
-     char* input = "propname.dat";
-     VEVPMFH_FP::readinput (matrprop, inclprop, input);
-     matrprop->E0 = (source.matrprop)->E0;
-     matrprop->nu = (source.matrprop)->nu;
+     inputfilename = source.inputfilename;
+     
+     VEVPMFH_FP::readinput (&matrprop, &inclprop, source.inputfilename);
+     matrprop.E0 = (source.matrprop).E0;
+     matrprop.nu = (source.matrprop).nu;
     
-     matrprop->yldstrs = (source.matrprop)->yldstrs;
-     matrprop->hardmodel = (source.matrprop)->hardmodel;
-     matrprop->hardmodulus1 = (source.matrprop)->hardmodulus1;
-     matrprop->hardmodulus2 = (source.matrprop)->hardmodulus2;
-     matrprop->hardexp = (source.matrprop)->hardexp;
+     matrprop.yldstrs = (source.matrprop).yldstrs;
+     matrprop.hardmodel = (source.matrprop).hardmodel;
+     matrprop.hardmodulus1 = (source.matrprop).hardmodulus1;
+     matrprop.hardmodulus2 = (source.matrprop).hardmodulus2;
+     matrprop.hardexp = (source.matrprop).hardexp;
      
-     matrprop->viscmodel = (source.matrprop)->viscmodel;
-     matrprop->viscexp1 = (source.matrprop)->viscexp1;
-     matrprop->viscexp2 = (source.matrprop)->viscexp2;
-     matrprop->visccoef = (source.matrprop)->visccoef;
+     matrprop.viscmodel = (source.matrprop).viscmodel;
+     matrprop.viscexp1 = (source.matrprop).viscexp1;
+     matrprop.viscexp2 = (source.matrprop).viscexp2;
+     matrprop.visccoef = (source.matrprop).visccoef;
      
-     matrprop->v = (source.matrprop)->v;
-     matrprop->ar = (source.matrprop)->ar;
+     matrprop.v = (source.matrprop).v;
+     matrprop.ar = (source.matrprop).ar;
      
      int i;
      for(i=0; i<3; i++){
-       matrprop->orient[i] = (source.matrprop)->orient[i];
+       matrprop.orient[i] = (source.matrprop).orient[i];
      }
      
-     matrprop->G0 = (source.matrprop)->G0;
-     matrprop->G_inf = (source.matrprop)->G_inf;
-     matrprop->nbrGi = (source.matrprop)->nbrGi;
+     matrprop.G0 = (source.matrprop).G0;
+     matrprop.G_inf = (source.matrprop).G_inf;
+     matrprop.nbrGi = (source.matrprop).nbrGi;
      
-     for(i=0; i< matrprop->nbrGi; i++){
-       matrprop->Gi[i] = (source.matrprop)->Gi[i];
-       matrprop->gi[i] = (source.matrprop)->gi[i];
+     for(i=0; i< matrprop.nbrGi; i++){
+       matrprop.Gi[i] = (source.matrprop).Gi[i];
+       matrprop.gi[i] = (source.matrprop).gi[i];
      }
      
-     matrprop->K0 = (source.matrprop)->K0;
-     matrprop->K_inf = (source.matrprop)->K_inf;
-     matrprop->nbrKj = (source.matrprop)->nbrKj;
+     matrprop.K0 = (source.matrprop).K0;
+     matrprop.K_inf = (source.matrprop).K_inf;
+     matrprop.nbrKj = (source.matrprop).nbrKj;
      
-     for(i=0; i<matrprop->nbrKj ; i++){
-       matrprop->Kj[i] = (source.matrprop)->Kj[i];
-       matrprop->kaj[i] = (source.matrprop)->kaj[i];
+     for(i=0; i<matrprop.nbrKj ; i++){
+       matrprop.Kj[i] = (source.matrprop).Kj[i];
+       matrprop.kaj[i] = (source.matrprop).kaj[i];
      }
-  }
+     
+     // Copy inclprop remaining 
+}
+
+
+
+
+
+mlawVEVPMFH& mlawVEVPMFH::operator=(const materialLaw &source){
 
-  mlawVEVPMFH& mlawVEVPMFH::operator=(const materialLaw &source)
-  {
      materialLaw::operator=(source);
      const mlawVEVPMFH* src =static_cast<const mlawVEVPMFH*>(&source);
      
      //copy
      _rho = src->_rho;
-     matrprop->E0 = (src->matrprop)->E0;
-     matrprop->nu = (src->matrprop)->nu;
+     matrprop.E0 = (src->matrprop).E0;
+     matrprop.nu = (src->matrprop).nu;
     
-     matrprop->yldstrs = (src->matrprop)->yldstrs;
-     matrprop->hardmodel = (src->matrprop)->hardmodel;
-     matrprop->hardmodulus1 = (src->matrprop)->hardmodulus1;
-     matrprop->hardmodulus2 = (src->matrprop)->hardmodulus2;
-     matrprop->hardexp = (src->matrprop)->hardexp;
+     matrprop.yldstrs = (src->matrprop).yldstrs;
+     matrprop.hardmodel = (src->matrprop).hardmodel;
+     matrprop.hardmodulus1 = (src->matrprop).hardmodulus1;
+     matrprop.hardmodulus2 = (src->matrprop).hardmodulus2;
+     matrprop.hardexp = (src->matrprop).hardexp;
      
-     matrprop->viscmodel = (src->matrprop)->viscmodel;
-     matrprop->viscexp1 = (src->matrprop)->viscexp1;
-     matrprop->viscexp2 = (src->matrprop)->viscexp2;
-     matrprop->visccoef = (src->matrprop)->visccoef;
+     matrprop.viscmodel = (src->matrprop).viscmodel;
+     matrprop.viscexp1 = (src->matrprop).viscexp1;
+     matrprop.viscexp2 = (src->matrprop).viscexp2;
+     matrprop.visccoef = (src->matrprop).visccoef;
      
-     matrprop->v = (src->matrprop)->v;
-     matrprop->ar = (src->matrprop)->ar;
+     matrprop.v = (src->matrprop).v;
+     matrprop.ar = (src->matrprop).ar;
      
      int i;
      for(i=0; i<3; i++){
-       matrprop->orient[i] = (src->matrprop)->orient[i];
+       matrprop.orient[i] = (src->matrprop).orient[i];
      }
      
-     matrprop->G0 = (src->matrprop)->G0;
-     matrprop->G_inf = (src->matrprop)->G_inf;
-     matrprop->nbrGi = (src->matrprop)->nbrGi;
+     matrprop.G0 = (src->matrprop).G0;
+     matrprop.G_inf = (src->matrprop).G_inf;
+     matrprop.nbrGi = (src->matrprop).nbrGi;
      
-     for(i=0; i< matrprop->nbrGi; i++){
-       matrprop->Gi[i] = (src->matrprop)->Gi[i];
-       matrprop->gi[i] = (src->matrprop)->gi[i];
+     for(i=0; i< matrprop.nbrGi; i++){
+       matrprop.Gi[i] = (src->matrprop).Gi[i];
+       matrprop.gi[i] = (src->matrprop).gi[i];
      }
      
-     matrprop->K0 = (src->matrprop)->K0;
-     matrprop->K_inf = (src->matrprop)->K_inf;
-     matrprop->nbrKj = (src->matrprop)->nbrKj;
+     matrprop.K0 = (src->matrprop).K0;
+     matrprop.K_inf = (src->matrprop).K_inf;
+     matrprop.nbrKj = (src->matrprop).nbrKj;
      
-     for(i=0; i<matrprop->nbrKj ; i++){
-       matrprop->Kj[i] = (src->matrprop)->Kj[i];
-       matrprop->kaj[i] = (src->matrprop)->kaj[i];
+     for(i=0; i<matrprop.nbrKj ; i++){
+       matrprop.Kj[i] = (src->matrprop).Kj[i];
+       matrprop.kaj[i] = (src->matrprop).kaj[i];
      }
-
+  
+  
+  // copy incl prop remaining 
   return *this;
 }
 
 
-  mlawVEVPMFH::~mlawVEVPMFH()
-  {
-    free(matrprop->Gi);
-    free(matrprop->gi);
-    free(matrprop->Kj);
-    free(matrprop->kaj);
-  }
+
+
+
+mlawVEVPMFH::~mlawVEVPMFH(){
+    free(matrprop.Gi);
+    free(matrprop.gi);
+    free(matrprop.Kj);
+    free(matrprop.kaj);
+}
+
+
+
   
   
-void mlawVEVPMFH::createIPState(IPStateBase* &ips,bool hasBodyForce,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
-{
+void mlawVEVPMFH::createIPState(IPStateBase* &ips,bool hasBodyForce,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const{
   IPVariable* ivi = new IPVEVPMFH();
   IPVariable* iv1 = new IPVEVPMFH();
   IPVariable* iv2 = new IPVEVPMFH();
@@ -147,68 +164,76 @@ void mlawVEVPMFH::createIPState(IPStateBase* &ips,bool hasBodyForce,const bool*
   ips = new IP3State(state_,ivi,iv1,iv2);
 }
 
-double mlawVEVPMFH::soundSpeed() const
-{
-  double nu = matrprop->nu;
-  double E = matrprop->E0;
+
+
+
+double mlawVEVPMFH::soundSpeed() const{
+  double nu = matrprop.nu;
+  double E = matrprop.E0;
   double mu = (E*nu)/((1+nu)*(1-2*nu));
   double factornu = (1.-nu)/((1.+nu)*(1.-2.*nu));
   
   return sqrt(E*factornu/_rho);
+
 }
 
 
+
 void mlawVEVPMFH::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P,const IPVariable *ipvprevi,
                                       IPVariable *ipvcuri,STensor43 &Tangent,const bool stiff,STensor43* elasticTangent, 
                                       const bool dTangent,
-                                      STensor63* dCalgdeps) const
-{
+                                      STensor63* dCalgdeps) const{
+                                      
+  
+  double **Calgo; 	VEVPMFH_FP::mallocmatrix(&Calgo,6,6);
+  double *dstrn; 	dstrn = (double*)malloc(6*sizeof(double));
+  double dt;
+  
   if(elasticTangent!=NULL) Msg::Error(" mlawVEVPMFH STensor43* elasticTangent not defined");
   const IPVEVPMFH *ipvprev = dynamic_cast<const IPVEVPMFH *> (ipvprevi);
   IPVEVPMFH *ipvcur = dynamic_cast<IPVEVPMFH *> (ipvcuri);
-  VEVPMFH_FP::STATE* statetn=ipvprev->_state;
-  VEVPMFH_FP::STATE* statetau=ipvcur->_state;
+  const VEVPMFH_FP::STATE* statetn= ipvprev->getSTATE();
+  VEVPMFH_FP::STATE* statetau= ipvcur->getSTATE();
   
+  // Set time step
+  dt = _timeStep;
+  //printf("\n time step in %g \n", dt);
   
-  double **Calgo; VEVPMFH_FP::mallocmatrix(&Calgo,6,6);
-  double *dstrn; dstrn = (double*)malloc(6*sizeof(double));
-  double dt = getTimeStep();
-  
+  // Compute logarithmic strain increment from deformation gradient F
   for(int i=0;i<3;i++){
-    dstrn[i]=Fn(i+1,i+1)-F0(i+1,i+1);
+    dstrn[i]=Fn(i,i)-F0(i,i);
+  }
+  
+  dstrn[3]=0.5*sqrt(2)*(Fn(0,1)-F0(0,1)+Fn(1,0)-F0(1,0));
+  dstrn[4]=0.5*sqrt(2)*(Fn(1,2)-F0(1,2)+Fn(2,1)-F0(2,1));
+  dstrn[5]=0.5*sqrt(2)*(Fn(0,2)-F0(0,2)+Fn(2,0)-F0(2,0));
+  
+  /*printf("\n dstrn:\n");
+  VEVPMFH_FP::printvect(dstrn,6);
+  
+  printf("\n Fn:\n");
+  for(int k=1; k<4;k++){
+    for(int n=1; n<4; n++){
+      printf("F(%i,%i)= %g\t", k, n, Fn(k,n));
+    }
+    printf("\n");
   }
   
-  dstrn[3]=0.5*sqrt(2)*(Fn(1,2)-F0(1,2)+Fn(2,1)-F0(2,1));
-  dstrn[4]=0.5*sqrt(2)*(Fn(2,3)-F0(2,3)+Fn(3,2)-F0(3,2));
-  dstrn[5]=0.5*sqrt(2)*(Fn(1,3)-F0(1,3)+Fn(3,1)-F0(3,1));
+  */
+  // Compute state at t = 0
+  if (dt == 0){
   
-  /** Call the constitutive box in VEVP small strains **/
-  /* If converged : return 1				**/
-  /* Else: return 0					**/
-  /*****************************************************/
-  int error= ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(matrprop, statetn, dt, dstrn, statetau, Calgo);
-
-  if(!error)
-  {
-    /** update viscous stresses after convergence **/
-    double *dstrnp; dstrnp = (double*)malloc(6*sizeof(double));
-    double *dstrn_ve; dstrn_ve = (double*)malloc(6*sizeof(double));
-    double *dstrn_ve_vol; dstrn_ve_vol = (double*)malloc(6*sizeof(double));
-    double *dstrn_ve_dev; dstrn_ve_dev = (double*)malloc(6*sizeof(double));
+    double *dstrs; 	dstrs = (double*)malloc(6*sizeof(double)); 
+    VEVPMFH_FP::compute_Hooke(matrprop.E0, matrprop.nu, Calgo);
+    VEVPMFH_FP::contraction42(Calgo, dstrn, dstrs);
+    VEVPMFH_FP::addtens2(1, statetn->strs, 1, dstrs, statetau->strs);
     
-    VEVPMFH_FP::addtens2(-1, statetn->pstrn, 1, statetau->pstrn, dstrnp); /*plastic strain increment */
-    VEVPMFH_FP::addtens2(-1, dstrnp, 1, dstrn, dstrn_ve); 	    /*VE strain increment*/
-
-    VEVPMFH_FP::decomposition_dstrn_ve(dstrn_ve, dstrn_ve_vol, dstrn_ve_dev);
-     
-    VEVPMFH_FP::compute_SigmaHj(matrprop, statetn, dt, dstrn_ve_vol, statetau->Sigma_Hj);
-    VEVPMFH_FP::compute_Si(matrprop, statetn, dt, dstrn_ve_dev, statetau->Si);
+    VEVPMFH_FP::addtens2(1, statetn->strn, 1, dstrn, statetau->strn);
+    VEVPMFH_FP::copystate(statetau, &(ipvcur->_state));
     
-    free(dstrnp); free(dstrn_ve); free(dstrn_ve_vol); free(dstrn_ve_dev);
       
     //convert str to P  :order xx yy zz xy xz yz
-    for(int i=0;i<3;i++)
-    {
+    for(int i=0;i<3;i++){
       P(i,i) = statetau->strs[i];
     }
     P(0,1) = statetau->strs[3]/ sqrt(2);
@@ -221,10 +246,6 @@ void mlawVEVPMFH::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P
     
     //convert MaterialTensor to Tangent  :order xx yy zz xy xz yz
    if(stiff){
-    //if(mat->isInLargeDeformation())
-    //{
-    //  Msg::Info("Material tensor does not exist in large deformation");
-    //}
     Tangent(0,0,0,0) = Calgo[0][0];
     Tangent(0,0,0,1) = Calgo[0][3]/sqrt(2);
     Tangent(0,0,0,2) = Calgo[0][5]/sqrt(2);
@@ -270,10 +291,10 @@ void mlawVEVPMFH::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P
     Tangent(1,0,0,2) = Calgo[3][5]/2;
     Tangent(1,0,1,0) = Calgo[3][3]/2;
     Tangent(1,0,1,1) = Calgo[3][1]/sqrt(2);
-    Tangent(1,0,1,2) = Calgo[3][4]/sqrt(2);
-    Tangent(1,0,2,0) = Calgo[3][5]/sqrt(2);
-    Tangent(1,0,2,1) = Calgo[3][4]/sqrt(2);
-    Tangent(1,0,2,2) = Calgo[3][2]/2;
+    Tangent(1,0,1,2) = Calgo[3][4]/2;
+    Tangent(1,0,2,0) = Calgo[3][5]/2;
+    Tangent(1,0,2,1) = Calgo[3][4]/2;
+    Tangent(1,0,2,2) = Calgo[3][2]/sqrt(2);
 
     Tangent(0,2,0,0) = Calgo[5][0]/sqrt(2);
     Tangent(0,2,0,1) = Calgo[5][3]/2;
@@ -314,41 +335,187 @@ void mlawVEVPMFH::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P
     Tangent(1,2,2,0) = Calgo[4][5]/2;
     Tangent(1,2,2,1) = Calgo[4][4]/2;
     Tangent(1,2,2,2) = Calgo[4][2]/sqrt(2);
-   }
+    
+    }
    
-    ipvcur->_elasticEne = 0.;
-    ipvcur->_plasticEne = 0.;
+  free(dstrs);
+  }
+  
+  else{
+    int error= ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(&matrprop, statetn, dt, dstrn, statetau, Calgo);
+
+    if(error){
+      /** update viscous stresses after convergence **/
+      double *dstrnp; dstrnp = (double*)malloc(6*sizeof(double));
+      double *dstrn_ve; dstrn_ve = (double*)malloc(6*sizeof(double));
+      double *dstrn_ve_vol; dstrn_ve_vol = (double*)malloc(6*sizeof(double));
+      double *dstrn_ve_dev; dstrn_ve_dev = (double*)malloc(6*sizeof(double));
+    
+      VEVPMFH_FP::addtens2(-1, statetn->pstrn, 1, statetau->pstrn, dstrnp); 
+      VEVPMFH_FP::addtens2(-1, dstrnp, 1, dstrn, dstrn_ve); 	    
+
+      VEVPMFH_FP::decomposition_dstrn_ve(dstrn_ve, dstrn_ve_vol, dstrn_ve_dev);
+     
+      VEVPMFH_FP::compute_SigmaHj(&matrprop, statetn, dt, dstrn_ve_vol, statetau->Sigma_Hj);
+      VEVPMFH_FP::compute_Si(&matrprop, statetn, dt, dstrn_ve_dev, statetau->Si);
+    
+      VEVPMFH_FP::addtens2(1, statetn->strn, 1, dstrn, statetau->strn);
+      VEVPMFH_FP::copystate(statetau, &(ipvcur->_state));
+    
+      free(dstrnp); free(dstrn_ve); free(dstrn_ve_vol); free(dstrn_ve_dev);
+      
+      //convert str to P  :order xx yy zz xy xz yz
+      for(int i=0;i<3;i++) {
+        P(i,i) = statetau->strs[i];
+      }
+      P(0,1) = statetau->strs[3]/ sqrt(2);
+      P(1,0) = P(0,1);
+      P(0,2) = statetau->strs[5]/ sqrt(2);
+      P(2,0) = P(0,2);
+      P(1,2) = statetau->strs[4]/ sqrt(2);
+      P(2,1) = P(1,2);
+    
+    
+     //convert MaterialTensor to Tangent  :order xx yy zz xy xz yz
+     if(stiff){
+      Tangent(0,0,0,0) = Calgo[0][0];
+      Tangent(0,0,0,1) = Calgo[0][3]/sqrt(2);
+      Tangent(0,0,0,2) = Calgo[0][5]/sqrt(2);
+      Tangent(0,0,1,0) = Calgo[0][3]/sqrt(2);
+      Tangent(0,0,1,1) = Calgo[0][1];
+      Tangent(0,0,1,2) = Calgo[0][4]/sqrt(2);
+      Tangent(0,0,2,0) = Calgo[0][5]/sqrt(2);
+      Tangent(0,0,2,1) = Calgo[0][4]/sqrt(2);
+      Tangent(0,0,2,2) = Calgo[0][2];
+
+      Tangent(1,1,0,0) = Calgo[1][0];  
+      Tangent(1,1,0,1) = Calgo[1][3]/sqrt(2);
+      Tangent(1,1,0,2) = Calgo[1][5]/sqrt(2);
+      Tangent(1,1,1,0) = Calgo[1][3]/sqrt(2);
+      Tangent(1,1,1,1) = Calgo[1][1];
+      Tangent(1,1,1,2) = Calgo[1][4]/sqrt(2);
+      Tangent(1,1,2,0) = Calgo[1][5]/sqrt(2);
+      Tangent(1,1,2,1) = Calgo[1][4]/sqrt(2);
+      Tangent(1,1,2,2) = Calgo[1][2];
+
+      Tangent(2,2,0,0) = Calgo[2][0];  
+      Tangent(2,2,0,1) = Calgo[2][3]/sqrt(2);
+      Tangent(2,2,0,2) = Calgo[2][5]/sqrt(2);
+      Tangent(2,2,1,0) = Calgo[2][3]/sqrt(2);
+      Tangent(2,2,1,1) = Calgo[2][1];
+      Tangent(2,2,1,2) = Calgo[2][4]/sqrt(2);
+      Tangent(2,2,2,0) = Calgo[2][5]/sqrt(2);
+      Tangent(2,2,2,1) = Calgo[2][4]/sqrt(2);
+      Tangent(2,2,2,2) = Calgo[2][2];
+
+      Tangent(0,1,0,0) = Calgo[3][0]/sqrt(2);
+      Tangent(0,1,0,1) = Calgo[3][3]/2;
+      Tangent(0,1,0,2) = Calgo[3][5]/2;
+      Tangent(0,1,1,0) = Calgo[3][3]/2;
+      Tangent(0,1,1,1) = Calgo[3][1]/sqrt(2);
+      Tangent(0,1,1,2) = Calgo[3][4]/2;
+      Tangent(0,1,2,0) = Calgo[3][5]/2;
+      Tangent(0,1,2,1) = Calgo[3][4]/2;
+      Tangent(0,1,2,2) = Calgo[3][2]/sqrt(2);
+
+      Tangent(1,0,0,0) = Calgo[3][0]/sqrt(2);
+      Tangent(1,0,0,1) = Calgo[3][3]/2;
+      Tangent(1,0,0,2) = Calgo[3][5]/2;
+      Tangent(1,0,1,0) = Calgo[3][3]/2;
+      Tangent(1,0,1,1) = Calgo[3][1]/sqrt(2);
+      Tangent(1,0,1,2) = Calgo[3][4]/2;
+      Tangent(1,0,2,0) = Calgo[3][5]/2;
+      Tangent(1,0,2,1) = Calgo[3][4]/2;
+      Tangent(1,0,2,2) = Calgo[3][2]/sqrt(2);
+
+      Tangent(0,2,0,0) = Calgo[5][0]/sqrt(2);
+      Tangent(0,2,0,1) = Calgo[5][3]/2;
+      Tangent(0,2,0,2) = Calgo[5][5]/2;
+      Tangent(0,2,1,0) = Calgo[5][3]/2;
+      Tangent(0,2,1,1) = Calgo[5][1]/sqrt(2);
+      Tangent(0,2,1,2) = Calgo[5][4]/2;
+      Tangent(0,2,2,0) = Calgo[5][5]/2;
+      Tangent(0,2,2,1) = Calgo[5][4]/2;
+      Tangent(0,2,2,2) = Calgo[5][2]/sqrt(2);
+
+      Tangent(2,0,0,0) = Calgo[5][0]/sqrt(2);
+      Tangent(2,0,0,1) = Calgo[5][3]/2;
+      Tangent(2,0,0,2) = Calgo[5][5]/2;
+      Tangent(2,0,1,0) = Calgo[5][3]/2;
+      Tangent(2,0,1,1) = Calgo[5][1]/sqrt(2);
+      Tangent(2,0,1,2) = Calgo[5][4]/2;
+      Tangent(2,0,2,0) = Calgo[5][5]/2;
+      Tangent(2,0,2,1) = Calgo[5][4]/2;
+      Tangent(2,0,2,2) = Calgo[5][2]/sqrt(2);
+
+      Tangent(2,1,0,0) = Calgo[4][0]/sqrt(2);
+      Tangent(2,1,0,1) = Calgo[4][3]/2;
+      Tangent(2,1,0,2) = Calgo[4][5]/2;
+      Tangent(2,1,1,0) = Calgo[4][3]/2;
+      Tangent(2,1,1,1) = Calgo[4][1]/sqrt(2);
+      Tangent(2,1,1,2) = Calgo[4][4]/2;
+      Tangent(2,1,2,0) = Calgo[4][5]/2;
+      Tangent(2,1,2,1) = Calgo[4][4]/2;
+      Tangent(2,1,2,2) = Calgo[4][2]/sqrt(2);
+
+      Tangent(1,2,0,0) = Calgo[4][0]/sqrt(2);
+      Tangent(1,2,0,1) = Calgo[4][3]/2;
+      Tangent(1,2,0,2) = Calgo[4][5]/2;
+      Tangent(1,2,1,0) = Calgo[4][3]/2;
+      Tangent(1,2,1,1) = Calgo[4][1]/sqrt(2);
+      Tangent(1,2,1,2) = Calgo[4][4]/2;
+      Tangent(1,2,2,0) = Calgo[4][5]/2;
+      Tangent(1,2,2,1) = Calgo[4][4]/2;
+      Tangent(1,2,2,2) = Calgo[4][2]/sqrt(2);
+     }
+   
+      ipvcur->_elasticEne = 0.;
+      ipvcur->_plasticEne = 0.;
+    
+    
+    }
     
     
-  }
-  else
-  {
     //constbox did not converge
-    P(0,0) = P(1,1) = P(2,2) = sqrt(-1.); // to exist NR and performed next step with a smaller incremenet
-  }
+    // to exist NR and performed next step with a smaller incremenet
+    else{
+      P(0,0) = P(1,1) = P(2,2) = sqrt(-1.); 
+    }
 
+  /*printf("\n Calgo:\n");
+  VEVPMFH_FP::printmatr(Calgo,6,6);*/
+  }
+  
+  /** Freeing **/
   VEVPMFH_FP::freematrix(Calgo,6);
   free(dstrn);
+  
 }
 
 
 
-const double mlawVEVPMFH::bulkModulus() const
-{
-  return matrprop->K0;
+
+
+const double mlawVEVPMFH::bulkModulus() const{
+  return matrprop.K0;
 }
-const double mlawVEVPMFH::shearModulus() const
-{
-  double nu = matrprop->nu;
-  double E = matrprop->E0;
+
+
+
+const double mlawVEVPMFH::shearModulus() const{
+  double nu = matrprop.nu;
+  double E = matrprop.E0;
   double mu = (E*nu)/((1+nu)*(1-2*nu));
   
   return mu ;
 }
 
-const double mlawVEVPMFH::poissonRatio() const
-{
-  return matrprop->nu;
+
+
+
+const double mlawVEVPMFH::poissonRatio() const{
+  return matrprop.nu;
+
 }
 
 
diff --git a/NonLinearSolver/materialLaw/mlawVEVPMFH.h b/NonLinearSolver/materialLaw/mlawVEVPMFH.h
index cd78dd248c03368e1f1dfc6bc33e5735c3dfe860..83dd5ea54bdf3d51c55a2ad752c05cbfbac0eed1 100644
--- a/NonLinearSolver/materialLaw/mlawVEVPMFH.h
+++ b/NonLinearSolver/materialLaw/mlawVEVPMFH.h
@@ -19,20 +19,19 @@
 #include "fonctions_pratiques.h"
 //#include "box_secant.h"
 
-class mlawVEVPMFH : public materialLaw
-{
-
- protected:
-  double _rho ; //maximal mu for anisotropic
-
-  // to be compatible with umat: use same variable
-
-  VEVPMFH_FP::MATPROP *matrprop, *inclprop;
+class mlawVEVPMFH : public materialLaw{
+  
+  protected:
+    double _rho; //maximal mu for anisotropic
 
- public:
-  mlawVEVPMFH(const int num, const double rho, const char *propName);
+    // to be compatible with umat: use same variable
+    VEVPMFH_FP::MATPROP matrprop;
+    VEVPMFH_FP::MATPROP inclprop;
+    const char* inputfilename;
+  public:
+    mlawVEVPMFH(const int num, const double rho, const char *propName);
 
- #ifndef SWIG
+  #ifndef SWIG
   mlawVEVPMFH(const mlawVEVPMFH &source);
   mlawVEVPMFH& operator=(const materialLaw &source);
   virtual ~mlawVEVPMFH();
@@ -45,14 +44,14 @@ class mlawVEVPMFH : public materialLaw
   virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
 
   virtual bool withEnergyDissipation() const {return true;};
-  virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; // this law is initialized so nothing to do
-  virtual double soundSpeed() const; // default but you can redefine it for your case
+  virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; 	// this law is initialized so nothing to do
+  virtual double soundSpeed() const; 						// default but you can redefine it for your case
   virtual const double bulkModulus() const;
   virtual const double shearModulus() const;
   virtual const double poissonRatio() const;
 
   // specific function
- public:
+  public:
   virtual void constitutive(
                             const STensor3& F0,         // initial deformation gradient (input @ time n)
                             const STensor3& Fn,         // updated deformation gradient (input @ time n+1)
@@ -70,7 +69,7 @@ class mlawVEVPMFH : public materialLaw
 
   /*int getNsdv() const { return nsdv;}*/
   
- #endif // SWIG
+  #endif // SWIG
 };
 
 #endif // MLAWNONLOCALDAMAGE_H_
diff --git a/dG3D/benchmarks/uniaxialTest_PBC/run.py b/dG3D/benchmarks/uniaxialTest_PBC/run.py
index 6fd6851b7e12cd3c27d59288fa793058d082e16c..729ecbabc89dcbc3f6bbb184eaabf2d788f5c156 100644
--- a/dG3D/benchmarks/uniaxialTest_PBC/run.py
+++ b/dG3D/benchmarks/uniaxialTest_PBC/run.py
@@ -3,6 +3,7 @@
 from gmshpy import *
 from dG3Dpy import*
 
+
 #script to launch PBC problem with a python script
 
 # material law
@@ -16,6 +17,8 @@ rho = 2.7e-9 # Bulk mass
 
 law1 = dG3DLinearElasticMaterialLaw(11,rho,E,nu)
 
+
+
 # geometry
 meshfile="rve.msh" # name of mesh file
 
diff --git a/dG3D/benchmarks/voidedRVE/voidedRVE.py b/dG3D/benchmarks/voidedRVE/voidedRVE.py
index 439169c0c073940eba76b81ad697c84b043a40b1..c4db68d2b4c97cf95206d608ea076c52ade6d8a6 100644
--- a/dG3D/benchmarks/voidedRVE/voidedRVE.py
+++ b/dG3D/benchmarks/voidedRVE/voidedRVE.py
@@ -1,7 +1,7 @@
 #coding-Utf-8-*-
 from gmshpy import *
-#from dG3DpyDebug import*
-from dG3Dpy import*
+from dG3DpyDebug import*
+#from dG3Dpy import*
 
 #script to launch composite problem with a python script
 
@@ -38,6 +38,8 @@ number=1
 # creation of law
 harden = PowerLawJ2IsotropicHardening(1,sy0EP, hEP, hexp,0.01)
 law1 = J2LinearDG3DMaterialLaw(lawep,rhoEP,youngEP,nuEP,harden)
+#BulkLawNum1 = 1
+#law1 = VEVPMFHDG3DMaterialLaw(BulkLawNum1, rho, "propname.dat")
 
 # creation of ElasticField
 Mfield = 12 # number of the field (physical number of Mtrix)
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 049902454ee960cbf8b5cd55a45112ce832b4296..1761e29feb045ec460ae5d7d7bfebc6cd8a82188 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -2178,6 +2178,87 @@ void localDamageIsotropicElasticityDG3DIPVariable::restart()
   return;
 }
 
+
+//====================================================VEVPMFHDG3DIPVariable==========================================================//
+//
+VEVPMFHDG3DIPVariable::VEVPMFHDG3DIPVariable(bool createBodyForceHO, const bool oninter) : dG3DIPVariable(createBodyForceHO, oninter)
+{
+  _VEVPipv = new IPVEVPMFH();
+}
+
+//VEVPMFHDG3DIPVariable::VEVPMFHDG3DIPVariable(const VEVPMFHDG3DIPVariable &source) : dG3DIPVariable(source), _tiipv(source._tiipv)
+VEVPMFHDG3DIPVariable::VEVPMFHDG3DIPVariable(const VEVPMFHDG3DIPVariable &source) : dG3DIPVariable(source)
+{
+  _VEVPipv = NULL;
+  if (source._VEVPipv){
+    _VEVPipv = new IPVEVPMFH();
+
+  }
+}
+
+
+VEVPMFHDG3DIPVariable& VEVPMFHDG3DIPVariable::operator=(const IPVariable &source)
+{
+  dG3DIPVariable::operator=(source);
+  const VEVPMFHDG3DIPVariable* src = dynamic_cast<const VEVPMFHDG3DIPVariable*>(&source);
+  /*if(src != NULL)
+    _tiipv.operator=(*dynamic_cast<const IPVariable*>(&( src->_tiipv)));
+  return *this;*/
+  
+  if(src != NULL){
+    if (src->_VEVPipv){
+      if (_VEVPipv){
+        _VEVPipv->operator=(*dynamic_cast<const IPVariable*>(src->_VEVPipv));
+      }
+      else{
+        _VEVPipv = dynamic_cast<IPVEVPMFH*>(src->_VEVPipv->clone());
+      }
+
+    }
+    else{
+      if (_VEVPipv) {
+        delete _VEVPipv; _VEVPipv = NULL;
+      }
+    }
+  }
+  return *this;
+}
+
+double VEVPMFHDG3DIPVariable::get(const int comp) const
+{
+  double val=_VEVPipv->get(comp);
+  if(val==0.)
+  {
+    return dG3DIPVariable::get(comp);
+  }
+  return val;
+}
+
+
+double VEVPMFHDG3DIPVariable::defoEnergy() const
+{
+  return getIPVEVPMFH()->defoEnergy();
+}
+
+double  VEVPMFHDG3DIPVariable::plasticEnergy() const
+{
+  return getIPVEVPMFH()->plasticEnergy();
+}
+
+void VEVPMFHDG3DIPVariable::restart()
+{
+  dG3DIPVariable::restart();
+  restartManager::restart(_VEVPipv);
+  return;
+}
+
+
+//
+
+
+
+
+
 //====================================================ViscoelasticDG3DIPVariable==========================================================//
 
 ViscoelasticDG3DIPVariable::ViscoelasticDG3DIPVariable(const mlawViscoelastic &_Vislaw, const bool createBodyForceHO, const bool oninter) : dG3DIPVariable(createBodyForceHO, oninter)
@@ -2670,51 +2751,10 @@ void AnisotropicStochDG3DIPVariable::restart()
   return;
 }
 
-//
-VEVPMFHDG3DIPVariable::VEVPMFHDG3DIPVariable(bool createBodyForceHO, const bool oninter) : dG3DIPVariable(createBodyForceHO, oninter)
-{
 
-}
-VEVPMFHDG3DIPVariable::VEVPMFHDG3DIPVariable(const VEVPMFHDG3DIPVariable &source) : dG3DIPVariable(source),
-_tiipv(source._tiipv)
-{
 
-}
-VEVPMFHDG3DIPVariable& VEVPMFHDG3DIPVariable::operator=(const IPVariable &source)
-{
-  dG3DIPVariable::operator=(source);
-  const VEVPMFHDG3DIPVariable* src = dynamic_cast<const VEVPMFHDG3DIPVariable*>(&source);
-  if(src != NULL)
-    _tiipv.operator=(*dynamic_cast<const IPVariable*>(&( src->_tiipv)));
-  return *this;
-}
-double VEVPMFHDG3DIPVariable::get(const int comp) const
-{
-  double val=_tiipv.get(comp);
-  if(val==0.)
-  {
-    return dG3DIPVariable::get(comp);
-  }
-  return val;
-}
-double VEVPMFHDG3DIPVariable::defoEnergy() const
-{
-  return getIPVEVPMFH()->defoEnergy();
-}
-double  VEVPMFHDG3DIPVariable::plasticEnergy() const
-{
-  return getIPVEVPMFH()->plasticEnergy();
-}
-
-void VEVPMFHDG3DIPVariable::restart()
-{
-  dG3DIPVariable::restart();
-  restartManager::restart(&_tiipv);
-  return;
-}
 
 
-//
 
 TFADG3DIPVariable::TFADG3DIPVariable(IPTFA* ipc, const bool createBodyForceHO, const bool oninter) : dG3DIPVariable(createBodyForceHO, oninter)
 {
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index f3bb61492b62e5bd1e056b18e3f5f2362a4d7087..f483e465fd2d2695c0428383ea7c96221ae9a5bf 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -2581,6 +2581,40 @@ class localDamageIsotropicElasticityDG3DIPVariable : public dG3DIPVariable
 
 
 
+
+
+class VEVPMFHDG3DIPVariable : public dG3DIPVariable // or store data in a different way
+{
+ protected:
+  IPVEVPMFH* _VEVPipv;
+ public:
+  VEVPMFHDG3DIPVariable(const bool createBodyForceHO, const bool oninter);
+  VEVPMFHDG3DIPVariable(const VEVPMFHDG3DIPVariable &source);
+  virtual VEVPMFHDG3DIPVariable& operator=(const IPVariable &source);
+  virtual ~VEVPMFHDG3DIPVariable(){if(_VEVPipv!=NULL) { delete _VEVPipv; _VEVPipv=NULL;}};
+
+  virtual void setLocation(const IPVariable::LOCATION loc) {
+    dG3DIPVariable::setLocation(loc);
+    if (_VEVPipv)
+       _VEVPipv->setLocation(loc);
+  };
+
+  virtual IPVariable* getInternalData() {return _VEVPipv;};
+  virtual const IPVariable* getInternalData()  const {return _VEVPipv;};
+
+ /* specific function */
+  IPVEVPMFH* getIPVEVPMFH(){return _VEVPipv;}
+  const IPVEVPMFH* getIPVEVPMFH() const{return _VEVPipv;}
+  virtual double get(const int i) const;
+  virtual double defoEnergy() const;
+  virtual double plasticEnergy() const;
+  virtual IPVariable* clone() const {return new VEVPMFHDG3DIPVariable(*this);};
+  virtual void restart();
+};
+
+
+
+
 class ViscoelasticDG3DIPVariable : public dG3DIPVariable
 {
  protected:
@@ -2866,32 +2900,6 @@ class AnisotropicStochDG3DIPVariable : public dG3DIPVariable // or store data in
 };
 
 
-class VEVPMFHDG3DIPVariable : public dG3DIPVariable // or store data in a different way
-{
- protected:
-  IPVEVPMFH _tiipv;
- public:
-  VEVPMFHDG3DIPVariable(const bool createBodyForceHO, const bool oninter);
-  VEVPMFHDG3DIPVariable(const VEVPMFHDG3DIPVariable &source);
-  virtual VEVPMFHDG3DIPVariable& operator=(const IPVariable &source);
-
-  virtual void setLocation(const IPVariable::LOCATION loc) {
-    dG3DIPVariable::setLocation(loc);
-    _tiipv.setLocation(loc);
-  };
-
-  virtual IPVariable* getInternalData() {return &_tiipv;};
-  virtual const IPVariable* getInternalData()  const {return &_tiipv;};
-
- /* specific function */
-  IPVEVPMFH* getIPVEVPMFH(){return &_tiipv;}
-  const IPVEVPMFH* getIPVEVPMFH() const{return &_tiipv;}
-  virtual double get(const int i) const;
-  virtual double defoEnergy() const;
-  virtual double plasticEnergy() const;
-  virtual IPVariable* clone() const {return new VEVPMFHDG3DIPVariable(*this);};
-  virtual void restart();
-};
 
 
 class TFADG3DIPVariable : public dG3DIPVariable // or store data in a different way
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index eeec69433dfe1269e93966d2a536b9b523babe13..fcb797284d9ca244997433be8646a363d43bf259 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -3098,6 +3098,155 @@ void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::stress(IPVariable* ipv, c
 }
 
 
+
+
+//=========================================================VEVPMFHDG3DMaterialLaw===============================================================//
+
+VEVPMFHDG3DMaterialLaw::VEVPMFHDG3DMaterialLaw(const int num, const double rho, const char *propName) :
+                               dG3DMaterialLaw(num,rho,true)
+{
+  _glaw = new mlawVEVPMFH(num,rho,propName);
+  /*double nu = _glaw->poissonRatio();
+  double mu = _glaw->shearModulus();
+  double E = 2.*mu*(1.+nu);
+  _rho = _glaw->density();
+  fillElasticStiffness(E, nu, elasticStiffness);*/
+
+}
+
+
+VEVPMFHDG3DMaterialLaw::~VEVPMFHDG3DMaterialLaw()
+{
+  if (_glaw !=NULL) delete _glaw;
+  _glaw = NULL;
+};
+
+
+VEVPMFHDG3DMaterialLaw::VEVPMFHDG3DMaterialLaw(const VEVPMFHDG3DMaterialLaw &source) :
+                                                             dG3DMaterialLaw(source)
+{
+	_glaw = NULL;
+	if (source._glaw != NULL){
+		_glaw = dynamic_cast<mlawVEVPMFH*>(source._glaw->clone());
+	}
+
+}
+
+void VEVPMFHDG3DMaterialLaw::setTime(const double t,const double dtime){
+  dG3DMaterialLaw::setTime(t,dtime);
+  _glaw->setTime(t,dtime);
+  return;
+}
+
+
+void VEVPMFHDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+  // check interface element
+  bool inter=true;
+  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+  if(iele==NULL) inter=false;
+
+
+  IPVariable* ipvi = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
+  IPVariable* ipv1 = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
+  IPVariable* ipv2 = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
+
+  if(ips != NULL) delete ips;
+  ips = new IP3State(state_,ipvi,ipv1,ipv2);
+}
+
+void VEVPMFHDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+  if(ipv !=NULL) delete ipv;
+  bool inter=true;
+  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+  if(iele == NULL){
+    inter=false;
+  }
+  
+  ipv = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
+
+}
+
+
+
+double VEVPMFHDG3DMaterialLaw::soundSpeed() const{return _glaw->soundSpeed();}
+
+void VEVPMFHDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const{
+  /* get ipvariable */
+  VEVPMFHDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
+  const VEVPMFHDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
+
+  FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+  if(ipvtmp !=NULL)
+  {
+    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipvtmp->getIPvBulk());
+    const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+   }
+  else
+  {
+    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipv);
+    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvp);
+  }
+  _glaw->checkInternalState(ipvcur->getInternalData(),ipvprev->getInternalData());
+
+};
+
+
+
+void VEVPMFHDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent)
+{
+  /* get ipvariable */
+  VEVPMFHDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
+  const VEVPMFHDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
+
+  FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+  if(ipvtmp !=NULL)
+  {
+    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipvtmp->getIPvBulk());
+    const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+   }
+   
+  else
+  {
+    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipv);
+    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvp);
+  }
+
+  /* strain xyz */
+  const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
+  const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
+
+  /* data for J2 law */
+  IPVEVPMFH* q1 = ipvcur->getIPVEVPMFH();
+  const IPVEVPMFH* q0 = ipvprev->getIPVEVPMFH();
+
+  /* compute stress */
+  _glaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
+                        stiff, NULL,dTangent,ipvcur->getPtrTodCmdFm());
+
+  ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+
+}
+
+
+double VEVPMFHDG3DMaterialLaw::scaleFactor() const{return _glaw->scaleFactor();};
+
+void VEVPMFHDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
+	dG3DMaterialLaw::setMacroSolver(sv);
+	if (_glaw != NULL){
+		_glaw->setMacroSolver(sv);
+	}
+};
+
+//
+
+
+
+
+
 //=========================================================ViscoelastiDG3DMaterialLaw===============================================================//
 std::vector<double> ViscoelasticDG3DMaterialLaw::read_file(const char* file_name, int file_size)
 {
@@ -6439,144 +6588,7 @@ void gursonUMATDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
 };
 
 
-//
-VEVPMFHDG3DMaterialLaw::VEVPMFHDG3DMaterialLaw(const int num, const double rho, const char *propName) :
-                               dG3DMaterialLaw(num,rho,true)
-{
-  _glaw = new mlawVEVPMFH(num,rho,propName);
-  double nu = _glaw->poissonRatio();
-  double mu = _glaw->shearModulus();
-  double E = 2.*mu*(1.+nu);
-  _rho = _glaw->density();
-  fillElasticStiffness(E, nu, elasticStiffness);
-
-
-}
-
-VEVPMFHDG3DMaterialLaw::~VEVPMFHDG3DMaterialLaw()
-{
-  if (_glaw !=NULL) delete _glaw;
-  _glaw = NULL;
-};
-
-
-VEVPMFHDG3DMaterialLaw::VEVPMFHDG3DMaterialLaw(const VEVPMFHDG3DMaterialLaw &source) :
-                                                             dG3DMaterialLaw(source)
-{
-	_glaw = NULL;
-	if (source._glaw != NULL){
-		_glaw = dynamic_cast<mlawVEVPMFH*>(source._glaw->clone());
-	}
-
-}
-
-void VEVPMFHDG3DMaterialLaw::setTime(const double t,const double dtime){
-  dG3DMaterialLaw::setTime(t,dtime);
-  _glaw->setTime(t,dtime);
-  return;
-}
-
-materialLaw::matname VEVPMFHDG3DMaterialLaw::getType() const {return _glaw->getType();}
-
-bool VEVPMFHDG3DMaterialLaw::withEnergyDissipation() const {return _glaw->withEnergyDissipation();};
-
-void VEVPMFHDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
-{
-  // check interface element
-  bool inter=true;
-  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
-  if(iele==NULL) inter=false;
-
-
-  IPVariable* ipvi = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
-  IPVariable* ipv1 = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
-  IPVariable* ipv2 = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
-
-  if(ips != NULL) delete ips;
-  ips = new IP3State(state_,ipvi,ipv1,ipv2);
-}
-
-void VEVPMFHDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
-{
-  if(ipv !=NULL) delete ipv;
-  bool inter=true;
-  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
-
-  ipv = new  VEVPMFHDG3DIPVariable(hasBodyForce,inter);
-
-}
-
-
-double VEVPMFHDG3DMaterialLaw::soundSpeed() const{return _glaw->soundSpeed();}
-
-void VEVPMFHDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const{
-  /* get ipvariable */
-  VEVPMFHDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
-  const VEVPMFHDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
-
-  FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
-  if(ipvtmp !=NULL)
-  {
-
-    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipvtmp->getIPvBulk());
-    const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
-    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvtmp2->getIPvBulk());
-   }
-  else
-  {
-    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipv);
-    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvp);
-  }
-  _glaw->checkInternalState(ipvcur->getInternalData(),ipvprev->getInternalData());
-
-};
-
-void VEVPMFHDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent)
-{
-  /* get ipvariable */
-  VEVPMFHDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
-  const VEVPMFHDG3DIPVariable* ipvprev; //= static_cast <const  nonLocalDamageDG3DIPVariable *> (ipvp);
-
-  FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
-  if(ipvtmp !=NULL)
-  {
-
-    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipvtmp->getIPvBulk());
-    const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
-    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvtmp2->getIPvBulk());
-   }
-  else
-  {
-    ipvcur = static_cast<VEVPMFHDG3DIPVariable*>(ipv);
-    ipvprev = static_cast<const VEVPMFHDG3DIPVariable*>(ipvp);
-  }
-
-  /* strain xyz */
-  const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
-  const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
-
-  /* data for J2 law */
-  IPVEVPMFH* q1 = ipvcur->getIPVEVPMFH();
-  const IPVEVPMFH* q0 = ipvprev->getIPVEVPMFH();
 
-  /* compute stress */
-  _glaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
-                        stiff, NULL,dTangent,ipvcur->getPtrTodCmdFm());
-
-  ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
-
-}
-
-double VEVPMFHDG3DMaterialLaw::scaleFactor() const{return _glaw->scaleFactor();};
-
-void VEVPMFHDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
-	dG3DMaterialLaw::setMacroSolver(sv);
-	if (_glaw != NULL){
-		_glaw->setMacroSolver(sv);
-	}
-};
-
-//
 LinearElecTherMechDG3DMaterialLaw::LinearElecTherMechDG3DMaterialLaw(const int num,const double rho,const double Ex, const double Ey, const double Ez, const double Vxy,
 			 const double Vxz, const double Vyz,const double MUxy, const double MUxz, const double MUyz, const double alpha, const double beta, const double gamma,
 			 const double t0,const double Kx,const double Ky,const double Kz,const double alphax,const double alphay,const double alphaz,const  double lx,const  double ly,
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index f0100f9143f1786fcb3cb7150dba33ecdac5a820..e53d1c3e6bf28887818190ec2e4a4282486f651c 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -829,6 +829,47 @@ class HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw  : public dG3DMaterialLaw
     #endif //SWIG
 };
 
+
+//======================================================VEVPDG3DMaterialLaw===============================================================//
+class VEVPMFHDG3DMaterialLaw : public dG3DMaterialLaw{
+  protected:
+    #ifndef SWIG
+    mlawVEVPMFH *_glaw;
+    #endif //SWIG
+
+  public:
+    VEVPMFHDG3DMaterialLaw(const int num, const double rho, const char *propName);
+    #ifndef SWIG
+    VEVPMFHDG3DMaterialLaw(const VEVPMFHDG3DMaterialLaw &source);
+    virtual ~VEVPMFHDG3DMaterialLaw();
+
+    // set the time of _cplaw
+    virtual void setTime(const double t,const double dtime);
+    virtual materialLaw::matname getType() const {return _glaw->getType();}
+    virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
+    virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
+    virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}
+    virtual double soundSpeed() const;// or change value ??
+    virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const;
+    virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
+    virtual double scaleFactor() const;
+    virtual materialLaw* clone() const{return new VEVPMFHDG3DMaterialLaw(*this);};
+    virtual bool withEnergyDissipation() const {return _glaw->withEnergyDissipation();};
+    virtual void setMacroSolver(const nonLinearMechSolver* sv);
+    virtual const materialLaw *getConstNonLinearSolverMaterialLaw() const
+    {
+      return _glaw;
+    }
+    virtual materialLaw *getNonLinearSolverMaterialLaw()
+    {
+      return _glaw;
+    }
+    #endif
+};
+
+
+
+
 //======================================================ViscoelasticDG3DMaterialLaw===============================================================//
 class ViscoelasticDG3DMaterialLaw : public dG3DMaterialLaw
 {
@@ -863,10 +904,10 @@ class ViscoelasticDG3DMaterialLaw : public dG3DMaterialLaw
   virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
   virtual const double bulkModulus() const {return _Vislaw.bulkModulus();}
   virtual std::vector<double> read_file(const char* file_name, int file_size); //read from files
-	virtual materialLaw* clone() const{return new ViscoelasticDG3DMaterialLaw(*this);};
-	virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const;
+  virtual materialLaw* clone() const{return new ViscoelasticDG3DMaterialLaw(*this);};
+  virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const;
   virtual bool withEnergyDissipation() const {return _Vislaw.withEnergyDissipation();};
-	virtual void setMacroSolver(const nonLinearMechSolver* sv){
+  virtual void setMacroSolver(const nonLinearMechSolver* sv){
 		dG3DMaterialLaw::setMacroSolver(sv);
 		_Vislaw.setMacroSolver(sv);
 	};
@@ -2265,41 +2306,7 @@ class gursonUMATDG3DMaterialLaw : public dG3DMaterialLaw{
     #endif
 };
 
-class VEVPMFHDG3DMaterialLaw : public dG3DMaterialLaw{
-  protected:
-    #ifndef SWIG
-    mlawVEVPMFH *_glaw;
-    #endif //SWIG
 
-  public:
-    VEVPMFHDG3DMaterialLaw(const int num, const double rho, const char *propName);
-    #ifndef SWIG
-    VEVPMFHDG3DMaterialLaw(const VEVPMFHDG3DMaterialLaw &source);
-    virtual ~VEVPMFHDG3DMaterialLaw();
-
-    // set the time of _cplaw
-    virtual void setTime(const double t,const double dtime);
-    virtual materialLaw::matname getType() const;
-    virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
-    virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
-    virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}
-    virtual double soundSpeed() const;// or change value ??
-    virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const;
-    virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
-    virtual double scaleFactor() const;
-    virtual materialLaw* clone() const{return new VEVPMFHDG3DMaterialLaw(*this);};
-    virtual bool withEnergyDissipation() const;
-    virtual void setMacroSolver(const nonLinearMechSolver* sv);
-    virtual const materialLaw *getConstNonLinearSolverMaterialLaw() const
-    {
-      return _glaw;
-    }
-    virtual materialLaw *getNonLinearSolverMaterialLaw()
-    {
-      return _glaw;
-    }
-    #endif
-};