Work on multi-level FETI-DP method.

AMDiS/src/parallel/DofComm.cc

@@ -18,22 +18,26 @@ namespace AMDiS {
 
   void DofComm::removeEmpty()
   {
-    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;
+    FUNCNAME("DofComm::removeEmpty()");
+
+    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;
       }
-      
-      if (dit->second.size() == 0)
-	data.erase(dit++);
-      else
-	++dit;
     }
   }
 
@@ -42,7 +46,7 @@ namespace AMDiS {
   {
     FUNCNAME("DofComm::Iterator::setNextFeMap()");
 
-    if (dataIter != dofComm.data.end()) {
+    if (dataIter != dofComm.data[0].end()) {
       TEST_EXIT_DBG(dataIter->second.size())("Should not happen!\n");
 
       feMapIter = dataIter->second.begin();

AMDiS/src/parallel/DofComm.h

@@ -40,29 +40,33 @@ namespace AMDiS {
     
     typedef map<const FiniteElemSpace*, DofContainer> FeMapType;
     typedef FeMapType::iterator FeMapIter;
-    // meshLevel: map[rank -> map[feSpace -> DofContainer]]
-    typedef vector<int, map<int, FeMapType> > DataType;
+    typedef map<int, FeMapType> DataType;
     typedef DataType::iterator DataIter;
+    // meshLevel: map[rank -> map[feSpace -> DofContainer]]
+    typedef vector<DataType> LevelDataType;
 
-    DofContainer& getDofContainer(int rank, const FiniteElemSpace *feSpace)
+    inline DofContainer& getDofContainer(int rank, 
+					 const FiniteElemSpace *feSpace, 
+					 int level = 0)
     {
-      return data[0][rank][feSpace];
+      return data[level][rank][feSpace];
     }
 
     void removeEmpty();
 
-    void clear()
+    void init(int nLevels = 1)
     {
       data.clear();
+      data.resize(nLevels);
     }
 
-    DataType& getData()
+    DataType& getData(int level = 0)
     {
-      return data;
+      return data[level];
     }
 
   protected:
-    DataType data;
+    LevelDataType data;
 
     friend class Iterator;
 

AMDiS/src/parallel/MeshDistributor.cc

@@ -87,7 +87,8 @@ namespace AMDiS {
       debugOutputDir(""),
       lastMeshChangeIndex(0),
       createBoundaryDofFlag(0),
-      sebastianMode(false)
+      sebastianMode(false),
+      boundaryDofInfo(1)
   {
     FUNCNAME("MeshDistributor::ParalleDomainBase()");
 
@@ -293,6 +294,11 @@ namespace AMDiS {
 	 it != mesh->getPeriodicAssociations().end(); ++it)
       const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
 
+
+    // If required, create hierarchical mesh level structure.
+    createMeshLevelStructure();
+
+
     updateLocalGlobalNumbering();
 
     // === In 3D we have to make some test, if the resulting mesh is valid.  ===
@@ -326,9 +332,6 @@ namespace AMDiS {
     ParallelDebug::testPeriodicBoundary(*this);
 #endif
 
-    // If required, create hierarchical mesh level structure.
-    createMeshLevelStructure();
-
     // === Global refinements. ===
     
     int globalRefinement = 0;
@@ -780,9 +783,6 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::createMeshLevelStructure()");
 
-    if (mpiSize != 16)
-      return;
-
     std::set<int> neighbours;
     for (InteriorBoundary::iterator it(rankIntBoundary); !it.end(); ++it)
       neighbours.insert(it.getRank());
@@ -792,6 +792,9 @@ namespace AMDiS {
 
     levelData.init(neighbours);
 
+    if (mpiSize != 16)
+      return;
+
     bool multiLevelTest = false;
     Parameters::get("parallel->multi level test", multiLevelTest);
     if (multiLevelTest) {
@@ -1871,10 +1874,13 @@ namespace AMDiS {
   {
     FUNCNAME("MeshDistributor::createBoundaryDofs()");
 
-    sendDofs.clear();
-    recvDofs.clear();
-
     int nLevels = levelData.getLevelNumber();
+    TEST_EXIT_DBG(nLevels >= 1)("Should not happen!\n");
+
+    sendDofs.init(nLevels);
+    recvDofs.init(nLevels);
+
+    boundaryDofInfo.resize(nLevels);
 
     for (unsigned int i = 0; i < feSpaces.size(); i++)
       for (int j = 0; j < nLevels; j++)
@@ -1887,42 +1893,43 @@ namespace AMDiS {
     FUNCNAME("MeshDistributor::createBoundaryDofs()");
 
     if (createBoundaryDofFlag.isSet(BOUNDARY_SUBOBJ_SORTED)) {
-      
-      TEST_EXIT(level == 0)
-	("This function does not support the usage of multi level structure!\n");
+
+      TEST_EXIT(level < boundaryDofInfo.size())("Should not happen!\n");
 
       // === Clear data. ===
       for (int geo = FACE; geo >= VERTEX; geo--)
-	boundaryDofInfo[feSpace].geoDofs[static_cast<GeoIndex>(geo)].clear();
+	boundaryDofInfo[level][feSpace].geoDofs[static_cast<GeoIndex>(geo)].clear();
 
       // === Create send DOFs. ===
       for (int geo = FACE; geo >= VERTEX; geo--) {
-	for (InteriorBoundary::iterator it(rankIntBoundary); !it.end(); ++it) {
+	for (InteriorBoundary::iterator it(rankIntBoundary, levelData, level); 
+	     !it.end(); ++it) {
 	  if (it->rankObj.subObj == geo) {
 	    DofContainer dofs;
 	    it->rankObj.el->getAllDofs(feSpace, it->rankObj, dofs);
 
-	    DofContainer& tmp = sendDofs.getDofContainer(it.getRank(), feSpace);
+	    DofContainer& tmp = sendDofs.getDofContainer(it.getRank(), feSpace, level);
 	    tmp.insert(tmp.end(), dofs.begin(), dofs.end());
 
 	    if (createBoundaryDofFlag.isSet(BOUNDARY_FILL_INFO_SEND_DOFS))
-	      boundaryDofInfo[feSpace].geoDofs[static_cast<GeoIndex>(geo)].insert(dofs.begin(), dofs.end());
+	      boundaryDofInfo[level][feSpace].geoDofs[static_cast<GeoIndex>(geo)].insert(dofs.begin(), dofs.end());
 	  }
 	}
       }
 
       // === Create recv DOFs. ===
       for (int geo = FACE; geo >= VERTEX; geo--) {
-	for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
+	for (InteriorBoundary::iterator it(otherIntBoundary, levelData, level); 
+	     !it.end(); ++it) {
 	  if (it->rankObj.subObj == geo) {
 	    DofContainer dofs;
 	    it->rankObj.el->getAllDofs(feSpace, it->rankObj, dofs);
 
-	    DofContainer& tmp = recvDofs.getDofContainer(it.getRank(), feSpace);
+	    DofContainer& tmp = recvDofs.getDofContainer(it.getRank(), feSpace, level);
 	    tmp.insert(tmp.end(), dofs.begin(), dofs.end());
 
 	    if (createBoundaryDofFlag.isSet(BOUNDARY_FILL_INFO_RECV_DOFS))
-	      boundaryDofInfo[feSpace].geoDofs[static_cast<GeoIndex>(geo)].insert(dofs.begin(), dofs.end());
+	      boundaryDofInfo[level][feSpace].geoDofs[static_cast<GeoIndex>(geo)].insert(dofs.begin(), dofs.end());
 	  }
 	}
       }
@@ -1930,12 +1937,12 @@ namespace AMDiS {
       for (InteriorBoundary::iterator it(rankIntBoundary, levelData, level);
 	   !it.end(); ++it)
 	it->rankObj.el->getAllDofs(feSpace, it->rankObj, 
-				   sendDofs.getDofContainer(it.getRank(), feSpace));
+				   sendDofs.getDofContainer(it.getRank(), feSpace, level));
       
       for (InteriorBoundary::iterator it(otherIntBoundary, levelData, level); 
 	   !it.end(); ++it)
 	it->rankObj.el->getAllDofs(feSpace, it->rankObj, 
-				   recvDofs.getDofContainer(it.getRank(), feSpace));
+				   recvDofs.getDofContainer(it.getRank(), feSpace, level));
     }
 
     // === Delete all empty DOF send and recv positions ===
@@ -1968,8 +1975,12 @@ namespace AMDiS {
     debug::createSortedDofs(mesh, elMap);   
 #endif
 
-    sendDofs.clear();
-    recvDofs.clear();
+    int nLevels = levelData.getLevelNumber();
+    TEST_EXIT_DBG(nLevels >= 1)("Should not happen!\n");
+
+    sendDofs.init(nLevels);
+    recvDofs.init(nLevels);
+    boundaryDofInfo.resize(nLevels);
 
     for (unsigned int i = 0; i < feSpaces.size(); i++)
       updateLocalGlobalNumbering(feSpaces[i]);
@@ -2116,7 +2127,7 @@ namespace AMDiS {
     FUNCNAME("MeshDistributor::createPeriodicMap()");
 
     // Clear all periodic DOF mappings calculated before. We do it from scratch.
-    periodicDofs.clear();
+    periodicDofs.init(levelData.getLevelNumber());
     periodicMap.clear();
 
     // If there are no periodic boundaries, return. Note that periodicDofs and

AMDiS/src/parallel/MeshDistributor.h

@@ -436,9 +436,15 @@ namespace AMDiS {
       createBoundaryDofFlag = flag;
     }
 
-    BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace)
+    BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace, int level = 0)
     {
-      return boundaryDofInfo[feSpace];
+      FUNCNAME("MeshDistributor::getBoundaryDofInfo()");
+
+      TEST_EXIT_DBG(level < static_cast<int>(boundaryDofInfo.size()))
+	("Wrong level number: %d, whereas array size is %d!\n", 
+	 level, boundaryDofInfo.size());
+
+      return boundaryDofInfo[level][feSpace];
     }
 
     void getAllBoundaryDofs(const FiniteElemSpace *feSpace,
@@ -738,7 +744,9 @@ namespace AMDiS {
 
     Flag createBoundaryDofFlag;
 
-    map<const FiniteElemSpace*, BoundaryDofInfo> boundaryDofInfo;
+    /// Stores on each mesh level for all FE spaces the information about
+    /// all boundary DOFs.
+    vector<map<const FiniteElemSpace*, BoundaryDofInfo> > boundaryDofInfo;
 
     MeshLevelData levelData;
 

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -56,7 +56,7 @@ namespace AMDiS {
 
     MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
     KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
-    MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
+   MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
     data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
     MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);