diff --git a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
index d3dd364feba69d8d678eb0fdf0c9157073a5a6a9..3dfe11b0f7f45be084c5a0ac3c7322b4296b8d7c 100644
--- a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
+++ b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
@@ -14,14 +14,14 @@
 nonLinearMicroBC::nonLinearMicroBC(const int tag, const int dim) : _tag(tag),_dim(dim),Fcur(1.0),
                     Fprev(1.),Gcur(0.0),Gprev(0.),
                     _time(1.0),_order(1),_timeDependency(true),FI(0),G(0.){
-		
+
 	_totalMechaDofNumber = 3; // by default , functionSpace generates 3 DOFs per node
 	_totalConDofNumber = 0;
 	_totalNonConDofNumber = 0;
 	// if 2D problem is considered, please formulated problem in X-Y plan
 	// two fist DOF components (following x,y) are constrained
 	// third DOF (following z) is constrained by usual DIRICHLET BC
-	for (int i=0; i< _dim; i++){ 
+	for (int i=0; i< _dim; i++){
 		_constrainedDofs.insert(i);
 		_constrainedMechaDofs.insert(i);
 	}
@@ -82,7 +82,7 @@ void nonLinearMicroBC::setBCPhysical(const int g1, const int g2, const int g3, c
 	if (g4>0) _physical.push_back(g4);
 	if (g5>0) _physical.push_back(g5);
 	if (g6>0) _physical.push_back(g6);
-	// all 
+	// all
 };
 
 
@@ -95,16 +95,16 @@ void nonLinearMicroBC::setDofPartition(const int mechaDofPerNode, const int conD
 	// an arbitrary DOF can be constrained by using two functions
 	// clearAllConstraintDofs()
 	// insertConstrainDof(comp)
-	
+
 	clearAllConstraintDofs();
-	
+
 	// only mechanical DOF from
 	_totalMechaDofNumber = mechaDofPerNode;
 	for (int i=0; i< std::min(_dim,_totalMechaDofNumber); i++){
 		_constrainedDofs.insert(i);
 		_constrainedMechaDofs.insert(i);
 	}
-	
+
 	// constitive extarDof
 	_totalConDofNumber = conDofPerNode;
 	for (int i=0; i< _totalConDofNumber; i++){
@@ -128,7 +128,7 @@ void nonLinearMicroBC::setDofPartition(const int mechaDofPerNode, const int conD
 		Gprev *= 0.;
 		G*= 0.;
 	}
-	
+
 	gradTcur.clear();
 	gradTprev.clear();
 	gradT.clear();
@@ -145,7 +145,7 @@ void nonLinearMicroBC::setDofPartition(const int mechaDofPerNode, const int conD
 		T.push_back(0.);
 		Tinit.push_back(0.);
 	}
-	
+
 	gradVcur.clear();
 	gradVprev.clear();
 	gradV.clear();
@@ -243,7 +243,7 @@ void nonLinearMicroBC::clearAllConstraintDofs(){
 };
 void nonLinearMicroBC::insertConstrainDof(const int dofNum){
 	_constrainedDofs.insert(dofNum);
-	// check what kind of DOF 
+	// check what kind of DOF
 	if (dofNum < _totalMechaDofNumber) {
 		_constrainedMechaDofs.insert(dofNum);
 	}
@@ -251,9 +251,9 @@ void nonLinearMicroBC::insertConstrainDof(const int dofNum){
 		_constrainedConDofs.insert(dofNum);
 	}
 	else if ((dofNum>= _totalMechaDofNumber+_totalConDofNumber) and (dofNum < _totalMechaDofNumber+_totalConDofNumber+_totalNonConDofNumber)){
-	
+
 	}
-	else 
+	else
 		Msg::Fatal("component %d cannot be constrained ",dofNum);
 };
 
@@ -264,7 +264,7 @@ nonLinearPeriodicBC::nonLinearPeriodicBC(const int tag, const int dim):
                         _ZaddNewVertices(false),_Zdegree(3),
                         _wM(nonLinearPeriodicBC::CEM),_forceInterp(false),
 												_isShifted(false),_shiftPBCNormal(0.,0.,0.){ }
-nonLinearPeriodicBC::nonLinearPeriodicBC(const nonLinearPeriodicBC& src):nonLinearMicroBC(src), 
+nonLinearPeriodicBC::nonLinearPeriodicBC(const nonLinearPeriodicBC& src):nonLinearMicroBC(src),
   _wM(src._wM),_Xdegree(src._Xdegree),_XaddNewVertices(src._XaddNewVertices),_Ydegree(src._Ydegree),
   _YaddNewVertices(src._YaddNewVertices),_Zdegree(src._Zdegree),_ZaddNewVertices(src._ZaddNewVertices),
   _forceInterp(src._forceInterp),_isShifted(src._isShifted),_shiftPBCNormal(src._shiftPBCNormal)
@@ -341,6 +341,7 @@ void nonLinearPeriodicBC::setShiftPBCNormal(const double n1, const double n2, co
 	_shiftPBCNormal[0] = n1;
 	_shiftPBCNormal[1] = n2;
 	_shiftPBCNormal[2] = n3;
+	printf("set Shifted PBC = [%f %f %f]\n",n1,n2,n3);
 	_shiftPBCNormal.normalize();
 };
 void nonLinearPeriodicBC::setShiftPBCNormal(const SVector3& n){
@@ -353,7 +354,7 @@ void nonLinearPeriodicBC::updateWithSolver(const nonLinearMechSolver* s){
 	printf("shitfed with lost of solution uniqueness normal : n = [%f %f %f]",_shiftPBCNormal[0],_shiftPBCNormal[1],_shiftPBCNormal[2]);
 };
 
-												
+
 nonLinearDisplacementBC::nonLinearDisplacementBC( const int tag, const int dim):
                     nonLinearMicroBC(tag,dim){};
 
@@ -365,12 +366,12 @@ nonLinearMinimalKinematicBC::nonLinearMinimalKinematicBC(const int tag, const in
 nonLinearMinimalKinematicBC::nonLinearMinimalKinematicBC(const nonLinearMinimalKinematicBC& src): nonLinearMicroBC(src){}
 
 nonLinearMixedBC::nonLinearMixedBC(const int tag, const int dim) : nonLinearMicroBC(tag,dim),_negativeFactor(1.){}
-nonLinearMixedBC::nonLinearMixedBC(const nonLinearMixedBC& src) : nonLinearMicroBC(src), 
+nonLinearMixedBC::nonLinearMixedBC(const nonLinearMixedBC& src) : nonLinearMicroBC(src),
   _kPhysical(src._kPhysical),_sPhysical(src._sPhysical),_fixVertices(src._fixVertices),
  _additionalSecondOrderStaticPhysical(src._additionalSecondOrderStaticPhysical),_negativeFactor(src._negativeFactor),
- _periodicNegative(src._periodicNegative),_periodicPositive(src._periodicPositive){} 
- 
-nonLinearOrthogonalMixedBCDirectionFollowing::nonLinearOrthogonalMixedBCDirectionFollowing(const int tag, const int dim): 
+ _periodicNegative(src._periodicNegative),_periodicPositive(src._periodicPositive){}
+
+nonLinearOrthogonalMixedBCDirectionFollowing::nonLinearOrthogonalMixedBCDirectionFollowing(const int tag, const int dim):
 nonLinearMicroBC(tag,dim){}
 
 nonLinearOrthogonalMixedBCDirectionFollowing::nonLinearOrthogonalMixedBCDirectionFollowing(const nonLinearOrthogonalMixedBCDirectionFollowing& src):
@@ -397,16 +398,16 @@ void nonLinearOrthogonalMixedBCDirectionFollowing::setSUBCDirection(const double
 void nonLinearOrthogonalMixedBCDirectionFollowing::clearDirs(){ _KUBCDirection.clear(); _SUBCDirection.clear();};
 
 
-nonLinearMixedBCDirectionFollowing::nonLinearMixedBCDirectionFollowing(const int tag, const int dim): 
+nonLinearMixedBCDirectionFollowing::nonLinearMixedBCDirectionFollowing(const int tag, const int dim):
 	nonLinearMicroBC(tag,dim){}
-	
+
 nonLinearMixedBCDirectionFollowing::nonLinearMixedBCDirectionFollowing(const nonLinearMixedBCDirectionFollowing& src):
 	nonLinearMicroBC(src),_KUBCDirection(src._KUBCDirection),_SUBCDirection(src._SUBCDirection),
 	_PBCNormalDirection(src._PBCNormalDirection), _PBCDirection(src._PBCDirection),
 	_rootPhysical(src._rootPhysical),
 	_interpolationDegree(src._interpolationDegree),
 	_fullConstrained(src._fullConstrained){}
-	
+
 void nonLinearMixedBCDirectionFollowing::setKUBCDirection(const SVector3& dir){
 	_KUBCDirection.push_back(dir);
 };
@@ -419,7 +420,7 @@ void nonLinearMixedBCDirectionFollowing::setPBCDirection(const SVector3& dir, co
 	_PBCNormalDirection.push_back(norm);
 	_fullConstrained.push_back(full);
 };
-	
+
 void nonLinearMixedBCDirectionFollowing::setKUBCDirection(const double nx, const double ny, const double nz){
 	SVector3 n(nx,ny,nz);
 	n.normalize();
@@ -436,11 +437,11 @@ void nonLinearMixedBCDirectionFollowing::setPBCDirection(const double nx, const
 	SVector3 n(nx,ny,nz);
 	n.normalize();
 	_PBCNormalDirection.push_back(n);
-	
+
 	SVector3 t(tx,ty,tz);
 	t.normalize();
 	_PBCDirection.push_back(t);
-	
+
 	_fullConstrained.push_back(fullConstrained);
 };
 
@@ -450,4 +451,4 @@ void nonLinearMixedBCDirectionFollowing::setRootPhysical(const int phy){
 
 void nonLinearMixedBCDirectionFollowing::setInterpolationDegree(const int deg){
 	_interpolationDegree.push_back(deg);
-};
\ No newline at end of file
+};
diff --git a/NonLinearSolver/materialLaw/elasticPotential.cpp b/NonLinearSolver/materialLaw/elasticPotential.cpp
index 3d306fc38b37159420e7f29913e215cc41e08630..573ef440acbc50c5673cd2fc15ba6804f8436e9d 100644
--- a/NonLinearSolver/materialLaw/elasticPotential.cpp
+++ b/NonLinearSolver/materialLaw/elasticPotential.cpp
@@ -9,40 +9,26 @@
 #include "elasticPotential.h"
 
 smallStrainLinearElasticPotential::smallStrainLinearElasticPotential(const int num, const double E, const double nu): elasticPotential(num), _E(E), _nu(nu){
-	
+
 	_K = _E/(3.*(1.-2.*_nu));
 	_mu = _E/(2.*(1.+_nu));
-	
-	_Cel(0,0,0,0) = (1.-_nu)*_E/(1.+_nu)/(1.-2.*_nu);
-	_Cel(1,1,1,1) = _Cel(0,0,0,0);
-	_Cel(2,2,2,2) = _Cel(0,0,0,0);
-
-	_Cel(0,0,1,1) = nu*_E/(1.+_nu)/(1.-2.*_nu);
-	_Cel(0,0,2,2) = _Cel(0,0,1,1);
-	_Cel(1,1,0,0) = _Cel(0,0,1,1);
-	_Cel(1,1,2,2) = _Cel(0,0,1,1);
-	_Cel(2,2,0,0) = _Cel(0,0,1,1);
-	_Cel(2,2,1,1) = _Cel(0,0,1,1);
-
-	_Cel(0,1,0,1) = _mu;
-	_Cel(0,1,1,0) = _mu;
-	_Cel(1,0,0,1) = _mu;
-	_Cel(1,0,1,0) = _mu;
-	
-	_Cel(0,2,0,2) = _mu;
-	_Cel(0,2,2,0) = _mu;
-	_Cel(2,0,0,2) = _mu;
-	_Cel(2,0,2,0) = _mu;
-
-	_Cel(1,2,1,2) = _mu;
-	_Cel(1,2,2,1) = _mu;
-	_Cel(2,1,1,2) = _mu;
-	_Cel(2,1,2,1) = _mu;
+	double lambda = _K- 2.*_mu/3.;
+	STensor3 I(1.);
+	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++){
+          _Cel(i,j,k,l) = lambda*I(i,j)*I(k,l)+_mu*(I(i,k)*I(j,l)+I(i,l)*I(j,k));
+        }
+      }
+    }
+	}
 };
 smallStrainLinearElasticPotential::smallStrainLinearElasticPotential(const smallStrainLinearElasticPotential& src) : elasticPotential(src),
 	_E(src._E),_nu(src._nu),_K(src._K),_mu(src._mu), _Cel(src._Cel){};
 
 double smallStrainLinearElasticPotential::get(const STensor3& F) const{
+  // small strain value
 	STensor3 E(0.);
 	for (int i=0; i<3; i++){
 		E(i,i) -= 1.;
@@ -85,20 +71,19 @@ void smallStrainLinearElasticPotential::constitutive(const STensor3& F, STensor3
 	}
 };
 void smallStrainLinearElasticPotential::getLocalValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const{
-	double Y = get(F);
-	double E = getYoungModulus();
-	
-	val = sqrt(2*Y/E);
+	double Y = get(F); // potential
+	double E = getYoungModulus(); // young modulus
+
+	val = sqrt(2*Y/E); // isotropic nonlocal strain
 	if (stiff){
 		this->constitutive(F,*DvalDF,false,NULL);
-		
 		if (val > 1e-10){
 			(*DvalDF)*= (1./(E*val));
 		}
 		else{
 			(*DvalDF)*= (1./(E*1e-10));
 		}
-		 
+
 	}
 };
 
@@ -106,46 +91,32 @@ biLogarithmicElasticPotential::biLogarithmicElasticPotential(const int num, cons
 	elasticPotential(num), _E(E), _nu(nu), _order(order), _ByPerturbation(pert),_tol(tol){
 	_K = _E/(3.*(1.-2.*_nu));
 	_mu = _E/(2.*(1.+_nu));
-	_Cel(0,0,0,0) = (1.-_nu)*_E/(1.+_nu)/(1.-2.*_nu);
-	_Cel(1,1,1,1) = _Cel(0,0,0,0);
-	_Cel(2,2,2,2) = _Cel(0,0,0,0);
-
-	_Cel(0,0,1,1) = nu*_E/(1.+_nu)/(1.-2.*_nu);
-	_Cel(0,0,2,2) = _Cel(0,0,1,1);
-	_Cel(1,1,0,0) = _Cel(0,0,1,1);
-	_Cel(1,1,2,2) = _Cel(0,0,1,1);
-	_Cel(2,2,0,0) = _Cel(0,0,1,1);
-	_Cel(2,2,1,1) = _Cel(0,0,1,1);
-
-	_Cel(0,1,0,1) = _mu;
-	_Cel(0,1,1,0) = _mu;
-	_Cel(1,0,0,1) = _mu;
-	_Cel(1,0,1,0) = _mu;
-	
-	_Cel(0,2,0,2) = _mu;
-	_Cel(0,2,2,0) = _mu;
-	_Cel(2,0,0,2) = _mu;
-	_Cel(2,0,2,0) = _mu;
-
-	_Cel(1,2,1,2) = _mu;
-	_Cel(1,2,2,1) = _mu;
-	_Cel(2,1,1,2) = _mu;
-	_Cel(2,1,2,1) = _mu;
+	double lambda = _K- 2.*_mu/3.;
+	STensor3 I(1.);
+	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++){
+          _Cel(i,j,k,l) = lambda*I(i,j)*I(k,l)+_mu*(I(i,k)*I(j,l)+I(i,l)*I(j,k));
+        }
+      }
+    }
+	}
 };
 
 biLogarithmicElasticPotential::biLogarithmicElasticPotential(const biLogarithmicElasticPotential& src) : elasticPotential(src),
 	_E(src._E),_nu(src._nu),_K(src._K),_mu(src._mu), _Cel(src._Cel),_order(src._order),
 	_ByPerturbation(src._ByPerturbation),_tol(src._tol){};
-	
+
 double biLogarithmicElasticPotential::get(const STensor3& F) const{
 	STensor3 C, E;
   for (int i=0; i<3; i++)
     for (int j=0; j<3; j++)
       C(i,j) = F(0,i)*F(0,j)+F(1,i)*F(1,j)+F(2,i)*F(2,j);
-			
+
 	STensorOperation::logSTensor3(C,_order,E);
 	E *= 0.5;
-	
+
 	double val = 0.;
 	for (int i=0; i<3;i++){
 		for (int j=0; j<3; j++){
@@ -170,7 +141,7 @@ void biLogarithmicElasticPotential::constitutive(const STensor3& F, STensor3& P,
 				P(i,j) = (potentialPlus - potential)/_tol;
 			}
 		}
-		
+
 		if (stiff){
 			(*L) *= 0.;
 			STensor3 Pp(0.), Fp(0.);
@@ -189,17 +160,22 @@ void biLogarithmicElasticPotential::constitutive(const STensor3& F, STensor3& P,
 		}
 	}
 	else{
-		
+
 		STensor3 C, E;
 		for (int i=0; i<3; i++)
 			for (int j=0; j<3; j++)
 				C(i,j) = F(0,i)*F(0,j)+F(1,i)*F(1,j)+F(2,i)*F(2,j);
-		
-		STensor43 dlnsqrtC(0.);
-		STensorOperation::logSTensor3(C,_order,E,&dlnsqrtC);
+
+		STensor43 dlnCdC(0.);
+		STensor63 ddlnCdC(0.);
+		if (stiff){
+      STensorOperation::logSTensor3(C,_order,E,&dlnCdC,&ddlnCdC);
+		}
+		else{
+      STensorOperation::logSTensor3(C,_order,E,&dlnCdC);
+		}
 		E *= 0.5;
-		dlnsqrtC *= 0.5;
-		
+
 		STensor3 sig(0.);
 		for (int i=0; i<3;i++){
 			for (int j=0; j<3; j++){
@@ -210,61 +186,196 @@ void biLogarithmicElasticPotential::constitutive(const STensor3& F, STensor3& P,
 				}
 			}
 		}
-		
-		for (int i=0; i<3; i++)
-			for (int j=0; j<3; j++)
-				P(i,j) = F(i,0)*sig(0,j)+F(i,1)*sig(1,j)+F(i,2)*sig(2,j);
-				
+
+		STensor3 PK2(0.);
+		STensorOperation::multSTensor43STensor3(dlnCdC,sig,PK2);
+
+		for (int i=0; i<3; i++){
+			for (int j=0; j<3; j++){
+				P(i,j) = F(i,0)*PK2(0,j)+F(i,1)*PK2(1,j)+F(i,2)*PK2(2,j);
+      }
+    }
+
+		if (stiff){
+      STensor43 dPK2dC(0.);
+      for (int k=0; k<3; k++){
+        for (int j=0; j<3; j++){
+          for (int m=0; m<3; m++){
+            for (int n=0; n<3; n++){
+              for (int q=0; q<3; q++){
+                for (int s=0; s<3; s++){
+                  dPK2dC(k,j,q,s) += sig(m,n)*ddlnCdC(m,n,k,j,q,s);
+                  for (int u=0; u<3; u++){
+                    for (int v=0; v<3; v++){
+                      dPK2dC(k,j,q,s) += 0.5*_Cel(m,n,u,v)*dlnCdC(m,n,k,j)*dlnCdC(u,v,q,s);
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+
+      (*L) *= 0.;
+      for (int i=0; i<3; i++){
+        for (int j=0; j<3; j++){
+          for (int p=0; p<3; p++){
+            (*L)(i,j,i,p) += PK2(p,j);
+            for (int q=0; q<3; q++){
+              for (int s=0; s<3; s++){
+                for (int k=0; k<3; k++){
+                  (*L)(i,j,p,q) += 2.*F(i,k)*dPK2dC(k,j,q,s)*F(p,s);
+                }
+              }
+            }
+          }
+        }
+      }
+		}
+
+	}
+
+};
+void biLogarithmicElasticPotential::getLocalValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const{
+	double Y = get(F);
+	double E = getYoungModulus();
+
+	val = sqrt(2*Y/E);
+	if (stiff){
+		(*DvalDF) *= 0.;
+		if (_ByPerturbation){
+			for (int i=0; i<3; i++){
+				for (int j=0; j<3; j++){
+					STensor3 Fp(F);
+					Fp(i,j) += _tol;
+					double valp;
+					getLocalValForNonLocalEquation(Fp,valp,false,NULL);
+					(*DvalDF)(i,j) = (valp-val)/_tol;
+				}
+			}
+		}
+		else{
+			this->constitutive(F,*DvalDF,false,NULL);
+			if (val > 1e-10){
+				(*DvalDF) *= (1./(E*val));
+			}
+			else{
+				(*DvalDF) *= (1./(E*1e-10));
+			}
+
+		}
+	}
+};
+
+NeoHookeanElasticPotential::NeoHookeanElasticPotential(const int num, const double E, const double nu, const bool pert, const double tol) :
+	elasticPotential(num), _E(E), _nu(nu), _ByPerturbation(pert),_tol(tol){
+	_K = _E/(3.*(1.-2.*_nu));
+	_mu = _E/(2.*(1.+_nu));
+	double lambda = _K- 2.*_mu/3.;
+	STensor3 I(1.);
+	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++){
+          _Cel(i,j,k,l) = lambda*I(i,j)*I(k,l)+_mu*(I(i,k)*I(j,l)+I(i,l)*I(j,k));
+        }
+      }
+    }
+	}
+};
+
+NeoHookeanElasticPotential::NeoHookeanElasticPotential(const NeoHookeanElasticPotential& src) : elasticPotential(src),
+	_E(src._E),_nu(src._nu),_K(src._K),_mu(src._mu), _Cel(src._Cel),
+	_ByPerturbation(src._ByPerturbation),_tol(src._tol){};
+
+double NeoHookeanElasticPotential::get(const STensor3& F) const{
+	double J = F.determinant();
+  double J2over3= pow(J,2./3.);
+  STensor3 C(0.);
+  for (int i=0; i<3; i++)
+    for (int j=0; j<3; j++)
+      C(i,j) = (F(0,i)*F(0,j)+F(1,i)*F(1,j)+F(2,i)*F(2,j))/J2over3;
+  
+  double val = 0.5*_K*log(J)*log(J) + 0.5*_mu*(C.trace()-3.);
+  return val;
+};
+void NeoHookeanElasticPotential::constitutive(const STensor3& F, STensor3& P, const bool stiff, STensor43* L) const{
+	if (_ByPerturbation){
+		double potential = get(F);
+		P *= 0.;
+		for (int i=0; i<3; i++){
+			for (int j=0; j<3; j++){
+				STensor3 Fp(F);
+				Fp(i,j) += _tol;
+				double potentialPlus = get(Fp);
+				P(i,j) = (potentialPlus - potential)/_tol;
+			}
+		}
+
 		if (stiff){
-			static STensor3 I(1.);
-			static STensor43 dCdF(0.);
-			dCdF *= 0.;
+			(*L) *= 0.;
+			STensor3 Pp(0.), Fp(0.);
 			for (int i=0; i<3; i++){
 				for (int j=0; j<3; j++){
-					for (int m=0; m<3; m++){
-						for (int n=0; n<3; n++){
-							dCdF(i,j,m,n) += I(i,n)*F(m,j)+F(m,i)*I(j,n);
+					Fp = F;
+					Fp(i,j) += _tol;
+					constitutive(Fp,Pp,false,NULL);
+					for (int p=0; p<3; p++){
+						for (int q=0; q<3; q++){
+							(*L)(p,q,i,j) = (Pp(p,q) - P(p,q))/_tol;
 						}
 					}
 				}
 			}
-		
-		
-			static STensor43 dlnsqrtCDCDF(0.);
-			dlnsqrtCDCDF *= 0.;
-			STensorOperation::multSTensor43(dlnsqrtC,dCdF,dlnsqrtCDCDF);
-			
-			static STensor43 DsigDF(0.);
-			DsigDF *= 0.;
-			STensorOperation::multSTensor43(_Cel,dlnsqrtCDCDF,DsigDF);
-			
+		}
+	}
+	else{
+    double J = F.determinant();
+    double J2over3= pow(J,2./3.);
+    double logJ = log(J);
+    STensor3 Finv = F.invert();
+    STensor3 C(0.);
+    for (int i=0; i<3; i++)
+      for (int j=0; j<3; j++)
+        C(i,j) = (F(0,i)*F(0,j)+F(1,i)*F(1,j)+F(2,i)*F(2,j));
+    
+    double trC = C.trace();
+    
+    for (int i=0; i<3; i++){
+      for (int j=0; j<3; j++){
+        P(i,j) = (_mu/J2over3)*(F(i,j) - Finv(j,i)*(trC/3.)) + _K*logJ*Finv(j,i);
+      }
+    }
+    
+    if (stiff){
 			(*L) *= 0.;
+			STensor3 Pp(0.), Fp(0.);
 			for (int i=0; i<3; i++){
 				for (int j=0; j<3; j++){
-					for (int m=0; m<3; m++){
-						for (int n=0; n<3; n++){
-							(*L)(i,j,m,n) += I(i,m)*sig(n,j);
-							for (int k=0; k<3; k++){
-								(*L)(i,j,m,n) += F(i,k)*DsigDF(k,j,m,n);
-							}
+					Fp = F;
+					Fp(i,j) += _tol;
+					constitutive(Fp,Pp,false,NULL);
+					for (int p=0; p<3; p++){
+						for (int q=0; q<3; q++){
+							(*L)(p,q,i,j) = (Pp(p,q) - P(p,q))/_tol;
 						}
 					}
 				}
 			}
 		}
-		
-	}
 	
+	}
+
 };
-void biLogarithmicElasticPotential::getLocalValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const{
+void NeoHookeanElasticPotential::getLocalValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const{
 	double Y = get(F);
 	double E = getYoungModulus();
 
 	val = sqrt(2*Y/E);
 	if (stiff){
 		(*DvalDF) *= 0.;
-		//if (_ByPerturbation){
-    if (1){
+		if (_ByPerturbation){
 			for (int i=0; i<3; i++){
 				for (int j=0; j<3; j++){
 					STensor3 Fp(F);
@@ -276,14 +387,15 @@ void biLogarithmicElasticPotential::getLocalValForNonLocalEquation(const STensor
 			}
 		}
 		else{
-			this->constitutive(F,*DvalDF,false,NULL);			
+			this->constitutive(F,*DvalDF,false,NULL);
 			if (val > 1e-10){
 				(*DvalDF) *= (1./(E*val));
 			}
 			else{
 				(*DvalDF) *= (1./(E*1e-10));
 			}
-			
+
 		}
 	}
 };
+
diff --git a/NonLinearSolver/materialLaw/elasticPotential.h b/NonLinearSolver/materialLaw/elasticPotential.h
index 8044cbbfbc095d0ba4b46f445303989d6d66662d..b9a4619ee425af8afa2f2d8223cae784e6bd496d 100644
--- a/NonLinearSolver/materialLaw/elasticPotential.h
+++ b/NonLinearSolver/materialLaw/elasticPotential.h
@@ -79,4 +79,28 @@ class biLogarithmicElasticPotential : public elasticPotential{
 		#endif //SWIG
 };
 
+class NeoHookeanElasticPotential : public elasticPotential{
+  protected:
+    #ifndef SWIG
+    double _K, _mu, _E, _nu;
+    STensor43 _Cel;
+		bool _ByPerturbation;
+		double _tol;
+    #endif //SWIG
+  public:
+    NeoHookeanElasticPotential(const int num, const double E, const double nu, const bool pert, const double tol =1e-8);
+		#ifndef SWIG
+		NeoHookeanElasticPotential(const NeoHookeanElasticPotential& src);
+		virtual ~NeoHookeanElasticPotential(){}	
+		virtual double getYoungModulus() const {return _E;};
+		virtual double getPoissonRatio() const {return _nu;};
+		virtual double get(const STensor3& F) const; // elastic potential
+		virtual void constitutive(const STensor3& F, STensor3& P, const bool sitff, STensor43* L) const; // constitutive law
+		virtual void getLocalValForNonLocalEquation(const STensor3& F, double& val, const bool stiff, STensor3* DvalDF) const;
+		virtual elasticPotential* clone() const {
+			return new NeoHookeanElasticPotential(*this);
+		}
+    #endif //SWIG
+};
+
 #endif //ELASTICPOTENTIAL_H_
\ No newline at end of file
diff --git a/NonLinearSolver/materialLaw/mlawHyperelasticDamage.cpp b/NonLinearSolver/materialLaw/mlawHyperelasticDamage.cpp
index d982ec4cd16f3c28078e20bfb5f4effdf4842a8f..7caf938392c544a84329698cd4e5b14997eacfaa 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelasticDamage.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelasticDamage.cpp
@@ -12,7 +12,7 @@
 mlawNonlocalDamageHyperelastic::mlawNonlocalDamageHyperelastic(const int num, const double rho, const elasticPotential& elaw, const CLengthLaw &cLLaw,const DamageLaw &damLaw):
 materialLaw(num,true), _rho(rho)
 {
-	_elasticLaw = elaw.clone(); 
+	_elasticLaw = elaw.clone();
 	_cLLaw = cLLaw.clone();
 	_damLaw = damLaw.clone();
 };
@@ -20,13 +20,13 @@ materialLaw(num,true), _rho(rho)
 mlawNonlocalDamageHyperelastic::mlawNonlocalDamageHyperelastic(const mlawNonlocalDamageHyperelastic& src): materialLaw(src){
 	_elasticLaw = NULL;
 	if (src._elasticLaw!=NULL) _elasticLaw = src._elasticLaw->clone();
-	
+
 	_cLLaw = NULL;
 	if (src._cLLaw != NULL) _cLLaw = src._cLLaw->clone();
-	
+
 	_damLaw = NULL;
 	if (src._damLaw != NULL) _damLaw = src._damLaw->clone();
-	
+
 };
 
 mlawNonlocalDamageHyperelastic::~mlawNonlocalDamageHyperelastic(){
@@ -63,7 +63,7 @@ void mlawNonlocalDamageHyperelastic::constitutive(
 	// compute current nonlocal lengths matrix
 	double p0 = q0->getLocalEffectiveStrains();
   _cLLaw->computeCL(p0, q->getRefToIPCLength());
-	
+
 	// elastic law
 	static STensor3 effP;
 	effP *= 0.;
@@ -71,14 +71,15 @@ void mlawNonlocalDamageHyperelastic::constitutive(
 		Msg::Fatal("elastic potential is NULL");
 	}
 	_elasticLaw->constitutive(F,effP,stiff,&dPdF); // compute
-	
+
 	// elastic potential
 	q->getRefToElasticEnergy() = _elasticLaw->get(F);
-	
-	// local strain	
+
+	// local strain
 	_elasticLaw->getLocalValForNonLocalEquation(F,q->getRefToLocalEffectiveStrains(),stiff,&dLocalVardF);
-	
-	if (q0->damageIsBlocked()){
+
+
+	if (q0->damageIsBlocked() or (q0->getDamage()> 0.99999)){
     q->blockDamage(true);
     q->setNonLocalToLocal(false);
     q->setCriticalDamage(q0->getCriticalDamage());
@@ -88,14 +89,14 @@ void mlawNonlocalDamageHyperelastic::constitutive(
     curDama.getRefToDDamageDp() = 0.;
     curDama.getRefToDDamageDFe() *= 0.;
     curDama.getRefToDeltaDamage() = 0;
-    curDama.getRefToMaximalP() = q0->getMaximalP();		
+    curDama.getRefToMaximalP() = q0->getMaximalP();
 		q->activeDamage(false);
   }
   else
   {
     q->blockDamage(false);
 		static STensor3 Fp(1.);
-		
+
     if (q0->getNonLocalToLocal()){
       q->setNonLocalToLocal(true);
       q->getRefToNonLocalEffectiveStrains() = q0->getNonLocalEffectiveStrains() + (q->getLocalEffectiveStrains() - q0->getLocalEffectiveStrains());
@@ -135,10 +136,11 @@ void mlawNonlocalDamageHyperelastic::constitutive(
   }
 
   double D = q->getDamage();
-	P = effP;
+
+  P = effP;
   P *= (1.-D);
 
-	// elastic energy
+  // elastic energy
   q->getRefToElasticEnergy() *= (1.-D);
 
 
@@ -149,7 +151,7 @@ void mlawNonlocalDamageHyperelastic::constitutive(
       dPdNonlocalVar *= 0.;
     }
     else{
-			dPdNonlocalVar = effP;
+      dPdNonlocalVar = effP;
       dPdNonlocalVar *= (-q->getDDamageDp());
     }
 
@@ -167,7 +169,7 @@ void mlawNonlocalDamageHyperelastic::constitutive(
 
     dLocalVardNonlocalVar *= 0.;
   }
-  
+
   q->getRefToIrreversibleEnergy() = q0->irreversibleEnergy() + q->defoEnergy()/(1.-q->getDamage()) * (D-q0->getDamage());
   if (stiff){
     if (q0->getNonLocalToLocal()){
@@ -188,5 +190,5 @@ void mlawNonlocalDamageHyperelastic::constitutive(
       q->getRefToDIrreversibleEnergyDNonLocalVariable() = q->defoEnergy()/(1.-q->getDamage()) * q->getDDamageDp();
     }
   }
-	
-};
\ No newline at end of file
+
+};
diff --git a/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp b/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
index f75702705da4f089a7b9f5b1f536bbfae7593b1e..9b9bb338cb0f928de52000f35c7d39d596526ec2 100644
--- a/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
+++ b/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
@@ -452,7 +452,7 @@ void pbcConstraintElementGroup::__init__(){
 			Msg::Info("Corner vertex = %d",v->getNum());
 		}
 		#endif //_DEBUG
-		
+
 		addControlPoints();
 		addOtherCornerPoints();
 	}
@@ -694,29 +694,29 @@ void pbcConstraintElementGroup::addControlPoints(){
     }
 
   }
-	
+
 	_periodic12[0] = _v2->x() - _v1->x();
 	_periodic12[1] = _v2->y() - _v1->y();
 	_periodic12[2] = _v2->z() - _v1->z();
-	
+
 	_periodic14[0] = _v4->x() - _v1->x();
 	_periodic14[1] = _v4->y() - _v1->y();
 	_periodic14[2] = _v4->z() - _v1->z();
-	
+
 	if (_solver->getMicroBC()->getDim() == 3){
 		_periodic15[0] = _v5->x() - _v1->x();
 		_periodic15[1] = _v5->y() - _v1->y();
 		_periodic15[2] = _v5->z() - _v1->z();
 	}
-	
+
 	_periodic12.print("periodic 1 --> 2");
 	_periodic14.print("periodic 1 --> 4");
-	
-	
+
+
 	if (_solver->getMicroBC()->getDim() == 3){
 		_periodic15.print("periodic 1 --> 5");
 	}
-	
+
 };
 
 void pbcConstraintElementGroup::addOtherCornerPoints(){
@@ -987,7 +987,7 @@ void pbcConstraintElementGroup::intersectionBetweenTwoFaces(const groupOfElement
 
 		}
 	}
-	
+
 	// add to grEle
 	if (line == NULL) line = new groupOfElements();
 	for (int i=0; i< allEdges.size(); i++){
@@ -997,7 +997,7 @@ void pbcConstraintElementGroup::intersectionBetweenTwoFaces(const groupOfElement
 	}
 };
 
-void pbcConstraintElementGroup::intersectionBetweenTwoLines(const groupOfElements* line1, const groupOfElements* line2, groupOfElements* & pt) const{	
+void pbcConstraintElementGroup::intersectionBetweenTwoLines(const groupOfElements* line1, const groupOfElements* line2, groupOfElements* & pt) const{
 	std::vector<MVertex*> allVertex;
 	for (groupOfElements::vertexContainer::iterator itL = line1->vbegin(); itL!= line1->vend(); itL++){
 		MVertex* v1 = *itL;
@@ -1009,7 +1009,7 @@ void pbcConstraintElementGroup::intersectionBetweenTwoLines(const groupOfElement
 			}
 		}
 	}
-	
+
 	if (pt == NULL) pt = new groupOfElements();
 	for (int i=0; i< allVertex.size(); i++){
 		MElement* ele = new MPoint(allVertex[i]);
@@ -1170,7 +1170,7 @@ void pbcConstraintElementGroup::periodicConditionForPeriodicMesh(Iterator posivt
 };
 
 template<class Iterator>
-void pbcConstraintElementGroup::periodicConditionForPeriodicMesh(const int comp, const double fact, 
+void pbcConstraintElementGroup::periodicConditionForPeriodicMesh(const int comp, const double fact,
 																							 Iterator posivtiveitbegin, Iterator positiveitend,
                                                Iterator negativeitbegin, Iterator negativeitend,
                                                std::set<MVertex*>& others, MVertex* vrootP, MVertex* vrootN){
@@ -1460,7 +1460,7 @@ void pbcConstraintElementGroup::lagrangePeriodicCondition(Iterator itbegin, Iter
 template<class Iterator>
 void pbcConstraintElementGroup::shiftedLagrangePeriodicCondition(Iterator itbegin, Iterator itend,
                                    std::vector<MVertex*>& vlist,std::set<MVertex*>& others, const SVector3& normal){
-																		 
+
 	for (Iterator it = itbegin; it != itend; it++){
     MVertex* v= *it;
     if (others.find(v) == others.end()){
@@ -1486,13 +1486,13 @@ void pbcConstraintElementGroup::directionalLagrangeFormCondition(Iterator itbegi
         break;
       }
     };
-		
+
     if (flag){
 			constraintElement* cele = new directionalFormLagrangeConstraintElement(_solver->getMicroBC(),_LagSpace,_MultSpace,v,vlist,pbcDir);
 			_allConstraint.insert(cele);
     };
-  };																			
-																	
+  };
+
 };
 
 template<class Iterator>
@@ -1522,14 +1522,14 @@ void pbcConstraintElementGroup::directionalLagrangePeriodicCondition(Iterator it
 				}
 			};
 		}
-		
+
     if (flag){
 			SPoint3 pt = v->point();
 			SPoint3 v1 = xlist[0]->point();
 			SPoint3 v2 = xlist[xsize-1]->point();
-			
+
 			SPoint3 pp = planeDirProject(pt,v1,v2,pbcNormal);
-			
+
 			SVector3 ppv1(pp,v1);
 			SVector3 ppv2(pp,v2);
 			if (dot(ppv1,ppv2) <=0.){
@@ -1540,7 +1540,7 @@ void pbcConstraintElementGroup::directionalLagrangePeriodicCondition(Iterator it
 				else{
 					cele = new directionalPBCLagrangeConstraintElement(_solver->getMicroBC(),_LagSpace,_MultSpace,v,xlist,pbcNormal);
 				}
-				 
+
 				_allConstraint.insert(cele);
 			}
 			else{
@@ -1554,7 +1554,7 @@ void pbcConstraintElementGroup::directionalLagrangePeriodicCondition(Iterator it
 				_allConstraint.insert(cele);
 			}
     };
-  };																	 
+  };
 };
 
 // for lagrange
@@ -1832,15 +1832,16 @@ void pbcConstraintElementGroup::createCubicSplineConstraintElementGroup(){
 
 void pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup(){
   nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
-	
+
 	const SVector3& normal = pbc->getShiftPBCNormal();
+	normal.print("normal");
 	if (normal.norm() < 1e-10){
 		printf("no shifted if normal is zero pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup \n");
 		this->createLagrangeConstraintElementGroup();
 		pbc->setShiftedBC(false);
 		return;
 	}
-	
+
 	int dim = _solver->getMicroBC()->getDim();
   const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
  // create base
@@ -1849,26 +1850,26 @@ void pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup(){
   if (dim ==2){
 		this->lagrangeFormCondition(l12.begin(),l12.end(),_xBase);
 		this->lagrangeFormCondition(l41.begin(),l41.end(),_yBase);
-		
+
 		SPoint3 p1 = _v1->point();
 		SPoint3 p2 = _v2->point();
 		SPoint3 p3 = _v3->point();
 		SPoint3 p4 = _v4->point();
 		SVector3 p1p2(p1,p2);
-		
+
 		SPoint3 p = planeDirProject(p3,p1,p2,normal);
 		SVector3 p1p3(p1,p3);
 		SVector3 p2p4(p2,p4);
-		
+
 		double phi0 = acos(dot(p1p2,p1p3)/(p1p2.norm()*p1p3.norm()));
 		double phi1 = acos(dot(p1p2,p2p4)/(p1p2.norm()*p2p4.norm()));
-		
+
 		SVector3 pp3(p,p3);
 		double phi = acos(dot(pp3,p1p2)/(pp3.norm()*p1p2.norm()));
-		
+
 		printf("phi = %f phi0 = %f phi1 = %f \n",phi*180./3.14159265359, phi0*180/3.14159265359,phi1*180./3.14159265359);
-		
-		if ((phi >=phi0) and (phi<=phi1)){				
+
+		if ((phi >=phi0) and (phi<=phi1)){
 			Msg::Info("shifted l34");
 			std::set<MVertex*> others;
 			this->lagrangePeriodicCondition(l23.begin(),l23.end(),_yBase,others);
@@ -1882,7 +1883,7 @@ void pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup(){
 			others.insert(_v3); // avoid duplicating constraint at v3
 			this->shiftedLagrangePeriodicCondition(l23.begin(),l23.end(),_yBase,others,normal);
 		}
-		
+
   }
 	else{
 		Msg::Fatal("createShiftedLagrangeConstraintElementGroup has not been implemented for 3D problems");
@@ -2062,7 +2063,7 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroup(){
 
 		}
 	};
-	
+
 	const std::vector<std::pair<int,int> >& sphysical = mixBC->getStaticPhysical();
 
 	for (int i=0; i<sphysical.size(); i++){
@@ -2087,7 +2088,7 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroup(){
 		}
 
 	}
-	
+
 	// periodic BC
 	const std::vector<std::pair<int,int> >& periodicNegative = mixBC->getNegativePhysical();
 	const std::vector<std::pair<int,int> >& periodicPositive = mixBC->getPositivePhysical();
@@ -2096,39 +2097,39 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroup(){
 		int phyPositive = periodicPositive[i].first;
 		int phyNegative = periodicNegative[i].first;
 		int comp = periodicPositive[i].second;
-		
+
 		Msg::Info("periodic BC negative = %d positive = %d comp = %d",phyNegative,phyPositive,comp);
-		
+
 		groupOfElements* gPositive = _solver->getMicroBC()->getBCGroup(phyPositive);
 		groupOfElements* gNegative = _solver->getMicroBC()->getBCGroup(phyNegative);
-		
+
 		std::set<MVertex*> nullSet;
 		periodicConditionForPeriodicMesh(comp,negFact,gPositive->vbegin(), gPositive->vend(),
 									gNegative->vbegin(), gNegative->vend(),nullSet,NULL,NULL);
-		
+
 	}
-	
+
 };
 
-void pbcConstraintElementGroup::createOrthogonalMixBCConstraintElementGroupDirectionFollowing(){	
+void pbcConstraintElementGroup::createOrthogonalMixBCConstraintElementGroupDirectionFollowing(){
 	nonLinearOrthogonalMixedBCDirectionFollowing* mixBC = dynamic_cast<nonLinearOrthogonalMixedBCDirectionFollowing*>(_solver->getMicroBC());
 	const std::vector<SVector3>& kubcDir = mixBC->getKUBCDirection();
 	const std::vector<SVector3>& subcDir = mixBC->getSUBCDirection();
 
 	const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
 	std::set<MVertex*> allVertex;
-	
+
 	for (int i=0; i< bgroup.size(); i++){
 		for (groupOfElements::vertexContainer::iterator it = bgroup[i]->vbegin(); it!= bgroup[i]->vend(); it++){
 			allVertex.insert(*it);
 		}
 	}
-	
+
 	for (int i=0; i< kubcDir.size(); i++){
 		const SVector3& dir = kubcDir[i];
 		for (std::set<MVertex*>::iterator it = allVertex.begin(); it!= allVertex.end(); it++){
 			MVertex* v = *it;
-			
+
 			if (it == allVertex.begin()){
 				MVertex* v = *(allVertex.begin());
 				constraintElement* ppc = new periodicMeshConstraint(_solver->getMicroBC(), _LagSpace,_LagSpace,_MultSpace,-1,v);
@@ -2140,8 +2141,8 @@ void pbcConstraintElementGroup::createOrthogonalMixBCConstraintElementGroupDirec
 			}
 		}
 	}
-	
-	int halfOfgsize = bgroup.size()/2;	
+
+	int halfOfgsize = bgroup.size()/2;
 	for (int i=0; i< subcDir.size(); i++){
 		const SVector3& dir = subcDir[i];
 		for (int j=0; j< halfOfgsize; j++){
@@ -2153,7 +2154,7 @@ void pbcConstraintElementGroup::createOrthogonalMixBCConstraintElementGroupDirec
 	}
 };
 
-void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollowing(){	
+void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollowing(){
 	nonLinearMixedBCDirectionFollowing* mixBC = dynamic_cast<nonLinearMixedBCDirectionFollowing*>(_solver->getMicroBC());
 	const std::vector<SVector3>& kubcDir = mixBC->getKUBCDirection();
 	const std::vector<SVector3>& subcDir = mixBC->getSUBCDirection();
@@ -2164,18 +2165,18 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollow
 
 	const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
 	std::set<MVertex*> allVertex;
-	
+
 	for (int i=0; i< bgroup.size(); i++){
 		for (groupOfElements::vertexContainer::iterator it = bgroup[i]->vbegin(); it!= bgroup[i]->vend(); it++){
 			allVertex.insert(*it);
 		}
 	}
-	
+
 	for (int i=0; i< kubcDir.size(); i++){
 		const SVector3& dir = kubcDir[i];
 		for (std::set<MVertex*>::iterator it = allVertex.begin(); it!= allVertex.end(); it++){
 			MVertex* v = *it;
-			
+
 			if ((pbcNormalDir.size() == 0) and (it == allVertex.begin())){
 				MVertex* v = *(allVertex.begin());
 				constraintElement* ppc = new periodicMeshConstraint(_solver->getMicroBC(), _LagSpace,_LagSpace,_MultSpace,-1,v);
@@ -2190,7 +2191,7 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollow
 	if (pbcNormalDir.size() > 0){
 		const std::vector<int>& rootPhysical = mixBC->getRootPhysical();
 		const std::vector<int>& interpDeg = mixBC->getInterpolationDegree();
-		
+
 		if (mixBC->getDim() == 2){
 			MVertex* v1(NULL), *v2(NULL), *v4(NULL);
 			for (std::set<MVertex*>::iterator it = allVertex.begin(); it!= allVertex.end(); it++){
@@ -2211,22 +2212,22 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollow
 			if ((v1==NULL) or (v2 == NULL) or (v4 ==NULL)){
 				Msg::Fatal("root physical vertces are not correctly defined");
 			}
-			
+
 			constraintElement* ppc = new periodicMeshConstraint(_solver->getMicroBC(), _LagSpace,_LagSpace,_MultSpace,-1,v2);
 			_allConstraint.insert(ppc);
-			
+
 			// create local basis
 			this->getInterpolatioPointsFromExistingVertices(interpDeg[0],v1,v4,bgroup[0]->vbegin(),bgroup[0]->vend(),_xBase);
 			this->getInterpolatioPointsFromExistingVertices(interpDeg[1],v1,v2,bgroup[1]->vbegin(),bgroup[1]->vend(),_yBase);
-			
+
 			for (int i=0; i< _xBase.size(); i++){
 				Msg::Info("xbase :  %d",_xBase[i]->getNum());
 			}
-			
+
 			for (int i=0; i< _yBase.size(); i++){
 				Msg::Info("ybase :  %d",_yBase[i]->getNum());
 			}
-			
+
 			std::set<MVertex*> allPositiveVertex;
 			for (groupOfElements::vertexContainer::iterator it = bgroup[2]->vbegin(); it!= bgroup[2]->vend(); it++){
 				allPositiveVertex.insert(*it);
@@ -2234,17 +2235,17 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollow
 			for (groupOfElements::vertexContainer::iterator it = bgroup[3]->vbegin(); it!= bgroup[3]->vend(); it++){
 				allPositiveVertex.insert(*it);
 			}
-			
-			
+
+
 			for (int i=0; i< pbcNormalDir.size(); i++){
 				const SVector3& normalPBC = pbcNormalDir[i];
 				const SVector3& dirPBC = pbcDir[i];
 				const bool& fullCon = pbcFullConstrainedFlag[i];
-				
+
 				if (fullCon){
 					this->lagrangeFormCondition(bgroup[0]->vbegin(),bgroup[0]->vend(),_xBase);
 					this->lagrangeFormCondition(bgroup[1]->vbegin(),bgroup[1]->vend(),_yBase);
-					this->directionalLagrangePeriodicCondition(allPositiveVertex.begin(), allPositiveVertex.end(),_xBase,_yBase,normalPBC,true);					
+					this->directionalLagrangePeriodicCondition(allPositiveVertex.begin(), allPositiveVertex.end(),_xBase,_yBase,normalPBC,true);
 				}
 				else{
 					this->directionalLagrangeFormCondition(bgroup[0]->vbegin(),bgroup[0]->vend(),_xBase,dirPBC);
@@ -2252,14 +2253,14 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollow
 					this->directionalLagrangePeriodicCondition(allPositiveVertex.begin(), allPositiveVertex.end(),_xBase,_yBase,normalPBC,false);
 				}
 			}
-		
+
 		}
 		else{
 			Msg::Fatal("this bc is not implemented for 3D problems");
 		}
 	}
-	
-	int halfOfgsize = bgroup.size()/2;	
+
+	int halfOfgsize = bgroup.size()/2;
 	for (int i=0; i< subcDir.size(); i++){
 		const SVector3& dir = subcDir[i];
 		for (int j=0; j< halfOfgsize; j++){
@@ -2358,7 +2359,7 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroup_DG(){
 		}
 
 	}
-	
+
 	if (mixBC->getPositivePhysical().size() > 0){
 		Msg::Error("periodic BC on components is not considered ");
 	}
@@ -2368,7 +2369,7 @@ void pbcConstraintElementGroup::createCornerConstraintElementGroupShiftedBPC(){
 	if (_solver->getMicroBC()->getOrder() == 2){
 		Msg::Fatal("shifted BC is implemented for first order BC only");
 	}
-	
+
 	if (_cVertices.size() > 0){
 		constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,-1,_v1);
 		_allConstraint.insert(cel);
@@ -2739,7 +2740,7 @@ void pbcConstraintElementGroup::createConstraintGroup(){
 			nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
 			if (pbc->isShifted() and pbc->getPBCMethod() != nonLinearPeriodicBC::LIM){
 				Msg::Warning("shifted has not implemented for this option");
-			} 
+			}
 			const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
 			if (bgroup.size()>0){
 				if  (pbc->getPBCMethod() == nonLinearPeriodicBC::CEM){
@@ -2775,7 +2776,7 @@ void pbcConstraintElementGroup::createConstraintGroup(){
 				else{
 					Msg::Fatal("this method is not implemented to impose PBC");
 				}
-				
+
 				if (pbc->isShifted()){
 					this->createCornerConstraintElementGroupShiftedBPC();
 				}
@@ -3194,14 +3195,14 @@ SPoint3 pbcConstraintElementGroup::getRootPoint(){
 				SPoint3 firstVn = (*it)->point();
 				it++;
 				SPoint3 secondVn = (*it)->point();
-				
+
 				SPoint3 p = project(oneVp,firstVn,secondVn);
-				_v1 = new MVertex(p.x(),p.y(),p.z());				
+				_v1 = new MVertex(p.x(),p.y(),p.z());
 			}
 			return _v1->point();
 		}
 		else if (_solver->getMicroBC()->getDim() == 3){
-			
+
 			MElement* e1 = *(bgroup[0]->begin());
       MElement* e2 = *(bgroup[2]->begin());
       MElement* e3 = *(bgroup[1]->begin());
@@ -3221,7 +3222,7 @@ SPoint3 pbcConstraintElementGroup::getRootPoint(){
       m.invert(invm);
       fullVector<double> f(3);
       invm.mult(b,f);
-			
+
 			SPoint3 p(f(0),f(1),f(2));
 			_v1 = new MVertex(p.x(),p.y(),p.z());
 			return p;
@@ -3230,6 +3231,6 @@ SPoint3 pbcConstraintElementGroup::getRootPoint(){
 			Msg::Fatal("getRootPoint has not been implemented for %d D problems",_solver->getMicroBC()->getDim());
 		}
 	}
-	
+
 };