diff --git a/NonLinearSolver/modelReduction/Tree.cpp b/NonLinearSolver/modelReduction/Tree.cpp
index 208855b6941784e020ac16b71e83d85fde1988cc..3d545a74ef9b2f32f81f9cb069df94f25ef0f5be 100644
--- a/NonLinearSolver/modelReduction/Tree.cpp
+++ b/NonLinearSolver/modelReduction/Tree.cpp
@@ -64,6 +64,7 @@ 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
@@ -917,7 +918,13 @@ void Tree::loadDataFromFile(std::string filename)
         node->depth = depth;
         node->location = location;
         node->childOrder = childOrder;
-        node->parent = nodeMap[std::pair<int,int>(depth-1,parentIdLoc)];
+        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->childs[childOrder] = node;
         node->matIndex = matIndex;
         if (numChild > 0)
@@ -1195,75 +1202,71 @@ bool Tree::treeMerging()
       for (std::map<int, std::vector<std::pair<TreeNode*, int> > >::iterator it = mapIndexes.begin(); it != mapIndexes.end(); it++)
       {
         int mIndex = it->first;
-        if (mIndex >=0)
+        std::vector<std::pair<TreeNode*, int> >& sameMat = it->second;
+        if ((mIndex >=0) && (sameMat.size()> 1))
         {
-          std::vector<std::pair<TreeNode*, int> >& sameMat = it->second;
-          if (sameMat.size()> 1)
+          std::vector<TreeNode*> nodes;
+          std::set<int> position;
+          for (int k=0; k< sameMat.size(); k++)
           {
-            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)
             {
-              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;
+              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);
           
-            std::vector<TreeNode*> newChilds;
-            std::vector<double> newNormals;
-            for (int j=0; j< parent->childs.size(); j++)
+          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())
             {
-              if (position.find(j) == position.end())
+              newChilds.push_back(parent->childs[j]);
+              if (j > 0)
               {
-                newChilds.push_back(parent->childs[j]);
-                if (j > 0)
+                if (numberVarsPerNormal == 1)
                 {
-                  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]);
-                  }              
+                  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;
-          }
+          };
           
+          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;
         }
       }
     };
@@ -1359,6 +1362,14 @@ bool Tree::removeZeroBranches(double tol)
             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++)
+            {
+              allAssociatedNodesToChild[j]->depth--;
+            };
+            
             Msg::Info("------------------");
             Msg::Info("node is replaced by child");
           }
@@ -1463,7 +1474,7 @@ void Tree::printTree() const
 {
   printf("Tree print \n");
   if (_root == NULL) return;
-  std::vector< TreeNode*> listPrev, listCur;
+  std::vector< TreeNode*> listPrev(0), listCur(0);
   listPrev.push_back(_root);
   int depth = 0;
   while (true)
@@ -1994,6 +2005,35 @@ 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();
diff --git a/NonLinearSolver/modelReduction/Tree.h b/NonLinearSolver/modelReduction/Tree.h
index b2a94811333a323b9c9cd4bef6896bf3af0ee6f3..7a5c631b980d7654f5ae67b0e274326e9bcc0264 100644
--- a/NonLinearSolver/modelReduction/Tree.h
+++ b/NonLinearSolver/modelReduction/Tree.h
@@ -93,6 +93,7 @@ class Tree
     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: