diff --git a/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp b/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp
index c9b2b94e06a46fa9ebbf8471a2af11548e2f5cfd..ca01520d3df82d5bd03a188f595938cdde69f4ef 100644
--- a/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp
+++ b/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp
@@ -468,7 +468,6 @@ void fullMatrixOperation::convertMat66ToMat99(const fullMatrix<double>& mat6x6,
     else if ((i==1 && j==2) || (i==2 && j==1)) return 5;
     {
       Msg::Error("index %d %d does not exst",i,j);
-			return -1;
     }
   };
   for (int i=0; i< 3; i++)
@@ -722,7 +721,7 @@ void BimaterialHomogenization::computeBimaterial(const fullMatrix<double>& C1, c
                         fullMatrix<double>* DCeffDnormal1,fullMatrix<double>* DCeffDnormal2, fullMatrix<double>* DCeffDnormal3) const
 {
   double tol = 1e-10;
-  if ((C1.norm() < tol*C2.norm()))
+  if (C1.norm() < tol*C2.norm())
   {
     // C1 is void
     Ceff  = C2;
@@ -743,7 +742,7 @@ void BimaterialHomogenization::computeBimaterial(const fullMatrix<double>& C1, c
       fullMatrixOperation::allocateAndMakeZero(*DCeffDnormal3,Ceff.size1(),Ceff.size2());
     };
   }
-  else if ((C2.norm() < tol*C1.norm()))
+  else if (C2.norm() < tol*C1.norm())
   {
     Ceff = C1;
     Ceff.scale(f1);
@@ -785,10 +784,6 @@ void BimaterialHomogenization::computeBimaterial(const fullMatrix<double>& C1, c
     bool ok = fullMatrixOperation::invertMatrix(A,invA,stiff,DinvADA);
     if (!ok)
     {
-      C1.print("C1");
-      C2.print("C2");
-      normal.print("normal");
-      printf("f1 = %e, f2 = %e\n",f1,f2);
       A.print("A");
     }
     //
@@ -927,30 +922,6 @@ void BimaterialHomogenization::computeBimaterial(const fullMatrix<double>& C1, c
   };
 };
 
-void BimaterialHomogenization::compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, 
-                        const std::vector<SVector3>& normal,
-                        fullMatrix<double>& Ceff, 
-                        bool stiff, 
-                        std::vector<hyperFullMatrix>& DCeffDcin, 
-                        std::vector<fullMatrix<double> >& DCeffDf, 
-                        std::vector<std::vector<fullMatrix<double> > >&DCeffDnormal) const
-{
-
-  std::vector<fullMatrix<double> >  Dvec;
-  compute(Cin,f,normal[0],Ceff,stiff,DCeffDcin,DCeffDf,Dvec);
-  if (stiff)
-  {
-    if (DCeffDnormal.size() != normal.size())
-    {
-      DCeffDnormal.resize(normal.size());
-    }
-    for (int i=0; i< normal.size(); i++)
-    {
-      DCeffDnormal[i] = Dvec;
-    }
-  }
-};
-
 
 
 void BimaterialHomogenization::compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, const SVector3& normal,
@@ -980,14 +951,9 @@ void BimaterialHomogenization::compute(const std::vector<fullMatrix<double> >& C
     {
       DCeffDnormal.resize(3);
     }    
-    computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal,Ceff,
-                    stiff,&DCeffDcin[0],&DCeffDcin[1],&DCeffDf[0],&DCeffDf[1],&DCeffDnormal[0],&DCeffDnormal[1],&DCeffDnormal[2]);
-  }
-  else
-  {
-    computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal,Ceff);
   };
-  
+  computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal,Ceff,
+                    stiff,&DCeffDcin[0],&DCeffDcin[1],&DCeffDf[0],&DCeffDf[1],&DCeffDnormal[0],&DCeffDnormal[1],&DCeffDnormal[2]);
 };
 
 void LaminateHomogenization::compute(const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& C3, 
@@ -1110,17 +1076,10 @@ void LaminateHomogenization::compute(const std::vector<fullMatrix<double> >& Cin
   }
   else if (numPhase == 2)
   {
-    if (stiff)
-    {
-      computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal,Ceff,
-                      stiff,&DCeffDcin[0],&DCeffDcin[1],
-                      &DCeffDf[0],&DCeffDf[1],&DCeffDnormal[0],
-                      &DCeffDnormal[1],&DCeffDnormal[2]);
-    }
-    else
-    {
-      computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal,Ceff);
-    }
+    computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal,Ceff,
+                    stiff,&DCeffDcin[0],&DCeffDcin[1],
+                    &DCeffDf[0],&DCeffDf[1],&DCeffDnormal[0],
+                    &DCeffDnormal[1],&DCeffDnormal[2]);
   }
   else
   {
@@ -1185,14 +1144,7 @@ void LaminateHomogenization::compute(const std::vector<fullMatrix<double> >& Cin
         };
       }
       // two phase
-      if (stiff)
-      {
-        computeBimaterial(CC[i],Cin[i+1],f0i,f1i,normal,CC[i+1],stiff,&A[i],&B[i],&F0[i],&F1[i],&N0[i],&N1[i],&N2[i]);
-      }
-      else
-      {
-        computeBimaterial(CC[i],Cin[i+1],f0i,f1i,normal,CC[i+1]);
-      }
+      computeBimaterial(CC[i],Cin[i+1],f0i,f1i,normal,CC[i+1],stiff,&A[i],&B[i],&F0[i],&F1[i],&N0[i],&N1[i],&N2[i]);
     };
     //
     Ceff = CC[numPhase-1];
@@ -1239,184 +1191,7 @@ void LaminateHomogenization::compute(const std::vector<fullMatrix<double> >& Cin
       DCeffDcin[0]=AA;
     }
   }
-};
-
-
-void LaminateHomogenization::compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, 
-                        const std::vector<SVector3>& normal,
-                        fullMatrix<double>& Ceff, 
-                        bool stiff, 
-                        std::vector<hyperFullMatrix>& DCeffDcin, 
-                        std::vector<fullMatrix<double> >& DCeffDf, 
-                        std::vector<std::vector<fullMatrix<double> > >&DCeffDnormal) const
-{
-  int numPhase = Cin.size();
-  if (stiff)
-  {
-    if (DCeffDcin.size() != numPhase)
-    {
-      DCeffDcin.resize(numPhase);
-    }
-    if (DCeffDf.size() != numPhase)
-    {
-      DCeffDf.resize(numPhase);
-    }
-    if (DCeffDnormal.size() != normal.size())
-    {
-      DCeffDnormal.resize(normal.size(), std::vector<fullMatrix<double> >(3));
-    }
-  };
-  
-  if (numPhase ==1) 
-  {
-    Ceff = Cin[0];
-    if (stiff)
-    {
-      DCeffDcin[0].allocate(Ceff.size1(),Cin[0].size1(),1.);
-      fullMatrixOperation::allocateAndMakeZero(DCeffDf[0],Cin[0].size1(),Cin[0].size2());
-      for (int j=0; j< normal.size(); j++)
-      {
-        fullMatrixOperation::allocateAndMakeZero(DCeffDnormal[j][0],Cin[0].size1(),Cin[0].size2());
-        fullMatrixOperation::allocateAndMakeZero(DCeffDnormal[j][1],Cin[0].size1(),Cin[0].size2());
-        fullMatrixOperation::allocateAndMakeZero(DCeffDnormal[j][2],Cin[0].size1(),Cin[0].size2());        
-      }
-    }
-  }
-  else if (numPhase == 2)
-  {
-    if (stiff)
-    {
-      computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal[0],Ceff,
-                    stiff,&DCeffDcin[0],&DCeffDcin[1],
-                    &DCeffDf[0],&DCeffDf[1],&DCeffDnormal[0][0],
-                    &DCeffDnormal[0][1],&DCeffDnormal[0][2]);
-    }
-    else
-    {
-      computeBimaterial(Cin[0],Cin[1],f[0],f[1],normal[0],Ceff);
-    }
-  }
-  else
-  {
-    //
-    auto getSumFaction = [](const std::vector<double>& ff, int i, fullVector<double>& DoutDf)
-    {
-      DoutDf.resize(ff.size());
-      DoutDf.setAll(0.);
-      double v= 0;
-      for (int k=0; k< i+1; k++)
-      {
-        v += ff[k];
-        DoutDf(k) = 1.;
-      }
-      return v;
-    };
-    
-    static std::vector<hyperFullMatrix> A, B;
-    static std::vector<fullMatrix<double> > CC, F0, F1, N0, N1, N2;
-    static fullMatrix<double> Df0iDf, Df1iDf;
-    
-    if (CC.size() != numPhase)
-    {
-      CC.resize(numPhase);
-    };
-    //
-    if (stiff)
-    {
-      // resallocate for static variable
-      if (A.size() != numPhase-1) A.resize(numPhase-1);
-      if (B.size() != numPhase-1) B.resize(numPhase-1);
-      if (F0.size() != numPhase-1) F0.resize(numPhase-1);
-      if (F1.size() != numPhase-1) F1.resize(numPhase-1);
-      if (N0.size() != numPhase-1) N0.resize(numPhase-1);
-      if (N1.size() != numPhase-1) N1.resize(numPhase-1);
-      if (N2.size() != numPhase-1) N2.resize(numPhase-1);
-      
-      fullMatrixOperation::allocateAndMakeZero(Df0iDf,numPhase-1,numPhase);
-      fullMatrixOperation::allocateAndMakeZero(Df1iDf,numPhase-1,numPhase);
-      
-      for (int j=0; j< numPhase; j++)
-      {
-        fullMatrixOperation::allocateAndMakeZero(DCeffDf[j],Cin[0].size1(),Cin[0].size2());
-      }
-    }
-    //
-    CC[0] = Cin[0];
-    for (int i=0; i< numPhase-1; i++)
-    {
-      static fullVector<double> DftolDf, DftolPrevDf;
-      double ftol = getSumFaction(f,i+1,DftolDf);
-      double ftolPrev = getSumFaction(f,i,DftolPrevDf);
-      
-      double f0i = ftolPrev/ftol;
-      double f1i = 1.-f0i;
-      if (stiff)
-      {     
-        for (int r=0; r< numPhase; r++)
-        {
-          Df0iDf(i,r) = (-ftolPrev/ftol/ftol)*DftolDf(r) + (1./ftol)*DftolPrevDf(r);
-          Df1iDf(i,r) = -Df0iDf(i,r);
-        };
-      }
-      // two phase
-      if (stiff)
-      {
-        computeBimaterial(CC[i],Cin[i+1],f0i,f1i,normal[i],CC[i+1],stiff,&A[i],&B[i],&F0[i],&F1[i],&N0[i],&N1[i],&N2[i]);
-      }
-      else
-      {
-        computeBimaterial(CC[i],Cin[i+1],f0i,f1i,normal[i],CC[i+1]);
-      }
-    };
-    //
-    Ceff = CC[numPhase-1];
-    //
-    if (stiff)
-    {
-      
-      DCeffDcin[numPhase-1] = B[numPhase-2];
-      //
-      for (int j=0; j< numPhase; j++)
-      {
-        DCeffDf[j].axpy(F0[numPhase-2],Df0iDf(numPhase-2,j));
-        DCeffDf[j].axpy(F1[numPhase-2],Df1iDf(numPhase-2,j));
-      }
-      //
-      DCeffDnormal[numPhase-2][0] = N0[numPhase-2];
-      DCeffDnormal[numPhase-2][1] = N1[numPhase-2];
-      DCeffDnormal[numPhase-2][2] = N2[numPhase-2];
-      
-      static hyperFullMatrix AA;
-      AA = A[numPhase-2];
-      for (int i=numPhase-2; i>0; i--)
-      {        
-        fullMatrixOperation::mult(AA,B[i-1],DCeffDcin[i],1.);
-        
-        static fullMatrix<double> AAF0, AAF1;
-        fullMatrixOperation::mult(AA,F0[i-1],AAF0);
-        fullMatrixOperation::mult(AA,F1[i-1],AAF1);
-        //
-        for (int j=0; j< numPhase; j++)
-        {
-          DCeffDf[j].axpy(AAF0,Df0iDf(i-1,j));
-          DCeffDf[j].axpy(AAF1,Df1iDf(i-1,j));
-        }
-        //
-        
-        fullMatrixOperation::mult(AA,N0[i-1],DCeffDnormal[i-1][0]);
-        fullMatrixOperation::mult(AA,N1[i-1],DCeffDnormal[i-1][1]);
-        fullMatrixOperation::mult(AA,N2[i-1],DCeffDnormal[i-1][2]);
-        //
-        hyperFullMatrix temp;
-        fullMatrixOperation::mult(AA,A[i-1],temp,1.);
-        AA = temp;
-      }
-      DCeffDcin[0]=AA;
-    }
-  }
-};
-                        
-                        
+}
 void LaminateHomogenization::printInfos() const
 {
   printf("using Voigt, Reuss with laminate \n");
@@ -1437,7 +1212,6 @@ LossMeasure* LossMeasure::createLossMeasure(const char what[])
   else
   {
     Msg::Error("%s has not been implemented",what);
-		return NULL;
   }
 };
 
@@ -1497,7 +1271,7 @@ double SquareRootError::get(const fullMatrix<double>& Cref, const fullMatrix<dou
   return sqrt(y/yref);
 };
 
-TrainingDeepMaterialNetwork::TrainingDeepMaterialNetwork(Tree& T, const Homogenization& ho): _T(&T), _homo(&ho),_lambda(0.), _sameNormal(true)
+TrainingDeepMaterialNetwork::TrainingDeepMaterialNetwork(Tree& T, const Homogenization& ho): _T(&T), _homo(&ho),_lambda(0.)
 {
   Tree::nodeContainer allNodes;
   _T->getAllNodes(allNodes);
@@ -1516,19 +1290,6 @@ TrainingDeepMaterialNetwork::~TrainingDeepMaterialNetwork()
   _offlineData.clear();
 }
 
-void TrainingDeepMaterialNetwork::sameNormal(bool fl)
-{
-  _sameNormal = fl;
-  if (fl)
-  {
-    Msg::Info("same normal is used");
-  }
-  else
-  {
-    Msg::Info("normal is different at different laminate intefaces");
-  }
-}
-
 void TrainingDeepMaterialNetwork::setPenalty(double l)
 {
   _lambda = l;
@@ -1559,7 +1320,7 @@ void TrainingDeepMaterialNetwork::setTrainingSample(int row, const fullMatrix<do
     fullMatrixOperation::reduceMat99ToMat66(C1,_XTrain[row][0]);
     fullMatrixOperation::reduceMat99ToMat66(C2,_XTrain[row][1]);
     fullMatrixOperation::reduceMat99ToMat66(Ceff,_YTrain[row]);
-  };
+  }
 };
 
 void TrainingDeepMaterialNetwork::setTrainingSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& Ceff)
@@ -1598,26 +1359,6 @@ void TrainingDeepMaterialNetwork::setTrainingSample(int row, const fullMatrix<do
   }
 };
 
-void TrainingDeepMaterialNetwork::setTrainingSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2,const fullMatrix<double>& C3, const fullMatrix<double>& C4, const fullMatrix<double>& Ceff)
-{
-  if (_homo->withPlaneStrain())
-  {
-    fullMatrixOperation::reducePlaneStrain(C1,_XTrain[row][0]);
-    fullMatrixOperation::reducePlaneStrain(C2,_XTrain[row][1]);
-    fullMatrixOperation::reducePlaneStrain(C3,_XTrain[row][2]);
-    fullMatrixOperation::reducePlaneStrain(C4,_XTrain[row][3]);
-    fullMatrixOperation::reducePlaneStrain(Ceff,_YTrain[row]);
-  }
-  else
-  {
-    fullMatrixOperation::reduceMat99ToMat66(C1,_XTrain[row][0]);
-    fullMatrixOperation::reduceMat99ToMat66(C2,_XTrain[row][1]);
-    fullMatrixOperation::reduceMat99ToMat66(C3,_XTrain[row][2]);
-    fullMatrixOperation::reduceMat99ToMat66(C4,_XTrain[row][3]);
-    fullMatrixOperation::reduceMat99ToMat66(Ceff,_YTrain[row]);
-  }
-};
-
 void TrainingDeepMaterialNetwork::setTestSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& Ceff)
 {
   static fullMatrix<double> C2;
@@ -1672,26 +1413,6 @@ void TrainingDeepMaterialNetwork::setTestSample(int row, const fullMatrix<double
   }
 };
 
-void TrainingDeepMaterialNetwork::setTestSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& C3, const fullMatrix<double>& C4, const fullMatrix<double>& Ceff)
-{
-  if (_homo->withPlaneStrain())
-  {
-    fullMatrixOperation::reducePlaneStrain(C1,_XTest[row][0]);
-    fullMatrixOperation::reducePlaneStrain(C2,_XTest[row][1]);
-    fullMatrixOperation::reducePlaneStrain(C3,_XTest[row][2]);
-    fullMatrixOperation::reducePlaneStrain(C4,_XTest[row][3]);
-    fullMatrixOperation::reducePlaneStrain(Ceff,_YTest[row]);
-  }
-  else
-  {
-    fullMatrixOperation::reduceMat99ToMat66(C1,_XTest[row][0]);
-    fullMatrixOperation::reduceMat99ToMat66(C2,_XTest[row][1]);
-    fullMatrixOperation::reduceMat99ToMat66(C3,_XTest[row][2]);
-    fullMatrixOperation::reduceMat99ToMat66(C4,_XTest[row][3]);
-    fullMatrixOperation::reduceMat99ToMat66(Ceff,_YTest[row]);
-  }
-};
-
 void TrainingDeepMaterialNetwork::testDataSize(int Nt, int numPhase)
 {
   printf("size of test set = %d numPhase = %d\n",Nt,numPhase);
@@ -1703,9 +1424,7 @@ void TrainingDeepMaterialNetwork::testDataSize(int Nt, int numPhase)
 };
 
 
-void TrainingDeepMaterialNetwork::train(double lr, int maxEpoch, std::string loss, std::string historyFileName, 
-                      std::string saveTreeFileName, int saveEpoch, int numberBatches, double lrmin,
-                      bool removeZeroContribution, int numStepRemove, double tolRemove)
+void TrainingDeepMaterialNetwork::train(double lr, int maxEpoch, std::string loss, std::string historyFileName, std::string saveTreeFileName, int saveEpoch, int numberBatches, double lrmin)
 {
   
   /*
@@ -1729,8 +1448,7 @@ void TrainingDeepMaterialNetwork::train(double lr, int maxEpoch, std::string los
   WBest = Wcur;
   Wprev = Wcur;
   g.resize(numDof,true);
-  double costfuncPrev = evaluateTrainingSet(*lossEval);
-  Msg::Info("costfuncPrev = %e",costfuncPrev);
+  double costfuncPrev = 0.;
   
   //OfflineData
   int epoch = 0;
@@ -1739,7 +1457,6 @@ void TrainingDeepMaterialNetwork::train(double lr, int maxEpoch, std::string los
   int Ntrain = _XTrain.size();
   double lrCur = lr;
   int numEpochIncrease = 0;
-  int numEpochLrMin = 0;
   int numEpochDecrease = 0;
   while (true)
   {
@@ -1808,23 +1525,17 @@ void TrainingDeepMaterialNetwork::train(double lr, int maxEpoch, std::string los
     double costfunc = evaluateTrainingSet(*lossEval);
     if (numberBatches==Ntrain)
     {
-      if (costfunc > costfuncPrev)
+      if ((costfunc > costfuncPrev) && (epoch > 0))
       {
         numEpochIncrease ++;
         if (numEpochIncrease == 1)
         {
           numEpochIncrease = 0;
-          Msg::Info("costfunc = %e, reduce learning rate from %.5e to %.5e",costfunc,lrCur,0.8*lrCur);
+          Msg::Info("reduce learning rate from %.5e to %.5e",lrCur,0.8*lrCur);
           lrCur *= 0.8;
           if (lrCur < lrmin)
           {
             lrCur = lrmin;
-            numEpochLrMin ++;
-            if (numEpochLrMin>10)
-            {
-              Msg::Info("maximal number of iterations %d at minimal learning step is reached !!",numEpochLrMin);
-              break;
-            }
           }
           numEpochDecrease = 0;
           Wcur = Wprev;
@@ -1854,12 +1565,9 @@ void TrainingDeepMaterialNetwork::train(double lr, int maxEpoch, std::string los
     //
     double errTrain = evaluateTrainingSet(*errorEval);
     double errTest = evaluateTestSet(*errorEval);
-    
-    
     Msg::Info("epoch = %d/%d (%.3fs) lrCur=%.5e loss-train =%.5e, loss-test =%.5e mare-train = %.5e mare-test = %.5e",epoch,maxEpoch,Cpu()-time,lrCur,costfunc,lossTest,errTrain,errTest);
     fprintf(f,"%d;%e;%e;%e;%e\n",epoch,costfunc,lossTest,errTrain,errTest);
     fflush(f);
-    
     if ((epoch > 0) && (epoch%saveEpoch==0))
     {
       if (saveTreeFileName.length() > 0)
@@ -1883,30 +1591,6 @@ void TrainingDeepMaterialNetwork::train(double lr, int maxEpoch, std::string los
       }
     };
     
-    if (removeZeroContribution)
-    {
-      
-      if ((epoch+1) %numStepRemove==0)
-      {
-        double costfuncPrev_pre = evaluateTrainingSet(*lossEval);
-        Msg::Info("tree removing");
-        bool removed =  _T->removeLeavesWithZeroContribution(tolRemove);
-        if (removed)
-        {
-          numberDofs(Wcur);
-          Msg::Info("number of Nodes = %d",_T->getNumberOfNodes());
-          Msg::Info("number of unknowns = %d",_unknown.size());
-          int numDof = Wcur.size();
-          WBest = Wcur;
-          Wprev = Wcur;
-          g.resize(numDof,true);
-        }
-        costfuncPrev = evaluateTrainingSet(*lossEval);
-        Msg::Info("pre-removing costfuncPrev = %e and after removing costfuncPrev = %e",costfuncPrev_pre,costfuncPrev);
-      }
-      
-    }
-    
     epoch++;
     if (epoch > maxEpoch)
     {
@@ -1932,22 +1616,9 @@ void TrainingDeepMaterialNetwork::getKeys(const TreeNode* node, std::vector<Dof>
   }
   else
   {
-    if (_sameNormal)
+    for (int i=0; i< node->direction.size(); i++)
     {
-      int numberNormal = node->childs.size()-1;
-      int totalNumberDirVars = node->direction.size();
-      int numberVarsPerNormal = totalNumberDirVars/numberNormal;
-      for (int i=0; i< numberVarsPerNormal; i++)
-      {
-        keys.push_back(Dof(node->location, Dof::createTypeWithTwoInts(i, node->depth+1)));
-      }
-    }
-    else
-    {
-      for (int i=0; i< node->direction.size(); i++)
-      {
-        keys.push_back(Dof(node->location, Dof::createTypeWithTwoInts(i, node->depth+1)));
-      }
+      keys.push_back(Dof(node->location, Dof::createTypeWithTwoInts(i, node->depth+1)));
     }
   }
 };
@@ -2067,25 +1738,9 @@ void TrainingDeepMaterialNetwork::updateFieldFromUnknown(const fullVector<double
     }
     else
     {
-      if (_sameNormal)
-      {
-        int numberNormal = node.childs.size()-1;
-        int totalNumberDirVars = node.direction.size();
-        int numberVarsPerNormal = totalNumberDirVars/numberNormal;
-        for (int j=0; j< numberNormal; j++)
-        {
-          for (int k=0; k< numberVarsPerNormal; k++)
-          {
-            node.direction[j*numberVarsPerNormal+k] = val(k);
-          }
-        }
-      }
-      else
+      for (int j=0; j< node.direction.size(); j++)
       {
-        for (int j=0; j< node.direction.size(); j++)
-        {
-          node.direction[j] = val(j);
-        }        
+        node.direction[j] = val(j);
       }
     };
   };
@@ -2108,6 +1763,58 @@ void TrainingDeepMaterialNetwork::updateFieldFromUnknown(const fullVector<double
   };
 };
 
+void TrainingDeepMaterialNetwork::getNormal(const TreeNode* node, SVector3& vec, bool stiff, std::vector<SVector3>* DnormalDunknown) const
+{
+  double pi = 3.14159265359;
+  int ncomp = node->direction.size(); 
+  if (ncomp == 1)
+  {
+    double theta = 2.*pi*node->direction[0];
+    vec[0] = cos(theta);
+    vec[1] = sin(theta);
+    vec[2] = 0.;
+    if (stiff)
+    {
+      if ((*DnormalDunknown).size() != 1)
+      {
+        (*DnormalDunknown).resize(1);
+      }
+      (*DnormalDunknown)[0][0] = -sin(theta)*2.*pi;
+      (*DnormalDunknown)[0][1] = cos(theta)*2.*pi;
+      (*DnormalDunknown)[0][2] = 0.;
+      
+    }
+  }
+  else if (ncomp == 2)
+  {
+    // cylindrical coordinates
+    double phi = 2.*pi*node->direction[0];
+    double theta = pi*node->direction[1];
+    vec[0] = sin(theta)*cos(phi);
+    vec[1] = sin(theta)*sin(phi);
+    vec[2] = cos(theta);
+    
+    if (stiff)
+    {
+      if ((*DnormalDunknown).size() != 2)
+      {
+        (*DnormalDunknown).resize(2);
+      }
+      
+      (*DnormalDunknown)[0][0] = -sin(theta)*sin(phi)*2*pi;
+      (*DnormalDunknown)[0][1] = sin(theta)*cos(phi)*2*pi;
+      (*DnormalDunknown)[0][2] = 0.;
+      
+      (*DnormalDunknown)[1][0] = cos(theta)*cos(phi)*pi;
+      (*DnormalDunknown)[1][1] = cos(theta)*sin(phi)*pi;
+      (*DnormalDunknown)[1][2] = -sin(theta)*pi;
+    }
+  }
+  else
+  {
+    printf("other dimension %d has not been implemented \n",ncomp+1);
+  }
+};
 
 void TrainingDeepMaterialNetwork::homogenizationProcess(bool stiff)
 {
@@ -2120,9 +1827,9 @@ void TrainingDeepMaterialNetwork::homogenizationProcess(bool stiff)
   static std::vector<double> fin;
   static std::vector<hyperFullMatrix> DCeffDCin;
   static std::vector<fullMatrix<double> > DCeffDfin;
-  static std::vector<std::vector<fullMatrix<double> > > DCeffDnormal;
-  static std::vector<SVector3> normal;
-  static std::vector<std::vector<SVector3> >DnormalDunknown;
+  static std::vector<fullMatrix<double> > DCeffDnormal;
+  static SVector3 normal;
+  static std::vector<SVector3> DnormalDunknown;
   //
   for (int i=0; i< maxDepth; i++)
   {
@@ -2176,10 +1883,9 @@ void TrainingDeepMaterialNetwork::homogenizationProcess(bool stiff)
           fin[j] = Wchild/W;
         };
         // direction parameters
-        node->getNormal(normal,stiff,&DnormalDunknown);
+        getNormal(node,normal,stiff,&DnormalDunknown);
         // homogenized behavior
         _homo->compute(Cin,fin,normal,data.Cel,stiff,DCeffDCin,DCeffDfin,DCeffDnormal);
-        
         if (stiff)
         {
           // update  data.DCeffDunknown
@@ -2198,31 +1904,15 @@ void TrainingDeepMaterialNetwork::homogenizationProcess(bool stiff)
           // keys inside node - rotation dofs
           std::vector<Dof> Rrotation;
           getKeys(node,Rrotation);
-          
-          int numberNormal = node->childs.size()-1;
-          int totalNumberDirVars = node->direction.size();
-          int numberVarsPerNormal = totalNumberDirVars/numberNormal;
-          
-          for (int j=0; j< numberNormal; j++)
+          // the first Dofs are rotation
+          for (int j=0; j< Rrotation.size(); j++)
           {
-            for (int k=0; k< numberVarsPerNormal; k++)
-            {  
-              if (_sameNormal)
-              {
-                fullMatrix<double>& DCeffDnormalDof = data.getDCeffDunknown(Rrotation[k],_homo->withPlaneStrain());
-                DCeffDnormalDof.axpy(DCeffDnormal[j][0],DnormalDunknown[j][k][0]); // comp 0
-                DCeffDnormalDof.axpy(DCeffDnormal[j][1],DnormalDunknown[j][k][1]); // comp 1
-                DCeffDnormalDof.axpy(DCeffDnormal[j][2],DnormalDunknown[j][k][2]); // comp 2
-              }
-              else
-              {
-                fullMatrix<double>& DCeffDnormalDof = data.getDCeffDunknown(Rrotation[j*numberVarsPerNormal+k],_homo->withPlaneStrain());
-                DCeffDnormalDof.axpy(DCeffDnormal[j][0],DnormalDunknown[j][k][0]); // comp 0
-                DCeffDnormalDof.axpy(DCeffDnormal[j][1],DnormalDunknown[j][k][1]); // comp 1
-                DCeffDnormalDof.axpy(DCeffDnormal[j][2],DnormalDunknown[j][k][2]); // comp 2
-              }
-            }
-          }          
+            fullMatrix<double>& DCeffDnormalDof = data.getDCeffDunknown(Rrotation[j],_homo->withPlaneStrain());
+            DCeffDnormalDof.axpy(DCeffDnormal[0],DnormalDunknown[j][0]); // comp 0
+            DCeffDnormalDof.axpy(DCeffDnormal[1],DnormalDunknown[j][1]); // comp 1
+            DCeffDnormalDof.axpy(DCeffDnormal[2],DnormalDunknown[j][2]); // comp 2
+          };
+          
           // for phase percentage
           std::vector<double> DWDWeight(leaves.size(),0.);
           for (int j=0; j< leaves.size(); j++)
@@ -2543,12 +2233,6 @@ DeepMaterialNetwork::DeepMaterialNetwork(const Tree& T): _T(&T), _isInitialized(
   #if defined(HAVE_PETSC) || defined(HAVE_GMM)
   _sys = NULL;
   #endif //
-  
-  if (_T->checkDuplicate())
-  {
-    Msg::Error("Duplicating node exists, the computation cannot continue !!!");
-    Msg::Exit(0);
-  }
 };
 
 
@@ -2567,6 +2251,33 @@ DeepMaterialNetwork::~DeepMaterialNetwork()
   }
 };
 
+SVector3 DeepMaterialNetwork::getNormal(const TreeNode* node) const
+{
+  double pi = 3.14159265359;
+  int ncomp = node->direction.size(); 
+  SVector3 vec;
+  if (ncomp == 1)
+  {
+    double theta = 2.*pi*node->direction[0];
+    vec[0] = cos(theta);
+    vec[1] = sin(theta);
+    vec[2] = 0.;
+  }
+  else if (ncomp == 2)
+  {
+    // cylindrical coordinates
+    double phi = 2.*pi*node->direction[0];
+    double theta = pi*node->direction[1];
+    vec[0] = sin(theta)*cos(phi);
+    vec[1] = sin(theta)*sin(phi);
+    vec[2] = cos(theta);
+  }
+  else
+  {
+    printf("other dimension %d has not been implemented \n",ncomp+1);
+  }
+  return vec;
+};
 
 double DeepMaterialNetwork::get(int comp)
 {
@@ -2614,7 +2325,7 @@ void DeepMaterialNetwork::initializeOnlineEvaluation()
     _T->getBackTrackingPath(_leaves[i],path);
     _backTrackingMap[_leaves[i]] = path;
   };
-  printf("size of void leaves = %ld\n",_voidLeaves.size());
+  printf("size of void leaves = %d\n",_voidLeaves.size());
   
   // leaf map
   printf("creating associated leaves map\n");
@@ -2705,7 +2416,7 @@ void DeepMaterialNetwork::initializeOnlineEvaluation()
     }
   };
   
-  printf("size of void _notUnknownRegularNodes = %ld\n",_notUnknownRegularNodes.size());
+  printf("size of void _notUnknownRegularNodes = %d\n",_notUnknownRegularNodes.size());
 };
 
 void DeepMaterialNetwork::getKeys(const TreeNode* node, std::vector<Dof>& keys) const
@@ -2781,8 +2492,7 @@ void DeepMaterialNetwork::getLinearTerm(const TreeNode* node, fullVector<double>
   {
     vec.setAll(0.);
   }
-  std::vector<SVector3> normal;
-  node->getNormal(normal);
+  SVector3 normal = getNormal(node);  
   for (int k=0; k< numChilds-1; k++)
   {
     const TreeNode* childLeft = node->childs[k];
@@ -2794,7 +2504,7 @@ void DeepMaterialNetwork::getLinearTerm(const TreeNode* node, fullVector<double>
     {
       for (int j=0; j<3; j++)
       {
-        vec(i+3*k) += (PL(i,j)-PR(i,j))*normal[k](j);
+        vec(i+3*k) += (PL(i,j)-PR(i,j))*normal(j);
       }
     }    
   }
@@ -2881,8 +2591,7 @@ void DeepMaterialNetwork::getTangentTerm(const TreeNode* node, fullMatrix<double
     mat.setAll(0.);
   }
   //
-  std::vector<SVector3> normal;
-  node->getNormal(normal);
+  SVector3 normal = getNormal(node);
   for (int k=0; k< numChilds-1; k++)
   {
     const TreeNode* childLeft = node->childs[k];
@@ -2910,7 +2619,7 @@ void DeepMaterialNetwork::getTangentTerm(const TreeNode* node, fullMatrix<double
             {
               for (int t=0; t<3; t++)
               {
-                mat(p+3*k,Tensor23::getIndex(q,r)) += L(p,t,q,r)*normal[k](t)*getWeight(no)/wL;
+                mat(p+3*k,Tensor23::getIndex(q,r)) += L(p,t,q,r)*normal(t)*getWeight(no)/wL;
               }
             }
           }
@@ -2935,7 +2644,7 @@ void DeepMaterialNetwork::getTangentTerm(const TreeNode* node, fullMatrix<double
             {
               for (int t=0; t<3; t++)
               {
-                mat(p+3*k,Tensor23::getIndex(q,r)) -= L(p,t,q,r)*normal[k](t)*getWeight(no)/wR;
+                mat(p+3*k,Tensor23::getIndex(q,r)) -= L(p,t,q,r)*normal(t)*getWeight(no)/wR;
               }
             }
           }
@@ -2973,8 +2682,7 @@ void DeepMaterialNetwork::getBilinearTerm(const TreeNode* node, Tree::nodeContai
   }
   
   //
-  std::vector<SVector3> normal;
-  node->getNormal(normal);
+  SVector3 normal = getNormal(node);
   for (int k=0; k< numChilds-1; k++)
   {
     const TreeNode* childLeft = node->childs[k];
@@ -3006,7 +2714,7 @@ void DeepMaterialNetwork::getBilinearTerm(const TreeNode* node, Tree::nodeContai
             {
               for (int t=0; t<3; t++)
               {
-                LDFlocal(p,q,r) += L(p,t,q,r)*normal[k](t)*getWeight(no)/wL;
+                LDFlocal(p,q,r) += L(p,t,q,r)*normal(t)*getWeight(no)/wL;
               }
             }
           }
@@ -3059,7 +2767,7 @@ void DeepMaterialNetwork::getBilinearTerm(const TreeNode* node, Tree::nodeContai
             {
               for (int t=0; t<3; t++)
               {
-                LDFlocal(p,q,r) -= L(p,t,q,r)*normal[k](t)*getWeight(no)/wR;
+                LDFlocal(p,q,r) -= L(p,t,q,r)*normal(t)*getWeight(no)/wR;
               }
             }
           }
@@ -3199,8 +2907,7 @@ void DeepMaterialNetwork::downscale(const STensor3& Fcur, STensor3& F,  const Tr
   {
     const TreeNode* nodeParent = path[i+1]; 
     const OnlineNodeData* dataNodeParent = getOnlineNodeData(nodeParent);
-    std::vector<SVector3> allNormal;
-    nodeParent->getNormal(allNormal);
+    SVector3 normal = getNormal(nodeParent);
     double wParent = getWeight(nodeParent);
     int numChilds = nodeParent->childs.size();
     const std::vector<SVector3>& acur = dataNodeParent->acur;
@@ -3210,10 +2917,9 @@ void DeepMaterialNetwork::downscale(const STensor3& Fcur, STensor3& F,  const Tr
     double wChild = getWeight(node);
     int childOrder = node->childOrder;
     //
-   //double f = wChild/wParent;
-    double fact = 1./wChild;
+    //double f = wChild/wParent;
+    //double fact = 1./f;
     
-    /*
     double fact = 1.;
     for (int ic=0; ic< numChilds; ic++)
     {
@@ -3221,81 +2927,54 @@ void DeepMaterialNetwork::downscale(const STensor3& Fcur, STensor3& F,  const Tr
       {
         fact *= (getWeight(nodeParent->childs[ic])/wParent);
       }
-    }*/
-    
-    std::vector<STensor33>& DFDa = DFcurDa[i];
+    }
     
     static STensor3 I(1.);
+    SVector3 aaStrain;
+    std::vector<STensor3> DaaStrainDa(numChilds-1,STensor3(0.));
     if (childOrder == 0)
     {
       // first phase
-      for (int r=0; r<3; r++)
-      {
-        for (int c=0; c<3; c++)
-        {
-          F(r,c) += acur[0](r)*allNormal[0](c)*fact;
-        }
-      }
-      
-      for (int r=0; r<3; r++)
-      {
-        for (int c=0; c<3; c++)
-        {
-          for (int t=0; t< 3; t++)
-          {
-            DFDa[0](r,c,t) = I(r,t)*allNormal[0](c)*fact;
-          }
-        }
-      };
-      
+      aaStrain = acur[0];
+      STensorOperation::diag(DaaStrainDa[0],1.);
     }
     else if (childOrder == numChilds-1)
     {
-      // last phase      
-      for (int r=0; r<3; r++)
-      {
-        for (int c=0; c<3; c++)
-        {
-          F(r,c) += -acur[numChilds-2](r)*allNormal[numChilds-2](c)*fact;
-        }
-      }
-      
-      for (int r=0; r<3; r++)
-      {
-        for (int c=0; c<3; c++)
-        {
-          for (int t=0; t< 3; t++)
-          {
-            DFDa[numChilds-2](r,c,t) = -I(r,t)*allNormal[numChilds-2](c)*fact;
-          }
-        }
-      };
-      
+      // last phase
+      aaStrain = acur[numChilds-2];
+      aaStrain *= (-1.);
+      STensorOperation::diag(DaaStrainDa[numChilds-2],-1.);
     }
     else
     {
-      //
-      for (int r=0; r<3; r++)
+      aaStrain = acur[childOrder] - acur[childOrder-1];
+      STensorOperation::diag(DaaStrainDa[childOrder],1.);
+      STensorOperation::diag(DaaStrainDa[childOrder-1],-1.);
+    }
+    
+    for (int r=0; r<3; r++)
+    {
+      for (int c=0; c<3; c++)
       {
-        for (int c=0; c<3; c++)
-        {
-          F(r,c) += (acur[childOrder](r)*allNormal[childOrder](c) -  acur[childOrder-1](r)*allNormal[childOrder-1](c))*fact;
-        }
+        F(r,c) += aaStrain(r)*normal(c)*fact;
       }
-      
+    }
+    
+    std::vector<STensor33>& DFDa = DFcurDa[i];
+    for (int k=0; k< numChilds-1; k++)
+    {
       for (int r=0; r<3; r++)
       {
         for (int c=0; c<3; c++)
         {
           for (int t=0; t< 3; t++)
           {
-            DFDa[childOrder](r,c,t) = I(r,t)*allNormal[childOrder](c)*fact;
-            DFDa[childOrder-1](r,c,t) = -I(r,t)*allNormal[childOrder-1](c)*fact;
+            DFDa[k](r,c,t) = DaaStrainDa[k](r,t)*normal(c)*fact;
           }
         }
       };
     }
-  }
+  };
 };
 
 void DeepMaterialNetwork::setVoidMatIndex(int index)
@@ -3620,7 +3299,7 @@ bool DeepMaterialNetwork::stress(const STensor3& F, STensor3& P, bool stiff, STe
       {
         assembleDresDu();
       }
-      bool succ = tangentSolve(L);
+      tangentSolve(L);
     }
     //
     return true;
@@ -3787,13 +3466,8 @@ bool DeepMaterialNetwork::tangentSolve(STensor43& L)
       }
     }
   }
-  bool ok = _sys->systemSolveMultipleRHS();
-  if (!ok)
-  {
-    Msg::Error("system multiple RHS is not sucessfully solved");
-    return false;
-  }
-
+  _sys->systemSolveMultipleRHS();
+  
   STensorOperation::zero(L);
   const TreeNode* root = _T->getRootNode();
   double Wroot = getWeight(root);
@@ -3871,7 +3545,6 @@ bool DeepMaterialNetwork::tangentSolve(STensor43& L)
     }
   };
   //L.print("homogenized tangent");
-  return true;
 }
 
 void DeepMaterialNetwork::systemNextStep()
@@ -3960,7 +3633,6 @@ bool DeepMaterialNetwork::systemSolve()
 bool DeepMaterialNetwork::tangentSolve(STensor43& L)
 {
   Msg::Error("tangent solver requires using PETSC ");
-  return false;
 }
 
 void DeepMaterialNetwork::systemNextStep()
diff --git a/NonLinearSolver/modelReduction/DeepMaterialNetworks.h b/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
index 51ce1bf3d9116b9ca7f414b932f11ad09c04786c..c1ccc4c0a85a42ac9441e24292dac84891782e02 100644
--- a/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
+++ b/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
@@ -17,7 +17,7 @@
 #include "highOrderTensor.h"
 #include "STensorOperations.h"
 #include "dofManager.h"
-#include "ipField.h"
+#include "staticDofManager.h"
 #if defined(HAVE_PETSC)
 template<class T>
 class nonLinearSystemPETSc;
@@ -90,14 +90,12 @@ class Homogenization
   public:
     #ifndef SWIG
     virtual ~Homogenization(){}
-    virtual void compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, 
-                        const std::vector<SVector3>& normal,
+    virtual void compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, const SVector3& normal,
                         fullMatrix<double>& Ceff, 
                         bool stiff, 
                         std::vector<hyperFullMatrix>& DCeffDcin, 
-                        std::vector<fullMatrix<double> >& DCeffDf, 
-                        std::vector<std::vector<fullMatrix<double> > >&DCeffDnormal) const = 0; 
-                 
+                        std::vector<fullMatrix<double>>& DCeffDf, 
+                        std::vector<fullMatrix<double>>&DCeffDnormal) const = 0; 
     virtual void printInfos() const = 0;
     virtual bool withPlaneStrain() const  = 0;
     #endif //SWIG
@@ -115,15 +113,6 @@ class BimaterialHomogenization : public Homogenization
                         fullMatrix<double>& DCeffDf1, fullMatrix<double>& DCeffDf2 ) const;
     #ifndef SWIG
     virtual ~BimaterialHomogenization(){}
-    
-    virtual void compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, 
-                        const std::vector<SVector3>& normal,
-                        fullMatrix<double>& Ceff, 
-                        bool stiff, 
-                        std::vector<hyperFullMatrix>& DCeffDcin, 
-                        std::vector<fullMatrix<double> >& DCeffDf, 
-                        std::vector<std::vector<fullMatrix<double> > >&DCeffDnormal) const; 
-    
     virtual void compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, const SVector3& normal,
                         fullMatrix<double>& Ceff, 
                         bool stiff, 
@@ -142,10 +131,10 @@ class BimaterialHomogenization : public Homogenization
           fullMatrix<double>& DNDn2) const;
       virtual void computeBimaterial(const fullMatrix<double>& C1, const fullMatrix<double>& C2, double f1, double f2, const SVector3& normal,
                         fullMatrix<double>& Ceff, 
-                        bool stiff=false, 
-                        hyperFullMatrix* DCeffDC1=NULL, hyperFullMatrix* DCeffDC2=NULL,
-                        fullMatrix<double>* DCeffDf1=NULL, fullMatrix<double>* DCeffDf2=NULL, 
-                        fullMatrix<double>* DCeffDnormal1=NULL,fullMatrix<double>* DCeffDnormal2=NULL, fullMatrix<double>* DCeffDnormal3=NULL) const;
+                        bool stiff, 
+                        hyperFullMatrix* DCeffDC1, hyperFullMatrix* DCeffDC2,
+                        fullMatrix<double>* DCeffDf1, fullMatrix<double>* DCeffDf2, 
+                        fullMatrix<double>* DCeffDnormal1,fullMatrix<double>* DCeffDnormal2, fullMatrix<double>* DCeffDnormal3) const;
     #endif //SWIG
 };
 
@@ -166,14 +155,6 @@ class LaminateHomogenization : public BimaterialHomogenization
                         std::vector<hyperFullMatrix>& DCeffDcin, 
                         std::vector<fullMatrix<double>>& DCeffDf, 
                         std::vector<fullMatrix<double>>&DCeffDnormal) const;
-                        
-    virtual void compute(const std::vector<fullMatrix<double> >& Cin, const std::vector<double>& f, 
-                        const std::vector<SVector3>& normal,
-                        fullMatrix<double>& Ceff, 
-                        bool stiff, 
-                        std::vector<hyperFullMatrix>& DCeffDcin, 
-                        std::vector<fullMatrix<double> >& DCeffDf, 
-                        std::vector<std::vector<fullMatrix<double> > >&DCeffDnormal) const; 
     
     virtual void printInfos() const;
     #endif //SWIG
@@ -265,34 +246,27 @@ class TrainingDeepMaterialNetwork
     std::vector<fullMatrix<double> > _YTrain;
     std::vector<std::vector<fullMatrix<double> > >_XTest;
     std::vector<fullMatrix<double> > _YTest;
-    bool _sameNormal;
     
   #endif //SWIG  
   public:
     TrainingDeepMaterialNetwork(Tree& T, const Homogenization& ho);
     ~TrainingDeepMaterialNetwork();
     void setPenalty(double l);
-    void sameNormal(bool fl);
     void trainingDataSize(int Ns, int numPhase);
     void setTrainingSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& Ceff);
     void setTrainingSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& Ceff);
     void setTrainingSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& C3, const fullMatrix<double>& Ceff);
-    void setTrainingSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& C3, const fullMatrix<double>& C4, const fullMatrix<double>& Ceff);
     void testDataSize(int Ns, int numPhase);
     void setTestSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& Ceff);
     void setTestSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& Ceff);
     void setTestSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& C3, const fullMatrix<double>& Ceff);
-    void setTestSample(int row, const fullMatrix<double>& C1, const fullMatrix<double>& C2, const fullMatrix<double>& C3, const fullMatrix<double>& C4, const fullMatrix<double>& Ceff);
     void train(double lr, int maxEpoch, 
               std::string loss = "mare", 
               std::string historyFileName="history.csv", 
               std::string saveTreeFileName="tree_trained.txt", 
               int saveEpoch=50, 
               int numberBatches=1,
-              double lrmin=1e-4,
-              bool removeZeroContribution = false,
-              int numStepRemove = 3,
-              double tolRemove=1e-4);
+              double lrmin=1e-4);
     // test homogenization procedure
     void testHomogenizationProcess(fullMatrix<double>& C1, const fullMatrix<double>& C2, fullMatrix<double>& Ceff);
     
@@ -304,6 +278,7 @@ class TrainingDeepMaterialNetwork
     int getDofNumber(Dof key) const;
     void updateFieldFromUnknown(const fullVector<double>& ufield);
     
+    void getNormal(const TreeNode* node, SVector3& vec, bool stiff, std::vector<SVector3>* DnormalDunknown) const;
     void homogenizationProcess(bool stiff);
     double computeError(const LossMeasure& loss, const std::vector<fullMatrix<double> >& Cin, const fullMatrix<double>& Ceff, bool stiff=false, fullVector<double>* g=NULL);
                  
@@ -367,6 +342,7 @@ class DeepMaterialNetwork
     //
     
   protected:
+    SVector3 getNormal(const TreeNode* node) const;
     const OnlineNodeData* getOnlineNodeData(const TreeNode* node) const;
     const Tree::nodeContainer* getAssociatedLeaves(const TreeNode* node) const;
     const Tree::nodeContainer* getBackTrackingPath(const TreeNode* leaf) const;
diff --git a/NonLinearSolver/modelReduction/Tree.cpp b/NonLinearSolver/modelReduction/Tree.cpp
index e123656bae15a5bfdfc577d41054c163c151a96d..6e7b591c20c00ccfc3846eb9627cd74585759d6e 100644
--- a/NonLinearSolver/modelReduction/Tree.cpp
+++ b/NonLinearSolver/modelReduction/Tree.cpp
@@ -21,8 +21,8 @@ direction(), weight(0.), childOrder(-1), af(NULL), matIndex(-1)
   
 }
 
-TreeNode::TreeNode(int k, int c, TreeNode* p, int o, const std::string affunc): depth(k), location(c), childs(0),parent(p),
-direction(0), weight(0.), childOrder(o), matIndex(-1)
+TreeNode::TreeNode(int k, int c, TreeNode* p, int numDirVars, int o, const std::string affunc): depth(k), location(c), childs(0),parent(p),
+direction(numDirVars,0), weight(0.), childOrder(o), matIndex(-1)
 {
   af = activationFunction::createActivationFunction(affunc.c_str());
 };
@@ -33,7 +33,7 @@ TreeNode::~TreeNode()
 
 void TreeNode::printData() const
 {
-  printf("Node: depth %d- loc %d  -mat index %d- child order %d - weight = %e - nb of direction vars %ld",depth,location,matIndex,childOrder,weight,direction.size());
+  printf("Node: depth %d- loc %d  -mat index %d- child order %d - weight = %e - direction",depth,location,matIndex,childOrder,weight);
   for (int i=0; i< direction.size(); i++)
   {
     printf(" %e",direction[i]);
@@ -49,7 +49,7 @@ void TreeNode::saveToFile(FILE* f) const
   {
     parentloc = parent->location;
   }
-  fprintf(f,"%d %d %d %d %d %ld %ld",depth,location, matIndex, childOrder,parentloc,childs.size(), direction.size());
+  fprintf(f,"%d %d %d %d %d %d %d",depth,location, matIndex, childOrder,parentloc,childs.size(), direction.size());
   //
   if (childs.size() > 0)
   {
@@ -64,89 +64,6 @@ void TreeNode::saveToFile(FILE* f) const
     fprintf(f," %f",direction[i]);
   }
   fprintf(f," %s\n",af->getName().c_str());
-  fflush(f);
-};
-
-void TreeNode::getNormal(std::vector<SVector3>& allVec, bool stiff, std::vector<std::vector<SVector3> >* allDnormalDunknown) const
-{
-  // number of normal
-  double pi = 3.14159265359;
-  int numberNormal = childs.size()-1;
-  int totalNumberDirVars = direction.size();
-  int numberVarsPerNormal = totalNumberDirVars/numberNormal;
-  // reallocated
-  if (numberNormal > 0)
-  {
-    if (allVec.size() != numberNormal)
-    {
-      allVec.resize(numberNormal);
-    }
-    
-    if (stiff)
-    {
-      if (allDnormalDunknown->size() != numberNormal)
-      {
-        allDnormalDunknown->resize(numberNormal,std::vector<SVector3>(numberVarsPerNormal,SVector3()));
-      }
-    }
-  }
-  
-  for (int i=0; i< numberNormal; i++)
-  {
-    SVector3& vec = allVec[i];
-    std::vector<SVector3>* DnormalDunknown = NULL;
-    if (stiff)
-    {
-      DnormalDunknown = &(*allDnormalDunknown)[i];
-    }
-    if (numberVarsPerNormal == 1)
-    {
-      double theta = 2.*pi*direction[i];
-      vec[0] = cos(theta);
-      vec[1] = sin(theta);
-      vec[2] = 0.;
-      if (stiff)
-      {
-        if ((*DnormalDunknown).size() != 1)
-        {
-          (*DnormalDunknown).resize(1);
-        }
-        (*DnormalDunknown)[0][0] = -sin(theta)*2.*pi;
-        (*DnormalDunknown)[0][1] = cos(theta)*2.*pi;
-        (*DnormalDunknown)[0][2] = 0.;
-        
-      }
-    }
-    else if (numberVarsPerNormal == 2)
-    {
-      // cylindrical coordinates
-      double phi = 2.*pi*direction[i*numberVarsPerNormal];
-      double theta = pi*direction[i*numberVarsPerNormal+1];
-      vec[0] = sin(theta)*cos(phi);
-      vec[1] = sin(theta)*sin(phi);
-      vec[2] = cos(theta);
-      if (stiff)
-      {
-        if ((*DnormalDunknown).size() != 2)
-        {
-          (*DnormalDunknown).resize(2);
-        }
-        
-        (*DnormalDunknown)[0][0] = -sin(theta)*sin(phi)*2*pi;
-        (*DnormalDunknown)[0][1] = sin(theta)*cos(phi)*2*pi;
-        (*DnormalDunknown)[0][2] = 0.;
-        
-        (*DnormalDunknown)[1][0] = cos(theta)*cos(phi)*pi;
-        (*DnormalDunknown)[1][1] = cos(theta)*sin(phi)*pi;
-        (*DnormalDunknown)[1][2] = -sin(theta)*pi;
-      }
-    }
-    else
-    {
-      Msg::Error("number of vars per normal must be 1 or 2, a value of %d is impossible",numberVarsPerNormal);
-      Msg::Exit(0);
-    }
-  }
 };
 
 
@@ -163,7 +80,8 @@ Tree::Tree(): _root(NULL)
 void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::string affunc)
 {
   clear();
-  _root = new TreeNode(0,0,NULL,0);
+  _root = new TreeNode(0,0,NULL,numDirVars,0);
+  _nodeMap[std::pair<int,int>(0,0)]=_root;
   
   TreeNode* np = _root;
   int depth = 1;
@@ -174,13 +92,13 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri
     numAllocatedNodes += 2;
     //
     np->childs.resize(2);
-    np->direction.resize(numDirVars);
-    np->childs[0] = new TreeNode(depth,0,np,0,affunc);
+    np->childs[0] = new TreeNode(depth,0,np,numDirVars,0,affunc);
+    _nodeMap[std::pair<int,int>(depth,0)]=np->childs[0];
     np->childs[0]->matIndex = matStart;
     
     if (numAllocatedNodes >= numNodes)
     {
-      np->childs[1] = new TreeNode(depth,1,np,1,affunc);
+      np->childs[1] = new TreeNode(depth,1,np,numDirVars,1,affunc);
       if (np->childs[0]->matIndex == 0)
       {
         np->childs[1]->matIndex = 1;
@@ -192,8 +110,9 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri
     }
     else
     {
-      np->childs[1] = new TreeNode(depth,1,np,1);
+      np->childs[1] = new TreeNode(depth,1,np,numDirVars,1);
     }
+    _nodeMap[std::pair<int,int>(depth,1)]=np->childs[1];
     if (matStart == 0) matStart = 1;
     else matStart= 0;
     depth++;
@@ -204,7 +123,9 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri
 void Tree::createSpecialTripleTree(int numNodes, int numDirVars, const std::string affunc)
 {
   clear();
-  _root = new TreeNode(0,0,NULL,0);  
+  _root = new TreeNode(0,0,NULL,numDirVars,0);
+  _nodeMap[std::pair<int,int>(0,0)]=_root;
+  
   TreeNode* np = _root;
   int depth = 1;
   int numAllocatedNodes = 1;
@@ -215,21 +136,22 @@ void Tree::createSpecialTripleTree(int numNodes, int numDirVars, const std::stri
     if (numAllocatedNodes < numNodes)
     {
       np->childs.resize(3);
-      np->direction.resize(2*numDirVars);
     }
     else
     {
       np->childs.resize(2);
-      np->direction.resize(numDirVars);
     }
-    np->childs[0] = new TreeNode(depth,0,np,0,affunc);
+    np->childs[0] = new TreeNode(depth,0,np,numDirVars,0,affunc);
+    _nodeMap[std::pair<int,int>(depth,0)]=np->childs[0];
     np->childs[0]->matIndex = 0;
-    np->childs[1] = new TreeNode(depth,1,np,1,affunc);
+    np->childs[1] = new TreeNode(depth,1,np,numDirVars,1,affunc);
+    _nodeMap[std::pair<int,int>(depth,1)]=np->childs[1];
     np->childs[1]->matIndex = 1;
     
     if (numAllocatedNodes < numNodes)
     {
-      np->childs[2] = new TreeNode(depth,2,np,2);
+      np->childs[2] = new TreeNode(depth,2,np,numDirVars,2);
+      _nodeMap[std::pair<int,int>(depth,2)]=np->childs[2];
       np = np->childs[2];
     }
     
@@ -237,419 +159,6 @@ void Tree::createSpecialTripleTree(int numNodes, int numDirVars, const std::stri
   }
 };
 
-void Tree::createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int leafPerNode, int numDirVars,const std::string affunc, bool randomChild, bool randomMat)
-{ 
-  auto makeGroup = [&](const std::vector<int>& input, std::vector<int>& outPut, bool first)
-  {
-    outPut.clear();
-    int numInPut = input.size();
-    int iter = 0;
-    while (true)
-    {
-      if (numInPut == 2 || numInPut == 3)
-      {
-        outPut.push_back(numInPut);
-        break;
-      }
-      else if (numInPut == 4)
-      {
-        outPut.push_back(2);
-        outPut.push_back(2);
-        break;
-      }
-            
-      int numChilds;
-      if (first)
-      {
-        numChilds = leafPerNode;
-        if (randomMat)
-        {
-          numChilds =  rand() % (leafPerNode-1) + 2; // number takes a random value from 2 to leafPerNode 
-        }
-      }
-      else
-      {
-        numChilds = numChildMax;
-        if (randomChild)
-        {
-          numChilds =  rand() % (numChildMax-1) + 2; // number takes a random value from 2 to numChildMax 
-        }
-      }
-      if (numInPut <= numChilds)
-      {
-        outPut.push_back(numInPut);
-        break;
-      }
-      int remaining = numInPut - iter - numChilds;
-      if (remaining > 1)
-      {
-        outPut.push_back(numChilds);
-        iter += numChilds;
-      }
-      else 
-      {
-        if (first)
-        {
-          if (numChilds+remaining > leafPerNode)
-          {
-            int num1 = leafPerNode;
-            int num2 =numChilds+remaining -leafPerNode;
-            while (num2 < 2)
-            {
-              num1 --;
-              num2 ++;
-            }
-            outPut.push_back(num1);
-            outPut.push_back(num2);
-            break;
-          }
-          else
-          {
-            outPut.push_back(numChilds+remaining);
-            break;
-          }
-        }
-        else
-        {
-          if (numChilds+remaining < 4)
-          {
-            outPut.push_back(numChilds+remaining);
-          }
-          else
-          {
-            outPut.push_back(2);
-            outPut.push_back(numChilds+remaining-2);
-          }
-          iter += numChilds+remaining;
-        }
-      };
-      if (iter >= numInPut) break;
-    }
-  };
-  
-  auto printVector = [](const std::vector<int>& vv)
-  {
-    printf("%ld:",vv.size());
-    for (int i=0; i< vv.size(); i++)
-    {
-      printf("%d ",vv[i]);
-    }
-    printf("\n");
-  }; 
-  
-  
-  std::vector<std::vector<int> > allListByDepth;
-  
-  // first list
-  allListByDepth.resize(1);
-  for (int i=0; i< nbLeaves; i++)
-  {
-    allListByDepth[0].push_back(1);
-  }
-  std::vector<int> intPut = allListByDepth[0];
-  printVector(intPut);
-  int iter=0;
-  while (true)
-  {
-    std::vector<int> outPut;
-    makeGroup(intPut,outPut,iter==0);
-    allListByDepth.push_back(outPut);
-    printVector(outPut);
-    intPut = outPut;
-    if (outPut.size() == 1)
-    {
-      break;
-    }
-    iter++;
-  };
-  
-  clear();
-  if (nbLeaves <0) return;
-  srand (time(NULL));
-  
-  int depthMax = allListByDepth.size()-1;
-  
-  std::vector<TreeNode*> listPrev, listCur;  
-  _root = new TreeNode(0,0,NULL,0);
-  
-  listPrev.push_back(_root);
-  int depth = 1;
-  while (true)
-  {
-    listCur.clear();
-    int inode=0;
-    
-    for (int i =0; i< listPrev.size(); i++)
-    {
-      TreeNode* np = listPrev[i];          
-      //
-      printVector(allListByDepth[depthMax-depth+1]);
-      int trueNbChilds = allListByDepth[depthMax-depth+1][i];
-      Msg::Info("i=%d/%d,trueNbChilds = %d",i,listPrev.size(),trueNbChilds);
-      np->childs.resize(trueNbChilds);
-      np->direction.resize((trueNbChilds-1)*numDirVars);
-      for (int j=0; j< trueNbChilds; j++)
-      {
-        if (depth == depthMax)
-        {
-          // leaf
-          np->childs[j] = new TreeNode(depth,inode,np,j,affunc);
-        }
-        else
-        {
-          np->childs[j] = new TreeNode(depth,inode,np,j);
-          listCur.push_back(np->childs[j]);
-        }
-        inode++;
-      }
-    }
-    depth ++;
-    listPrev = listCur;
-    if (listCur.size() ==0)
-    {
-      break;
-    }
-  };
-  std::vector<TreeNode*> allLeaves;
-  getRefToAllLeaves(allLeaves);
-  std::set<TreeNode*> setNodes;
-  for (int i=0; i< allLeaves.size(); i++)
-  {
-    TreeNode* leaf = allLeaves[i];
-    if (setNodes.find(leaf) == setNodes.end())
-    {
-      TreeNode* parent = leaf->parent;
-      std::vector<TreeNode*> ll;
-      for (int j=0; j< parent->childs.size(); j++)
-      {
-        TreeNode* nj = parent->childs[j];
-        if (nj->childs.size() == 0)
-        {
-          ll.push_back(nj);
-          setNodes.insert(nj);
-        }
-      }
-      if (ll.size() == leafPerNode)
-      {
-        for (int j=0; j< ll.size(); j++)
-        {
-          ll[j]->matIndex  = j;
-        }
-      }
-      else
-      {
-        std::vector<int> v(leafPerNode,0);
-        for (int j=0; j< leafPerNode; j++)
-        {
-          v[j] = j;
-        }
-        for (int j=0; j< ll.size(); j++)
-        {
-          int pos =  rand() % (v.size()); // random from 0 to v.size() -1
-          TreeNode* otherNode = ll[j];
-          otherNode->matIndex = v[pos];
-          v.erase(v.begin() + pos);
-        }
-      }
-    }
-  }
-  
-    
-};
-
-void Tree::createRandomTree(int nbLeaves, int numChildMax, int leafPerNode, int numDirVars,const std::string affunc, bool randomSubdivion, bool randomChild, bool sameNbChildLayer)
-{
-  
-  auto randomDistribution = [&randomSubdivion](int inSize, int numPart, int leafPerNode, std::vector<int>& out)
-  {
-    out.clear();
-    bool suc = true;
-    if (inSize <= numPart)
-    {
-      if (inSize > leafPerNode)
-      {
-        out.resize(2);
-        out[0] = leafPerNode;
-        out[1] = inSize - leafPerNode;
-      }
-      else
-      {
-        out.resize(inSize);
-        for (int i=0; i< inSize; i++)
-        {
-          out[i] = 1;
-        }
-      }
-    }
-    else
-    {
-      
-      int maxIter = 0;
-      while (maxIter <10)
-      {
-        out.resize(numPart,0);
-        int remaining = inSize;
-        int remainingPart = numPart;
-        int numOneSize = 0;
-        for (int i=0; i< numPart-1; i++)
-        {
-          if (randomSubdivion)
-          {
-            out[i] = rand()%(remaining-remainingPart)+1;
-          }
-          else
-          {
-            out[i] = inSize/numPart;
-          }
-          remainingPart--;
-          remaining -= out[i];    
-          if (out[i] == 1) numOneSize++;
-        };
-        out[numPart-1] = remaining;
-        if (remaining == 1) numOneSize++;
-        if (remaining == 0 || numOneSize > leafPerNode)
-        {
-          Msg::Error("re distribution remaining = %d  numOneSize = %d",remaining,numOneSize);
-          numPart --;
-          suc = false;
-        }
-        else
-        {
-          suc = true;
-          break;
-        }
-        maxIter++;
-      }
-    }
-    
-    Msg::Info("partition %d into %d parts: ",inSize,out.size());
-    for (int i=0; i< out.size(); i++)
-    {
-      Msg::Info("num part %d = %d",i,out[i]);
-    }
-    if (!suc)
-    {
-      Msg::Exit(0);
-    }
-    
-  };
-  
-  clear();
-  if (nbLeaves <0) return;
-  srand (time(NULL));
-
-  std::vector<TreeNode*> listPrev, listCur;
-  std::vector<int> listCurLeaves, listPrevLeaves;
-
-  _root = new TreeNode(0,0,NULL,0);
-  listPrev.push_back(_root);
-  listPrevLeaves.push_back(nbLeaves);
-  
-  int depth = 1;
-  while (true)
-  {
-    listCur.clear();
-    listCurLeaves.clear();
-    int inode=0;
-    int numChilds= numChildMax;
-    if (sameNbChildLayer)
-    {
-      if (randomChild)
-      {
-        numChilds =  rand() % (numChildMax-1) + 2;
-      }
-    }
-    
-    for (int i =0; i< listPrev.size(); i++)
-    {
-      TreeNode* np = listPrev[i];
-      int indexes = listPrevLeaves[i];
-      
-      if (!sameNbChildLayer)
-      {
-        if (randomChild)
-        {
-          numChilds =  rand() % (numChildMax-1) + 2;
-        }
-      }
-      
-      std::vector<int > partIndexes;
-      randomDistribution(indexes, numChilds, leafPerNode, partIndexes);
-      
-      //
-      int trueNbChilds = partIndexes.size();
-      Msg::Info("trueNbChilds = %d",trueNbChilds);
-      np->childs.resize(trueNbChilds);
-      np->direction.resize((trueNbChilds-1)*numDirVars);
-      for (int j=0; j< trueNbChilds; j++)
-      {
-        if (partIndexes[j] == 1)
-        {
-          // leaf
-          np->childs[j] = new TreeNode(depth,inode,np,j,affunc);
-        }
-        else if (partIndexes[j] > 1)
-        {
-          np->childs[j] = new TreeNode(depth,inode,np,j);
-          listCurLeaves.push_back(partIndexes[j]);
-          listCur.push_back(np->childs[j]);
-        }
-        inode++;
-      }
-    }
-    depth ++;
-    listPrev = listCur;
-    listPrevLeaves = listCurLeaves;
-    if (listCur.size() ==0)
-    {
-      break;
-    }
-  };
-  std::vector<TreeNode*> allLeaves;
-  getRefToAllLeaves(allLeaves);
-  std::set<TreeNode*> setNodes;
-  for (int i=0; i< allLeaves.size(); i++)
-  {
-    TreeNode* leaf = allLeaves[i];
-    if (setNodes.find(leaf) == setNodes.end())
-    {
-      TreeNode* parent = leaf->parent;
-      std::vector<TreeNode*> ll;
-      for (int j=0; j< parent->childs.size(); j++)
-      {
-        TreeNode* nj = parent->childs[j];
-        if (nj->childs.size() == 0)
-        {
-          ll.push_back(nj);
-          setNodes.insert(nj);
-        }
-      }
-      if (ll.size() == leafPerNode)
-      {
-        for (int j=0; j< ll.size(); j++)
-        {
-          ll[j]->matIndex  = j;
-        }
-      }
-      else
-      {
-        std::vector<int> v(leafPerNode,0);
-        for (int j=0; j< leafPerNode; j++)
-        {
-          v[j] = j;
-        }
-        for (int j=0; j< ll.size(); j++)
-        {
-          int pos =  rand() % (v.size()); // random from 0 to v.size() -1
-          TreeNode* otherNode = ll[j];
-          otherNode->matIndex = v[pos];
-          v.erase(v.begin() + pos);
-        }
-      }
-    }
-  }
-};
-
 void Tree::createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int numDirVars,const std::string affunc)
 {
   clear();
@@ -661,7 +170,8 @@ void Tree::createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int n
   {
     if (depth == 0)
     {
-      _root = new TreeNode(0,0,NULL,0);
+      _root = new TreeNode(0,0,NULL,numDirVars,0);
+      _nodeMap[std::pair<int,int>(0,0)]=_root;
       listPrev.push_back(_root);
     }
     else
@@ -682,17 +192,17 @@ void Tree::createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int n
         }
         //Msg::Info("numchild = %d",numChilds);
         np->childs.resize(numChilds);
-        np->direction.resize((numChilds-1)*numDirVars);
         for (int j=0; j< numChilds; j++)
         {
           if (depth == maxDepth)
           {
-            np->childs[j] = new TreeNode(depth,inode,np,j,affunc);
+            np->childs[j] = new TreeNode(depth,inode,np,numDirVars,j,affunc);
           }
           else
           {
-            np->childs[j] = new TreeNode(depth,inode,np,j);
+            np->childs[j] = new TreeNode(depth,inode,np,numDirVars,j);
           }
+          _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
           listCur.push_back(np->childs[j]);
           inode++;
         }
@@ -721,7 +231,8 @@ void Tree::createPerfectTree(int maxDepth, int numChildRegular, int numChildLeaf
   {
     if (depth == 0)
     {
-      _root = new TreeNode(0,0,NULL,0);
+      _root = new TreeNode(0,0,NULL,numDirVars,0);
+      _nodeMap[std::pair<int,int>(0,0)]=_root;
       listPrev.push_back(_root);
     }
     else
@@ -737,17 +248,17 @@ void Tree::createPerfectTree(int maxDepth, int numChildRegular, int numChildLeaf
           numChilds = numChildLeaf;
         }
         np->childs.resize(numChilds);
-        np->direction.resize((numChilds-1)*numDirVars);
         for (int j=0; j< numChilds; j++)
         {
           if (depth == maxDepth)
           {
-            np->childs[j] = new TreeNode(depth,inode,np,j,affunc);
+            np->childs[j] = new TreeNode(depth,inode,np,numDirVars,j,affunc);
           }
           else
           {
-            np->childs[j] = new TreeNode(depth,inode,np,j);
+            np->childs[j] = new TreeNode(depth,inode,np,numDirVars,j);
           }
+          _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
           listCur.push_back(np->childs[j]);
           inode++;
         }
@@ -777,7 +288,8 @@ void Tree::createPerfectTree(int maxDepth, int numChilds, int numDirVars,const s
   {
     if (depth == 0)
     {
-      _root = new TreeNode(0,0,NULL,0);
+      _root = new TreeNode(0,0,NULL,numDirVars,0);
+      _nodeMap[std::pair<int,int>(0,0)]=_root;
       listPrev.push_back(_root);
     }
     else
@@ -788,17 +300,17 @@ void Tree::createPerfectTree(int maxDepth, int numChilds, int numDirVars,const s
       {
         TreeNode* np = listPrev[i];
         np->childs.resize(numChilds);
-        np->direction.resize((numChilds-1)*numDirVars);
         for (int j=0; j< numChilds; j++)
         {
           if (depth == maxDepth)
           {
-            np->childs[j] = new TreeNode(depth,inode,np,j,affunc);
+            np->childs[j] = new TreeNode(depth,inode,np,numDirVars,j,affunc);
           }
           else
           {
-            np->childs[j] = new TreeNode(depth,inode,np,j);
+            np->childs[j] = new TreeNode(depth,inode,np,numDirVars,j);
           }
+          _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j];
           listCur.push_back(np->childs[j]);
           inode++;
         }
@@ -854,7 +366,7 @@ void Tree::loadDataFromFile(std::string filename)
 {
   // clear tree
   clear();
-  std::map<std::pair<int,int>, TreeNode* > nodeMap;
+  _nodeMap.clear();
   printf("loading data from file %s\n",filename.c_str());
   //
   FILE* fp = fopen(filename.c_str(),"r");
@@ -868,8 +380,8 @@ void Tree::loadDataFromFile(std::string filename)
       int depth, location, matIndex, childOrder, parentIdLoc, numChild, numDir;
       if (fscanf(fp, "%d %d %d %d %d %d %d", &depth, &location, &matIndex, &childOrder, &parentIdLoc, &numChild, &numDir)==7)
       {
-        printf("--------depth=%d, location=%d, childOrder=%d, parentIdLoc=%d, numChild=%d, numDir=%d \n",
-                depth, location, childOrder, parentIdLoc, numChild, numDir);
+        //printf("depth=%d, location=%d, childOrder=%d, parentIdLoc=%d, numChild=%d, numDir=%d \n",
+        //        depth, location, childOrder, parentIdLoc, numChild, numDir);
       }
       else
       {
@@ -878,21 +390,21 @@ void Tree::loadDataFromFile(std::string filename)
       std::vector<int> allChilds(numChild);
       for (int i=0; i< numChild; i++)
       {
-        int fl=fscanf(fp, "%d",&allChilds[i]);
-        printf("child %d location %d\n",i,allChilds[i]);
+        fscanf(fp, "%d",&allChilds[i]);
+        //printf("child %d location %d\n",i,allChilds[i]);
       }
       double weight;
-      int fl=fscanf(fp, "%lf",&weight);
-      printf("weight = %f\n",weight);
+      fscanf(fp, "%lf",&weight);
+      //printf("weight = %f\n",weight);
       std::vector<double> direction(numDir);
       for (int i=0; i< numDir; i++)
       {
-        fl=fscanf(fp, "%lf",&direction[i]);
-        printf("direction %d value %f\n",i,direction[i]);
+        fscanf(fp, "%lf",&direction[i]);
+        //printf("direction %d value %f\n",i,direction[i]);
       }
       char what[256];
-      fl=fscanf(fp, "%s", what);
-      printf("aftype = %s\n",what);
+      fscanf(fp, "%s", what);
+      //printf("aftype = %s\n",what);
       
       // create node
       if (depth == 0)
@@ -915,7 +427,7 @@ void Tree::loadDataFromFile(std::string filename)
           node->childs.clear();
         }
         _root = node;
-        nodeMap[std::pair<int,int>(depth,location)]=node;
+        _nodeMap[std::pair<int,int>(depth,location)]=node;
       }
       else
       {
@@ -926,13 +438,7 @@ void Tree::loadDataFromFile(std::string filename)
         node->depth = depth;
         node->location = location;
         node->childOrder = childOrder;
-        std::pair<int,int> twoNum(depth-1,parentIdLoc);
-        if (nodeMap.find(twoNum) == nodeMap.end())
-        {
-          Msg::Error("node %d %d is not found !!!",depth-1,parentIdLoc);
-          Msg::Exit(0);
-        }
-        node->parent = nodeMap[twoNum];
+        node->parent = _nodeMap[std::pair<int,int>(depth-1,parentIdLoc)];
         node->parent->childs[childOrder] = node;
         node->matIndex = matIndex;
         if (numChild > 0)
@@ -943,9 +449,9 @@ void Tree::loadDataFromFile(std::string filename)
         {
           node->childs.clear();
         }
-        printf("node %d %d:\n",depth,location);
-        node->parent->printData();
-        nodeMap[std::pair<int,int>(depth,location)]=node;
+        //printf("node %d %d:\n",depth,location);
+        //node->parent->printData();
+        _nodeMap[std::pair<int,int>(depth,location)]=node;
       }
     }
     fclose(fp);
@@ -966,74 +472,21 @@ int Tree::getNumberOfNodes() const
   return allNodes.size();
 };
 
-bool Tree::checkDuplicate() const
+bool Tree::removeLeavesWithZeroContribution()
 {
   bool ok = false;
-  nodeContainer allNodes;
-  getAllNodes(allNodes);
-  std::set<int> nodeId;
-  for (int i=0; i< allNodes.size(); i++)
-  {
-    const TreeNode* node = allNodes[i];
-    std::set<int>::const_iterator itF = nodeId.find(node->getNodeId());
-    if (itF == nodeId.end())
-    {
-      nodeId.insert(node->getNodeId());
-    }
-    else
-    {
-      Msg::Error("node duplicating:");
-      node->printData();
-      ok = true;
-    }
-  }
-  if (!ok)
-  {
-    Msg::Info("no duplicatint node is found !!!");
-  }
-  return ok;
-}
-
-bool Tree::removeZeroLeaves(double tol, bool iteration)
-{
-  bool ok = false;
-  int maxInter = 100;
-  int it = 0;
-  while (true)
+  std::vector<TreeNode*> allLeaves;
+  getRefToAllLeaves(allLeaves);
+  for (int i=0; i< allLeaves.size(); i++)
   {
-    it ++;
-    if (it == maxInter)
-    {
-      Msg::Error("removing duplicating leaves: maximal number of iterations (%d) reaches", maxInter);
-      Msg::Exit(0);
-    }
-    std::vector<TreeNode*> allLeaves;
-    getRefToAllLeaves(allLeaves);
-    std::vector<TreeNode*> allLeavesWillRemove;
-    for (int i=0; i< allLeaves.size(); i++)
-    {
-      TreeNode* leaf = allLeaves[i];
-      if (leaf->af->getVal(leaf->weight) <tol*_root->weight)
-      {
-        allLeavesWillRemove.push_back(leaf);
-      }
-    }
-    if (allLeavesWillRemove.size() == 0)
-    {
-      Msg::Info("done removing !!! ");
-      return ok; 
-    }
-    else
-    {
-      Msg::Info("iter = %d /%d",it,maxInter);
-      ok = true;
-    }
-    for (int i=0; i< allLeavesWillRemove.size(); i++)
+    TreeNode* leaf = allLeaves[i];
+    TreeNode* parent = leaf->parent;
+    double tol = 1e-6;
+    if (leaf->af->getVal(leaf->weight) <tol*_root->weight)
     {
-      TreeNode* leaf = allLeavesWillRemove[i];
-      TreeNode* parent = leaf->parent;
       if (parent->childs.size() ==2)
       {
+        ok = true;
         TreeNode* grandParent = parent->parent;
         TreeNode* otherChild = NULL;
         if (leaf->childOrder == 0)
@@ -1044,406 +497,58 @@ bool Tree::removeZeroLeaves(double tol, bool iteration)
         {
           otherChild =  parent->childs[0];
         }
+        Msg::Info("otherChild weight = %e parent weight = %e",otherChild->weight,parent->weight);
         
         otherChild->depth = parent->depth;
         otherChild->location = parent->location;
         otherChild->parent = parent->parent;
         otherChild->childOrder = parent->childOrder;
         grandParent->childs[parent->childOrder] = otherChild;
-        
-        // decrease all associated nodes by 1
-        std::vector<TreeNode*> allAssociatedNodesToChild;
-        getRefToAssociatedNodes(otherChild,allAssociatedNodesToChild);
-        for (int j=0; j< allAssociatedNodesToChild.size(); j++)
-        {
-          std::vector<int> loc;
-          getMaxLocationByDepth(loc);
-          allAssociatedNodesToChild[j]->depth--;
-          allAssociatedNodesToChild[j]->location=loc[allAssociatedNodesToChild[j]->depth]+1;
-        };
-        
-        Msg::Info("child weight = %e otherChild weight = %e parent weight = %e",leaf->af->getVal(leaf->weight),otherChild->af->getVal(otherChild->weight),parent->weight);
       }
-      else if (parent->childs.size() > 2)
+      else
       {
         Msg::Info("----------------------------------");
         parent->printData();
-        //
-        std::vector<TreeNode*> newChilds;
-        std::vector<double> newDirection;
-        //
-        int numberNormal = parent->childs.size()-1;
-        int totalNumberDirVars = parent->direction.size();
-        int numberVarsPerNormal = totalNumberDirVars/numberNormal;
-        
-        for (int j=0; j< parent->childs.size(); j++)
-        {
-          if (parent->childs[j]->location != leaf->location)
-          {
-            newChilds.push_back(parent->childs[j]);
-            if (newChilds.size() > 1)
-            {
-              if (numberVarsPerNormal == 1)
-              {
-                newDirection.push_back(parent->direction[j-1]);
-              }
-              else if (numberVarsPerNormal == 2)
-              {
-                newDirection.push_back(parent->direction[(j-1)*numberVarsPerNormal]);
-                newDirection.push_back(parent->direction[(j-1)*numberVarsPerNormal+1]);
-              }
-              else
-              {
-                Msg::Error("numberVarsPerNormal = %d is not correct !!!",numberVarsPerNormal);
-                Msg::Exit(0);
-              }
-            }
-          }
-        }
-        Msg::Info("\nSTART removing ");
-        parent->printData();
-        parent->childs = newChilds;
-        for (int j=0; j <parent->childs.size(); j++)
+        for (int i=0; i< parent->childs.size(); i++)
         {
-          parent->childs[j]->childOrder = j;
+          parent->childs[i]->printData();
         }
-        parent->direction = newDirection;   
-        Msg::Info("------------------");
-        parent->printData();     
-        Msg::Info("DONE removing\n");
-      }
-      else
-      {
-        Msg::Error("num childs must be larger than 2");
-        Msg::Exit(0);
-      }
-    };
-    
-    if (!iteration)
-    {
-      break;
-    }
-  };
-  return false;
-};
-
-bool Tree::treeMerging()
-{
-  bool ok = false;
-  int itMerge = 0;
-  int maxInter = 100;
-  while (true)
-  {
-    itMerge ++;
-    if (itMerge == maxInter)
-    {
-      Msg::Error("maximal number of iterations (%d) reaches", maxInter);
-      Msg::Exit(0);
-    }
-    std::vector<TreeNode*> allLeaves;
-    getRefToAllLeaves(allLeaves);
-    std::set<TreeNode*> allParents;
-    std::set<TreeNode*> allParentsOneChild;
-    for (int i=0; i< allLeaves.size(); i++)
-    {
-      TreeNode* leaf = allLeaves[i];
-      TreeNode* parent = leaf->parent;
-      if (parent->childs.size() > 1)
-      {
-        std::map<int, std::vector<TreeNode*> > mapIndexes;
-        for (int j=0; j< parent->childs.size(); j++)
-        {
-          std::vector<TreeNode*>& cc = mapIndexes[parent->childs[j]->matIndex];      
-          cc.push_back(parent->childs[j]);
-        }
-        for (std::map<int, std::vector<TreeNode*> >::iterator it = mapIndexes.begin(); it != mapIndexes.end(); it++)
-        {
-          int mIndex = it->first;
-          if (mIndex >=0)
-          {
-            std::vector<TreeNode*>& sameMat = it->second;
-            if (sameMat.size()> 1)
-            {
-              allParents.insert(parent);
-            }
-          }
-        }
-      }
-      else if (parent->childs.size() == 1)
-      {
-        allParentsOneChild.insert(parent);
-      }
-    }
-    
-    if (allParents.size() == 0 && allParentsOneChild.size() == 0)
-    {
-      Msg::Info("done tree merging");
-      return ok;
-    }
-    else
-    {
-      Msg::Info("tree merging %d / %d",itMerge,maxInter);
-      ok = true;
-    }
-    
-    for (std::set<TreeNode*>::iterator itp = allParentsOneChild.begin(); itp != allParentsOneChild.end(); itp++)
-    {
-      TreeNode* parent = *itp;
-      TreeNode* node = parent->childs[0];
-      TreeNode* grandParent = parent->parent;
-      //
-      node->depth = parent->depth;
-      node->location = parent->location;
-      node->parent = parent->parent;
-      node->childOrder = parent->childOrder;
-      grandParent->childs[parent->childOrder] = node;
-      
-      //
-      // decrease all associated nodes by 1
-      std::vector<TreeNode*> allAssociatedNodesToChild;
-      getRefToAssociatedNodes(node,allAssociatedNodesToChild);
-      for (int j=0; j< allAssociatedNodesToChild.size(); j++)
-      {
-        std::vector<int> loc;
-        getMaxLocationByDepth(loc);
-        allAssociatedNodesToChild[j]->depth--;
-        allAssociatedNodesToChild[j]->location=loc[allAssociatedNodesToChild[j]->depth]+1;
-      };
-    };
-    
-    for (std::set<TreeNode*>::iterator itp = allParents.begin(); itp != allParents.end(); itp++)
-    {
-      TreeNode* parent = *itp;
-      std::map<int, std::vector<std::pair<TreeNode*, int> > > mapIndexes;
-      for (int j=0; j< parent->childs.size(); j++)
-      {
-        std::vector<std::pair<TreeNode*, int > >& cc = mapIndexes[parent->childs[j]->matIndex];      
-        cc.push_back(std::pair<TreeNode*, int>(parent->childs[j],j));
-      }
-      for (std::map<int, std::vector<std::pair<TreeNode*, int> > >::iterator it = mapIndexes.begin(); it != mapIndexes.end(); it++)
-      {
-        int mIndex = it->first;
-        std::vector<std::pair<TreeNode*, int> >& sameMat = it->second;
-        if ((mIndex >=0) && (sameMat.size()> 1))
-        {
-          std::vector<TreeNode*> nodes;
-          std::set<int> position;
-          for (int k=0; k< sameMat.size(); k++)
-          {
-            nodes.push_back(sameMat[k].first);
-            if (k > 0)
-            {
-              position.insert(sameMat[k].second);
-            }
-          };
-          
-          double totalW = 0.;
-          for (int k=0; k< nodes.size(); k++)
-          {
-            totalW += nodes[k]->af->getVal(nodes[k]->weight);
-          };
-          nodes[0]->weight = nodes[0]->af->getReciprocalVal(totalW);
-          
-          int numberNormal = parent->childs.size()-1;
-          int totalNumberDirVars = parent->direction.size();
-          int numberVarsPerNormal = totalNumberDirVars/numberNormal;
         
-          std::vector<TreeNode*> newChilds;
-          std::vector<double> newNormals;
-          for (int j=0; j< parent->childs.size(); j++)
-          {
-            if (position.find(j) == position.end())
-            {
-              newChilds.push_back(parent->childs[j]);
-              if (j > 0)
-              {
-                if (numberVarsPerNormal == 1)
-                {
-                  newNormals.push_back(parent->direction[j-1]);
-                }
-                else if (numberVarsPerNormal == 2)
-                {
-                  newNormals.push_back(parent->direction[(j-1)*numberVarsPerNormal]);
-                  newNormals.push_back(parent->direction[(j-1)*numberVarsPerNormal+1]);
-                }              
-              }
-            }
-          };
-          
-          Msg::Info("\nSTART merging ");
-          parent->printData();
-          printf("merge loc (number of phases = %ld): %d (mat %d) ",parent->childs.size(),sameMat[0].second, (sameMat[0].first)->matIndex);
-          for (std::set<int>::iterator itpos = position.begin(); itpos != position.end(); itpos++)
-          {
-            printf(" %d (mat %d)",*itpos,parent->childs[*itpos]->matIndex);
-          }
-          printf("\n----------------------\n");
-          parent->childs = newChilds;
-          for (int j=0; j <parent->childs.size(); j++)
-          {
-            parent->childs[j]->childOrder = j;
-          }
-          parent->direction = newNormals;
-          parent->printData();
-          Msg::Info("DONE merging\n");
-          //
-          break;
-        }
-      }
-    };
-  }
-  return false;
-};
-
-bool Tree::removeZeroBranches(double tol)
-{
-  int depth = 0;
-  bool ok = false;
-  while (true)
-  {
-    std::map<int,std::vector<TreeNode*> > allNodesByDepth;
-    getAllNodesByDepth(allNodesByDepth);
-    std::map<int,std::vector<TreeNode*> >::iterator itF= allNodesByDepth.find(depth);
-    if ( itF!= allNodesByDepth.end())
-    {
-      std::vector<TreeNode*>& depthNodes = itF->second;
-      for (int i=0; i<depthNodes.size(); i++)
-      {
-        TreeNode* node = depthNodes[i];
-        double Wt = getNonNegativeWeight(node);
-        //
-        bool removed = false;
         std::vector<TreeNode*> newChilds;
-        std::vector<double> newDirection;
-        std::vector<int> removeLoc;
-        //
-        int numberNormal = node->childs.size()-1;
-        int totalNumberDirVars = node->direction.size();
-        int numberVarsPerNormal = totalNumberDirVars/numberNormal;
-        
-        for (int j=0; j< node->childs.size(); j++)
-        {
-          double Wchild = getNonNegativeWeight(node->childs[j]);
-          if (Wchild < tol* Wt)
-          {
-            removed = true;
-            removeLoc.push_back(j);
-          }
-          else
-          {
-            newChilds.push_back(node->childs[j]);
-            if (newChilds.size() > 1)
-            {
-              if (numberVarsPerNormal == 1)
-              {
-                newDirection.push_back(node->direction[j-1]);
-              }
-              else if (numberVarsPerNormal == 2)
-              {
-                newDirection.push_back(node->direction[(j-1)*numberVarsPerNormal]);
-                newDirection.push_back(node->direction[(j-1)*numberVarsPerNormal+1]);
-              }
-              else
-              {
-                Msg::Error("numberVarsPerNormal = %d is not correct !!!",numberVarsPerNormal);
-                Msg::Exit(0);
-              }
-            }
-          }
-        }
-        if (removed)
+        bool found = false;
+        for (int i=0; i< parent->childs.size(); i++)
         {
-          Msg::Info("\nSTART removing ");
-          node->printData();
-          ok = true;
-          printf("pos (number of phases= %ld): ",node->childs.size());
-          for (int j=0; j< removeLoc.size(); j++)
+          if (parent->childs[i] == leaf)
           {
-            printf("%d (f=%e) ",removeLoc[j],getNonNegativeWeight(node->childs[removeLoc[j]])/Wt);
+            found = true;
           }
-          printf("\n");
-          if (newChilds.size() >1)
+          if (!found)
           {
-            node->childs = newChilds;
-            for (int j=0; j <node->childs.size(); j++)
-            {
-              node->childs[j]->childOrder = j;
-            }
-            node->direction = newDirection;   
-            Msg::Info("------------------");
-            node->printData();     
+            newChilds.push_back(parent->childs[i]);
           }
           else
           {
-            TreeNode* grandParent = node->parent;
-            TreeNode* child = newChilds[0];
-            //
-            child->depth = node->depth;
-            child->location = node->location;
-            child->parent = node->parent;
-            child->childOrder = node->childOrder;
-            grandParent->childs[node->childOrder] = child;
-            // decrease all associated nodes by 1
-            std::vector<TreeNode*> allAssociatedNodesToChild;
-            getRefToAssociatedNodes(child,allAssociatedNodesToChild);
-            for (int j=0; j< allAssociatedNodesToChild.size(); j++)
+            for (int j=i+1; j<parent->childs.size(); j++)
             {
-              std::vector<int> loc;
-              getMaxLocationByDepth(loc);
-              allAssociatedNodesToChild[j]->depth--;
-              allAssociatedNodesToChild[j]->location=loc[allAssociatedNodesToChild[j]->depth]+1;
+              TreeNode* n = parent->childs[j];
+              n->childOrder = j-1;
+              newChilds.push_back(n);
             };
-            
-            Msg::Info("------------------");
-            Msg::Info("node is replaced by child");
+            break;
           }
-          Msg::Info("DONE removing\n");
+        };
+        parent->childs = newChilds;
+        Msg::Info("after removing");
+        parent->printData();
+        for (int i=0; i< parent->childs.size(); i++)
+        {
+          parent->childs[i]->printData();
         }
+        Msg::Info("----------------------------------");
       }
     }
-    else
-    {
-      break;
-    }
-    depth++;
-  };
-  return ok;
-};
-
-bool Tree::removeLeavesWithZeroContribution(double tol)
-{
-  int iterRemove = 0;
-  int maxIterRemove = 100;
-  bool ok = false;
-  while (true)
-  {
-    iterRemove ++;
-    if (iterRemove == maxIterRemove)
-    {
-      Msg::Error("removeLeavesWithZeroContribution: maximal number of iterations (%d) reaches", maxIterRemove);
-      Msg::Exit(0);
-    }
-    //bool okLeaf = removeZeroBranches(tol);
-    bool okLeaf = removeZeroLeaves(tol,true);
-    bool okMerge = treeMerging();
-    // check
-    if (okLeaf || okLeaf)
-    {
-      ok = true;
-    }
-    
-    if ((!okLeaf)&& (!okMerge))
-    {
-      Msg::Info("done removing");
-      return ok;
-    }
-    else
-    {
-      Msg::Info("removing %d / %d",iterRemove,maxIterRemove);
-    }
-    
   }
-  return false;
+  return ok;
 };
 
 double Tree::getPhaseFraction(int i) const
@@ -1459,21 +564,7 @@ double Tree::getPhaseFraction(int i) const
     total += node->af->getVal(node->weight);
   }
   return mapFraction[i]/total;
-};
-
-void Tree::printNodeFraction() const
-{
-  nodeContainer allLeaves;
-  getAllLeaves(allLeaves);
-  double total = 0;
-  for (int i=0; i< allLeaves.size(); i++)
-  {
-    const TreeNode* node = allLeaves[i];
-    Msg::Info("fraction node %d = %e",node->getNodeId(), node->af->getVal(node->weight)/_root->weight);
-    total += node->af->getVal(node->weight);
-  }
-  Msg::Info("Wroot = %e, total = %e",_root->weight,total);
-};
+}
 
 void Tree::printPhaseFraction() const
 {
@@ -1495,27 +586,11 @@ void Tree::printPhaseFraction() const
   allPhase.print("all phase fraction:"); 
 };
 
-void Tree::printLeafFraction() const
-{
-  nodeContainer allLeaves;
-  getAllLeaves(allLeaves);
-  double total = getNonNegativeWeight(_root);
-  printf("print fraction at leaf:\n");
-  int ord=0;
-  for (int i=0; i< allLeaves.size(); i++)
-  {
-    const TreeNode* leaf = allLeaves[i];
-    printf("loc %d: %d %d, fraction = %.5e\n",ord,leaf->depth,leaf->location,getNonNegativeWeight(leaf)/total);
-    ord++;
-  }
-  printf("done printing:\n"); 
-};
-
 void Tree::printTree() const
 {
   printf("Tree print \n");
   if (_root == NULL) return;
-  std::vector< TreeNode*> listPrev(0), listCur(0);
+  std::vector< TreeNode*> listPrev, listCur;
   listPrev.push_back(_root);
   int depth = 0;
   while (true)
@@ -1542,7 +617,7 @@ void Tree::printTree() const
   
   nodeContainer allLeaves;
   getAllLeaves(allLeaves);
-  printf("all leaves = %ld\n",allLeaves.size());
+  printf("all leaves = %d\n",allLeaves.size());
   for (int i=0; i< allLeaves.size(); i++)
   {
     allLeaves[i]->printData();
@@ -1550,140 +625,6 @@ void Tree::printTree() const
   printf("\n");
 };
 
-void Tree::printTreeInteraction(const std::string fname, bool colorMat, int dir) const
-{
-  printf("Tree print \n");
-  if (_root == NULL) return;
-  std::vector< TreeNode*> listPrev, listCur;
-  listPrev.push_back(_root);
-  nodeContainer allLeaves;
-  getAssociatedLeavesForNode(_root,allLeaves);
-  printf("*|%ld|*\n",allLeaves.size());
-  while (true)
-  {
-    listCur.clear();
-    for (int i =0; i< listPrev.size(); i++)
-    {
-      printf("*|");
-      TreeNode& np = *(listPrev[i]);
-      for (int j=0; j< np.childs.size(); j++)
-      {
-        TreeNode* nj = np.childs[j];
-        allLeaves.clear();
-        getAssociatedLeavesForNode(nj,allLeaves);
-        if (allLeaves.size() == 1)
-        {
-          if (j==0)
-          {
-            printf("%ld(%d)",allLeaves.size(),nj->matIndex);
-          }
-          else
-          {
-            printf("-%ld(%d)",allLeaves.size(),nj->matIndex);
-          }
-        }
-        else
-        {
-          if (j==0)
-          {
-            printf("%ld",allLeaves.size());
-          }
-          else
-          {
-            printf("-%ld",allLeaves.size());
-          }
-          listCur.push_back(nj);
-        }
-      }
-      printf("|*");
-    }
-    printf("\n");
-    listPrev = listCur;
-    if (listCur.size() == 0)
-    {
-      break;
-    }
-  };
-  
-  
-  std::map<const TreeNode*, std::string> nameMap;
-  nodeContainer allNodes;
-  getAllNodes(allNodes);
-  for (int i=0; i< allNodes.size(); i++)
-  {
-    nameMap[allNodes[i]] = "n00"+std::to_string(i+2);
-  };
-  
-  std::string txt ="##Command to get the layout: \"dot  -Teps "+fname+" > treeInteraction.eps\"\n";
-  txt += "graph \"\" \n{node [fontsize=60];\ngraph[fontsize=50,style=invis,ratio=1];\nsubgraph cluster01\n{";
-  
-  
-  for (std::map<const TreeNode*, std::string>::iterator it = nameMap.begin(); it != nameMap.end(); it++)
-  {
-    allLeaves.clear();
-    const TreeNode* node = it->first;
-    getAssociatedLeavesForNode(node,allLeaves);
-    std::string label;
-    if (node->childs.size() > 0)
-    {
-      label = " [label=\""+ std::to_string(allLeaves.size())+"\"];\n";
-    }
-    else
-    {
-      if (colorMat)
-      {
-        if (node->matIndex == 0)
-          label = " [label=\""+ std::to_string(allLeaves.size()) +"\",color=red,style=filled];\n";
-        else if (node->matIndex == 1)
-          label = " [label=\""+ std::to_string(allLeaves.size()) +"\",color=green,style=filled];\n";
-        else if (node->matIndex == 2)
-          label = " [label=\""+ std::to_string(allLeaves.size()) +"\",color=cyan,style=filled];\n";
-        else if (node->matIndex == 3)
-          label = " [label=\""+ std::to_string(allLeaves.size()) +"\",color=orange,style=filled];\n";
-        else
-        {
-          label = " [label=\""+ std::to_string(allLeaves.size()) +"\",color=blue,style=filled];\n";
-        }
-      }
-      else
-      {
-        label = " [label=\""+ std::to_string(allLeaves.size()) +"\",color=azure4,style=filled];\n";
-      }
-    }
-    txt += it->second+label;
-  }
-  // connection
-  txt += nameMap[_root]+"\n";
-  for (std::map<const TreeNode*, std::string>::iterator it = nameMap.begin(); it != nameMap.end(); it++)
-  {
-    const TreeNode* node = it->first;
-    for (int i=0; i< node->childs.size(); i++)
-    {
-      const TreeNode* conn = node->childs[i];
-      if (dir == 0)
-      {
-        if (conn->childs.size() > 0)
-          txt += nameMap[node] + " -- " + nameMap[conn]+"[dir=forward,arrowsize=2];\n";
-        else
-          txt += nameMap[node] + " -- " + nameMap[conn]+"[dir=forward,arrowsize=2];\n";
-      }
-      else
-      {
-        if (conn->childs.size() > 0)
-          txt += nameMap[node] + " -- " + nameMap[conn]+"[dir=back,arrowsize=2];\n";
-        else
-          txt += nameMap[node] + " -- " + nameMap[conn]+"[dir=back,arrowsize=2];\n";
-      }
-    }
-  }
- 
-  txt += "}\n}";
-  FILE* ff = fopen(fname.c_str(),"w");
-  fprintf(ff,"%s",txt.c_str());
-  fclose(ff); 
-  
-};
-
 void Tree::getAllLeaves(nodeContainer& allLeaves) const
 {
   allLeaves.clear();
@@ -1742,11 +683,11 @@ void Tree::getRefToAllLeaves(std::vector<TreeNode*>& allLeaves)
   };
 };
 
-void Tree::getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes)
+void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const
 {
   allNodes.clear();
   if (_root == NULL) return;
-  std::vector<TreeNode*> listPrev, listCur;
+  nodeContainer listPrev, listCur;
   listPrev.push_back(_root);
   //
   allNodes[_root->depth].push_back(_root);
@@ -1755,7 +696,7 @@ void Tree::getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes)
     listCur.clear();    
     for (int i =0; i< listPrev.size(); i++)
     {
-      TreeNode* np = listPrev[i];
+      const TreeNode* np = listPrev[i];
       for (int j=0; j< np->childs.size(); j++)
       {
         listCur.push_back(np->childs[j]);
@@ -1770,61 +711,24 @@ void Tree::getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes)
   };
 };
 
-void Tree::getMaxLocationByDepth(std::vector<int>& loc) const
+const TreeNode* Tree::getNode(int depth, int pos) const
 {
-  std::map<int,nodeContainer > allNodes;
-  getAllNodesByDepth(allNodes);
-  loc.resize(allNodes.size(),0);
-  for (std::map<int,nodeContainer >::const_iterator it = allNodes.begin(); it != allNodes.end(); it++)
+  std::map<std::pair<int,int>, TreeNode* >::const_iterator it = _nodeMap.find(std::pair<int,int>(depth,pos));
+  if (it == _nodeMap.end())
   {
-    const nodeContainer& depthNode = it->second;
-    for (int j=0; j< depthNode.size(); j++)
-    {
-      const TreeNode* node = depthNode[j];
-      if (node->location >loc[it->first])
-      {
-        loc[it->first]= node->location;
-      }
-    }
+    printf("node in depth %d pos %d is not found \n",depth,pos);
+    return NULL;
   }
+  return it->second;
 };
 
-void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const
-{
-  allNodes.clear();
-  if (_root == NULL) return;
-  nodeContainer listPrev, listCur;
-  listPrev.push_back(_root);
-  //
-  allNodes[_root->depth].push_back(_root);
-  while (true)
-  {
-    listCur.clear();    
-    for (int i =0; i< listPrev.size(); i++)
-    {
-      const TreeNode* np = listPrev[i];
-      for (int j=0; j< np->childs.size(); j++)
-      {
-        listCur.push_back(np->childs[j]);
-        allNodes[np->childs[j]->depth].push_back(np->childs[j]);
-      }
-    }
-    listPrev = listCur;
-    if (listCur.size() == 0)
-    {
-      break;
-    }
-  };
-};
-
-
 void Tree::printBackTrackingPath(const TreeNode* node) const
 {
   nodeContainer allNodes;
   getBackTrackingPath(node,allNodes);
   printf("-----------------------------\n");
   node->printData();
-  printf("backtracking laves to node: %ld \n",allNodes.size());
+  printf("backtracking laves to node: %d \n",allNodes.size());
   for (int i=0; i< allNodes.size(); i++)
   {
     allNodes[i]->printData();
@@ -1839,7 +743,7 @@ void Tree::printAssociatedLeavesForNode(const TreeNode* node) const
   
   printf("-----------------------------\n");
   node->printData();
-  printf("associated laves to node: %ld \n",allLeaves.size());
+  printf("associated laves to node: %d \n",allLeaves.size());
   for (int i=0; i< allLeaves.size(); i++)
   {
     allLeaves[i]->printData();
@@ -1859,54 +763,6 @@ void Tree::getBackTrackingPath(const TreeNode* node, nodeContainer& allNodes) co
   };
 };
 
-void Tree::getAssociatedLeavesForNode(TreeNode* node, std::vector<TreeNode*>& allLeaves)
-{
-  allLeaves.clear();
-  std::vector<TreeNode*> listPrev, listCur;
-  listPrev.push_back(node);
-  while (true)
-  {
-    listCur.clear();    
-    for (int i =0; i< listPrev.size(); i++)
-    {
-      TreeNode* np = listPrev[i];
-      if (np->childs.size() == 0)
-      {
-        allLeaves.push_back(np);
-      }
-      for (int j=0; j< np->childs.size(); j++)
-      {
-        listCur.push_back(np->childs[j]);
-      }
-    }
-    listPrev = listCur;
-    if (listCur.size() == 0)
-    {
-      break;
-    }
-  };
-};
-
-double Tree::getNonNegativeWeight(const TreeNode* node) const
-{
-  if (node->childs.size() == 0)
-  {
-    return node->af->getVal(node->weight);
-  }
-  else
-  {
-    nodeContainer allLeaves;
-    getAssociatedLeavesForNode(node,allLeaves);
-    double W = 0;
-    for (int i=0; i< allLeaves.size(); i++)
-    {
-      const TreeNode* leaf = allLeaves[i];
-      W += leaf->af->getVal(leaf->weight);
-    };
-    return W;
-  }
-}
-
 void Tree::getAssociatedLeavesForNode(const TreeNode* node, nodeContainer& allLeaves) const
 {
   allLeaves.clear();
@@ -1989,12 +845,13 @@ void Tree::assignMaterialLaws(int numPhases)
   }
 }
 
+
 void Tree::printAllRegularNodes() const
 {
   nodeContainer allNodes;
   getAllRegularNodes(allNodes);
   printf("-----------------------------\n");
-  printf("all regular nodes: %ld \n",allNodes.size());
+  printf("all regular nodes: %d \n",allNodes.size());
   for (int i=0; i< allNodes.size(); i++)
   {
     allNodes[i]->printData();
@@ -2065,35 +922,6 @@ void Tree::getAllNodes(nodeContainer& allNodes) const
   };
 };
 
-void Tree::getRefToAssociatedNodes(TreeNode* node, std::vector<TreeNode*>& allNodes)
-{
-  allNodes.clear();
-  std::vector<TreeNode*> listPrev, listCur;
-  listPrev.push_back(node);
-  //
-  while (true)
-  {
-    listCur.clear();    
-    for (int i =0; i< listPrev.size(); i++)
-    {
-      TreeNode* np = listPrev[i];
-      for (int j=0; j< np->childs.size(); j++)
-      {
-        listCur.push_back(np->childs[j]);
-      }
-    }
-    listPrev = listCur;
-    if (listCur.size() == 0)
-    {
-      break;
-    }
-    else
-    {
-      allNodes.insert(allNodes.end(),listCur.begin(),listCur.end());
-    }
-  };
-}
-
 void Tree::getRefToAllNodes(std::vector<TreeNode*>& allNodes)
 {
   allNodes.clear();
@@ -2138,73 +966,7 @@ void Tree::clear()
   _root = NULL;
 };
 
-#if 0
-void Tree::initialize(bool rand, bool normalizedWeight)
-{
-  // random weight for leaves
-  unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
-  std::default_random_engine generator(seed);
-  std::uniform_real_distribution<double> distribution(0., 1.); 
-  
-  auto getVal = [&](double min, double max)
-  {
-    double x = 0.5;
-    if (rand)
-    {
-      x= distribution(generator);
-    }
-    return (max-min)*x + min;
-  };
-  
-  double minWeight = 0.2;
-  double maxWeight = 0.8;
-  std::vector<TreeNode*> allLeaves;
-  getRefToAllLeaves(allLeaves);
-  double toltalWeightLeaves = 0.;
-  for (int i=0; i< allLeaves.size(); i++)
-  {
-    TreeNode* n = allLeaves[i];
-    n->weight = getVal(minWeight,maxWeight);
-    toltalWeightLeaves += n->af->getVal(n->weight);
-  }
-  
-  for (int i=0; i< allLeaves.size(); i++)
-  {
-    TreeNode* n = allLeaves[i];
-    if (normalizedWeight)
-    {
-      double vv = (n->af->getVal(n->weight))/toltalWeightLeaves;
-      //n->weight = n->af->getReciprocalVal(vv);   
-      n->weight = vv;
-    }
-    // propagate data
-    double w = n->weight;
-    TreeNode* parent = n->parent;
-    while (parent != NULL)
-    {
-      parent->weight += n->af->getVal(w);
-      parent = parent->parent;
-    };
-  };
-  //printf("pi = %e\n",pi);
-  double minAlpha = 0.;
-  double maxAlpha = 1;
-  //
-  std::vector<TreeNode*> allNodes;
-  getRefToAllNodes(allNodes);
-  printf("num of nodes = %ld\n",allNodes.size());
-  for (int i=0; i< allNodes.size(); i++)
-  {
-    TreeNode& node = *(allNodes[i]);
-    for (int j=0; j<node.direction.size(); j++)
-    {
-      node.direction[j] = getVal(minAlpha,maxAlpha);   
-    }
-  };
-};
-
-#else
-void Tree::initialize(bool rand,bool normalizedWeight)
+void Tree::initialize(bool rand)
 {
   // random weight for leaves
   unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
@@ -2238,11 +1000,10 @@ void Tree::initialize(bool rand,bool normalizedWeight)
     n->weight *= (1./toltalWeightLeaves);
     // propagate data
     double w = n->weight;
-    TreeNode* pp = n->parent;
-    while (pp != NULL)
+    while (n->parent != NULL)
     {
-      pp->weight += n->af->getVal(w);
-      pp = pp->parent;
+      n->parent->weight += n->af->getVal(w);
+      n = n->parent;
     };
   };
     
@@ -2263,4 +1024,3 @@ void Tree::initialize(bool rand,bool normalizedWeight)
   };
 };
 
-#endif //
\ No newline at end of file
diff --git a/NonLinearSolver/modelReduction/Tree.h b/NonLinearSolver/modelReduction/Tree.h
index 2f775afd5068e61beaa40ecd4409930e1fae1493..fbe574db22991f4f566461b9f881199bee1b07b5 100644
--- a/NonLinearSolver/modelReduction/Tree.h
+++ b/NonLinearSolver/modelReduction/Tree.h
@@ -13,15 +13,13 @@
 #ifndef SWIG
 #include "ANNUtils.h"
 #include <vector>
-#include "SVector3.h"
 #endif //SWIG
 
 class TreeNode
 {
   public:
-    #ifndef SWIG
     activationFunction* af; // to make sure weight is posible when upscale 
-    std::vector<double> direction; // direction data = (nuChilds-1) in 2D and (nuChilds-1)*2 in 3D 
+    std::vector<double> direction; // direction data
     double weight; 
     int depth; // in layer 0 - root, 
     int location; // location in vector of childs of parent node
@@ -29,16 +27,12 @@ class TreeNode
     TreeNode* parent; // parent node (NULL for root node)
     std::vector<TreeNode*> childs; // child nodes (NULL for last nodes)
     int matIndex; // specify material number -1 by default
-    #endif //SWIG
     TreeNode();
-    TreeNode(int k, int c, TreeNode* p, int o, const std::string affunc="relu");
+    TreeNode(int k, int c, TreeNode* p, int numDirVars, int o, const std::string affunc="relu");
     virtual ~TreeNode();
     int getNodeId() const{return location*1000+depth;};
     void printData() const;
-    #ifndef SWIG
     void saveToFile(FILE* file) const;
-    void getNormal(std::vector<SVector3>& vec, bool stiff=false, std::vector<std::vector<SVector3> >* DnormalDunknown=NULL) const;
-    #endif //SWIG
 };
 
 class Tree
@@ -52,19 +46,16 @@ class Tree
     void createPerfectTree(int maxDepth, int numChild, int numDirVars, const std::string affunc="relu");
     void createPerfectTree(int maxDepth, int numChildRegular, int numChildLeaf, int numDirVars, const std::string affunc="relu", int lastRepeat=1);
     void createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int numDirVars,const std::string affunc="relu");
-    void createRandomTree(int nbLeaves, int numChildMax, int leafPerNode, int numDirVars,const std::string affunc="relu", bool randomSubdivion=true, bool randomChild=true, bool sameNbChildLayer=true);
-    void createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int leafPerNode, int numDirVars,const std::string affunc="relu", bool randomChild=true, bool randomMat=true);
     void createSpecialBinaryTree(int numNodes, int numDirVars,const std::string affunc="relu");
     void createSpecialTripleTree(int numNodes, int numDirVars,const std::string affunc="relu");
     void printTree() const;
-    void printTreeInteraction(const std::string fname="treeInteraction.txt", bool colorMat = true, int dir=1) const;
     void clear();
-    void initialize(bool rand=true, bool normalizedWeight=true);
+    void initialize(bool rand=true);
     const TreeNode* getRootNode() const {return _root;};
-    
-    void assignMaterialLaws(int numPhases);
+    const TreeNode* getNode(int depth, int pos) const;
     
     void printAllRegularNodes() const;
+    void assignMaterialLaws(int numPhases);
     
     void printAssociatedLeavesForNode(const TreeNode* node) const;
     void printBackTrackingPath(const TreeNode* node) const;
@@ -72,36 +63,26 @@ class Tree
     void saveDataToFile(std::string filename) const;
     void loadDataFromFile(std::string filename);
     int getNumberOfNodes() const;
-    bool checkDuplicate() const;
     double getPhaseFraction(int i) const;
     void printPhaseFraction() const;
-    void printLeafFraction() const;
-    void printNodeFraction() const;
-    
-    bool removeLeavesWithZeroContribution(double tol=1e-6);
-    bool removeZeroBranches(double tol);
-    bool removeZeroLeaves(double tol, bool iteration=true);
-    bool treeMerging();
     
+    bool removeLeavesWithZeroContribution();
     #ifndef SWIG
     void getAllRegularNodes(nodeContainer& allRegNodes) const;
     void getAllLeaves(nodeContainer& allLeaves) const;
     void getAllNodes(nodeContainer& allNodes) const;
     void getAllNodesByDepth(std::map<int,nodeContainer>& allNodes) const;
-    void getMaxLocationByDepth(std::vector<int>& loc) const;
-    void getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes);
     void getAssociatedLeavesForNode(const TreeNode* node, nodeContainer& allLeaves) const;
-    void getAssociatedLeavesForNode(TreeNode* node, std::vector<TreeNode*>& allLeaves);
     void getBackTrackingPath(const TreeNode* node, nodeContainer& path) const;
-    double getNonNegativeWeight(const TreeNode* node) const;
+    
     void getRefToAllLeaves(std::vector<TreeNode*>& allLeaves);
     void getRefToAllNodes(std::vector<TreeNode*>& allNodes);  
-    void getRefToAssociatedNodes(TreeNode* node, std::vector<TreeNode*>& allNodes);
     #endif //SWIG
   
   protected:
     #ifndef SWIG
     TreeNode* _root;
+    std::map<std::pair<int,int>, TreeNode* > _nodeMap;
     #endif //SWIG
 };
 #endif //_TREE_H_