FETI-DP matrices are assembled now.

AMDiS/src/DOFAdmin.h

@@ -262,7 +262,7 @@ namespace AMDiS {
      */
     int getDOFIndex();
 
-    /// Frees index dof. Used by Mesh::getDof()
+    /// Frees index DOF. Used by Mesh::getDof()
     void freeDofIndex(int dof);
 
     ///
@@ -287,10 +287,10 @@ namespace AMDiS {
     /// allocated size of managed vectors and matrices
     int size;
 
-    /// number of used dof indices
+    /// number of used DOF indices
     int usedCount;
 
-    /// number of FREED dof indices (NOT size-sizeUsed)
+    /// number of FREED DOF indices (NOT size - sizeUsed)
     int holeCount;
 
     /// > max. index of a used entry
@@ -302,7 +302,7 @@ namespace AMDiS {
      */
     DimVec<int> nDof;
 
-    /// Dofs from previous DOFAdmins
+    /// DOFs from previous DOFAdmins
     DimVec<int> nPreDof;
 
     /// List of all managed DOFIndexed objects.

AMDiS/src/parallel/MeshDistributor.cc

@@ -1685,42 +1685,35 @@ namespace AMDiS {
     recvDofs.clear();
 
     if (createBoundaryDofFlag.isSet(BOUNDARY_SUBOBJ_SORTED)) {
-      DofContainer dofs;
+      for (int geo = FACE; geo >= VERTEX; geo--)
+	boundaryDofInfo.geoDofs[static_cast<GeoIndex>(geo)].clear();
 
       for (int geo = FACE; geo >= VERTEX; geo--) {
-	std::set<const DegreeOfFreedom*> &dofSet =
-	  boundaryDofInfo.geoDofs[static_cast<GeoIndex>(geo)];
-	dofSet.clear();
-
 	for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
 	  if (it->rankObj.subObj == geo) {
-	    dofs.clear();
+	    DofContainer dofs;
 	    it->rankObj.el->getAllDofs(feSpace, it->rankObj, dofs);
 
 	    DofContainer& tmp = sendDofs[it.getRank()];
 	    tmp.insert(tmp.end(), dofs.begin(), dofs.end());
 
 	    if (createBoundaryDofFlag.isSet(BOUNDARY_FILL_INFO_SEND_DOFS))
-	      dofSet.insert(dofs.begin(), dofs.end());
+	      boundaryDofInfo.geoDofs[static_cast<GeoIndex>(geo)].insert(dofs.begin(), dofs.end());
 	  }
 	}
       }
 
       for (int geo = FACE; geo >= VERTEX; geo--) {
-	std::set<const DegreeOfFreedom*> &dofSet =
-	  boundaryDofInfo.geoDofs[static_cast<GeoIndex>(geo)];
-	dofSet.clear();
-
 	for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
 	  if (it->rankObj.subObj == geo) {
-	    dofs.clear();
+	    DofContainer dofs;
 	    it->rankObj.el->getAllDofs(feSpace, it->rankObj, dofs);
 
 	    DofContainer& tmp = recvDofs[it.getRank()];
 	    tmp.insert(tmp.end(), dofs.begin(), dofs.end());
 
 	    if (createBoundaryDofFlag.isSet(BOUNDARY_FILL_INFO_RECV_DOFS))
-	      dofSet.insert(dofs.begin(), dofs.end());
+	      boundaryDofInfo.geoDofs[static_cast<GeoIndex>(geo)].insert(dofs.begin(), dofs.end());
 	  }
 	}
       }

AMDiS/src/parallel/ParallelTypes.h

@@ -33,11 +33,17 @@ namespace AMDiS {
 
   using namespace std;
 
+  /// Defines set of DOF indices.
+  typedef std::set<DegreeOfFreedom> DofIndexSet;
+
   /// Defines a mapping type from DOFs to rank numbers.
   typedef map<const DegreeOfFreedom*, int> DofToRank;
   
   /// Defines a mapping type from DOFs to a set of rank numbers.
   typedef map<const DegreeOfFreedom*, std::set<int> > DofToPartitions;
+
+  /// Defines a mapping type from DOF indices to a set of rank numbers.
+  typedef map<DegreeOfFreedom, std::set<int> > DofIndexToPartitions;
   
   /// Defines a mapping type from rank numbers to sets of DOFs.
   typedef map<int, DofContainer> RankToDofContainer;

AMDiS/src/parallel/PetscSolver.h

@@ -46,7 +46,8 @@ namespace AMDiS {
   {
   public:
     PetscSolver()
-      : meshDistributor(NULL)
+      : meshDistributor(NULL),
+	mpiRank(-1)
     {}
 
     virtual ~PetscSolver() {}
@@ -54,6 +55,7 @@ namespace AMDiS {
     void setMeshDistributor(MeshDistributor *m)
     {
       meshDistributor = m;
+      mpiRank = meshDistributor->getMpiRank();
     }
 
     /** \brief
@@ -81,6 +83,8 @@ namespace AMDiS {
 
   protected:
     MeshDistributor *meshDistributor;
+    
+    int mpiRank;
 
     /// Petsc's matrix structure.
     Mat petscMatrix;

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -20,33 +20,49 @@ namespace AMDiS {
 
 
 #ifdef HAVE_PETSC_DEV 
-  void PetscSolverFeti::updateDofData(int nComponents)
+  void PetscSolverFeti::updateDofData()
   {
     FUNCNAME("PetscSolverFeti::updateDofData()");
+
+    TEST_EXIT(meshDistributor->getMesh()->getDim() == 2)
+      ("Works for 2D problems only!");
+
+    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
+      ("Works for linear basis functions only!\n");
    
-    definePrimals();
+    createPrimals();
+
+    createDuals();
+
+    createLagrange();
+
+    createIndexB();
   }
 
 
-  void PetscSolverFeti::definePrimals()
+  void PetscSolverFeti::createPrimals()
   {
-    FUNCNAME("PetscSolverFeti::definePrimals()");
+    FUNCNAME("PetscSolverFeti::createPrimals()");  
 
-    TEST_EXIT_DBG(meshDistributor->getMesh()->getDim() == 2)
-      ("Works for 2D problems only!");
+    primals.clear();
+    DofContainerSet& vertices = 
+      meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX];
+    TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n");
+    for (DofContainerSet::iterator it = vertices.begin(); 
+	 it != vertices.end(); ++it)
+      primals.insert(**it);
 
-    primals = meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX];
     globalPrimalIndex.clear();
-   
-    int nRankPrimals = 0;
-    for (DofContainerSet::iterator it = primals.begin(); 
-	 it != primals.end(); ++it)
-      if (meshDistributor->getIsRankDof(**it)) {
-	globalPrimalIndex[**it] = nRankPrimals;
+    nRankPrimals = 0;
+    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
+      if (meshDistributor->getIsRankDof(*it)) {
+	globalPrimalIndex[*it] = nRankPrimals;
 	nRankPrimals++;
       }
 
-    int nOverallPrimals = 0, rStartPrimals = 0;
+
+    nOverallPrimals = 0;
+    int rStartPrimals = 0;
     mpi::getDofNumbering(meshDistributor->getMpiComm(),
 			 nRankPrimals, rStartPrimals, nOverallPrimals);
 
@@ -56,18 +72,480 @@ namespace AMDiS {
 
     MSG_DBG("nRankPrimals = %d   nOverallPrimals = %d\n",
 	    nRankPrimals, nOverallPrimals);
+
+    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
+    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
+    for (RankToDofContainer::iterator it = sendDofs.begin();
+	 it != sendDofs.end(); ++it)
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt)
+	if (globalPrimalIndex.count(**dofIt))
+	  stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]);
+    stdMpi.updateSendDataSize();
+
+    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
+    for (RankToDofContainer::iterator it = recvDofs.begin();
+	 it != recvDofs.end(); ++it) {
+      bool recvFromRank = false;
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt)
+	if (primals.count(**dofIt) && 
+	    meshDistributor->getIsRankDof(**dofIt) == false) {
+	  recvFromRank = true;
+	  break;
+	}
+
+      if (recvFromRank) 
+	stdMpi.recv(it->first);
+    }
+    stdMpi.startCommunication();
+
+    for (RankToDofContainer::iterator it = recvDofs.begin();
+	 it != recvDofs.end(); ++it) {
+      int i = 0;
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt) {
+	if (primals.count(**dofIt) && 
+	    meshDistributor->getIsRankDof(**dofIt) == false)
+	  globalPrimalIndex[**dofIt] = stdMpi.getRecvData(it->first)[i++];
+      }
+    }
+
+    TEST_EXIT_DBG(primals.size() == globalPrimalIndex.size())
+      ("Number of primals %d, but number of global primals on this rank is %d!\n",
+       primals.size(), globalPrimalIndex.size());
+
+
+    TEST_EXIT_DBG(nOverallPrimals > 0)
+      ("There are zero primal nodes in domain!\n");
+  }
+
+
+  void PetscSolverFeti::createDuals()
+  {
+    FUNCNAME("PetscSolverFeti::createDuals()");
+    
+    // === Create for each dual node the set of ranks that contain this ===
+    // === node (denoted by W(x_j)).                                    ===
+
+    boundaryDofRanks.clear();
+
+    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
+    for (RankToDofContainer::iterator it = sendDofs.begin();
+	 it != sendDofs.end(); ++it) {
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt) {
+	// If DOF is not primal, i.e., its a dual node
+	if (primals.count(**dofIt) == 0) {
+	  boundaryDofRanks[**dofIt].insert(mpiRank);
+	  boundaryDofRanks[**dofIt].insert(it->first);
+	}
+      }
+    }
+
+    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
+    for (RankToDofContainer::iterator it = sendDofs.begin();
+	 it != sendDofs.end(); ++it)
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt)
+	if (primals.count(**dofIt) == 0)
+	  stdMpi.getSendData(it->first).push_back(boundaryDofRanks[**dofIt]);
+
+    stdMpi.updateSendDataSize();
+
+    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
+    for (RankToDofContainer::iterator it = recvDofs.begin();
+	 it != recvDofs.end(); ++it) {
+      bool recvFromRank = false;
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt)
+	if (primals.count(**dofIt) == 0) {
+	  recvFromRank = true;
+	  break;
+	}
+
+      if (recvFromRank)
+	stdMpi.recv(it->first);
+    }
+    stdMpi.startCommunication();
+
+    for (RankToDofContainer::iterator it = recvDofs.begin();
+	 it != recvDofs.end(); ++it) {
+      int i = 0;
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt)	
+	if (primals.count(**dofIt) == 0)
+	  boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++];	      
+    }
+
+
+    // === Create global index of the dual nodes on each rank. ===
+
+    duals.clear();
+    globalDualIndex.clear();
+
+    int nRankAllDofs = meshDistributor->getFeSpace()->getAdmin()->getUsedDofs();
+    nRankB = nRankAllDofs - primals.size();
+    nOverallB = 0;
+    rStartB = 0;
+    mpi::getDofNumbering(meshDistributor->getMpiComm(),
+			 nRankB, rStartB, nOverallB);
+    DofContainer allBoundaryDofs;
+    meshDistributor->getAllBoundaryDofs(allBoundaryDofs);
+    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();
+
+    int nRankDuals = 0;
+    for (DofContainer::iterator it = allBoundaryDofs.begin();
+	 it != allBoundaryDofs.end(); ++it) {
+      if (primals.count(**it) == 0) {
+	duals.insert(**it);
+	globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals;
+	nRankDuals++;
+      }
+    }
+
+    int nOverallDuals = nRankDuals;
+    mpi::globalAdd(nOverallDuals);
+
+    MSG_DBG("nRankDuals = %d   nOverallDuals = %d\n",
+	    nRankDuals, nOverallDuals);
+  }
+
+  
+  void PetscSolverFeti::createLagrange()
+  {
+    FUNCNAME("PetscSolverFeti::createLagrange()");
+
+    nRankLagrange = 0;
+    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
+      if (meshDistributor->getIsRankDof(*it)) {
+	dofFirstLagrange[*it] = nRankLagrange;
+	int degree = boundaryDofRanks[*it].size();
+	nRankLagrange += (degree * (degree - 1)) / 2;
+      }
+    }
+
+    nOverallLagrange = 0;
+    int rStartLagrange = 0;
+    mpi::getDofNumbering(meshDistributor->getMpiComm(),
+			 nRankLagrange, rStartLagrange, nOverallLagrange);
+
+    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it)
+      if (meshDistributor->getIsRankDof(*it))
+	dofFirstLagrange[*it] += rStartLagrange;
+
+    MSG_DBG("nRankLagrange = %d  nOverallLagrange = %d\n",
+	    nRankLagrange, nOverallLagrange);
+
+
+    // === ===
+
+    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
+    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
+    for (RankToDofContainer::iterator it = sendDofs.begin();
+	 it != sendDofs.end(); ++it)
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt) {
+	if (primals.count(**dofIt) == 0) {
+	  TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n");
+	  stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]);
+	}
+      }
+    stdMpi.updateSendDataSize();
+
+    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
+    for (RankToDofContainer::iterator it = recvDofs.begin();
+	 it != recvDofs.end(); ++it) {
+      bool recvData = false;
+      for (DofContainer::iterator dofIt = it->second.begin();
+	   dofIt != it->second.end(); ++dofIt)
+	if (primals.count(**dofIt) == 0) {
+	  recvData = true;
+	  break;
+	}
+	  
+      if (recvData)
+	stdMpi.recv(it->first);
+    }
+
+    stdMpi.startCommunication();
+
+    for (RankToDofContainer::iterator it = recvDofs.begin();
+	 it != recvDofs.end(); ++it) {
+      int counter = 0;
+      for (unsigned int i = 0; i < it->second.size(); i++)
+	if (primals.count(*(it->second[i])) == 0)
+	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
+    }
+     
+  }
+
+
+  void PetscSolverFeti::createIndexB()
+  {
+    FUNCNAME("PetscSolverFeti::createIndeB()");
+
+    globalIndexB.clear();
+    DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin();
+    
+    for (int i = 0; i < admin->getUsedSize(); i++)
+      if (admin->isDofFree(i) == false && primals.count(i) == 0)
+	if (duals.count(i) == 0 && primals.count(i) == 0)
+	  globalIndexB[i] = -1;
+
+    int nInterior = 0;
+    for (DofMapping::iterator it = globalIndexB.begin(); 
+	 it != globalIndexB.end(); ++it) {
+      it->second = nInterior + rStartB;
+      nInterior++;
+    }
+
+    TEST_EXIT_DBG(nInterior + primals.size() + duals.size() == 
+		  static_cast<unsigned int>(admin->getUsedDofs()))
+      ("Should not happen!\n");
+
+    for (DofIndexSet::iterator it = duals.begin();
+	 it != duals.end(); ++it)
+      globalIndexB[*it] = globalDualIndex[*it];
   }
 
 
-  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
+  void PetscSolverFeti::createMatLagrange(int nComponents)
+  {
+    FUNCNAME("PetscSolverFeti::createMatLagrange()");
+
+    MatCreateMPIAIJ(PETSC_COMM_WORLD,
+		    nRankLagrange * nComponents, nRankB * nComponents,
+		    nOverallLagrange * nComponents, nOverallB * nComponents,
+		    2, PETSC_NULL, 2, PETSC_NULL,
+		    &mat_lagrange);
+
+    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
+      TEST_EXIT_DBG(dofFirstLagrange.count(*it))("Should not happen!\n");
+      TEST_EXIT_DBG(boundaryDofRanks.count(*it))("Should not happen!\n");
+
+      int index = dofFirstLagrange[*it];
+      vector<int> W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end());
+      int degree = W.size();
+
+      TEST_EXIT_DBG(globalDualIndex.count(*it))("Should not happen!\n");
+      int dualCol = globalDualIndex[*it];
+
+      for (int i = 0; i < degree; i++) {
+	for (int j = i + 1; j < degree; j++) {
+	  if (W[i] == mpiRank || W[j] == mpiRank) {	      
+	    for (int k = 0; k < nComponents; k++) {
+	      int rowIndex = index * nComponents + k;
+	      int colIndex = dualCol * nComponents + k;
+	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
+	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
+			  INSERT_VALUES);
+	    }
+	  }
+
+	  index++;
+	}
+      }
+    }
+
+    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
+    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
+  }
+
+
+  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat, 
+					SystemVector *vec)
   {
     FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
+
+    updateDofData();
+
+    int nComponents = vec->getSize();
+
+    MatCreateMPIAIJ(PETSC_COMM_WORLD,
+		    nRankB * nComponents, nRankB * nComponents,
+		    nOverallB * nComponents, nOverallB * nComponents,
+		    100, PETSC_NULL, 100, PETSC_NULL, 
+		    &mat_b_b);
+
+    MatCreateMPIAIJ(PETSC_COMM_WORLD,
+		    nRankPrimals * nComponents, nRankPrimals * nComponents,
+		    nOverallPrimals * nComponents, nOverallPrimals * nComponents,
+		    10, PETSC_NULL, 10, PETSC_NULL, 
+		    &mat_primal_primal);
+
+    MatCreateMPIAIJ(PETSC_COMM_WORLD,
+		    nRankB * nComponents, nRankPrimals * nComponents,
+		    nOverallB * nComponents, nOverallPrimals * nComponents,
+		    100, PETSC_NULL, 100, PETSC_NULL, 
+		    &mat_b_primal);
+
+    MatCreateMPIAIJ(PETSC_COMM_WORLD,
+		    nRankPrimals * nComponents, nRankB * nComponents,
+		    nOverallPrimals * nComponents, nOverallB * nComponents,
+		    100, PETSC_NULL, 100, PETSC_NULL, 
+		    &mat_primal_b);
+
+    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
+    namespace traits = mtl::traits;
+    typedef DOFMatrix::base_matrix_type Matrix;
+
+    typedef traits::range_generator<row, Matrix>::type cursor_type;
+    typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+
+    vector<int> cols, colsOther;
+    vector<double> values, valuesOther;
+    cols.reserve(300