diff --git a/NonLinearSolver/internalPoints/ipVEVPMFH.cpp b/NonLinearSolver/internalPoints/ipVEVPMFH.cpp
index 8887d1596c3e3f9d563e231cea8eac7362adab8f..f73dffe7d1085b771a2d6a91273f78aa1647a1ca 100644
--- a/NonLinearSolver/internalPoints/ipVEVPMFH.cpp
+++ b/NonLinearSolver/internalPoints/ipVEVPMFH.cpp
@@ -1,34 +1,28 @@
-//
-// Description: storing class for non local damage
-//
-//
-// Author:  <L. Noels>, (C) 2011
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
+/*********************************************************************************************/
+/** Copyright (C) 2022 - See COPYING file that comes with this distribution	                **/
+/** Author: Mohamed HADDAD (MEMA-iMMC-UCLouvain)							                              **/
+/** Contact: mohamed.haddad@uclouvain.be					                                      		**/
+/**											                                                                  	**/
+/** Description: storing class of IP variabales for VE-VP material law (small strains)		  **/
+/** VE-VP material law : See B.Miled et I.Doghri 2011                                       **/
+/**********************************************************************************************/
+
 #include "ipVEVPMFH.h"
 #include "ipField.h"
 #include "restartManager.h"
-#include "fonctions_pratiques.h"
-
 
 
 double IPVEVPMFH::defoEnergy() const{
   return _elasticEne;
 }
 
-
-
 double IPVEVPMFH::plasticEnergy() const{
   return _plasticEne;
 }
 
-
-
 double IPVEVPMFH::get(int comp) const{
-  
-  if(comp == IPField::SIG_XX)			// Cauchy stress field 
+  // Cauchy stress 
+  if(comp == IPField::SIG_XX)			
 	  return _state.strs[0];
   else if(comp == IPField::SIG_YY)
 	  return _state.strs[1];
@@ -41,14 +35,18 @@ double IPVEVPMFH::get(int comp) const{
   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);
+  // Von Mises equivelent stress
+  else if(comp == IPField::SVM){	
+     double svm= VEVPMFHSPACE::vonmises(_state.strs);
      return svm;
   }
-
-  else if(comp == IPField::PLASTICSTRAIN) 	// Accumulated plasticity 
+  
+  // Plastic true strain  
+  else if(comp == IPField::PLASTICSTRAIN) 	
      return _state.p;
-  else if(comp == IPField::STRAIN_XX) 	// True strain field
+
+  // Total true strain    
+  else if(comp == IPField::STRAIN_XX) 	
 	  return _state.strn[0];
   else if(comp == IPField::STRAIN_YY) 
 	  return _state.strn[1];
@@ -58,15 +56,12 @@ double IPVEVPMFH::get(int comp) const{
 	  return _state.strn[3]/sqrt(2);
   else if(comp == IPField::STRAIN_XZ)
 	  return _state.strn[5]/sqrt(2);
-  else if(comp == IPField::STRAIN_XZ)
+  else if(comp == IPField::STRAIN_YZ)
 	  return _state.strn[4]/sqrt(2);
   else
      return 0.;
 }
 
-
-
-
 void IPVEVPMFH::restart(){
   
   IPVariableMechanics::restart();
@@ -86,4 +81,3 @@ void IPVEVPMFH::restart(){
   return;
 }
 
-
diff --git a/NonLinearSolver/internalPoints/ipVEVPMFH.h b/NonLinearSolver/internalPoints/ipVEVPMFH.h
index 5c615a4642bacc31ca966ea67e9dc8e8c0860dfa..63113fabefa85a97514987c7f530a1672f62a78c 100644
--- a/NonLinearSolver/internalPoints/ipVEVPMFH.h
+++ b/NonLinearSolver/internalPoints/ipVEVPMFH.h
@@ -1,92 +1,78 @@
-//
-// Description: storing class for non local damage law it is a "virtual implementation" you have to define an ipvariable
-//              in your project which derive from this one. (all data in this ipvarible have name beggining by _nld...)
-//              so don't do the same in your project...
-// Author:  <L. NOELS>, (C) 2011
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
+/*********************************************************************************************/
+/** Copyright (C) 2022 - See COPYING file that comes with this distribution	                **/
+/** Author: Mohamed HADDAD (MEMA-iMMC-UCLouvain)							                              **/
+/** Contact: mohamed.haddad@uclouvain.be					                                      		**/
+/**											                                                                  	**/
+/** Description: storing class of IP variabales for VE-VP material law (small strains)		  **/
+/** VE-VP material law : See B.Miled et I.Doghri 2011                                       **/
+/**********************************************************************************************/
 
 #ifndef IPVEVPMFH_H_
 #define IPVEVPMFH_H_
 #include "ipvariable.h"
 #include "STensor3.h"
 #include "fullMatrix.h"
-#include "fonctions_pratiques.h"
+
+#include "UsefulFunctionsVEVPMFH.h"
 
 class IPVEVPMFH : public IPVariableMechanics{
 
-  public: // All data public to avoid the creation of function to access
-  // state variables
-    VEVPMFH_FP::STATE _state;
+  public: 
+    VEVPMFHSPACE::STATE _state;               // Structure containing all state variables (defined in UsefulFunctionsVEVPMFH.h)
     double _elasticEne;
     double _plasticEne;
   
-
   protected:
-  // parthFollowing
     double _irreversibleEnergy;
-    STensor3 _DirreversibleEnergyDF;    // Derivatives in terms of deformation gradient
+    STensor3 _DirreversibleEnergyDF;    
+
 
   public:
     IPVEVPMFH() : IPVariableMechanics(), _elasticEne(0.),_plasticEne(0.), _irreversibleEnergy(0),  _DirreversibleEnergyDF(0.){ 
-      //create and initialize
-       VEVPMFH_FP::mallocstate(&_state); 
-    }
-
+       VEVPMFHSPACE::mallocstate(&_state);    //create and initialize state structure
+    };
 
     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;
+        VEVPMFHSPACE::mallocstate(&_state);   //create and state structure
+        VEVPMFHSPACE::copystate(&(source._state), &_state); //copy state structure
+	      _elasticEne=source._elasticEne;
+ 	      _plasticEne=source._plasticEne;
         _irreversibleEnergy = source._irreversibleEnergy;
         _DirreversibleEnergyDF = source._DirreversibleEnergyDF;
-    }
-
+    };
 
     IPVEVPMFH &operator = (const IPVariable &_source){
       IPVariableMechanics::operator=(_source);
       const IPVEVPMFH* source=static_cast<const IPVEVPMFH *> (&_source);
       
-      //copy 
-      VEVPMFH_FP::copystate(&(source->_state), &_state);
+      VEVPMFHSPACE::copystate(&(source->_state), &_state);
       _elasticEne=source->_elasticEne;
       _plasticEne=source->_plasticEne;
       _irreversibleEnergy = source->_irreversibleEnergy;
       _DirreversibleEnergyDF = source->_DirreversibleEnergyDF; 
       
       return *this; 
-    }
-  
-  
+    };
   
     virtual ~IPVEVPMFH(){
-      //Destructor 
-      VEVPMFH_FP::freestate(_state);
-    }
-  
+      VEVPMFHSPACE::freestate(_state);     //Destructor of state structure
+    };
   
     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();
 
-    // 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 VEVPMFHSPACE::STATE* getSTATE() {return &_state;}
+    virtual const VEVPMFHSPACE::STATE* getSTATE() const {return &_state;}
+   
 };
 
-#endif // IPVEVPMFH_H_
+#endif // IPVEVPMFH_H_
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/CMakeLists.txt b/NonLinearSolver/materialLaw/CMakeLists.txt
index 109a929bb9d2dd4eed2eb0dcefa06b8d0fe87e0b..341ace0e973297e1d0654f9ac325875f6b289b94 100644
--- a/NonLinearSolver/materialLaw/CMakeLists.txt
+++ b/NonLinearSolver/materialLaw/CMakeLists.txt
@@ -31,6 +31,7 @@ set(SRC
   mlawTransverseIsoCurvature.cpp
   mlawTransverseIsoYarnB.cpp
   mlawNonLocalDamage.cpp
+  UsefulFunctionsVEVPMFH.cpp
   mlawVEVPMFH.cpp
   mlawNonLocalDamage_Stoch.cpp
   mlawVUMATinterface.cpp
@@ -64,8 +65,6 @@ set(SRC
   nonLocalDamageLaw.cpp
   mlawJ2SmallStrains.cpp
   numericalMaterial.cpp
-  fonctions_pratiques.cpp
-  box_secant.cpp
   mlawJ2ThermoMechanics.cpp
   mlawJ2FullyCoupledThermoMechanics.cpp
   mlawNonLocalDamageJ2FullyCoupledThermoMechanics.cpp
diff --git a/NonLinearSolver/materialLaw/fonctions_pratiques.cpp b/NonLinearSolver/materialLaw/UsefulFunctionsVEVPMFH.cpp
similarity index 73%
rename from NonLinearSolver/materialLaw/fonctions_pratiques.cpp
rename to NonLinearSolver/materialLaw/UsefulFunctionsVEVPMFH.cpp
index d99a2d05d56fd98bf76418f32abeb125dd468b92..4f11fe5405e54343e7f69a4d1d979db862964dbb 100644
--- a/NonLinearSolver/materialLaw/fonctions_pratiques.cpp
+++ b/NonLinearSolver/materialLaw/UsefulFunctionsVEVPMFH.cpp
@@ -1,18 +1,18 @@
-#ifndef FONCTIONS_PRATIQUES_C
-#define FONCTIONS_PRATIQUES_C 1
+/*************************************************************************************************/
+/** Copyright (C) 2022 - See COPYING file that comes with this distribution	                    **/
+/** Author: Mohamed HADDAD (MEMA-iMMC-UCLouvain)							                    **/
+/** Contact: mohamed.haddad@uclouvain.be					                                    **/
+/**											                                                    **/
+/** Description: storing functions needed for VE-VP material law (small strains)    		    **/
+/** VE-VP material law : See B.Miled et I.Doghri 2011                                           **/
+/*************************************************************************************************/
 
 #include <math.h>
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
 
-
-#include "fonctions_pratiques.h"
-#include "box_secant.h"
-
-
-
-
+#include "UsefulFunctionsVEVPMFH.h"
 
 
 //Description: Convert a matrix(3*3) to a vector(6*1) (Issam Doghri's convention)
@@ -20,7 +20,7 @@
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::convertmatrvect (double **matr, double *vect){
+void VEVPMFHSPACE::convertmatrvect (double **matr, double *vect){
   
   /** Declaration **/
   int i;
@@ -35,13 +35,12 @@ void VEVPMFH_FP::convertmatrvect (double **matr, double *vect){
 }
 
 
-
 //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){
+void VEVPMFHSPACE::convertvectmatr (double *vect, double **matr){
   
   /** Declaration **/
   int i;
@@ -59,14 +58,12 @@ void VEVPMFH_FP::convertvectmatr (double *vect, double **matr){
 }
 
 
-
-
 //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){
+void VEVPMFHSPACE::addtens4 (double a, double **A, double b, double **B, double **res){
   
   /** Declaration **/
   int i,j;
@@ -98,7 +95,7 @@ void VEVPMFH_FP::addtens4 (double a, double **A, double b, double **B, double **
 //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){
+void VEVPMFHSPACE::addtens2 (double a, double *A, double b, double *B, double *res){
   
   /** Declaration **/
   int i;
@@ -125,7 +122,7 @@ void VEVPMFH_FP::addtens2 (double a, double *A, double b, double *B, double *res
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::contraction22 (const double *TA, const double *TB){
+double VEVPMFHSPACE::contraction22 (const double *TA, const double *TB){
 
   /** Declaration **/
   int i;
@@ -145,7 +142,7 @@ double VEVPMFH_FP::contraction22 (const double *TA, const double *TB){
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::contraction42 (double **T4, double *T2, double *Sol){
+void VEVPMFHSPACE::contraction42 (double **T4, double *T2, double *Sol){
 
   /** Declaration **/
   int i,j;
@@ -168,7 +165,7 @@ void VEVPMFH_FP::contraction42 (double **T4, double *T2, double *Sol){
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::contraction44 (double **TA, double **TB, double **Tres){
+void VEVPMFHSPACE::contraction44 (double **TA, double **TB, double **Tres){
   
   /** Declaration **/
   int i,j,k;
@@ -194,7 +191,7 @@ void VEVPMFH_FP::contraction44 (double **TA, double **TB, double **Tres){
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::dcontr44 (double **TA, double **TB){
+double VEVPMFHSPACE::dcontr44 (double **TA, double **TB){
   
   /** Declaration **/
   int i,j;
@@ -220,14 +217,14 @@ double VEVPMFH_FP::dcontr44 (double **TA, double **TB){
 //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){
+void VEVPMFHSPACE::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);
+  VEVPMFHSPACE::contraction42(B,C,temp);
+  VEVPMFHSPACE::contraction42(A, temp, Tres);
   
 }
 
@@ -238,7 +235,7 @@ void VEVPMFH_FP::contraction442 (double **A, double **B, double *C, double *Tres
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::contr444_reduced (double **A, double **B, double **C){
+double VEVPMFHSPACE::contr444_reduced (double **A, double **B, double **C){
 
   /** Declaration **/
   int i;
@@ -257,10 +254,10 @@ double VEVPMFH_FP::contr444_reduced (double **A, double **B, double **C){
   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);
+  VEVPMFHSPACE::contraction42(A, I2, temp1);
+  VEVPMFHSPACE::contraction42(C, I2, temp2);
+  VEVPMFHSPACE::contraction42(B, temp2, temp3);
+  sol=VEVPMFHSPACE::contraction22(temp1, temp3);
   
   /** Freeing **/
   free(temp1);
@@ -278,21 +275,21 @@ double VEVPMFH_FP::contr444_reduced (double **A, double **B, double **C){
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::contraction444 (double **A, double **B, double **C){
+double VEVPMFHSPACE::contraction444 (double **A, double **B, double **C){
   
   /** Declaration **/
   double **temp;
   double sol;
 
   /** Allocation & Init. **/
-  VEVPMFH_FP::mallocmatrix(&temp, 6, 6);
+  VEVPMFHSPACE::mallocmatrix(&temp, 6, 6);
 
   /**** Start ****/
-  VEVPMFH_FP::contraction44(A, B, temp);
-  sol=VEVPMFH_FP::dcontr44(temp, C);
+  VEVPMFHSPACE::contraction44(A, B, temp);
+  sol=VEVPMFHSPACE::dcontr44(temp, C);
 
   /** Freeing **/
-  VEVPMFH_FP::freematrix(temp, 6);
+  VEVPMFHSPACE::freematrix(temp, 6);
 
   return(sol);
   
@@ -305,7 +302,7 @@ double VEVPMFH_FP::contraction444 (double **A, double **B, double **C){
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::produit_tensoriel (double *TA, double *TB, double **Tres){
+void VEVPMFHSPACE::produit_tensoriel (double *TA, double *TB, double **Tres){
   
   /** Declaration **/
   int i,j;
@@ -328,7 +325,7 @@ void VEVPMFH_FP::produit_tensoriel (double *TA, double *TB, double **Tres){
 //Created: april 2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::max(double a, double b){
+double VEVPMFHSPACE::max(double a, double b){
 
   if (a>b) return(a);
   else return (b);
@@ -336,13 +333,6 @@ double VEVPMFH_FP::max(double a, double b){
 }
 
 
-
-
-
-
-
-
-
 //Description: DECOMPOSITION LU D'UNE MATRICE CAREE D'ORDRE QUELCONQUE
 //Author: Olivier Pierard
 //Created: april 2003
@@ -355,7 +345,7 @@ double VEVPMFH_FP::max(double a, double b){
 
 #define TINY 1.0e-14
 
-void VEVPMFH_FP::ludcmp(double **a, int n, int *indx, double *d, int *error){
+void VEVPMFHSPACE::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
@@ -442,7 +432,7 @@ If one row is null, the function is automatically stopped*/
 //Created: april 2003
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::lubksb (double **a, int n, int *indx, double b[]){
+void VEVPMFHSPACE::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 
@@ -486,24 +476,24 @@ with many zero elements, do it is efficient for use in matrix inversion. */
 //Created: april 2003
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::inverse (double **a, double **y, int *error, int n){
+void VEVPMFHSPACE::inverse (double **a, double **y, int *error, int n){
 
   double d;
   double *col;
   int i,j,*indx;
   double **aoriginal; 
   
-  VEVPMFH_FP::mallocmatrix(&aoriginal, n, n);
+  VEVPMFHSPACE::mallocmatrix(&aoriginal, n, n);
   col = (double*)malloc(n*sizeof(double));
   indx = (int*)malloc (n*sizeof(int));
 
-  for(i=0;i<n;i++) {     /* On mémorise la matrice de départ à VEVPMFH_FP::inverser */
+  for(i=0;i<n;i++) {     /* On mémorise la matrice de départ à VEVPMFHSPACE::inverser */
     for(j=0;j<n;j++) {
       aoriginal[i][j]=a[i][j];
     }
   }
 
-  VEVPMFH_FP::ludcmp(a,n,indx,&d,error);
+  VEVPMFHSPACE::ludcmp(a,n,indx,&d,error);
   if (*error != 0) {
     return;
   }
@@ -513,7 +503,7 @@ void VEVPMFH_FP::inverse (double **a, double **y, int *error, int n){
       col[i]=0.0;
     }
     col[j]=1.0;
-    VEVPMFH_FP::lubksb(a,n,indx,col);
+    VEVPMFHSPACE::lubksb(a,n,indx,col);
     for (i=0;i<n;i++) {
       y[i][j]=col[i];
     }
@@ -526,7 +516,7 @@ void VEVPMFH_FP::inverse (double **a, double **y, int *error, int n){
 
   free(col);
   free(indx);
-  VEVPMFH_FP::freematrix(aoriginal, n);
+  VEVPMFHSPACE::freematrix(aoriginal, n);
 }
 
 
@@ -540,7 +530,7 @@ void VEVPMFH_FP::inverse (double **a, double **y, int *error, int n){
 //** INPUT : A fourth order tensor (6*6)
 //** OUTPUT : An isotropic fourth order tensor
 
-void VEVPMFH_FP::extractiso (double **tensani, double **tensiso)
+void VEVPMFHSPACE::extractiso (double **tensani, double **tensiso)
 {
   double deltaij[6], deltakl[6];
   double I[6][6]={0};
@@ -557,8 +547,8 @@ void VEVPMFH_FP::extractiso (double **tensani, double **tensiso)
     deltakl[i+3] = 0;
   }
 
-  VEVPMFH_FP::mallocmatrix(&I_vol,6,6);
-  VEVPMFH_FP::produit_tensoriel(deltaij, deltakl, I_vol);
+  VEVPMFHSPACE::mallocmatrix(&I_vol,6,6);
+ VEVPMFHSPACE::produit_tensoriel(deltaij, deltakl, I_vol);
 
   for (i=0;i<6;i++) {
     for (j=0;j<6;j++) {
@@ -590,7 +580,7 @@ void VEVPMFH_FP::extractiso (double **tensani, double **tensiso)
     }
   }
 
-  VEVPMFH_FP::freematrix(I_vol,6);
+  VEVPMFHSPACE::freematrix(I_vol,6);
 }
 
 
@@ -601,7 +591,7 @@ void VEVPMFH_FP::extractiso (double **tensani, double **tensiso)
 //Created: 03/04/2020
 //Copyright: See COPYING file that comes with this distribution
  
-void VEVPMFH_FP::mallocmatrix (double ***A, int m, int n){
+void VEVPMFHSPACE::mallocmatrix (double ***A, int m, int n){
   
   /** Dec **/
   int i,j;
@@ -630,7 +620,7 @@ void VEVPMFH_FP::mallocmatrix (double ***A, int m, int n){
 //Created: 01/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::freematrix (double **A, int m){
+void VEVPMFHSPACE::freematrix (double **A, int m){
   
   /** Dec **/
   int i;
@@ -653,7 +643,7 @@ void VEVPMFH_FP::freematrix (double **A, int m){
 //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){
+void VEVPMFHSPACE::mallocmatrix3D (double ****A, int p, int m, int n){
   
   /** Dec **/
   int i,j,k;
@@ -686,7 +676,7 @@ void VEVPMFH_FP::mallocmatrix3D (double ****A, int p, int m, int n){
 //Created: 03/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::freematrix3D (double ***A, int p, int m){
+void VEVPMFHSPACE::freematrix3D (double ***A, int p, int m){
 
   /** Dec **/
   int k, i;
@@ -713,7 +703,7 @@ void VEVPMFH_FP::freematrix3D (double ***A, int p, int m){
 //Created: 03/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::compute_Hooke(double E, double nu, double **C){
+void VEVPMFHSPACE::compute_Hooke(double E, double nu, double **C){
   
   /** Dec **/
   double lambda;
@@ -749,15 +739,12 @@ void VEVPMFH_FP::compute_Hooke(double E, double nu, double **C){
 
 
 
-
-
-
 //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){
+void VEVPMFHSPACE::printvect(double *vect, int n){
   
   /** Dec **/
   int i;
@@ -772,13 +759,12 @@ void VEVPMFH_FP::printvect(double *vect, int n){
 
 
 
-
 //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){
+void VEVPMFHSPACE::fprintvect(FILE *pfile, double *vect, int n){
 
   /** Dec **/
   int i;
@@ -793,14 +779,12 @@ void VEVPMFH_FP::fprintvect(FILE *pfile, double *vect, 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){
+void VEVPMFHSPACE::printmatr(double **matr4, int m, int n){
   
   /** Dec **/
   int i, j;
@@ -824,7 +808,7 @@ void VEVPMFH_FP::printmatr(double **matr4, int m, int n){
 //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){
+void VEVPMFHSPACE::fprintmatr(FILE *pfile, double **matr4, int m, int n){
 
   /** Dec **/
   int i, j;
@@ -848,7 +832,7 @@ void VEVPMFH_FP::fprintmatr(FILE *pfile, double **matr4, int m, int n){
 //Created: 05/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::trace(double *tens2){
+double VEVPMFHSPACE::trace(double *tens2){
 
   /** Dec **/
   int i;
@@ -871,7 +855,7 @@ double VEVPMFH_FP::trace(double *tens2){
 //Created: 05/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::dev(double *tens2, double *dev_tens2){
+void VEVPMFHSPACE::dev(double *tens2, double *dev_tens2){
 
   /** Dec **/
   int i;
@@ -902,7 +886,7 @@ void VEVPMFH_FP::dev(double *tens2, double *dev_tens2){
 //Created: 05/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::vonmises(double *tens2){
+double VEVPMFHSPACE::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));
 
@@ -917,7 +901,7 @@ double VEVPMFH_FP::vonmises(double *tens2){
 //Created: 05/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-double VEVPMFH_FP::eqstrn(double *tens2){
+double VEVPMFHSPACE::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));
 
 }
@@ -931,7 +915,7 @@ double VEVPMFH_FP::eqstrn(double *tens2){
 //Created: 05/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-int VEVPMFH_FP::isnull(double *A, int size){
+int VEVPMFHSPACE::isnull(double *A, int size){
 
   /** Dec & Init. **/
   int i, nil=1;
@@ -961,7 +945,7 @@ int VEVPMFH_FP::isnull(double *A, int size){
 //**Input: A: a fourth order tensor
 //**Output: res: 1 if A is isotropic; 0 otherwise 
 
-int VEVPMFH_FP::isisotropic(double**A){
+int VEVPMFHSPACE::isisotropic(double**A){
 
   /** Dec & Init. **/ 
   double lambda, mu;
@@ -1007,7 +991,7 @@ int VEVPMFH_FP::isisotropic(double**A){
 //Created: 06/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::copyvect(double *A, double *B, int size){
+void VEVPMFHSPACE::copyvect(double *A, double *B, int size){
 
   /** Dec **/
   int i;
@@ -1027,7 +1011,7 @@ void VEVPMFH_FP::copyvect(double *A, double *B, int size){
 //Created: November 2020 
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::compute_I(double **I){
+void VEVPMFHSPACE::compute_I(double **I){
   
   /** Dec **/
   int i,j;
@@ -1051,7 +1035,7 @@ void VEVPMFH_FP::compute_I(double **I){
 //Created: November 2020 
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::compute_Ivol(double **I_vol){
+void VEVPMFHSPACE::compute_Ivol(double **I_vol){
   
   /** Dec **/
   int i,j;
@@ -1071,7 +1055,7 @@ void VEVPMFH_FP::compute_Ivol(double **I_vol){
   }
   
 
-  VEVPMFH_FP::produit_tensoriel(deltaij, deltakl, I_vol);
+  VEVPMFHSPACE::produit_tensoriel(deltaij, deltakl, I_vol);
 
     for (i=0;i<6;i++) {
       for (j=0;j<6;j++) {
@@ -1095,25 +1079,25 @@ void VEVPMFH_FP::compute_Ivol(double **I_vol){
 //Created: 06/04/2020
 //Copyright: See COPYING file that comes with this distribution
 
-void VEVPMFH_FP::compute_Idev(double **I_dev){
+void VEVPMFHSPACE::compute_Idev(double **I_dev){
   
   /** Declaration **/
   double **I, **Ivol;
 
   /** Allocation **/
-  VEVPMFH_FP::mallocmatrix(&I,6,6);
-  VEVPMFH_FP::mallocmatrix(&Ivol,6,6);
+  VEVPMFHSPACE::mallocmatrix(&I,6,6);
+  VEVPMFHSPACE::mallocmatrix(&Ivol,6,6);
 
   /**** Start ****/
-  VEVPMFH_FP::compute_I(I);
-  VEVPMFH_FP::compute_Ivol(Ivol);
+  VEVPMFHSPACE::compute_I(I);
+  VEVPMFHSPACE::compute_Ivol(Ivol);
 
-  VEVPMFH_FP::addtens4(1,I,-1,Ivol,I_dev);
+  VEVPMFHSPACE::addtens4(1,I,-1,Ivol,I_dev);
 
 
   /** Freeing **/
-  VEVPMFH_FP::freematrix(I,6);
-  VEVPMFH_FP::freematrix(Ivol,6);
+  VEVPMFHSPACE::freematrix(I,6);
+  VEVPMFHSPACE::freematrix(Ivol,6);
 
 }
 
@@ -1125,7 +1109,7 @@ void VEVPMFH_FP::compute_Idev(double **I_dev){
 //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){
+void VEVPMFHSPACE::mallocMATPROP(VEVPMFHSPACE::MATPROP* matr, VEVPMFHSPACE::MATPROP* incl, const char* inputfilename){
   FILE *inputfile;
   char c;
   int i;
@@ -1183,15 +1167,12 @@ void VEVPMFH_FP::mallocMATPROP(VEVPMFH_FP::MATPROP* matr, VEVPMFH_FP::MATPROP* i
 
 
 
-
-
-
 //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){
+void VEVPMFHSPACE::readinput (VEVPMFHSPACE::MATPROP *matr, VEVPMFHSPACE::MATPROP *incl, const char* inputfilename){
 
   /** Dec **/
   FILE *inputfile;
@@ -1442,11 +1423,11 @@ void VEVPMFH_FP::readinput (VEVPMFH_FP::MATPROP *matr, VEVPMFH_FP::MATPROP *incl
 //**OUTPUT: 	-R(P)	
 //	   	-dR/dp*/							
 
-void VEVPMFH_FP::compute_R (double p, const VEVPMFH_FP::MATPROP *mat, double *Rp, double *dRdp){
+void VEVPMFHSPACE::compute_R (double p, const VEVPMFHSPACE::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);
+    printf("Invalid parameter(s) for the hardening law in function VEVPMFHSPACE::compute_R:\np: %g, hardmod1: %g, hardmod2: %g, hardexp: %g\n", p, mat->hardmodulus1 , mat->hardmodulus2, mat->hardexp);
     exit(0);
   }
 
@@ -1516,7 +1497,7 @@ void VEVPMFH_FP::compute_R (double p, const VEVPMFH_FP::MATPROP *mat, double *Rp
        
     
   default :
-    printf("parameter chhard invalid in function VEVPMFH_FP::compute_R\n");
+    printf("parameter chhard invalid in function VEVPMFHSPACE::compute_R\n");
     exit(0);
   
   // End switch
@@ -1542,13 +1523,13 @@ void VEVPMFH_FP::compute_R (double p, const VEVPMFH_FP::MATPROP *mat, double *Rp
 //             dgvdstrseq : Derivative of gv w.r.t. strseq
 //             dgvdp : Derivative of gv w.r.t. p */
 
-void VEVPMFH_FP::computevisc (double strseq, double p, const VEVPMFH_FP::MATPROP *mat, double *gv, double *dgvdstrseq, double *dgvdp){
+void VEVPMFHSPACE::computevisc (double strseq, double p, const VEVPMFHSPACE::MATPROP *mat, double *gv, double *dgvdstrseq, double *dgvdp){
   
   /** Dec **/
   double Rp, dRdp, f;
 
   /**** Start ****/
-  VEVPMFH_FP::compute_R(p, mat, &Rp, &dRdp);
+  VEVPMFHSPACE::compute_R(p, mat, &Rp, &dRdp);
   f = strseq - (mat->yldstrs) - Rp;
 
 
@@ -1634,7 +1615,7 @@ void VEVPMFH_FP::computevisc (double strseq, double p, const VEVPMFH_FP::MATPROP
 
 
   default :
-    printf("parameter choice_viscoplastic invalid in function VEVPMFH_FP::computevisc in fonctions_affine.c\n");
+    printf("parameter choice_viscoplastic invalid in function VEVPMFHSPACE::computevisc in fonctions_affine.c\n");
     exit(0);
     
   //End of switch  
@@ -1651,7 +1632,7 @@ void VEVPMFH_FP::computevisc (double strseq, double p, const VEVPMFH_FP::MATPROP
 //Created: 07/04/2020
 //Copyright: See COPYING file that comes with this distribution 
 
-void VEVPMFH_FP::mallocstate(VEVPMFH_FP::STATE *st) {
+void VEVPMFHSPACE::mallocstate(VEVPMFHSPACE::STATE *st) {
 
   /** Dec **/
   int i, j;
@@ -1663,8 +1644,8 @@ void VEVPMFH_FP::mallocstate(VEVPMFH_FP::STATE *st) {
   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);
+  VEVPMFHSPACE::mallocmatrix(&st->Si, VEVPMFHSPACE::max_nbrGi, 6);
+  VEVPMFHSPACE::mallocmatrix(&st->Sigma_Hj, VEVPMFHSPACE::max_nbrKj, 6);
   
 }
   
@@ -1679,7 +1660,7 @@ void VEVPMFH_FP::mallocstate(VEVPMFH_FP::STATE *st) {
 //Created: 07/04/2020
 //Copyright: See COPYING file that comes with this distribution 
 
-void VEVPMFH_FP::freestate(VEVPMFH_FP::STATE st) {
+void VEVPMFHSPACE::freestate(VEVPMFHSPACE::STATE st) {
 
   /** Dec **/
   int i;
@@ -1688,8 +1669,8 @@ void VEVPMFH_FP::freestate(VEVPMFH_FP::STATE st) {
   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);
+  VEVPMFHSPACE::freematrix(st.Si,VEVPMFHSPACE::max_nbrGi);
+  VEVPMFHSPACE::freematrix(st.Sigma_Hj,VEVPMFHSPACE::max_nbrKj);
 
 }
 
@@ -1703,7 +1684,7 @@ void VEVPMFH_FP::freestate(VEVPMFH_FP::STATE st) {
 //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) {
+void VEVPMFHSPACE::copystate(const VEVPMFHSPACE::STATE *A, VEVPMFHSPACE::STATE *B) {
 
   /** Dec **/
   int i,j;
@@ -1717,14 +1698,14 @@ void VEVPMFH_FP::copystate(const VEVPMFH_FP::STATE *A, VEVPMFH_FP::STATE *B) {
   }
 
 
-   for (i=0; i<VEVPMFH_FP::max_nbrGi; i++) {
+   for (i=0; i<VEVPMFHSPACE::max_nbrGi; i++) {
     for (j=0; j<6; j++) {
       B->Si[i][j] = A->Si[i][j];
      }
     }
   
 
-  for (i=0; i<VEVPMFH_FP::max_nbrKj; i++) {
+  for (i=0; i<VEVPMFHSPACE::max_nbrKj; i++) {
     for (j=0; j<6; j++) {
       B->Sigma_Hj[i][j] = A->Sigma_Hj[i][j];
      }
@@ -1742,7 +1723,7 @@ void VEVPMFH_FP::copystate(const VEVPMFH_FP::STATE *A, VEVPMFH_FP::STATE *B) {
 //Created: 07/04/2020
 //Copyright: See COPYING file that comes with this distribution 
 
-void VEVPMFH_FP::fprintstate(FILE *pfile, VEVPMFH_FP::STATE *A) {
+void VEVPMFHSPACE::fprintstate(FILE *pfile, VEVPMFHSPACE::STATE *A) {
 
   /** Dec **/
   int i;
@@ -1753,11 +1734,11 @@ void VEVPMFH_FP::fprintstate(FILE *pfile, VEVPMFH_FP::STATE *A) {
 
   /**** Start ****/
   fprintf(pfile, "Strain:");
-  VEVPMFH_FP::fprintvect(pfile,A->strn,6);
+  VEVPMFHSPACE::fprintvect(pfile,A->strn,6);
   fprintf(pfile, "Stress:");
-  VEVPMFH_FP::fprintvect(pfile,A->strs,6);
+  VEVPMFHSPACE::fprintvect(pfile,A->strs,6);
   fprintf(pfile, "Plastic strain:");
-  VEVPMFH_FP::fprintvect(pfile, A->pstrn, 6);
+  VEVPMFHSPACE::fprintvect(pfile, A->pstrn, 6);
   fprintf(pfile, "p :%f \n",temp);
   fprintf(pfile, "\n");
   
@@ -1765,15 +1746,12 @@ void VEVPMFH_FP::fprintstate(FILE *pfile, VEVPMFH_FP::STATE *A) {
 
 
 
-
-
-
 //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) {
+void VEVPMFHSPACE::printstate(VEVPMFHSPACE::STATE *A) {
 
   /** Dec **/
   int i;
@@ -1784,19 +1762,19 @@ void VEVPMFH_FP::printstate(VEVPMFH_FP::STATE *A) {
   
   /**** Start ****/
   printf("\nStrain:");
-  VEVPMFH_FP::printvect(A->strn,6);
+  VEVPMFHSPACE::printvect(A->strn,6);
   printf("Stress:");
-  VEVPMFH_FP::printvect(A->strs,6);
+  VEVPMFHSPACE::printvect(A->strs,6);
   printf("Plastic strain:");
-  VEVPMFH_FP::printvect(A->pstrn, 6);
+  VEVPMFHSPACE::printvect(A->pstrn, 6);
 
   printf("p : %g \n ",temp);
 
   printf("Si:"); 
-  printmatr(A->Si,VEVPMFH_FP::max_nbrGi,6);
+  printmatr(A->Si,VEVPMFHSPACE::max_nbrGi,6);
 
   printf("sigma_Hj:"); 
-  VEVPMFH_FP::printmatr(A->Sigma_Hj,VEVPMFH_FP::max_nbrKj,6);
+  VEVPMFHSPACE::printmatr(A->Sigma_Hj,VEVPMFHSPACE::max_nbrKj,6);
 
   printf("\n");
 
@@ -1804,33 +1782,28 @@ void VEVPMFH_FP::printstate(VEVPMFH_FP::STATE *A) {
 
 
 
-
-
 //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){
+void fileprintresults(FILE *pfile, double time, VEVPMFHSPACE::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->strs[3]/sqrt(2), st->strs[4]/sqrt(2),st->strs[5]/sqrt(2), st->p, VEVPMFHSPACE::vonmises(st->strs), VEVPMFHSPACE::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]); 
   
 }
 
 
 
-
-
-
 //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) {
+void VEVPMFHSPACE::fprintmatprop(FILE *pfile, VEVPMFHSPACE::MATPROP *A) {
 
   /** Dec **/
   int i;
@@ -1885,16 +1858,12 @@ void VEVPMFH_FP::fprintmatprop(FILE *pfile, VEVPMFH_FP::MATPROP *A) {
 
 
 
-
-
-
-
 //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) {
+void VEVPMFHSPACE::printmatprop(VEVPMFHSPACE::MATPROP *A) {
 
   /**** Start ****/
   printf("\n E(t=0): %g\n",A->E0);
@@ -1916,17 +1885,17 @@ void VEVPMFH_FP::printmatprop(VEVPMFH_FP::MATPROP *A) {
   printf("G_inf: %e\n",A->G_inf);
   printf("number of Gi: %d\n",A->nbrGi);
   printf("Gi:");
-  VEVPMFH_FP::printvect(A->Gi, A->nbrGi);
+  VEVPMFHSPACE::printvect(A->Gi, A->nbrGi);
   printf("gi:");
-  VEVPMFH_FP::printvect(A->gi, A->nbrGi);
+  VEVPMFHSPACE::printvect(A->gi, A->nbrGi);
 
   printf("K0: %e\n",A->K0);
   printf("K_inf: %e\n",A->K_inf);
   printf("number of Kj: %d\n",A->nbrKj);
   printf("Kj:");
-  VEVPMFH_FP::printvect(A->Kj, A->nbrKj);
+  VEVPMFHSPACE::printvect(A->Kj, A->nbrKj);
   printf("kj:");
-  VEVPMFH_FP::printvect(A->kaj, A->nbrKj);
+  VEVPMFHSPACE::printvect(A->kaj, A->nbrKj);
 
 }
 
@@ -1940,7 +1909,7 @@ void VEVPMFH_FP::printmatprop(VEVPMFH_FP::MATPROP *A) {
 //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){
+double VEVPMFHSPACE::compute_Gtilde(const VEVPMFHSPACE::MATPROP *material_prop, double dt){
 
   /** Dec. & Init. **/ 
   int i;
@@ -1957,16 +1926,12 @@ double VEVPMFH_FP::compute_Gtilde(const VEVPMFH_FP::MATPROP *material_prop, doub
 
 
 
-
-
-
-
 //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){
+double VEVPMFHSPACE::compute_Ktilde(const VEVPMFHSPACE::MATPROP *material_prop, double dt){
 
   int i;
   double result=material_prop->K_inf;
@@ -1981,94 +1946,88 @@ double VEVPMFH_FP::compute_Ktilde(const VEVPMFH_FP::MATPROP *material_prop, doub
 
 
 
-
-
 //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){
+void VEVPMFHSPACE::compute_Etilde(const VEVPMFHSPACE::MATPROP *material_prop, double dt, double **E_tilde){
   
   /** Dec. **/
   double Gtilde, Ktilde;
   double **Idev, **Ivol;
 
   /** Alloc. & Init. **/
-  VEVPMFH_FP::mallocmatrix(&Idev,6,6);
-  VEVPMFH_FP::mallocmatrix(&Ivol,6,6);
+  VEVPMFHSPACE::mallocmatrix(&Idev,6,6);
+  VEVPMFHSPACE::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);
+  VEVPMFHSPACE::compute_Idev(Idev);
+  VEVPMFHSPACE::compute_Ivol(Ivol);
 
-  VEVPMFH_FP::addtens4(2*Gtilde,Idev,3*Ktilde,Ivol,E_tilde);
+  VEVPMFHSPACE::addtens4(2*Gtilde,Idev,3*Ktilde,Ivol,E_tilde);
   
   /** Freeing **/
-  VEVPMFH_FP::freematrix(Idev,6);
-  VEVPMFH_FP::freematrix(Ivol,6); 
+  VEVPMFHSPACE::freematrix(Idev,6);
+  VEVPMFHSPACE::freematrix(Ivol,6); 
 
 }
 
 
 
-
-
-
 //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){
+void VEVPMFHSPACE::compute_Einf(const VEVPMFHSPACE::MATPROP *material_prop, double **E_inf){
 
   double **Idev, **Ivol;
 
-  VEVPMFH_FP::mallocmatrix(&Idev,6,6);
-  VEVPMFH_FP::mallocmatrix(&Ivol,6,6);
+  VEVPMFHSPACE::mallocmatrix(&Idev,6,6);
+  VEVPMFHSPACE::mallocmatrix(&Ivol,6,6);
 
-  VEVPMFH_FP::compute_Idev(Idev);
-  VEVPMFH_FP::compute_Ivol(Ivol);
+  VEVPMFHSPACE::compute_Idev(Idev);
+  VEVPMFHSPACE::compute_Ivol(Ivol);
 
 
-  VEVPMFH_FP::addtens4(2*(material_prop->G_inf),Idev,3*(material_prop->K_inf),Ivol,E_inf);
+  VEVPMFHSPACE::addtens4(2*(material_prop->G_inf),Idev,3*(material_prop->K_inf),Ivol,E_inf);
 
-  VEVPMFH_FP::freematrix(Idev,6);
-  VEVPMFH_FP::freematrix(Ivol,6); 
+  VEVPMFHSPACE::freematrix(Idev,6);
+  VEVPMFHSPACE::freematrix(Ivol,6); 
 
 }
 
 
 
 
-
 //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){
+void VEVPMFHSPACE::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);
+  VEVPMFHSPACE::mallocmatrix(&Idev,6,6);
+  VEVPMFHSPACE::mallocmatrix(&Ivol,6,6);
 
   /**** Start ****/
-  VEVPMFH_FP::compute_Ivol(Ivol);
-  VEVPMFH_FP::compute_Idev(Idev);
+  VEVPMFHSPACE::compute_Ivol(Ivol);
+  VEVPMFHSPACE::compute_Idev(Idev);
  
-  VEVPMFH_FP::contraction42(Ivol, dstrn_ve, dstrn_ve_volpart);
-  VEVPMFH_FP::contraction42(Idev, dstrn_ve, dstrn_ve_devpart);
+  VEVPMFHSPACE::contraction42(Ivol, dstrn_ve, dstrn_ve_volpart);
+  VEVPMFHSPACE::contraction42(Idev, dstrn_ve, dstrn_ve_devpart);
 
   /** Freeing **/
-  VEVPMFH_FP::freematrix(Idev,6);
-  VEVPMFH_FP::freematrix(Ivol,6); 
+  VEVPMFHSPACE::freematrix(Idev,6);
+  VEVPMFHSPACE::freematrix(Ivol,6); 
 
 }
 
@@ -2082,7 +2041,7 @@ void VEVPMFH_FP::decomposition_dstrn_ve(double *dstrn_ve, double *dstrn_ve_volpa
 //Created: November 2020
 //Copyright: See COPYING file that comes with this distribution 
 
-void VEVPMFH_FP::compute_Si(const VEVPMFH_FP::MATPROP *material_prop, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn_ve, double **Si){
+void VEVPMFHSPACE::compute_Si(const VEVPMFHSPACE::MATPROP *material_prop, const VEVPMFHSPACE::STATE *tn, double dt, double *dstrn_ve, double **Si){
 
    /** Dec. **/
    int i,j;
@@ -2094,7 +2053,7 @@ void VEVPMFH_FP::compute_Si(const VEVPMFH_FP::MATPROP *material_prop, const VEVP
 
    /**** Start ****/
    // comute dev and vol parts of the strain 
-   VEVPMFH_FP::decomposition_dstrn_ve(dstrn_ve,dstrn_ve_vol,dstrn_ve_dev);
+   VEVPMFHSPACE::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++){
@@ -2120,7 +2079,7 @@ void VEVPMFH_FP::compute_Si(const VEVPMFH_FP::MATPROP *material_prop, const VEVP
 //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){
+void VEVPMFHSPACE::compute_SigmaHj(const VEVPMFHSPACE::MATPROP *material_prop, const VEVPMFHSPACE::STATE *tn, double dt, double *dstrn_ve, double **Sigma_Hj){
 
    /** Dec. **/
    int i,j;
@@ -2131,7 +2090,7 @@ void VEVPMFH_FP::compute_SigmaHj(const VEVPMFH_FP::MATPROP *material_prop, const
    dstrn_ve_dev= (double*) calloc(6,sizeof(double));
 
    /**** Start ****/
-   VEVPMFH_FP::decomposition_dstrn_ve(dstrn_ve,dstrn_ve_vol,dstrn_ve_dev);
+   VEVPMFHSPACE::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++){
@@ -2156,7 +2115,7 @@ void VEVPMFH_FP::compute_SigmaHj(const VEVPMFH_FP::MATPROP *material_prop, const
 //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){
+void VEVPMFHSPACE::compute_sum_Siexp(const VEVPMFHSPACE::MATPROP *material_prop, double dt, double **Si, double *sum_Siexp){
   
   /** Dec. **/
   int i,j;
@@ -2185,7 +2144,7 @@ void VEVPMFH_FP::compute_sum_Siexp(const VEVPMFH_FP::MATPROP *material_prop, dou
 //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){
+void VEVPMFHSPACE::compute_sum_SigmaHjexp(const VEVPMFHSPACE::MATPROP *material_prop, double dt, double **Sigma_Hj, double *sum_SigmaHjexp){
 
   /** Dec. **/
   int i,j;
@@ -2204,68 +2163,212 @@ void VEVPMFH_FP::compute_sum_SigmaHjexp(const VEVPMFH_FP::MATPROP *material_prop
 
 
 
-
-
-
-//Description: COMPUTE Calgo USING PERTURBATION METHOD
+//Description: Constitutive box for VE-VP homogeneous isotropic phase (see B.Miled & I.Doghri, 2011, IJSS.)
 //Author: Mohamed Haddad 
-//Created: September 2022
+//Created: November 2020
 //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;
+//	** 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 VEVPMFHSPACE::ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(const VEVPMFHSPACE::MATPROP *material, const VEVPMFHSPACE::STATE *tn, double dt, double *dstrn, VEVPMFHSPACE::STATE *tau, double **Calgo){
+
+  /**	Declaration	**/
+  /***********************/
+  double **E_tilde;
+  double **E_inf;
+  double *sum_Siexp_tn;
+  double *sum_SigmaHjexp_tn;
+  double *strntn_ve;
+  double *part1; 
+  double *part2; 
+  double *strs_pred;
+  double strseq_pred, f_pred = 0;
+  double R, dRdp;
+  double *S_pred;
+  double *N_pred;
+  double strseq_tau, ptau=0, G_tld;
+  double *strstau;
+  double *pstrntau; 
+  double **tens4;
+  double gv, dgvdstrseq, dgvdp;
+  double kp, kstrs, res, tol = 1.e-8;
+  int m, i;
+  double cp;
+  double *Stau;
+  double trace_strs_pred;
+  double dp, hv;
+  double **NcrossN; 
+  double **dNdstrs; 
+  double **Idev;
+  
+  
+  /**	Memory allocation	**/
+  /******************************/
+  VEVPMFHSPACE::mallocmatrix(&E_tilde, 6, 6); 
+  VEVPMFHSPACE::mallocmatrix(&E_inf, 6, 6);
+  strntn_ve=(double*)calloc(6,sizeof(double));
+  sum_Siexp_tn=(double*)calloc(6,sizeof(double));
+  sum_SigmaHjexp_tn=(double*)calloc(6,sizeof(double));
+  part1 = (double*)calloc(6,sizeof(double));
+  part2 = (double*)calloc(6,sizeof(double));
+  strs_pred = (double*)calloc(6,sizeof(double));
+  strstau = (double*)calloc(6,sizeof(double));
+  Stau = (double*)calloc(6,sizeof(double));
+  pstrntau = (double*)calloc(6,sizeof(double));
+  S_pred = (double*)calloc(6,sizeof(double));
+  N_pred = (double*)calloc(6,sizeof(double));
+  VEVPMFHSPACE::mallocmatrix(&tens4, 6, 6); 
+  VEVPMFHSPACE::mallocmatrix(&NcrossN, 6, 6); 
+  VEVPMFHSPACE::mallocmatrix(&dNdstrs, 6, 6); 
+  VEVPMFHSPACE::mallocmatrix(&Idev, 6, 6); 
+   
+   
+  /**	Begin code	**/
+  /***********************/
+  VEVPMFHSPACE::compute_Etilde(material,dt,E_tilde);
+  VEVPMFHSPACE::compute_Einf(material, E_inf);
+  VEVPMFHSPACE::compute_sum_Siexp(material, dt, tn->Si, sum_Siexp_tn);
+  VEVPMFHSPACE::compute_sum_SigmaHjexp(material, dt, tn->Sigma_Hj, sum_SigmaHjexp_tn);
   
-  int cvarstrn, i;
+  VEVPMFHSPACE::addtens2(1, tn->strn, -1, tn->pstrn, strntn_ve);
+  VEVPMFHSPACE::contraction42(E_inf, strntn_ve,part1);
+  VEVPMFHSPACE::contraction42(E_tilde, dstrn, part2);
   
-  double *vardstrn, *dstrnpert1, *dstrnpert2;
-  double **Calgopert1, **Calgopert2;
+  VEVPMFHSPACE::addtens2(1, part1, 1, part2, part1);
+  VEVPMFHSPACE::addtens2(1, part1, 1, sum_Siexp_tn, part1);
+  VEVPMFHSPACE::addtens2(1, part1, 1, sum_SigmaHjexp_tn, strs_pred);
   
-  /** Alloc. & Init. **/
-  VEVPMFH_FP::mallocstate(&matrtaupert1);
-  VEVPMFH_FP::mallocstate(&matrtaupert2);
+  strseq_pred = VEVPMFHSPACE::vonmises(strs_pred);
+  VEVPMFHSPACE::compute_R (tn->p, material, &R, &dRdp);
+  /** yield function  **/
+  f_pred = strseq_pred - R - (material->yldstrs);
   
-  vardstrn=(double*)calloc(6,sizeof(double));
-  dstrnpert1=(double*)calloc(6,sizeof(double));
-  dstrnpert2=(double*)calloc(6,sizeof(double));
-  VEVPMFH_FP::mallocmatrix(&Calgopert1,6,6);
-  VEVPMFH_FP::mallocmatrix(&Calgopert2,6,6);
   
-  /**** 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*/ 
+  /** VE Predictor **/
+  /******************/
+  if (f_pred <= 0){
+    VEVPMFHSPACE::copyvect(strs_pred, tau->strs, 6);
+    tau->p = tn->p;
+    VEVPMFHSPACE::copyvect(tn->pstrn, tau->pstrn, 6);
+    VEVPMFHSPACE::addtens4(1, E_tilde, 0, E_tilde, Calgo);  
+  }
+  
+  /** VP Corrector  **/
+  /*******************/
+  else{
+  
+    /* compute N_pred */
+     VEVPMFHSPACE::dev(strs_pred, S_pred);
+     VEVPMFHSPACE::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 = VEVPMFHSPACE::compute_Gtilde(material,dt);
+
+    VEVPMFHSPACE::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 = VEVPMFHSPACE::max(fabs(kp),fabs(kstrs/(material->E0)));
+
+
+
+    for (m=1; res > tol; m++) { 
     
-    // Perurbation of dstrn increment 
-    addtens2(1, dstrn, 1, vardstrn, dstrnpert1);                            
-    // Compute the new mechanical state after perturbation 
-    ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(matrprop, matrtn, dt, dstrnpert1, &matrtaupert1, Calgopert1);
+      /* 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));
+
+
+      /*update accumulated plastic strain p and eq_stress */
+      ptau += cp;
+      strseq_tau += -kstrs-(3 * G_tld * cp);
+
+      VEVPMFHSPACE::compute_R (ptau, material, &R, &dRdp);
+      VEVPMFHSPACE::computevisc (strseq_tau, ptau, material, &gv, &dgvdstrseq, &dgvdp);
     
-    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]);
-    }
+      /* 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 = VEVPMFHSPACE::max(fabs(kp),fabs(kstrs/(material->E0)));
 
-  }   
-  
-  /** 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); 
+
+      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("fabs(kp):%g;  fabs(kstrs/(material->yldstrs))=%g; res:%g; tol:%g\n", fabs(kp),fabs(kstrs/(material->yldstrs)), res, tol); 
+        return (0);
+      }    
+   }
   
-}
+ 
 
+    VEVPMFHSPACE::addtens2((2.0*strseq_tau)/3.0, N_pred, 0, N_pred, Stau);
+    trace_strs_pred = VEVPMFHSPACE::trace(strs_pred);     
+   
+   
+    /* Update stress*/
+    for (i=0;i<3;i++) {
+      (tau->strs)[i] = Stau[i] + (1.0/3.0)*trace_strs_pred;
+      (tau->strs)[i+3] = Stau[i+3];
+    }
 
+    /* Update p and plastic strain*/
+    tau->p = ptau;
+    VEVPMFHSPACE::addtens2((ptau-(tn->p)), N_pred, 1, tn->pstrn, tau->pstrn);
+    
+    
+    /* Update Calgo **/
+    VEVPMFHSPACE::computevisc (strseq_tau, ptau, material, &gv, &dgvdstrseq, &dgvdp);
+    dp = ptau - (tn->p);
+    
+    hv = (1 / (dt * dgvdstrseq)) + (3*G_tld) - (dgvdp/dgvdstrseq);
+    VEVPMFHSPACE::produit_tensoriel(N_pred, N_pred, NcrossN);
+    VEVPMFHSPACE::compute_Idev(Idev);
+    VEVPMFHSPACE::addtens4((3/(2*strseq_tau)), Idev, (-1/strseq_tau), NcrossN, dNdstrs);  
+    VEVPMFHSPACE::addtens4(1, E_tilde, -(pow((2*G_tld),2)/hv), NcrossN, tens4);
+    VEVPMFHSPACE::addtens4(1, tens4, -(pow((2*G_tld),2)*((strseq_tau*dp)/(strseq_tau + (3*G_tld*dp)))), dNdstrs, Calgo);
+  }
+  
+
+  /** Freeing **/ 
+  VEVPMFHSPACE::freematrix(E_tilde,6); 
+  VEVPMFHSPACE::freematrix(E_inf,6);
+  free(strntn_ve);
+  free(sum_Siexp_tn); 
+  free(sum_SigmaHjexp_tn);
+  free(part1); 
+  free(part2);
+  free(strs_pred); 
+  free(strstau); 
+  free(pstrntau);
+  free(S_pred); 
+  free(N_pred);
+  free(Stau);
+  VEVPMFHSPACE::freematrix(tens4,6); 
+  VEVPMFHSPACE::freematrix(NcrossN,6);  
+  VEVPMFHSPACE::freematrix(dNdstrs,6); 
+  VEVPMFHSPACE::freematrix(Idev,6);
+   
+  return(1) ;
+
+}
 /** END **/
 /**********/
 
 
-#endif
   
-  
-
diff --git a/NonLinearSolver/materialLaw/fonctions_pratiques.h b/NonLinearSolver/materialLaw/UsefulFunctionsVEVPMFH.h
similarity index 54%
rename from NonLinearSolver/materialLaw/fonctions_pratiques.h
rename to NonLinearSolver/materialLaw/UsefulFunctionsVEVPMFH.h
index 746840ce5a43e4110d13a27810568c42ed758345..63f5341fd840722fab96784c543c97988c6a7c3c 100644
--- a/NonLinearSolver/materialLaw/fonctions_pratiques.h
+++ b/NonLinearSolver/materialLaw/UsefulFunctionsVEVPMFH.h
@@ -1,68 +1,66 @@
-#ifndef FONCTIONS_PRATIQUES_H
-#define FONCTIONS_PRATIQUES_H 1
-
-
-
-
-namespace VEVPMFH_FP{
+/*************************************************************************************************/
+/** Copyright (C) 2022 - See COPYING file that comes with this distribution	                    **/
+/** Author: Mohamed HADDAD (MEMA-iMMC-UCLouvain)							                    **/
+/** Contact: mohamed.haddad@uclouvain.be					                                    **/
+/**											                                                    **/
+/** Description: namespace for storing functions needed for VE-VP material law (small strains)  **/
+/** VE-VP material law : See B.Miled et I.Doghri 2011                                           **/
+/*************************************************************************************************/
+#ifndef USEFULFUNCTIONSVEVPMFH_H
+#define USEFULFUNCTIONSVEVPMFH_H 1
+
+namespace VEVPMFHSPACE{
 
   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 *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) */
+    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 E0;              /* Initial relaxation modulus */
+    double nu;              /* Poisson's ratio */
+    double yldstrs;         /* Yield stress */
+    int hardmodel;          /* Hardening model type */
     double hardmodulus1;
     double hardmodulus2;
     double hardexp;
-    int viscmodel;
+    int viscmodel;          /* Viscous model type */
     double visccoef;
     double viscexp1;
     double viscexp2;
-    double v;
-    double ar;
-    double orient[3];
+    double v;               /* Volume fraction of the phase */
+    double ar;              /* Aspect ratio */
+    double orient[3];       /* Orientation of inclusions  */
     
- 
     // VE properties 
     double G0;
     double G_inf;
-    int nbrGi;  		/* number of shear relaxation times */
-    double *Gi; 		/* shear relaxation weights */ 
-    double *gi; 		/* shear relaxation times */
+    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 */ 
+    int nbrKj;   		        /* number of bulk relaxation times */ 
+    double *Kj;  		        /* bulk relaxation weights */
+    double *kaj; 		        /* bulk relaxation times */ 
   }MATPROP;
 
 
@@ -89,7 +87,7 @@ namespace VEVPMFH_FP{
   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 mallocMATPROP(MATPROP* matr, MATPROP* incl, const char* inputfilename);
   
   void compute_Hooke(double E, double nu, double **C);
 
@@ -132,26 +130,23 @@ namespace VEVPMFH_FP{
   void printmatprop(MATPROP *A);
   void fprintmatprop(FILE *pfile, MATPROP *A);
 
-  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 compute_Gtilde(const MATPROP *material_prop, double dt);
+  double compute_Ktilde(const MATPROP *material_prop, double dt);
+  void compute_Etilde(const MATPROP *material_prop, double dt, double **E_tilde);
+  void compute_Einf(const MATPROP *material_prop, double **E_inf);
   
   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);
+  void compute_Si(const MATPROP *material_prop, const STATE *tn, double dt, double *dstrn_ve, double **Si);
+  void compute_SigmaHj(const MATPROP *material_prop, const 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_Siexp(const MATPROP *material_prop, double dt, double **Si, double *sum_Siexp);
+  void compute_sum_SigmaHjexp(const 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);
 
-/** End namespace **/
-}
+  int ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(const MATPROP *material, const STATE *tn, double dt, double *dstrn, STATE *tau, double **Calgo);
 
+}
  
-#endif
+#endif
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/box_secant.cpp b/NonLinearSolver/materialLaw/box_secant.cpp
deleted file mode 100644
index 19ec142a1ce6fb1c776a5cfbb2add41ff7a1206f..0000000000000000000000000000000000000000
--- a/NonLinearSolver/materialLaw/box_secant.cpp
+++ /dev/null
@@ -1,222 +0,0 @@
-#ifndef BOX_AFFINE_C
-#define BOX_AFFINE_C 1
-
-#include <math.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-#include "box_secant.h"
-
-
-//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(const VEVPMFH_FP::MATPROP *material, const VEVPMFH_FP::STATE *tn, double dt, double *dstrn, STATE *tau, double **Calgo){
-
-  /**	Declaration	**/
-  /***********************/
-  double **E_tilde;
-  double **E_inf;
-  double *sum_Siexp_tn;
-  double *sum_SigmaHjexp_tn;
-  double *strntn_ve;
-  double *part1; 
-  double *part2; 
-  double *strs_pred;
-  double strseq_pred, f_pred = 0;
-  double R, dRdp;
-  double *S_pred;
-  double *N_pred;
-  double strseq_tau, ptau=0, G_tld;
-  double *strstau;
-  double *pstrntau; 
-  double **tens4;
-  double gv, dgvdstrseq, dgvdp;
-  double kp, kstrs, res, tol = 1.e-8;
-  int m, i;
-  double cp;
-  double *Stau;
-  double trace_strs_pred;
-  double dp, hv;
-  double **NcrossN; 
-  double **dNdstrs; 
-  double **Idev;
-  
-  
-  /**	Memory allocation	**/
-  /******************************/
-  mallocmatrix(&E_tilde, 6, 6); 
-  mallocmatrix(&E_inf, 6, 6);
-  strntn_ve=(double*)calloc(6,sizeof(double));
-  sum_Siexp_tn=(double*)calloc(6,sizeof(double));
-  sum_SigmaHjexp_tn=(double*)calloc(6,sizeof(double));
-  part1 = (double*)calloc(6,sizeof(double));
-  part2 = (double*)calloc(6,sizeof(double));
-  strs_pred = (double*)calloc(6,sizeof(double));
-  strstau = (double*)calloc(6,sizeof(double));
-  Stau = (double*)calloc(6,sizeof(double));
-  pstrntau = (double*)calloc(6,sizeof(double));
-  S_pred = (double*)calloc(6,sizeof(double));
-  N_pred = (double*)calloc(6,sizeof(double));
-  mallocmatrix(&tens4, 6, 6); 
-  mallocmatrix(&NcrossN, 6, 6); 
-  mallocmatrix(&dNdstrs, 6, 6); 
-  mallocmatrix(&Idev, 6, 6); 
-   
-   
-  /**	Begin code	**/
-  /***********************/
-  
-  compute_Etilde(material,dt,E_tilde);
-  compute_Einf(material, E_inf);
-  compute_sum_Siexp(material, dt, tn->Si, sum_Siexp_tn);
-  compute_sum_SigmaHjexp(material, dt, tn->Sigma_Hj, sum_SigmaHjexp_tn);
-  
-  addtens2(1, tn->strn, -1, tn->pstrn, strntn_ve);
-  contraction42(E_inf, strntn_ve,part1);
-  contraction42(E_tilde, dstrn, part2);
-  
-  addtens2(1, part1, 1, part2, part1);
-  addtens2(1, part1, 1, sum_Siexp_tn, part1);
-  addtens2(1, part1, 1, sum_SigmaHjexp_tn, strs_pred);
-  
-  strseq_pred = vonmises(strs_pred);
-  compute_R (tn->p, material, &R, &dRdp);
-  /** yield function  **/
-  f_pred = strseq_pred - R - (material->yldstrs);
-  
-  
-  /** VE Predictor **/
-  /******************/
-  if (f_pred <= 0){
-    copyvect(strs_pred, tau->strs, 6);
-    tau->p = tn->p;
-    copyvect(tn->pstrn, tau->pstrn, 6);
-    addtens4(1, E_tilde, 0, E_tilde, Calgo);  
-  }
-  
-  /** VP Corrector  **/
-  /*******************/
-  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_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));
-
-
-      /*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("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.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++) {
-      (tau->strs)[i] = Stau[i] + (1.0/3.0)*trace_strs_pred;
-      (tau->strs)[i+3] = Stau[i+3];
-    }
-
-    /* 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_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_tau*dp)/(strseq_tau + (3*G_tld*dp)))), dNdstrs, Calgo);
-  }
-  
-
-  /** Freeing **/ 
-  freematrix(E_tilde,6); 
-  freematrix(E_inf,6);
-  free(strntn_ve);
-  free(sum_Siexp_tn); 
-  free(sum_SigmaHjexp_tn);
-  free(part1); 
-  free(part2);
-  free(strs_pred); 
-  free(strstau); 
-  free(pstrntau);
-  free(S_pred); 
-  free(N_pred);
-  free(Stau);
-  freematrix(tens4,6); 
-  freematrix(NcrossN,6);  
-  freematrix(dNdstrs,6); 
-  freematrix(Idev,6);
-   
-  return(1) ;
-
-}
-
-#endif
-
-
diff --git a/NonLinearSolver/materialLaw/box_secant.h b/NonLinearSolver/materialLaw/box_secant.h
deleted file mode 100644
index c788822d1eda62e408a4509ebdbef6ad470a442a..0000000000000000000000000000000000000000
--- a/NonLinearSolver/materialLaw/box_secant.h
+++ /dev/null
@@ -1,11 +0,0 @@
-#ifndef BOX_SECANT_H
-#define BOX_SECANT_H 1
-
-#include "fonctions_pratiques.h"
-
-using namespace VEVPMFH_FP;
-
-
-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/mlawVEVPMFH.cpp b/NonLinearSolver/materialLaw/mlawVEVPMFH.cpp
index 76ed95bba2d0f0ab946822d300d05ad83a919e74..ae2a62eebff9104f58610772bc89179f23b29d16 100644
--- a/NonLinearSolver/materialLaw/mlawVEVPMFH.cpp
+++ b/NonLinearSolver/materialLaw/mlawVEVPMFH.cpp
@@ -1,46 +1,44 @@
-//
-// C++ Interface: material law
-//
-// Description: non-local damage elasto-plastic law
-//
-//
-// Author:  <L. NOELS>, (C) 2011
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
+/*********************************************************************************************/
+/** Copyright (C) 2022 - See COPYING file that comes with this distribution	                **/
+/** Author: Mohamed HADDAD (MEMA-iMMC-UCLouvain)							                              **/
+/** Contact: mohamed.haddad@uclouvain.be					                                      		**/
+/**											                                                                  	**/
+/** C++ Interface: material law                                                             **/
+/** Description: ViscoElastic-ViscoPlastic material law (small strains)             				**/
+/** See B.Miled et I.Doghri 2011                                                            **/
+/**********************************************************************************************/
+
 #include <math.h>
 #include "MInterfaceElement.h"
 #include <fstream>
-#include "mlawVEVPMFH.h"
 #include "STensorOperations.h"
-#include "fonctions_pratiques.h"
-#include "box_secant.h"
 #include <cstring>
-
 #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);
+#include "mlawVEVPMFH.h"
+#include "UsefulFunctionsVEVPMFH.h"
+
+mlawVEVPMFH::mlawVEVPMFH(const int num, const double rho, const char *propName) : materialLaw(num,true), _rho(rho)
+{
         
-        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);*/
-}
-  
-  
+    std::ifstream in(propName);
+      if(!in)
+      {
+          Msg::Error("Cannot open the file %s! Maybe is missing",propName);
+          std::exit(EXIT_FAILURE);
+      }
+      
+    inputfilename = propName;
+    VEVPMFHSPACE::readinput(&matrprop, &inclprop, propName);
+};
   
   
-mlawVEVPMFH::mlawVEVPMFH(const mlawVEVPMFH &source) : materialLaw(source), _rho(source._rho){
-
-     // create and copy
+mlawVEVPMFH::mlawVEVPMFH(const mlawVEVPMFH &source) : materialLaw(source), _rho(source._rho)
+{
      _rho = source._rho;
      inputfilename = source.inputfilename;
      
-     VEVPMFH_FP::readinput (&matrprop, &inclprop, source.inputfilename);
+     VEVPMFHSPACE::readinput (&matrprop, &inclprop, source.inputfilename);
      matrprop.E0 = (source.matrprop).E0;
      matrprop.nu = (source.matrprop).nu;
     
@@ -81,19 +79,16 @@ mlawVEVPMFH::mlawVEVPMFH(const mlawVEVPMFH &source) : materialLaw(source), _rho(
        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
+     //copy material structures
      _rho = src->_rho;
      matrprop.E0 = (src->matrprop).E0;
      matrprop.nu = (src->matrprop).nu;
@@ -135,39 +130,35 @@ mlawVEVPMFH& mlawVEVPMFH::operator=(const materialLaw &source){
        matrprop.kaj[i] = (src->matrprop).kaj[i];
      }
   
-  
-  // copy incl prop remaining 
-  return *this;
-}
-
-
 
+    return *this;
+};
 
 
-mlawVEVPMFH::~mlawVEVPMFH(){
+mlawVEVPMFH::~mlawVEVPMFH()
+{
+    //Destructor
     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();
      
   if(ips != NULL) delete ips;
   ips = new IP3State(state_,ivi,iv1,iv2);
-}
-
-
+};
 
 
-double mlawVEVPMFH::soundSpeed() const{
+double mlawVEVPMFH::soundSpeed() const
+{
   double nu = matrprop.nu;
   double E = matrprop.E0;
   double mu = (E*nu)/((1+nu)*(1-2*nu));
@@ -175,7 +166,7 @@ double mlawVEVPMFH::soundSpeed() const{
   
   return sqrt(E*factornu/_rho);
 
-}
+};
 
 
 
@@ -185,21 +176,20 @@ void mlawVEVPMFH::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P
                                       STensor63* dCalgdeps) const{
                                       
   
-  double **Calgo; 	VEVPMFH_FP::mallocmatrix(&Calgo,6,6);
+  double **Calgo; 	VEVPMFHSPACE::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);
-  const VEVPMFH_FP::STATE* statetn= ipvprev->getSTATE();
-  VEVPMFH_FP::STATE* statetau= ipvcur->getSTATE();
+  const VEVPMFHSPACE::STATE* statetn= ipvprev->getSTATE();
+  VEVPMFHSPACE::STATE* statetau= ipvcur->getSTATE();
   
-  // Set time step
+  //time step
   dt = _timeStep;
-  //printf("\n time step in %g \n", dt);
   
-  // Compute logarithmic strain increment from deformation gradient F
+  //true strain from the deformation gradient
   for(int i=0;i<3;i++){
     dstrn[i]=Fn(i,i)-F0(i,i);
   }
@@ -208,31 +198,18 @@ void mlawVEVPMFH::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P
   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");
-  }
-  
-  */
-  // Compute state at t = 0
+  //case t=0
   if (dt == 0){
-  
+    //strs(t=0)
     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);
+    VEVPMFHSPACE::compute_Hooke(matrprop.E0, matrprop.nu, Calgo);
+    VEVPMFHSPACE::contraction42(Calgo, dstrn, dstrs);
+    VEVPMFHSPACE::addtens2(1, statetn->strs, 1, dstrs, statetau->strs);
     
-    VEVPMFH_FP::addtens2(1, statetn->strn, 1, dstrn, statetau->strn);
-    VEVPMFH_FP::copystate(statetau, &(ipvcur->_state));
+    VEVPMFHSPACE::addtens2(1, statetn->strn, 1, dstrn, statetau->strn);
+    VEVPMFHSPACE::copystate(statetau, &(ipvcur->_state));
     
-      
-    //convert str to P  :order xx yy zz xy xz yz
+    //convert strs(t=0) to P(t=0)  
     for(int i=0;i<3;i++){
       P(i,i) = statetau->strs[i];
     }
@@ -243,141 +220,8 @@ void mlawVEVPMFH::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P
     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
+    //convert calgo(t=0) to Tangent(t=0) 
    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);
-    
-    }
-   
-  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);
@@ -467,55 +311,183 @@ 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.;
+    //freeing variables at t=0
+    free(dstrs);
+  }
+  
+  //case t!=0
+  else{
+    //update the stress, p, Calgo
+    int error= VEVPMFHSPACE::ConstitutiveBoxViscoElasticViscoPlasticSmallStrains(&matrprop, statetn, dt, dstrn, statetau, Calgo);
+
+    if(error){//converged 
+      //update viscous stresses  
+      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));
     
+      VEVPMFHSPACE::addtens2(-1, statetn->pstrn, 1, statetau->pstrn, dstrnp); 
+      VEVPMFHSPACE::addtens2(-1, dstrnp, 1, dstrn, dstrn_ve); 	    
+
+      VEVPMFHSPACE::decomposition_dstrn_ve(dstrn_ve, dstrn_ve_vol, dstrn_ve_dev);
+     
+      VEVPMFHSPACE::compute_SigmaHj(&matrprop, statetn, dt, dstrn_ve_vol, statetau->Sigma_Hj);
+      VEVPMFHSPACE::compute_Si(&matrprop, statetn, dt, dstrn_ve_dev, statetau->Si);
     
-    }
+      VEVPMFHSPACE::addtens2(1, statetn->strn, 1, dstrn, statetau->strn);
+      VEVPMFHSPACE::copystate(statetau, &(ipvcur->_state));
+      
+      //convert strs to P 
+      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 Calgo to Tangent 
+     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);
+      }
+
+      // compute VE & VP energies
+      //double* ve_strn; ve_strn=(double*)malloc(6*sizeof(double));
+      //VEVPMFHSPACE::addtens2(1,statetau->strn,-1,statetau->pstrn,ve_strn);
+      //ipvcur->_elasticEne=0.5*VEVPMFHSPACE::contraction22(statetau->strs,ve_strn);
+      //ipvcur->_plasticEne=0.5*VEVPMFHSPACE::contraction22(statetau->strs,statetau->pstrn);
+      //free(ve_strn);
+
+      free(dstrnp); 
+      free(dstrn_ve); 
+      free(dstrn_ve_vol); 
+      free(dstrn_ve_dev);
+        
+    }
     
-    //constbox did not converge
-    // to exist NR and performed next step with a smaller incremenet
-    else{
+    else{//no convergence of the constitutive box
+      // to exist NR and performed next step with a smaller incremenet
       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);
+  //freeing
+  VEVPMFHSPACE::freematrix(Calgo,6);
   free(dstrn);
   
 }
 
 
-
-
-
 const double mlawVEVPMFH::bulkModulus() const{
   return matrprop.K0;
 }
 
-
+const double mlawVEVPMFH::YoungModulus() const{
+  return 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 ;
+  return matrprop.G0 ;
 }
 
 
-
-
 const double mlawVEVPMFH::poissonRatio() const{
   return matrprop.nu;
 
 }
 
-
diff --git a/NonLinearSolver/materialLaw/mlawVEVPMFH.h b/NonLinearSolver/materialLaw/mlawVEVPMFH.h
index 83dd5ea54bdf3d51c55a2ad752c05cbfbac0eed1..91b530de4401bde2bec3354d0889686feabe342a 100644
--- a/NonLinearSolver/materialLaw/mlawVEVPMFH.h
+++ b/NonLinearSolver/materialLaw/mlawVEVPMFH.h
@@ -1,33 +1,31 @@
-//
-// C++ Interface: material law
-//
-// Description: non local damage elasto-plastic law it is a "virtual implementation" you have to define a law
-//              in your project which derive from this law. you have to do THE SAME for your IPVariable
-//              which have to derive from IPJ2linear (all data in this ipvarible have name beggining by _nld...
-//
-// Author:  <L. Noels>, (C) 2011
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
+/*********************************************************************************************/
+/** Copyright (C) 2022 - See COPYING file that comes with this distribution	                **/
+/** Author: Mohamed HADDAD (MEMA-iMMC-UCLouvain)							                              **/
+/** Contact: mohamed.haddad@uclouvain.be					                                      		**/
+/**											                                                                  	**/
+/** C++ Interface: material law                                                             **/
+/** Description: ViscoElastic-ViscoPlastic material law (small strains)             				**/
+/** See B.Miled et I.Doghri 2011                                                            **/
+/**********************************************************************************************/
+
+
 #ifndef MLAWVEVPMFH_H_
 #define MLAWVEVPMFH_H_
 #include "mlaw.h"
 #include "STensor3.h"
 #include "STensor43.h"
+
 #include "ipVEVPMFH.h"
-#include "fonctions_pratiques.h"
-//#include "box_secant.h"
+#include "UsefulFunctionsVEVPMFH.h"
 
 class mlawVEVPMFH : public materialLaw{
   
   protected:
-    double _rho; //maximal mu for anisotropic
-
-    // to be compatible with umat: use same variable
-    VEVPMFH_FP::MATPROP matrprop;
-    VEVPMFH_FP::MATPROP inclprop;
-    const char* inputfilename;
+    double _rho;                                              // maximal mu for anisotropic
+    VEVPMFHSPACE::MATPROP matrprop;                           // structures to store matrix & inclusions VE-VP properties
+    VEVPMFHSPACE::MATPROP inclprop;
+    const char* inputfilename;                                // input file
+  
   public:
     mlawVEVPMFH(const int num, const double rho, const char *propName);
 
@@ -37,15 +35,15 @@ class mlawVEVPMFH : public materialLaw{
   virtual ~mlawVEVPMFH();
   virtual materialLaw* clone() const {return new mlawVEVPMFH(*this);};
   virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{}; // do nothing
-  
-  // function of materialLaw
-  virtual matname getType() const{return materialLaw::vevpmfh;}
+
+  virtual matname getType() const{return materialLaw::vevpmfh;}       // function of 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 double soundSpeed() const; 	
+  virtual const double YoungModulus() 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;
@@ -53,23 +51,19 @@ class mlawVEVPMFH : public materialLaw{
   // specific function
   public:
   virtual void constitutive(
-                            const STensor3& F0,         // initial deformation gradient (input @ time n)
-                            const STensor3& Fn,         // updated deformation gradient (input @ time n+1)
-                            STensor3 &P,                // updated 1st Piola-Kirchhoff stress tensor (output)
-                                                        // contains the initial values on input
-                            const IPVariable *q0,       // array of initial internal variable
-                            IPVariable *q1,             // updated array of internal variable (in ipvcur on output),
-                            STensor43 &Tangent,         // constitutive tangents (output)
-                            const bool stiff,            // if true compute the tangents
+                            const STensor3& F0,                     // initial deformation gradient (input @ time n)
+                            const STensor3& Fn,                     // updated deformation gradient (input @ time n+1)
+                            STensor3 &P,                            // updated 1st Piola-Kirchhoff stress tensor (output)
+                                                                    // contains the initial values on input
+                            const IPVariable *q0,                   // array of initial internal variable
+                            IPVariable *q1,                         // updated array of internal variable (in ipvcur on output),
+                            STensor43 &Tangent,                     // constitutive tangents (output)
+                            const bool stiff,                       // if true compute the tangents
                             STensor43* elasticTangent = NULL, 
                             const bool dTangent =false,
                             STensor63* dCalgdeps = NULL) const;
-
-
-
-  /*int getNsdv() const { return nsdv;}*/
   
   #endif // SWIG
 };
 
-#endif // MLAWNONLOCALDAMAGE_H_
+#endif 
diff --git a/NonLinearSolver/modelReduction/DeepMaterialNetworks.h b/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
index 93f08583e1b59b562413076e92ff2d0e38582529..e639391090ec99f10693738b4f0ba2ee41a29585 100644
--- a/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
+++ b/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
@@ -160,6 +160,7 @@ class DeepMaterialNetwork
     int getTag() const{return _tag;}
     bool stress(const STensor3& F, STensor3& P, bool stiff, STensor43& L, int maxIter, double tol, double absTol, 
                 bool messageActive=false, bool reuseStiffnessMatrix=false);
+    virtual void setTime(double tcur,double timeStep)=0;
     bool initialElasticTangentByStress(fullMatrix<double>& Lmat, 
                             bool stiff=false, 
                             const CoefficientReduction* coeffReduc=NULL,
diff --git a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp
index 89dcea0ad3eebd147c97ebf0e2382f34dae234bf..bf6809c4c1d58a8fb7e34c885e2e0394ad122bd3 100644
--- a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp
+++ b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp
@@ -840,7 +840,7 @@ void TrainingDeepMaterialNetworkNonLinear::setMaximumStrainStep(double val)
 
 int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strainPath, PathSTensor3& stressPath, bool messageActive, bool stiff, 
                                                       bool extractPhaseDefo, PathSTensor3* phaseDefoPath,
-                                                      bool extractIPComp, std::vector<PathDouble>* IpCompPath)
+                                                      bool extractIPComp, std::vector<PathDouble>* IpCompPath, PathDouble* timeSteps)
 {
   double tstart = Cpu();
   //
@@ -903,6 +903,8 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
       }
     }
   }
+
+  double tn=0;
   for (int ip = 0; ip <pathLength; ip++)
   {    
     static STensor3 F0, dF;
@@ -920,7 +922,13 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
     dF -= F0;
     //
     double strainStep = sqrt(STensorOperation::doubledot(dF,dF));
+    
     double timeStep = 1.;
+    if(timeSteps!=NULL)
+    {
+        timeStep=timeSteps->getValue(ip);
+    }
+
     if (strainStep > _maximumStrainStep)
     {
       int numsubstep =1+trunc(strainStep/_maximumStrainStep);
@@ -928,12 +936,15 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
       {
         Msg::Info("using %d substeps",numsubstep);
       }
-      timeStep = 1./double(numsubstep);
+      timeStep = timeStep/double(numsubstep);
     };
     //
     double initialTimeStep = timeStep;
-    double tprev =0.;
-    double tcur = 0;
+    double tprev=tn;
+    double tcur=tn;
+    double ttau=tn+timeStep;
+    //Msg::Info("tn=%g, time step=%g, tn+1=%g, tprev=%g, tcur=%g",tn,timeStep,ttau,tprev,tcur);
+
     STensor3& Pcur = stressPath.getValue(ip);
     std::vector<STensor3>* DPcurDW = NULL;
     std::vector<STensor3>* DPcurDAlpha = NULL;
@@ -971,18 +982,21 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
     int numberOfReducedSteps=0;
     
     int iter = 0;
-    while (fabs(tcur-1.) > 1e-6)    
+    double timeStep0 = timeStep;
+    while (fabs(tcur-ttau) > 1e-6)    
     {
       tcur = tprev + timeStep;
+      tn=tcur;
       //  
       static STensor3 Fcur;
       Fcur = F0;
-      Fcur.daxpy(dF,tcur);
+      Fcur.daxpy(dF,timeStep/timeStep0);
       //
       if (messageActive)
       {
         Msg::Info("iteration %d time = %f: ",iter, tcur);
       }
+      _DMN->setTime(tcur,timeStep);
       bool ok = _DMN->stress(Fcur,Pcur,false,L,maxIter,tol,absTol,messageActive);
       if (!ok)
       {
@@ -1006,6 +1020,7 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
       else
       {
         // everything is good
+        Msg::Info("tn=%g, time step=%g, tn+1=%g",tn,timeStep,ttau);
         tprev = tcur;
         if (timeStep < initialTimeStep)
         {
@@ -1021,9 +1036,9 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
           }
         }
         
-        if (tcur+timeStep > 1.)
+        if (tcur+timeStep > ttau)
         {
-          timeStep = 1.-tprev;
+          timeStep = ttau-tprev;
         };
         _DMN->nextStep();
         //
@@ -2026,11 +2041,11 @@ void TrainingDeepMaterialNetworkNonLinear::reinitializeModel()
   _interaction->reInitialize();
 };
 
-int evaluatePath(DeepMaterialNetwork* dmn, double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath)
+int evaluatePath(DeepMaterialNetwork* dmn, double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath, PathDouble* timeSteps)
 {
   CoefficientReductionVoid reduction;
   TrainingDeepMaterialNetworkNonLinear offTraining(*dmn,reduction);
   offTraining.setMaximumStrainStep(maximumStrainStep);
-  int nb = offTraining.evaluatePath(strainPath,stressPath,true,false);
+  int nb = offTraining.evaluatePath(strainPath,stressPath,false,false,false,NULL,false,NULL,timeSteps);
   return nb;
 };
\ No newline at end of file
diff --git a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h
index 4ec712ae83015e854bdbb8cef0803f88b2763272..69b6289a1d574313f441e5b9d21e6dd3cbfc1d8e 100644
--- a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h
+++ b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h
@@ -247,7 +247,7 @@ class TrainingDeepMaterialNetworkNonLinear : public TrainingArbitraryDeepMateria
 	 */
     int evaluatePath(const PathSTensor3& strainPath, PathSTensor3& stressPath, bool messageActive, bool stiff=false, 
                           bool extractPhaseDefo=false, PathSTensor3* phaseDefoPath=NULL,
-                          bool extractIPComp = false, std::vector<PathDouble>* IpCompPath=NULL);
+                          bool extractIPComp = false, std::vector<PathDouble>* IpCompPath=NULL, PathDouble* timeSteps=NULL);
     void computeErrorVec(std::vector<fullVector<double> >& allErrors, const std::vector<LossPathMeasure*>& lossMeasures, const PathSTensor3& Fref, const PathSTensor3& Pref, 
                                 bool extractPhaseDefo=false, PathSTensor3* phaseDefoPathRef=NULL,
                                 bool extractIPComp = false, std::vector<PathDouble>* IpCompPathRef=NULL);
@@ -283,6 +283,6 @@ class TrainingDeepMaterialNetworkNonLinear : public TrainingArbitraryDeepMateria
 };
 
 
-int evaluatePath(DeepMaterialNetwork* dmn, double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath);
+int evaluatePath(DeepMaterialNetwork* dmn, double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath, PathDouble* timeSteps=NULL);
 
 #endif //_TRAININGDEEPMATERIALNETWORKSNONLINEAR_H_ 
\ No newline at end of file
diff --git a/dG3D/src/dG3DDeepMaterialNetworks.cpp b/dG3D/src/dG3DDeepMaterialNetworks.cpp
index 16be25cffe834d852c4dc8bffcfbcb25f89985fe..a77e3918a17198827f1e44beb8b07aaa6d01127e 100644
--- a/dG3D/src/dG3DDeepMaterialNetworks.cpp
+++ b/dG3D/src/dG3DDeepMaterialNetworks.cpp
@@ -191,6 +191,14 @@ void dG3DDeepMaterialNetwork::addLaw(int matIndex, materialLaw& l1)
   _mlaw[matIndex] = &l1;
 };
 
+void dG3DDeepMaterialNetwork::setTime(double tcur,double timeStep)
+{
+  for (std::map<int,materialLaw*>::iterator it = _mlaw.begin(); it != _mlaw.end(); it++)
+  {
+    it->second->setTime(tcur,timeStep);
+  }
+};
+
 void dG3DDeepMaterialNetwork::setMacroSolver(const nonLinearMechSolver* sv)
 {
   for (std::map<int,materialLaw*>::iterator it = _mlaw.begin(); it != _mlaw.end(); it++)
diff --git a/dG3D/src/dG3DDeepMaterialNetworks.h b/dG3D/src/dG3DDeepMaterialNetworks.h
index 2bc4458513c5f3c2c87a803db22781c5564f53a2..9d7e713d0476441358a27a0d0c1cb0c2b3b2d901 100644
--- a/dG3D/src/dG3DDeepMaterialNetworks.h
+++ b/dG3D/src/dG3DDeepMaterialNetworks.h
@@ -79,6 +79,7 @@ class dG3DDeepMaterialNetwork : public DeepMaterialNetwork
     dG3DDeepMaterialNetwork(const std::string fname);
     virtual ~dG3DDeepMaterialNetwork();
     void addLaw(int matIndex, materialLaw& l1);
+    void setTime(double tcur,double timeStep);
   #ifndef SWIG
   protected:
     virtual void setMacroSolver(const nonLinearMechSolver* sv);
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index cd7f56808858eb9ad7782cf44aef4ebeea4c415d..ddc830c7edf6862749bd90996014216a78ce327d 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -2980,14 +2980,14 @@ 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)
 {
@@ -2999,6 +2999,7 @@ VEVPMFHDG3DIPVariable::VEVPMFHDG3DIPVariable(const VEVPMFHDG3DIPVariable &source
 }
 
 
+
 VEVPMFHDG3DIPVariable& VEVPMFHDG3DIPVariable::operator=(const IPVariable &source)
 {
   dG3DIPVariable::operator=(source);
@@ -3026,6 +3027,8 @@ VEVPMFHDG3DIPVariable& VEVPMFHDG3DIPVariable::operator=(const IPVariable &source
   return *this;
 }
 
+
+
 double VEVPMFHDG3DIPVariable::get(const int comp) const
 {
   double val=_VEVPipv->get(comp);
@@ -3037,30 +3040,28 @@ double VEVPMFHDG3DIPVariable::get(const int comp) const
 }
 
 
+
+
 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)
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 44eef4a5f6e7f5a33d4b73ee034a99d02ff1d824..6b195fb68d656556caa7dc696e8fba91410efd88 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -2696,10 +2696,7 @@ class localDamageIsotropicElasticityDG3DIPVariable : public dG3DIPVariable
 };
 
 
-
-
-
-class VEVPMFHDG3DIPVariable : public dG3DIPVariable // or store data in a different way
+class VEVPMFHDG3DIPVariable : public dG3DIPVariable
 {
  protected:
   IPVEVPMFH* _VEVPipv;
@@ -2730,7 +2727,6 @@ class VEVPMFHDG3DIPVariable : public dG3DIPVariable // or store data in a differ
 
 
 
-
 class ViscoelasticDG3DIPVariable : public dG3DIPVariable
 {
  protected:
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index bb1f96e5a58056f33e28d35c213f9be63d353855..e733835895991618fac38a692dc13ffdc65cac1b 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -3983,20 +3983,15 @@ 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);*/
-
+  fillElasticStiffness(_glaw->YoungModulus(), _glaw->poissonRatio(), elasticStiffness);
 }
 
 
+
 VEVPMFHDG3DMaterialLaw::~VEVPMFHDG3DMaterialLaw()
 {
   if (_glaw !=NULL) delete _glaw;
@@ -4004,8 +3999,8 @@ VEVPMFHDG3DMaterialLaw::~VEVPMFHDG3DMaterialLaw()
 };
 
 
-VEVPMFHDG3DMaterialLaw::VEVPMFHDG3DMaterialLaw(const VEVPMFHDG3DMaterialLaw &source) :
-                                                             dG3DMaterialLaw(source)
+
+VEVPMFHDG3DMaterialLaw::VEVPMFHDG3DMaterialLaw(const VEVPMFHDG3DMaterialLaw &source):dG3DMaterialLaw(source)
 {
 	_glaw = NULL;
 	if (source._glaw != NULL){
@@ -4014,6 +4009,8 @@ VEVPMFHDG3DMaterialLaw::VEVPMFHDG3DMaterialLaw(const VEVPMFHDG3DMaterialLaw &sou
 
 }
 
+
+
 void VEVPMFHDG3DMaterialLaw::setTime(const double t,const double dtime){
   dG3DMaterialLaw::setTime(t,dtime);
   _glaw->setTime(t,dtime);
@@ -4021,6 +4018,7 @@ void VEVPMFHDG3DMaterialLaw::setTime(const double t,const double dtime){
 }
 
 
+
 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
@@ -4037,6 +4035,8 @@ void VEVPMFHDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce,
   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;
@@ -4052,12 +4052,16 @@ void VEVPMFHDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForc
 
 
 
-double VEVPMFHDG3DMaterialLaw::soundSpeed() const{return _glaw->soundSpeed();}
+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);
+  VEVPMFHDG3DIPVariable* ipvcur; 
+  const VEVPMFHDG3DIPVariable* ipvprev; 
 
   FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
   if(ipvtmp !=NULL)
@@ -4073,15 +4077,14 @@ void VEVPMFHDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariabl
   }
   _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);
+  VEVPMFHDG3DIPVariable* ipvcur; 
+  const VEVPMFHDG3DIPVariable* ipvprev;
 
   FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
   if(ipvtmp !=NULL)
@@ -4097,15 +4100,12 @@ void VEVPMFHDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, con
     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());
 
@@ -4114,19 +4114,18 @@ void VEVPMFHDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, con
 }
 
 
-double VEVPMFHDG3DMaterialLaw::scaleFactor() const{return _glaw->scaleFactor();};
+double VEVPMFHDG3DMaterialLaw::scaleFactor() const{
+  return _glaw->scaleFactor();
+
+}
+
 
 void VEVPMFHDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
 	dG3DMaterialLaw::setMacroSolver(sv);
 	if (_glaw != NULL){
 		_glaw->setMacroSolver(sv);
 	}
-};
-
-//
-
-
-
+}
 
 
 //=========================================================ViscoelastiDG3DMaterialLaw===============================================================//
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 20aaacf2d10764416e83b10c9744b6b24dd45424..10b0952020f907f7e9feb023d7e2f1f39cf9ff53 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -1033,9 +1033,6 @@ class VEVPMFHDG3DMaterialLaw : public dG3DMaterialLaw{
     #endif
 };
 
-
-
-
 //======================================================ViscoelasticDG3DMaterialLaw===============================================================//
 class ViscoelasticDG3DMaterialLaw : public dG3DMaterialLaw
 {