Work on multi level methods, NEVER, NEVER, NEVER USE THIS RELEASE OTHERWISE YOU WILL BURN IN HELL.

AMDiS/src/Operator.cc

@@ -190,7 +190,7 @@ namespace AMDiS {
 
 
   void Operator::addTerm(FirstOrderTerm *term,
-           FirstOrderType type)
+			 FirstOrderType type)
   {
     addFirstOrderTerm(term, type);
   }

AMDiS/src/est/ResidualEstimator.cc

@@ -172,13 +172,13 @@ namespace AMDiS {
     mesh->getDofIndexCoords(coords);
 
     InteriorBoundary &intBoundary = 
-      MeshDistributor::globalMeshDistributor->getIntBoundary();
+      MeshDistributor::globalMeshDistributor->getIntBoundary(0);
 
     ElInfo *elInfo = mesh->createNewElInfo();
     elInfo->setFillFlag(Mesh::FILL_COORDS);
 
-    StdMpi<vector<double> > stdMpiDet(MeshDistributor::globalMeshDistributor->getMpiComm());
-    StdMpi<vector<vector<WorldVector<double> > > > stdMpiGrdUh(MeshDistributor::globalMeshDistributor->getMpiComm());
+    StdMpi<vector<double> > stdMpiDet(MeshDistributor::globalMeshDistributor->getMpiComm(0));
+    StdMpi<vector<vector<WorldVector<double> > > > stdMpiGrdUh(MeshDistributor::globalMeshDistributor->getMpiComm(0));
 
     RankToBoundMap allBounds = intBoundary.getOther();
     allBounds.insert(intBoundary.getOwn().begin(), intBoundary.getOwn().end());

AMDiS/src/parallel/DofComm.cc

@@ -16,32 +16,22 @@
 #include "FiniteElemSpace.h"
 #include "Debug.h"
 #include "ElementDofIterator.h"
+#include "DOFVector.h"
 
 namespace AMDiS {
 
   using namespace std;
 
 
-  void DofComm::init(int level, 
-		     MeshLevelData &ld,
-		     vector<const FiniteElemSpace*> &fe)
+  void DofComm::init(vector<const FiniteElemSpace*> &fe)
   {
     FUNCNAME("DofComm::init()");
     
-    meshLevel = level;
-    levelData = &ld;
     feSpaces = fe;
     
-    nLevel = levelData->getLevelNumber() - meshLevel;
-    TEST_EXIT_DBG(nLevel >= 1)("Should not happen!\n");
-
     sendDofs.clear();
     recvDofs.clear();
-    periodicDofs.clear();
-    
-    sendDofs.resize(nLevel);
-    recvDofs.resize(nLevel);
-    periodicDofs.resize(nLevel);
+    periodicDofs.clear();   
   }
 
 
@@ -55,40 +45,36 @@ namespace AMDiS {
 
 
   void DofComm::createContainer(RankToBoundMap &boundary,
-				LevelDataType &data)
+				DataType &data)
   {
     FUNCNAME("DofComm::createContainer()");
 
     // === Fill data. ===
 
     for (unsigned int i = 0; i < feSpaces.size(); i++)
-      for (int level = 0; level < nLevel; level++)
-		for (InteriorBoundary::iterator it(boundary, level); !it.end(); ++it){
-		  it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, 
-					     data[level][it.getRank()][feSpaces[i]]);
-		}
+      for (InteriorBoundary::iterator it(boundary); !it.end(); ++it)
+	it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, 
+				   data[it.getRank()][feSpaces[i]]);
 
 
     // === Remove empty data containers. ===
 
-    for (unsigned int i = 0; i < data.size(); i++) {
-      	DataIter dit = data[i].begin();
-      	while (dit != data[i].end()) {
-			FeMapIter it = dit->second.begin();
-			while (it != dit->second.end()) {
-			  if (it->second.size() == 0) {
-			    const FiniteElemSpace *fe = it->first;
-			    ++it;
-			    dit->second.erase(fe);
-			  } else
-			    ++it;
-			}
-	
-			if (dit->second.size() == 0)
-			  data[i].erase(dit++);
-			else
-			  ++dit;
-    	}
+    DataIter dit = data.begin();
+    while (dit != data.end()) {
+      FeMapIter it = dit->second.begin();
+      while (it != dit->second.end()) {
+	if (it->second.size() == 0) {
+	  const FiniteElemSpace *fe = it->first;
+	  ++it;
+	  dit->second.erase(fe);
+	} else
+	  ++it;
+      }
+      
+      if (dit->second.size() == 0)
+	data.erase(dit++);
+      else
+	++dit;
     }
   }
 
@@ -101,20 +87,16 @@ namespace AMDiS {
   }
 
   
-  int DofComm::getNumberDofs(LevelDataType &data, 
-			     int level, 
+  int DofComm::getNumberDofs(DataType &data, 
 			     const FiniteElemSpace *feSpace,
 			     bool countDouble)
   {
     FUNCNAME("DofComm::getNumberDofs()");
 
-    TEST_EXIT_DBG(level < data.size())("Should not happen!\n");
-
     DofContainerSet dofSet;
     DofContainer dofVec;
 
-    for (DataIter rankIt = data[level].begin(); 
-	 rankIt != data[level].end(); ++rankIt)
+    for (DataIter rankIt = data.begin(); rankIt != data.end(); ++rankIt)
       for (FeMapIter feIt = rankIt->second.begin();
 	   feIt != rankIt->second.end(); ++feIt)
 	if (feIt->first == feSpace)
@@ -133,7 +115,7 @@ namespace AMDiS {
   {
     FUNCNAME("DofComm::Iterator::setNextFeMap()");
 
-    if (dataIter != data[traverseLevel].end()) {
+    if (dataIter != data.end()) {
       TEST_EXIT_DBG(dataIter->second.size())("Should not happen!\n");
 
       feMapIter = dataIter->second.begin();
@@ -158,4 +140,25 @@ namespace AMDiS {
     return true;
   }
 
+
+  void MultiLevelDofComm::init(MeshLevelData &levelData,
+			       vector<const FiniteElemSpace*> &fe)
+  {
+    FUNCNAME("MultiLevelDofComm::init()");
+
+    int nLevel = levelData.getNumberOfLevels();
+    for (int level = 0; level < nLevel; level++)
+      levelDofComm[level].init(fe);
+  }
+
+  
+  void MultiLevelDofComm::create(MultiLevelInteriorBoundary &boundary)
+  {
+    FUNCNAME("MultiLevelDofComm::create()");
+
+    for (map<int, DofComm>::iterator it = levelDofComm.begin();
+	 it != levelDofComm.end(); ++it)
+      it->second.create(boundary[it->first]);
+  }
+
 }

AMDiS/src/parallel/DofComm.h

@@ -37,38 +37,28 @@ namespace AMDiS {
   {
   public:
     DofComm() 
-      : recvDofs(1),
-	sendDofs(1),
-	periodicDofs(0),
-	meshLevel(-1),
-	nLevel(0),
-	levelData(NULL)
-    {}
+      {}
     
     typedef map<const FiniteElemSpace*, DofContainer> FeMapType;
     typedef FeMapType::iterator FeMapIter;
     typedef map<int, FeMapType> DataType;
     typedef DataType::iterator DataIter;
-    // meshLevel: map[rank -> map[feSpace -> DofContainer]]
-    typedef vector<DataType> LevelDataType;
 
-    void init(int level, 
-	      MeshLevelData &levelData, 
-	      vector<const FiniteElemSpace*> &fe);
+    void init(vector<const FiniteElemSpace*> &fe);
 
     void create(InteriorBoundary &boundary);
 
-    LevelDataType& getSendDofs()
+    DataType& getSendDofs()
     {
       return sendDofs;
     }
 
-    LevelDataType& getRecvDofs()
+    DataType& getRecvDofs()
     {
       return recvDofs;
     }
 
-    LevelDataType& getPeriodicDofs()
+    DataType& getPeriodicDofs()
     {
       return periodicDofs;
     }
@@ -76,35 +66,29 @@ namespace AMDiS {
     // Writes all data of this object to an output stream.
     void serialize(ostream &out);
 
-    int getNumberDofs(LevelDataType &data, 
-		      int level, 
+    int getNumberDofs(DataType &data, 
 		      const FiniteElemSpace *feSpace,
 		      bool countDouble = false);
 
   protected:
-    void createContainer(RankToBoundMap &boundary, LevelDataType &data);
+    void createContainer(RankToBoundMap &boundary, 
+			 DataType &data);
 
   protected:
     /// This map contains for each rank the list of DOFs the current rank must 
     /// end to exchange solution DOFs at the interior boundaries.
-    LevelDataType sendDofs;
+    DataType sendDofs;
 
     /// This map contains on each rank the list of DOFs from which the current 
     /// rank will receive DOF values (i.e., this are all DOFs at an interior 
     /// boundary). The DOF indices are given in rank's local numbering.
-    LevelDataType recvDofs;
+    DataType recvDofs;
 
     /// This map contains on each rank a list of DOFs along the interior bound-
     /// aries to communicate with other ranks. The DOF indices are given in rank's
     /// local numbering. Periodic boundaries within one subdomain are not 
     /// considered here. 
-    LevelDataType periodicDofs;
-
-    int meshLevel;
-
-    int nLevel;
-
-    MeshLevelData *levelData;
+    DataType periodicDofs;
 
     vector<const FiniteElemSpace*> feSpaces;
 
@@ -114,24 +98,11 @@ namespace AMDiS {
     class Iterator
     {
     public:
-      Iterator(LevelDataType &d,
-	       const FiniteElemSpace *fe = NULL)
-	: data(d),
-	  dofCounter(-1),
-	  traverseFeSpace(fe),
-	  traverseLevel(0),
-	  removedDof(false)
-      {
-	goFirst();
-      }
-
-      Iterator(LevelDataType &d,
-	       int level,
+      Iterator(DataType &d,
 	       const FiniteElemSpace *fe = NULL)
 	: data(d),
 	  dofCounter(-1),
 	  traverseFeSpace(fe),
-	  traverseLevel(level),
 	  removedDof(false)
       {
 	goFirst();
@@ -139,7 +110,7 @@ namespace AMDiS {
 
       inline bool end()
       {
-	return (dataIter == data[traverseLevel].end());
+	return (dataIter == data.end());
       }
       
       inline void nextRank()
@@ -242,7 +213,7 @@ namespace AMDiS {
     protected:
       void goFirst()
       {
-	dataIter = data[traverseLevel].begin();
+	dataIter = data.begin();
 
 	while (setNextFeMap() == false)
 	  ++dataIter;
@@ -251,7 +222,7 @@ namespace AMDiS {
       bool setNextFeMap();
 
     protected:
-      LevelDataType &data;
+      DataType &data;
       
       DofComm::DataIter dataIter;
       
@@ -263,12 +234,26 @@ namespace AMDiS {
 
       const FiniteElemSpace *traverseFeSpace;
 
-      int traverseLevel;
-
       bool removedDof;
     };
 
+  };
+
+  
+  class MultiLevelDofComm {
+  public:
+    void init(MeshLevelData &levelData,
+	      vector<const FiniteElemSpace*> &fe);
+
+    void create(MultiLevelInteriorBoundary &boundary);
+
+    inline DofComm& operator[](int level) 
+    {
+      return levelDofComm[level];
+    }
 
+  private:
+    map<int, DofComm> levelDofComm;
   };
 
 }

AMDiS/src/parallel/ElementObjectDatabase.cc

@@ -722,11 +722,13 @@ namespace AMDiS {
   int ElementObjectDatabase::getOwner(vector<ElementObjectData>& objData, 
 				      int level)
   {
+    FUNCNAME("ElementObjectDatabase::getOwner()");
+
     int owner = -1;
 
     std::set<int> &levelRanks = levelData->getLevelRanks(level);
-    bool allRanks = (levelRanks.size() == 1 && *(levelRanks.begin()) == -1);
-    
+    bool allRanks = (level == 0);
+   
     for (vector<ElementObjectData>::iterator it = objData.begin();
 	 it != objData.end(); ++it) {
       int elRank = (*macroElementRankMap)[it->elIndex];
@@ -744,7 +746,7 @@ namespace AMDiS {
   {
     FUNCNAME("ElementObjectDatabase::getIterateMaxLevel()");
 
-    int nLevel = levelData->getLevelNumber();
+    int nLevel = levelData->getNumberOfLevels();
     TEST_EXIT_DBG(nLevel > 0)("Should not happen!\n");
 
     vector<std::set<int> > ranksInLevel;

AMDiS/src/parallel/ElementObjectDatabase.h

@@ -143,6 +143,12 @@ namespace AMDiS {
     /// All data from the database is dropped. 
     void clear();
 
+    /// Resets iteration;
+    bool resetIterator()
+    {
+      iterGeoPos = CENTER;
+    }
+
     /** \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
@@ -151,7 +157,7 @@ namespace AMDiS {
      * \param[in]  pos   Must be either VERTEX, EDGE or FACE and defines the
      *                   elements that should be traversed.
      */
-    bool iterate(GeoIndex pos)
+    inline bool iterate(GeoIndex pos)
     {
       // CENTER marks the variable "iterGeoPos" to be in an undefined state. I.e.,
       // there is no iteration that is actually running.

AMDiS/src/parallel/InteriorBoundary.cc

@@ -35,19 +35,24 @@ namespace AMDiS {
 
     Mesh *mesh = elObjDb.getMesh();
     TEST_EXIT_DBG(mesh)("Should not happen!\n");
-    TEST_EXIT_DBG(level < levelData.getLevelNumber())("Should not happen!\n");
+    TEST_EXIT_DBG(level < levelData.getNumberOfLevels())
+      ("Should not happen!\n");
 
     MPI::Intracomm mpiComm = levelData.getMpiComm(level);
+
+    if (mpiComm == MPI::COMM_SELF)
+      return;
+
     int levelMpiRank = mpiComm.Get_rank();
     int globalMpiRank = MPI::COMM_WORLD.Get_rank();
     std::set<int> levelRanks = levelData.getLevelRanks(level);
-    bool allRanks = (levelRanks.size() == 1 && *(levelRanks.begin()) == -1);
 
     // === Create interior boundary data structure. ===
     
     for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) {
       GeoIndex geoIndex = INDEX_OF_DIM(geoPos, mesh->getDim());
 
+      elObjDb.resetIterator();
       while (elObjDb.iterate(geoIndex)) {
 	flat_map<int, ElementObjectData>& objData = elObjDb.getIterateData();
 
@@ -58,17 +63,22 @@ namespace AMDiS {
 	// Test, if the boundary object defines an interior boundary within the
 	// ranks of the MPI group. If not, go to next element.
 	bool boundaryWithinMpiGroup = false;
-	if (allRanks) {
+	if (levelData.getNumberOfLevels() == 1) {
 	  boundaryWithinMpiGroup = true;
 	} else {
+	  TEST_EXIT_DBG(level + 1 < levelData.getNumberOfLevels())
+	    ("Level = %d    but number of levels = %d\n",
+	     level, levelData.getNumberOfLevels());
+
 	  for (flat_map<int, ElementObjectData>::iterator it = objData.begin();
 	       it != objData.end(); ++it) {
-	    if (it->first != globalMpiRank && levelRanks.count(it->first)) {
+	    if (!levelData.rankInLevel(it->first, level + 1)) {
 	      boundaryWithinMpiGroup = true;
 	      break;
 	    }
 	  }
 	}
+	
 	if (!boundaryWithinMpiGroup)
 	  continue;
 
@@ -100,7 +110,7 @@ namespace AMDiS {
 	    if (it2->first == globalMpiRank)
 	      continue;
 
-	    if (!allRanks && levelRanks.count(it2->first) == 0)
+	    if (!levelData.rankInLevel(it2->first, level))
 	      continue;
 
 	    bound.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
@@ -111,7 +121,8 @@ namespace AMDiS {
 	    
 	    bound.type = INTERIOR;
 
-	    AtomicBoundary& b = getNewOwn(it2->first);
+	    int rankInLevel = levelData.mapRank(it2->first, 0, level);
+	    AtomicBoundary& b = getNewOwn(rankInLevel);
 	    b = bound;
 	    if (geoIndex == EDGE)
 	      b.neighObj.reverseMode = 
@@ -134,8 +145,9 @@ namespace AMDiS {
 	  bound.neighObj.ithObj = ownerBoundEl.ithObject;
 	  
 	  bound.type = INTERIOR;
-	  
-	  AtomicBoundary& b = getNewOther(owner);
+
+	  int rankInLevel = levelData.mapRank(owner, 0, level);	  
+	  AtomicBoundary& b = getNewOther(rankInLevel);
 	  b = bound;	    
 	  if (geoIndex == EDGE)
 	    b.rankObj.reverseMode =
@@ -224,23 +236,26 @@ namespace AMDiS {
 	  // interior database.
 
 	  if (owner == globalMpiRank) {
-	    AtomicBoundary& b = getNewOwn(otherElementRank);	  
+	    int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
+	    AtomicBoundary& b = getNewOwn(rankInLevel);	  
 	    b = bound;	  
 	  } else {
-	    ElementObjectData& ownerBoundEl = objData[owner];    
+	    ElementObjectData& ownerBoundEl = objData[owner];
 	    bound.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
 	    bound.neighObj.elIndex = ownerBoundEl.elIndex;
 	    bound.neighObj.elType = -1;
 	    bound.neighObj.ithObj = ownerBoundEl.ithObject;
 
-	    AtomicBoundary& b = getNewOther(owner);
+	    int rankInLevel = levelData.mapRank(owner, 0, level);
+	    AtomicBoundary& b = getNewOther(rankInLevel);
 	    b = bound;
 	  }
 
 	} else {
 	  bound.type = it->second;
 	  
-	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
+	  int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
+	  AtomicBoundary& b = getNewPeriodic(rankInLevel);
 	  b = bound;	  
 	}
       }
@@ -293,10 +308,12 @@ namespace AMDiS {
 	  }
 
 	  if (owner == globalMpiRank) {
-	    AtomicBoundary& b = getNewOwn(otherElementRank); 
+	    int rankInLevel = levelData.mapRank(owner, 0, level);
+	    AtomicBoundary& b = getNewOwn(rankInLevel); 
 	    b = bound;	  
 	  } else {
-	    AtomicBoundary& b = getNewOther(owner);
+	    int rankInLevel = levelData.mapRank(owner, 0, level);
+	    AtomicBoundary& b = getNewOther(rankInLevel);
 	    b = bound;
 	    
 	    ElementObjectData& ownerBoundEl = objData[owner];    
@@ -309,8 +326,9 @@ namespace AMDiS {
 	  }
 	} else {
 	  bound.type = it->second;
-	  
-	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
+
+	  int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
+	  AtomicBoundary& b = getNewPeriodic(rankInLevel);
 	  b = bound;
 	  
 	  if (globalMpiRank > otherElementRank)
@@ -369,10 +387,12 @@ namespace AMDiS {
 	  ElementObjectData& rankBoundEl = objData[globalMpiRank];
 
 	  if (owner == globalMpiRank) {
-	    AtomicBoundary& b = getNewOwn(otherElementRank);
+	    int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
+	    AtomicBoundary& b = getNewOwn(rankInLevel);
 	    b = bound;	  
 	  } else {
-	    AtomicBoundary& b = getNewOther(owner);
+	    int rankInLevel = levelData.mapRank(owner, 0, level);
+	    AtomicBoundary& b = getNewOther(rankInLevel);
 	    b = bound;
 	    
 	    ElementObjectData& ownerBoundEl = objData[owner];    
@@ -386,7 +406,8 @@ namespace AMDiS {
 	} else {
 	  bound.type = it->second;
 	  
-	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
+	  int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
+	  AtomicBoundary& b = getNewPeriodic(rankInLevel);
 	  b = bound;
 	  
 	  if (globalMpiRank > otherElementRank)
@@ -406,18 +427,8 @@ namespace AMDiS {
     // === share the bounday.                                                 ===
 
     StdMpi<vector<AtomicBoundary> > stdMpi(mpiComm);
-    if (level == 0) {
-      stdMpi.send(own);
-      stdMpi.recv(other);
-    } else {
-      for (RankToBoundMap::iterator rankIt = own.begin();	 
-	   rankIt != own.end(); ++rankIt)
-	stdMpi.send(levelData.mapRank(rankIt->first, 0, level), rankIt->second);
-      for (RankToBoundMap::iterator rankIt = other.begin();
-	   rankIt != other.end(); ++rankIt)
-	stdMpi.recv(levelData.mapRank(rankIt->first, 0, level));
-    }
-
+    stdMpi.send(own);
+    stdMpi.recv(other);
     stdMpi.startCommunication();
 
 
@@ -430,8 +441,6 @@ namespace AMDiS {
 	 rankIt != other.end(); ++rankIt) {
 
       int rank = rankIt->first;
-      if (level > 0)
-	rank = levelData.mapRank(rank, 0, level);
 
       // === We have received from "rank" the ordered list of element       ===
       // === indices. Now, we have to sort the corresponding list in this   ===
@@ -819,4 +828,31 @@ namespace AMDiS {
     }
   }
 
+
+  void MultiLevelInteriorBoundary::create(MeshLevelData &levelData,
+					  ElementObjectDatabase &elObjDb)
+  {
+    FUNCNAME(" MultLevelInteriorBoundary::create()");
+
+    levelIntBound.clear();
+
+    int nLevel = levelData.getNumberOfLevels();
+    for (int level = 0; level < nLevel; level++)
+      levelIntBound[level].create(levelData, level, elObjDb);
+  }
+
+
+  void MultiLevelInteriorBoundary::serialize(ostream &out)
+  {
+    FUNCNAME("MultiLevelInteriorBoundary::serialize()");
+    ERROR_EXIT("Not yet implemented!\n");
+  }
+
+
+  void MultiLevelInteriorBoundary::deserialize(istream &in, Mesh *mesh)