Added debuging functionality for face coarse space in FETI-DP solver.

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -734,19 +734,13 @@ namespace AMDiS {
   }
 
 
-  void PetscSolverFeti::createMatAugmentedLagrange()
+  std::set<BoundaryObject> PetscSolverFeti::getAugmentedLagrange()
   {
-    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");
-
-    if (!augmentedLagrange)
-      return;
+    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
 
     int nEdges = 0;
     int nFaces = 0;
 
-    double wtime = MPI::Wtime();
-
-    nOverallEdges = 0;
     InteriorBoundary &intBound = meshDistributor->getIntBoundary();
     std::set<BoundaryObject> allEdges;
     for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
@@ -779,17 +773,31 @@ namespace AMDiS {
       }
     }
 
+    mpi::globalAdd(nEdges);
+    mpi::globalAdd(nFaces);
+
+    MSG("NEDGES = %d  NFACES = %d\n", nEdges, nFaces);   
+
+    return allEdges;
+  }
+
+
+  void PetscSolverFeti::createMatAugmentedLagrange()
+  {
+    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");
+
+    if (!augmentedLagrange)
+      return;
+
+    double wtime = MPI::Wtime();
+
+    std::set<BoundaryObject> allEdges = getAugmentedLagrange();
     nRankEdges = allEdges.size();    
     int rStartEdges = 0;
     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();
     nOverallEdges *= componentSpaces.size();

AMDiS/src/parallel/PetscSolverFeti.h

@@ -131,6 +131,8 @@ namespace AMDiS {
     /// to \ref mat_lagrange.
     void createMatLagrange();
 
+    std::set<BoundaryObject> getAugmentedLagrange();
+
     void createMatAugmentedLagrange();
 
     bool testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace);

AMDiS/src/parallel/PetscSolverFetiDebug.cc

@@ -636,9 +636,125 @@ namespace AMDiS {
 
          
     // === Write face primals. ===
-    
-    {
-      
+
+    if (mesh->getDim() == 3) {
+      StdMpi<vector<double> > stdMpi(MPI::COMM_WORLD);
+      vector<double> sendData;
+
+      std::set<BoundaryObject> coarseSpace = feti.getAugmentedLagrange();
+      for (std::set<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]);
+	}
+      }
+
+      if (MPI::COMM_WORLD.Get_rank() == 0) {
+	TEST_EXIT(sendData.size() == 0)("Should not happen!\n");
+      }
+
+      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> > faceNodes;
+	
+	for (int i = 1; i < MPI::COMM_WORLD.Get_size(); i++) {
+	  vector<double> &recvData = stdMpi.getRecvData(i);
+	  
+	  TEST_EXIT(recvData.size() % 9 == 0)
+	    ("Wrong number of coordinates received!\n");
+
+	  int nNodes = recvData.size() / 3;
+	  int counter = 0;
+	  for (int j = 0; j < nNodes; j++) {
+	    WorldVector<double> c;
+	    c[0] = recvData[counter++];
+	    c[1] = recvData[counter++];
+	    c[2] = recvData[counter++];
+	    faceNodes.push_back(c);
+	  }
+	}
+
+	TEST_EXIT(faceNodes.size() % 3 == 0)("Should not happen!\n");
+
+	int nElements = faceNodes.size() / 3;
+
+	ofstream file;
+	file.open("faces.vtu");
+
+	file << "<?xml version=\"1.0\"?>\n";
+	file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
+	file << "  <UnstructuredGrid>\n";
+	file << "    <Piece NumberOfPoints=\"" << faceNodes.size()
+	     << "\" NumberOfCells=\"" << nElements << "\">\n";
+	file << "      <Points>\n";
+	file << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
+
+	for (int i = 0; i < static_cast<int>(faceNodes.size()); i++) 
+	  file << " " << faceNodes[i][0] << " " << faceNodes[i][1] << " " << faceNodes[i][2] << "\n";
+	
+	file << "        </DataArray>\n";
+	file << "      </Points>\n";
+	file << "      <Cells>\n";
+	file << "        <DataArray type=\"Int32\" Name=\"offsets\">\n";
+
+	for (int i = 0; i < nElements; i++)
+	  file << " " << (i + 1) * 3 << "\n";
+	
+	file << "        </DataArray>\n";
+	file << "        <DataArray type=\"UInt8\" Name=\"types\">\n";
+	
+	for (int i = 0; i < nElements; i++)
+	  file << " 5\n";
+	
+	file << "        </DataArray>\n";
+	file << "        <DataArray type=\"Int32\" Name=\"connectivity\">\n";
+	
+	{
+	  int counter = 0;
+	  for (int i = 0; i < nElements; i++) {
+	    file << " " << counter << " " << counter + 1 << " " << counter + 2 << "\n";
+	    counter += 3;
+	  }
+	}
+	
+	file << "        </DataArray>\n";    
+	file << "      </Cells>\n";
+	file << "      <PointData>\n";
+	
+	file << "        <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">\n";
+
+	for (int i = 0; i < static_cast<int>(faceNodes.size()); i++) 
+	  file << " 0.0\n";
+	
+	file << "        </DataArray>\n";
+	
+	file << "      </PointData>\n";
+	file << "    </Piece>\n";
+	file << "  </UnstructuredGrid>\n";
+	file << "</VTKFile>\n";
+	
+	file.close();
+      }
     }
       
   }