diff --git a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
index 88dd0aba71b752382b8b2fec93b468137b58d146..026967d09880d02209113ae2ddbbd15c573f444e 100644
--- a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
+++ b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
@@ -348,7 +348,7 @@ void nonLinearMicroBC::setPathFunctionDeformationGradient(int compi, int compj,
    _pathFunctionDefoGradientF[index] = fct.clone();
 };
 
-void nonLinearMicroBC::setPathFunctionDeformationGradient(int compi, int compj, int compk, const scalarFunction& fct){
+void nonLinearMicroBC::setPathFunctionGradientOfDeformationGradient(int compi, int compj, int compk, const scalarFunction& fct){
    int index=Tensor33::getIndex(compi,compj,compk);
    delete _pathFunctionDefoGradientG[index];
    _pathFunctionDefoGradientG[index] = fct.clone();
diff --git a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
index 639f81c6ef169448cb1f40e0ef4c5b7410a0e985..3f3c5e0e7697e051cb365877b6ace38e95308917 100644
--- a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
+++ b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
@@ -95,7 +95,7 @@ public:
     
     void setPathFunctionDeformationGradient(const scalarFunction& fct);
     void setPathFunctionDeformationGradient(int compi, int compj, const scalarFunction& fct);
-    void setPathFunctionDeformationGradient(int compi, int compj, int compk, const scalarFunction& fct);
+    void setPathFunctionGradientOfDeformationGradient(int compi, int compj, int compk, const scalarFunction& fct);
 		
 		// mechanics
 		void setDeformationGradient(const int i,const int j, const double val);
@@ -128,7 +128,7 @@ public:
     const STensor3& getPreviousDeformationGradient() const {return Fprev;};
     const STensor3& getCurrentDeformationGradient() const{return Fcur;};
     const STensor33& getPreviousDeformationGradientGradient() const {return Gprev;};
-    const STensor33& getCurrentDeformationGradientGradient() const{return G;};
+    const STensor33& getCurrentDeformationGradientGradient() const{return Gcur;};
 		
 		void setDeformationGradient(const STensor3& F);
 		void setGradientOfDeformationGradient(const STensor33& G);
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 495b8b7425f5fbd97991b41bf44f1f9274dc64fd..7a2d41a87d5117c413785dddbb9be4ae72c47b73 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -5847,7 +5847,7 @@ void nonLinearMechSolver::computeBodyForceVector(){
     this->getMicroBC()->setTime(1.); 
     G1 = this->getMicroBC()->getSecondOrderKinematicalVariable();
     this->getMicroBC()->setTime(curtime);
-       
+           
     G1 -=G0;
     for (int idom=0; idom<domainVector.size(); idom++){
         partDomain* dom = domainVector[idom];
@@ -9373,7 +9373,7 @@ double nonLinearMechSolver::computeRightHandSide()
     if (_systemType == nonLinearMechSolver::MULT_ELIM){
       _pAl->assembleConstraintResidualToSystem();
       this->computeStiffMatrix();
-      if(computedUdF() and _pathFollowing){
+      if(computedUdF() and _pathFollowing and usePresentModuli()){
          this->computeBodyForceVector();
       }
     }
@@ -9616,7 +9616,7 @@ int nonLinearMechSolver::NewtonRaphson(const int numstep){
 
   // compute ipvariable
   if(computedUdF()){
-      _ipf->setDeformationGradientGradient(this->getMicroBC()->getCurrentDeformationGradientGradient(), IPStateBase::current);
+      _ipf->setDeformationGradientGradient(this->getMicroBC()->getSecondOrderKinematicalVariable(), IPStateBase::current);
   }
   _ipf->compute1state(IPStateBase::current, true);
   // Compute Right Hand Side (Fext-Fint)
@@ -13227,8 +13227,6 @@ void nonLinearMechSolver::extractdUdFByInSystemCondensation(homogenizedData* hom
   t = Cpu() -t;
   if (_multiscaleFlag == false)
     Msg::Info("Done performing in-system condensation (%f s)",t);
-    
-  //_ipf->setdFmdFM(homoData->getdUnknowndF(), IPStateBase::current);
 };
 
 void nonLinearMechSolver::extractAverageStress(homogenizedData* homoData){
@@ -14588,6 +14586,11 @@ void nonLinearMechSolver::nextStep(const double curtime, const int step)
   }
   _ipf->nextStep(curtime); //ipvariable
   pAssembler->nextStep();
+  // update for bodyforce 
+  if(computedUdF()){
+     this->extractdUdFByInSystemCondensation(this->getHomogenizationState(IPStateBase::current));
+  }
+  
   // contact
   for(contactContainer::iterator it = _allContact.begin(); it!=_allContact.end(); ++it){
     contactDomain *cdom = *it;
@@ -15207,7 +15210,6 @@ double nonLinearMechSolver::solveMicroSolverSNL(){
   //
   _sucessMicroSolve = false;
   bool willFinish = false;
-  bool dUdF_Flag = false;
   
   while(!_timeManager->reachEndTime())
   {
@@ -15227,24 +15229,21 @@ double nonLinearMechSolver::solveMicroSolverSNL(){
     pAssembler->setTimeStep(dt);
     _ipf->setTime(curtime,dt,ii+1);
     if(computedUdF()){
-      if(dUdF_Flag){
-         homogenizedData* homoData = this->getHomogenizationState(IPStateBase::current);
-         this->extractdUdFByInSystemCondensation(homoData);
-         _ipf->setdFmdFM(homoData->getdUnknowndF(), IPStateBase::current);
-     }
-     else{
-         dUdF_Flag = true;
-     } 
+       _ipf->setdFmdFM(this->getHomogenizationState(IPStateBase::current)->getdUnknowndF(), IPStateBase::current);
+       
     }
     
     // Solve one step by NR scheme
     int niteNR = 0;
     if (_pathFollowing)
     {
-
       pathSys->setPathFollowingIncrement(dt);
       _pAl->updateConstraint(_ipf);
       
+      if(computedUdF() and !usePresentModuli()){
+         this->computeBodyForceVector();
+      }
+   
       niteNR = microNewtonRaphsonPathFollowing(ii+1);
 
         
@@ -15416,7 +15415,7 @@ double nonLinearMechSolver::solveMicroSolverStaticLinear(){
   this->setTimeForBC(curtime);
   this->setTimeForLaw(curtime,timestep,step);
   if(computedUdF()){
-     _ipf->setDeformationGradientGradient(this->getMicroBC()->getCurrentDeformationGradientGradient(), IPStateBase::current);
+     _ipf->setDeformationGradientGradient(this->getMicroBC()->getSecondOrderKinematicalVariable(), IPStateBase::current);
   }
 
   if (_systemType == DISP_ELIM)
@@ -15446,7 +15445,7 @@ double nonLinearMechSolver::solveMicroSolverForwardEuler(){
   //
   this->setTimeForBC(1.);
   if(computedUdF()){
-      _ipf->setDeformationGradientGradient(this->getMicroBC()->getCurrentDeformationGradientGradient(), IPStateBase::current);
+      _ipf->setDeformationGradientGradient(this->getMicroBC()->getSecondOrderKinematicalVariable(), IPStateBase::current);
   }
   _pAl->assembleRightHandSide();
   //
@@ -15478,7 +15477,7 @@ double nonLinearMechSolver::solveMicroSolverForwardEuler(){
     // update ipvariable
 
     if(computedUdF()){
-      _ipf->setDeformationGradientGradient(this->getMicroBC()->getCurrentDeformationGradientGradient(), IPStateBase::current);
+      _ipf->setDeformationGradientGradient(this->getMicroBC()->getSecondOrderKinematicalVariable(), IPStateBase::current);
     }
     _ipf->compute1state(IPStateBase::current,true);
 
@@ -15491,7 +15490,7 @@ double nonLinearMechSolver::solveMicroSolverForwardEuler(){
     {
       this->nextStep(control,lastIter+1);
     }
-		lsys->zeroMatrix();
+    lsys->zeroMatrix();
     elsys->nextStep();
     _timeManager->computeTimeStepForNextSolving(1);
   }
@@ -15622,7 +15621,7 @@ int nonLinearMechSolver::microNewtonRaphson(const int numstep){
   // compute ipvariable
 
   if(computedUdF()){
-      _ipf->setDeformationGradientGradient(this->getMicroBC()->getCurrentDeformationGradientGradient(), IPStateBase::current);
+      _ipf->setDeformationGradientGradient(this->getMicroBC()->getSecondOrderKinematicalVariable(), IPStateBase::current);
   }
   _ipf->compute1state(IPStateBase::current, true);
   double normFinfInit = this->computeRightHandSide();
diff --git a/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h b/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
index 8cc34fb380dd635f092efd524c27ce1fb174bc21..b1d1279eec94150dff847dd32e6ff416dceaff58 100644
--- a/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
+++ b/NonLinearSolver/pathFollowing/pbcPathFollowingSystemPETSc.h
@@ -234,14 +234,7 @@ class pbcPathFollowingSystemPETSC : public pbcNonLinearSystemPETSc<scalar>,
       _try(VecScale(_qeff,this->_scale));
       
       if (_isAllocatedBodyForceVector){
-       //  if(_iter==0){
-       //    _try(MatMultAdd(this->_Q,_BF,_qeff,_qeff));
-       //  }
-       //  else{
-           PetscScalar inv_dt = 1./_pseudoTimeIncrement;
-           _try(VecScale(_BF, inv_dt));
-           _try(MatMultAdd(this->_Q,_BF,_qeff,_qeff));
-      //   }
+         _try(MatMultAdd(this->_Q,_BF,_qeff,_qeff));
        }
       
       _try(VecAssemblyBegin(_qeff));
diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp
index e5f71e42bc5f6a1bce875a67c84de1833cdf2915..1d5c72e6e91353f5a556a19ef7631267d92d3d3f 100644
--- a/dG3D/src/dG3DDomain.cpp
+++ b/dG3D/src/dG3DDomain.cpp
@@ -2606,10 +2606,9 @@ void dG3DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichStat
   bool needDTangent=false;
   bool usePresentModuli = false;
   
-  if (hasBodyForceForHO()){
-         needDTangent=true;
-    if(usePresentModuliHO())
-         usePresentModuli = true;
+  if (hasBodyForceForHO() and usePresentModuliHO()){
+      needDTangent= true;
+      usePresentModuli = true;
   }      
   if(_evalStiff) { this->computeStrain(e,npts_bulk,GP,vips,ws,disp,useBarF);};
   for(int j=0;j<npts_bulk;j++)
@@ -2626,7 +2625,7 @@ void dG3DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichStat
       else
       	mlaw->stress(ipv,ipvprev, stiff,true,needDTangent);
       	
-      if(needDTangent){
+      if(hasBodyForceForHO()){
         if(usePresentModuli){  
           ipv->setValueForBodyForce(ipv->getConstRefToTangentModuli(),usePresentModuli); 
         }
@@ -2693,7 +2692,7 @@ void dG3DDomain::computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichStat
         // next step
         iter++;
       }
-      if(needDTangent){
+      if(hasBodyForceForHO()){
         if(usePresentModuli){
             ipv->setValueForBodyForce(ipv->getConstRefToTangentModuli(),usePresentModuli); 
         }    
@@ -3725,9 +3724,8 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
     bool needDTangent=false;
     bool usePresentModuli = false;
   
-    if (hasBodyForceForHO()){
+    if (hasBodyForceForHO() and usePresentModuliHO()){
          needDTangent=true;
-      if(usePresentModuliHO())
          usePresentModuli = true;
     } 
     AllIPState::ipstateElementContainer *vips = aips->getIPstate(ele->getNum());
@@ -3902,7 +3900,7 @@ void dG3DDomain::computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,co
       else{
         Msg::Error("This substepping method %d has not been implemented",_subSteppingMethod);
       }
-      if(needDTangent){
+      if(hasBodyForceForHO()){
         if(usePresentModuli){  
           ipvp->setValueForBodyForce(ipvp->getConstRefToTangentModuli(),usePresentModuli);
           ipvm->setValueForBodyForce(ipvm->getConstRefToTangentModuli(),usePresentModuli);
@@ -5522,9 +5520,8 @@ void dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufiel
   bool needDTangent=false;
   bool usePresentModuli = false;
   
-  if (hasBodyForceForHO()){
-         needDTangent=true;
-    if(usePresentModuliHO())
+  if (hasBodyForceForHO() and usePresentModuliHO()){
+         needDTangent=true;    
          usePresentModuli = true;
   } 
   for (std::set<int>::iterator it = _domainIPBulk.begin(); it != _domainIPBulk.end(); it++){
@@ -5547,7 +5544,19 @@ void dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufiel
       mlaw->stress3DTo2D(ipv,ipvprev,stiff,false,needDTangent);
     else
       mlaw->stress(ipv,ipvprev,stiff,false,needDTangent);
+      
+    if(hasBodyForceForHO()){
+      if(usePresentModuli){  
+        ipv->setValueForBodyForce(ipv->getConstRefToTangentModuli(),usePresentModuli); 
+      }
+      else{
+        ipv->setValueForBodyForce(ipvprev->getConstRefToTangentModuli(),usePresentModuli);  
+      } 
+    }  
   };
+  
+  
+  
 
   for (std::set<int>::iterator it = _domainIPInterfaceMinus.begin(); it != _domainIPInterfaceMinus.end(); it++){
     int num = *it;
@@ -5565,10 +5574,18 @@ void dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufiel
     IPVariable *ipv = ips->getState(ws);
     const IPVariable *ipvprev = ips->getState(IPStateBase::previous);
     dG3DMaterialLaw *mlawMinus = dynamic_cast<dG3DMaterialLaw*>(this->getMaterialLawMinus());
-		if (_planeStressState)
-			mlawMinus->stress3DTo2D(ipv,ipvprev,stiff,true,needDTangent);
-		else
+    if (_planeStressState)
+	mlawMinus->stress3DTo2D(ipv,ipvprev,stiff,true,needDTangent);
+    else
     	mlawMinus->stress(ipv,ipvprev,stiff,true,needDTangent);
+    if(hasBodyForceForHO()){
+      if(usePresentModuli){  
+        ipv->setValueForBodyForce(ipv->getConstRefToTangentModuli(),usePresentModuli); 
+      }
+      else{
+        ipv->setValueForBodyForce(ipvprev->getConstRefToTangentModuli(),usePresentModuli);  
+      } 
+    }  
   };
 
 
@@ -5588,10 +5605,18 @@ void dG3DDomain::computeIPVariableMPI(AllIPState *aips,const unknownField *ufiel
     IPVariable *ipv = ips->getState(ws);
     const IPVariable *ipvprev = ips->getState(IPStateBase::previous);
     dG3DMaterialLaw *mlawPlus= dynamic_cast<dG3DMaterialLaw*>(this->getMaterialLawPlus());
-		if (_planeStressState)
-			mlawPlus->stress3DTo2D(ipv,ipvprev,stiff,true,needDTangent);
-		else
-    	mlawPlus->stress(ipv,ipvprev,stiff,true,needDTangent);
+    if (_planeStressState)
+	mlawPlus->stress3DTo2D(ipv,ipvprev,stiff,true,needDTangent);
+    else
+       mlawPlus->stress(ipv,ipvprev,stiff,true,needDTangent);
+    if(hasBodyForceForHO()){
+      if(usePresentModuli){  
+        ipv->setValueForBodyForce(ipv->getConstRefToTangentModuli(),usePresentModuli); 
+      }
+      else{
+        ipv->setValueForBodyForce(ipvprev->getConstRefToTangentModuli(),usePresentModuli);  
+      } 
+    }  
   }
 
 
diff --git a/dG3D/src/dG3DHomogenizedTangentTerms.cpp b/dG3D/src/dG3DHomogenizedTangentTerms.cpp
index 79eb95eef755a0c69b8d436c8c4a606baf0cd5e0..eba072b5d973c9793cab1405ce41d45937fef6b5 100644
--- a/dG3D/src/dG3DHomogenizedTangentTerms.cpp
+++ b/dG3D/src/dG3DHomogenizedTangentTerms.cpp
@@ -454,18 +454,14 @@ void dG3DBodyForceTermPathFollowing::get(MElement *ele,int npts,IntPt *GP,fullVe
     double &detJ = ipv->getJacobianDeterminant(ele,GP[i]);
 		
     const double ratio = detJ*weight;
-  //  SVector3 Bm;
-  //  ipv->computeBodyForce(G, Bm);
-    const SVector3 &Bm0 = ipvprev->getConstRefToBm();
-    const SVector3 &Bm = ipv->getConstRefToBm();
-  
-    
+    SVector3 Bm;
+    ipv->computeBodyForce(G, Bm);
       
     if (_dom->getHomogenizationOrder() == 2 and ipv->hasBodyForceForHO())
     {
        for (int j=0; j<nbFF; j++)
           for (int k=0; k<3; k++){
-              bm(j+k*nbFF)+= (Bm(k)-Bm0(k))*Vals[j+k*nbFF]*ratio;
+              bm(j+k*nbFF)+= Bm(k)*Vals[j+k*nbFF]*ratio;
           }           
     }	
   }
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 016944b880bd11b77ff359cd264bc8b04bd42e05..0e1ff081a7c469c44b5d242361b9422eaa977116 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -1176,7 +1176,6 @@ void dG3DIPVariable::computeBodyForce(const STensor33& G, SVector3& B) const{
 
 
 void dG3DIPVariable::setValueForBodyForce(const STensor43& K, const bool usePresentModuli){
-  // STensor43& K = elasticTangentModuli;
 
  if(hasBodyForceForHO()){
    static STensor33 tmp1;