Fixed problem for repartitioning with periodic boundary conditions.

AMDiS/src/parallel/ElementObjectData.h

@@ -90,13 +90,13 @@ namespace AMDiS {
 
   /** \brief
    * This class is a database of element objects. An element object is either a
-   * vertex, edge or the face of a specific element. This database is used to store
-   * all objects of all elements of a mesh. The information is stored in a way that
-   * makes it possible to identify all elements, which have a given vertex, edge or
-   * face in common. If is is known which element is owned by which rank in parallel
-   * computations, it is thus possible to get all interior boundaries on object 
-   * level. This is required, because two elements may share a common vertex without
-   * beging neighbours in the definition of AMDiS.
+   * vertex, edge or the face of a specific element. This database is used to
+   * store all objects of all elements of a mesh. The information is stored in a
+   * way that makes it possible to identify all elements, which have a given
+   * vertex, edge or face in common. If is is known which element is owned by 
+   * which rank in parallel computations, it is thus possible to get all interior
+   * boundaries on object level. This is required, because two elements may share
+   * a common vertex without beging neighbours in the definition of AMDiS.
    */
   class ElementObjects {
   public:
@@ -114,9 +114,9 @@ namespace AMDiS {
 
 
     /** \brief
-     * Adds an element to the object database. If the element is part of a periodic
-     * boundary, all information about subobjects of the element on this boundary
-     * are collected.
+     * Adds an element to the object database. If the element is part of a
+     * periodic boundary, all information about subobjects of the element on
+     * this boundary are collected.
      *
      * \param[in]  elInfo    ElInfo object of the element. 
      */
@@ -136,15 +136,16 @@ namespace AMDiS {
 
 
     /** \brief
-     * Create for a filled object database the membership information for all element
-     * objects. An object is owned by a rank, if the rank has the heighest rank
-     * number of all ranks where the object is part of.
+     * Create for a filled object database the membership information for all
+     * element objects. An object is owned by a rank, if the rank has the
+     * heighest rank number of all ranks where the object is part of.
      *
-     * \param[in]  macroElementRankMap   Maps to each macro element of the mesh the
-     *                                   rank that owns this macro element.
+     * \param[in]  macroElementRankMap   Maps to each macro element of the mesh
+     *                                   the rank that owns this macro element.
      */
     void createRankData(map<int, int>& macroElementRankMap);
 
+
     /** \brief
      * Creates on all boundaries the reverse mode flag.
      *
@@ -162,12 +163,12 @@ namespace AMDiS {
 
 
     /** \brief
-     * Iterates over all elements for one geometrical index, i.e., over all vertices,
-     * edges or faces in the mesh. The function returns true, if the result is valid.
-     * Otherwise the iterator is at the end position.
+     * Iterates over all elements for one geometrical index, i.e., over all
+     * vertices, edges or faces in the mesh. The function returns true, if the
+     * result is valid. Otherwise the iterator is at the end position.
      *
-     * \param[in]  pos   Must be either VERTEX, EDGE or FACE and defines the elements
-     *                   that should be traversed.
+     * \param[in]  pos   Must be either VERTEX, EDGE or FACE and defines the
+     *                   elements that should be traversed.
      */
     bool iterate(GeoIndex pos)
     {
@@ -419,7 +420,8 @@ namespace AMDiS {
       return periodicFaces;
     }
 
-    inline bool getEdgeReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+    inline bool getEdgeReverseMode(ElementObjectData &obj0, 
+				   ElementObjectData &obj1)
     {
       if (mesh->getDim() == 2)
 	return true;
@@ -430,7 +432,8 @@ namespace AMDiS {
       return edgeReverseMode[make_pair(obj0, obj1)];
     }
 
-    inline bool getFaceReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+    inline bool getFaceReverseMode(ElementObjectData &obj0, 
+				   ElementObjectData &obj1)
     {
       TEST_EXIT_DBG(faceReverseMode.count(make_pair(obj0, obj1)))
 	("Should not happen!\n");
@@ -540,12 +543,12 @@ namespace AMDiS {
     /// objects that have this vertex DOF in common.
     map<DegreeOfFreedom, map<int, ElementObjectData> > vertexInRank;
 
-    /// Defines to each edge a map that maps to each rank number the element objects
-    /// that have this edge in common.
+    /// Defines to each edge a map that maps to each rank number the element 
+    /// objects that have this edge in common.
     map<DofEdge, map<int, ElementObjectData> > edgeInRank;
 
-    /// Defines to each face a map that maps to each rank number the element objects
-    /// that have this face in common.
+    /// Defines to each face a map that maps to each rank number the element 
+    /// objects that have this face in common.
     map<DofFace, map<int, ElementObjectData> > faceInRank;
 
 

AMDiS/src/parallel/MeshDistributor.cc

@@ -94,16 +94,10 @@ namespace AMDiS {
     mpiSize = MPI::COMM_WORLD.Get_size();
     mpiComm = MPI::COMM_WORLD;
 
-    int tmp = 0;
-    Parameters::get(name + "->repartitioning", tmp);
-    repartitioningAllowed = (tmp > 0);
-
+    Parameters::get(name + "->repartitioning", repartitioningAllowed);
     Parameters::get(name + "->debug output dir", debugOutputDir);
     Parameters::get(name + "->repartition ith change", repartitionIthChange);
-
-    tmp = 0;
-    Parameters::get(name + "->log main rank", tmp);
-    Msg::outputMainRank = (tmp > 0);
+    Parameters::get(name + "->log main rank", Msg::outputMainRank);
 
     string partStr = "parmetis";
     Parameters::get(name + "->partitioner", partStr);
@@ -125,7 +119,7 @@ namespace AMDiS {
     if (partStr == "simple")
       partitioner = new SimplePartitioner(&mpiComm);
 
-    tmp = 0;
+    int tmp = 0;
     Parameters::get(name + "->box partitioning", tmp);
     partitioner->setBoxPartitioning(static_cast<bool>(tmp));
 
@@ -874,7 +868,7 @@ namespace AMDiS {
       repartitionMesh();
       nMeshChangesAfterLastRepartitioning = 0;
     } else {
-      MSG_DBG("Repartitioning not tried because tryRepartitioning = %d   reparttitioningAllowed = %d   nMeshChange =%d    repartitionIthChange = %d\n",
+      MSG_DBG("Repartitioning not tried because tryRepartitioning = %d   repartitioningAllowed = %d   nMeshChange = %d   repartitionIthChange = %d\n",
 	      tryRepartition, repartitioningAllowed, 
 	      nMeshChangesAfterLastRepartitioning, repartitionIthChange);
     }
@@ -896,15 +890,18 @@ namespace AMDiS {
     if (mpiRank == 0) {
       int nOverallDofs = 0;
       int maxDofs = numeric_limits<int>::min();
+      int minDofs = numeric_limits<int>::max();
       for (int i = 0; i < mpiSize; i++) {
 	nOverallDofs += nDofsInRank[i];
 	maxDofs = std::max(maxDofs, nDofsInRank[i]);
+	minDofs = std::min(minDofs, nDofsInRank[i]);
       }      
-      int avrgDofs = nOverallDofs / mpiSize;
-      double imbalance = 
-	(static_cast<double>(maxDofs - avrgDofs) /  avrgDofs) * 100.0;
+//       int avrgDofs = nOverallDofs / mpiSize;
+//       double imbalance0 = 
+// 	(static_cast<double>(maxDofs - avrgDofs) /  avrgDofs) * 100.0;
+      double imbalance1 = (static_cast<double>(maxDofs) / minDofs - 1.0) * 100.0;
 
-      MSG("Imbalancing factor: %.1f\%\n", imbalance);
+      MSG("Imbalancing factor: %.1f\%\n", imbalance1);
     }
   }
 
@@ -1157,12 +1154,18 @@ namespace AMDiS {
     // === Run mesh partitioner to calculate a new mesh partitioning.  ===
 
     partitioner->setLocalGlobalDofMap(&(dofFeData[feSpaces[0]].mapDofToGlobal));
-    bool partitioningSucceed = partitioner->partition(elemWeights, ADAPTIVE_REPART);
+    bool partitioningSucceed = 
+      partitioner->partition(elemWeights, ADAPTIVE_REPART);
     if (!partitioningSucceed) {
       MSG("Mesh partitioner created empty partition!\n");
       return;
     }
 
+    // In the case the partitioner does not create a new mesh partition, return
+    // without and changes.
+    if (!partitioner->meshChanged())
+      return;
+
 
     // === Create map that maps macro element indices to pointers to the ===
     // === macro elements.                                               ===
@@ -2023,13 +2026,17 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::createPeriodicMap()");
 
-    if (periodicBoundary.boundary.size() == 0)
-      return;
-
     // Clear all periodic DOF mappings calculated before. We do it from scratch.
     periodicDofs.clear();
     periodicMap.clear();
 
+    // If there are no periodic boundaries, return. Note that periodicDofs and
+    // periodicMap must be still cleared before: if we do repartitioning and
+    // there were periodic boundaries in subdomain before and after repartitioning
+    // there are no more periodic boundaries.
+    if (periodicBoundary.boundary.size() == 0)
+      return;
+
     for (unsigned int i = 0; i < feSpaces.size(); i++)
       createPeriodicMap(feSpaces[i]);
   }

AMDiS/src/parallel/MeshPartitioner.h

@@ -29,13 +29,13 @@
 
 #include "AMDiS_fwd.h"
 #include "Mesh.h"
+#include "parallel/MpiHelper.h"
 
-namespace AMDiS {
 
+namespace AMDiS {
 
   using namespace std;
 
-
   enum PartitionMode {
     INITIAL = 0,          // initial partitioning of a unpartitioned mesh 
     ADAPTIVE_REPART = 1,  // repartitioning of a adaptively refined mesh
@@ -130,6 +130,16 @@ namespace AMDiS {
     {
       return sendElements;
     }
+    
+    /// After mesh repartition this function returns true if the mesh must be 
+    /// redistributed on at least one rank.
+    bool meshChanged()
+    {
+      int nChanges = recvElements.size() + sendElements.size();
+      mpi::globalAdd(nChanges);
+
+      return static_cast<bool>(nChanges);
+    }
 
   protected:
     /// Pointer to the MPI communicator the mesh partitioner should make use of.
@@ -170,6 +180,8 @@ namespace AMDiS {
     /// partitiong mode) the rank number the element belongs to.
     map<int, int> partitionMap;
 
+    /// After mesh repartitioning these maps stores which elements are communicated
+    /// from this rank to other ranks.
     map<int, vector<int> > recvElements, sendElements;
   };
 }

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -398,6 +398,8 @@ namespace AMDiS {
     MSG("nRankDuals = %d   nOverallDuals = %d\n",
 	dualDofMap[feSpace].nRankDofs, 
 	dualDofMap[feSpace].nOverallDofs);
+
+    nLocalDuals = dualDofMap[feSpace].size();
   }
 
   
@@ -409,26 +411,26 @@ namespace AMDiS {
     // === appropriate number of Lagrange constraints.                      ===
 
     const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-    dofFirstLagrange.addFeSpace(feSpace);
+    lagrangeMap.addFeSpace(feSpace);
 
     int nRankLagrange = 0;
     DofMapping& dualMap = dualDofMap[feSpace].getMap();
     for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
       if (meshDistributor->getIsRankDof(feSpace, it->first)) {
-	dofFirstLagrange[feSpace].insert(it->first, nRankLagrange);
+	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
 	int degree = boundaryDofRanks[it->first].size();
 	nRankLagrange += (degree * (degree - 1)) / 2;
       }
     }
-    dofFirstLagrange[feSpace].nRankDofs = nRankLagrange;
-    dofFirstLagrange[feSpace].update();
+    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
+    lagrangeMap[feSpace].update();
     
     MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
-	dofFirstLagrange[feSpace].nRankDofs, 
-	dofFirstLagrange[feSpace].nOverallDofs);
+	lagrangeMap[feSpace].nRankDofs, 
+	lagrangeMap[feSpace].nOverallDofs);
 
 
-    // === Communicate dofFirstLagrange to all other ranks. ===
+    // === Communicate lagrangeMap to all other ranks. ===
 
     StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
 
@@ -436,9 +438,7 @@ namespace AMDiS {
 	 !it.end(); it.nextRank()) {
       for (; !it.endDofIter(); it.nextDof())
 	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
-	  TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it.getDofIndex()))
-	    ("Should not happen!\n");
-	  DegreeOfFreedom d = dofFirstLagrange[feSpace][it.getDofIndex()];
+	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
 	  stdMpi.getSendData(it.getRank()).push_back(d);
 	}
     }
@@ -466,7 +466,7 @@ namespace AMDiS {
       for (; !it.endDofIter(); it.nextDof()) {
 	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
 	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
-	  dofFirstLagrange[feSpace].insert(it.getDofIndex(), d); 
+	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
 	}
       }
     }
@@ -475,36 +475,29 @@ namespace AMDiS {
 
   void PetscSolverFeti::createIndexB()
   {
-    FUNCNAME("PetscSolverFeti::createIndeB()");
+    FUNCNAME("PetscSolverFeti::createIndexB()");
 
     const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-    localIndexB.clear();
+    localDofMap.addFeSpace(feSpace);
     DOFAdmin* admin = feSpace->getAdmin();
 
     // === To ensure that all interior node on each rank are listen first in ===
     // === the global index of all B nodes, insert all interior nodes first, ===
     // === without defining a correct index.                                 ===
-    
-    for (int i = 0; i < admin->getUsedSize(); i++)
+
+    nLocalInterior = 0;    
+    for (int i = 0; i < admin->getUsedSize(); i++) {
       if (admin->isDofFree(i) == false && 
 	  primalDofMap[feSpace].isSet(i) == false && 
-	  dualDofMap[feSpace].isSet(i) == false)
-	localIndexB[i] = -1;
-
-
-    // === Get correct index for all interior nodes. ===
-
-    nLocalInterior = 0;
-    for (DofMapping::iterator it = localIndexB.begin(); 
-	 it != localIndexB.end(); ++it) {
-      it->second = nLocalInterior;
-      nLocalInterior++;
+	  dualDofMap[feSpace].isSet(i) == false) {
+	localDofMap[feSpace].insertRankDof(i);
+	nLocalInterior++;
+      }
     }
-    nLocalDuals = dualDofMap[feSpace].size();
 
     TEST_EXIT_DBG(nLocalInterior + 
 		  primalDofMap[feSpace].size() + 
-		  dualDofMap[feSpace].size() == 
+		  nLocalDuals == 
 		  static_cast<unsigned int>(admin->getUsedDofs()))
       ("Should not happen!\n");
 
@@ -513,7 +506,8 @@ namespace AMDiS {
 
     for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
 	 it != dualDofMap[feSpace].getMap().end(); ++it)
-      localIndexB[it->first] = it->second - rStartB;
+      localDofMap[feSpace].insert(it->first, it->second - rStartB);
+      //      localDofMap[feSpace].insertRankDof(it->first);
   }
 
 
@@ -526,9 +520,9 @@ namespace AMDiS {
     // === Create distributed matrix for Lagrange constraints. ===
 
     MatCreateMPIAIJ(PETSC_COMM_WORLD,
-		    dofFirstLagrange[feSpace].nRankDofs * nComponents, 
+		    lagrangeMap[feSpace].nRankDofs * nComponents, 
 		    nRankB * nComponents,
-		    dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
+		    lagrangeMap[feSpace].nOverallDofs * nComponents, 
 		    nOverallB * nComponents,
 		    2, PETSC_NULL, 2, PETSC_NULL,
 		    &mat_lagrange);
@@ -542,13 +536,11 @@ namespace AMDiS {
 
     DofMapping &dualMap = dualDofMap[feSpace].getMap();
     for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
-      TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it->first))
-	("Should not happen!\n");
       TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
 	("Should not happen!\n");
 
       // Global index of the first Lagrange constriant for this node.
-      int index = dofFirstLagrange[feSpace][it->first];
+      int index = lagrangeMap[feSpace][it->first];
       // Copy set of all ranks that contain this dual node.
       vector<int> W(boundaryDofRanks[it->first].begin(), 
 		    boundaryDofRanks[it->first].end());
@@ -742,18 +734,18 @@ namespace AMDiS {
 		 nRankB * nComponents, nOverallB * nComponents,
 		 &(fetiData.tmp_vec_b));
     VecCreateMPI(PETSC_COMM_WORLD,
-		 dofFirstLagrange[feSpace].nRankDofs * nComponents, 
-		 dofFirstLagrange[feSpace].nOverallDofs * nComponents,
+		 lagrangeMap[feSpace].nRankDofs * nComponents, 
+		 lagrangeMap[feSpace].nOverallDofs * nComponents,
 		 &(fetiData.tmp_vec_lagrange));
     VecCreateMPI(PETSC_COMM_WORLD, 
 		 primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
 		 &(fetiData.tmp_vec_primal));
 
     MatCreateShell(PETSC_COMM_WORLD,
-		   dofFirstLagrange[feSpace].nRankDofs * nComponents, 
-		   dofFirstLagrange[feSpace].nRankDofs * nComponents,
-		   dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
-		   dofFirstLagrange[feSpace].nOverallDofs * nComponents,
+		   lagrangeMap[feSpace].nRankDofs * nComponents, 
+		   lagrangeMap[feSpace].nRankDofs * nComponents,
+		   lagrangeMap[feSpace].nOverallDofs * nComponents, 
+		   lagrangeMap[feSpace].nOverallDofs * nComponents,
 		   &fetiData, &mat_feti);
     MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);
 
@@ -942,8 +934,8 @@ namespace AMDiS {
     for (int i = 0; i < nComponents; i++) {
       DOFVector<double>& dofVec = *(vec.getDOFVector(i));
 
-      for (DofMapping::iterator it = localIndexB.begin();
-	   it != localIndexB.end(); ++it) {
+      for (DofMapping::iterator it = localDofMap[feSpace].getMap().begin();
+	   it != localDofMap[feSpace].getMap().end(); ++it) {
 	int petscIndex = it->second * nComponents + i;
 	dofVec[it->first] = localSolB[petscIndex];
       }
@@ -973,7 +965,8 @@ namespace AMDiS {
 
     primalDofMap.setMpiComm(mpiComm);
     dualDofMap.setMpiComm(mpiComm);
-    dofFirstLagrange.setMpiComm(mpiComm);
+    lagrangeMap.setMpiComm(mpiComm);
+    localDofMap.setMpiComm(mpiComm);