And even more work....

AMDiS/src/CoarseningManager3d.cc

@@ -85,17 +85,18 @@ namespace AMDiS {
 
 
 #if HAVE_PARALLEL_DOMAIN_AMDIS
-    Element *otherEl = NULL;
-    int otherEdge = -1;      
-    FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge);
-    
+    vector<FixRefinementPatch::EdgeInEl> refineEdges;
+    FixRefinementPatch::getOtherEl(stack, refineEdges);
+
     // === If the refinement edge must be fixed, add also the other part of this ===
     // === edge to the refinement patch.                                         ===
     
-    if (otherEl) {
+    if (refineEdges.size()) {
       // TODO: Remove these two lines and make something more meaningful!!
       el->setMark(0);
       return;
+
+      Element *otherEl = refineEdges[0].first;
       
       TraverseStack stack2;
       ElInfo *elInfo2 = 

AMDiS/src/CoarseningManager3d.h

@@ -46,31 +46,25 @@ namespace AMDiS {
     /// Implements \ref CoarseningManager::coarsenFunction
     void coarsenFunction(ElInfo *el_info);
 
-    /** \brief
-     * Coarsens a single Tetrahedron of the coarsening patch. DOFs
-     * in the interior of the element are removed; DOFs for higher order
-     * at the boundary or the coarsening patch still belong to 
-     * the parent. Do not remove them form the mesh!!!                  
-     */
+    /// Coarsens a single Tetrahedron of the coarsening patch. DOFs
+    /// in the interior of the element are removed; DOFs for higher order
+    /// at the boundary or the coarsening patch still belong to 
+    /// the parent. Do not remove them form the mesh!!!                  
     void coarsenTetrahedron(RCNeighbourList &coarsenList, int index);
 
-    /** \brief
-     * Gets the patch for coarsening starting on element    
-     * el_info->el in direction of neighbour [3-dir]; returns 1 if a boundary
-     * reached and 0 if we come back to the starting element.         
-     * We complete the loop also in the case of a incompatible 
-     * coarsening patch since then all marks of patch elements are reset by    
-     * coarsenPatch() and this minimizes calls of traverseNeighbour();     
-     * if we reach a boundary while looping around the edge we loop back to    
-     * the starting element before we return                        
-     */         
+    /// Gets the patch for coarsening starting on element    
+    /// el_info->el in direction of neighbour [3-dir]; returns 1 if a boundary
+    /// reached and 0 if we come back to the starting element.         
+    /// We complete the loop also in the case of a incompatible 
+    /// coarsening patch since then all marks of patch elements are reset by    
+    /// coarsenPatch() and this minimizes calls of traverseNeighbour();     
+    /// if we reach a boundary while looping around the edge we loop back to    
+    /// the starting element before we return                        
     bool getCoarsenPatch(ElInfo* el_info, DegreeOfFreedom *edge[2],
 			 int dir, RCNeighbourList &coarsenList, int *n_neigh);
 
-    /** \brief
-     * First rebuild the DOFs on the parents then do restriction
-     * of data (if possible) and finally coarsen the patch elements
-     */
+    /// First rebuild the DOFs on the parents then do restriction
+    /// of data (if possible) and finally coarsen the patch elements
     void coarsenPatch(RCNeighbourList &coarsenList, int n_neigh, int bound);
     
   };

AMDiS/src/ProblemStat.cc

@@ -253,8 +253,10 @@ namespace AMDiS {
       }
     }
 
+#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
     if (!solver) 
       WARNING("no solver created\n");
+#endif
 
     // === create estimator ===
     if (initFlag.isSet(INIT_ESTIMATOR))

AMDiS/src/RefinementManager3d.cc

@@ -655,9 +655,14 @@ namespace AMDiS {
 
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    Element *otherEl = NULL;
-    int otherEdge = -1;
-    FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge);
+    vector<FixRefinementPatch::EdgeInEl> refineEdges;
+    FixRefinementPatch::getOtherEl(stack, refineEdges);
+    if (refineEdges.size()) {
+      MSG("FIX REFINEMENT PATH ON ELEMENT %d %d: %d additional edges\n", 
+	  elInfo->getElement()->getIndex(),
+	  elInfo->getMacroElement()->getIndex(),
+	  refineEdges.size());
+    }
 #endif
 
     // === Traverse and refine the refinement patch. ====
@@ -669,10 +674,20 @@ namespace AMDiS {
     }
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    // === If the refinement edge must be fixed, add also the other part of this ===
-    // === edge to the refinement patch.                                         ===
+    // === If the refinement edge must be fixed, add also the other part  ===
+    // === of this edge to the refinement patch.                          ===
+
+    for (int edgeIndex = 0; 
+	 edgeIndex < static_cast<unsigned int>(refineEdges.size()); edgeIndex++) {
+      MSG("  IN REF FIX: %d %d\n", refineEdges[edgeIndex].first->getIndex(), refineEdges[edgeIndex].second);
+    }
+
+    for (int edgeIndex = 0; 
+	 edgeIndex < static_cast<unsigned int>(refineEdges.size()); edgeIndex++) {
 
-    if (otherEl) {
+      MSG("  MAKE -> %d\n", edgeIndex);
+
+      Element *otherEl = refineEdges[edgeIndex].first;   
       TraverseStack stack2;
       ElInfo *elInfo2 = 
 	stack2.traverseFirstOneMacro(mesh, otherEl->getIndex(), -1, 
@@ -811,15 +826,14 @@ namespace AMDiS {
 
 
   void FixRefinementPatch::getOtherEl(TraverseStack *stack, 
-				      Element **otherEl, 
-				      int &otherEdge)
+				      vector<EdgeInEl>& refineEdges)
   {
     FUNCNAME("FixRefinementPatch::getOtherEl()");
 
     if (!FixRefinementPatch::connectedEdges.empty()) {
       // === Get stack of current traverse. ===
-      std::vector<ElInfo*> elInfos;
-      std::vector<int> infos;      
+      vector<ElInfo*> elInfos;
+      vector<int> infos;      
       int stackUsed = stack->getStackData(elInfos, infos);         
       int checkIndex = stackUsed;
       int localEdgeNo = 0;
@@ -849,28 +863,26 @@ namespace AMDiS {
 
 	checkIndex--;
       }
-
+      
      
       // If the refinement edge is part of an edge on the macro level, we must 
       // check if the refinement edge is part of an edge that must be fixed.
       if (localEdgeNo >= 0) {
+	int macroElIndex = elInfos[checkIndex]->getElement()->getIndex();
+
 	TEST_EXIT_DBG(elInfos[checkIndex]->getLevel() == 0)
 	  ("Should not happen!\n");
-	TEST_EXIT_DBG(elInfos[checkIndex]->getElement()->getIndex() ==
-		      elInfo->getMacroElement()->getIndex())
+	TEST_EXIT_DBG(macroElIndex == elInfo->getMacroElement()->getIndex())
 	  ("Should not happen!\n");
 	TEST_EXIT_DBG(localEdgeNo <= 5)("Should not happen!\n");
-	
-	for (unsigned int i = 0; i < FixRefinementPatch::connectedEdges.size(); i++) {
-	  if (FixRefinementPatch::connectedEdges[i].first.first->getIndex() ==
-	      elInfos[checkIndex]->getElement()->getIndex() &&
-	      FixRefinementPatch::connectedEdges[i].first.second ==
-	      localEdgeNo) {
-	    // Okay, we have found that this edge must be fixed.
-	    *otherEl = FixRefinementPatch::connectedEdges[i].second.first;
-	    otherEdge = FixRefinementPatch::connectedEdges[i].second.second;
-	    
-	    break;
+
+	refineEdges.clear();
+
+	for (int i = 0; i < static_cast<int>(connectedEdges.size()); i++) {
+	  if (connectedEdges[i].first.first->getIndex() == macroElIndex &&
+	      connectedEdges[i].first.second == localEdgeNo) {
+	    // We have found that this edge must be fixed.
+	    refineEdges.push_back(connectedEdges[i].second);
 	  }
 	}
       }

AMDiS/src/RefinementManager3d.h

@@ -86,7 +86,8 @@ namespace AMDiS {
 
     static ConnectedEdges connectedEdges;
 
-    static void getOtherEl(TraverseStack *stack, Element **otherEl, int &otherEdge);
+    static void getOtherEl(TraverseStack *stack, 
+			   vector<EdgeInEl> &refineEdges);
   };
 
 }

AMDiS/src/io/ArhReader.cc

@@ -102,9 +102,50 @@ namespace AMDiS {
   {
     FUNCNAME("ArhReader::readMeta()");
 
+    
+    // === Read the meta arh file. ===
+
+    string arhPrefix = "";
+    map<int, int> elInRank;
+    map<int, int> elCodeSize;
+    readMetaData(filename, elInRank, elCodeSize, arhPrefix);
+    
+
+    // === Check which arh files must be read by current rank. ===
+
+    // Set of all file indices which should be read to restore all macro elements.
+    std::set<int> readArhFiles;
+
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
+    while (elInfo) {
+      int macroElIndex = elInfo->getElement()->getIndex();
+      TEST_EXIT(elInRank.count(macroElIndex))("Should not happen!\n");
+      readArhFiles.insert(elInRank[macroElIndex]);
+      elInfo = stack.traverseNext(elInfo);
+    } 
+
+
+    // === Read the individual arh files. ===
+
     boost::filesystem::path p(filename.c_str());
     boost::filesystem::path directory = p.parent_path();
 
+    for (std::set<int>::iterator it = readArhFiles.begin();
+	 it != readArhFiles.end(); ++it) {
+      string arhFilename = 
+	directory.native() + "/" + arhPrefix + "-p" + boost::lexical_cast<string>(*it) + "-.arh";
+      MSG("ARH file read from: %s\n", arhFilename.c_str());
+      readFile(arhFilename, mesh, vecs);
+    }
+  }
+
+
+  int ArhReader::readMetaData(string filename,
+			      map<int, int> &elInRank,
+			      map<int, int> &elCodeSize,
+			      string &arhPrefix)
+  {
     ifstream file;
     file.open(filename.c_str());
     TEST_EXIT(file.is_open())
@@ -112,13 +153,12 @@ namespace AMDiS {
 
     string readStr = "";
     file >> readStr;
-    string arhPrefix = "";
+    arhPrefix = "";
     file >> arhPrefix;
     int nProc;
     file >> nProc;
 
     // Maps to each macro element index the arh file index it is stored in.
-    map<int, int> macroInProc;
     for (int i = 0; i < nProc; i++) {
       int tmp;
       file >> tmp;
@@ -126,34 +166,36 @@ namespace AMDiS {
       int nMacroEl;
       file >> nMacroEl;
       for (int j = 0; j < nMacroEl; j++) {
-	int elIndex;
+	int elIndex, codeSize;
 	file >> elIndex;
-	macroInProc[elIndex] = i;
+	file >> codeSize;
+	elInRank[elIndex] = i;
+	elCodeSize[elIndex] = codeSize;
       }
     }
 
     file.close();
 
+    return nProc;
+  }
 
-    // Set of all file indices which should be read to restore all macro elements.
-    std::set<int> readArhFiles;
 
-    TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
-    while (elInfo) {
-      int macroElIndex = elInfo->getElement()->getIndex();
-      TEST_EXIT(macroInProc.count(macroElIndex))("Should not happen!\n");
-      readArhFiles.insert(macroInProc[macroElIndex]);
-      elInfo = stack.traverseNext(elInfo);
-    } 
+  int ArhReader::readMetaData(string filename)
+  {
+    FUNCNAME("ArhReader::readMetaData()");
 
-    for (std::set<int>::iterator it = readArhFiles.begin();
-	 it != readArhFiles.end(); ++it) {
-      string arhFilename = 
-	directory.native() + "/" + arhPrefix + "-p" + boost::lexical_cast<string>(*it) + "-.arh";
-      MSG("READ FILE: %s\n", arhFilename.c_str());
-      readFile(arhFilename, mesh, vecs);
-    }
+    ifstream file;
+    file.open(filename.c_str());
+    TEST_EXIT(file.is_open())
+      ("Cannot open arh meta file \"%s\"\n", filename.c_str());
+
+    string readStr = "";
+    file >> readStr;
+    file >> readStr;
+    int nProc;
+    file >> nProc;
+    file.close();
+    return nProc;
   }
 
 

AMDiS/src/io/ArhReader.h

@@ -64,6 +64,22 @@ namespace AMDiS {
 			 Mesh *mesh,
 			 vector<DOFVector<double>*> vecs);
 
+    static int readMetaData(string filename,
+			    map<int, int> &elInRank,
+			    map<int, int> &elCodeSize,
+			    string &arhPrefix);
+
+    static int readMetaData(string filename,
+			    map<int, int> &elInRank,
+			    map<int, int> &elCodeSize)
+    {
+      string tmp;
+      return readMetaData(filename, elInRank, elCodeSize, tmp);
+    }
+
+    /// Returns just the number of subdomains a meta ARH file is defined for.
+    static int readMetaData(string filename);
+
     static void readFromMemoryBlock(vector<char> &data, Mesh *mesh,
 				    DOFVector<double>* vec0 = NULL,
 				    DOFVector<double>* vec1 = NULL,

AMDiS/src/parallel/MeshDistributor.cc

@@ -34,6 +34,7 @@
 #include "io/MacroInfo.h"
 #include "io/MacroWriter.h"
 #include "io/VtkWriter.h"
+#include "io/ArhReader.h"
 #include "Mesh.h"
 #include "Traverse.h"
 #include "ElInfo.h"
@@ -205,26 +206,9 @@ namespace AMDiS {
       return;
     }
 
-
-    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is
-    // supported only for macro meshes, so it will not work yet if the mesh is
-    // already refined in some way.
-    testForMacroMesh();
-
-    // For later mesh repartitioning, we need to store some information about
-    // the macro mesh.
-    createMacroElementInfo();
-
-    // create an initial partitioning of the mesh
-    partitioner->createInitialPartitioning();
-
-    // set the element weights, which are 1 at the very first begin
-    setInitialElementWeights();
-
-    // and now partition the mesh    
-    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
-    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
-    partitioner->createPartitionMap(partitionMap);
+    
+    // Create a very first pariitioning of the mesh.
+    createInitialPartitioning();
 
 
 #if (DEBUG != 0)    
@@ -357,6 +341,87 @@ namespace AMDiS {
   }
 
 
+  void MeshDistributor::createInitialPartitioning()
+  {
+    FUNCNAME("MeshDistributor::createInitialPartitioning()");
+
+    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is
+    // supported only for macro meshes, so it will not work yet if the mesh is
+    // already refined in some way.
+    testForMacroMesh();
+
+    // For later mesh repartitioning, we need to store some information about
+    // the macro mesh.
+    createMacroElementInfo();
+
+    // create an initial partitioning of the mesh
+    partitioner->createInitialPartitioning();
+
+    // set the element weights, which are 1 at the very first begin
+    setInitialElementWeights();
+
+    string filename = "";
+    Parameters::get("parallel->partitioner->read meta arh", filename);
+    if (filename == "" || ArhReader::readMetaData(filename) != mpiSize) {
+      // and now partition the mesh    
+      bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
+      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
+      partitioner->createPartitionMap(partitionMap);
+    }
+  }
+
+
+  void MeshDistributor::setInitialElementWeights() 
+  {
+    FUNCNAME("MeshDistributor::setInitialElementWeights()");
+
+    elemWeights.clear();
+      
+    string filename = "";
+    Parameters::get(mesh->getName() + "->macro weights", filename);
+    if (filename != "") {
+      MSG("Read macro weights from %s\n", filename.c_str());
+
+      ifstream infile;
+      infile.open(filename.c_str(), ifstream::in);
+      while (!infile.eof()) {
+	int elNum, elWeight;
+	infile >> elNum;
+	if (infile.eof())
+	  break;
+	infile >> elWeight;
+
+	elemWeights[elNum] = elWeight;
+      }
+
+      infile.close();
+
+      return;
+    }
+
+    filename = "";
+    Parameters::get("parallel->partitioner->read meta arh", filename);
+    if (filename != "") {
+      map<int, int> arhElInRank;
+      map<int, int> arhElCodeSize;
+
+      int nProc = ArhReader::readMetaData(filename, arhElInRank, arhElCodeSize);
+      for (map<int, int>::iterator it = arhElCodeSize.begin(); 
+	   it != arhElCodeSize.end(); ++it)
+	elemWeights[it->first] = it->second;
+
+      return;
+    }
+
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    while (elInfo) {
+      elemWeights[elInfo->getElement()->getIndex()] = 1.0;
+      elInfo = stack.traverseNext(elInfo);    
+    }
+  }
+
+
   void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
   {
     FUNCNAME("MeshDistributor::addProblemStat()");
@@ -1232,41 +1297,6 @@ namespace AMDiS {
   }
 
 
-  void MeshDistributor::setInitialElementWeights() 
-  {
-    FUNCNAME("MeshDistributor::setInitialElementWeights()");
-
-    elemWeights.clear();
-      
-    string filename = "";
-    Parameters::get(mesh->getName() + "->macro weights", filename);
-    if (filename != "") {
-      MSG("Read macro weights from %s\n", filename.c_str());
-
-      ifstream infile;
-      infile.open(filename.c_str(), ifstream::in);
-      while (!infile.eof()) {
-	int elNum, elWeight;
-	infile >> elNum;
-	if (infile.eof())
-	  break;
-	infile >> elWeight;
-
-	elemWeights[elNum] = elWeight;
-      }
-
-      infile.close();
-    } else {           
-      TraverseStack stack;
-      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
-      while (elInfo) {
-	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
-	elInfo = stack.traverseNext(elInfo);
-      }
-    }
-  }
-
-
   void MeshDistributor::repartitionMesh()
   {
     FUNCNAME("MeshDistributor::repartitionMesh()");

AMDiS/src/parallel/MeshDistributor.h

@@ -128,10 +128,6 @@ namespace AMDiS {
     /// the program if it finds a non macro element in the mesh.
     void testForMacroMesh();
 
-    /// Set for each element on the partitioning level the number of 
-    /// leaf elements.
-    void setInitialElementWeights();
-
     inline string getName() 
     { 
       return name; 
@@ -333,12 +329,21 @@ namespace AMDiS {
 				    const FiniteElemSpace *feSpace);
 
   protected:
+    /// Creates an initial paritioning of the mesh.
+    void createInitialPartitioning();
+
+    /// Set for each element on the partitioning level the number of 
+    /// leaf elements.
+    void setInitialElementWeights();
+
+    ///
     void addProblemStat(ProblemStatSeq *probStat);
 
     /// Determines the interior boundaries, i.e. boundaries between ranks, and
     /// stores all information about them in \ref interiorBoundary.
     void createInteriorBoundary(bool firstCall);
 
+    ///
     void createBoundaryDofs();
 
     /// Removes all macro elements from the mesh that are not part of ranks 

AMDiS/src/parallel/MeshPartitioner.cc

@@ -11,6 +11,7 @@
 
 
 #include "parallel/MeshPartitioner.h"
+#include "io/ArhReader.h"
 #include "Mesh.h"
 #include "Traverse.h"
 #include "Serializer.h"
@@ -26,6 +27,21 @@ namespace AMDiS {
     int nLeaves = mesh->getNumberOfLeaves();
     int elPerRank = nLeaves / mpiSize;
 
+
+    // === Check for reading ARH meta file to create initial partitioning. ===
+    
+    map<int, int> arhElInRank;
+    map<int, int> arhElCodeSize;
+    string arhMetaFile = "";
+    Parameters::get("parallel->partitioner->read meta arh", arhMetaFile);
+    bool partitioningArhBased = (arhMetaFile != "");
+    if (partitioningArhBased) {
+      int nProc = ArhReader::readMetaData(arhMetaFile, arhElInRank, arhElCodeSize);
+      if (nProc != mpiSize)
+	partitioningArhBased = false;
+    }
+    
+
     // === Create initial partitioning of the AMDiS mesh. ===
 
     elementInRank.clear();
@@ -53,10 +69,15 @@ namespace AMDiS {
       if (!boxPartitioning) {
 	// In standard mode assign to each macro element an arbitrary but unique
 	// rank number.
-	
-	int elInRank = std::min(elIndex / elPerRank, mpiSize - 1);
+
+	int elInRank = 0;
+	if (!partitioningArhBased)
+	  elInRank = std::min(elIndex / elPerRank, mpiSize - 1);
+ 	else
+ 	  elInRank = arhElInRank[elIndex];
+
 	elementInRank[elIndex] = (elInRank == mpiRank);
-	partitionMap[elIndex] = elInRank;	
+	partitionMap[elIndex] = elInRank;
       } else {
 	// In box partitioning mode, we do the assignment of boxes to ranks later.