Commit 4873445c authored by Thomas Witkowski's avatar Thomas Witkowski

Fixed problem when using multiple periodic boundary conditions in parallel domain decomposition.

parent 6d3687b7
......@@ -371,6 +371,35 @@ namespace AMDiS {
/// Returns whether Element has sideElem as one of its sides.
virtual bool hasSide(Element *sideElem) const = 0;
/** \brief
* Traverses an edge of a given element (this includes also all children of the
* element having the same edge). All vertex dofs alonge this edge are assembled
* and put together to a list.
*
* \param[in] feSpace FE space which is used to get the dofs.
* \param[in] ithEdge Defines the edge on which all the vertex dofs
* are assembled.
* \param[out] dofs List of dofs, where the result is stored.
* \param[in] parentVertices If true, also the two vertices of the parent
* element are put into the result list.
*/
virtual void getVertexDofs(FiniteElemSpace* feSpace, int ithEdge,
DofContainer& dofs, bool parentVertices = 0) const = 0;
/** \brief
* Traverses an edge of a given element (this includes also all children of the
* element having the same edge). All non vertex dofs alonge this edge are
* assembled and put together to a list.
*
* \param[in] feSpace FE space which is used to get the dofs.
* \param[in] ithEdge Defines the edge on which all the non vertex
* dofs are assembled.
* \param[out] dofs All dofs are put to this dof list.
*/
virtual void getNonVertexDofs(FiniteElemSpace* feSpace, int ithEdge,
DofContainer& dofs) const = 0;
/** \} */
// ===== other public methods =================================================
......
......@@ -6,7 +6,7 @@
namespace AMDiS {
void ElementDofIterator::reset(Element* el)
void ElementDofIterator::reset(const Element* el)
{
FUNCNAME("ElementDofIterator::reset()");
......
......@@ -54,7 +54,7 @@ namespace AMDiS {
{}
/// Start a new traverse with the given element.
void reset(Element* el);
void reset(const Element* el);
/// Go to next dof. Returns false, if there is dof anymore.
bool next();
......@@ -116,7 +116,7 @@ namespace AMDiS {
int* orderPosition;
Element* element;
const Element* element;
/// Current position (i.e., vertex, edge, face) of the traverse.
int pos;
......
......@@ -74,6 +74,9 @@ namespace AMDiS {
/// datatype for degrees of freedom
typedef signed int DegreeOfFreedom;
/// Defines type for a vector of DOF pointers.
typedef std::vector<const DegreeOfFreedom*> DofContainer;
/// returns the GeoIndex of d for dimension dim.
#define INDEX_OF_DIM(d, dim) (static_cast<GeoIndex>((d == dim) ? CENTER : d + 1))
......
......@@ -21,13 +21,13 @@ namespace AMDiS {
for (int i = 0; i < boundSize; i++) {
AtomicBoundary &bound = (it->second)[i];
SerUtil::serialize(out, bound.rankObject.elIndex);
SerUtil::serialize(out, bound.rankObject.subObjAtBoundary);
SerUtil::serialize(out, bound.rankObject.ithObjAtBoundary);
SerUtil::serialize(out, bound.rankObj.elIndex);
SerUtil::serialize(out, bound.rankObj.subObj);
SerUtil::serialize(out, bound.rankObj.ithObj);
SerUtil::serialize(out, bound.neighbourObject.elIndex);
SerUtil::serialize(out, bound.neighbourObject.subObjAtBoundary);
SerUtil::serialize(out, bound.neighbourObject.ithObjAtBoundary);
SerUtil::serialize(out, bound.neighObj.elIndex);
SerUtil::serialize(out, bound.neighObj.subObj);
SerUtil::serialize(out, bound.neighObj.ithObj);
}
}
}
......@@ -47,16 +47,16 @@ namespace AMDiS {
for (int i = 0; i < boundSize; i++) {
AtomicBoundary &bound = boundary[rank][i];
SerUtil::deserialize(in, bound.rankObject.elIndex);
SerUtil::deserialize(in, bound.rankObject.subObjAtBoundary);
SerUtil::deserialize(in, bound.rankObject.ithObjAtBoundary);
SerUtil::deserialize(in, bound.rankObj.elIndex);
SerUtil::deserialize(in, bound.rankObj.subObj);
SerUtil::deserialize(in, bound.rankObj.ithObj);
SerUtil::deserialize(in, bound.neighbourObject.elIndex);
SerUtil::deserialize(in, bound.neighbourObject.subObjAtBoundary);
SerUtil::deserialize(in, bound.neighbourObject.ithObjAtBoundary);
SerUtil::deserialize(in, bound.neighObj.elIndex);
SerUtil::deserialize(in, bound.neighObj.subObj);
SerUtil::deserialize(in, bound.neighObj.ithObj);
bound.rankObject.el = elIndexMap[bound.rankObject.elIndex];
bound.neighbourObject.el = NULL;
bound.rankObj.el = elIndexMap[bound.rankObj.elIndex];
bound.neighObj.el = NULL;
}
}
}
......
......@@ -43,18 +43,18 @@ namespace AMDiS {
* 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).
*/
GeoIndex subObjAtBoundary;
GeoIndex subObj;
/** \brief
* Defines which of vertix, edge or face of the macro element is part of the
* boundary.
*
* Example: If the macro element is a triangle, than \ref subObjAtBoundary may
* be either 1 (vertex) or 2 (edge). Assume its the last one. So this variable
* defines which of the three possible edges of the triangle is at the interior
* Example: If the macro element is a triangle, than \ref subObj may be either
* 1 (vertex) or 2 (edge). Assume its the last one. So this variable defines
* which of the three possible edges of the triangle is at the interior
* boundary.
*/
int ithObjAtBoundary;
int ithObj;
};
/** \brief
......@@ -63,10 +63,10 @@ namespace AMDiS {
*/
struct AtomicBoundary {
/// The rank's part of the boundary.
BoundaryObject rankObject;
BoundaryObject rankObj;
/// The object on the other side of the boundary.
BoundaryObject neighbourObject;
BoundaryObject neighObj;
};
/** \brief
......
......@@ -124,19 +124,19 @@ namespace AMDiS {
return false;
}
/// implements Element::isLine. Returns true because this element is a Line
/// Returns true because this element is a Line.
inline bool isLine() const
{
return true;
}
/// implements Element::isTriangle. Returns false because this element is a Line
/// Returns false because this element is a Line.
inline bool isTriangle() const
{
return false;
}
/// implements Element::isTetrahedron. Returns false because this element is a Line
/// Returns false because this element is a Line
inline bool isTetrahedron() const
{
return false;
......@@ -147,6 +147,22 @@ namespace AMDiS {
return "Line";
}
void getVertexDofs(FiniteElemSpace* feSpace, int ithEdge,
DofContainer& dofs, bool parentVertices = 0) const
{
FUNCNAME("Line::getVertexDofs()");
ERROR_EXIT("Not yet implemented!\n");
}
void getNonVertexDofs(FiniteElemSpace* feSpace, int ithEdge,
DofContainer& dofs) const
{
FUNCNAME("Line::getNonVertexDofs()");
ERROR_EXIT("Not yet implemented!\n");
}
protected:
/** \brief
* vertexOfEdge[i][j] is the local number of the j-th vertex of the i-th
......
......@@ -153,6 +153,7 @@ namespace AMDiS {
createPeriodicMap();
}
// exit(0);
#if (DEBUG != 0)
DbgTestCommonDofs(true);
......@@ -982,16 +983,16 @@ namespace AMDiS {
// === Create information about the boundary between the two elements. ===
AtomicBoundary bound;
bound.rankObject.el = element;
bound.rankObject.elIndex = element->getIndex();
bound.rankObject.subObjAtBoundary = EDGE;
bound.rankObject.ithObjAtBoundary = i;
bound.rankObj.el = element;
bound.rankObj.elIndex = element->getIndex();
bound.rankObj.subObj = EDGE;
bound.rankObj.ithObj = 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 = elInfo->getSideOfNeighbour(i);
bound.neighObj.el = NULL;
bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
bound.neighObj.subObj = EDGE;
bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
// i == 2 => getSideOfNeighbour(i) == 2
TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
......@@ -1049,8 +1050,8 @@ namespace AMDiS {
int nSendInt = rankIt->second.size();
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;
buffer[i * 2] = (rankIt->second)[i].rankObj.elIndex;
buffer[i * 2 + 1] = (rankIt->second)[i].rankObj.ithObj;
}
sendBuffers.push_back(buffer);
......@@ -1086,12 +1087,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.elIndex != recvBuffers[i][j * 2] ||
(rankIt->second)[j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
if ((rankIt->second)[j].neighObj.elIndex != recvBuffers[i][j * 2] ||
(rankIt->second)[j].neighObj.ithObj != recvBuffers[i][j * 2 + 1]) {
int k = j + 1;
for (; k < static_cast<int>(rankIt->second.size()); k++)
if ((rankIt->second)[k].neighbourObject.elIndex == recvBuffers[i][j * 2] &&
(rankIt->second)[k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
if ((rankIt->second)[k].neighObj.elIndex == recvBuffers[i][j * 2] &&
(rankIt->second)[k].neighObj.ithObj == recvBuffers[i][j * 2 + 1])
break;
// The element must always be found, because the list is just in another order.
......@@ -1130,8 +1131,8 @@ namespace AMDiS {
int nSendInt = rankIt->second.size();
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;
buffer[i * 2] = (rankIt->second)[i].rankObj.elIndex;
buffer[i * 2 + 1] = (rankIt->second)[i].rankObj.ithObj;
}
sendBuffers.push_back(buffer);
......@@ -1157,13 +1158,13 @@ namespace AMDiS {
if (rankIt->first > mpiRank) {
for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++)
if (periodicBoundary.boundary[rankIt->first][j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
periodicBoundary.boundary[rankIt->first][j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
if (periodicBoundary.boundary[rankIt->first][j].neighObj.elIndex != recvBuffers[i][j * 2] ||
periodicBoundary.boundary[rankIt->first][j].neighObj.ithObj != recvBuffers[i][j * 2 + 1]) {
int k = j + 1;
for (; k < static_cast<int>(rankIt->second.size()); k++)
if (periodicBoundary.boundary[rankIt->first][k].neighbourObject.elIndex == recvBuffers[i][j * 2] &&
periodicBoundary.boundary[rankIt->first][k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
if (periodicBoundary.boundary[rankIt->first][k].neighObj.elIndex == recvBuffers[i][j * 2] &&
periodicBoundary.boundary[rankIt->first][k].neighObj.ithObj == recvBuffers[i][j * 2 + 1])
break;
// The element must always be found, because the list is just in another order.
......@@ -1505,10 +1506,9 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
DofContainer dofs;
addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs);
addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs);
boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
TEST_EXIT_DBG(find(dofsToSend.begin(), dofsToSend.end(), dofs[i]) == dofsToSend.end())
("Should not happen!\n");
......@@ -1528,10 +1528,8 @@ namespace AMDiS {
boundIt != it->second.end(); ++boundIt) {
DofContainer dofs;
addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs);
addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs);
boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
for (int i = static_cast<int>(dofs.size()) - 1; i >= 0; i--) {
TEST_EXIT_DBG(find(dofsToRecv.begin(), dofsToRecv.end(), dofs[i]) == dofsToRecv.end())
......@@ -1666,121 +1664,7 @@ namespace AMDiS {
mapLocalToDofIndex[dofIt->second] = *(dofIt->first);
}
void ParallelDomainBase::addVertexDofs(Element *el, int ithEdge, DofContainer& dofs,
bool parentVertices)
{
FUNCNAME("ParallelDomainBase::addVertexDofs()");
TEST_EXIT_DBG(mesh->getDim() == 2)("Not yet implemented for 3D!\n");
if (parentVertices) {
switch (ithEdge) {
case 0:
dofs.push_back(el->getDOF(1));
break;
case 1:
dofs.push_back(el->getDOF(0));
break;
case 2:
dofs.push_back(el->getDOF(0));
break;
default:
ERROR_EXIT("Should never happen!\n");
}
}
switch (ithEdge) {
case 0:
if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
addVertexDofs(el->getSecondChild()->getFirstChild(), 0, dofs);
dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
addVertexDofs(el->getSecondChild()->getSecondChild(), 1, dofs);
}
break;
case 1:
if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
addVertexDofs(el->getFirstChild()->getFirstChild(), 0, dofs);
dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
addVertexDofs(el->getFirstChild()->getSecondChild(), 1, dofs);
}
break;
case 2:
if (el->getFirstChild()) {
addVertexDofs(el->getFirstChild(), 0, dofs);
dofs.push_back(el->getFirstChild()->getDOF(2));
addVertexDofs(el->getSecondChild(), 1, dofs);
}
break;
default:
ERROR_EXIT("Should never happen!\n");
}
if (parentVertices) {
switch (ithEdge) {
case 0:
dofs.push_back(el->getDOF(2));
break;
case 1:
dofs.push_back(el->getDOF(2));
break;
case 2:
dofs.push_back(el->getDOF(1));
break;
default:
ERROR_EXIT("Should never happen!\n");
}
}
}
void ParallelDomainBase::addEdgeDofs(Element *el, int ithEdge, DofContainer& dofs)
{
FUNCNAME("ParallelDomainBase::addEdgeDofs()");
TEST_EXIT_DBG(mesh->getDim() == 2)("Not yet implemented for 3D!\n");
bool addThisEdge = false;
switch (ithEdge) {
case 0:
if (el->getSecondChild())
addEdgeDofs(el->getSecondChild(), 2, dofs);
else
addThisEdge = true;
break;
case 1:
if (el->getFirstChild())
addEdgeDofs(el->getFirstChild(), 2, dofs);
else
addThisEdge = true;
break;
case 2:
if (el->getFirstChild()) {
addEdgeDofs(el->getFirstChild(), 0, dofs);
addEdgeDofs(el->getSecondChild(), 1, dofs);
} else {
addThisEdge = true;
}
break;
default:
ERROR_EXIT("Should never happen!\n");
}
if (addThisEdge) {
ElementDofIterator elDofIter(feSpace, true);
elDofIter.reset(el);
do {
if (elDofIter.getCurrentPos() == 1 &&
elDofIter.getCurrentElementPos() == ithEdge)
dofs.push_back(elDofIter.getDofPtr());
} while(elDofIter.next());
}
}
void ParallelDomainBase::createDofMemberInfo(DofToPartitions& partitionDofs,
DofContainer& rankOwnedDofs,
DofContainer& rankAllDofs,
......@@ -1876,7 +1760,7 @@ namespace AMDiS {
// Clear all periodic dof mappings calculated before. We do it from scratch.
periodicDof.clear();
MPI::Request request[periodicBoundary.boundary.size() * 2];
MPI::Request request[min(static_cast<int>(periodicBoundary.boundary.size() * 2), 4)];
int requestCounter = 0;
std::vector<int*> sendBuffers, recvBuffers;
......@@ -1895,10 +1779,8 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs, true);
addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
dofs);
boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs, true);
boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
}
// Send the global indices to the rank on the other side.
......@@ -1922,14 +1804,16 @@ namespace AMDiS {
}
MPI::Request::Waitall(requestCounter, request);
MPI::Request::Waitall(requestCounter, request);
// === The rank has received the dofs from the rank on the other side of ===
// === the boundary. Now it can use them to create the mapping between ===
// === the periodic dofs in this rank and the corresponding periodic ===
// === dofs from the other ranks. ===
std::map<DegreeOfFreedom, std::set<int> > dofFromRank;
int i = 0;
for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
it != periodicBoundary.boundary.end(); ++it) {
......@@ -1939,39 +1823,74 @@ namespace AMDiS {
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
DofContainer tmpdofs;
addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
tmpdofs);
addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
tmpdofs, true);
boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, tmpdofs);
boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, tmpdofs, true);
for (int j = static_cast<int>(tmpdofs.size()) - 1; j >= 0; j--)
dofs.push_back(tmpdofs[j]);
}
// Added the received dofs to the mapping.
for (int j = 0; j < static_cast<int>(dofs.size()); j++)
periodicDof[mapLocalGlobalDOFs[*(dofs[j])]].insert(recvBuffers[i][j]);
for (int j = 0; j < static_cast<int>(dofs.size()); j++) {
int globalDofIndex = mapLocalGlobalDOFs[*(dofs[j])];
periodicDof[globalDofIndex].insert(recvBuffers[i][j]);
dofFromRank[globalDofIndex].insert(it->first);
}
delete [] sendBuffers[i];
delete [] recvBuffers[i];
i++;
}
// TODO search for 2 dof mappings.
switch (mpiRank) {
case 0:
periodicDof[0].insert(3136);
break;
case 1:
periodicDof[1024].insert(2080);
break;
case 2:
periodicDof[2080].insert(1024);
break;
case 3:
periodicDof[3136].insert(0);
break;
}
sendBuffers.clear();
recvBuffers.clear();
TEST_EXIT_DBG(mesh->getDim() == 2)
("Periodic boundary corner problem must be generalized to 3d!\n");
requestCounter = 0;
for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin();
it != dofFromRank.end(); ++it) {
if (it->second.size() == 2) {
TEST_EXIT_DBG(periodicDof[it->first].size() == 2)("Missing periodic dof!\n");
int *sendbuf = new int[2];
sendbuf[0] = *(periodicDof[it->first].begin());
sendbuf[1] = *(++(periodicDof[it->first].begin()));
request[requestCounter++] =
mpiComm.Isend(sendbuf, 2, MPI_INT, *(it->second.begin()), 0);
request[requestCounter++] =
mpiComm.Isend(sendbuf, 2, MPI_INT, *(++(it->second.begin())), 0);
sendBuffers.push_back(sendbuf);
int *recvbuf1 = new int[2];
int *recvbuf2 = new int[2];
request[requestCounter++] =
mpiComm.Irecv(recvbuf1, 2, MPI_INT, *(it->second.begin()), 0);
request[requestCounter++] =
mpiComm.Irecv(recvbuf2, 2, MPI_INT, *(++(it->second.begin())), 0);
recvBuffers.push_back(recvbuf1);
recvBuffers.push_back(recvbuf2);
}
}
MPI::Request::Waitall(requestCounter, request);
i = 0;
for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin();
it != dofFromRank.end(); ++it) {
if (it->second.size() == 2) {
for (int k = 0; k < 2; k++)
for (int j = 0; j < 2; j++)
if (recvBuffers[i + k][j] != it->first)
periodicDof[it->first].insert(recvBuffers[i + k][j]);
i++;
}
}
}
......@@ -2133,7 +2052,7 @@ namespace AMDiS {
int nSendInt = rankIt->second.size();
int* buffer = new int[nSendInt];
for (int i = 0; i < nSendInt; i++)
buffer[i] = (rankIt->second)[i].rankObject.elIndex;
buffer[i] = (rankIt->second)[i].rankObj.elIndex;
sendBuffers.push_back(buffer);
......@@ -2166,7 +2085,7 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++) {
int elIndex1 = recvBuffers[bufCounter][i];
int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighbourObject.elIndex;
int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighObj.elIndex;
TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n");
}
......@@ -2371,6 +2290,8 @@ namespace AMDiS {
void ParallelDomainBase::printMapPeriodic(int rank)
{
FUNCNAME("ParallelDomainBase::printMapPeriodic()");
if (rank == -1 || mpiRank == rank) {
std::cout << "====== DOF MAP PERIODIC ====== " << std::endl;
......@@ -2382,12 +2303,14 @@ namespace AMDiS {
std::cout << *dofit << " ";
std::cout << std::endl;
DegreeOfFreedom localdof;
DegreeOfFreedom localdof = -1;
for (DofMapping::iterator dofIt = mapLocalGlobalDOFs.begin();
dofIt != mapLocalGlobalDOFs.end(); ++dofIt)
if (dofIt->second == it->first)
localdof = dofIt->first;
TEST_EXIT(localdof != -1)("There is something wrong!\n");
WorldVector<double> coords;
mesh->getDofIndexCoords(localdof, feSpace, coords);
coords.print();
......