diff --git a/NonLinearSolver/internalPoints/ipNonLocalDamage.cpp b/NonLinearSolver/internalPoints/ipNonLocalDamage.cpp
index dffbaca036351e3595bfe651a2f2ed21273ec1fc..31a5be851d4c3fb21584391919e84b405cdb8dc9 100644
--- a/NonLinearSolver/internalPoints/ipNonLocalDamage.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLocalDamage.cpp
@@ -1522,6 +1522,9 @@ void IPNonLocalDamage::restart()
   restartManager::restart(pos_Msy0);
   restartManager::restart(pos_Mhmod1);
   restartManager::restart(pos_Mhmod2);
+  restartManager::restart(pos_Malpha_DP);
+  restartManager::restart(pos_Mm_DP);
+  restartManager::restart(pos_Mnup);
   restartManager::restart(pos_Mhexp);
   restartManager::restart(pos_DamParm1);
   restartManager::restart(pos_DamParm2);
@@ -1663,6 +1666,21 @@ void IPNonLocalDamage::fillValuesToCommunicateMPI(double *arrayMPI)const
        arrayMPI[i]=_nldStatev[pos_Mhmod2];
        i++;
      }
+     if(pos_Malpha_DP !=0) 
+     {
+       arrayMPI[i]=_nldStatev[pos_Malpha_DP];
+       i++;
+     }
+     if(pos_Mm_DP !=0) 
+     {
+       arrayMPI[i]=_nldStatev[pos_Mm_DP];
+       i++;
+     }
+     if(pos_Mnup !=0) 
+     {
+       arrayMPI[i]=_nldStatev[pos_Mnup];
+       i++;
+     }
      if(pos_Mhexp !=0)  
      {
        arrayMPI[i]=_nldStatev[pos_Mhexp];
@@ -1750,6 +1768,21 @@ void IPNonLocalDamage::getValuesFromMPI(const double *arrayMPI)
        _nldStatev[pos_Mhmod2]=arrayMPI[i];
        i++;
      }
+     if(pos_Malpha_DP !=0) 
+     {
+       _nldStatev[pos_Malpha_DP]=arrayMPI[i];
+       i++;
+     }
+     if(pos_Mm_DP !=0) 
+     {
+       _nldStatev[pos_Mm_DP]=arrayMPI[i];
+       i++;
+     }
+     if(pos_Mnup !=0) 
+     {
+       _nldStatev[pos_Mnup]=arrayMPI[i];
+       i++;
+     }
      if(pos_Mhexp !=0)  
      {
        _nldStatev[pos_Mhexp]=arrayMPI[i];
diff --git a/NonLinearSolver/internalPoints/ipNonLocalDamage.h b/NonLinearSolver/internalPoints/ipNonLocalDamage.h
index d9fdca4235ac40e31fee5d05d935e621e4d2743d..1431995b86c0540e4f95282ba62204930330b64e 100644
--- a/NonLinearSolver/internalPoints/ipNonLocalDamage.h
+++ b/NonLinearSolver/internalPoints/ipNonLocalDamage.h
@@ -46,6 +46,9 @@ class IPNonLocalDamage : public IPVariableMechanics
   int pos_Msy0;
   int pos_Mhmod1;
   int pos_Mhmod2;
+  int pos_Malpha_DP;
+  int pos_Mm_DP;
+  int pos_Mnup;
   int pos_Mhexp;
   int pos_DamParm1;
   int pos_DamParm2;
@@ -187,7 +190,7 @@ class IPNonLocalDamage : public IPVariableMechanics
                             _dissipationBlocked(false),  
                             _pos_maxD(-1), _pos_inc_maxD(-1), _pos_mtx_maxD(-1), pos_vfi(0),
                             pos_euler(0), pos_aspR(0), pos_ME(0), pos_Mnu(0), pos_Msy0(0), pos_Mhmod1(0),
-                            pos_Mhmod2(0), pos_Mhexp(0), pos_DamParm1(0), pos_DamParm2(0), pos_DamParm3(0), 
+                            pos_Mhmod2(0), pos_Malpha_DP(0), pos_Mm_DP(0), pos_Mnup(0), pos_Mhexp(0), pos_DamParm1(0), pos_DamParm2(0), pos_DamParm3(0), 
                             pos_INCDamParm1(0), pos_INCDamParm2(0), pos_INCDamParm3(0), Randnum(0),
                             _irreversibleEnergy(0),  _DirreversibleEnergyDF(0.),  _DirreversibleEnergyDNonLocalVariableMatrix(0.),
                             _DirreversibleEnergyDNonLocalVariableFiber(0.),
@@ -269,6 +272,9 @@ class IPNonLocalDamage : public IPVariableMechanics
      pos_Msy0 = source.pos_Msy0;
      pos_Mhmod1 = source.pos_Mhmod1;
      pos_Mhmod2 = source.pos_Mhmod2;
+     pos_Malpha_DP = source.pos_Malpha_DP;
+     pos_Mm_DP = source.pos_Mm_DP;
+     pos_Mnup = source.pos_Mnup;
      pos_Mhexp = source.pos_Mhexp;
      pos_DamParm1 = source.pos_DamParm1;
      pos_DamParm2 = source.pos_DamParm2;
@@ -413,6 +419,9 @@ class IPNonLocalDamage : public IPVariableMechanics
      pos_Msy0 = source->pos_Msy0;
      pos_Mhmod1 = source->pos_Mhmod1;
      pos_Mhmod2 = source->pos_Mhmod2;
+     pos_Malpha_DP = source->pos_Malpha_DP;
+     pos_Mm_DP = source->pos_Mm_DP;
+     pos_Mnup = source->pos_Mnup;
      pos_Mhexp = source->pos_Mhexp;
      pos_DamParm1 = source->pos_DamParm1;
      pos_DamParm2 = source->pos_DamParm2;
@@ -594,6 +603,9 @@ class IPNonLocalDamage : public IPVariableMechanics
   virtual int getPos_Msy0() const { return pos_Msy0;};
   virtual int getPos_Mhmod1() const { return pos_Mhmod1;};
   virtual int getPos_Mhmod2() const { return pos_Mhmod2;};
+  virtual int getPos_Malpha_DP() const { return pos_Malpha_DP;};
+  virtual int getPos_Mm_DP() const { return pos_Mm_DP;};
+  virtual int getPos_Mnup() const { return pos_Mnup;};
   virtual int getPos_Mhexp() const { return pos_Mhexp;};
   virtual int getPos_DamParm1() const { return pos_DamParm1;};
   virtual int getPos_DamParm2() const { return pos_DamParm2;};
@@ -611,6 +623,9 @@ class IPNonLocalDamage : public IPVariableMechanics
   virtual void setPos_Msy0( int p ) {  pos_Msy0=p;};
   virtual void setPos_Mhmod1( int p ) {  pos_Mhmod1=p;};
   virtual void setPos_Mhmod2( int p ) {  pos_Mhmod2=p;};
+  virtual void setPos_Malpha_DP( int p ) {  pos_Malpha_DP=p;};
+  virtual void setPos_Mm_DP( int p ) {  pos_Mm_DP=p;};
+  virtual void setPos_Mnup( int p ) {  pos_Mnup=p;};
   virtual void setPos_Mhexp( int p ) {  pos_Mhexp=p;};
   virtual void setPos_DamParm1( int p ) {  pos_DamParm1=p;};
   virtual void setPos_DamParm2( int p ) {  pos_DamParm2=p;};
diff --git a/NonLinearSolver/materialLaw/damage.h b/NonLinearSolver/materialLaw/damage.h
index ccd86915800a9a2955a2943303e040df081c7490..93f1310d7d0d367b52b02bbfc183c74689d1188e 100755
--- a/NonLinearSolver/materialLaw/damage.h
+++ b/NonLinearSolver/materialLaw/damage.h
@@ -14,11 +14,13 @@ class Damage{
 		int ndamprops;
 		int nsdv;
                 int pos_maxD;
+                bool localDamage;
 		
 	public:
 		
 		Damage(double* props, int iddam){
 			dammodel = (int) props[iddam];
+			localDamage = false;
 		}
 		virtual ~Damage(){ ; }
 		inline int get_dammodel(){  return dammodel; }
@@ -43,7 +45,10 @@ class Damage{
                 virtual double getMaxD(double *statev, int pos_dam) const {return statev[pos_dam+get_pos_maxD()];}
                 virtual void setMaxD(double *statev, int pos_dam, double Dmax) { statev[pos_dam+get_pos_maxD()]=Dmax;}
                 virtual int getNbNonLocalVariables() const {return 1;}
-
+#ifdef NONLOCALGMSH
+                virtual void setLocalDamage() {localDamage=true;}
+                virtual void setNonLocalDamage() {localDamage=false;}
+#endif
 };
 
 //Linear elastic material
diff --git a/NonLinearSolver/materialLaw/material.cpp b/NonLinearSolver/materialLaw/material.cpp
index b27c085ba012417893d36d11642ccda13d569378..2ba2fd8a4d5d9f1314b412180415dbc11efc9425 100644
--- a/NonLinearSolver/materialLaw/material.cpp
+++ b/NonLinearSolver/materialLaw/material.cpp
@@ -1416,7 +1416,7 @@ void EVP_Material::unload_step(double* dstrn,double* strs_n, double* statev_n, i
 
 //********************************************************************************************
 // COMPOSITE - VOIGT  
-//PROPS: MODEL, IDDAM, Nph, IDMAT phase 1, IDTOP phase 1, ... , 
+//PROPS: MODEL, IDDAM, NGrain, IDMAT phase 1, IDTOP phase 1, ... , 
 //	 	...	IDMAT phase N, IDTOP phase N
 //  (material properties of each phase), (volume fraction of each phase)
 //STATEV: strn, strs, dstrn, (statev of each phase)   
@@ -1430,13 +1430,13 @@ VT_Material::VT_Material(double* props, int idmat):Material(props, idmat){
 	nsdv=18;
 	iddam = (int) props[idmat+1];
         idclength = (int)props[idmat+2];
-	Nph = (int) props[idmat+3];
-	mat = (Material**) malloc(Nph*sizeof(Material*));	
-	mallocvector(&vf,Nph);
+        NGrain = (int) props[idmat+3];
+	mat = (Material**) malloc(NGrain*sizeof(Material*));	
+	mallocvector(&vf,NGrain);
         
 	k=idmat+4;
 
-	for(i=0;i<Nph;i++){
+	for(i=0;i<NGrain;i++){
 		idmat_i = (int)props[k++];
 		idtop_i = (int)props[k++];
 		
@@ -1516,7 +1516,7 @@ double VT_Material::getPlasticEnergy (double* statev)
 VT_Material::~VT_Material(){
 
 	int i;
-	for (i=0;i<Nph;i++){
+	for (i=0;i<NGrain;i++){
 		if(mat[i]){
 			delete mat[i];
 		}
@@ -1535,9 +1535,9 @@ void VT_Material::print(){
 	printf("Voigt\n");
 	if(dam){ dam->print();}
 
-	printf("Number of phases: %d\n", Nph);
+	printf("Number of phases: %d\n", NGrain);
 	
-	for(i=0;i<Nph;i++){
+	for(i=0;i<NGrain;i++){
 		printf("phase %d, volume fraction: %lf\n",i+1,vf[i]);
 		mat[i]->print();
 	}
@@ -1644,7 +1644,7 @@ double LAM2Ply_Material::getPlasticEnergy (double* statev)
 LAM2Ply_Material::~LAM2Ply_Material(){
 
 	int i;
-	for (i=0;i<Nph;i++){
+	for (i=0;i<NPly;i++){
 		if(mat[i]){
 			delete mat[i];
 		}
@@ -1663,9 +1663,9 @@ void LAM2Ply_Material::print(){
 	printf("LAM2Ply\n");
 	if(dam){ dam->print();}
 
-	printf("Number of phases: %d\n", Nph);
+	printf("Number of phases: %d\n", NPly);
 	
-	for(i=0;i<Nph;i++){
+	for(i=0;i<NPly;i++){
 		printf("phase %d, volume fraction: %lf\n",i+1,vf[i]);
 		mat[i]->print();
 	}
@@ -1845,6 +1845,23 @@ void MT_Material::setEuler(double e1, double e2, double e3)
 }
 #endif
 
+int MT_Material::get_pos_Fd_bar()const
+{
+   if(icl_mat->get_pos_Fd_bar()!= 0) return idsdv_i + icl_mat->get_pos_Fd_bar();
+   else return 0;
+}
+int MT_Material::get_pos_dFd_bar()const
+{
+   if(icl_mat->get_pos_dFd_bar()!= 0) return idsdv_i + icl_mat->get_pos_dFd_bar();
+   else return 0;
+}
+int MT_Material::get_pos_locFd()const
+{
+   if(icl_mat->get_pos_locFd()!= 0) return idsdv_i + icl_mat->get_pos_locFd();
+   else return 0;
+}
+
+
 //destructor
 MT_Material::~MT_Material(){
 
@@ -2259,6 +2276,27 @@ void MTSecF_Material::unload_step(double* dstrn, double* strs_n, double* statev_
 
 }
 
+#ifdef NONLOCALGMSH
+
+void MTSecF_Material:: get_elOp_mtx(double* statev, double** Cel){
+       mtx_mat->get_elOp(&statev[idsdv_m],Cel);
+}
+
+void MTSecF_Material::get_elOp_icl(double* statev, double** Cel){
+       icl_mat->get_elOp(&statev[idsdv_i],Cel);
+}
+
+void MTSecF_Material:: get_elOp_mtx(double** Cel){
+       mtx_mat->get_elOp(Cel);
+}
+
+void MTSecF_Material::get_elOp_icl(double** Cel){
+       icl_mat->get_elOp(Cel);
+}
+
+#endif
+
+
 
 //**************************************************************************************************************************************************
 // COMPOSITE - MORI-TANAKA   
@@ -2835,6 +2873,21 @@ int MTSecF_Stoch_Material::get_pos_hmod2()const
    if(mtx_mat->get_pos_hmod2()!= 0) return idsdv_m + mtx_mat->get_pos_hmod2();
    else return 0;
 }
+int MTSecF_Stoch_Material::get_pos_alpha_DP()const
+{
+   if(mtx_mat->get_pos_alpha_DP()!= 0) return idsdv_m + mtx_mat->get_pos_alpha_DP();
+   else return 0;
+}
+int MTSecF_Stoch_Material::get_pos_m_DP()const
+{
+   if(mtx_mat->get_pos_m_DP()!= 0) return idsdv_m + mtx_mat->get_pos_m_DP();
+   else return 0;
+}
+int MTSecF_Stoch_Material::get_pos_nup()const
+{
+   if(mtx_mat->get_pos_nup()!= 0) return idsdv_m + mtx_mat->get_pos_nup();
+   else return 0;
+}
 int MTSecF_Stoch_Material::get_pos_hp0()const 
 {
    if(mtx_mat->get_pos_hp0()!= 0) return idsdv_m + mtx_mat->get_pos_hp0();
diff --git a/NonLinearSolver/materialLaw/material.h b/NonLinearSolver/materialLaw/material.h
index bd8555c79b62663b03b94d24f48a9ee73c024a14..19d059f43cf5ff0b788fe4c9a8f485fd30f6ba22 100644
--- a/NonLinearSolver/materialLaw/material.h
+++ b/NonLinearSolver/materialLaw/material.h
@@ -6,6 +6,7 @@
 #include "clength.h"
 #include "lccfunctions.h"
 #include "elasticlcc.h"
+#include <vector>
 
 
 #ifdef NONLOCALGMSH
@@ -22,6 +23,8 @@ class Material{
 		Damage* dam;
 		int pos_dam;	//position of damage state variables in statev vector
 
+		bool localDamage;
+
                 int pos_p;
 
                 Clength* clength;
@@ -34,6 +37,7 @@ class Material{
 		Material(double* props, int idmat){
 			constmodel = (int) props[idmat];
                         evaluateStiffness = true;
+                        localDamage = false;
 		}
 		virtual ~Material(){ ; }
 		inline int get_constmodel(){  return constmodel; }
@@ -117,6 +121,9 @@ double E, double nu, int htype, double sy0, double hexp, double hmod1, double hm
 		virtual int get_pos_sy0() const {  return 0; }  
 		virtual int get_pos_hmod1() const {  return 0; } 
 		virtual int get_pos_hmod2() const {  return 0; } 
+		virtual int get_pos_alpha_DP() const {  return 0; } 
+		virtual int get_pos_m_DP() const {  return 0; } 
+		virtual int get_pos_nup() const {  return 0; } 
 		virtual int get_pos_hp0() const {  return 0; } 
 		virtual int get_pos_hexp() const {  return 0; } 
 		virtual int get_pos_DamParm1() const {  return 0; } 
@@ -137,6 +144,13 @@ double E, double nu, int htype, double sy0, double hexp, double hmod1, double hm
                 virtual double get_mtx_MaxD(double *statev) const {return -1.; }
                 virtual void set_mtx_MaxD(double *statev, double Dmax) {}
                 
+                virtual std::vector<double> &getEuler(std::vector<double> &_euler) {return _euler;}
+                virtual int getNPatch() const {return -1;}
+                virtual int get_mat_NPatch(int matNb) const {return -1;}
+		virtual int get_mat_constmodel(int matNb) const {return -1;}
+		virtual int get_submat_constmodel(int matNb, int submatNb) const {return -1;}
+
+
 #ifdef NONLOCALGMSH
 		virtual int get_pos_currentplasticstrainfornonlocal() const = 0;
 		virtual int get_pos_effectiveplasticstrainfornonlocal() const = 0;
@@ -160,54 +174,103 @@ double E, double nu, int htype, double sy0, double hexp, double hmod1, double hm
                 virtual double get_vfI();
                 
                 virtual void setEuler(double e1, double e2, double e3);
-#endif
-		
-    		// some new function
-		virtual int getNPatch() const {return -1;}
-		virtual int get_mat_NPatch(int matNb) const {return -1;}
-		virtual int get_mat_constmodel(int matNb) const {return -1;}
-		virtual int get_submat_constmodel(int matNb, int submatNb) const {return -1;}
-		virtual void populateMaxD(double *statev, double Dmax_Inc, double Dmax_Mtx) const{};
+                
+		virtual void setLocalDamage() {
+			localDamage=true;
+			if(dam) dam->setLocalDamage();
+		}
 
+		virtual void setNonLocalDamage() {
+			localDamage=false;
+			if(dam) dam->setNonLocalDamage();
+		}
+		
+		virtual int get_pos_mat_inc_maxD(int matNb) const {return -1;}
+		virtual int get_pos_mat_mtx_maxD(int matNb) const {return -1;}
+                virtual void populateMaxD(double *statev, double Dmax_Inc, double Dmax_Mtx) const {}
+
+
+		// for VT or LAM_2PLY case (VM, LVM or VLM)
+                virtual int get_pos_mat_strain(int matNb) const {return -1;}
+                virtual int get_pos_mat_stress(int matNb) const {return -1;}
+                virtual int get_pos_mat_Dam(int matNb) const {return -1;}
+                virtual int get_pos_mat_mtx_strain(int matNb) const {return -1;}
+                virtual int get_pos_mat_inc_strain(int matNb) const {return -1;}
+                virtual int get_pos_mat_mtx_stress(int matNb) const {return -1;}
+                virtual int get_pos_mat_inc_stress(int matNb) const {return -1;}
+                virtual int get_pos_mat_mtx_Dam(int matNb) const {return -1;}
+                virtual int get_pos_mat_inc_Dam(int matNb) const {return -1;}
+                
+                // for VT case (VM or VLM)
+                // for VLM case
+		// laminate euler angles
 		virtual std::vector<double> &getEulerVTLam(std::vector<double> &euler_lam, int grainNb) {return euler_lam;}
-
+		// for each VT grain
+		// strain
 		virtual int getPosStrnVTM0(int grainNb) const {return -1;}
 		virtual int getPosStrnVTMT(int grainNb) const {return -1;}
 		virtual int getPosStrnVTLam(int grainNb) const {return -1;}
+		// stress
 		virtual int getPosStrsVTM0(int grainNb) const {return -1.;}
 		virtual int getPosStrsVTMT(int grainNb) const {return -1;}
 		virtual int getPosStrsVTLam(int grainNb) const {return -1;}
+		// damage
                 virtual int getPosDamVTM0(int grainNb) const {return -1;}
+                // for each laminate ply in VT grains
+		// strain
 		virtual int getPosStrnVTLamM0(int grainNb) const {return -1;}
 		virtual int getPosStrnVTLamMT(int grainNb) const {return -1;}
+		// stress
 		virtual int getPosStrsVTLamM0(int grainNb) const {return -1;}
 		virtual int getPosStrsVTLamMT(int grainNb) const {return -1;}
+		// damage
 		virtual int getPosDamVTLamM0(int grainNb) const {return -1;}
+		// for each material phase
+		// strain
 		virtual int getPosStrnVTMTMtx(int grainNb) const {return -1;}
 		virtual int getPosStrnVTMTInc(int grainNb) const {return -1;}
 		virtual int getPosStrnVTLamMTMtx(int grainNb) const {return -1;}
 		virtual int getPosStrnVTLamMTInc(int grainNb) const {return -1;}
+		// stress
 		virtual int getPosStrsVTMTMtx(int grainNb) const {return -1;}
 		virtual int getPosStrsVTMTInc(int grainNb) const {return -1;}
 		virtual int getPosStrsVTLamMTMtx(int grainNb) const {return -1;}
 		virtual int getPosStrsVTLamMTInc(int grainNb) const {return -1;}
+		// damage
 		virtual int getPosDamVTMTMtx(int grainNb) const {return -1;}
 		virtual int getPosDamVTMTInc(int grainNb) const {return -1;}
 		virtual int getPosDamVTLamMTMtx(int grainNb) const {return -1;}
 		virtual int getPosDamVTLamMTInc(int grainNb) const {return -1;}
+                
+                
+                // for LAM_2PLY case (LVM)
+                // for each laminate ply
+                // strain
 		virtual int getPosStrnLamM0(int laminateNb) const {return -1;}
 		virtual int getPosStrnLamVT(int laminateNb) const {return -1;}
+		// stress
 		virtual int getPosStrsLamM0(int laminateNb) const {return -1;}
 		virtual int getPosStrsLamVT(int laminateNb) const {return -1;}
+		// damage
 		virtual int getPosDamLamM0(int laminateNb) const {return -1;}
+		// for each MT grain in VT ply
+		// strain
 		virtual int getPosStrnLamVTMT(int grainNb) const {return -1;}
+		// stress
 		virtual int getPosStrsLamVTMT(int grainNb) const {return -1;}
+		// for each material phase
+		// strain
 		virtual int getPosStrnLamVTMTMtx(int grainNb) const {return -1;}
 		virtual int getPosStrnLamVTMTInc(int grainNb) const {return -1;}
+		// stress
 		virtual int getPosStrsLamVTMTMtx(int grainNb) const {return -1;}
 		virtual int getPosStrsLamVTMTInc(int grainNb) const {return -1;}
+		// damage
 		virtual int getPosDamLamVTMTMtx(int grainNb) const {return -1;}
 		virtual int getPosDamLamVTMTInc(int grainNb) const {return -1;}
+#endif
+
+
 };
 
 //Linear elastic material
@@ -318,10 +381,10 @@ double* dnu, double &dnudp_bar, double alpha, double* dpdE, double* dstrsdp_bar,
 class EP_Stoch_Material : public Material{
 
 	private:
-		double _E,_nu,_sy0;
+		double _E,_nu,_sy0, _alpha_DP, _m_DP, _nup;
 		int htype;
 		double _hmod1, _hmod2, _hp0, _hexp;
-                int pos_E, pos_nu, pos_sy0, alpha_DP, m_DP, nup;
+                int pos_E, pos_nu, pos_sy0, pos_alpha_DP, pos_m_DP, pos_nup;
                 int pos_hmod1, pos_hmod2, pos_hp0, pos_hexp;
                 int pos_DamParm1, pos_DamParm2, pos_DamParm3;
 		int pos_pstrn, pos_dp, pos_twomu, pos_eqstrs;
@@ -354,6 +417,10 @@ class EP_Stoch_Material : public Material{
 		virtual int get_pos_E() const {  return pos_E; } 
 		virtual int get_pos_nu() const {  return pos_nu; } 
 		virtual int get_pos_sy0() const {  return pos_sy0; } 
+
+		virtual int get_pos_alpha_DP() const {  return pos_alpha_DP; } 
+		virtual int get_pos_m_DP () const {  return pos_m_DP ; } 
+		virtual int get_pos_nup() const {  return pos_nup; } 
  
 		virtual int get_pos_hmod1() const {  return pos_hmod1; } 
 		virtual int get_pos_hmod2() const {  return pos_hmod2; } 
@@ -434,10 +501,40 @@ double* dnu, double &dnudp_bar, double alpha, double* dpdE, double* dstrsdp_bar,
 class VT_Material : public Material{
 
 	private:
-		int Nph;					//number of phases in the composite
+		int NGrain;					//number of grains
 		Material** mat;		//array of pointers on Materials
-		double* vf;				//array of volume fraction of each phase
-			
+		double* vf;				//array of volume fraction of each grain
+				
+#ifdef NONLOCALGMSH
+		// for VT case (VM or VLM)
+		// for VLM case
+		// laminate euler angles
+		std::vector<double> euler_vt_lam;
+		// for each VT grain
+		// strain
+		std::vector<int> pos_strn_vt_m0, pos_strn_vt_mt, pos_strn_vt_lam;
+		// stress
+		std::vector<int> pos_strs_vt_m0, pos_strs_vt_mt, pos_strs_vt_lam;
+		// damage
+		std::vector<int> pos_dam_vt_m0;
+		// for each laminate ply in VT grains
+		// strain
+		std::vector<int> pos_strn_vt_lam_m0, pos_strn_vt_lam_mt;
+		// stress
+		std::vector<int> pos_strs_vt_lam_m0, pos_strs_vt_lam_mt;
+		// damage
+		std::vector<int> pos_dam_vt_lam_m0;
+		// for each material phase
+		// strain
+		std::vector<int> pos_strn_vt_mt_mtx, pos_strn_vt_mt_inc;
+		std::vector<int> pos_strn_vt_lam_mt_mtx, pos_strn_vt_lam_mt_inc;
+		// stress
+		std::vector<int> pos_strs_vt_mt_mtx, pos_strs_vt_mt_inc;
+		std::vector<int> pos_strs_vt_lam_mt_mtx, pos_strs_vt_lam_mt_inc;
+		// damage
+		std::vector<int> pos_dam_vt_mt_mtx, pos_dam_vt_mt_inc;
+		std::vector<int> pos_dam_vt_lam_mt_mtx, pos_dam_vt_lam_mt_inc;
+#endif
 
 	public:
 		
@@ -450,12 +547,257 @@ class VT_Material : public Material{
 		virtual void get_elOp(double** Cel);
 		virtual void get_elOD(double* statev, double** Cel);
                 virtual void unload_step(double* dstrn, double* strs_n, double* statev, int kinc, int kstep);
+
+                virtual int getNPatch() const {return NGrain;}
+                virtual int get_mat_NPatch(int grainNb) const {return mat[grainNb]->getNPatch();}                
+                virtual int get_mat_constmodel(int grainNb) const {return mat[grainNb]->get_constmodel();}
+                virtual int get_submat_constmodel(int grainNb, int laminateNb) const {return mat[grainNb]->get_mat_constmodel(laminateNb);}
+                             
 #ifdef NONLOCALGMSH
 		virtual int get_pos_currentplasticstrainfornonlocal() const;
 		virtual int get_pos_effectiveplasticstrainfornonlocal() const;
 		virtual int get_pos_damagefornonlocal() const;
                 virtual double getElasticEnergy (double* statev) ;
 		virtual double getPlasticEnergy (double* statev) ;
+		
+		/*virtual void setLocalDamage()
+		{
+                    Material::setLocalDamage();
+                    for(i=0;i<NGrain;i++){
+		        mat[i]->setLocalDamage();
+                    }
+                }*/
+                
+                virtual int get_pos_mat_inc_maxD(int grainNb) const {             
+                	int idsdv = 18;  // position in statev where the statev of the first phase start
+                	for(int i=0;i<grainNb;i++) {
+                		idsdv += mat[i]->get_nsdv();
+                	}
+                	
+			switch(mat[grainNb]->get_constmodel()){		                
+				case MTSecF:  // MFH Fisrst order secant
+					if(mat[grainNb]->get_pos_inc_maxD()>-1)  // MT-Inc phase
+						return (idsdv + mat[grainNb]->get_pos_inc_maxD());
+					return -1;
+					break;
+					
+				default:
+					return -1;
+					break;
+			}        
+                }
+                
+                virtual int get_pos_mat_mtx_maxD(int grainNb) const {             
+                	int idsdv = 18;  // position in statev where the statev of the first phase start
+                	for(int i=0;i<grainNb;i++) {
+                		idsdv += mat[i]->get_nsdv();
+                	}
+                	
+			switch(mat[grainNb]->get_constmodel()){		                
+				case EP:  // pure matrix
+		                	if(mat[grainNb]->get_pos_maxD()>-1)  // EP case
+		                		return (idsdv + mat[grainNb]->get_pos_maxD());
+		                	return -1;
+		                	break;
+				
+				case MTSecF:  // MFH Fisrst order secant
+					if(mat[grainNb]->get_pos_mtx_maxD()>-1)  // MT-Inc phase
+						return (idsdv + mat[grainNb]->get_pos_mtx_maxD());
+					return -1;
+					break;
+					
+				default:
+					return -1;
+					break;
+			}        
+                }
+            
+		// here we only work for the initial value, will not work to block damage evolution in the grain on the fly
+                virtual void populateMaxD(double *statev, double Dmax_Inc, double Dmax_Mtx) const {                    
+                	int idsdv = 18;  // position in statev where the statev of the first phase start
+                	for(int i=0;i<NGrain;i++) {
+                        	switch(mat[i]->get_constmodel()){
+                        		case EP:  // pure matrix
+                        			if (mat[i]->get_pos_maxD()>-1)
+		                		{
+		                    			statev[idsdv + mat[i]->get_pos_maxD()] = Dmax_Mtx;
+		                		}
+		                		break;
+		                		
+		                	case MTSecF:  // MFH First order secant
+		                		if(mat[i]->get_pos_inc_maxD()>-1)  // MT-Inc phase
+		                		{
+		                    			statev[idsdv + mat[i]->get_pos_inc_maxD()] = Dmax_Inc;
+		                		}
+		                		if(mat[i]->get_pos_mtx_maxD()>-1)  // MT-Mtx phase
+		                		{
+		                    			statev[idsdv + mat[i]->get_pos_mtx_maxD()] = Dmax_Mtx;
+		                		}
+		                		break;
+		                		
+		                	case LAM_2PLY:
+		                		for(int j=0;j<mat[i]->getNPatch();j++)
+		                		{
+		                			if(mat[i]->get_pos_mat_inc_maxD(j)>-1)  // LAM_2PLY-Inc phase
+		                			{
+		                    				statev[idsdv + mat[i]->get_pos_mat_inc_maxD(j)] = Dmax_Inc;
+		                			}
+		                			if(mat[i]->get_pos_mat_mtx_maxD(j)>-1)  // LAM_2PLY-Mtx phase
+		                			{
+		                    				statev[idsdv + mat[i]->get_pos_mat_mtx_maxD(j)] = Dmax_Mtx;
+		                			}
+		                		}
+		                		break;
+		                }
+		                
+		                idsdv += mat[i]->get_nsdv();
+		        }
+                }
+                
+                // for VT case (VM or VLM)
+                // for each VT grain
+                // strain
+                virtual int get_pos_mat_strain (int grainNb) const
+                {
+                    int idsdv = 18;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<grainNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[grainNb]->get_pos_strn();
+		}
+		// stress
+		virtual int get_pos_mat_stress (int grainNb) const
+                {
+                    int idsdv = 18;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<grainNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[grainNb]->get_pos_strs();
+		}
+		// damage
+		virtual int get_pos_mat_Dam (int grainNb) const
+                {
+                    if(mat[grainNb]->get_pos_damagefornonlocal()>-1){
+                        int idsdv = 18;  // position in statev where the statev of the first phase start
+                        for(int i=0;i<grainNb;i++){
+                            idsdv += mat[i]->get_nsdv();
+                        }
+		        return idsdv + mat[grainNb]->get_pos_damagefornonlocal();
+		    }
+		    else{
+		        return -1;
+		    }
+		}
+		// for each material phase
+		// strain
+		virtual int get_pos_mat_mtx_strain (int grainNb) const
+		{
+		    int idsdv = 18;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<grainNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[grainNb]->get_pos_mtx_strain();
+		}
+		virtual int get_pos_mat_inc_strain (int grainNb) const
+		{
+		    int idsdv = 18;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<grainNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    } 
+		    return idsdv + mat[grainNb]->get_pos_inc_strain();
+		}
+		// stress
+		virtual int get_pos_mat_mtx_stress (int grainNb) const
+		{
+		    int idsdv = 18;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<grainNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[grainNb]->get_pos_mtx_stress();
+		}
+		virtual int get_pos_mat_inc_stress (int grainNb) const
+		{
+		    int idsdv = 18;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<grainNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[grainNb]->get_pos_inc_stress();
+		}
+		// damage
+		virtual int get_pos_mat_mtx_Dam (int grainNb) const
+		{
+		    if(mat[grainNb]->get_pos_mtx_Dam()>-1){
+                        int idsdv = 18;  // position in statev where the statev of the first phase start
+                        for(int i=0;i<grainNb;i++){
+                            idsdv += mat[i]->get_nsdv();
+                        }
+		        return idsdv + mat[grainNb]->get_pos_mtx_Dam();
+		    }
+		    else{
+		        return -1;
+		    }
+		}
+		virtual int get_pos_mat_inc_Dam (int grainNb) const
+		{
+		    if(mat[grainNb]->get_pos_inc_Dam()>-1){
+                        int idsdv = 18;  // position in statev where the statev of the first phase start
+                        for(int i=0;i<grainNb;i++){
+                            idsdv += mat[i]->get_nsdv();
+                        }
+		        return idsdv + mat[grainNb]->get_pos_inc_Dam();
+		    }
+		    else{
+		        return -1;
+		    }
+		}
+
+		// for VT case (VM or VLM)
+		// for VLM case
+		// laminate euler angles
+		virtual std::vector<double> &getEulerVTLam(std::vector<double> &euler_lam, int grainNb)
+		{
+			euler_lam.clear();
+			euler_lam.emplace_back(euler_vt_lam[0+grainNb*3]);
+			euler_lam.emplace_back(euler_vt_lam[1+grainNb*3]);
+			euler_lam.emplace_back(euler_vt_lam[2+grainNb*3]);
+			return euler_lam;
+		}
+		// for each VT grain
+		// strain
+		virtual int getPosStrnVTM0(int grainNb) const {return (int)pos_strn_vt_m0[grainNb];}
+		virtual int getPosStrnVTMT(int grainNb) const {return (int)pos_strn_vt_mt[grainNb];}
+		virtual int getPosStrnVTLam(int grainNb) const {return (int)pos_strn_vt_lam[grainNb];}
+		// stress
+		virtual int getPosStrsVTM0(int grainNb) const {return (int)pos_strs_vt_m0[grainNb];}
+		virtual int getPosStrsVTMT(int grainNb) const {return (int)pos_strs_vt_mt[grainNb];}
+		virtual int getPosStrsVTLam(int grainNb) const {return (int)pos_strs_vt_lam[grainNb];}
+		// damage
+		virtual int getPosDamVTM0(int grainNb) const {return (int)pos_dam_vt_m0[grainNb];}
+		// for each laminate ply in VT grains
+		// strain
+		virtual int getPosStrnVTLamM0(int grainNb) const {return (int)pos_strn_vt_lam_m0[grainNb];}
+		virtual int getPosStrnVTLamMT(int grainNb) const {return (int)pos_strn_vt_lam_mt[grainNb];}
+		// stress
+		virtual int getPosStrsVTLamM0(int grainNb) const {return (int)pos_strs_vt_lam_m0[grainNb];}
+		virtual int getPosStrsVTLamMT(int grainNb) const {return (int)pos_strs_vt_lam_mt[grainNb];}
+		// damage
+		virtual int getPosDamVTLamM0(int grainNb) const {return (int)pos_dam_vt_lam_m0[grainNb];}
+                // for each material phase
+                // strain
+                virtual int getPosStrnVTMTMtx(int grainNb) const {return (int)pos_strn_vt_mt_mtx[grainNb];}
+                virtual int getPosStrnVTMTInc(int grainNb) const {return (int)pos_strn_vt_mt_inc[grainNb];}
+                virtual int getPosStrnVTLamMTMtx(int grainNb) const {return (int)pos_strn_vt_lam_mt_mtx[grainNb];}
+                virtual int getPosStrnVTLamMTInc(int grainNb) const {return (int)pos_strn_vt_lam_mt_inc[grainNb];}
+                // stress
+                virtual int getPosStrsVTMTMtx(int grainNb) const {return (int)pos_strs_vt_mt_mtx[grainNb];}
+                virtual int getPosStrsVTMTInc(int grainNb) const {return (int)pos_strs_vt_mt_inc[grainNb];}
+                virtual int getPosStrsVTLamMTMtx(int grainNb) const {return (int)pos_strs_vt_lam_mt_mtx[grainNb];}
+                virtual int getPosStrsVTLamMTInc(int grainNb) const {return (int)pos_strs_vt_lam_mt_inc[grainNb];}
+                // damage
+                virtual int getPosDamVTMTMtx(int grainNb) const {return (int)pos_dam_vt_mt_mtx[grainNb];}
+                virtual int getPosDamVTMTInc(int grainNb) const {return (int)pos_dam_vt_mt_inc[grainNb];}
+                virtual int getPosDamVTLamMTMtx(int grainNb) const {return (int)pos_dam_vt_lam_mt_mtx[grainNb];}
+                virtual int getPosDamVTLamMTInc(int grainNb) const {return (int)pos_dam_vt_lam_mt_inc[grainNb];}
 #endif
 
 };  
@@ -464,12 +806,35 @@ class VT_Material : public Material{
 class LAM2Ply_Material : public Material{
 
 	private:
-		int Nph;					//number of phases in the composite
+		int NPly=2;					//number of plies
 		Material** mat;		//array of pointers on Materials
-		double* vf;				//array of volume fraction of each phase
+		double* vf;				//array of volume fraction of each ply
                 int idsdv_A, idsdv_B;
                 double* euler;
-                int pos_C;  
+                int pos_C;
+                
+#ifdef NONLOCALGMSH
+		// for LAM_2PLY case (LVM)
+		// for each laminate ply
+		// strain
+		std::vector<int> pos_strn_lam_m0, pos_strn_lam_vt;
+		// stress
+		std::vector<int> pos_strs_lam_m0, pos_strs_lam_vt;
+		// damage
+		std::vector<int> pos_dam_lam_m0;
+		// for each MT grain in VT ply
+		// strain
+		std::vector<int> pos_strn_lam_vt_mt;
+		// stress
+		std::vector<int> pos_strs_lam_vt_mt;
+		// for each material phase
+		// strain
+		std::vector<int> pos_strn_lam_vt_mt_mtx, pos_strn_lam_vt_mt_inc;
+		// stress
+		std::vector<int> pos_strs_lam_vt_mt_mtx, pos_strs_lam_vt_mt_inc;
+		// damage
+		std::vector<int> pos_dam_lam_vt_mt_mtx, pos_dam_lam_vt_mt_inc;
+#endif                
 	
 	public:
 		
@@ -482,13 +847,292 @@ class LAM2Ply_Material : public Material{
 		virtual void get_elOp(double** Cel);
 		virtual void get_elOD(double* statev, double** Cel);
                 virtual void unload_step(double* dstrn, double* strs_n, double* statev, int kinc, int kstep);
-                virtual int get_pos_C() const {	return pos_C; } 
+                virtual int get_pos_C() const {	return pos_C; }
+
+		virtual std::vector<double> &getEuler(std::vector<double> &_euler)
+		{
+			_euler.clear();
+			_euler.emplace_back(euler[0]);
+			_euler.emplace_back(euler[1]);
+			_euler.emplace_back(euler[2]);
+			return _euler;
+		}
+
+		virtual int getNPatch() const {return NPly;}
+		virtual int get_mat_NPatch(int laminateNb) const {return mat[laminateNb]->getNPatch();}
+		virtual int get_mat_constmodel(int laminateNb) const {return mat[laminateNb]->get_constmodel();}
+		virtual int get_submat_constmodel(int laminateNb, int grainNb) const {return mat[laminateNb]->get_mat_constmodel(grainNb);}
+
 #ifdef NONLOCALGMSH
 		virtual int get_pos_currentplasticstrainfornonlocal() const;
 		virtual int get_pos_effectiveplasticstrainfornonlocal() const;
 		virtual int get_pos_damagefornonlocal() const;
                 virtual double getElasticEnergy (double* statev) ;
 		virtual double getPlasticEnergy (double* statev) ;
+
+		/*virtual void setLocalDamage()
+		{
+                    Material::setLocalDamage();
+                    for(i=0;i<NPly;i++){
+		        mat[i]->setLocalDamage();
+                    }
+                }*/
+                
+                virtual int get_pos_mat_inc_maxD(int laminateNb) const {             
+                	int idsdv = 39;  // position in statev where the statev of the first phase start
+                	for(int i=0;i<laminateNb;i++) {
+                		idsdv += mat[i]->get_nsdv();
+                	}
+                	
+			switch(mat[laminateNb]->get_constmodel()){		                
+				case MTSecF:  // MFH Fisrst order secant
+					if(mat[laminateNb]->get_pos_inc_maxD()>-1)  // MT-Inc phase
+						return (idsdv + mat[laminateNb]->get_pos_inc_maxD());
+					return -1;
+					break;
+				
+				default:
+					return -1;
+					break;
+				
+			}        
+                }
+                
+                virtual int get_pos_mat_mtx_maxD(int laminateNb) const {             
+                	int idsdv = 39;  // position in statev where the statev of the first phase start
+                	for(int i=0;i<laminateNb;i++) {
+                		idsdv += mat[i]->get_nsdv();
+                	}
+                	
+			switch(mat[laminateNb]->get_constmodel()){		                
+				case EP:  // pure matrix
+		                	if(mat[laminateNb]->get_pos_maxD()>-1)  // EP case
+						return (idsdv + mat[laminateNb]->get_pos_maxD());
+		                	return -1;
+		                	break;
+				
+				case MTSecF:  // MFH Fisrst order secant
+					if(mat[laminateNb]->get_pos_mtx_maxD()>-1)  // MT-Inc phase
+						return (idsdv + mat[laminateNb]->get_pos_mtx_maxD());
+					return -1;
+					break;
+					
+				default:
+					return -1;
+					break;
+			}        
+                }
+                
+                // here we only work for the initial value, will not work to block damage evolution in the grain on the fly
+                virtual void populateMaxD(double *statev, double Dmax_Inc, double Dmax_Mtx) const {                    
+                	int idsdv = 39;  // position in statev where the statev of the first phase start
+                	for(int i=0;i<NPly;i++) {
+                        	switch(mat[i]->get_constmodel()){
+                        		case EP:  // pure matrix
+                        			if (mat[i]->get_pos_maxD()>-1)
+		                		{
+		                    			statev[idsdv + mat[i]->get_pos_maxD()] = Dmax_Mtx;
+		                		}
+		                		break;
+		                		
+		                	case MTSecF:  // MFH First order secant
+		                		if(mat[i]->get_pos_inc_maxD()>-1)  // MT-Inc phase
+		                		{
+		                    			statev[idsdv + mat[i]->get_pos_inc_maxD()] = Dmax_Inc;
+		                		}
+		                		if(mat[i]->get_pos_mtx_maxD()>-1)  // MT-Mtx phase
+		                		{
+		                    			statev[idsdv + mat[i]->get_pos_mtx_maxD()] = Dmax_Mtx;
+		                		}
+		                		break;
+		                		
+		                	case VT:
+		                		for(int j=0;j<mat[i]->getNPatch();j++)
+		                		{
+		                			if(mat[i]->get_pos_mat_inc_maxD(j)>-1)  // LAM_2PLY-Inc phase
+		                			{
+		                    				statev[idsdv + mat[i]->get_pos_mat_inc_maxD(j)] = Dmax_Inc;
+		                			}
+		                			if(mat[i]->get_pos_mat_mtx_maxD(j)>-1)  // LAM_2PLY-Mtx phase
+		                			{
+		                    				statev[idsdv + mat[i]->get_pos_mat_mtx_maxD(j)] = Dmax_Mtx;
+		                			}
+		                		}
+		                		break;
+		                }
+		                
+		                idsdv += mat[i]->get_nsdv();
+		        }
+                }
+
+		/*virtual std::vector<int> &get_allPos_mtx_maxD() const {
+                    std::vector<int> allPos_mtx_maxD;
+                    
+                    int idsdv = 39;  // position in statev where the statev of the first phase start                    
+		    for(int i=0;i<NPly;i++) {
+		        switch(mat[i]->get_constmodel()){
+		            case EP:  // pure matrix
+		                if (mat[i]->get_pos_maxD()>-1)  // EP case
+		                {
+		                    allPos_mtx_maxD.emplace_back(idsdv + mat[i]->get_pos_maxD());
+		                }
+		                break;
+		                
+		            case MTSecF:  // MFH Fisrst order secant
+		                if(mat[i]->get_pos_mtx_maxD()>-1)  // MT-Mtx phase
+		                {
+		                    allPos_mtx_maxD.emplace_back(idsdv + mat[i]->get_pos_mtx_maxD());
+		                }
+		                break;
+		                
+		            case VT:
+		                std::vector<int> allPos_mat_mtx_maxD = mat[i]->get_allPos_mtx_maxD();
+		                for(int j=0;j<allPos_mat_mtx_maxD.size();j++)  // LAM_2PLY-Mtx phase
+		                {
+		                    allPos_mtx_maxD.emplace_back(idsdv + allPos_mat_mtx_maxD[j]);
+		                }
+		                break;		        
+		        }
+		        		        
+		        idsdv += mat[i]->get_nsdv();
+		    }
+		    
+		    return allPos_mtx_maxD;
+                }
+
+                virtual void populateMaxD(double *statev, double Dmax_Inc, double Dmax_Mtx) const {
+                    std::vector<int> allPos_inc_maxD = this->get_allPos_inc_maxD();
+		    for(int i=0;i<allPos_inc_maxD.size();i++) {
+		        statev[allPos_inc_maxD[i]] = Dmax_Inc;
+		    }
+		    std::vector<int> allPos_mtx_maxD = this->get_allPos_mtx_maxD();
+		    for(int i=0;i<allPos_mtx_maxD.size();i++) {
+		        statev[allPos_mtx_maxD[i]] = Dmax_Mtx;
+		    }
+                }*/
+
+		// for LAM_2PLY case (LVM)
+                // for each laminate ply
+                // strain
+                virtual int get_pos_mat_strain (int laminateNb) const
+                {
+                    int idsdv = 39;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<laminateNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[laminateNb]->get_pos_strn();
+		}
+		// stress
+		virtual int get_pos_mat_stress (int laminateNb) const
+                {
+                    int idsdv = 39;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<laminateNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[laminateNb]->get_pos_strs();
+		}
+		// damage
+		virtual int get_pos_mat_Dam (int laminateNb) const
+                {
+                    if(mat[laminateNb]->get_pos_damagefornonlocal()>-1){
+                        int idsdv = 39;  // position in statev where the statev of the first phase start
+                        for(int i=0;i<laminateNb;i++){
+                            idsdv += mat[i]->get_nsdv();
+                        }
+		        return idsdv + mat[laminateNb]->get_pos_damagefornonlocal();
+		    }
+		    else{
+		        return -1;
+		    }
+		}
+		// for each material phase
+		// strain
+		virtual int get_pos_mat_mtx_strain (int laminateNb) const
+                {
+                    int idsdv = 39;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<laminateNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[laminateNb]->get_pos_mtx_strain();
+		}
+		virtual int get_pos_mat_inc_strain (int laminateNb) const
+                {
+                    int idsdv = 39;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<laminateNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[laminateNb]->get_pos_inc_strain();
+		}
+		// stress
+		virtual int get_pos_mat_mtx_stress (int laminateNb) const
+                {
+                    int idsdv = 39;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<laminateNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[laminateNb]->get_pos_mtx_stress();
+		}
+		virtual int get_pos_mat_inc_stress (int laminateNb) const
+                {
+                    int idsdv = 39;  // position in statev where the statev of the first phase start
+                    for(int i=0;i<laminateNb;i++){
+                        idsdv += mat[i]->get_nsdv();
+                    }
+		    return idsdv + mat[laminateNb]->get_pos_inc_stress();
+		}
+		// damage
+		virtual int get_pos_mat_mtx_Dam (int laminateNb) const
+                {
+                    if(mat[laminateNb]->get_pos_mtx_Dam()>-1){
+                        int idsdv = 39;  // position in statev where the statev of the first phase start
+                        for(int i=0;i<laminateNb;i++){
+                            idsdv += mat[i]->get_nsdv();
+                        }
+		        return idsdv + mat[laminateNb]->get_pos_mtx_Dam();
+		    }
+		    else{
+		        return -1;
+		    }
+		}
+		virtual int get_pos_mat_inc_Dam (int laminateNb) const
+                {
+                    if(mat[laminateNb]->get_pos_inc_Dam()>-1){
+                        int idsdv = 39;  // position in statev where the statev of the first phase start
+                        for(int i=0;i<laminateNb;i++){
+                            idsdv += mat[i]->get_nsdv();
+                        }
+		        return idsdv + mat[laminateNb]->get_pos_inc_Dam();
+		    }
+		    else{
+		        return -1;
+		    }
+		}
+		
+		// for LAM_2PLY case (LVM)
+		// for each laminate ply
+		// strain
+		virtual int getPosStrnLamM0(int laminateNb) const {return (int)pos_strn_lam_m0[laminateNb];}
+		virtual int getPosStrnLamVT(int laminateNb) const {return (int)pos_strn_lam_vt[laminateNb];}
+		// stress
+		virtual int getPosStrsLamM0(int laminateNb) const {return (int)pos_strs_lam_m0[laminateNb];}
+		virtual int getPosStrsLamVT(int laminateNb) const {return (int)pos_strs_lam_vt[laminateNb];}
+		// damage
+		virtual int getPosDamLamM0(int laminateNb) const {return (int)pos_dam_lam_m0[laminateNb];}
+		// for each MT grain in VT ply
+		// strain
+		virtual int getPosStrnLamVTMT(int grainNb) const {return (int)pos_strn_lam_vt_mt[grainNb];}
+		// stress
+		virtual int getPosStrsLamVTMT(int grainNb) const {return (int)pos_strs_lam_vt_mt[grainNb];}
+		// for each material phase
+		// strain
+		virtual int getPosStrnLamVTMTMtx(int grainNb) const {return (int)pos_strn_lam_vt_mt_mtx[grainNb];}
+		virtual int getPosStrnLamVTMTInc(int grainNb) const {return (int)pos_strn_lam_vt_mt_inc[grainNb];}
+		// stress
+		virtual int getPosStrsLamVTMTMtx(int grainNb) const {return (int)pos_strs_lam_vt_mt_mtx[grainNb];}
+		virtual int getPosStrsLamVTMTInc(int grainNb) const {return (int)pos_strs_lam_vt_mt_inc[grainNb];}
+		// damage
+		virtual int getPosDamLamVTMTMtx(int grainNb) const {return (int)pos_dam_lam_vt_mt_mtx[grainNb];}
+		virtual int getPosDamLamVTMTInc(int grainNb) const {return (int)pos_dam_lam_vt_mt_inc[grainNb];}
 #endif
 
 };  
@@ -518,37 +1162,39 @@ class MT_Material : public Material{
 		virtual void get_elOp(double** Cel);
 		virtual void get_elOD(double* statev, double** Cel);
                 virtual void unload_step(double* dstrn, double* strs_n, double* statev, int kinc, int kstep);
+
                 virtual int get_pos_inc_maxD() const {
-                   if(icl_mat->get_pos_maxD()>=0)	
+                   if(icl_mat->get_pos_maxD()>-1)
 	              return icl_mat->get_pos_maxD()+idsdv_i;
                    return -1;
                 }
                 virtual double get_inc_MaxD(double *statev) const {
-                   if(icl_mat->get_pos_maxD()>=0)	
+                   if(icl_mat->get_pos_maxD()>-1)
 	              return statev[get_pos_inc_maxD()];
                    return -1.;
                 }
-
                 virtual void set_inc_MaxD(double *statev, double Dmax)
                 {
-                   if(icl_mat->get_pos_maxD()>=0)	
+                   if(icl_mat->get_pos_maxD()>-1)
 	              statev[get_pos_inc_maxD()]=Dmax;
                 }
+
                 virtual int get_pos_mtx_maxD() const {
-                   if(mtx_mat->get_pos_maxD()>=0)	
+                   if(mtx_mat->get_pos_maxD()>-1)
 	              return mtx_mat->get_pos_maxD()+idsdv_m;
                    return -1;
                 }
                 virtual double get_mtx_MaxD(double *statev) const {
-                   if(mtx_mat->get_pos_maxD()>=0)	
+                   if(mtx_mat->get_pos_maxD()>-1)
 	              return statev[get_pos_mtx_maxD()];
                    return -1.;
                 }
                 virtual void set_mtx_MaxD(double *statev, double Dmax)
                 {
-                   if(mtx_mat->get_pos_maxD()>=0)	
+                   if(mtx_mat->get_pos_maxD()>-1)
 	              statev[get_pos_mtx_maxD()]=Dmax;
                 }
+
 		virtual int getNbNonLocalVariables() const 
                 {
                    int nb=0;
@@ -556,6 +1202,10 @@ class MT_Material : public Material{
                    nb=nb+icl_mat->getNbNonLocalVariables(); 
                    return nb;
                 }
+                virtual int get_pos_Fd_bar()const;
+                virtual int get_pos_dFd_bar()const;
+                virtual int get_pos_locFd() const;
+
 
 #ifdef NONLOCALGMSH
 		virtual int get_pos_currentplasticstrainfornonlocal() const;
@@ -572,9 +1222,12 @@ class MT_Material : public Material{
 		
 		virtual void setEuler(double e1, double e2, double e3);
 
-
-
-
+		/*virtual void setLocalDamage()
+		{
+                    Material::setLocalDamage();
+                    mtx_mat->setLocalDamage();
+                    inc_mat->setLocalDamage();
+                }*/
 #endif
 };
 
@@ -588,9 +1241,18 @@ class MTSecF_Material : public MT_Material{
 		
 		MTSecF_Material(double* props, int idmat);
 		~MTSecF_Material();
+		
 
 		virtual int constbox(double* dstrn, double* strs_n,  double* strs, double* statev_n, double* statev, double** Cref, double*** dCref, double* tau, double** dtau, double** Calgo, double alpha, double* dpdE, double* strs_dDdp_bar, double *Sp_bar, double** c_g, double* dFd_d_bar, double* dstrs_dFd_bar, double *dFd_dE, double** c_gF,  double* dpdFd, double* dFddp, int kinc, int kstep, double dt);
                 virtual void unload_step(double* dstrn, double* strs_n, double* statev, int kinc, int kstep);
+                
+#ifdef NONLOCALGMSH
+		virtual void get_elOp_mtx(double* statev, double** Cel);
+		virtual void get_elOp_icl(double* statev, double** Cel);
+		virtual void get_elOp_mtx(double** Cel);
+		virtual void get_elOp_icl(double** Cel);
+                inline double get_vfI() {return vf_i;}
+#endif
 };
 
 //Composite material homogenized with Mori-Tanaka scheme with incremental secant method, second moment for the elastic predictor
@@ -649,6 +1311,9 @@ class MTSecF_Stoch_Material : public MT_Material{
 		virtual int get_pos_sy0()const;   
 		virtual int get_pos_hmod1()const;  
 		virtual int get_pos_hmod2()const;  
+		virtual int get_pos_alpha_DP ()const; 
+		virtual int get_pos_m_DP ()const;
+		virtual int get_pos_nup ()const;
 		virtual int get_pos_hp0()const;  
 		virtual int get_pos_hexp()const;  
 		virtual int get_pos_DamParm1()const;  
@@ -822,8 +1487,8 @@ class ANEL_Material : public Material{
                 double dDeldD;
                 double v12, v13, v31, E1, E3;
                 int pos_euler;
-                
-    
+
+
 	public:
 		
 		ANEL_Material(double* props, int idmat);
@@ -851,6 +1516,7 @@ class ANEL_Material : public Material{
 		virtual int get_pos_locFd() const {  return pos_locFd; }
                 virtual int get_pos_euler() const {  return pos_euler; }
                 
+                
 #ifdef NONLOCALGMSH
 		virtual int get_pos_currentplasticstrainfornonlocal() const;
 		virtual int get_pos_effectiveplasticstrainfornonlocal() const;
@@ -861,7 +1527,7 @@ class ANEL_Material : public Material{
 		virtual void setEuler(double e1, double e2, double e3);
 #endif
 		
-}; 
+};
 
 #ifdef NONLOCALGMSH
 #else
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.cpp
index d321f89867d632c4ac849148e1dd90dcab7db47c..788988aff97f0cab7a701b2006811de26502c2cf 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.cpp
@@ -27,6 +27,9 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
         pos_Msy0 = mat->get_pos_sy0();
         pos_Mhmod1 = mat->get_pos_hmod1();
         pos_Mhmod2 = mat->get_pos_hmod2();
+        pos_Malpha_DP = mat->get_pos_alpha_DP();
+        pos_Mm_DP = mat->get_pos_m_DP();
+        pos_Mnup = mat->get_pos_nup();
         pos_Mhexp = mat->get_pos_hexp();  
         pos_DamParm1 = mat->get_pos_DamParm1();
         pos_DamParm2 = mat->get_pos_DamParm2();
@@ -43,6 +46,9 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
         if(pos_Msy0 !=0) Randnum +=1;
         if(pos_Mhmod1 !=0) Randnum +=1;
         if(pos_Mhmod2 !=0) Randnum +=1;
+        if(pos_Malpha_DP !=0) Randnum +=1;
+        if(pos_Mm_DP !=0) Randnum +=1;
+        if(pos_Mnup !=0) Randnum +=1;
         if(pos_Mhexp !=0) Randnum +=1; 
         if(pos_euler !=0) Randnum +=1; 
         if(pos_DamParm1 !=0) Randnum +=1; 
@@ -87,6 +93,9 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
     pos_Msy0 = mat->get_pos_sy0();
     pos_Mhmod1 = mat->get_pos_hmod1();
     pos_Mhmod2 = mat->get_pos_hmod2();
+    pos_Malpha_DP = mat->get_pos_alpha_DP();
+    pos_Mm_DP = mat->get_pos_m_DP();
+    pos_Mnup = mat->get_pos_nup();
     pos_Mhexp = mat->get_pos_hexp();  
     pos_DamParm1 = mat->get_pos_DamParm1();
     pos_DamParm2 = mat->get_pos_DamParm2();
@@ -103,6 +112,9 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
     if(pos_Msy0 !=0) Randnum +=1;
     if(pos_Mhmod1 !=0) Randnum +=1;
     if(pos_Mhmod2 !=0) Randnum +=1;
+    if(pos_Malpha_DP !=0) Randnum +=1;
+    if(pos_Mm_DP !=0) Randnum +=1;
+    if(pos_Mnup !=0) Randnum +=1;
     if(pos_Mhexp !=0) Randnum +=1; 
     if(pos_euler !=0) Randnum +=1; 
     if(pos_DamParm1 !=0) Randnum +=1; 
@@ -115,7 +127,7 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
     //allocate memory if random variable are required
     // read random property variables from 
     if ( Randnum != 0){
-       double Rprop[12];
+       double Rprop[15];
        int nxyz[3];
        int k;
 
@@ -132,6 +144,9 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
              if(pos_Msy0 !=0) resizeFlag = _sy0_Mat.resize(nxyz[0], nxyz[1], true);
              if(pos_Mhmod1 !=0) resizeFlag = _hmod1_Mat.resize(nxyz[0], nxyz[1], true);
              if(pos_Mhmod2 !=0) resizeFlag = _hmod2_Mat.resize(nxyz[0], nxyz[1], true);
+             if(pos_Malpha_DP !=0) resizeFlag = _alpha_DP_Mat.resize(nxyz[0], nxyz[1], true);
+             if(pos_Mm_DP !=0) resizeFlag = _m_DP_Mat.resize(nxyz[0], nxyz[1], true);
+             if(pos_Mnup !=0) resizeFlag = _nup_Mat.resize(nxyz[0], nxyz[1], true);
              if(pos_Mhexp !=0) resizeFlag = _hexp_Mat.resize(nxyz[0], nxyz[1], true);
              if(pos_euler !=0) resizeFlag = _euler_Mat.resize(nxyz[0], nxyz[1], true);
              if(pos_DamParm1 !=0) resizeFlag = _dam_Parm1_Mat.resize(nxyz[0], nxyz[1], true);
@@ -143,8 +158,7 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
              if (resizeFlag){
                 for(int i=0; i<nxyz[0]; i++){
                    for(int j=0; j<nxyz[1]; j++){
-                      okf = fscanf(Props, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", &Rprop[0], &Rprop[1],&Rprop[2],
-                          &Rprop[3],&Rprop[4],&Rprop[5],&Rprop[6],&Rprop[7],&Rprop[8],&Rprop[9],&Rprop[10],&Rprop[11]);
+                      okf = fscanf(Props, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", &Rprop[0],&Rprop[1],&Rprop[2],&Rprop[3],&Rprop[4],&Rprop[5],&Rprop[6],&Rprop[7],&Rprop[8],&Rprop[9],&Rprop[10],&Rprop[11],&Rprop[12],&Rprop[13],&Rprop[14]);
                       k = 0;
                       if(pos_vfi !=0){
                          _VfMat.set(i, j, Rprop[k++]);
@@ -159,6 +173,9 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
                       if(pos_Msy0 !=0) _sy0_Mat.set(i, j, Rprop[k++]);
                       if(pos_Mhmod1 !=0) _hmod1_Mat.set(i, j, Rprop[k++]);
                       if(pos_Mhmod2 !=0) _hmod2_Mat.set(i, j, Rprop[k++]);
+                      if(pos_Malpha_DP !=0) _alpha_DP_Mat.set(i, j, Rprop[k++]);
+                      if(pos_Mm_DP !=0) _m_DP_Mat.set(i, j, Rprop[k++]);
+                      if(pos_Mnup !=0) _nup_Mat.set(i, j, Rprop[k++]);
                       if(pos_Mhexp !=0) _hexp_Mat.set(i, j, Rprop[k++]);
                       if(pos_euler !=0) _euler_Mat.set(i, j, Rprop[k++]);
                       if(pos_DamParm1 !=0) _dam_Parm1_Mat.set(i, j, Rprop[k++]);
@@ -286,12 +303,12 @@ mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const int num, const double r
 mlawNonLocalDamage_Stoch::mlawNonLocalDamage_Stoch(const mlawNonLocalDamage_Stoch &source):                                        
            mlawNonLocalDamage(source), pos_vfi(source.pos_vfi), pos_euler(source.pos_euler), pos_aspR(source.pos_aspR), 
            pos_ME(source.pos_ME), pos_Mnu(source.pos_Mnu), pos_Msy0(source.pos_Msy0), pos_Mhmod1(source.pos_Mhmod1), 
-           pos_Mhmod2(source.pos_Mhmod2), pos_Mhexp(source.pos_Mhexp), pos_DamParm1(source.pos_DamParm1), pos_DamParm2(source.pos_DamParm2), pos_DamParm3(source.pos_DamParm3),  
+           pos_Mhmod2(source.pos_Mhmod2), pos_Malpha_DP(source.pos_Malpha_DP), pos_Mm_DP(source.pos_Mm_DP), pos_Mnup(source.pos_Mnup), pos_Mhexp(source.pos_Mhexp), pos_DamParm1(source.pos_DamParm1), pos_DamParm2(source.pos_DamParm2), pos_DamParm3(source.pos_DamParm3),  
            pos_INCDamParm1(source.pos_INCDamParm1), pos_INCDamParm2(source.pos_INCDamParm2), pos_INCDamParm3(source.pos_INCDamParm3), Randnum(source.Randnum), _intpl(source._intpl), 
            _dx(source._dx), _dy(source._dy), _dz(source._dz), _OrigX(source._OrigX), _OrigY(source._OrigY), 
            _OrigZ(source._OrigZ), _Ldomain(source._Ldomain), _VfMat(source._VfMat), _aspRMat(source._aspRMat), 
            _E_Mat(source._E_Mat), _nu_Mat(source._nu_Mat), _sy0_Mat(source._sy0_Mat), _hmod1_Mat(source._hmod1_Mat),
-           _hmod2_Mat(source._hmod2_Mat), _hexp_Mat(source._hexp_Mat),_euler_Mat(source._euler_Mat),
+           _hmod2_Mat(source._hmod2_Mat), _alpha_DP_Mat(source._alpha_DP_Mat), _m_DP_Mat(source._m_DP_Mat),_nup_Mat(source._nup_Mat),_hexp_Mat(source._hexp_Mat),_euler_Mat(source._euler_Mat),
            _dam_Parm1_Mat(source._dam_Parm1_Mat),_dam_Parm2_Mat(source._dam_Parm2_Mat) ,_dam_Parm3_Mat(source._dam_Parm3_Mat){
 
      Reuler = NULL;
@@ -317,6 +334,9 @@ mlawNonLocalDamage_Stoch& mlawNonLocalDamage_Stoch::operator=(const materialLaw
      pos_Msy0 = src->pos_Msy0;
      pos_Mhmod1 = src->pos_Mhmod1;
      pos_Mhmod2 = src->pos_Mhmod2;
+     pos_Malpha_DP = src->pos_Malpha_DP;
+     pos_Mm_DP = src->pos_Mm_DP;
+     pos_Mnup = src->pos_Mnup;
      pos_Mhexp = src->pos_Mhexp;
      pos_DamParm1 = src->pos_DamParm1;
      pos_DamParm2 = src->pos_DamParm2;
@@ -346,6 +366,12 @@ mlawNonLocalDamage_Stoch& mlawNonLocalDamage_Stoch::operator=(const materialLaw
      _sy0_Mat = src->_sy0_Mat;
      _hmod1_Mat = src->_hmod1_Mat;
      _hmod2_Mat = src->_hmod2_Mat;
+
+     _alpha_DP_Mat = src->_alpha_DP_Mat;
+     _m_DP_Mat = src->_m_DP_Mat;
+     _nup_Mat = src->_nup_Mat;
+
+
      _hexp_Mat = src->_hexp_Mat;
      _euler_Mat = src->_euler_Mat;
 
@@ -435,6 +461,9 @@ void mlawNonLocalDamage_Stoch::createIPState(const SVector3 &GaussP, IPNonLocalD
       ivi->setPos_Msy0(pos_Msy0);
       ivi->setPos_Mhmod1(pos_Mhmod1);
       ivi->setPos_Mhmod2(pos_Mhmod2);
+      ivi->setPos_Malpha_DP(pos_Malpha_DP);
+      ivi->setPos_Mm_DP(pos_Mm_DP);
+      ivi->setPos_Mnup(pos_Mnup);
       ivi->setPos_Mhexp(pos_Mhexp);
       ivi->setPos_DamParm1(pos_DamParm1);
       ivi->setPos_DamParm2(pos_DamParm2);
@@ -452,6 +481,9 @@ void mlawNonLocalDamage_Stoch::createIPState(const SVector3 &GaussP, IPNonLocalD
       iv1->setPos_Msy0(pos_Msy0);
       iv1->setPos_Mhmod1(pos_Mhmod1);
       iv1->setPos_Mhmod2(pos_Mhmod2);
+      iv1->setPos_Malpha_DP(pos_Malpha_DP);
+      iv1->setPos_Mm_DP(pos_Mm_DP);
+      iv1->setPos_Mnup(pos_Mnup);
       iv1->setPos_Mhexp(pos_Mhexp);
       iv1->setPos_DamParm1(pos_DamParm1);
       iv1->setPos_DamParm2(pos_DamParm2);
@@ -469,6 +501,9 @@ void mlawNonLocalDamage_Stoch::createIPState(const SVector3 &GaussP, IPNonLocalD
       iv2->setPos_Msy0(pos_Msy0);
       iv2->setPos_Mhmod1(pos_Mhmod1);
       iv2->setPos_Mhmod2(pos_Mhmod2);
+      iv2->setPos_Malpha_DP(pos_Malpha_DP);
+      iv2->setPos_Mm_DP(pos_Mm_DP);
+      iv2->setPos_Mnup(pos_Mnup);
       iv2->setPos_Mhexp(pos_Mhexp);
       iv2->setPos_DamParm1(pos_DamParm1);
       iv2->setPos_DamParm2(pos_DamParm2);
@@ -493,6 +528,9 @@ void mlawNonLocalDamage_Stoch::createIPState(const SVector3 &GaussP, IPNonLocalD
       if(pos_Msy0 !=0) stv1[pos_Msy0] = Rprop[k++];
       if(pos_Mhmod1 !=0) stv1[pos_Mhmod1] = Rprop[k++];
       if(pos_Mhmod2 !=0) stv1[pos_Mhmod2] = Rprop[k++];
+      if(pos_Malpha_DP !=0) stv1[pos_Malpha_DP] = Rprop[k++];
+      if(pos_Mm_DP !=0) stv1[pos_Mm_DP] = Rprop[k++];
+      if(pos_Mnup !=0) stv1[pos_Mnup] = Rprop[k++];
       if(pos_Mhexp !=0) stv1[pos_Mhexp] = Rprop[k++];
       if(pos_DamParm1 !=0) stv1[pos_DamParm1] = Rprop[k++];
       if(pos_DamParm2 !=0) stv1[pos_DamParm2] = Rprop[k++];
@@ -518,6 +556,9 @@ void mlawNonLocalDamage_Stoch::createIPState(const SVector3 &GaussP, IPNonLocalD
        if(pos_Msy0 !=0) stv1[pos_Msy0] = _sy0_Mat(nx_int,ny_int);
        if(pos_Mhmod1 !=0) stv1[pos_Mhmod1] = _hmod1_Mat(nx_int,ny_int);
        if(pos_Mhmod2 !=0) stv1[pos_Mhmod2] = _hmod2_Mat(nx_int,ny_int);
+       if(pos_Malpha_DP !=0) stv1[pos_Malpha_DP] = _alpha_DP_Mat(nx_int,ny_int);
+       if(pos_Mm_DP !=0) stv1[pos_Mm_DP] = _m_DP_Mat(nx_int,ny_int);
+       if(pos_Mnup !=0) stv1[pos_Mnup] = _nup_Mat(nx_int,ny_int);
        if(pos_Mhexp !=0)  stv1[pos_Mhexp] = _hexp_Mat(nx_int,ny_int);
        if(pos_DamParm1 !=0) stv1[pos_DamParm1] = _dam_Parm1_Mat(nx_int,ny_int);
        if(pos_DamParm2 !=0) stv1[pos_DamParm2] = _dam_Parm2_Mat(nx_int,ny_int);
@@ -865,6 +906,9 @@ void mlawNonLocalDamage_Stoch::createIPVariable(const SVector3 &GaussP, IPNonLoc
       ipv->setPos_Msy0(pos_Msy0);
       ipv->setPos_Mhmod1(pos_Mhmod1);
       ipv->setPos_Mhmod2(pos_Mhmod2);
+      ipv->setPos_Malpha_DP(pos_Malpha_DP);
+      ipv->setPos_Mm_DP(pos_Mm_DP);
+      ipv->setPos_Mnup(pos_Mnup);
       ipv->setPos_Mhexp(pos_Mhexp);
       ipv->setPos_DamParm1(pos_DamParm1);
       ipv->setPos_DamParm2(pos_DamParm2);
@@ -888,6 +932,9 @@ void mlawNonLocalDamage_Stoch::createIPVariable(const SVector3 &GaussP, IPNonLoc
       if(pos_Msy0 !=0) stv1[pos_Msy0] = Rprop[k++];
       if(pos_Mhmod1 !=0) stv1[pos_Mhmod1] = Rprop[k++];
       if(pos_Mhmod2 !=0) stv1[pos_Mhmod2] = Rprop[k++];
+      if(pos_Malpha_DP !=0) stv1[pos_Malpha_DP] = Rprop[k++];
+      if(pos_Mm_DP !=0) stv1[pos_Mm_DP] = Rprop[k++];
+      if(pos_Mnup !=0) stv1[pos_Mnup] = Rprop[k++];
       if(pos_Mhexp !=0) stv1[pos_Mhexp] = Rprop[k++];
       if(pos_DamParm1 !=0) stv1[pos_DamParm1] = Rprop[k++];
       if(pos_DamParm2 !=0) stv1[pos_DamParm2] = Rprop[k++];
@@ -912,6 +959,9 @@ void mlawNonLocalDamage_Stoch::createIPVariable(const SVector3 &GaussP, IPNonLoc
        if(pos_Msy0 !=0) stv1[pos_Msy0] = _sy0_Mat(nx_int,ny_int);
        if(pos_Mhmod1 !=0) stv1[pos_Mhmod1] = _hmod1_Mat(nx_int,ny_int);
        if(pos_Mhmod2 !=0) stv1[pos_Mhmod2] = _hmod2_Mat(nx_int,ny_int);
+       if(pos_Malpha_DP !=0) stv1[pos_Malpha_DP] = _alpha_DP_Mat(nx_int,ny_int);
+       if(pos_Mm_DP !=0) stv1[pos_Mm_DP] = _m_DP_Mat(nx_int,ny_int);
+       if(pos_Mnup !=0) stv1[pos_Mnup] = _nup_Mat(nx_int,ny_int);
        if(pos_Mhexp !=0)  stv1[pos_Mhexp] = _hexp_Mat(nx_int,ny_int);
        if(pos_DamParm1 !=0) stv1[pos_DamParm1] = _dam_Parm1_Mat(nx_int,ny_int);
        if(pos_DamParm2 !=0) stv1[pos_DamParm2] = _dam_Parm2_Mat(nx_int,ny_int);
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.h b/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.h
index d543daf0127427f7960bc9c582a3eaa7facc8456..c78a68b518d5d52b92ae09d64a53c5534fc35d0e 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.h
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamage_Stoch.h
@@ -33,6 +33,9 @@ class mlawNonLocalDamage_Stoch : public mlawNonLocalDamage
   int pos_Msy0;
   int pos_Mhmod1;
   int pos_Mhmod2;
+  int pos_Malpha_DP;
+  int pos_Mm_DP;
+  int pos_Mnup;
   int pos_Mhexp;
   int pos_DamParm1;
   int pos_DamParm2;
@@ -63,6 +66,9 @@ class mlawNonLocalDamage_Stoch : public mlawNonLocalDamage
   fullMatrix<double> _E_Mat;
   fullMatrix<double> _nu_Mat;
   fullMatrix<double> _sy0_Mat;
+  fullMatrix<double> _alpha_DP_Mat;
+  fullMatrix<double> _m_DP_Mat;
+  fullMatrix<double> _nup_Mat;
 
   fullMatrix<double> _hmod1_Mat;
   fullMatrix<double> _hmod2_Mat;