Work on FETI-DP for Stokes like problems.

AMDiS/src/io/MacroReader.cc

@@ -61,12 +61,12 @@ namespace AMDiS {
       int *verticesEl1 = new int[dim];
       int *verticesEl2 = new int[dim];
       int mode = -1; // 0: drop dofs, 1: associate dofs
-      int result;
+      int result = 0;
       BoundaryType boundaryType;
       MacroReader::PeriodicMap periodicMap;
 
-      fscanf(file, "%*s %d", &n);
-      fscanf(file, "%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s");
+      result = fscanf(file, "%*s %d", &n);
+      result = fscanf(file, "%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s");
    
       for (int i = 0; i < n; i++) {
 	std::map<int, int> vertexMapEl1;

AMDiS/src/parallel/MeshDistributor.cc

@@ -1692,13 +1692,13 @@ namespace AMDiS {
     int nLevels = levelData.getLevelNumber();
     TEST_EXIT_DBG(nLevels >= 1)("Should not happen!\n");
 
-    dofMap.init(levelData, feSpaces, feSpaces, true, true);
+    dofMap.init(levelData, feSpaces, feSpaces);
     dofMap.setMpiComm(levelData.getMpiComm(0), 0);
     dofMap.setDofComm(dofComm);
     dofMap.clear();
 
     if (nLevels > 1) {
-      dofMapSd.init(levelData, feSpaces, feSpaces, true, true);
+      dofMapSd.init(levelData, feSpaces, feSpaces);
       dofMapSd.setMpiComm(levelData.getMpiComm(1), 1);
       dofMapSd.setDofComm(dofCommSd);
       dofMapSd.clear();

AMDiS/src/parallel/ParallelDofMapping.cc

@@ -71,7 +71,7 @@ namespace AMDiS {
     if (needGlobalMapping) {
       computeGlobalMapping();
       
-      if (hasNonLocalDofs)
+      if (isNonLocal)
 	computeNonLocalIndices();
     }
   }
@@ -161,15 +161,14 @@ namespace AMDiS {
   void ParallelDofMapping::init(MeshLevelData &ldata,
 				vector<const FiniteElemSpace*> &fe,
 				vector<const FiniteElemSpace*> &uniqueFe,
-				bool needGlobalMapping,
-				bool bNonLocalDofs)
+				bool b)
   {
     FUNCNAME("ParallelDofMapping::init()");
 
     levelData = &ldata;
     feSpaces = fe;
     feSpacesUnique = uniqueFe;
-    hasNonLocalDofs = bNonLocalDofs;
+    isNonLocal = b;
     
 
     // === Init the mapping for all different FE spaces. ===
@@ -177,8 +176,8 @@ namespace AMDiS {
     for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
 	 it != feSpacesUnique.end(); ++it) {
       addFeSpace(*it);
-      data[*it].setNeedGlobalMapping(needGlobalMapping);
-      data[*it].setNonLocalDofs(hasNonLocalDofs);
+      data[*it].setNeedGlobalMapping(isNonLocal);
+      data[*it].setNonLocal(isNonLocal);
     }
   }
 
@@ -310,7 +309,6 @@ namespace AMDiS {
     rStartDofs = computeStartDofs();
     
     // And finally, compute the matrix indices.
-    if (needMatIndex)
     computeMatIndex(needMatIndexFromGlobal);
   }
   
@@ -363,7 +361,7 @@ namespace AMDiS {
       offset += data[feSpaces[i]].nRankDofs;
 	
       // If there are no non local DOFs, continue with the next FE space.
-      if (!hasNonLocalDofs)
+      if (!isNonLocal)
 	continue;
       
       TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");

AMDiS/src/parallel/ParallelDofMapping.h

@@ -115,7 +115,7 @@ namespace AMDiS {
 	dofComm(NULL),
 	feSpace(NULL),
 	needGlobalMapping(false),
-	hasNonLocalDofs(false)
+	isNonLocal(false)
     {
       clear();
     }
@@ -230,9 +230,9 @@ namespace AMDiS {
 
     /// Informs the mapping whether the mapping will include DOFs that are not
     /// owned by the rank.
-    void setNonLocalDofs(bool b)
+    void setNonLocal(bool b)
     {
-      hasNonLocalDofs = b;
+      isNonLocal = b;
     }
 
     /// Informs the mapping whether a global index must be computed.
@@ -287,7 +287,7 @@ namespace AMDiS {
     /// Is true if there are DOFs in at least one subdomain that are not owned
     /// by the rank. If the value is false, each rank contains only DOFs that
     /// are also owned by this rank.
-    bool hasNonLocalDofs;
+    bool isNonLocal;
 
   public:
     /// 
@@ -306,8 +306,7 @@ namespace AMDiS {
     ParallelDofMapping() 
       : levelData(NULL),
 	dofComm(NULL),
-	hasNonLocalDofs(false),
-	needMatIndex(false),
+	isNonLocal(true),
 	needMatIndexFromGlobal(false),
 	feSpaces(0),
 	nRankDofs(1),
@@ -329,16 +328,24 @@ namespace AMDiS {
      * \param[in]  uniqueFe           Unique list of FE spaces. Thus, two
      *                                arbitrary elements of this list are always
      *                                different.
-     * \param[in]  needGlobalMapping  If true, the mapping computes also a global
-     *                                index for the DOFs.
-     * \param[in]  bNonLocalDofs      If true, at least one rank's mapping con-
+     * \param[in]  isNonLocal         If true, at least one rank's mapping con-
      *                                taines DOFs that are not owend by the rank.
      */
     void init(MeshLevelData& mld,
 	      vector<const FiniteElemSpace*> &fe,
 	      vector<const FiniteElemSpace*> &uniqueFe,
-	      bool needGlobalMapping,
-	      bool bNonLocalDofs);
+	      bool isNonLocal = true);
+
+    /// In the case of having only one FE space, this init function can be used.
+    void init(MeshLevelData& mld,
+	      const FiniteElemSpace *feSpace,
+	      bool isNonLocal = true)
+    {
+      vector<const FiniteElemSpace*> feSpaces;
+      feSpaces.push_back(feSpace);
+      init(mld, feSpaces, feSpaces, isNonLocal);
+    }
+
 
     void setMpiComm(MPI::Intracomm &m, int l);
     
@@ -359,9 +366,10 @@ namespace AMDiS {
       return *dofComm;
     }
 
-    void setComputeMatIndex(bool b, bool global = false)
+    /// Changes the computation of matrix indices based of either local or
+    /// global DOF indices, see \ref needMatIndexFromGlobal
+    void setComputeMatIndex(bool global)
     {
-      needMatIndex = b;
       needMatIndexFromGlobal = global;
     }
 
@@ -494,11 +502,7 @@ namespace AMDiS {
     /// Is true if there are DOFs in at least one subdomain that are not owned
     /// by the rank. If the value is false, each rank contains only DOFs that
     /// are also owned by this rank.
-    bool hasNonLocalDofs;
-
-    /// If true, matrix indeces for the stored DOFs are computed, see
-    /// \ref computeMatIndex.
-    bool needMatIndex;
+    bool isNonLocal;
 
     /// If matrix indices should be computed, this variable defines if the
     /// mapping from DOF indices to matrix row indices is defined on local

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -263,6 +263,7 @@ namespace AMDiS {
       ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1);
 
     MeshLevelData& levelData = meshDistributor->getMeshLevelData();
+    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
 
     if (subdomain == NULL) {
       subdomain = new PetscSolverGlobalMatrix();
@@ -278,42 +279,19 @@ namespace AMDiS {
       }
     }
 
-    primalDofMap.init(levelData,
-		      feSpaces, meshDistributor->getFeSpaces(), 
-		      true, true);
-    primalDofMap.setComputeMatIndex(true);
-    
-    dualDofMap.init(levelData, 
-		    feSpaces, meshDistributor->getFeSpaces(),
-		    false, false);
-    dualDofMap.setComputeMatIndex(true);
-
-    if (meshLevel == 0) {
-      localDofMap.init(levelData, 
-		       feSpaces, meshDistributor->getFeSpaces(),
-		       false, false);
-      localDofMap.setComputeMatIndex(true);
-    } else {
-      localDofMap.init(levelData, 
-		       feSpaces, meshDistributor->getFeSpaces(),
-		       true, true);
-      localDofMap.setComputeMatIndex(true);
-    }
-
-    lagrangeMap.init(levelData, 
-		     feSpaces, meshDistributor->getFeSpaces(),
-		     true, true);
-    lagrangeMap.setComputeMatIndex(true);
+    primalDofMap.init(levelData, feSpaces, uniqueFe);   
+    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
+    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
+    lagrangeMap.init(levelData, feSpaces, uniqueFe);
 
+    if (fullInterface != NULL)
+      interfaceDofMap.init(levelData, fullInterface);
 
     if (fetiPreconditioner != FETI_NONE) {
       TEST_EXIT(meshLevel == 0)
 	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");
 
-      interiorDofMap.init(levelData, 
-			  feSpaces,  meshDistributor->getFeSpaces(),
-			  false, false);
-      interiorDofMap.setComputeMatIndex(true);
+      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
     }
   }
 
@@ -352,6 +330,12 @@ namespace AMDiS {
     else
       localDofMap.setDofComm(meshDistributor->getDofCommSd());
 
+    if (fullInterface != NULL) {
+      interfaceDofMap.clear();
+      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
+      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
+    }
+
     for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
       const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
 
@@ -359,6 +343,8 @@ namespace AMDiS {
 
       createDuals(feSpace);
 
+      createInterfaceNodes(feSpace);
+
       createIndexB(feSpace);
     }
 
@@ -368,6 +354,9 @@ namespace AMDiS {
     if (fetiPreconditioner != FETI_NONE)
       interiorDofMap.update();
 
+    if (fullInterface != NULL)
+      interfaceDofMap.update();
+
     for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
       const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
       createLagrange(feSpace);
@@ -435,6 +424,9 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverFeti::createPrimals()");  
 
+    if (feSpace == fullInterface)
+      return;
+
     // === Define all vertices on the interior boundaries of the macro mesh ===
     // === to be primal variables.                                          ===
 
@@ -477,6 +469,9 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverFeti::createDuals()");
 
+    if (feSpace == fullInterface)
+      return;
+
     // === Create global index of the dual nodes on each rank. ===
 
     DofContainer allBoundaryDofs;
@@ -494,10 +489,32 @@ namespace AMDiS {
   }
 
   
+  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
+  {
+    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");
+
+    if (feSpace != fullInterface)
+      return;
+
+    DofContainer allBoundaryDofs;
+    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
+
+    for (DofContainer::iterator it = allBoundaryDofs.begin();
+	 it != allBoundaryDofs.end(); ++it)
+      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
+	interfaceDofMap[feSpace].insertRankDof(*it);
+      else
+	interfaceDofMap[feSpace].insertNonRankDof(*it);
+  }
+
+
   void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
   {
     FUNCNAME("PetscSolverFeti::createLagrange()");
 
+    if (feSpace == fullInterface)
+      return;
+
     boundaryDofRanks[feSpace].clear();
 
     // Stores for all rank owned communication DOFs, if the counterpart is
@@ -645,7 +662,8 @@ namespace AMDiS {
     for (int i = 0; i < admin->getUsedSize(); i++) {
       if (admin->isDofFree(i) == false && 
 	  isPrimal(feSpace, i) == false &&
-	  isDual(feSpace, i) == false) {
+	  isDual(feSpace, i) == false &&
+	  isInterface(feSpace, i) == false) {
 	if (meshLevel == 0) {
 	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
 

AMDiS/src/parallel/PetscSolverFeti.h

@@ -108,9 +108,12 @@ namespace AMDiS {
     void createPrimals(const FiniteElemSpace *feSpace);
 
     /// Defines the set of dual variables and creates the global index of
-    // dual variables.
+    /// dual variables.
     void createDuals(const FiniteElemSpace *feSpace);
 
+    /// 
+    void createInterfaceNodes(const FiniteElemSpace *feSpace);
+
     /// Create Lagrange multiplier variables corresponding to the dual 
     /// variables.
     void createLagrange(const FiniteElemSpace *feSpace);
@@ -192,6 +195,13 @@ namespace AMDiS {
       return dualDofMap[feSpace].isSet(dof);
     }
 
+    /// Checks whether a given DOF in a give FE space is an interface DOF.
+    inline bool isInterface(const FiniteElemSpace *feSpace, 
+			    DegreeOfFreedom dof)
+    {
+      return interfaceDofMap[feSpace].isSet(dof);
+    }
+
   protected:
     /// Mapping from primal DOF indices to a global index of primals.
     ParallelDofMapping primalDofMap;
@@ -199,6 +209,11 @@ namespace AMDiS {
     /// Mapping from dual DOF indices to a global index of duals.
     ParallelDofMapping dualDofMap;
 
+    /// Mapping from interface DOF indices to a global index of interface 
+    /// nodes. This is mainly used for Stokes-like solvers, where the pressure
+    /// interface nodes are neither primal nor dual.
+    ParallelDofMapping interfaceDofMap;
+
     /// Index for each non primal DOF to the global index of B variables (thus,
     /// all pure local variables).
     ParallelDofMapping localDofMap;
@@ -271,6 +286,8 @@ namespace AMDiS {
 
     int nGlobalOverallInterior;
 
+    const FiniteElemSpace* fullInterface;
+
     bool printTimings;
   };
 

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -598,7 +598,7 @@ namespace AMDiS {
     if (!nnzInterior.dnnz || recvAllValues != 0 || alwaysCreateNnzStructure) {
       vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
       
-      interiorMap->setComputeMatIndex(true, !localMatrix);
+      interiorMap->setComputeMatIndex(!localMatrix);
       interiorMap->update(feSpaces);
 
       nnzInterior.clear();