From 36a3ffe1e2c3f3f26a3ce9515205e8c863dfee7e Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 26 Nov 2010 15:19:08 +0000
Subject: [PATCH] First working repartitioning of parallel distributed mesh.

---
 AMDiS/CMakeLists.txt                   |   1 +
 AMDiS/bin/Makefile.am                  |   1 +
 AMDiS/bin/Makefile.in                  |  12 ++
 AMDiS/src/DataCollector.cc             |   4 +-
 AMDiS/src/Element.cc                   | 103 ++++++++---------
 AMDiS/src/Element.h                    |  12 +-
 AMDiS/src/Mesh.cc                      |   5 +-
 AMDiS/src/Mesh.h                       |   5 +-
 AMDiS/src/parallel/MeshDistributor.cc  | 124 ++------------------
 AMDiS/src/parallel/MeshDistributor.h   |  19 ++--
 AMDiS/src/parallel/MeshManipulation.cc | 150 +++++++++++++++++++++++++
 AMDiS/src/parallel/MeshManipulation.h  |  52 +++++++++
 12 files changed, 302 insertions(+), 186 deletions(-)
 create mode 100644 AMDiS/src/parallel/MeshManipulation.cc
 create mode 100644 AMDiS/src/parallel/MeshManipulation.h

diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 07cddfbe..5a0b88c5 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -180,6 +180,7 @@ if(ENABLE_PARALLEL_DOMAIN)
 	SET(PARALLEL_DOMAIN_AMDIS_SRC
                 ${SOURCE_DIR}/parallel/ParMetisPartitioner.cc
 		${SOURCE_DIR}/parallel/MeshDistributor.cc 
+		${SOURCE_DIR}/parallel/MeshManipulation.cc
 		${SOURCE_DIR}/parallel/StdMpi.cc
 		${SOURCE_DIR}/parallel/ParallelDebug.cc
 		${SOURCE_DIR}/parallel/MpiHelper.cc
diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am
index d4427386..dfd523bd 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -18,6 +18,7 @@ if USE_PARALLEL_DOMAIN_AMDIS
   $(SOURCE_DIR)/parallel/StdMpi.h $(SOURCE_DIR)/parallel/StdMpi.cc \
   $(SOURCE_DIR)/parallel/ParMetisPartitioner.h $(SOURCE_DIR)/parallel/ParMetisPartitioner.cc \
   $(SOURCE_DIR)/parallel/MeshDistributor.h $(SOURCE_DIR)/parallel/MeshDistributor.cc \
+  $(SOURCE_DIR)/parallel/MeshManipulation.h $(SOURCE_DIR)/parallel/MeshManipulation.cc \
   $(SOURCE_DIR)/parallel/ParallelDebug.h $(SOURCE_DIR)/parallel/ParallelDebug.cc \
   $(SOURCE_DIR)/parallel/ParallelProblemStatBase.h \
   $(SOURCE_DIR)/parallel/PetscSolver.h $(SOURCE_DIR)/parallel/PetscSolver.cc \
diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in
index 3f02633f..2afe484d 100644
--- a/AMDiS/bin/Makefile.in
+++ b/AMDiS/bin/Makefile.in
@@ -38,6 +38,7 @@ host_triplet = @host@
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@  $(SOURCE_DIR)/parallel/StdMpi.h $(SOURCE_DIR)/parallel/StdMpi.cc \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@  $(SOURCE_DIR)/parallel/ParMetisPartitioner.h $(SOURCE_DIR)/parallel/ParMetisPartitioner.cc \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@  $(SOURCE_DIR)/parallel/MeshDistributor.h $(SOURCE_DIR)/parallel/MeshDistributor.cc \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@  $(SOURCE_DIR)/parallel/MeshManipulation.h $(SOURCE_DIR)/parallel/MeshManipulation.cc \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@  $(SOURCE_DIR)/parallel/ParallelDebug.h $(SOURCE_DIR)/parallel/ParallelDebug.cc \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@  $(SOURCE_DIR)/parallel/ParallelProblemStatBase.h \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@  $(SOURCE_DIR)/parallel/PetscSolver.h $(SOURCE_DIR)/parallel/PetscSolver.cc \
@@ -96,6 +97,8 @@ am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
 	$(SOURCE_DIR)/parallel/ParMetisPartitioner.cc \
 	$(SOURCE_DIR)/parallel/MeshDistributor.h \
 	$(SOURCE_DIR)/parallel/MeshDistributor.cc \
+	$(SOURCE_DIR)/parallel/MeshManipulation.h \
+	$(SOURCE_DIR)/parallel/MeshManipulation.cc \
 	$(SOURCE_DIR)/parallel/ParallelDebug.h \
 	$(SOURCE_DIR)/parallel/ParallelDebug.cc \
 	$(SOURCE_DIR)/parallel/ParallelProblemStatBase.h \
@@ -260,6 +263,7 @@ am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__objects_1 = libamdis_la-StdMpi.lo \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@	libamdis_la-ParMetisPartitioner.lo \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@	libamdis_la-MeshDistributor.lo \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@	libamdis_la-MeshManipulation.lo \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@	libamdis_la-ParallelDebug.lo \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@	libamdis_la-PetscSolver.lo \
 @USE_PARALLEL_DOMAIN_AMDIS_TRUE@	libamdis_la-MpiHelper.lo \
@@ -820,6 +824,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Marker.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Mesh.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MeshDistributor.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MeshManipulation.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MeshStructure.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MpiHelper.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-NonLinUpdater.Plo@am__quote@
@@ -926,6 +931,13 @@ libamdis_la-MeshDistributor.lo: $(SOURCE_DIR)/parallel/MeshDistributor.cc
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-MeshDistributor.lo `test -f '$(SOURCE_DIR)/parallel/MeshDistributor.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/MeshDistributor.cc
 
+libamdis_la-MeshManipulation.lo: $(SOURCE_DIR)/parallel/MeshManipulation.cc
+@am__fastdepCXX_TRUE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-MeshManipulation.lo -MD -MP -MF $(DEPDIR)/libamdis_la-MeshManipulation.Tpo -c -o libamdis_la-MeshManipulation.lo `test -f '$(SOURCE_DIR)/parallel/MeshManipulation.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/MeshManipulation.cc
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/libamdis_la-MeshManipulation.Tpo $(DEPDIR)/libamdis_la-MeshManipulation.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='$(SOURCE_DIR)/parallel/MeshManipulation.cc' object='libamdis_la-MeshManipulation.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-MeshManipulation.lo `test -f '$(SOURCE_DIR)/parallel/MeshManipulation.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/MeshManipulation.cc
+
 libamdis_la-ParallelDebug.lo: $(SOURCE_DIR)/parallel/ParallelDebug.cc
 @am__fastdepCXX_TRUE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ParallelDebug.lo -MD -MP -MF $(DEPDIR)/libamdis_la-ParallelDebug.Tpo -c -o libamdis_la-ParallelDebug.lo `test -f '$(SOURCE_DIR)/parallel/ParallelDebug.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ParallelDebug.cc
 @am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/libamdis_la-ParallelDebug.Tpo $(DEPDIR)/libamdis_la-ParallelDebug.Plo
diff --git a/AMDiS/src/DataCollector.cc b/AMDiS/src/DataCollector.cc
index 325be21f..1194027f 100644
--- a/AMDiS/src/DataCollector.cc
+++ b/AMDiS/src/DataCollector.cc
@@ -208,8 +208,8 @@ namespace AMDiS {
     DegreeOfFreedom vertexDOF;
     WorldVector<double> vertexCoords;
 
-    MSG("ELEMENT: %d\n", elInfo->getElement()->getIndex());
-    MSG("DOFs: %d %d %d\n", dof[0][0], dof[1][0], dof[2][0]);
+    // MSG("ELEMENT: %d\n", elInfo->getElement()->getIndex());
+    // MSG("DOFs: %d %d %d\n", dof[0][0], dof[1][0], dof[2][0]);
 
     // create ElementInfo
     ElementInfo elementInfo(dim);
diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc
index 73307c61..7bfdffd4 100644
--- a/AMDiS/src/Element.cc
+++ b/AMDiS/src/Element.cc
@@ -198,107 +198,78 @@ namespace AMDiS {
   /*  should be used only at the end of dof_compress()!!!!!                   */
   /****************************************************************************/
 
-  /* CHANGE_DOFS_1 changes old dofs to NEGATIVE new dofs */
-
-#define CHANGE_DOFS_1(el)						\
-  ldof = el->dof[n0 + i] + nd0;						\
-    for (j = 0; j < nd; j++) {						\
-      if ((k = ldof[j]) >= 0) {						\
-	/* do it only once! (dofs are visited more than once) */	\
-	ldof[j] = - admin->getMesh()->newDOF[k] - 1;			\
-      } }
-
-  /* CHANGE_DOFS_2 changes NEGATIVE new dofs to POSITIVE */
-
-#define CHANGE_DOFS_2(el)						\
-  ldof = el->dof[n0+i] + nd0;						\
-    for (j = 0; j < nd; j++) {						\
-      if ((k = ldof[j]) < 0) {						\
-	/* do it only once! (dofs are visited more than once) */	\
-	ldof[j] = - k - 1;						\
-      } }
-
-  void Element::newDOFFct1(const DOFAdmin* admin)
+  void Element::newDofFct1(const DOFAdmin* admin)
   {
-    int j, k, n0, nd, nd0;
-    DegreeOfFreedom *ldof;  
-    int vertices = mesh->getGeo(VERTEX);
-    int edges = mesh->getGeo(EDGE); 
-    int faces = mesh->getGeo(FACE);
+    int n0, nd, nd0;
 
     if ((nd = admin->getNumberOfDofs(VERTEX)))  {
+      int vertices = mesh->getGeo(VERTEX);
       nd0 = admin->getNumberOfPreDOFs(VERTEX);
       n0 = admin->getMesh()->getNode(VERTEX);
-      for (int i = 0; i < vertices; i++) {
-	CHANGE_DOFS_1(this);
-      }
+      for (int i = 0; i < vertices; i++)
+	changeDofs1(admin, n0, nd0, nd, i);
     }
 
     if (mesh->getDim() > 1) {
       if ((nd = admin->getNumberOfDofs(EDGE)))  {
+	int edges = mesh->getGeo(EDGE); 
 	nd0 = admin->getNumberOfPreDOFs(EDGE);
 	n0 = admin->getMesh()->getNode(EDGE);
-	for (int i = 0; i < edges; i++) {
-	  CHANGE_DOFS_1(this);
-	}
+	for (int i = 0; i < edges; i++)
+	  changeDofs1(admin, n0, nd0, nd, i);
       }
     }
 
     if (mesh->getDim() == 3) {
       if ((nd = admin->getNumberOfDofs(FACE)))  {
+	int faces = mesh->getGeo(FACE);
 	nd0 = admin->getNumberOfPreDOFs(FACE);
 	n0 = admin->getMesh()->getNode(FACE);
-	for (int i = 0; i < faces; i++) {
-	  CHANGE_DOFS_1(this);
-	}
+	for (int i = 0; i < faces; i++)
+	  changeDofs1(admin, n0, nd0, nd, i);
       }
     }
 
     if ((nd = admin->getNumberOfDofs(CENTER)))  {
       nd0 = admin->getNumberOfPreDOFs(CENTER);
-      n0 = admin->getMesh()->getNode(CENTER);
-      int i = 0;          /* only one center */
-      CHANGE_DOFS_1(this);
+      n0 = admin->getMesh()->getNode(CENTER);      
+      changeDofs1(admin, n0, nd0, nd, 0);
     }
   }
 
 
-  void Element::newDOFFct2(const DOFAdmin* admin)
+  void Element::newDofFct2(const DOFAdmin* admin)
   {
-    int j, k, n0, nd0;
-    DegreeOfFreedom  *ldof;
-    int vertices = mesh->getGeo(VERTEX);
-    int edges = mesh->getGeo(EDGE); 
-    int faces = mesh->getGeo(FACE);
+    int n0, nd0;
 
     int nd = admin->getNumberOfDofs(VERTEX);
     if (nd) {
+      int vertices = mesh->getGeo(VERTEX);
       nd0 = admin->getNumberOfPreDOFs(VERTEX);
       n0 = admin->getMesh()->getNode(VERTEX);
-      for (int i = 0; i < vertices; i++) {
-	CHANGE_DOFS_2(this);
-      }
+      for (int i = 0; i < vertices; i++)
+	changeDofs2(n0, nd0, nd, i);
     }
 
     if (mesh->getDim() > 1) {
       nd = admin->getNumberOfDofs(EDGE);
       if (nd) {
+	int edges = mesh->getGeo(EDGE); 
 	nd0 = admin->getNumberOfPreDOFs(EDGE);
 	n0 = admin->getMesh()->getNode(EDGE);
-	for (int i = 0; i < edges; i++) {
-	  CHANGE_DOFS_2(this);
-	}
+	for (int i = 0; i < edges; i++)
+	  changeDofs2(n0, nd0, nd, i);
       }
     }
 
     if (mesh->getDim() == 3) {
       nd = admin->getNumberOfDofs(FACE);
       if (nd) {
+	int faces = mesh->getGeo(FACE);
 	nd0 = admin->getNumberOfPreDOFs(FACE);
 	n0 = admin->getMesh()->getNode(FACE);
-	for (int i = 0; i < faces; i++) {
-	  CHANGE_DOFS_2(this);
-	}
+	for (int i = 0; i < faces; i++)
+	  changeDofs2(n0, nd0, nd, i);
       }
     }
 
@@ -307,13 +278,31 @@ namespace AMDiS {
       nd0 = admin->getNumberOfPreDOFs(CENTER);
       n0 = admin->getMesh()->getNode(CENTER);
       // only one center
-      int i = 0;   
-      CHANGE_DOFS_2(this);
+      changeDofs2(n0, nd0, nd, 0);	  
     }
   }
 
-#undef CHANGE_DOFS_1
-#undef CHANGE_DOFS_2
+
+  void Element::changeDofs1(const DOFAdmin* admin, int n0, int nd0, int nd, int pos)
+  {
+    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
+    for (int j = 0; j < nd; j++) {
+      int k = ldof[j];
+      if (k >= 0)
+	ldof[j] = -admin->getMesh()->newDOF[k] - 1;
+    }
+  }
+  
+  
+  void Element::changeDofs2(int n0, int nd0, int nd, int pos)
+  {
+    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
+    for (int j = 0; j < nd; j++) {
+      int k = ldof[j];
+      if (k < 0)
+	ldof[j] = -k - 1;      
+    }
+  }
 
 
   /****************************************************************************/
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index 8707c810..9c70dd2c 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -22,12 +22,12 @@
 #ifndef AMDIS_ELEMENT_H
 #define AMDIS_ELEMENT_H
 
+#include "AMDiS_fwd.h"
 #include "Global.h"
 #include "RefinementManager.h"
 #include "Serializable.h"
 #include "ElementData.h"
 #include "LeafData.h"
-#include "AMDiS_fwd.h"
 
 namespace AMDiS {
 
@@ -521,10 +521,16 @@ namespace AMDiS {
     }
 
     /// Used by friend class Mesh while dofCompress
-    void newDOFFct1(const DOFAdmin*);
+    void newDofFct1(const DOFAdmin*);
 
     /// Used by friend class Mesh while dofCompress
-    void newDOFFct2(const DOFAdmin*);
+    void newDofFct2(const DOFAdmin*);
+
+    /// Changes old dofs to negative new dofs
+    void changeDofs1(const DOFAdmin* admin, int n0, int nd0, int nd, int pos);
+
+    /// Changes negative new dofs to positive
+    void changeDofs2(int n0, int nd0, int nd, int pos);
 
   protected:
     /** \brief
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 170237c6..10001b80 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -472,13 +472,13 @@ namespace AMDiS {
       TraverseStack stack;
       ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
       while (elInfo) {
-	elInfo->getElement()->newDOFFct1(compressAdmin);
+	elInfo->getElement()->newDofFct1(compressAdmin);
 	elInfo = stack.traverseNext(elInfo);
       }
 
       elInfo = stack.traverseFirst(this, -1, fill_flag);
       while (elInfo) {
-	elInfo->getElement()->newDOFFct2(compressAdmin);
+	elInfo->getElement()->newDofFct2(compressAdmin);
 	elInfo = stack.traverseNext(elInfo);
       }
 
@@ -613,6 +613,7 @@ namespace AMDiS {
     }
 
     delete [] dof;
+    dof = NULL;
   }
 
 
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index d7071db7..639b5f60 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -880,8 +880,9 @@ namespace AMDiS {
     friend class MacroWriter;
     friend class MacroElement;
     friend class Element;
-    friend void Element::newDOFFct1(const DOFAdmin*);
-    friend void Element::newDOFFct2(const DOFAdmin*);
+    friend void Element::newDofFct1(const DOFAdmin*);
+    friend void Element::newDofFct2(const DOFAdmin*);
+    friend void Element::changeDofs1(const DOFAdmin*, int, int, int, int);
   };
 
 }
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 8e8c8efc..f1804c2e 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -7,6 +7,7 @@
 #include <boost/filesystem.hpp>
 
 #include "parallel/MeshDistributor.h"
+#include "parallel/MeshManipulation.h"
 #include "parallel/ParallelDebug.h"
 #include "parallel/StdMpi.h"
 #include "parallel/ParMetisPartitioner.h"
@@ -54,6 +55,7 @@ namespace AMDiS {
       rstart(0),
       deserialized(false),
       writeSerializationFile(false),
+      repartitioningAllowed(false),
       lastMeshChangeIndex(0),
       macroElementStructureConsisten(false)
   {
@@ -62,6 +64,10 @@ namespace AMDiS {
     mpiRank = MPI::COMM_WORLD.Get_rank();
     mpiSize = MPI::COMM_WORLD.Get_size();
     mpiComm = MPI::COMM_WORLD;
+
+    int tmp = 0;
+    GET_PARAMETER(0, name + "->repartitioning", "%d", &tmp);
+    repartitioningAllowed = (tmp > 0);
   }
 
 
@@ -566,7 +572,8 @@ namespace AMDiS {
 
     // === The mesh has changed, so check if it is required to repartition the mesh. ===
 
-    repartitionMesh();
+    if (repartitioningAllowed)
+      repartitionMesh();
   }
 
   
@@ -981,8 +988,6 @@ namespace AMDiS {
     if (repartitioning == 0)
       return;
 
-    return;
-
     DOFVector<double> tmpa(feSpace, "tmp");
     tmpa.set(mpiRank);
     VtkWriter::writeFile(tmpa, "before-repartition.vtu");
@@ -1040,8 +1045,6 @@ namespace AMDiS {
     for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
 	 it != newMacroEl.end(); ++it) {
       MacroElement *mel = *it;
-      //      MSG("NEW MACRO EL: %d\n", mel->getIndex());
-
       for (int i = 0; i < mesh->getGeo(NEIGH); i++)
 	mel->setNeighbour(i, NULL);
 
@@ -1059,7 +1062,6 @@ namespace AMDiS {
 	 it != partitioner->getSendElements().end(); ++it) {
       for (std::vector<int>::iterator elIt = it->second.begin();
 	   elIt != it->second.end(); ++elIt) {
-	//	MSG("MAKE EL-STRUC FOR EL %d TO RANK %d\n", *elIt, it->first);
 	MeshStructure elCode;
 	elCode.init(mesh, *elIt);
 	sendCodes[it->first].push_back(elCode);
@@ -1083,16 +1085,10 @@ namespace AMDiS {
 	("Should not happen!\n");
 
       for (std::vector<int>::iterator elIt = it->second.begin();
-	   elIt != it->second.end(); ++elIt) {
-	//	MSG("RECV EL-STRUC FOR EL %d FROM RANK %d\n", *elIt, it->first);
-
+	   elIt != it->second.end(); ++elIt)
 	recvCodes[i++].fitMeshToStructure(mesh, refineManager, false, *elIt);
-      }
     }
 
-    MSG("DAS IST FERTIG!\n");
-
-    
 
     // === Set correct neighbour information on macro elements. ===
 
@@ -1120,120 +1116,24 @@ namespace AMDiS {
 	   elIt != it->second.end(); ++elIt) {
 	TEST_EXIT_DBG(elIndexMap.count(*elIt) == 1)
 	  ("Could not find macro element %d\n", *elIt);
-	//	MSG("DELETE MACRO %d\n", (*elIt));
+
 	deleteMacroElements.insert(elIndexMap[*elIt]);
       }
     }
 
-    //    MSG("REMOVE MACRO ELEMENT!\n");
-    
     mesh->removeMacroElements(deleteMacroElements, feSpace);
 
 
     // === Remove double DOFs. ===
 
-    std::map<int, MacroElement*> leafInMacroEl;
-    
-    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
-    while (elInfo) {
-      leafInMacroEl[elInfo->getElement()->getIndex()] = elInfo->getMacroElement();
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    std::map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
-
-    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
-    while (elInfo) {
-      Element *el = elInfo->getElement();
-
-      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
-	Element *neigh = elInfo->getNeighbour(i);
-
-	if (!neigh)
-	  continue;      
-
-	if (leafInMacroEl[el->getIndex()] != 
-	    leafInMacroEl[neigh->getIndex()] && 
-	    newMacroEl.count(leafInMacroEl[el->getIndex()]) == 0 &&
-	    newMacroEl.count(leafInMacroEl[neigh->getIndex()]) == 1) {
-
-	  //	  MSG("EL: %d   NEIGH:  %d\n", el->getIndex(), neigh->getIndex());
-
-	  int elEdgeIndex = i;
-	  int neighEdgeIndex = elInfo->getOppVertex(i);
-
-	  DofEdge elEdge = el->getEdge(elEdgeIndex);
-	  DofEdge neighEdge = neigh->getEdge(neighEdgeIndex);	  
-
-	  if (elEdge != neighEdge) {
-	    int elEdgeIndex0 = el->getVertexOfEdge(elEdgeIndex, 0);
-	    int elEdgeIndex1 = el->getVertexOfEdge(elEdgeIndex, 1);
-
-	    int neighEdgeIndex0 = neigh->getVertexOfEdge(neighEdgeIndex, 0);
-	    int neighEdgeIndex1 = neigh->getVertexOfEdge(neighEdgeIndex, 1);
-
-	    DegreeOfFreedom elEdgeDof0 = el->getDof(elEdgeIndex0, 0);
-	    DegreeOfFreedom elEdgeDof1 = el->getDof(elEdgeIndex1, 0);
-
-	    DegreeOfFreedom neighEdgeDof0 = neigh->getDof(neighEdgeIndex0, 0);
-	    DegreeOfFreedom neighEdgeDof1 = neigh->getDof(neighEdgeIndex1, 0);
-
-	    if ((elEdgeDof0 < elEdgeDof1 && neighEdgeDof0 > neighEdgeDof1) ||
-		(elEdgeDof0 > elEdgeDof1 && neighEdgeDof0 < neighEdgeDof1)) {
-	      std::swap(neighEdgeIndex0, neighEdgeIndex1);
-	      std::swap(neighEdgeDof0, neighEdgeDof1);
-	    }
-
-	    if (neighEdgeDof0 != elEdgeDof0) {
-	      mapDelDofs[neigh->getDof(neighEdgeIndex0)] = el->getDof(elEdgeIndex0);
-	      //	      MSG("REPLACE DOF %d by %d\n", neigh->getDof(neighEdgeIndex0, 0), el->getDof(elEdgeIndex0, 0));
-	    }
-
-	    if (neighEdgeDof1 != elEdgeDof1) {
-	      mapDelDofs[neigh->getDof(neighEdgeIndex1)] = el->getDof(elEdgeIndex1);
-	      //	      MSG("REPLACE DOF %d by %d\n", neigh->getDof(neighEdgeIndex1, 0), el->getDof(elEdgeIndex1, 0));
-	    }
-	  }
-	}
-      }
-
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-
-    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
-    while (elInfo) {
-      for (int i = 0; i < mesh->getGeo(VERTEX); i++)
-	if (mapDelDofs.count(elInfo->getElement()->getDof(i)) == 1) {
-	  // MSG("SET DOF %d TO DOF %d in EL %d\n",
-	  //     elInfo->getElement()->getDof(i, 0),
-	  //     *(mapDelDofs[elInfo->getElement()->getDof(i)]),
-	  //     elInfo->getElement()->getIndex());
-
-	  elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)]));
-	}
-
-      elInfo = stack.traverseNext(elInfo);
-    }
-
-    for (std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
-	 it != mapDelDofs.end(); ++it) {
-      //      MSG("DOF TO DEL: %d\n", (*(it->first)));
-
-      mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), VERTEX);
-    }
-
+    MeshManipulation meshManipulation(mesh);
+    meshManipulation.deleteDoubleDofs(newMacroEl);
     mesh->dofCompress();
 
-    MSG("WRITE OUTPUT\n");
-
     DOFVector<double> tmp(feSpace, "tmp");
     tmp.set(mpiRank);
     VtkWriter::writeFile(tmp, "after-repartition.vtu");
 
-
-    MSG("DONE!\n");
-
     MSG("USED-SIZE B: %d\n", mesh->getDofAdmin(0).getUsedDofs());
     exit(0);
 
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index e3007237..6c3ccec9 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -107,13 +107,6 @@ namespace AMDiS {
     /// Set for each element on the partitioning level the number of leaf elements.
     void setInitialElementWeights();
 
-    /** \brief
-     * Checks if is required to repartition the mesh. If this is the case, a new
-     * partition will be created and the mesh will be redistributed between the
-     * ranks.
-     */
-    void repartitionMesh();
-
     inline virtual std::string getName() 
     { 
       return name; 
@@ -314,7 +307,14 @@ namespace AMDiS {
      *                       should be checked.
      */
     bool checkAndAdaptBoundary(RankToBoundMap &allBound);
-    
+  
+    /** \brief
+     * Checks if is required to repartition the mesh. If this is the case, a new
+     * partition will be created and the mesh will be redistributed between the
+     * ranks.
+     */
+    void repartitionMesh();
+
     /** \brief
      * This functions create a Paraview file with the macro mesh where the elements
      * are colored by the partition they are part of. This function can be used for
@@ -576,6 +576,9 @@ namespace AMDiS {
     /// Denotes whether there exists a filewriter for this object.
     bool writeSerializationFile;
 
+    /// If true, it is possible to repartition the mesh during computations.
+    bool repartitioningAllowed;
+
     /** \brief
      * Stores the mesh change index. This is used to recognize changes in the mesh 
      * structure (e.g. through refinement or coarsening managers).
diff --git a/AMDiS/src/parallel/MeshManipulation.cc b/AMDiS/src/parallel/MeshManipulation.cc
new file mode 100644
index 00000000..4c04faac
--- /dev/null
+++ b/AMDiS/src/parallel/MeshManipulation.cc
@@ -0,0 +1,150 @@
+#include "parallel/MeshManipulation.h"
+#include "Mesh.h"
+#include "Traverse.h"
+
+namespace AMDiS {
+
+  void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl)
+  {
+    std::map<int, MacroElement*> leafInMacroEl;
+    
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    while (elInfo) {
+      leafInMacroEl[elInfo->getElement()->getIndex()] = elInfo->getMacroElement();
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    deleteDoubleDofs(leafInMacroEl, newMacroEl, 0);
+    deleteDoubleDofs(leafInMacroEl, newMacroEl, 1);
+  }
+
+
+  void MeshManipulation::deleteDoubleDofs(std::map<int, MacroElement*>& leafInMacroEl,
+					  std::set<MacroElement*>& newMacroEl,
+					  int mode)
+  {
+    FUNCNAME("MeshManipulation::deleteDoubleDofs()");
+
+    std::map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
+    std::map<const DegreeOfFreedom*, int> mapDofsMacro;
+
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
+    while (elInfo) {
+      Element *el = elInfo->getElement();
+
+      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
+	Element *neigh = elInfo->getNeighbour(i);
+
+	if (!neigh)
+	  continue;      
+
+	int elMacroIndex = leafInMacroEl[el->getIndex()]->getIndex();
+	int neighMacroIndex = leafInMacroEl[neigh->getIndex()]->getIndex();
+
+	if (elMacroIndex == neighMacroIndex)
+	  continue;
+
+	if ((mode == 0 && 
+	     newMacroEl.count(leafInMacroEl[el->getIndex()]) == 1 &&
+	     newMacroEl.count(leafInMacroEl[neigh->getIndex()]) == 1 &&
+	     elMacroIndex > neighMacroIndex) ||
+	    (mode == 1 &&
+	     newMacroEl.count(leafInMacroEl[el->getIndex()]) == 0 &&
+	     newMacroEl.count(leafInMacroEl[neigh->getIndex()]) == 1)) {
+
+	  //	  MSG("EL: %d   NEIGH:  %d\n", el->getIndex(), neigh->getIndex());
+
+	  if (el->getEdge(i) != neigh->getEdge(elInfo->getOppVertex(i))) {
+	    std::vector<int> elEdgeIndex(2);
+	    elEdgeIndex[0] = el->getVertexOfEdge(i, 0);
+	    elEdgeIndex[1] = el->getVertexOfEdge(i, 1);
+
+	    std::vector<int> neighEdgeIndex(2);
+	    neighEdgeIndex[0] = neigh->getVertexOfEdge(elInfo->getOppVertex(i), 0);
+	    neighEdgeIndex[1] = neigh->getVertexOfEdge(elInfo->getOppVertex(i), 1);
+
+	    std::vector<DegreeOfFreedom> elEdgeDof(2);
+	    elEdgeDof[0] = el->getDof(elEdgeIndex[0], 0);
+	    elEdgeDof[1] = el->getDof(elEdgeIndex[1], 0);
+
+	    std::vector<DegreeOfFreedom> neighEdgeDof(2);
+	    neighEdgeDof[0] = neigh->getDof(neighEdgeIndex[0], 0);
+	    neighEdgeDof[1] = neigh->getDof(neighEdgeIndex[1], 0);
+
+	    //	    MSG("TEST DOFs %d  %d    <->    %d  %d\n", elEdgeDof0, elEdgeDof1, neighEdgeDof0, neighEdgeDof1);
+
+	    if ((elEdgeDof[0] < elEdgeDof[1] && neighEdgeDof[0] > neighEdgeDof[1]) ||
+		(elEdgeDof[0] > elEdgeDof[1] && neighEdgeDof[0] < neighEdgeDof[1])) {
+	      std::swap(neighEdgeIndex[0], neighEdgeIndex[1]);
+	      std::swap(neighEdgeDof[0], neighEdgeDof[1]);
+	    }
+
+	    for (int j = 0; j < 2; j++) {
+
+	      if (neighEdgeDof[j] != elEdgeDof[j]) {
+		const DegreeOfFreedom *delDof = neigh->getDof(neighEdgeIndex[j]);
+		const DegreeOfFreedom *replaceDof = el->getDof(elEdgeIndex[j]);
+
+		if (mapDelDofs.count(delDof) == 1 && mapDelDofs[delDof] != replaceDof) {
+		  if (mapDofsMacro[mapDelDofs[delDof]] > elMacroIndex) {
+		    mapDelDofs[replaceDof] = mapDelDofs[delDof];
+		  } else {
+		    mapDelDofs[mapDelDofs[delDof]] = replaceDof;
+		    mapDelDofs[delDof] = replaceDof;
+		    mapDofsMacro[replaceDof] = elMacroIndex;
+		  }
+		} else {
+		  mapDelDofs[delDof] = replaceDof;
+		  mapDofsMacro[replaceDof] = elMacroIndex;
+		}
+	      }
+	    }	      
+	  }
+	}
+      }
+
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    bool changed = false;
+    do {
+      changed = false;
+      for (std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
+	   it != mapDelDofs.end(); ++it) {
+	std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second);
+	if (findIt != mapDelDofs.end()) {
+	  it->second = findIt->second;   
+	  changed = true;
+	}
+      } 
+    } while (changed);
+
+
+    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+    while (elInfo) {
+      for (int i = 0; i < mesh->getGeo(VERTEX); i++)
+	if (mapDelDofs.count(elInfo->getElement()->getDof(i)) == 1) {
+	  // MSG("SET DOF %d TO DOF %d in EL %d\n",
+	  //      elInfo->getElement()->getDof(i, 0),
+	  //      *(mapDelDofs[elInfo->getElement()->getDof(i)]),
+	  //      elInfo->getElement()->getIndex());
+
+	  elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)]));
+	}      
+
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    for (std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
+	 it != mapDelDofs.end(); ++it) {
+      //      MSG("DOF TO DEL: %d\n", (*(it->first)));
+
+      mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), VERTEX);
+    }
+
+    MSG("FINISHED DELLING!\n");	
+  }
+
+}
diff --git a/AMDiS/src/parallel/MeshManipulation.h b/AMDiS/src/parallel/MeshManipulation.h
new file mode 100644
index 00000000..8fb4741d
--- /dev/null
+++ b/AMDiS/src/parallel/MeshManipulation.h
@@ -0,0 +1,52 @@
+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ============================================================================
+// ==                                                                        ==
+// ==  TU Dresden                                                            ==
+// ==                                                                        ==
+// ==  Institut für Wissenschaftliches Rechnen                               ==
+// ==  Zellescher Weg 12-14                                                  ==
+// ==  01069 Dresden                                                         ==
+// ==  germany                                                               ==
+// ==                                                                        ==
+// ============================================================================
+// ==                                                                        ==
+// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
+// ==                                                                        ==
+// ============================================================================
+
+/** \file MeshManipulation.h */
+
+#ifndef AMDIS_MESHMANIPULATION_H
+#define AMDIS_MESHMANIPULATION_H
+
+#include <set>
+#include <map>
+
+#include "AMDiS_fwd.h"
+
+namespace AMDiS {
+
+  class MeshManipulation {
+  public:
+    MeshManipulation(Mesh *m)
+      : mesh(m)
+    {}
+
+    void deleteDoubleDofs(std::set<MacroElement*>& newMacroEl);
+
+  private:
+    void deleteDoubleDofs(std::map<int, MacroElement*>& leafInMacroEl,
+			  std::set<MacroElement*>& newMacroEl,
+			  int mode);
+
+  private:
+    Mesh *mesh;
+  };
+
+}
+
+#endif
+
-- 
GitLab