diff --git a/AMDiS/src/BoundaryObject.h b/AMDiS/src/BoundaryObject.h
index ec273251d1dfdf47685900332d343a87fbe237b4..cc9edbc37b14cef95d59e6b588872c35793fd393 100644
--- a/AMDiS/src/BoundaryObject.h
+++ b/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;
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
index 725c2602e543e44702a5b4a2f993c60c37c56ee8..c9b08ce3069af739a86f8a59caf1f236cb2e6538 100644
--- a/AMDiS/src/Debug.cc
+++ b/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])];
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index a1ad385f81d9a42567f9039b29a824dd14f11bcb..5042fbdd9a166f0e4c78641f2ad1249cf50df9f2 100644
--- a/AMDiS/src/Mesh.cc
+++ b/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);
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index c44a3a6923a8b5c8700e9ba57015bf2367c6595e..be2eff48b386f1db863cb2709578df461c94ab52 100644
--- a/AMDiS/src/Mesh.h
+++ b/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.
diff --git a/AMDiS/src/est/ResidualEstimator.cc b/AMDiS/src/est/ResidualEstimator.cc
index 5e6f4a8314d2b220e99251942054332eb9b2e99c..c1995787c343189ca66f6ba357ed7d117ade54c4 100644
--- a/AMDiS/src/est/ResidualEstimator.cc
+++ b/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();
diff --git a/AMDiS/src/parallel/BddcMlSolver.cc b/AMDiS/src/parallel/BddcMlSolver.cc
index d08605d30d52f2d057ad0c77bb1fbe627b7e3bd0..841897d5719d0143992d23c04cf6538b15adebb5 100644
--- a/AMDiS/src/parallel/BddcMlSolver.cc
+++ b/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++)
diff --git a/AMDiS/src/parallel/ElementObjectDatabase.h b/AMDiS/src/parallel/ElementObjectDatabase.h
index ad8513995f9d2a76e523d8c22c18efca253e0ec9..523d78daa897d4565c1874b554b606f2685d95c6 100644
--- a/AMDiS/src/parallel/ElementObjectDatabase.h
+++ b/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)
{
diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc
index 8eee4cd50aa51db9693ada8ca137f2d213df952d..cdfcbcec13861ae62306ec9210cf86bfaadb48b8 100644
--- a/AMDiS/src/parallel/InteriorBoundary.cc
+++ b/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);
diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h
index 594c295eb8cf21fd245203c1d3ed7de6288bf547..d763459c3c0311ef734d944b77157847a67f728e 100644
--- a/AMDiS/src/parallel/InteriorBoundary.h
+++ b/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);
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 0b416addc2ad635b8593264fec1d21f4646cbb3a..84d5f278b076979c23cc6a9cab098e45d7fadfed 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/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");
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index 3022c26cc316c5ce651cf52b59b5563d80c71c3d..49d7ed092f4daf3d16ab4345663c8e7cdb4773fe 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -309,7 +309,7 @@ namespace AMDiS {
int level,
DofContainer& dofs);
- const ElementObjectDatabase& getElementObjectDb()
+ ElementObjectDatabase& getElementObjectDb()
{
return elObjDb;
}
diff --git a/AMDiS/src/parallel/MeshPartitioner.cc b/AMDiS/src/parallel/MeshPartitioner.cc
index a873e0bc1e8b1a8a3d5f784d1f475374ab276aeb..a501eed52b1b284438b5bf3d76faa14a1bfaec7f 100644
--- a/AMDiS/src/parallel/MeshPartitioner.cc
+++ b/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 {
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index 241ddb148e6492893e7fd417ba39dec6fa9a6900..bc8b88bbf36e646850734460e2230c2892745ddb 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/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;
diff --git a/AMDiS/src/parallel/ParallelTypes.h b/AMDiS/src/parallel/ParallelTypes.h
index b0d796aee88e60ee5dc773d69ad66fea2259fdb0..9f246f29f22431e41e156434ad8ec92a8b95db5c 100644
--- a/AMDiS/src/parallel/ParallelTypes.h
+++ b/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;
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 9379f352489d205f95f98e4de9f0874f88928003..00fc57d20ee37ed3bd4285e3c8b5bf56b856077f 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/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);
diff --git a/AMDiS/src/parallel/PetscSolverFetiDebug.cc b/AMDiS/src/parallel/PetscSolverFetiDebug.cc
index ad575e210fd4f95e2dc5a367c734b09876ac509f..09db80e7389f75ece02113ac751199c59dae1bd8 100644
--- a/AMDiS/src/parallel/PetscSolverFetiDebug.cc
+++ b/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. ===
+
+ {
+
+ }
+
+ }
+
}
diff --git a/AMDiS/src/parallel/PetscSolverFetiDebug.h b/AMDiS/src/parallel/PetscSolverFetiDebug.h
index e151158e98980df37f8be2760f9ca8215b008337..145e2d210c530b95c170b3f6624d0f0f374e1fc2 100644
--- a/AMDiS/src/parallel/PetscSolverFetiDebug.h
+++ b/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);
};
}