diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 00fc57d20ee37ed3bd4285e3c8b5bf56b856077f..7190ff4dca1d5d2e289bd446a2ba115b2cc208e1 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/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();
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 7b73c2e58c0d595fc964d41f17963e72b1032149..9e308d8dc34e2d5db69a0443d19ff50f15e4eb79 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/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);
diff --git a/AMDiS/src/parallel/PetscSolverFetiDebug.cc b/AMDiS/src/parallel/PetscSolverFetiDebug.cc
index 09db80e7389f75ece02113ac751199c59dae1bd8..dc7ec4953728212829a747bcbb10963ab80bbf04 100644
--- a/AMDiS/src/parallel/PetscSolverFetiDebug.cc
+++ b/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();
+ }
}
}