diff --git a/NonLinearSolver/materialLaw/elasticPotential.cpp b/NonLinearSolver/materialLaw/elasticPotential.cpp
index 926a87280aa1f18ecaa6333244e7c68941104edc..62fad959e354bb7eb7fd52a84d93808d34050efe 100644
--- a/NonLinearSolver/materialLaw/elasticPotential.cpp
+++ b/NonLinearSolver/materialLaw/elasticPotential.cpp
@@ -118,8 +118,7 @@ largeStrainAnisotropicPotential::largeStrainAnisotropicPotential(const int num,
 	const double Vzy= Vyz*Ez/Ey  ;
 	const double D=( 1-Vxy*Vyx-Vzy*Vyz-Vxz*Vzx-2*Vxy*Vyz*Vzx ) / ( Ex*Ey*Ez );
 
- 	static STensor43 ElasticityTensor;
-  STensorOperation::zero(ElasticityTensor);
+ 	STensor43 ElasticityTensor(0.);
         
  	ElasticityTensor(0,0,0,0)=( 1-Vyz*Vzy ) / (Ey*Ez*D );
  	ElasticityTensor(1,1,1,1)=( 1-Vxz*Vzx ) / (Ex*Ez*D );
@@ -160,7 +159,7 @@ largeStrainAnisotropicPotential::largeStrainAnisotropicPotential(const int num,
 	s1c2 = s1*c2;
 	c1c2 = c1*c2;
   
-  static STensor3 R; //3x3 rotation matrix
+  STensor3 R; //3x3 rotation matrix
 	R(0,0) = c3*c1 - s1c2*s3;
 	R(0,1) = c3*s1 + c1c2*s3;
 	R(0,2) = s2*s3;
@@ -324,7 +323,7 @@ biLogarithmicElasticPotential::biLogarithmicElasticPotential(const int num, cons
   _K = _E/(3.*(1.-2.*_nu));
   _mu = _E/(2.*(1.+_nu));
   double lambda = _K- 2.*_mu/3.;
-  static const STensor3 I(1.);
+  const STensor3 I(1.);
   for (int i=0; i<3; i++){
     for (int j=0; j<3; j++){
       for (int k=0; k<3; k++){
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 892dad69813571ab6803b568fa022fcfbde725d9..ecd3c59c9e4c97f261e0149745a6844a9e91a089 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -130,7 +130,8 @@
 																																			_GModelIsRotated(false), _checkFailureOnset(false),
 																																			_damageVolumeCenter(0.,0.,0.),_homogenizedCrackSurface(0.),
 																																			_microBCOld(NULL),_microFailureBC(NULL), _failureBCIsSwitched(false),
-                                                                      _pathFollowingSwitchCriterion(0.), _failureBasedOnPreviousState(true)
+                                                                      _pathFollowingSwitchCriterion(0.), _failureBasedOnPreviousState(true),
+                                                                      _currentTime(0.),_timeStep(0.)
 {
   // check parallelization of dofManager
  #if defined(HAVE_MPI)
@@ -11056,6 +11057,11 @@ double nonLinearMechSolver::solveMicroSolverForwardEuler(){
 	};
 };
 
+void nonLinearMechSolver::setTime(const double ctime,const double dtime){
+  _timeStep = dtime; 
+  _currentTime = ctime;
+}
+
 double nonLinearMechSolver::microSolveStaticLinear(){
   double t= Cpu();
   initMicroSolver();
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.h b/NonLinearSolver/nlsolver/nonLinearMechSolver.h
index 4963da0fd24896c693c9cb0d99aecdc2d0894fd2..f2ad270b26a02434b48aa0e734e3a3ad3fd6361d 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.h
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.h
@@ -707,6 +707,10 @@ class nonLinearMechSolver
   **/
  protected:
  #ifndef SWIG
+  // time and time step
+  double _timeStep; // for law which works on increment. (Use only in so no getTimeStep function)
+  double _currentTime; // To save results vs time
+	
   //for micro flag
   SYSTEM_TYPE _systemType;
   CONTROL_TYPE _controlType;
@@ -858,6 +862,7 @@ class nonLinearMechSolver
  #endif // SWIG
  public:
  #ifndef SWIG
+  void setTime(const double ctime,const double dtime);
   // RVE analysis newton raphson
   double microSolveStaticLinear();
   double microSolveSNL();
diff --git a/dG3D/src/dG3DMultiscaleMaterialLaw.cpp b/dG3D/src/dG3DMultiscaleMaterialLaw.cpp
index 27e4c9d5c5d1b33332ce888dda776ca647bf9378..82a0b0f2b6fae25a2a4a45b91c161d8e4adc403c 100644
--- a/dG3D/src/dG3DMultiscaleMaterialLaw.cpp
+++ b/dG3D/src/dG3DMultiscaleMaterialLaw.cpp
@@ -94,14 +94,15 @@ double dG3DElasticMultiscaleMaterialLaw::soundSpeed() const
 
 dG3DMultiscaleMaterialLaw::dG3DMultiscaleMaterialLaw(const int lnum, const int tag)
                             : dG3DMultiscaleMaterialLawBase(lnum,tag),
-														_blockDamageAfterFailureOnset(false),_lostSolutionUniquenssTolerance(0.),
+														_blockDamageAfterFailureOnset(false),_lostSolutionUniquenssTolerance(-1.e10),
 														_failureCrFollowingInterfaceNormal(true){};
 
 void dG3DMultiscaleMaterialLaw::initLaws(const std::map<int,materialLaw*> &maplaw){
   if (!_initialized){
     this->fillElasticStiffness(_rho,elasticStiffness);
 
-		if (_blockDamageAfterFailureOnset and (_lostSolutionUniquenssTolerance<=0.)){
+		if (_blockDamageAfterFailureOnset and (_lostSolutionUniquenssTolerance<0.)){
+      // find a smallest tolerance
 			_lostSolutionUniquenssTolerance = 1e10;
 			for (std::map<int,materialLaw*>::const_iterator it = maplaw.begin(); it != maplaw.end(); it++){
 				const GeneralBulkFollwedCohesive3DLaw* cohLaw = dynamic_cast<const GeneralBulkFollwedCohesive3DLaw*>(it->second);
@@ -112,8 +113,11 @@ void dG3DMultiscaleMaterialLaw::initLaws(const std::map<int,materialLaw*> &mapla
 				}
 			}
 			printf("tolerance of lost of ellipticity at bulk ip %f \n",_lostSolutionUniquenssTolerance);
-
 		}
+    
+    if (_lostSolutionUniquenssTolerance < 0){
+      _lostSolutionUniquenssTolerance = 0.;
+    }
     _initialized = true;
   }
 };
@@ -143,6 +147,11 @@ void dG3DMultiscaleMaterialLaw::createIPState(const bool isSolve, IPStateBase* &
 	IPVariable* ipv2 = new dG3DMultiscaleIPVariable(oninter);
 	if(ips != NULL) delete ips;
 	ips = new IP3State(state,ipvi,ipv1,ipv2);
+  
+  // if a cohesive law is added, ipstate at interface is created based on combined fracture law
+  // otherwise, only blocked failure will be considered 
+  // _blockDamageAfterFailureOnset is active for bulk homogenization
+  // and interface in case of fullDG formulation
 
   if (isSolve){
     int el = ele->getNum();
@@ -201,14 +210,18 @@ void dG3DMultiscaleMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
   solver->tangentAveragingFlag(stiff);
   // solve micro problem
   double time = solver->microSolve();
+  
   // get homogenized properties form solved micro problem
 	mipv->getRefToFirstPiolaKirchhoffStress() = solver->getHomogenizationState(IPStateBase::current)->getHomogenizedStress();
+  
 	mipv->getRefToDefoEnergy() = solver->getHomogenizationState(IPStateBase::current)->getDeformationEnergy();
 	mipv->getRefToPlasticEnergy() = solver->getHomogenizationState(IPStateBase::current)->getPlasticEnergy();
 	mipv->getRefToIrreversibleEnergy() = solver->getHomogenizationState(IPStateBase::current)->getIrreversibleEnergy();
 	
 	if (solver->checkFailureOnset()){
+    // if check failure onset
 		if (mipv->isInterface()){
+      // at interface
 			if (_failureCrFollowingInterfaceNormal){
 				solver->checkFailureOnset_withNormalVector();
 			}
@@ -217,12 +230,13 @@ void dG3DMultiscaleMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
 			}
 		}
 		else{
+      // at bulk element
 			solver->checkFailureOnset_withoutNormalVector();
 		}
-    mipv->getRefToLostEllipticityCriterion() = solver->getHomogenizationState(IPStateBase::current)->getLostOfSolutionUniquenessCriterion();
 	}
-	
+  mipv->getRefToLostEllipticityCriterion() = solver->getHomogenizationState(IPStateBase::current)->getLostOfSolutionUniquenessCriterion();
 	mipv->getRefToSolverBrokenFlag() =  solver->getHomogenizationState(IPStateBase::current)->getBrokenFlag();
+  
 	if (!(this->getMacroSolver()->withFailureBasedOnPreviousState())){ 
 		// when checking failure on the fly
 		// if current state is broken
@@ -274,7 +288,7 @@ void MultiscaleFractureByCohesive3DLaw::initLaws(const std::map<int,materialLaw*
 
     _nummat = dynamic_cast<numericalMaterialBase*>(bulkLaw);
 
-    if (_nummat == NULL) Msg::Error("bulk law is not a multiscale one");
+    if (_nummat == NULL) Msg::Error("bulk law is not a multiscale one in MultiscaleFractureByCohesive3DLaw::initLaws");
 
   }
 };
@@ -326,14 +340,18 @@ void MultiscaleFractureByCohesive3DLaw::createIPState(const bool isSolve, IPStat
 			 // if bulk law is a multiscale law
 			int el = ele->getNum();
 			nonLinearMechSolver* solver   = this->createMicroSolver(el,gpt);
+      
 			const GeneralMultiscaleBulkFollwedCohesive3DLaw* interfaceLaw = dynamic_cast<const GeneralMultiscaleBulkFollwedCohesive3DLaw*>(_mfrac);
 			if (interfaceLaw != NULL){
+        solver->setCheckFailureOnset(true); // set check failure onset
 				solver->setLostSolutionUniquenssTolerance(interfaceLaw->getLostSolutionUniquenssTolerance());
+        
+        // stress location
 				if (interfaceLaw->getLocationForStressExtraction()>=0){
 					solver->setLocationForStressExtraction(interfaceLaw->getLocationForStressExtraction());
 				}
-				solver->setCheckFailureOnset(true); // set check failure onset
-					
+				
+        // if extracting cohesive law
 				if (interfaceLaw->isDamageToCohesiveJump()){
 					printf("rank %d : extract cohesive law at element %d ip %d \n",Msg::GetCommRank(),el,gpt);
 					solver->setExtractCohesiveLawFromMicroDamage(true);
@@ -361,6 +379,7 @@ void MultiscaleFractureByCohesive3DLaw::createIPState(const bool isSolve, IPStat
 					solver->rotateModel(n,t,b);
 				}
 				
+        // update microbc with local basis
 				nonLinearMicroBC* mbc  = solver->getMicroBC();
 				if (mbc->getType() == nonLinearMicroBC::OrthogonalDirectionalMixedBC){
 					printf("set pbc %d n = %f %f %f t= %f %f %f",Msg::GetCommRank(),n[0],n[1],n[2],t[0],t[1],t[2]);
@@ -412,6 +431,9 @@ void MultiscaleFractureByCohesive3DLaw::createIPState(const bool isSolve, IPStat
 					
 				}
 			}
+      else{
+        Msg::Fatal("GeneralMultiscaleBulkFollwedCohesive3DLaw must be used in MultiscaleFractureByCohesive3DLaw::createIPState");
+      }
 			solver->initMicroSolver();
 			ips->setSolver(solver);
     }
@@ -608,12 +630,10 @@ void dG3DMultiscaleCohesiveLaw::transferInterfaceDataToBulk(IPVariable* ipv, con
     static SVector3 refN;
     refN = fipv->getConstRefToReferenceOutwardNormal();
     refN.normalize();
-    static SVector3 curN;
-    curN = fipvprev->getConstRefToCurrentOutwardNormal();
-    curN.normalize();
-    
-    double beta = (homoData->getActiveDissipationVolume() + _voidPartInLocalizationBand)/(_surfaceReductionRatio*solver->getRVEVolume());
+
+    double beta = 0.;
     if (!_withNormalPart){
+      beta = (homoData->getActiveDissipationVolume() + _voidPartInLocalizationBand)/(_surfaceReductionRatio*solver->getRVEVolume());
       static SVector3 normdF;
       for (int i=0; i<3; i++){
         normdF(i) = 0.;
@@ -697,6 +717,7 @@ void dG3DMultiscaleCohesiveLaw::stress(IPVariable* ipv, const IPVariable* ipvpre
     static SVector3 refN;
     refN = fipv->getConstRefToReferenceOutwardNormal();
     refN.normalize();
+    
     static SVector3 curN;
     curN = fipvprev->getConstRefToCurrentOutwardNormal();
     curN.normalize();
@@ -937,7 +958,6 @@ void gradientEnhanceddG3DMultiscaleCohesiveLaw::stress(IPVariable* ipv, const IP
     curN.normalize();
 
 
-
     const STensor3& P = fipv->getConstRefToFirstPiolaKirchhoffStress();
     SVector3& T = fipv->getRefToInterfaceForce();
     T *= 0.;
diff --git a/dG3D/src/dG3DMultiscaleMaterialLaw.h b/dG3D/src/dG3DMultiscaleMaterialLaw.h
index 8a91dc7a6a8c0cbe84ebf8e47724a4c43155779f..3f46732b06aa43146950f3900da9d45321cf2aa1 100644
--- a/dG3D/src/dG3DMultiscaleMaterialLaw.h
+++ b/dG3D/src/dG3DMultiscaleMaterialLaw.h
@@ -73,9 +73,9 @@ class dG3DElasticMultiscaleMaterialLaw : public dG3DMultiscaleMaterialLawBase{
 
 class dG3DMultiscaleMaterialLaw : public dG3DMultiscaleMaterialLawBase{
   protected:
-		bool _blockDamageAfterFailureOnset;
-		double _lostSolutionUniquenssTolerance;
-		bool _failureCrFollowingInterfaceNormal;
+		bool _blockDamageAfterFailureOnset; // true if failure is blocked after failure onset
+		double _lostSolutionUniquenssTolerance; 
+		bool _failureCrFollowingInterfaceNormal; // true if failure is checked based on interface normal
 
   public:
     dG3DMultiscaleMaterialLaw(const int lnum, const int tag);
@@ -222,7 +222,7 @@ class dG3DMultiscaleCohesiveLaw : public GeneralMultiscaleBulkFollwedCohesive3DL
   protected:
     bool _withNormalPart;
     bool _fixedCrackFaceFlag; // introduce a fixed crack length over RVE
-		bool _fixedCrackFace;
+		double _fixedCrackFace;
   public:
     dG3DMultiscaleCohesiveLaw(const int num, const double surfaceRatio=1., const bool normPart = true);
     void setFixedCrackFace(const double crackFace) {_fixedCrackFaceFlag = true; _fixedCrackFace = crackFace;};