Some work done even today.

AMDiS/src/BoundaryObject.h

@@ -50,7 +50,7 @@ namespace AMDiS {
 				   BoundaryObject &obj1,
 				   const FiniteElemSpace *feSpace,
 				   BoundaryType boundary);
-
+    
     bool operator==(const BoundaryObject& other) const;
     
     bool operator!=(const BoundaryObject& other) const;

AMDiS/src/Debug.cc

@@ -531,7 +531,7 @@ namespace AMDiS {
       FUNCNAME("debug::writeCoordsFile()");
 
       DOFVector<WorldVector<double> > coords(feSpace, "tmp");
-      feSpace->getMesh()->getDofIndexCoords(feSpace, coords);
+      feSpace->getMesh()->getDofIndexCoords(coords);
 
       std::ofstream file;
       file.open(filename.c_str());
@@ -647,7 +647,7 @@ namespace AMDiS {
     void importDofVectorByCoords(DOFVector<double>* vec, std::string filename)
     { 
       DOFVector<WorldVector<double> > coords(vec->getFeSpace(), "dofCoords");
-      vec->getFeSpace()->getMesh()->getDofIndexCoords(vec->getFeSpace(), coords);
+      vec->getFeSpace()->getMesh()->getDofIndexCoords(coords);
       int dim = vec->getFeSpace()->getMesh()->getDim();
 
       std::ifstream file;
@@ -690,7 +690,7 @@ namespace AMDiS {
 				 std::string filename)
     {
       DOFVector<WorldVector<double> > coords(vec->getFeSpace(), "dofCoords");
-      vec->getFeSpace()->getMesh()->getDofIndexCoords(vec->getFeSpace(), coords);
+      vec->getFeSpace()->getMesh()->getDofIndexCoords(coords);
       int dim = vec->getFeSpace()->getMesh()->getDim();
 
       std::ofstream file;
@@ -817,7 +817,7 @@ namespace AMDiS {
 	 dofs0.size(), dofs1.size());
 
       DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
-      feSpace->getMesh()->getDofIndexCoords(feSpace, coords);
+      feSpace->getMesh()->getDofIndexCoords(coords);
 
       for (unsigned int i = 0; i < dofs0.size(); i++) {
 	WorldVector<double> tmp = coords[*(dofs0[i])];

AMDiS/src/Mesh.cc

@@ -943,11 +943,11 @@ namespace AMDiS {
   }
 
 
-  void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace,
-			       DOFVector<WorldVector<double> >& coords)
+  void Mesh::getDofIndexCoords(DOFVector<WorldVector<double> >& coords)
   {
     FUNCNAME("Mesh::getDofIndexCoords()");
 
+    const FiniteElemSpace *feSpace = coords.getFeSpace();
     const BasisFunction* basFcts = feSpace->getBasisFcts();
     int nBasFcts = basFcts->getNumber();
     vector<DegreeOfFreedom> dofVec(nBasFcts);

AMDiS/src/Mesh.h

@@ -487,15 +487,13 @@ namespace AMDiS {
 			   WorldVector<double>& coords);
 
     /** \brief
-     * Traverse the whole mesh and stores to each DOF of the given finite
-     * element space the coordinates in a given DOFVector. Works in the same
-     * way as the function \ref getDofIndexCoords defined above.
+     * Traverse the whole mesh and stores to each DOF the coordinates in a given
+     * DOFVector. Works in the same way as the function \ref getDofIndexCoords 
+     * defined above.
      *
-     * @param[in]   feSpace   The FE space to be used for the search.
      * @param[out]  coords    DOF vector that stores the coordinates to each DOF.
      */
-    void getDofIndexCoords(const FiniteElemSpace* feSpace,
-			   DOFVector<WorldVector<double> >& coords);
+    void getDofIndexCoords(DOFVector<WorldVector<double> >& coords);
 
     /** \brief
      * Traverse the mesh and get all DOFs in this mesh for a given FE space.

AMDiS/src/est/ResidualEstimator.cc

@@ -169,7 +169,7 @@ namespace AMDiS {
       return;
 
     DOFVector<WorldVector<double> > coords(uh[0]->getFeSpace(), "tmp");
-    mesh->getDofIndexCoords(uh[0]->getFeSpace(), coords);
+    mesh->getDofIndexCoords(coords);
 
     InteriorBoundary &intBoundary = 
       MeshDistributor::globalMeshDistributor->getIntBoundary();

AMDiS/src/parallel/BddcMlSolver.cc

@@ -179,7 +179,7 @@ namespace AMDiS {
     
     {
       DOFVector<WorldVector<double> > coordDofs(feSpace, "tmp");
-      mesh->getDofIndexCoords(feSpace, coordDofs);
+      mesh->getDofIndexCoords(coordDofs);
 
       for (int i = 0; i < lxyz2; i++)
 	for (int j = 0; j < nnods; j++)

AMDiS/src/parallel/ElementObjectDatabase.h

@@ -345,6 +345,24 @@ namespace AMDiS {
       return faceInRank[face];
     }
 
+    /// Get degree of a DOF, thus the number of ranks which contain it.
+    inline int getDegree(DegreeOfFreedom dof)
+    {
+      return vertexInRank[dof].size();
+    }
+
+    /// Get degree of an edge, thus the number of ranks which contain it.
+    inline int getDegree(DofEdge edge)
+    {
+      return edgeInRank[edge].size();
+    }
+
+    /// Get degree of a face, thus the number of ranks which contain it.
+    inline int getDegree(DofFace face)
+    {
+      return faceInRank[face].size();
+    }
+
     /// Returns to an element object data the appropriate vertex DOF.
     DegreeOfFreedom getVertexLocalMap(ElementObjectData &data)
     {

AMDiS/src/parallel/InteriorBoundary.cc

@@ -532,6 +532,24 @@ namespace AMDiS {
   }
 
 
+  int InteriorBoundary::getDegreeOwn(BoundaryObject &bObj)
+  {
+    int counter = 0;
+
+    for (RankToBoundMap::iterator it = own.begin(); it != own.end(); ++it) {
+      for (vector<AtomicBoundary>::iterator bIt = it->second.begin(); 
+	   bIt != it->second.end(); ++bIt) {
+	if (bIt->rankObj == bObj) {
+	  counter++;
+	  break;
+	}
+      }
+    }
+
+    return counter;
+  }
+
+
   void InteriorBoundary::serialize(ostream &out)
   {
     serialize(out, own);

AMDiS/src/parallel/InteriorBoundary.h

@@ -79,6 +79,8 @@ namespace AMDiS {
       return static_cast<bool>(periodic.size());
     }
 
+    int getDegreeOwn(BoundaryObject &bObj);
+
     /// Writes this object to a file.
     void serialize(ostream &out);
 

AMDiS/src/parallel/MeshDistributor.cc

@@ -238,7 +238,7 @@ namespace AMDiS {
 #else
       int writePartMesh = 0;
 #endif
-      Parameters::get("dbg->write part mesh", writePartMesh);
+      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
 
       if (writePartMesh > 0) {
 	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");

AMDiS/src/parallel/MeshDistributor.h

@@ -309,7 +309,7 @@ namespace AMDiS {
 			    int level,
 			    DofContainer& dofs);
 
-    const ElementObjectDatabase& getElementObjectDb() 
+    ElementObjectDatabase& getElementObjectDb() 
     {
       return elObjDb;
     }

AMDiS/src/parallel/MeshPartitioner.cc

@@ -55,12 +55,6 @@ namespace AMDiS {
 	// rank number.
 	
 	int elInRank = std::min(elIndex / elPerRank, mpiSize - 1);
-//  	int elInRank = 0;
-//  	if (elIndex <= 11 || elIndex >= 28)
-//  	  elInRank = 0;
-//  	else
-//  	  elInRank = 1;
-
 	elementInRank[elIndex] = (elInRank == mpiRank);
 	partitionMap[elIndex] = elInRank;	
       } else {

AMDiS/src/parallel/ParallelDebug.cc

@@ -353,7 +353,7 @@ namespace AMDiS {
     RankToCoords recvCoords;
 
     DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
-    pdb.mesh->getDofIndexCoords(feSpace, coords);
+    pdb.mesh->getDofIndexCoords(coords);
 
     for (DofComm::Iterator it(pdb.dofComm.getSendDofs(), feSpace); 
 	 !it.end(); it.nextRank())
@@ -474,7 +474,7 @@ namespace AMDiS {
     const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1];
 
     DOFVector<WorldVector<double> > coords(feSpace, "tmp");
-    pdb.mesh->getDofIndexCoords(feSpace, coords);
+    pdb.mesh->getDofIndexCoords(coords);
 
     typedef map<WorldVector<double>, int> CoordsIndexMap;
     CoordsIndexMap coordsToIndex;
@@ -772,7 +772,7 @@ namespace AMDiS {
     filename << prefix << "-" << MPI::COMM_WORLD.Get_rank() << "." << postfix;
 
     DOFVector<WorldVector<double> > coords(feSpace, "tmp");
-    mesh->getDofIndexCoords(feSpace, coords);
+    mesh->getDofIndexCoords(coords);
 
     typedef map<int, vector<DegreeOfFreedom> > ElDofMap;
     ElDofMap elDofMap;

AMDiS/src/parallel/ParallelTypes.h

@@ -73,7 +73,6 @@ namespace AMDiS {
     int local, global;
   };
 
-  //  typedef std::map<DegreeOfFreedom, MultiIndex> DofMap;
   typedef boost::container::flat_map<DegreeOfFreedom, MultiIndex> DofMap;
  
   typedef vector<MeshStructure> MeshCodeVec;

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -279,6 +279,12 @@ namespace AMDiS {
       MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
 	  timeCounter);
     }
+
+
+    bool writePrimals = false;
+    Parameters::get("parallel->debug->write primals", writePrimals);
+    if (writePrimals)
+      PetscSolverFetiDebug::writePrimalFiles(*this);
   }
 
 
@@ -709,6 +715,9 @@ namespace AMDiS {
 
   bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
   {
+    if (meshDistributor->getIntBoundary().getDegreeOwn(edge) != 3)
+      return false;
+
     Element *el = edge.el;
     int i0 = el->getVertexOfEdge(edge.ithObj, 0);
     int i1 = el->getVertexOfEdge(edge.ithObj, 1);
@@ -732,6 +741,9 @@ namespace AMDiS {
     if (!augmentedLagrange)
       return;
 
+    int nEdges = 0;
+    int nFaces = 0;
+
     double wtime = MPI::Wtime();
 
     nOverallEdges = 0;
@@ -755,8 +767,15 @@ namespace AMDiS {
 	  }
 	}
 
-	if (!dirichletOnlyEdge)
+	if (!dirichletOnlyEdge) {
+	  if (it->rankObj.subObj == EDGE)
+	    nEdges++;
+	  
+	  if (it->rankObj.subObj == FACE)
+	    nFaces++;
+
 	  allEdges.insert(it->rankObj);
+	}
       }
     }
 
@@ -765,6 +784,11 @@ namespace AMDiS {
     mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);
 
     MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);
+
+    mpi::globalAdd(nEdges);
+    mpi::globalAdd(nFaces);
+
+    MSG("NEDGES = %d  NFACES = %d\n", nEdges, nFaces);
     
     nRankEdges *= componentSpaces.size();
     rStartEdges *= componentSpaces.size();
@@ -884,7 +908,7 @@ namespace AMDiS {
 	createMatExplicitAugmentedSchurPrimal();
 
 
-      // === Create KSP solver object and set appropriate solver options. ====
+      // === Create KSP solver object and set appropriate solver options. ===
 
       KSPCreate(mpiCommGlobal, &ksp_schur_primal);
       KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
@@ -1005,9 +1029,7 @@ namespace AMDiS {
       
       // matTmp = inv(A_BB) trans(J) trans(Q)
       Mat qT, jTqT;
-      MatTranspose(mat_augmented_lagrange, MAT_INITIAL_MATRIX, &qT);
-      //      Mat jT;
-      //      MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &jT);
+      MatTranspose(mat_augmented_lagrange, MAT_INITIAL_MATRIX, &qT);     
       MatTransposeMatMult(mat_lagrange, qT, MAT_INITIAL_MATRIX, PETSC_DEFAULT, 
 			  &jTqT);
       petsc_helper::blockMatMatSolve(subdomain->getSolver(), jTqT, matTmp);

AMDiS/src/parallel/PetscSolverFetiDebug.cc

@@ -565,4 +565,82 @@ namespace AMDiS {
 
     return nZeroRows;
   }
+
+
+  void PetscSolverFetiDebug::writePrimalFiles(PetscSolverFeti &feti)
+  {
+    FUNCNAME("PetscSolverFetiDebug::writePrimalFiles()");
+
+    Mesh *mesh = feti.componentSpaces[0]->getMesh();
+    DOFVector<WorldVector<double> > coords(feti.componentSpaces[0], "coords");
+    mesh->getDofIndexCoords(coords);
+
+
+    // === Write vertex primals. ===
+
+    {
+      StdMpi<vector<double> > stdMpi(MPI::COMM_WORLD);
+      vector<double> sendData;
+      
+      DofMap &primals = feti.primalDofMap[0].getMap();
+      for (DofMap::iterator it = primals.begin();it != primals.end(); ++it) {
+	if (feti.primalDofMap[0].isRankDof(it->first)) {
+	  for (int i = 0; i < mesh->getDim(); i++)
+	    sendData.push_back(coords[it->first][i]);
+	}
+      }
+      
+      if (MPI::COMM_WORLD.Get_rank() == 0) {
+	for (int i = 1; i < MPI::COMM_WORLD.Get_size(); i++)
+	  stdMpi.recv(i);
+      } else {
+	stdMpi.send(0, sendData);
+      }
+
+      stdMpi.startCommunication();
+      
+      if (MPI::COMM_WORLD.Get_rank() == 0) {
+	vector<WorldVector<double> > allPrimals;
+	
+	for (int i = 1; i < MPI::COMM_WORLD.Get_size(); i++) {
+	  vector<double> &recvData = stdMpi.getRecvData(i);
+	  
+	  TEST_EXIT(recvData.size() % mesh->getDim() == 0)
+	    ("Wrong number of coordinates received!\n");
+	  
+	  int nCoords = recvData.size() / mesh->getDim();
+	  int counter = 0;
+	  for (int j = 0; j < nCoords; j++) {
+	    WorldVector<double> c;
+	    for (int k = 0; k < mesh->getDim(); k++)
+	      c[k] = recvData[counter++];
+	    allPrimals.push_back(c);
+	  }
+	}
+	
+	ofstream file;
+	file.open("primals.xyz");
+	file << allPrimals.size() << "\n";
+	file << "Primals\n";
+	for (int i = 0; i < static_cast<int>(allPrimals.size()); i++) {
+	  file << "P ";
+	  if (mesh->getDim() == 2)
+	    file << allPrimals[i][0] << " " << allPrimals[i][1] << " 0.0";
+	  else
+	    file << allPrimals[i][0] << " " << allPrimals[i][1] << " " << allPrimals[i][2];
+	  file << "\n";
+	}
+	file.close();
+      }
+    }
+
+         
+    // === Write face primals. ===
+    
+    {
+      
+    }
+      
+  }
+  
 }

AMDiS/src/parallel/PetscSolverFetiDebug.h

@@ -137,6 +137,9 @@ namespace AMDiS {
      * \return     int    Number of zero rows in matrix;
      */
     static int testZeroRows(Mat mat);
+
+    /// Write files with information about the primal coarner nodes.
+    static void writePrimalFiles(PetscSolverFeti &feti);
   };
 
 }