Changed coarse space creation for FETI-DP method.

AMDiS/src/AMDiS.cc

@@ -22,7 +22,7 @@ namespace AMDiS {
   using namespace std;
 
 #if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
-  mtl::par::environment* mtl_environment=NULL;
+  mtl::par::environment* mtl_environment = NULL;
 #endif
 
   void init(int argc, char **argv, std::string initFileName)
@@ -37,8 +37,6 @@ namespace AMDiS {
     mpi::startRand();
 #endif
 
-    Parameters::get("parallel->log main rank", Msg::outputMainRank);
-
 #ifdef HAVE_ZOLTAN
     float zoltanVersion = 0.0;
     Zoltan_Initialize(argc, argv, &zoltanVersion);
@@ -98,6 +96,8 @@ namespace AMDiS {
       Parameters::init(initFileName);
     }
 
+    Parameters::get("parallel->log main rank", Msg::outputMainRank);
+
     // reset parameters from command line
     bool ignoreCommandline = false;
     Parameters::get("ignore commandline options", ignoreCommandline);

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -734,51 +734,78 @@ namespace AMDiS {
   }
 
 
-  std::set<BoundaryObject> PetscSolverFeti::getAugmentedLagrange()
+  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseEdges()
   {
     FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
 
-    int nEdges = 0;
-    int nFaces = 0;
-
     InteriorBoundary &intBound = meshDistributor->getIntBoundary();
     std::set<BoundaryObject> allEdges;
-    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
-      if ((it->rankObj.subObj == FACE ||
-	   (it->rankObj.subObj == EDGE && 
-	    testWirebasketEdge(it->rankObj, feSpaces[0]))) && 
-	  allEdges.count(it->rankObj) == 0) {
-	bool dirichletOnlyEdge = true;
-
-	DofContainer edgeDofs;
-	it->rankObj.el->getAllDofs(feSpaces[0], it->rankObj, edgeDofs);
-	
-	for (DofContainer::iterator dit = edgeDofs.begin();
-	     dit != edgeDofs.end(); ++dit) {
-	  if (dirichletRows[0].count(**dit) == 0) {
-	    dirichletOnlyEdge = false;
-	    break;
-	  }
-	}
+    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
+      if (it->rankObj.subObj == EDGE && 
+	  testWirebasketEdge(it->rankObj, feSpaces[0]) &&
+	  allEdges.count(it->rankObj) == 0)
+	allEdges.insert(it->rankObj);
+
+    int nEmptyEdges = 0;
+    vector<vector<BoundaryObject> > data;
+    for (std::set<BoundaryObject>::iterator it = allEdges.begin();
+	 it != allEdges.end(); ++it) {
+      DofContainer edgeDofs;
+      it->el->getAllDofs(feSpaces[0], *it, edgeDofs);
+      if (edgeDofs.size() == 0) {
+	nEmptyEdges++;
+      } else {
+	vector<BoundaryObject> oneBoundary;
+	oneBoundary.push_back(*it);      
+	data.push_back(oneBoundary);
+      }
+    }
 
-	if (!dirichletOnlyEdge) {
-	  if (it->rankObj.subObj == EDGE)
-	    nEdges++;
-	  
-	  if (it->rankObj.subObj == FACE)
-	    nFaces++;
+    int nEdges = allEdges.size();
+    mpi::globalAdd(nEdges);
+    mpi::globalAdd(nEmptyEdges);
 
-	  allEdges.insert(it->rankObj);
-	}
+    MSG("Coarse space edges: %d (empty: %d)\n", nEdges, nEmptyEdges);
+    
+    return data;
+  }
+
+
+  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseFaces()
+  {
+    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
+
+    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
+    map<int, std::set<BoundaryObject> > allFaces;
+    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
+      if (it->rankObj.subObj == FACE)
+	allFaces[it.getRank()].insert(it->rankObj);
+
+    int nEmptyFaces = 0;
+    vector<vector<BoundaryObject> > data;
+    for (map<int, std::set<BoundaryObject> >::iterator it = allFaces.begin();
+	 it != allFaces.end(); ++it) {
+      vector<BoundaryObject> oneBoundary;
+      for (std::set<BoundaryObject>::iterator bIt = it->second.begin();
+	   bIt != it->second.end(); ++bIt) {
+	DofContainer faceDofs;
+	bIt->el->getAllDofs(feSpaces[0], *bIt, faceDofs);
+	if (faceDofs.size())
+	  oneBoundary.push_back(*bIt);
       }
+      if (oneBoundary.size()) 
+	data.push_back(oneBoundary);
+      else 
+	nEmptyFaces++;
     }
 
-    mpi::globalAdd(nEdges);
+    int nFaces = allFaces.size();
     mpi::globalAdd(nFaces);
+    mpi::globalAdd(nEmptyFaces);
 
-    MSG("NEDGES = %d  NFACES = %d\n", nEdges, nFaces);   
-
-    return allEdges;
+    MSG("Coarse space faces: %d (empty: %d)\n", nFaces, nEmptyFaces);
+   
+    return data;
   }
 
 
@@ -791,8 +818,11 @@ namespace AMDiS {
 
     double wtime = MPI::Wtime();
 
-    std::set<BoundaryObject> allEdges = getAugmentedLagrange();
-    nRankEdges = allEdges.size();    
+    vector<vector<BoundaryObject> > allEdges = getCoarseEdges();
+    vector<vector<BoundaryObject> > allFaces = getCoarseFaces(); 
+    allEdges.insert(allEdges.end(), allFaces.begin(), allFaces.end());
+
+    nRankEdges = allEdges.size();
     int rStartEdges = 0;
     mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);
 
@@ -810,27 +840,29 @@ namespace AMDiS {
     MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 
     int rowCounter = rStartEdges;
-    for (std::set<BoundaryObject>::iterator edgeIt = allEdges.begin(); 
-	 edgeIt != allEdges.end(); ++edgeIt) {
+    for (vector<vector<BoundaryObject> >::iterator it = allEdges.begin(); 
+	 it != allEdges.end(); ++it) {
       for (int component = 0; component < componentSpaces.size(); component++) {
-	DofContainer edgeDofs;
-	edgeIt->el->getAllDofs(componentSpaces[component], *edgeIt, edgeDofs);
-	
-	for (DofContainer::iterator it = edgeDofs.begin();
-	     it != edgeDofs.end(); ++it) {
-	  TEST_EXIT_DBG(isPrimal(component, **it) == false)
-	    ("Should not be primal!\n");
+	for (vector<BoundaryObject>::iterator edgeIt = it->begin();
+	     edgeIt != it->end(); ++edgeIt) {
 	  
-	  if (dirichletRows[component].count(**it))
-	    continue;
+	  DofContainer edgeDofs;
+	  edgeIt->el->getAllDofs(componentSpaces[component], *edgeIt, edgeDofs);
+	  
+	  TEST_EXIT(edgeDofs.size())("Should not happen!\n");
 	  
-	  int col = lagrangeMap.getMatIndex(component, **it);
-	  double value = 1.0;
-	  MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
+	  for (DofContainer::iterator it = edgeDofs.begin();
+	       it != edgeDofs.end(); ++it) {
+	    TEST_EXIT(isPrimal(component, **it) == false)
+	      ("Should not be primal!\n");
+	    
+	    int col = lagrangeMap.getMatIndex(component, **it);
+	    double value = 1.0;
+	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
+	  }
 	}
-	
-	rowCounter++;
-      }      
+	rowCounter++;	
+      }
     }
 
     MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);

AMDiS/src/parallel/PetscSolverFeti.h

@@ -131,7 +131,9 @@ namespace AMDiS {
     /// to \ref mat_lagrange.
     void createMatLagrange();
 
-    std::set<BoundaryObject> getAugmentedLagrange();
+    vector<vector<BoundaryObject> > getCoarseEdges();
+
+    vector<vector<BoundaryObject> > getCoarseFaces();
 
     void createMatAugmentedLagrange();
 

AMDiS/src/parallel/PetscSolverFetiDebug.cc

@@ -641,23 +641,26 @@ namespace AMDiS {
       StdMpi<vector<double> > stdMpi(MPI::COMM_WORLD);
       vector<double> sendData;
 
-      std::set<BoundaryObject> coarseSpace = feti.getAugmentedLagrange();
-      for (std::set<BoundaryObject>::iterator it = coarseSpace.begin();
+      vector<vector<BoundaryObject> > coarseSpace = feti.getCoarseFaces();
+      for (vector<vector<BoundaryObject> >::iterator it = coarseSpace.begin();
 	   it != coarseSpace.end(); ++it) {
-	if (it->subObj == FACE) {
-	  DofFace face = it->el->getFace(it->ithObj);
-	  WorldVector<double> c0 = coords[face.get<0>()];
-	  WorldVector<double> c1 = coords[face.get<1>()];
-	  WorldVector<double> c2 = coords[face.get<2>()];
-	  sendData.push_back(c0[0]);
-	  sendData.push_back(c0[1]);
-	  sendData.push_back(c0[2]);
-	  sendData.push_back(c1[0]);
-	  sendData.push_back(c1[1]);
-	  sendData.push_back(c1[2]);
-	  sendData.push_back(c2[0]);
-	  sendData.push_back(c2[1]);
-	  sendData.push_back(c2[2]);
+	for (vector<BoundaryObject>::iterator edgeIt = it->begin();
+	     edgeIt != it->end(); ++edgeIt) {
+	  if (edgeIt->subObj == FACE) {
+	    DofFace face = edgeIt->el->getFace(edgeIt->ithObj);
+	    WorldVector<double> c0 = coords[face.get<0>()];
+	    WorldVector<double> c1 = coords[face.get<1>()];
+	    WorldVector<double> c2 = coords[face.get<2>()];
+	    sendData.push_back(c0[0]);
+	    sendData.push_back(c0[1]);
+	    sendData.push_back(c0[2]);
+	    sendData.push_back(c1[0]);
+	    sendData.push_back(c1[1]);
+	    sendData.push_back(c1[2]);
+	    sendData.push_back(c2[0]);
+	    sendData.push_back(c2[1]);
+	    sendData.push_back(c2[2]);
+	  }
 	}
       }