diff --git a/NonLinearSolver/materialLaw/mlaw.h b/NonLinearSolver/materialLaw/mlaw.h
index d440e62ad6564f0ef2d6b6a55c84a982076529a0..7da8c7a5640b10b36e389a72e1afd8a859b1b08d 100644
--- a/NonLinearSolver/materialLaw/mlaw.h
+++ b/NonLinearSolver/materialLaw/mlaw.h
@@ -33,7 +33,7 @@ class materialLaw{
                  hyperelastic, powerYieldLaw, powerYieldLawWithFailure,nonLocalDamagePowerYieldHyper,
                  localDamagePowerYieldHyperWithFailure,nonLocalDamagePowerYieldHyperWithFailure,ElecSMP,
                  ThermalConducter,AnIsotropicTherMech, localDamageJ2Hyper,linearElastic,nonLocalDamageGursonThermoMechanics,
-                 localDamageJ2SmallStrain,nonLocalDamageJ2SmallStrain,cluster,tfa,ANN, DMN, torchANN, NonlocalDamageTorchANN, LinearElecMagTherMech, LinearElecMagInductor, hyperviscoelastic,
+                 localDamageJ2SmallStrain,nonLocalDamageJ2SmallStrain,cluster,tfa,ANN, DMN, torchANN, NonlocalDamageTorchANN, StochDMN, LinearElecMagTherMech, LinearElecMagInductor, hyperviscoelastic,
                  GenericThermoMechanics, ElecMagGenericThermoMechanics, ElecMagInductor,
                  nonlineartvm,nonlinearTVE,nonlinearTVP,vevpmfh, VEVPUMAT, IMDEACPUMAT, Hill48,nonlinearTVEnonlinearTVP,
                  MultipleLaws};
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index c5c3875ee33b5b80b2d4da250943bc42333ba580..92e67aa01818038d541a5c85b54fe8ae7dd63ea7 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -2154,6 +2154,134 @@ double ANNBasedDG3DIPVariable::get(const int comp) const
   }
 }
 
+
+StochDMNDG3DIPVariable::StochDMNDG3DIPVariable(const int NRoot,const int NInterface, const bool createBodyForceHO, const bool oninter): 
+    dG3DIPVariable(createBodyForceHO, oninter), _NRoot(NRoot), _NInterface(NInterface), _IPVector(NRoot,NULL){
+    bool resizeFlag;
+    int entryFP = 9*_NRoot;
+    int entry_a = 3*_NInterface;
+    resizeFlag = Fvec.resize(entryFP, true); 
+    resizeFlag = Pvec.resize(entryFP, true); 
+    resizeFlag = a.resize(entry_a, true);     
+ }   
+    
+StochDMNDG3DIPVariable::StochDMNDG3DIPVariable(const StochDMNDG3DIPVariable& src): 
+  dG3DIPVariable(src), _NRoot(src._NRoot), Fvec(src.Fvec), Pvec(src.Pvec), a(src.a),_IPVector(src._NRoot,NULL){
+    for (int j=0; j< src._NRoot; j++){
+      if (src._IPVector[j] != NULL){
+        _IPVector[j] = src._IPVector[j]->clone();
+      }
+    }
+  }
+
+
+StochDMNDG3DIPVariable& StochDMNDG3DIPVariable::operator =(const IPVariable& src)
+{
+  dG3DIPVariable::operator=(src);
+  const StochDMNDG3DIPVariable* psrc = dynamic_cast<const StochDMNDG3DIPVariable*>(&src);
+  if (psrc!=NULL)
+  {
+    _NRoot = psrc->_NRoot;
+    Fvec = psrc->Fvec;
+    Pvec = psrc->Pvec;
+    a = psrc->a;
+    for (int j=0; j< _NRoot; j++)
+    {
+      if ((_IPVector[j]!=NULL) and (psrc->_IPVector[j] != NULL))
+      {
+        _IPVector[j]->operator =(*(psrc->_IPVector[j]));
+      }
+      else
+      {
+        Msg::Error("cannot perform equal operator in StochDMNDG3DIPVariable::operator =");
+      }
+    }
+  }
+  return *this;
+}
+
+StochDMNDG3DIPVariable::~StochDMNDG3DIPVariable()
+{
+  for (int j=0; j< _NRoot; j++)
+  {
+    if (_IPVector[j] != NULL)
+    {
+      delete _IPVector[j];
+    }
+  }
+}
+
+void StochDMNDG3DIPVariable::addIPv(int loc, IPVariable* ipv)
+{
+  _IPVector[loc] = ipv;
+}
+
+double StochDMNDG3DIPVariable::get(const int comp) const
+{
+  if (comp == IPField::DEFO_ENERGY ||
+      comp == IPField::PLASTIC_ENERGY || 
+      comp == IPField::DAMAGE_ENERGY ||
+      comp == IPField::DAMAGE )
+  {
+    double v = 0;
+    for (int j=0; j< _NRoot; j++)
+    {
+      v += _IPVector[j]->get(comp);
+    }
+    return v;
+  }
+  else
+  {
+    return dG3DIPVariable::get(comp);
+  }
+};
+
+
+void StochDMNDG3DIPVariable::restart()
+{
+  dG3DIPVariable::restart();
+  restartManager::restart(_NRoot);
+  restartManager::restart(Fvec.getDataPtr(),9*_NRoot);
+  restartManager::restart(Pvec.getDataPtr(),9*_NRoot);
+  restartManager::restart(a.getDataPtr(),3*_NInterface);
+  restartManager::restart(_IPVector);
+  for (int j=0; j< _NRoot; j++){
+    _IPVector[j]->restart();
+  }
+};
+
+
+
+double StochDMNDG3DIPVariable::defoEnergy() const
+{
+  double v = 0;
+  for (int j=0; j< _NRoot; j++)
+  {
+    v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->defoEnergy();
+  }
+  return v;
+}
+double StochDMNDG3DIPVariable::plasticEnergy() const
+{
+  double v = 0;
+  for (int j=0; j< _NRoot; j++)
+  {
+    v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->plasticEnergy();
+  }
+  return v;
+}
+double StochDMNDG3DIPVariable::damageEnergy() const
+{
+  double v = 0;
+  for (int j=0; j< _NRoot; j++)
+  {
+    v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->damageEnergy();
+  }
+  return v;
+}
+
+
+
 //torchANNBasedDG3DIPVariable::torchANNBasedDG3DIPVariable(const int n, const double initial_h, const bool createBodyForceHO, const bool oninter):
 //                 dG3DIPVariable(createBodyForceHO, oninter), _n(n), _initial_h(initial_h), restart_internalVars(1,n), _kinematicVariables(1,6)
 //{
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index ebb74105c74d4d6a8459d9653d74dcc5d8e9da1d..e8c05c2b14f7e34c42ea0958e053d454747b8d6f 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -2393,6 +2393,48 @@ class ANNBasedDG3DIPVariable :public dG3DIPVariable
     virtual void restart();
 };
 
+
+class StochDMNDG3DIPVariable : public dG3DIPVariable
+{
+  protected:
+    int _NRoot, _NInterface;
+    std::vector<IPVariable*> _IPVector;
+    fullVector<double> Fvec; 
+    fullVector<double> Pvec;
+    fullVector<double> a;
+
+ public:      
+    StochDMNDG3DIPVariable(const int NRoot,const int NInterface, const bool createBodyForceHO, const bool oninter);
+    StochDMNDG3DIPVariable(const StochDMNDG3DIPVariable& src);
+    virtual StochDMNDG3DIPVariable& operator =(const IPVariable& src);
+    virtual ~StochDMNDG3DIPVariable();
+
+    void InitializeFluctuationVector() {a.setAll(0.0);};
+    virtual double defoEnergy() const;
+    virtual double plasticEnergy() const;
+    virtual double damageEnergy() const;
+    
+    void addIPv(int i, IPVariable* ipv);
+    IPVariable* getIPv(int i) {return _IPVector[i];};
+    const IPVariable* getIPv(int i) const {return _IPVector[i];};
+    
+    virtual double get(const int comp) const;
+    virtual fullVector<double>& getRefToDeformationGradientVect() {return Fvec;};
+    virtual const fullVector<double>&  getConstRefToDeformationGradientVect() const {return Fvec;}; 
+    virtual fullVector<double>& getRefToFirstPiolaKirchhoffStressVect() {return Pvec;};
+    virtual const fullVector<double>& getConstRefToFirstPiolaKirchhoffStressVect() const {return Pvec;};
+    virtual fullVector<double>& getRefToInterfaceFluctuationVect() {return a;};
+    virtual const fullVector<double>& getConstRefToInterfaceFluctuationVect() const {return a;}; 
+    
+    
+    
+    virtual IPVariable* getInternalData() {return NULL;};
+    virtual const IPVariable* getInternalData()  const {return NULL;};
+    virtual IPVariable* clone() const {return new StochDMNDG3DIPVariable(*this);};
+    virtual void restart();
+};
+
+
 class torchANNBasedDG3DIPVariable :public dG3DIPVariable
 {
   protected:
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 6ed9c64645d0c5bec195a60db135b56d76f3c012..880ee1733b6f08cc9174c0043c631f41d9819d00 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -1368,7 +1368,7 @@ void DMNBasedDG3DMaterialLaw::initialIPVariable(IPVariable* ipv, bool stiff)
 
 void DMNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent)
 {
-  DMNBasedDG3DIPVariable* ipvcur = static_cast<DMNBasedDG3DIPVariable*>(ipv);;
+  DMNBasedDG3DIPVariable* ipvcur = static_cast<DMNBasedDG3DIPVariable*>(ipv);
   const DMNBasedDG3DIPVariable* ipvprev = static_cast<const DMNBasedDG3DIPVariable*>(ipvp);
 
   dG3DDeepMaterialNetwork* deepMN = dynamic_cast<dG3DDeepMaterialNetwork*>(ipvcur->getDeepMaterialNetwork());
@@ -3428,6 +3428,1502 @@ double torchANNBasedDG3DMaterialLaw::soundSpeed() const
   return 0.;
 };
 
+
+
+StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile,const bool porous,const double tol):
+  dG3DMaterialLaw(num,rho,false),_porous(porous),_tol(tol){
+          
+  read_parameter(ParaFile);
+  fill_Matrices(_Vf);         
+}
+    
+StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const StochDMNDG3DMaterialLaw &src): dG3DMaterialLaw(src), _porous(src._porous), _TotalLevel(src._TotalLevel), _NRoot(src._NRoot),
+       _NInterface(src._NInterface), _N_Node(src._N_Node), _NPara_Wt(src._NPara_Wt), _NPara_Norm(src._NPara_Norm), _Dim(src._Dim), _Vf(src._Vf), _tol(src._tol),
+       _Nv(src._Nv), _NTV(src._NTV), _V(src._V), _Para_Wt(src._Para_Wt), _Para_Norm(src._Para_Norm), _VA(src._VA), _VI(src._VI), _mapLaw(src._mapLaw){ }
+    
+StochDMNDG3DMaterialLaw::~StochDMNDG3DMaterialLaw(){};
+     
+          
+void StochDMNDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const{
+      // check interface element
+  bool inter=true;
+  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+  if(iele==NULL) inter=false;
+    IPVariable* ipvi = new  StochDMNDG3DIPVariable(_NRoot, _NInterface, hasBodyForce,inter);
+    IPVariable* ipv1 = new  StochDMNDG3DIPVariable(_NRoot, _NInterface, hasBodyForce,inter);
+    IPVariable* ipv2 = new  StochDMNDG3DIPVariable(_NRoot, _NInterface, hasBodyForce,inter);
+        
+    StochDMNDG3DIPVariable* ipvi_root = static_cast<StochDMNDG3DIPVariable*>(ipvi);
+    StochDMNDG3DIPVariable* ipv1_root = static_cast<StochDMNDG3DIPVariable*>(ipv1);
+    StochDMNDG3DIPVariable* ipv2_root = static_cast<StochDMNDG3DIPVariable*>(ipv2);
+    for (int i=0; i< _NRoot; i++){
+      int mat_IP = 0;
+      IPStateBase* ips_i = NULL;
+      if(_Modified_DMN){
+        if(i%3 == 1) mat_IP = 1;  
+      }   
+      else{
+         mat_IP = i%2;  
+      }    
+      if(_porous) mat_IP=0; 
+      _mapLaw[mat_IP]->createIPState(ips_i, hasBodyForce, state_, ele, nbFF_, GP, gpt);
+      std::vector<IPVariable*> ip_all;
+      ips_i->getAllIPVariable(ip_all);
+      ipvi_root->addIPv(i, ip_all[0]);
+      ipv1_root->addIPv(i, ip_all[1]);
+      ipv2_root->addIPv(i, ip_all[2]);
+    } 
+    if(ips != NULL) delete ips;
+      ips = new IP3State(state_,ipvi,ipv1,ipv2);
+}
+       
+   
+void StochDMNDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const {
+  if(ipv !=NULL) delete ipv;
+  bool inter=true;
+  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+  if(iele == NULL){
+    inter=false;
+  }
+  ipv = new  StochDMNDG3DIPVariable(_NRoot, _NInterface, hasBodyForce,inter);
+  StochDMNDG3DIPVariable* ipv_root = static_cast<StochDMNDG3DIPVariable*>(ipv);
+  for (int i=0; i< _NRoot; i++){
+    IPVariable* ipv_i = NULL;
+    int mat_IP =0;
+    if(_Modified_DMN){
+      if(i%3 == 1) mat_IP = 1;  
+    }   
+    else{
+       mat_IP = i%2;  
+    }  
+    if(_porous) mat_IP=0; 
+    _mapLaw[mat_IP]->createIPVariable(ipv_i, hasBodyForce, ele, nbFF_, GP, gpt);
+    ipv_root->addIPv(i, ipv_i);
+  }      
+};
+    
+void StochDMNDG3DMaterialLaw::initialIPVariable(IPVariable* ipv, bool stiff){
+  StochDMNDG3DIPVariable* ipv_all = dynamic_cast<StochDMNDG3DIPVariable*>(ipv);
+  ipv_all->InitializeFluctuationVector();
+  if (ipv_all == NULL){
+    Msg::Error("StochDMNDG3DIPVariable must be used in StochDMNDG3DMaterialLaw::initialIPVariable");
+  }
+  for (int i=0; i< _NRoot; i++){
+    IPVariable* ipv_i = ipv_all->getIPv(i);  
+    int mat_IP =0;      
+    if(_Modified_DMN){
+      if(i%3 == 1) mat_IP = 1;  
+    }   
+    else{
+        mat_IP = i%2;  
+    }   
+    if(_porous) mat_IP=0; 
+    _mapLaw[mat_IP]->initialIPVariable(ipv_i, stiff);
+  }
+}  
+        
+void StochDMNDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{
+  StochDMNDG3DIPVariable* ipv_all = dynamic_cast<StochDMNDG3DIPVariable*>(ipv);
+  const StochDMNDG3DIPVariable* ipvprev_all = dynamic_cast<const StochDMNDG3DIPVariable*>(ipvprev);
+  if (ipv_all == NULL){
+    Msg::Error("StochDMNDG3DIPVariable must be used in StochDMNDG3DMaterialLaw::checkInternalState");
+  }
+  for (int i=0; i< _NRoot; i++){
+    IPVariable* ipv_i = ipv_all->getIPv(i);
+    const IPVariable* ipvprev_i = ipvprev_all->getIPv(i);
+    int mat_IP = 0;
+    if(_Modified_DMN){
+      if(i%3 == 1) mat_IP = 1;  
+    }  
+    else{
+       mat_IP = i%2;  
+    }  
+    if(_porous) mat_IP=0; 
+    _mapLaw[mat_IP]->checkInternalState(ipv_i, ipvprev_i);
+  }
+}
+
+void StochDMNDG3DMaterialLaw::pos_kids(const int pos, const int level, int& pos0, int& pos1){
+  if(!_Modified_DMN or level<_TotalLevel-2 ){
+      pos0 = 2*pos+1;
+      pos1 = 2*pos+2;
+  }
+  else if(level ==_TotalLevel-2){
+    pos0 = pos + pow(2, level);
+    int k = pos - (pow(2, level)-1);
+    pos1 =  int(3*pow(2, level)-1) + k*3+2;
+  }
+  else if(level==_TotalLevel-1){
+    int k = pos - (pow(2, level)-1);
+    pos0 = int(3*pow(2, level-1)-1) +k*3;
+    pos1 = pos0+1;
+  }
+}
+
+     
+void StochDMNDG3DMaterialLaw::TwoPhasesInteract(const STensor43& C0, const STensor43& C1, fullMatrix<double> &mat, const int pos, const int level){
+  double n0 = _Para_Norm[pos][0];
+  double n1 = _Para_Norm[pos][1];
+  double n2 = _Para_Norm[pos][2];
+  double tmp1[6][6], tmp2[6][6], Norm[6][3], NormT[3][6];
+  double S_ba = 1.0;
+  int pos0=0;
+  int pos1=0; 
+  if(level ==_TotalLevel-1 and _porous){
+    S_ba = _Para_Wt[pos][1];
+  }
+       
+  for(int i=0; i<3; i++){
+    for(int j=0; j<3; j++){
+      Norm[i][j] = 0.0;
+      Norm[3+i][j] = 0.0;
+      NormT[i][j] = 0.0;
+      NormT[i][3+j] = 0.0;
+    }
+  }        
+  Norm[0][0] = n0;        // different N and NT are used in order to be consistent with function "fromSTensor43ToFullMatrix()"
+  Norm[1][1] = n1;
+  Norm[2][2] = n2;
+  Norm[3][0] = n1/2.0;
+  Norm[3][1] = n0/2.0;
+  Norm[4][0] = n2/2.0;
+  Norm[4][2] = n0/2.0;
+  Norm[5][1] = n2/2.0;
+  Norm[5][2] = n1/2.0;  
+       
+  NormT[0][0] = n0;
+  NormT[1][1] = n1;
+  NormT[2][2] = n2;
+  NormT[0][3] = n1;
+  NormT[0][4] = n2;
+  NormT[1][3] = n0;
+  NormT[1][5] = n2;
+  NormT[2][4] = n0;       
+  NormT[2][5] = n1;    
+       
+  fullMatrix<double> mat0(6,6), mat1(6,6), K(3,3);  
+  pos_kids(pos, level, pos0, pos1); 
+  double v_A = _VA[pos0];
+  double v_B = _VA[pos1]; 
+         
+  if(level == _TotalLevel-1){
+    STensorOperation::fromSTensor43ToFullMatrix(C0, mat0);
+    STensorOperation::fromSTensor43ToFullMatrix(C1, mat1);       
+  }
+  else if(_Modified_DMN and level == _TotalLevel-2){
+    TwoPhasesInteract(C0, C1, mat0, pos0, level+1);
+    STensorOperation::fromSTensor43ToFullMatrix(C0, mat1);
+  } 
+  else{
+    TwoPhasesInteract(C0, C1, mat0, pos0, level+1);
+    TwoPhasesInteract(C0, C1, mat1, pos1, level+1);       
+  }
+       
+
+  for(int i=0; i<6;i++){
+    for(int j=0; j<6;j++){
+      tmp1[i][j] = mat0(i,j)/v_A + S_ba*S_ba*mat1(i,j)/v_B;      
+      tmp2[i][j] = S_ba*mat1(i,j) - mat0(i,j);      
+    }  
+  }
+       
+  for(int m=0; m<3;m++){
+    for(int n=0; n<3;n++) {
+      for(int i=0; i<6;i++) {
+        for(int j=0; j<6;j++) { 
+          K(m,n) += NormT[m][i]*tmp1[i][j]*Norm[j][n];
+        }
+      }
+    }
+  }         
+  K.invertInPlace();
+  for(int m=0; m<6;m++){
+    for(int n=0; n<6;n++) {
+      tmp1[m][n] = 0.0;
+      for(int i=0; i<3;i++) {
+        for(int j=0; j<3;j++) { 
+          tmp1[m][n] += Norm[m][i]*K(i,j)*NormT[j][n];
+        }
+      }
+    }
+  }
+       
+  for(int m=0; m<6;m++){
+    for(int n=0; n<6; n++) {
+      mat(m,n) = mat0(m,n)*v_A + mat1(m,n)*v_B;
+      for(int i=0; i<6;i++) {
+        for(int j=0; j<6;j++) { 
+          mat(m,n) -= tmp2[m][i]*tmp1[i][j]*tmp2[j][n];        
+        }  
+      }
+    }
+  }   
+}
+  // end of TwoPhasesInteract    
+      
+void StochDMNDG3DMaterialLaw::initLaws(const std::map<int,materialLaw*> &maplaw){
+  if (!_initialized){
+    STensorOperation::zero(elasticStiffness);
+    fullMatrix<double> mat(6,6);
+    std::vector<STensor43> C;    
+    for (std::vector<dG3DMaterialLaw*>::iterator ite = _mapLaw.begin(); ite!= _mapLaw.end(); ite++){
+      (*ite)->initLaws(maplaw);
+      elasticStiffness = (*ite)->elasticStiffness;
+      C.push_back(elasticStiffness);
+    }
+    if(_porous) C.push_back(elasticStiffness);
+    TwoPhasesInteract(C[0], C[1],  mat, 0, 0);
+    
+    STensorOperation::fromFullMatrixToSTensor43(mat, elasticStiffness);
+   // elasticStiffness.print("elasticStiffness StochDMNDG3DMaterialLaw");
+    _initialized = true;
+   // mat.print("initialize material tensor");
+  }  
+}  
+
+void StochDMNDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
+   dG3DMaterialLaw::setMacroSolver(sv);
+    for (std::vector<dG3DMaterialLaw*>::iterator ite = _mapLaw.begin(); ite!= _mapLaw.end(); ite++){
+      (*ite)->setMacroSolver(sv);
+    }
+}
+         
+void StochDMNDG3DMaterialLaw::fill_NLocal(const int pos, double N[][3]){
+  double n0 = _Para_Norm[pos][0];
+  double n1 = _Para_Norm[pos][1];
+  double n2 = _Para_Norm[pos][2];
+   
+  for(int i=0; i<9; i++){
+    for(int j=0; j<3; j++){
+      N[i][j]=0.0;
+    }
+  }    
+  N[0][0]= n0;
+  N[1][0]= n1;
+  N[2][0]= n2;
+  N[3][1]= n0;
+  N[4][1]= n1;
+  N[5][1]= n2;
+  N[6][2]= n0;
+  N[7][2]= n1;
+  N[8][2]= n2;
+}   
+
+
+void StochDMNDG3DMaterialLaw::fill_W_WI(const int pos, const int level, double W[][2]){
+  double VI = _VI[pos];
+  double Wtmp0 = _Para_Wt[pos][0];
+  double Wtmp1 = _Para_Wt[pos][1];
+  double WA = Wtmp0*(1.0-VI)+ Wtmp1*VI;
+  
+  for(int i=0; i<2; i++){
+    for(int j=0; j<2; j++){
+      W[i][j] = 0.0;
+    }
+  }      
+  W[0][0]= WA;
+  W[0][1]= 1.0-WA;
+  W[1][0]= (Wtmp1*VI)/WA;
+  W[1][1]= ((1.0-Wtmp1)*VI)/(1.0-WA);      
+} 
+
+     
+// fill the matrices which keep the topological information
+void  StochDMNDG3DMaterialLaw::fill_Matrices(const double VI){    
+  double Norm[9][3];
+  double W_WI[2][2];
+  int pos0=0;
+  int pos1=0; 
+  int k, node;
+  int Col_start, Col_start0,Col_start1, Row_start,Row_start0,  Row_start1;
+  double tmp;
+  int col = 9;    
+      
+  for(int i=0; i <_N_Node; i++){
+    _VA.push_back(1.0);  
+  }    
+  // fill matrix for equilibrium of interface
+  fullMatrix<double> _NT(3*_NInterface, col*(_N_Node-1)); 
+  _VI.push_back(VI);
+  int l = -1;
+  for(int pos=0; pos <_NInterface; pos++){  
+    if(pos == int(pow(2,l+1)-1)) l +=1;
+    fill_NLocal(pos,Norm);
+    pos_kids(pos, l, pos0, pos1);
+    Row_start = pos*3;
+    Col_start0 = (pos0-1)*col;
+    Col_start1 = (pos1-1)*col;
+    
+    double S_ba = 1.0;    
+    if(l ==_TotalLevel-1 and _porous) S_ba = _Para_Wt[pos][1];
+    // information of interface norm
+    for(int nc=0; nc<col;nc++){
+      for(int nr=0; nr<3;nr++){
+        _NT.set(Row_start+nr, Col_start0+nc, Norm[nc][nr]);
+        _NT.set(Row_start+nr, Col_start1+nc,-S_ba*Norm[nc][nr]);
+      }
+    } 
+    
+    if(l ==_TotalLevel-1){
+      if(_porous){
+        _VA[pos0] = (1.0-_VI[pos])*_Para_Wt[pos][0];
+        _VA[pos1] = (1.0-_VI[pos])*(1.0-_Para_Wt[pos][0]);   
+      }
+      else{      
+        _VA[pos0] = 1.0-_VI[pos];
+        _VA[pos1] = _VI[pos];
+      }  
+    }
+    else if(_Modified_DMN and l ==_TotalLevel-2){
+      fill_W_WI(pos,l,W_WI);
+      _VA[pos0] = W_WI[0][0];
+      _VA[pos1] = W_WI[0][1];
+      _VI.push_back(W_WI[1][0]);
+    }
+    else{
+      fill_W_WI(pos,l,W_WI);
+      _VA[pos0] = W_WI[0][0];
+      _VA[pos1] = W_WI[0][1];
+      _VI.push_back(W_WI[1][0]);
+      _VI.push_back(W_WI[1][1]);
+    }
+    //information of homogenous strain to local strain 
+    Col_start = pos*3;
+    if(!_Modified_DMN or l< _TotalLevel-2){
+      k = int(pow(2,l));
+      node = _NRoot/k/2;
+      int i = pos-int(pow(2,l))+1;
+      for(int r=0; r<node; r++){
+        Row_start0 = (i*2*node+r)*col;
+        Row_start1 = (i*2*node+r+node)*col;
+        for(int nr=0; nr<col;nr++){
+          for(int nc=0; nc<3;nc++){
+            _Nv.set(Row_start0+nr, Col_start+nc,Norm[nr][nc]/_VA[pos0]);
+            _Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/_VA[pos1]);
+          }
+        }
+      }   
+    }
+    else if(l == _TotalLevel-2){
+      node = 3;
+      int i = pos-int(pow(2,l))+1;
+      for(int r=0; r<2; r++){
+        Row_start0 = (i*node+r)*col;
+        for(int nr=0; nr<col;nr++){
+          for(int nc=0; nc<3;nc++){
+            _Nv.set(Row_start0+nr, Col_start+nc,Norm[nr][nc]/_VA[pos0]);
+          }
+        }
+      }  
+      Row_start1 = (i*node+2)*col;
+      for(int nr=0; nr<col;nr++){
+        for(int nc=0; nc<3;nc++){
+          _Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/_VA[pos1]);
+        }
+      }
+    }  
+    else if(l == _TotalLevel-1){
+      node = 3;
+      int i = pos-int(pow(2,l))+1;
+      Row_start0 = (i*node)*col;
+      Row_start1 = (i*node+1)*col;
+      for(int nr=0; nr<col;nr++){
+        for(int nc=0; nc<3;nc++){
+          _Nv.set(Row_start0+nr, Col_start+nc,Norm[nr][nc]/_VA[pos0]);
+          _Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/_VA[pos1]);
+        }
+      }
+    }  
+   }
+   
+ // information of voulme fraction to stress on node   
+  Row_start = col*(_N_Node-1-_NRoot);  
+  for(int nc=0; nc<col*_NRoot; nc++){
+    _V.set(Row_start+nc,nc,1.0);
+  }  
+    
+  l = _TotalLevel-1;
+  for(int pos = _NInterface-1; pos > 0; pos--){  
+    if(pos == int(pow(2,l)-2)) l -=1; 
+    pos_kids(pos, l, pos0, pos1);
+    Row_start = (pos-1)*col;
+    Row_start0 = (pos0-1)*col;
+    Row_start1 = (pos1-1)*col;
+    for(int nr=0; nr<col;nr++){
+      for(int nc=0; nc<col*_NRoot ; nc++){ 
+        tmp = _VA[pos0]*_V(Row_start0+nr,nc) + _VA[pos1]*_V(Row_start1+nr,nc);
+        _V.set(Row_start+nr,nc, tmp);
+      }  
+    }  
+  }
+  _NT.mult(_V, _NTV);
+      
+  for(int i=0; i<col*_NRoot;i++){
+    _I(i, i%col) = 1.0;
+  }   
+}     
+ // end fill_matrices     
+         
+// read parameter from input file  
+void StochDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
+    
+  double sin0, cos0, sin1,cos1; 
+  int Model; 
+  FILE *Para = fopen(ParaFile, "r");
+  if ( Para != NULL ){ 
+    int okf = fscanf(Para, "%d\n", &Model);
+    _Modified_DMN = static_cast<bool>(Model);
+    
+    okf = fscanf(Para, "%d %d\n", &_Dim, &_TotalLevel);  
+    
+    bool resizeFlag;
+    if(_Modified_DMN){
+       _NRoot = int(pow(2,(_TotalLevel-2))*3);   
+    }
+    else{    
+      _NRoot = int(pow(2,_TotalLevel)); 
+    }  
+    _NInterface = _NRoot-1;
+    _N_Node = _NRoot + _NInterface;
+      
+    int row = 9*_NRoot;
+    int col = 3*_NInterface;
+    resizeFlag = _Nv.resize(row, col, true);
+        
+    resizeFlag = _I.resize(row, 9, true);
+        
+    row = 9*(_N_Node-1);
+    col = 9*_NRoot;
+    resizeFlag = _V.resize(row, col, true);
+        
+    row = 3*_NInterface;        
+    resizeFlag = _NTV.resize(row, col, true);             
+        
+    okf = fscanf(Para, "%lf\n", &_Vf);
+        
+    _NPara_Norm = _NInterface;
+        
+    double ParaN[2];
+    double TwoPi(2.0*3.14159265359);        
+    SPoint3 Para_Norm;
+    SPoint2 Para_Wt;
+        
+    if(_Dim==2){
+      ParaN[1] = 0.0;
+      sin1 = sin(TwoPi*ParaN[1]);
+      cos1 = cos(TwoPi*ParaN[1]);
+      for(int i=0;i<_NInterface;i++){          
+        okf = fscanf(Para, "%lf\n", &ParaN[0]);
+        sin0 = sin(TwoPi*ParaN[0]);
+        cos0 = cos(TwoPi*ParaN[0]);
+        Para_Norm[0] = cos1*sin0;
+        Para_Norm[1] = cos1*cos0;
+        Para_Norm[2] = sin1;
+        _Para_Norm.push_back(Para_Norm); 
+      }  
+    }
+        
+    else if(_Dim==3){          
+      for(int i=0;i<_NInterface;i++){          
+        okf = fscanf(Para, "%lf %lf\n", &ParaN[0],  &ParaN[1]);
+        sin0 = sin(TwoPi*ParaN[0]);
+        cos0 = cos(TwoPi*ParaN[0]);
+        sin1 = sin(TwoPi*ParaN[1]);
+        cos1 = cos(TwoPi*ParaN[1]);
+        Para_Norm[0] = cos1*sin0;
+        Para_Norm[1] = cos1*cos0;
+        Para_Norm[2] = sin1;
+        _Para_Norm.push_back(Para_Norm); 
+      }  
+    }       
+        
+    if(_porous){
+      _NPara_Wt = _NInterface;
+    }  
+    else{  
+      _NPara_Wt = int(pow(2,(_TotalLevel-1)))-1;
+    }
+    for(int i=0;i<_NPara_Wt;i++){          
+      okf = fscanf(Para, "%lf %lf\n", &Para_Wt[0],  &Para_Wt[1]);
+      _Para_Wt.push_back(Para_Wt); 
+    }
+  }         
+  Msg::Info("Finish reading parameter file, Model = %d, level = %d, Dim = %d, NRoot = %d",Model, _TotalLevel, _Dim, _NRoot);
+}
+   // end of read_parameter     
+        
+       
+void StochDMNDG3DMaterialLaw::addLaw(dG3DMaterialLaw* law){
+  _mapLaw.push_back(law);
+};  
+      
+      
+void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff, const bool checkfrac, const bool dTangent){
+   /* get ipvariable */
+  StochDMNDG3DIPVariable* ipvcur = static_cast<StochDMNDG3DIPVariable*>(ipv);;
+  const StochDMNDG3DIPVariable* ipv_prev = static_cast<const StochDMNDG3DIPVariable*>(ipvprev);
+  const STensor3& F = ipvcur->getConstRefToDeformationGradient(); 
+    
+  STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
+  STensor43& L = ipvcur->getRefToTangentModuli();
+  STensorOperation::zero(P);
+  STensorOperation::zero(L);
+   
+  static double C_hom[9][9]; 
+  static double P_hom[9];
+    
+  fullVector<double>& Fvec = ipvcur->getRefToDeformationGradientVect(); 
+  fullVector<double>& Pvec = ipvcur->getRefToFirstPiolaKirchhoffStressVect();
+  fullVector<double>& a_cur = ipvcur->getRefToInterfaceFluctuationVect();
+  Pvec.setAll(0.0);  
+ //  fullVector<double> P_node(9*(_N_Node-1));
+  static fullVector<double> Res_node_t(3*_NInterface);
+  static fullVector<double> Res_node(3*_NInterface);
+  static fullVector<double> a_step(3*_NInterface);
+  static fullMatrix<double> C(9*_NRoot,9*_NRoot);
+  static fullMatrix<double> NTVC(3*_NInterface, 9*_NRoot);
+  static fullMatrix<double> Jacobian_inv(3*_NInterface, 3*_NInterface); 
+  static fullMatrix<double> Jacobian(3*_NInterface, 3*_NInterface);  
+  static fullMatrix<double> Modify_Jacobian(3*_NInterface, 3*_NInterface); 
+  static fullMatrix<double> tmp(9*_NRoot, 9); 
+  static fullMatrix<double> dRdF(3*_NInterface, 9); 
+  static fullMatrix<double> dadF(3*_NInterface, 9); 
+    
+  Res_node_t.setAll(0.0);
+  Res_node.setAll(0.0);
+  a_step.setAll(0.0);
+  C.setAll(0.0);
+  NTVC.setAll(0.0); 
+  Jacobian.setAll(0.0);
+  dRdF.setAll(0.0);
+  dadF.setAll(0.0);
+  tmp.setAll(0.0);
+  int ite = 0;
+  double r = 1.0;
+  bool last = false;
+  int pos_start = _N_Node-_NRoot;
+  int max_iter= 15;
+  bool stiff_loc = true;
+ 
+  
+  while( (r>_tol and ite < max_iter) or last){  
+  //  if((1+ite)%5 == 0) stiff_loc=false;
+    ite +=1;
+    _Nv.mult(a_cur, Fvec);
+    for(int i=0; i<_NRoot; i++){       
+      if (_VA[pos_start+i] > 1e-6){     
+        dG3DIPVariableBase* ipv_i = dynamic_cast<dG3DIPVariableBase*>(ipvcur->getIPv(i));
+        const dG3DIPVariableBase* ipvprev_i = dynamic_cast<const dG3DIPVariableBase*>(ipv_prev->getIPv(i));
+        
+        STensor3& Floc = ipv_i->getRefToDeformationGradient();   
+        for(int m=0; m<3; m++){
+          for(int n=0; n<3; n++){  
+            Floc(m,n) = F(m,n) + Fvec(i*9 + m*3+n);
+          } 
+        }
+        int mat_IP = 0;
+        if(_Modified_DMN){
+          if(i%3 == 1) mat_IP = 1;  
+        }   
+        else{
+           mat_IP = i%2;  
+        }
+                
+        if(_porous) mat_IP=0;        
+        _mapLaw[mat_IP]->stress(ipv_i, ipvprev_i, stiff_loc, checkfrac, dTangent);
+        const STensor3& P_i = ipv_i->getConstRefToFirstPiolaKirchhoffStress();
+        
+        for(int m=0; m<3; m++){
+            for(int n=0; n<3; n++){     
+              Pvec(i*9 + m*3+n) = P_i(m,n);
+            }
+        }      
+        
+        if(stiff_loc){
+          const STensor43& L_i = ipv_i->getConstRefToTangentModuli();        
+          for(int m=0; m<3; m++){
+            for(int n=0; n<3; n++){                   
+              for(int k=0; k<3; k++){
+                for(int l=0; l<3; l++){  
+                  C(i*9+m*3+n, i*9+k*3+l) = L_i(m,n,k,l);
+                }
+              }
+            }
+          } 
+        }    
+      }     
+    }   
+    
+    _NTV.mult(Pvec, Res_node);        
+    r = Res_node.norm()/_NInterface;
+    
+    if(r<_tol){
+      if(last or stiff_loc){
+        _NTV.mult(C, NTVC);
+        NTVC.mult(_Nv, Jacobian);
+        Jacobian.invertInPlace();
+        break;
+      }
+      else{
+        last = true;
+        stiff_loc = true;
+      }  
+    } 
+            
+    else{         
+      if(stiff_loc){         
+        _NTV.mult(C, NTVC);
+        NTVC.mult(_Nv, Jacobian);
+        Jacobian.invertInPlace();
+      }
+      else{
+        _NTV.mult(Pvec, Res_node);
+        Res_node_t.scale(-1.0);
+        Res_node_t.axpy(Res_node, 1.0);
+            
+        double Rho = 1.0/(Res_node_t*a_step);
+        Modify_Jacobian.setAll(0.0);
+        for(int i =0; i<3*_NInterface; i++){
+          Modify_Jacobian(i,i) = 1.0;
+          for(int j =0; j<3*_NInterface; j++){
+            Modify_Jacobian(i,j) -= a_step(i)*Res_node_t(j)*Rho;
+          }
+        }
+        Modify_Jacobian.mult(Jacobian,Jacobian_inv);
+        Modify_Jacobian.transposeInPlace();
+        Jacobian_inv.mult(Modify_Jacobian,Jacobian);
+      
+        for(int i =0; i<3*_NInterface; i++){
+          for(int j =0; j<3*_NInterface; j++){
+            Jacobian(i,j) += a_step(i)*a_step(j)*Rho;
+          }
+        } 
+        stiff_loc = true;   
+      }
+       
+      Res_node_t = Res_node; 
+      Jacobian.mult(Res_node, a_step);
+      a_step.scale(-1.0);
+      a_cur.axpy(a_step, 1.0); 
+    }  
+  }  
+  Jacobian.scale(-1.0);  
+    
+  if(r>_tol){  
+    Msg::Error("DMN residual (= %lf) did n't converge to tolenerce after %d iteration",r,ite); 
+    return;  
+  }
+  else{
+    for(int m=0;m<9;m++){
+      P_hom[m] = 0.0;
+      for(int i=0;i<_NRoot;i++){
+        P_hom[m] += (_VA[1]*_V(m, 9*i+m) + _VA[2]*_V(9+m, 9*i+m)) *Pvec(9*i+m);
+      }
+    }
+    for(int i=0;i<3;i++){
+      for(int j=0;j<3;j++){
+        P(i,j) = P_hom[3*i+j];
+      }
+    }    
+      
+    if(stiff){    
+      NTVC.mult(_I, dRdF);
+      Jacobian.mult(dRdF, dadF);
+      _Nv.mult(dadF, tmp);
+      tmp.add(_I);
+          
+      for(int m=0;m<9;m++){
+        for(int n=0;n<9;n++){
+          C_hom[m][n] = 0.0;
+          for(int i=0;i<_NRoot;i++){
+            for(int j=0;j< 9*_NRoot;j++){
+              C_hom[m][n] += (_VA[1]*_V(m,m+9*i) + _VA[2]*_V(9+m,m+9*i))*C(m+9*i,j)*tmp(j,n);
+            }
+          }
+        }
+      }                     
+        
+      for(int i=0;i<3;i++){
+        for(int j=0;j<3;j++){
+          for(int k=0;k<3;k++){
+            for(int l=0;l<3;l++){  
+              L(i,j,k,l) = C_hom[i*3+j][k*3+l];
+            }
+          }
+        }
+      }        
+    }
+  }
+  ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+} 
+        
+  // end of stress      
+        
+      
+void StochDMNDG3DMaterialLaw::setTime(const double t,const double dtime){
+  dG3DMaterialLaw::setTime(t,dtime);
+  for (std::vector<dG3DMaterialLaw*>::iterator it = _mapLaw.begin(); it != _mapLaw.end(); it++) {
+    (*it)->setTime(t,dtime);
+  }
+}
+
+/* back up for original 2 phases interaction DMN
+
+StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile,const bool porous,const double tol):
+  dG3DMaterialLaw(num,rho,false),_porous(porous),_tol(tol){
+          
+  read_parameter(ParaFile);
+  fill_Matrices(_Vf);         
+}
+    
+StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const StochDMNDG3DMaterialLaw &src): dG3DMaterialLaw(src), _porous(src._porous), _TotalLevel(src._TotalLevel), _NRoot(src._NRoot),
+       _NInterface(src._NInterface), _N_Node(src._N_Node), _NPara_Wt(src._NPara_Wt), _NPara_Norm(src._NPara_Norm), _Dim(src._Dim), _Vf(src._Vf), _tol(src._tol),
+       _Nv(src._Nv), _NTV(src._NTV), _V(src._V), _Para_Wt(src._Para_Wt), _Para_Norm(src._Para_Norm), _VA(src._VA), _VI(src._VI), _mapLaw(src._mapLaw){ }
+    
+StochDMNDG3DMaterialLaw::~StochDMNDG3DMaterialLaw(){};
+     
+          
+void StochDMNDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const{
+      // check interface element
+  bool inter=true;
+  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+  if(iele==NULL) inter=false;
+    IPVariable* ipvi = new  StochDMNDG3DIPVariable(_NRoot, hasBodyForce,inter);
+    IPVariable* ipv1 = new  StochDMNDG3DIPVariable(_NRoot, hasBodyForce,inter);
+    IPVariable* ipv2 = new  StochDMNDG3DIPVariable(_NRoot, hasBodyForce,inter);
+       
+    StochDMNDG3DIPVariable* ipvi_root = static_cast<StochDMNDG3DIPVariable*>(ipvi);
+    StochDMNDG3DIPVariable* ipv1_root = static_cast<StochDMNDG3DIPVariable*>(ipv1);
+    StochDMNDG3DIPVariable* ipv2_root = static_cast<StochDMNDG3DIPVariable*>(ipv2);
+    for (int i=0; i< _NRoot; i++){
+      IPStateBase* ips_i = NULL;
+      int mat_IP = i%2;
+      if(_porous) mat_IP=0; 
+      _mapLaw[mat_IP]->createIPState(ips_i, hasBodyForce, state_, ele, nbFF_, GP, gpt);
+      std::vector<IPVariable*> ip_all;
+      ips_i->getAllIPVariable(ip_all);
+      ipvi_root->addIPv(i, ip_all[0]);
+      ipv1_root->addIPv(i, ip_all[1]);
+      ipv2_root->addIPv(i, ip_all[2]);
+    } 
+    if(ips != NULL) delete ips;
+      ips = new IP3State(state_,ipvi,ipv1,ipv2);
+}
+       
+   
+void StochDMNDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const {
+  if(ipv !=NULL) delete ipv;
+  bool inter=true;
+  const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+  if(iele == NULL){
+    inter=false;
+  }
+  ipv = new  StochDMNDG3DIPVariable(_NRoot, hasBodyForce,inter);
+  StochDMNDG3DIPVariable* ipv_root = static_cast<StochDMNDG3DIPVariable*>(ipv);
+  for (int i=0; i< _NRoot; i++){
+    IPVariable* ipv_i = NULL;
+    int mat_IP = i%2;
+    if(_porous) mat_IP=0; 
+    _mapLaw[mat_IP]->createIPVariable(ipv_i, hasBodyForce, ele, nbFF_, GP, gpt);
+    ipv_root->addIPv(i, ipv_i);
+  }      
+};
+    
+void StochDMNDG3DMaterialLaw::initialIPVariable(IPVariable* ipv, bool stiff){
+  StochDMNDG3DIPVariable* ipv_all = dynamic_cast<StochDMNDG3DIPVariable*>(ipv);
+  ipv_all->InitializeFluctuationVector();
+  if (ipv_all == NULL){
+    Msg::Error("StochDMNDG3DIPVariable must be used in StochDMNDG3DMaterialLaw::initialIPVariable");
+  }
+  for (int i=0; i< _NRoot; i++){
+    IPVariable* ipv_i = ipv_all->getIPv(i);        
+    int mat_IP = i%2;
+    if(_porous) mat_IP=0; 
+    _mapLaw[mat_IP]->initialIPVariable(ipv_i, stiff);
+  }
+}  
+        
+void StochDMNDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{
+  StochDMNDG3DIPVariable* ipv_all = dynamic_cast<StochDMNDG3DIPVariable*>(ipv);
+  const StochDMNDG3DIPVariable* ipvprev_all = dynamic_cast<const StochDMNDG3DIPVariable*>(ipvprev);
+  if (ipv_all == NULL){
+    Msg::Error("StochDMNDG3DIPVariable must be used in StochDMNDG3DMaterialLaw::checkInternalState");
+  }
+  for (int i=0; i< _NRoot; i++){
+    IPVariable* ipv_i = ipv_all->getIPv(i);
+    const IPVariable* ipvprev_i = ipvprev_all->getIPv(i);
+    int mat_IP = i%2;
+    if(_porous) mat_IP=0; 
+    _mapLaw[mat_IP]->checkInternalState(ipv_i, ipvprev_i);
+  }
+}
+     
+void StochDMNDG3DMaterialLaw::TwoPhasesInteract(const STensor43& C0, const STensor43& C1, fullMatrix<double> &mat, const int pos, const int level){
+  double n0 = _Para_Norm[pos][0];
+  double n1 = _Para_Norm[pos][1];
+  double n2 = _Para_Norm[pos][2];
+  double tmp1[6][6], tmp2[6][6], Norm[6][3], NormT[3][6];
+  double S_ba = 1.0;
+  if(level ==_TotalLevel-1 and _porous){
+    S_ba = _Para_Wt[pos][1];
+  }
+       
+  for(int i=0; i<3; i++){
+    for(int j=0; j<3; j++){
+      Norm[i][j] = 0.0;
+      Norm[3+i][j] = 0.0;
+      NormT[i][j] = 0.0;
+      NormT[i][3+j] = 0.0;
+    }
+  }        
+  Norm[0][0] = n0;        // different N and NT are used in order to be consistent with function "fromSTensor43ToFullMatrix()"
+  Norm[1][1] = n1;
+  Norm[2][2] = n2;
+  Norm[3][0] = n1/2.0;
+  Norm[3][1] = n0/2.0;
+  Norm[4][0] = n2/2.0;
+  Norm[4][2] = n0/2.0;
+  Norm[5][1] = n2/2.0;
+  Norm[5][2] = n1/2.0;  
+       
+  NormT[0][0] = n0;
+  NormT[1][1] = n1;
+  NormT[2][2] = n2;
+  NormT[0][3] = n1;
+  NormT[0][4] = n2;
+  NormT[1][3] = n0;
+  NormT[1][5] = n2;
+  NormT[2][4] = n0;       
+  NormT[2][5] = n1;    
+       
+  fullMatrix<double> mat0(6,6), mat1(6,6), K(3,3);       
+  int pos0 = 2*pos+1;
+  int pos1 = 2*pos+2;
+       
+  if(level == _TotalLevel-1){
+    STensorOperation::fromSTensor43ToFullMatrix(C0, mat0);
+    STensorOperation::fromSTensor43ToFullMatrix(C1, mat1);       
+  }
+  else{
+    TwoPhasesInteract(C0, C1, mat0, pos0, level+1);
+    TwoPhasesInteract(C0, C1, mat1, pos1, level+1);       
+  }
+       
+  double v_A = _VA[pos0];
+  double v_B = _VA[pos1];
+  for(int i=0; i<6;i++){
+    for(int j=0; j<6;j++){
+      tmp1[i][j] = mat0(i,j)/v_A + S_ba*S_ba*mat1(i,j)/v_B;      
+      tmp2[i][j] = S_ba*mat1(i,j) - mat0(i,j);      
+    }  
+  }
+       
+  for(int m=0; m<3;m++){
+    for(int n=0; n<3;n++) {
+      for(int i=0; i<6;i++) {
+        for(int j=0; j<6;j++) { 
+          K(m,n) += NormT[m][i]*tmp1[i][j]*Norm[j][n];
+        }
+      }
+    }
+  }         
+  K.invertInPlace();
+  for(int m=0; m<6;m++){
+    for(int n=0; n<6;n++) {
+      tmp1[m][n] = 0.0;
+      for(int i=0; i<3;i++) {
+        for(int j=0; j<3;j++) { 
+          tmp1[m][n] += Norm[m][i]*K(i,j)*NormT[j][n];
+        }
+      }
+    }
+  }
+       
+  for(int m=0; m<6;m++){
+    for(int n=0; n<6; n++) {
+      mat(m,n) = mat0(m,n)*v_A + mat1(m,n)*v_B;
+      for(int i=0; i<6;i++) {
+        for(int j=0; j<6;j++) { 
+          mat(m,n) -= tmp2[m][i]*tmp1[i][j]*tmp2[j][n];        
+        }  
+      }
+    }
+  }        
+}
+  // end of TwoPhasesInteract    
+      
+void StochDMNDG3DMaterialLaw::initLaws(const std::map<int,materialLaw*> &maplaw){
+  if (!_initialized){
+    STensorOperation::zero(elasticStiffness);
+    fullMatrix<double> mat(6,6);
+    std::vector<STensor43> C;    
+    for (std::vector<dG3DMaterialLaw*>::iterator ite = _mapLaw.begin(); ite!= _mapLaw.end(); ite++){
+      (*ite)->initLaws(maplaw);
+      elasticStiffness = (*ite)->elasticStiffness;
+      C.push_back(elasticStiffness);
+    }
+    if(_porous) C.push_back(elasticStiffness);
+    TwoPhasesInteract(C[0], C[1],  mat, 0, 0);
+    
+    STensorOperation::fromFullMatrixToSTensor43(mat, elasticStiffness);
+   // elasticStiffness.print("elasticStiffness StochDMNDG3DMaterialLaw");
+    _initialized = true;
+   // mat.print("initialize material tensor");
+  }  
+}  
+
+void StochDMNDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
+   dG3DMaterialLaw::setMacroSolver(sv);
+    for (std::vector<dG3DMaterialLaw*>::iterator ite = _mapLaw.begin(); ite!= _mapLaw.end(); ite++){
+      (*ite)->setMacroSolver(sv);
+    }
+}
+         
+void StochDMNDG3DMaterialLaw::fill_NLocal(const int pos, double N[][3]){
+  double n0 = _Para_Norm[pos][0];
+  double n1 = _Para_Norm[pos][1];
+  double n2 = _Para_Norm[pos][2];
+   
+  for(int i=0; i<9; i++){
+    for(int j=0; j<3; j++){
+      N[i][j]=0.0;
+    }
+  }    
+  N[0][0]= n0;
+  N[1][0]= n1;
+  N[2][0]= n2;
+  N[3][1]= n0;
+  N[4][1]= n1;
+  N[5][1]= n2;
+  N[6][2]= n0;
+  N[7][2]= n1;
+  N[8][2]= n2;
+}   
+
+/*void StochDMNDG3DMaterialLaw::fill_W_WI(const int pos,  double W[][2]){
+  double VI = _VI[pos];
+  double WA = _Para_Wt[pos][0];
+    
+  for(int i=0; i<2; i++){
+    for(int j=0; j<2; j++){
+      W[i][j] = 0.0;
+    }
+  }      
+  W[0][0]= WA;
+  W[0][1]= 1.0-WA;   
+}  
+
+void  StochDMNDG3DMaterialLaw::fill_Matrices(const double VI){    
+  double Norm[9][3];
+  double W_WI[2][2];
+  int pos, pos0, pos1; 
+  int k, node;
+  int Col_start, Col_start0,Col_start1, Row_start,Row_start0,  Row_start1;
+  double tmp;
+  int col = 9;
+      
+  fullMatrix<double> _NT(3*_NInterface, col*(_N_Node-1));     
+  _VA.push_back(1.0);
+  for(int l=0; l<_TotalLevel; l++){
+    k = int(pow(2,l));
+    node = _NRoot/k/2;
+    double S_ba = 1.0;    
+    for(int i=0; i<k ; i++){
+      pos = int(pow(2,l))-1+i;
+      if(l ==_TotalLevel-1 and _porous){
+        S_ba = _Para_Wt[pos][1];
+      }
+      fill_NLocal(pos,Norm);
+      pos0 = 2*pos+1;
+      pos1 = 2*pos+2;
+      Row_start = pos*3;
+      Col_start0 = (pos0-1)*col;
+      Col_start1 = (pos1-1)*col;
+      
+      // information of interface norm
+      for(int nc=0; nc<col;nc++){
+        for(int nr=0; nr<3;nr++){
+          _NT.set(Row_start+nr, Col_start0+nc,Norm[nc][nr]);
+          _NT.set(Row_start+nr, Col_start1+nc,-S_ba*Norm[nc][nr]);
+        }
+      }          
+          
+      if(l ==_TotalLevel-1 and _porous){
+          W_WI[0][0] = (1.0-_VI[pos])*_Para_Wt[pos][0];
+          W_WI[0][1] = (1.0-_VI[pos])*(1.0-_Para_Wt[pos][0]);
+          _VA.push_back(W_WI[0][0]);
+          _VA.push_back(W_WI[0][1]);         
+      }
+      else{
+        fill_W_WI(pos,W_WI);
+        _VA.push_back(W_WI[0][0]);
+        _VA.push_back(W_WI[0][1]);
+      }
+        
+      //information of homogenous strain to local strain 
+      Col_start = pos*3;
+      for(int r=0; r<node; r++){
+        Row_start0 = (i*2*node+r)*col;
+        Row_start1 = (i*2*node+r+node)*col;
+        for(int nr=0; nr<col;nr++){
+          for(int nc=0; nc<3;nc++){
+            _Nv.set(Row_start0+nr, Col_start+nc,Norm[nr][nc]/W_WI[0][0]);
+            _Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/W_WI[0][1]);
+          }
+        }   
+      }  
+    }
+  }  
+      
+  // information of voulme fraction 
+  Row_start = col*(_N_Node-1-_NRoot);  
+  for(int nr=0; nr<col*_NRoot; nr++){
+    _V.set(Row_start+nr,nr,1.0);
+  }  
+      
+  for(int l=_TotalLevel-1; l>0; l--){
+    k = int(pow(2,l));
+    for(int i=0; i<k ; i++){
+      pos = k-1+i;
+      pos0 = 2*pos+1;
+      pos1 = 2*pos+2; 
+      Row_start = (pos-1)*col;
+      Row_start0 = (pos0-1)*col;
+      Row_start1 = (pos1-1)*col;
+      for(int nr=0; nr<col;nr++){
+        for(int nc=0; nc<col*_NRoot ; nc++){ 
+          tmp = _VA[pos0]*_V(Row_start0+nr,nc) + _VA[pos1]*_V(Row_start1+nr,nc);
+          _V.set(Row_start+nr,nc, tmp);
+        }  
+      }
+    }
+  }
+  _NT.mult(_V, _NTV);
+      
+  for(int i=0; i<col*_NRoot;i++){
+    _I(i, i%col) = 1.0;
+  }       
+}     */
+
+
+/*void StochDMNDG3DMaterialLaw::fill_W_WI(const int pos,  double W[][2]){
+  double VI = _VI[pos];
+  double Wtmp0 = _Para_Wt[pos][0];
+  double Wtmp1 = _Para_Wt[pos][1];
+  double WA = Wtmp0*(1.0-VI)+ Wtmp1*VI;
+  
+  for(int i=0; i<2; i++){
+    for(int j=0; j<2; j++){
+      W[i][j] = 0.0;
+    }
+  }      
+  W[0][0]= WA;
+  W[0][1]= 1.0-WA;
+  W[1][0]= (Wtmp1*VI)/WA;
+  W[1][1]= ((1.0-Wtmp1)*VI)/(1.0-WA);      
+} 
+
+     
+// fill the matrices which keep the topological information
+void  StochDMNDG3DMaterialLaw::fill_Matrices(const double VI){    
+  double Norm[9][3];
+  double W_WI[2][2];
+  int pos, pos0, pos1; 
+  int k, node;
+  int Col_start, Col_start0,Col_start1, Row_start,Row_start0,  Row_start1;
+  double tmp;
+  int col = 9;
+      
+  fullMatrix<double> _NT(3*_NInterface, col*(_N_Node-1));     
+  _VI.push_back(VI);
+  _VA.push_back(1.0);
+  for(int l=0; l<_TotalLevel; l++){
+    k = int(pow(2,l));
+    node = _NRoot/k/2;
+    double S_ba = 1.0;    
+    for(int i=0; i<k ; i++){
+      pos = int(pow(2,l))-1+i;
+      if(l ==_TotalLevel-1 and _porous){
+        S_ba = _Para_Wt[pos][1];
+      }
+      fill_NLocal(pos,Norm);
+      pos0 = 2*pos+1;
+      pos1 = 2*pos+2;
+      Row_start = pos*3;
+      Col_start0 = (pos0-1)*col;
+      Col_start1 = (pos1-1)*col;
+      
+      // information of interface norm
+      for(int nc=0; nc<col;nc++){
+        for(int nr=0; nr<3;nr++){
+          _NT.set(Row_start+nr, Col_start0+nc,Norm[nc][nr]);
+          _NT.set(Row_start+nr, Col_start1+nc,-S_ba*Norm[nc][nr]);
+        }
+      }          
+          
+      if(l ==_TotalLevel-1){
+        if(_porous){
+          W_WI[0][0] = (1.0-_VI[pos])*_Para_Wt[pos][0];
+          W_WI[0][1] = (1.0-_VI[pos])*(1.0-_Para_Wt[pos][0]);
+          _VA.push_back(W_WI[0][0]);
+          _VA.push_back(W_WI[0][1]);     
+        }
+        else{      
+          W_WI[0][0] = 1.0-_VI[pos];
+          W_WI[0][1] = _VI[pos];
+          _VA.push_back(W_WI[0][0]);
+          _VA.push_back(W_WI[0][1]);
+        }  
+      }
+      else{
+        fill_W_WI(pos,W_WI);
+        _VA.push_back(W_WI[0][0]);
+        _VA.push_back(W_WI[0][1]);
+        _VI.push_back(W_WI[1][0]);
+        _VI.push_back(W_WI[1][1]);
+      }
+        
+      //information of homogenous strain to local strain 
+      Col_start = pos*3;
+      for(int r=0; r<node; r++){
+        Row_start0 = (i*2*node+r)*col;
+        Row_start1 = (i*2*node+r+node)*col;
+        for(int nr=0; nr<col;nr++){
+          for(int nc=0; nc<3;nc++){
+            _Nv.set(Row_start0+nr, Col_start+nc,Norm[nr][nc]/W_WI[0][0]);
+            _Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm[nr][nc]/W_WI[0][1]);
+          }
+        }   
+      }  
+    }
+  }  
+      
+  // information of voulme fraction 
+  Row_start = col*(_N_Node-1-_NRoot);  
+  for(int nr=0; nr<col*_NRoot; nr++){
+    _V.set(Row_start+nr,nr,1.0);
+  }  
+      
+  for(int l=_TotalLevel-1; l>0; l--){
+    k = int(pow(2,l));
+    for(int i=0; i<k ; i++){
+      pos = k-1+i;
+      pos0 = 2*pos+1;
+      pos1 = 2*pos+2; 
+      Row_start = (pos-1)*col;
+      Row_start0 = (pos0-1)*col;
+      Row_start1 = (pos1-1)*col;
+      for(int nr=0; nr<col;nr++){
+        for(int nc=0; nc<col*_NRoot ; nc++){ 
+          tmp = _VA[pos0]*_V(Row_start0+nr,nc) + _VA[pos1]*_V(Row_start1+nr,nc);
+          _V.set(Row_start+nr,nc, tmp);
+        }  
+      }
+    }
+  }
+  _NT.mult(_V, _NTV);
+      
+  for(int i=0; i<col*_NRoot;i++){
+    _I(i, i%col) = 1.0;
+  }       
+}     
+ // end fill_matrices     
+         
+// read parameter from input file  
+void StochDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
+    
+  double sin0, cos0, sin1,cos1;  
+  FILE *Para = fopen(ParaFile, "r");
+  if ( Para != NULL ){ 
+    int okf = fscanf(Para, "%d %d\n", &_Dim, &_TotalLevel);
+       
+    bool resizeFlag;
+    _NRoot = int(pow(2,_TotalLevel)); 
+    _NInterface = _NRoot-1;
+    _N_Node = _NRoot + _NInterface;
+      
+    int row = 9*_NRoot;
+    int col = 3*_NInterface;
+    resizeFlag = _Nv.resize(row, col, true);
+        
+    resizeFlag = _I.resize(row, 9, true);
+        
+    row = 9*(_N_Node-1);
+    col = 9*_NRoot;
+    resizeFlag = _V.resize(row, col, true);
+        
+    row = 3*_NInterface;        
+    resizeFlag = _NTV.resize(row, col, true);             
+        
+    okf = fscanf(Para, "%lf\n", &_Vf);
+        
+    _NPara_Norm = _NInterface;
+        
+    double ParaN[2];
+    double TwoPi(2.0*3.14159265359);        
+    SPoint3 Para_Norm;
+    SPoint2 Para_Wt;
+        
+    if(_Dim==2){
+      ParaN[1] = 0.0;
+      sin1 = sin(TwoPi*ParaN[1]);
+      cos1 = cos(TwoPi*ParaN[1]);
+      for(int i=0;i<_NInterface;i++){          
+        okf = fscanf(Para, "%lf\n", &ParaN[0]);
+        sin0 = sin(TwoPi*ParaN[0]);
+        cos0 = cos(TwoPi*ParaN[0]);
+        Para_Norm[0] = cos1*sin0;
+        Para_Norm[1] = cos1*cos0;
+        Para_Norm[2] = sin1;
+        _Para_Norm.push_back(Para_Norm); 
+      }  
+    }
+        
+    else if(_Dim==3){          
+      for(int i=0;i<_NInterface;i++){          
+        okf = fscanf(Para, "%lf %lf\n", &ParaN[0],  &ParaN[1]);
+        sin0 = sin(TwoPi*ParaN[0]);
+        cos0 = cos(TwoPi*ParaN[0]);
+        sin1 = sin(TwoPi*ParaN[1]);
+        cos1 = cos(TwoPi*ParaN[1]);
+        Para_Norm[0] = cos1*sin0;
+        Para_Norm[1] = cos1*cos0;
+        Para_Norm[2] = sin1;
+        _Para_Norm.push_back(Para_Norm); 
+      }  
+    }       
+        
+    if(_porous){
+      _NPara_Wt = _NInterface;
+    }  
+    else{  
+      _NPara_Wt = int(pow(2,(_TotalLevel-1)))-1;
+    }
+    for(int i=0;i<_NPara_Wt;i++){          
+      okf = fscanf(Para, "%lf %lf\n", &Para_Wt[0],  &Para_Wt[1]);
+      _Para_Wt.push_back(Para_Wt); 
+    }
+  }  
+  
+  Msg::Info("Finish reading parameter file, level = %d, Dim = %d, NRoot = %d",_TotalLevel, _Dim, _NRoot);
+}
+   // end of read_parameter     
+        
+       
+void StochDMNDG3DMaterialLaw::addLaw(dG3DMaterialLaw* law){
+  _mapLaw.push_back(law);
+};  
+      
+      
+void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff, const bool checkfrac, const bool dTangent){
+   // get ipvariable 
+  StochDMNDG3DIPVariable* ipvcur = static_cast<StochDMNDG3DIPVariable*>(ipv);;
+  const StochDMNDG3DIPVariable* ipv_prev = static_cast<const StochDMNDG3DIPVariable*>(ipvprev);
+  const STensor3& F = ipvcur->getConstRefToDeformationGradient(); 
+    
+  STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress();
+  STensor43& L = ipvcur->getRefToTangentModuli();
+  STensorOperation::zero(P);
+  STensorOperation::zero(L);
+   
+  static double C_hom[9][9]; 
+  static double P_hom[9];
+    
+  fullVector<double>& Fvec = ipvcur->getRefToDeformationGradientVect(); 
+  fullVector<double>& Pvec = ipvcur->getRefToFirstPiolaKirchhoffStressVect();
+  fullVector<double>& a_cur = ipvcur->getRefToInterfaceFluctuationVect();
+  Pvec.setAll(0.0);
+    
+ //  fullVector<double> P_node(9*(_N_Node-1));
+  static fullVector<double> Res_node_t(3*_NInterface);
+  static fullVector<double> Res_node(3*_NInterface);
+  static fullVector<double> a_step(3*_NInterface);
+  static fullMatrix<double> C(9*_NRoot,9*_NRoot);
+  static fullMatrix<double> NTVC(3*_NInterface, 9*_NRoot);
+  static fullMatrix<double> Jacobian_inv(3*_NInterface, 3*_NInterface); 
+  static fullMatrix<double> Jacobian(3*_NInterface, 3*_NInterface);  
+  static fullMatrix<double> Modify_Jacobian(3*_NInterface, 3*_NInterface); 
+  static fullMatrix<double> tmp(9*_NRoot, 9); 
+  static fullMatrix<double> dRdF(3*_NInterface, 9); 
+  static fullMatrix<double> dadF(3*_NInterface, 9); 
+    
+  Res_node_t.setAll(0.0);
+  Res_node.setAll(0.0);
+  a_step.setAll(0.0);
+  C.setAll(0.0);
+  NTVC.setAll(0.0); 
+  Jacobian.setAll(0.0);
+  dRdF.setAll(0.0);
+  dadF.setAll(0.0);
+  tmp.setAll(0.0);
+  int ite = 0;
+  double r = 1.0;
+  bool last = false;
+  int pos_start = int(pow(2,_TotalLevel))-1;
+  int max_iter= 15;
+  bool stiff_loc = true;
+  
+  while( (r>_tol and ite < max_iter) or last){  
+  //  if((1+ite)%5 == 0) stiff_loc=false;
+    ite +=1;
+    _Nv.mult(a_cur, Fvec);
+    for(int i=0; i<_NRoot; i++){       
+      if (_VA[pos_start+i] > 1e-6){     
+        dG3DIPVariableBase* ipv_i = dynamic_cast<dG3DIPVariableBase*>(ipvcur->getIPv(i));
+        const dG3DIPVariableBase* ipvprev_i = dynamic_cast<const dG3DIPVariableBase*>(ipv_prev->getIPv(i));
+        
+        STensor3& Floc = ipv_i->getRefToDeformationGradient();   
+        for(int m=0; m<3; m++){
+          for(int n=0; n<3; n++){  
+            Floc(m,n) = F(m,n) + Fvec(i*9 + m*3+n);
+          } 
+        }
+        int mat_IP = i%2;
+        if(_porous) mat_IP=0;        
+        _mapLaw[mat_IP]->stress(ipv_i, ipvprev_i, stiff_loc, checkfrac, dTangent);
+        const STensor3& P_i = ipv_i->getConstRefToFirstPiolaKirchhoffStress();
+        
+        for(int m=0; m<3; m++){
+            for(int n=0; n<3; n++){     
+              Pvec(i*9 + m*3+n) = P_i(m,n);
+            }
+        }      
+        
+        if(stiff_loc){
+          const STensor43& L_i = ipv_i->getConstRefToTangentModuli();        
+          for(int m=0; m<3; m++){
+            for(int n=0; n<3; n++){                   
+              for(int k=0; k<3; k++){
+                for(int l=0; l<3; l++){  
+                  C(i*9+m*3+n, i*9+k*3+l) = L_i(m,n,k,l);
+                }
+              }
+            }
+          } 
+        }    
+      }     
+    }   
+    
+    _NTV.mult(Pvec, Res_node);        
+    r = Res_node.norm()/_NInterface;
+    
+    if(r<_tol){
+      if(last or stiff_loc){
+        _NTV.mult(C, NTVC);
+        NTVC.mult(_Nv, Jacobian);
+        Jacobian.invertInPlace();
+        break;
+      }
+      else{
+        last = true;
+        stiff_loc = true;
+      }  
+    } 
+            
+    else{         
+      if(stiff_loc){         
+        _NTV.mult(C, NTVC);
+        NTVC.mult(_Nv, Jacobian);
+        Jacobian.invertInPlace();
+      }
+      else{
+        _NTV.mult(Pvec, Res_node);
+        Res_node_t.scale(-1.0);
+        Res_node_t.axpy(Res_node, 1.0);
+            
+        double Rho = 1.0/(Res_node_t*a_step);
+        Modify_Jacobian.setAll(0.0);
+        for(int i =0; i<3*_NInterface; i++){
+          Modify_Jacobian(i,i) = 1.0;
+          for(int j =0; j<3*_NInterface; j++){
+            Modify_Jacobian(i,j) -= a_step(i)*Res_node_t(j)*Rho;
+          }
+        }
+        Modify_Jacobian.mult(Jacobian,Jacobian_inv);
+        Modify_Jacobian.transposeInPlace();
+        Jacobian_inv.mult(Modify_Jacobian,Jacobian);
+      
+        for(int i =0; i<3*_NInterface; i++){
+          for(int j =0; j<3*_NInterface; j++){
+            Jacobian(i,j) += a_step(i)*a_step(j)*Rho;
+          }
+        } 
+        stiff_loc = true;   
+      }
+       
+      Res_node_t = Res_node; 
+      Jacobian.mult(Res_node, a_step);
+      a_step.scale(-1.0);
+      a_cur.axpy(a_step, 1.0); 
+    }  
+  }  
+  Jacobian.scale(-1.0);  
+    
+  if(r>_tol){  
+    Msg::Error("DMN residual (= %lf) did n't converge to tolenerce after %d iteration",r,ite); 
+    return;  
+  }
+  else{
+    for(int m=0;m<9;m++){
+      P_hom[m] = 0.0;
+      for(int i=0;i<_NRoot;i++){
+        P_hom[m] += (_VA[1]*_V(m, 9*i+m) + _VA[2]*_V(9+m, 9*i+m)) *Pvec(9*i+m);
+      }
+    }
+    for(int i=0;i<3;i++){
+      for(int j=0;j<3;j++){
+        P(i,j) = P_hom[3*i+j];
+      }
+    }    
+      
+    if(stiff){    
+      NTVC.mult(_I, dRdF);
+      Jacobian.mult(dRdF, dadF);
+      _Nv.mult(dadF, tmp);
+      tmp.add(_I);
+          
+      for(int m=0;m<9;m++){
+        for(int n=0;n<9;n++){
+          C_hom[m][n] = 0.0;
+          for(int i=0;i<_NRoot;i++){
+            for(int j=0;j< 9*_NRoot;j++){
+              C_hom[m][n] += (_VA[1]*_V(m,m+9*i) + _VA[2]*_V(9+m,m+9*i))*C(m+9*i,j)*tmp(j,n);
+            }
+          }
+        }
+      }                     
+        
+      for(int i=0;i<3;i++){
+        for(int j=0;j<3;j++){
+          for(int k=0;k<3;k++){
+            for(int l=0;l<3;l++){  
+              L(i,j,k,l) = C_hom[i*3+j][k*3+l];
+            }
+          }
+        }
+      }        
+    }
+  }
+  ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+} 
+        
+  // end of stress      
+        
+      
+void StochDMNDG3DMaterialLaw::setTime(const double t,const double dtime){
+  dG3DMaterialLaw::setTime(t,dtime);
+  for (std::vector<dG3DMaterialLaw*>::iterator it = _mapLaw.begin(); it != _mapLaw.end(); it++) {
+    (*it)->setTime(t,dtime);
+  }
+}
+
+*/
+
+
+
+
+
+
 //
 J2SmallStrainDG3DMaterialLaw::J2SmallStrainDG3DMaterialLaw(const int num,const double rho,
                                    double E,const double nu,const double sy0,const double h,
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 9c805dd5dcec27e7224417f05bd3457e45bbe432..fb0421e33e6c4cc11add62bf1f7f1e99afab6ef7 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -716,6 +716,80 @@ class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{
 };
 
 
+class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
+  protected:
+    #ifndef SWIG
+    bool _porous;
+    bool _Modified_DMN;
+    int _TotalLevel;
+    int _NRoot;
+    int _NInterface;
+    int _N_Node;
+    int _NPara_Wt;
+    int _NPara_Norm;  
+    int _Dim;
+    double _Vf;
+    double _tol;
+    
+    std::vector<SPoint2> _Para_Wt;
+    std::vector<SPoint3> _Para_Norm;    
+    
+     
+    std::vector<dG3DMaterialLaw*> _mapLaw;
+    fullMatrix<double> _Nv;
+    fullMatrix<double> _NTV;
+    fullMatrix<double> _V, _I;
+    std::vector<double> _VA;
+    std::vector<double> _VI;  
+        
+    void fill_NLocal(const int pos, double N[][3]);
+    void fill_W_WI(const int pos, const int level, double W[][2]);
+   // void fill_W_WI(const int pos,  double W[][2]);
+    void fill_Matrices(const double VI);
+    void read_parameter(const char *ParaFile);
+     
+    #endif //SWIG
+  public:
+    StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile,const bool porous,const double tol=1e-6);
+    void addLaw(dG3DMaterialLaw* law);
+    #ifndef SWIG
+    StochDMNDG3DMaterialLaw(const StochDMNDG3DMaterialLaw &source);
+    virtual ~StochDMNDG3DMaterialLaw();
+
+    virtual void setTime(const double t,const double dtime);
+    virtual materialLaw::matname getType() const {return materialLaw::StochDMN;} 
+        
+    virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
+    virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
+    virtual void initialIPVariable(IPVariable* ipv, bool stiff);
+    
+    virtual void pos_kids(const int pos, const int level, int& pos0, int& pos1);
+    virtual void TwoPhasesInteract(const STensor43& C0, const STensor43& C1, fullMatrix<double> &mat, const int pos, const int level);
+    virtual void initLaws(const std::map<int,materialLaw*> &maplaw);
+    virtual double scaleFactor() const {return 1.;};
+    virtual double soundSpeed() const {return 0.;};
+    virtual materialLaw* clone() const {return new StochDMNDG3DMaterialLaw(*this);};
+    virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
+   
+    virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const;
+    virtual bool withEnergyDissipation() const {return false;};
+    virtual const materialLaw *getConstNonLinearSolverMaterialLaw() const
+    {
+      Msg::Error("getConstNonLinearSolverMaterialLaw in  StochDMNDG3DMaterialLaw does not exist");
+      return NULL;
+    }
+    virtual materialLaw *getNonLinearSolverMaterialLaw()
+    {
+      Msg::Error("getNonLinearSolverMaterialLaw in  StochDMNDG3DMaterialLaw does not exist");
+      return NULL;
+    }
+    virtual void setMacroSolver(const nonLinearMechSolver* sv);
+
+    #endif //SWIG
+ };   
+
+
+
 class J2SmallStrainDG3DMaterialLaw : public dG3DMaterialLaw
 {
  protected: