diff --git a/AMDiS/src/AdaptInstationary.cc b/AMDiS/src/AdaptInstationary.cc
index 3ed06ba55e2b3e2ef967c14d72c7cce3b6951f09..e4571ccf51b6dad94ddda7bb5c15664b193e622c 100644
--- a/AMDiS/src/AdaptInstationary.cc
+++ b/AMDiS/src/AdaptInstationary.cc
@@ -108,7 +108,7 @@ namespace AMDiS {
INFO(info_,6)("time = %e, try timestep = %e\n",
adaptInfo->getTime(), adaptInfo->getTimestep());
-
+
problemIteration_->oneIteration(adaptInfo, NO_ADAPTION);
adaptInfo->incTimestepIteration();
diff --git a/AMDiS/src/ElementFileWriter.cc b/AMDiS/src/ElementFileWriter.cc
index 35884163991c4f480ce8908e44fa6df2f266cddd..8710c719638639794359a322ffc496b57010ad2f 100644
--- a/AMDiS/src/ElementFileWriter.cc
+++ b/AMDiS/src/ElementFileWriter.cc
@@ -6,7 +6,7 @@
namespace AMDiS {
- ElementFileWriter::ElementFileWriter(const std::string& name_,
+ ElementFileWriter::ElementFileWriter(std::string name_,
const FiniteElemSpace *feSpace_,
std::map<int, double> &mapvec)
: name(name_),
@@ -95,14 +95,14 @@ namespace AMDiS {
void ElementFileWriter::writeFile(std::map<int, double> &vec,
const FiniteElemSpace *feSpace,
- const std::string& filename)
+ std::string filename)
{
ElementFileWriter efw("", feSpace, vec);
efw.writeVtkValues(filename);
}
- void ElementFileWriter::writeTecPlotValues(const std::string &filename)
+ void ElementFileWriter::writeTecPlotValues(std::string filename)
{
FUNCNAME("ElementFileWriter::writeTecPlotValues()");
std::ofstream fout(filename.c_str());
@@ -154,10 +154,9 @@ namespace AMDiS {
// Write coordinates of all element vertices and element value.
for (int i = 0; i <= dim; ++i) {
-
- for (int j = 0; j < dim; ++j) {
+ for (int j = 0; j < dim; ++j)
fout << elInfo->getCoord(i)[j] << " ";
- }
+
fout << val << "\n";
}
@@ -182,7 +181,7 @@ namespace AMDiS {
fout.close();
}
- void ElementFileWriter::writeMeshDatValues(const std::string &filename, double time)
+ void ElementFileWriter::writeMeshDatValues(std::string filename, double time)
{
FUNCNAME("ElementFileWriter::writeMeshDatValues()");
std::ofstream fout(filename.c_str());
@@ -281,7 +280,7 @@ namespace AMDiS {
fout.close();
}
- void ElementFileWriter::writeVtkValues(const std::string &filename)
+ void ElementFileWriter::writeVtkValues(std::string filename)
{
FUNCNAME("ElementFileWriter::writeVtkValues()");
std::ofstream fout(filename.c_str());
@@ -289,6 +288,7 @@ namespace AMDiS {
TEST_EXIT(fout)("Could not open file %s !\n", filename.c_str());
int dim = mesh->getDim();
+ int dimOfWorld = Global::getGeo(WORLD);
int vertices = mesh->getGeo(VERTEX);
int nElements = mesh->getNumberOfLeaves();
double val;
@@ -312,12 +312,12 @@ namespace AMDiS {
while (elInfo) {
// Write coordinates of all element vertices.
for (int i = 0; i <= dim; i++) {
- for (int j = 0; j < dim; j++) {
+ for (int j = 0; j < dimOfWorld; j++)
fout << elInfo->getCoord(i)[j] << " ";
- }
- for (int j = dim; j < 3; j++) {
+
+ for (int j = dimOfWorld; j < 3; j++)
fout << "0.0";
- }
+
fout << "\n";
}
@@ -329,9 +329,8 @@ namespace AMDiS {
fout << " <Cells>" << std::endl;
fout << " <DataArray type=\"Int32\" Name=\"offsets\">" << std::endl;
- for (int i = 0; i < nElements; i++) {
- fout << " " << (i + 1) * vertices << std::endl;
- }
+ for (int i = 0; i < nElements; i++)
+ fout << " " << (i + 1) * vertices << std::endl;
fout << " </DataArray>" << std::endl;
diff --git a/AMDiS/src/ElementFileWriter.h b/AMDiS/src/ElementFileWriter.h
index 67cda223329101255c9012c100cdabe74c45182e..74e16582a6a507d8a95158483969bba22ac5b82b 100644
--- a/AMDiS/src/ElementFileWriter.h
+++ b/AMDiS/src/ElementFileWriter.h
@@ -21,7 +21,7 @@ namespace AMDiS {
{
public:
/// Constructor.
- ElementFileWriter(const std::string& name,
+ ElementFileWriter(std::string name,
const FiniteElemSpace *feSpace,
std::map<int, double> &vec);
@@ -34,17 +34,17 @@ namespace AMDiS {
/// Simple writing procedure for one element map.
static void writeFile(std::map<int, double> &vec,
const FiniteElemSpace *feSpace,
- const std::string& filename);
+ std::string filename);
protected:
/// Writes element data in tecplot format.
- void writeTecPlotValues(const std::string &filename);
+ void writeTecPlotValues(std::string filename);
/// Writes element data in AMDiS format (1 file !).
- void writeMeshDatValues(const std::string &filename, double time);
+ void writeMeshDatValues(std::string filename, double time);
/// Writes element data in VTK format.
- void writeVtkValues(const std::string &filename);
+ void writeVtkValues(std::string filename);
protected:
/// Name.
diff --git a/AMDiS/src/InteriorBoundary.h b/AMDiS/src/InteriorBoundary.h
index 077d6d40ea20933c9506582cd2fcdb734e8ec2bb..6de459da1f08f88d1b8ba6011014aca9e830af53 100644
--- a/AMDiS/src/InteriorBoundary.h
+++ b/AMDiS/src/InteriorBoundary.h
@@ -35,6 +35,9 @@ namespace AMDiS {
/// The macro element to which the boundary element corresponds to.
Element* el;
+ /// Index of the macro element.
+ int elIndex;
+
/** \brief
* Defines the geometrical object at the boundary. It must be "a part" of the
* macro element \ref el, i.e., either 1 (a vertex), 2 (an edge) or 3 (a face).
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 017476fbb127dbe1dd26d273606fc1464f70cef7..e13fbca3dba84857a7f9902ada58b8b60cf2e3e3 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -48,12 +48,11 @@ namespace AMDiS {
Lagrange::~Lagrange()
{
- for (int i = 0; i < static_cast<int>(bary->size()); i++) {
+ for (int i = 0; i < static_cast<int>(bary->size()); i++)
if ((*bary)[i]) {
delete (*bary)[i];
(*bary)[i] = NULL;
}
- }
}
Lagrange* Lagrange::getLagrange(int dim, int degree)
@@ -70,13 +69,12 @@ namespace AMDiS {
void Lagrange::clear()
{
- std::list<Lagrange*>::iterator it;
- for (it = allBasFcts.begin(); it != allBasFcts.end(); it++) {
+ for (std::list<Lagrange*>::iterator it = allBasFcts.begin();
+ it != allBasFcts.end(); it++)
if (*it) {
delete *it;
*it = NULL;
}
- }
}
void Lagrange::setFunctionPointer()
diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc
index 50ef534fb489b3431a9f0e34ee8db31686c5d287..21bbe419ccb30839bb7f6662df0490088ebb432f 100644
--- a/AMDiS/src/ParallelDomainBase.cc
+++ b/AMDiS/src/ParallelDomainBase.cc
@@ -15,6 +15,7 @@
#include "ElementDofIterator.h"
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
+#include "ElementFileWriter.h"
#include "petscksp.h"
@@ -22,7 +23,7 @@ namespace AMDiS {
PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
{
- if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
+ // if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
std::cout << " Iteration " << iter << ": " << rnorm << std::endl;
return 0;
@@ -66,6 +67,23 @@ namespace AMDiS {
if (mpiSize <= 1)
return;
+#if 0
+ if (mpiRank == 0) {
+ std::map<int, double> vec;
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
+ while (elInfo) {
+ vec[elInfo->getElement()->getIndex()] = static_cast<double>(elInfo->getElement()->getIndex());
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ ElementFileWriter::writeFile(vec, feSpace, "test.vtu");
+ }
+#endif
+
+
// Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
// only for macro meshes, so it will not work yet if the mesh is already refined
// in some way.
@@ -112,6 +130,41 @@ namespace AMDiS {
updateDofAdmins();
+#if 0
+ if (mpiRank == 0) {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ if (elInfo->getElement()->getIndex() == 4) {
+ WorldVector<double> x;
+ mesh->getDofIndexCoords(elInfo->getElement()->getDOF(0), feSpace, x);
+ std::cout << "FOUND!" << std::endl;
+ x.print();
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+
+ if (mpiRank == 1) {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ if (elInfo->getElement()->getIndex() == 3) {
+ WorldVector<double> x;
+ mesh->getDofIndexCoords(elInfo->getElement()->getDOF(2), feSpace, x);
+ std::cout << "FOUND!" << std::endl;
+ x.print();
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+ exit(0);
+#endif
+
// === Global refinements. ===
@@ -137,7 +190,7 @@ namespace AMDiS {
#if (DEBUG != 0)
- DbgTestCommonDofs();
+ DbgTestCommonDofs(true);
#endif
@@ -293,7 +346,7 @@ namespace AMDiS {
KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
KSPGetPC(ksp, &pc);
// PCSetType(pc, PCNONE);
- PCSetType(pc, PCJACOBI);
+ PCSetType(pc, PCILU);
KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetType(ksp, KSPBCGS);
//KSPSetType(ksp, KSPCG);
@@ -555,6 +608,24 @@ namespace AMDiS {
bool isRankDOF2 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF2) != rankDOFs.end());
bool ranksBoundary = isRankDOF1 || isRankDOF2;
+#if 0
+ if (mpiRank == 3 && ranksBoundary &&
+ partitionVec[elInfo->getNeighbour(i)->getIndex()] == 2) {
+ std::cout << "ADD MY BOUND " << element->getIndex() << "/" << i
+ << " with "
+ << elInfo->getNeighbour(i)->getIndex() << "/"
+ << elInfo->getSideOfNeighbour(i) << std::endl;
+ }
+
+ if (mpiRank == 2 && !ranksBoundary &&
+ partitionVec[elInfo->getNeighbour(i)->getIndex()] == 3) {
+ std::cout << "ADD OT BOUND " << element->getIndex() << "/" << i
+ << " with "
+ << elInfo->getNeighbour(i)->getIndex() << "/"
+ << elInfo->getSideOfNeighbour(i) << std::endl;
+ }
+#endif
+
// === And add the part of the interior boundary. ===
AtomicBoundary& bound =
@@ -563,11 +634,19 @@ namespace AMDiS {
otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));
bound.rankObject.el = element;
+ bound.rankObject.elIndex = element->getIndex();
bound.rankObject.subObjAtBoundary = EDGE;
bound.rankObject.ithObjAtBoundary = i;
- bound.neighbourObject.el = elInfo->getNeighbour(i);
+ // Do not set a pointer to the element, because if will be deleted from
+ // mesh after partitioning and the pointer would become invalid.
+ bound.neighbourObject.el = NULL;
+ bound.neighbourObject.elIndex = elInfo->getNeighbour(i)->getIndex();
bound.neighbourObject.subObjAtBoundary = EDGE;
- bound.neighbourObject.ithObjAtBoundary = -1;
+ bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);
+
+ // i == 2 => getSideOfNeighbour(i) == 2
+ TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
+ ("Should not happen!\n");
}
}
}
@@ -589,23 +668,25 @@ namespace AMDiS {
for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
rankIt != myIntBoundary.boundary.end(); ++rankIt) {
int nSendInt = rankIt->second.size();
- int* buffer = new int[nSendInt];
- for (int i = 0; i < nSendInt; i++)
- buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
+ int* buffer = new int[nSendInt * 2];
+ for (int i = 0; i < nSendInt; i++) {
+ buffer[i * 2] = (rankIt->second)[i].rankObject.elIndex;
+ buffer[i * 2 + 1] = (rankIt->second)[i].rankObject.ithObjAtBoundary;
+ }
sendBuffers.push_back(buffer);
request[requestCounter++] =
- mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
+ mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
}
for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
int nRecvInt = rankIt->second.size();
- int *buffer = new int[nRecvInt];
+ int *buffer = new int[nRecvInt * 2];
recvBuffers.push_back(buffer);
request[requestCounter++] =
- mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
+ mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
}
@@ -622,10 +703,12 @@ namespace AMDiS {
for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
// If the expected object is not at place, search for it.
- if ((rankIt->second)[j].neighbourObject.el->getIndex() != recvBuffers[i][j]) {
+ if ((rankIt->second)[j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
+ (rankIt->second)[j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
int k = j + 1;
for (; k < static_cast<int>(rankIt->second.size()); k++)
- if ((rankIt->second)[k].neighbourObject.el->getIndex() == recvBuffers[i][j])
+ if ((rankIt->second)[k].neighbourObject.elIndex == recvBuffers[i][j * 2] &&
+ (rankIt->second)[k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
break;
// The element must always be found, because the list is just in another order.
@@ -790,6 +873,11 @@ namespace AMDiS {
sendBuffers[i][c++] = *(dofIt->first);
sendBuffers[i][c++] = dofIt->second;
+#if 0
+ if (mpiRank == 3 && sendIt->first == 2)
+ std::cout << "SEND DOF: " << dofIt->first << std::endl;
+#endif
+
sendDofs[sendIt->first].push_back(dofIt->first);
}
@@ -852,6 +940,11 @@ namespace AMDiS {
if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) {
dofChanged[dofIt->first] = true;
+#if 0
+ if (mpiRank == 2 && recvIt->first == 3)
+ std::cout << "RECV DOF: " << dofIt->first << std::endl;
+#endif
+
recvDofs[recvIt->first].push_back(dofIt->first);
rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
@@ -897,8 +990,15 @@ namespace AMDiS {
}
DofContainer rankAllDofs;
- for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt)
+ for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt) {
+ /* if (mpiRank == 1) {
+ std::cout << "COORDs of dof = " << **dofIt << std::endl;
+ WorldVector<double> x;
+ mesh->getDofIndexCoords(*dofIt, feSpace, x);
+ x.print();
+ }*/
rankAllDofs.push_back(*dofIt);
+ }
sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
DofContainer rankDOFs = rankAllDofs;
@@ -915,6 +1015,16 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
+#if 0
+ if (mpiRank == 3 && it->first == 2)
+ std::cout << "GO ON MY BOUND: " << boundIt->rankObject.elIndex
+ << "/" << boundIt->rankObject.ithObjAtBoundary << " with "
+ << boundIt->neighbourObject.elIndex << "/"
+ << boundIt->neighbourObject.ithObjAtBoundary
+ << std::endl;
+#endif
+
+
DofContainer dofs;
DofContainer &dofsToSend = sendDofs[it->first];
@@ -935,22 +1045,32 @@ namespace AMDiS {
ERROR_EXIT("Should never happen!\n");
}
- for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) {
- if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end())
- dofsToSend.push_back(*dofIt);
+#if 0
+ if (mpiRank == 3 && it->first == 2) {
+ WorldVector<double> x;
+ mesh->getDofIndexCoords(dofs[0], feSpace, x);
+ x.print();
+ mesh->getDofIndexCoords(dofs[1], feSpace, x);
+ x.print();
}
+#endif
+
+
+ for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt)
+ if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end())
+ dofsToSend.push_back(*dofIt);
dofs.clear();
addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs);
addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs);
- for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
+ for (int i = 0; i < static_cast<int>(dofs.size()); i++)
dofsToSend.push_back(dofs[i]);
- }
-
}
- }
+ }
+
+
for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
it != otherIntBoundary.boundary.end(); ++it) {
@@ -958,17 +1078,37 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
+#if 0
+ if (mpiRank == 2 && it->first == 3)
+ std::cout << "GO ON OT BOUND: " << boundIt->rankObject.elIndex
+ << "/" << boundIt->rankObject.ithObjAtBoundary << " with "
+ << boundIt->neighbourObject.elIndex << "/"
+ << boundIt->neighbourObject.ithObjAtBoundary
+ << std::endl;
+#endif
+
+
DofContainer dofs;
DofContainer &dofsToRecv = recvDofs[it->first];
switch (boundIt->rankObject.ithObjAtBoundary) {
case 0:
- dofs.push_back(boundIt->rankObject.el->getDOF(1));
- dofs.push_back(boundIt->rankObject.el->getDOF(2));
+ if (boundIt->neighbourObject.ithObjAtBoundary == 0) {
+ dofs.push_back(boundIt->rankObject.el->getDOF(2));
+ dofs.push_back(boundIt->rankObject.el->getDOF(1));
+ } else {
+ dofs.push_back(boundIt->rankObject.el->getDOF(1));
+ dofs.push_back(boundIt->rankObject.el->getDOF(2));
+ }
break;
case 1:
- dofs.push_back(boundIt->rankObject.el->getDOF(0));
- dofs.push_back(boundIt->rankObject.el->getDOF(2));
+ if (boundIt->neighbourObject.ithObjAtBoundary == 0) {
+ dofs.push_back(boundIt->rankObject.el->getDOF(0));
+ dofs.push_back(boundIt->rankObject.el->getDOF(2));
+ } else {
+ dofs.push_back(boundIt->rankObject.el->getDOF(2));
+ dofs.push_back(boundIt->rankObject.el->getDOF(0));
+ }
break;
case 2:
dofs.push_back(boundIt->rankObject.el->getDOF(1));
@@ -978,6 +1118,16 @@ namespace AMDiS {
ERROR_EXIT("Should never happen!\n");
}
+#if 0
+ if (mpiRank == 2 && it->first == 3) {
+ WorldVector<double> x;
+ mesh->getDofIndexCoords(dofs[0], feSpace, x);
+ x.print();
+ mesh->getDofIndexCoords(dofs[1], feSpace, x);
+ x.print();
+ }
+#endif
+
for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) {
DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), *dofIt);
if (eraseIt != rankDOFs.end())
@@ -1002,6 +1152,7 @@ namespace AMDiS {
dofsToRecv.push_back(dofs[i]);
}
+
}
}
@@ -1344,10 +1495,25 @@ namespace AMDiS {
for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
rankIt != myIntBoundary.boundary.end(); ++rankIt) {
+
int nSendInt = rankIt->second.size();
int* buffer = new int[nSendInt];
- for (int i = 0; i < nSendInt; i++)
- buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
+ for (int i = 0; i < nSendInt; i++) {
+
+#if 0
+ if (mpiRank == 3 && rankIt->first == 2) {
+ std::cout << "MY BOUND: " << (rankIt->second)[i].rankObject.elIndex
+ << "/" << (rankIt->second)[i].rankObject.ithObjAtBoundary
+ << " with neig "
+ << (rankIt->second)[i].neighbourObject.elIndex
+ << "/" << (rankIt->second)[i].neighbourObject.ithObjAtBoundary
+ << std::endl;
+ }
+#endif
+
+
+ buffer[i] = (rankIt->second)[i].rankObject.elIndex;
+ }
sendBuffers.push_back(buffer);
request[requestCounter++] =
@@ -1378,8 +1544,22 @@ namespace AMDiS {
("Boundaries does not fit together!\n");
for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++) {
+
+#if 0
+ if (mpiRank == 2 && rankIt->first == 3) {
+ std::cout << "OT BOUND: " << (rankIt->second)[i].rankObject.elIndex
+ << "/" << (rankIt->second)[i].rankObject.ithObjAtBoundary
+ << " with neig "
+ << (rankIt->second)[i].neighbourObject.elIndex
+ << "/" << (rankIt->second)[i].neighbourObject.ithObjAtBoundary
+ << std::endl;
+ }
+
+#endif
+
+
int elIndex1 = recvBuffers[bufCounter][i];
- int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighbourObject.el->getIndex();
+ int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighbourObject.elIndex;
TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n");
}
@@ -1393,6 +1573,27 @@ namespace AMDiS {
{
FUNCNAME("ParallelDomainBase::DbgTestCommonDofs()");
+
+#if 0
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ Element *el = elInfo->getElement();
+ for (int i = 0; i < 3; i++) {
+ if (mpiRank == 2 && el->getDOF(i, 0) == 3)
+ std::cout << "RANK 2 EL = " << el->getIndex() << " i = " << i << std::endl;
+
+ if (mpiRank == 3 && el->getDOF(i, 0) == 4)
+ std::cout << "RANK 3 EL = " << el->getIndex() << " i = " << i << std::endl;
+
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+#endif
+
// Maps to each neighbour rank an array of WorldVectors. This array contain the
// coordinates of all dofs this rank shares on the interior boundary with the
// neighbour rank. A rank sends the coordinates to another rank, if it owns the
@@ -1411,6 +1612,14 @@ namespace AMDiS {
dofIt != it->second.end(); ++dofIt) {
bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
TEST_EXIT(b)("Cannot find DOF in mesh!\n");
+
+#if 0
+ if (mpiRank == 3 && it->first == 2) {
+ std::cout << "SEND COORDS: " << *dofIt << " = " << **dofIt << std::endl;
+ coords.print();
+ }
+#endif
+
sendCoords[it->first].push_back(coords);
}
}
@@ -1422,6 +1631,14 @@ namespace AMDiS {
bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
TEST_EXIT(b)("Cannot find DOF in mesh!\n");
recvCoords[it->first].push_back(coords);
+
+#if 0
+ if (mpiRank == 2 && it->first == 3) {
+ std::cout << "RECV COORDS: " << *dofIt << " = " << **dofIt << std::endl;
+ coords.print();
+ }
+#endif
+
}
}
@@ -1470,24 +1687,25 @@ namespace AMDiS {
std::map<int, double*> sendCoordsBuffer, recvCoordsBuffer;
requestCounter = 0;
- for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) {
- sendCoordsBuffer[it->first] = new double[it->second.size() * 2];
+ int dimOfWorld = Global::getGeo(WORLD);
- for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
- sendCoordsBuffer[it->first][i * 2] = (it->second)[i][0];
- sendCoordsBuffer[it->first][i * 2 + 1] = (it->second)[i][1];
- }
+ for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) {
+ sendCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld];
+ for (int i = 0; i < static_cast<int>(it->second.size()); i++)
+ for (int j = 0; j < dimOfWorld; j++)
+ sendCoordsBuffer[it->first][i * dimOfWorld + j] = (it->second)[i][j];
+
request[requestCounter++] = mpiComm.Isend(sendCoordsBuffer[it->first],
- it->second.size() * 2,
+ it->second.size() * dimOfWorld,
MPI_DOUBLE, it->first, 0);
}
for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) {
- recvCoordsBuffer[it->first] = new double[it->second.size() * 2];
+ recvCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld];
request[requestCounter++] = mpiComm.Irecv(recvCoordsBuffer[it->first],
- it->second.size() * 2,
+ it->second.size() * dimOfWorld,
MPI_DOUBLE, it->first, 0);
}
@@ -1500,24 +1718,27 @@ namespace AMDiS {
for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) {
for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
-
- if (printCoords) {
- std::cout << "[DBG] i = " << i << std::endl;
+ for (int j = 0; j < dimOfWorld; j++) {
+ bool c = fabs((it->second)[i][j] -
+ recvCoordsBuffer[it->first][i * dimOfWorld + j]) < eps;
- std::cout.precision(20);
- std::cout << "[DBG] "
- << "Rank " << mpiRank << " from rank " << it->first
- << " expect coords ("
- << (it->second)[i][0] << " , " << (it->second)[i][1]
- << ") received coords ("
- << recvCoordsBuffer[it->first][i * 2] << " , "
- << recvCoordsBuffer[it->first][i * 2 + 1] << ")" << std::endl;
- }
-
- bool c0 = fabs((it->second)[i][0] - recvCoordsBuffer[it->first][i * 2]) < eps;
- bool c1 = fabs((it->second)[i][1] - recvCoordsBuffer[it->first][i * 2 + 1]) < eps;
+ if (printCoords && !c) {
+ std::cout << "[DBG] i = " << i << std::endl;
+
+ std::cout.precision(5);
+ std::cout << "[DBG] "
+ << "Rank " << mpiRank << " from rank " << it->first
+ << " expect coords (";
+ for (int k = 0; k < dimOfWorld; k++)
+ std::cout << (it->second)[i][k] << " , ";
+ std::cout << ") received coords (";
+ for (int k = 0; k < dimOfWorld; k++)
+ std::cout << recvCoordsBuffer[it->first][i * dimOfWorld + k] << " , ";
+ std::cout << ")" << std::endl;
+ }
- TEST_EXIT(c0 && c1)("Wrong DOFs in rank %d!\n", mpiRank);
+ TEST_EXIT(c)("Wrong DOFs in rank %d!\n", mpiRank);
+ }
}
delete [] recvCoordsBuffer[it->first];