diff --git a/AMDiS/src/InteriorBoundary.cc b/AMDiS/src/InteriorBoundary.cc
index 999c30514084237d9492d63faabd675fead29ee7..4da79698d67505f8bfce5c8748aca88fa1e618df 100644
--- a/AMDiS/src/InteriorBoundary.cc
+++ b/AMDiS/src/InteriorBoundary.cc
@@ -1,4 +1,11 @@
#include "InteriorBoundary.h"
namespace AMDiS {
+
+ AtomicBoundary& InteriorBoundary::getNewAtomicBoundary(int rank)
+ {
+ boundary[rank].resize(boundary[rank].size() + 1);
+ return boundary[rank][boundary[rank].size() - 1];
+ }
+
}
diff --git a/AMDiS/src/InteriorBoundary.h b/AMDiS/src/InteriorBoundary.h
index 43eef37030002e151e09bcfca0d4d1b9eadc9ec9..59f5065cb4317ded13b8324f0eab347c3604738a 100644
--- a/AMDiS/src/InteriorBoundary.h
+++ b/AMDiS/src/InteriorBoundary.h
@@ -23,6 +23,7 @@
#define AMDIS_INTERIORBOUNDARY_H
#include <vector>
+#include <map>
#include "MacroElement.h"
@@ -32,7 +33,7 @@ namespace AMDiS {
struct BoundaryObject {
/// The macro element to which the boundary element corresponds to.
- MacroElement& el;
+ Element* el;
/** \brief
* Defines the geometrical object at the boundary. It must be "a part" of the
@@ -58,13 +59,10 @@ namespace AMDiS {
*/
struct AtomicBoundary {
/// The rank's part of the boundary.
- BoundaryObject& rankObject;
+ BoundaryObject rankObject;
/// The object on the other side of the boundary.
- BoundaryObject& neighbourObject;
-
- /// The rank number of the process at the other side of the boundary.
- int neighbourRank;
+ BoundaryObject neighbourObject;
};
/** \brief
@@ -73,8 +71,12 @@ namespace AMDiS {
*/
class InteriorBoundary {
public:
- protected:
- std::vector<AtomicBoundary> boundary;
+ InteriorBoundary() {}
+
+ AtomicBoundary& getNewAtomicBoundary(int rank);
+
+ public:
+ std::map<int, std::vector<AtomicBoundary> > boundary;
};
}
diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc
index 473bc60f116b7c477e0a479d51d6d0c89fa5e696..5b86db888072ca3cc26d6f48313d5dff5e2b3b09 100644
--- a/AMDiS/src/ParallelDomainProblem.cc
+++ b/AMDiS/src/ParallelDomainProblem.cc
@@ -45,31 +45,26 @@ namespace AMDiS {
// and now partition the mesh
partitionMesh(adaptInfo);
- std::map<int, std::set<int> > partitionDofs;
+ /// === Determine to each dof the set of partitions the dof belongs to. ===
+
+ std::map<const DegreeOfFreedom*, std::set<int> > partitionDofs;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- int nLeaves = 0;
while (elInfo) {
Element *element = elInfo->getElement();
- // Hidde elements which are not part of ranks partition.
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
- if (partitionData->getPartitionStatus() != IN) {
- } else {
- }
-
// Determine to each dof the partition(s) it corresponds to.
for (int i = 0; i < 3; i++)
- partitionDofs[element->getDOF(i, 0)].insert(partitionVec[element->getIndex()]);
+ partitionDofs[element->getDOF(i)].insert(partitionVec[element->getIndex()]);
elInfo = stack.traverseNext(elInfo);
}
- std::vector<int> rankDofs;
- for (std::map<int, std::set<int> >::iterator it = partitionDofs.begin();
+ /// === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
+
+ std::vector<const DegreeOfFreedom*> rankDofs;
+ for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDofs.begin();
it != partitionDofs.end();
++it) {
for (std::set<int>::iterator itpart1 = it->second.begin();
@@ -79,30 +74,64 @@ namespace AMDiS {
if (it->second.size() == 1) {
rankDofs.push_back(it->first);
} else {
+ // This dof is at the ranks boundary. It is owned by the rank only if
+ // the rank number is the highest of all ranks containing this dof.
+
bool insert = true;
+ int highestRank = mpiRank;
for (std::set<int>::iterator itpart2 = it->second.begin();
itpart2 != it->second.end();
++itpart2) {
- if (*itpart2 > mpiRank) {
+ if (*itpart2 > mpiRank)
insert = false;
- break;
- }
+
+ if (*itpart2 > highestRank)
+ highestRank = *itpart2;
}
- if (insert) {
+
+ if (insert)
rankDofs.push_back(it->first);
- }
+
+ boundaryDofs[it->first] = highestRank;
}
}
}
}
- if (mpiRank == 1) {
- std::cout << "RANKS dofs = ";
- for (int i = 0; i < rankDofs.size(); i++)
- std::cout << rankDofs[i] << " ";
- std::cout << std::endl;
+
+ // === Create interior boundary information ===
+
+ elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
+ while (elInfo) {
+ Element *element = elInfo->getElement();
+
+ // Hidde elements which are not part of ranks partition.
+ PartitionElementData *partitionData =
+ dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
+ if (partitionData->getPartitionStatus() == IN) {
+ for (int i = 0; i < 3; i++) {
+ if (!elInfo->getNeighbour(i))
+ continue;
+
+ PartitionElementData *neighbourPartitionData =
+ dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED));
+ if (neighbourPartitionData->getPartitionStatus() == OUT) {
+ AtomicBoundary& bound = interiorBoundary.
+ getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]);
+ bound.rankObject.el = element;
+ bound.rankObject.subObjAtBoundary = EDGE;
+ bound.rankObject.ithObjAtBoundary = i;
+ bound.neighbourObject.el = elInfo->getNeighbour(i);
+ bound.neighbourObject.subObjAtBoundary = EDGE;
+ bound.neighbourObject.ithObjAtBoundary = -1;
+ }
+ }
+ }
+
+ elInfo = stack.traverseNext(elInfo);
}
+
// === Remove all macro elements that are not part of the rank partition. ===
std::vector<MacroElement*> macrosToRemove;
@@ -124,15 +153,15 @@ namespace AMDiS {
int *gOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));
int *lOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));
- for (std::vector<int>::iterator it = rankDofs.begin();
+ for (std::vector<const DegreeOfFreedom*>::iterator it = rankDofs.begin();
it != rankDofs.end(); ++it) {
- gOrder[nRankDOFs++] = *it;
+ gOrder[nRankDOFs++] = (*it)[0];
}
int rstart = 0;
MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
rstart -= nRankDOFs;
-
+
for (int i = 0; i < nRankDOFs; i++) {
lOrder[i] = rstart + i;
}
@@ -141,6 +170,103 @@ namespace AMDiS {
free(gOrder);
free(lOrder);
+
+ /// === Create information which dof indices must be send and which must be received. ===
+
+ std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> > sendNewDofs;
+ std::map<int, std::vector<DegreeOfFreedom> > recvNewDofs;
+
+ for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDofs.begin();
+ it != boundaryDofs.end();
+ ++it) {
+ if (it->second == mpiRank) {
+ int oldDofIndex = (it->first)[0];
+ int newDofIndex = 0;
+ for (int i = 0; i < static_cast<int>(rankDofs.size()); i++) {
+ if (rankDofs[i] == it->first) {
+ newDofIndex = rstart + i;
+ break;
+ }
+ }
+
+ for (std::set<int>::iterator itRanks = partitionDofs[it->first].begin();
+ itRanks != partitionDofs[it->first].end();
+ ++itRanks) {
+ if (*itRanks != mpiRank) {
+ sendNewDofs[*itRanks][oldDofIndex] = newDofIndex;
+ }
+ }
+ } else {
+ recvNewDofs[it->second].push_back((it->first)[0]);
+ }
+ }
+
+ /// === Send and receive the dof indices at boundary. ===
+
+ std::vector<int*> sendBuffers(sendNewDofs.size());
+ std::vector<int*> recvBuffers(recvNewDofs.size());
+
+ int i = 0;
+ for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
+ sendIt != sendNewDofs.end();
+ ++sendIt, i++) {
+ sendBuffers[i] = new int[sendIt->second.size() * 2];
+ int c = 0;
+ for (std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator dofIt = sendIt->second.begin();
+ dofIt != sendIt->second.end();
+ ++dofIt, c += 2) {
+ sendBuffers[i][c] = dofIt->first;
+ sendBuffers[i][c + 1] = dofIt->second;
+ }
+
+ mpiComm.Isend(sendBuffers[i], sendIt->second.size() * 2, MPI_INT, sendIt->first, 0);
+ }
+
+ i = 0;
+ for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
+ recvIt != recvNewDofs.end();
+ ++recvIt, i++) {
+ recvBuffers[i] = new int[recvIt->second.size() * 2];
+
+ mpiComm.Irecv(recvBuffers[i], recvIt->second.size() * 2, MPI_INT, recvIt->first, 0);
+ }
+
+
+ mpiComm.Barrier();
+
+ /// === Change dof indices at boundary from other ranks. ===
+
+ i = 0;
+ for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
+ recvIt != recvNewDofs.end();
+ ++recvIt, i++) {
+
+ for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) {
+ for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDofs.begin();
+ dofIt != boundaryDofs.end();
+ ++dofIt) {
+ if ((dofIt->first)[0] == recvBuffers[i][j * 2]) {
+ const_cast<DegreeOfFreedom*>(dofIt->first)[0] = recvBuffers[i][j * 2 + 1];
+ break;
+ }
+ }
+ }
+
+ delete [] recvBuffers[i];
+ }
+
+ i = 0;
+ for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
+ sendIt != sendNewDofs.end();
+ ++sendIt, i++) {
+ delete [] sendBuffers[i];
+ }
+
+ /// === Change dof indices for rank partition. ===
+
+ for (int i = 0; i < static_cast<int>(rankDofs.size()); i++) {
+ const_cast<DegreeOfFreedom*>(rankDofs[i])[0] = rstart + i;
+ }
}
void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h
index ac369330bff2824c79dd01e916f802e2871166fd..fe206322637bc65cabfb0d469108f47b5d92d037 100644
--- a/AMDiS/src/ParallelDomainProblem.h
+++ b/AMDiS/src/ParallelDomainProblem.h
@@ -29,6 +29,7 @@
#include "ProblemIterationInterface.h"
#include "FiniteElemSpace.h"
#include "AdaptInfo.h"
+#include "InteriorBoundary.h"
#include "petscao.h"
#include "mpi.h"
@@ -156,6 +157,14 @@ namespace AMDiS {
/// Number of DOFs in the rank mesh.
int nRankDOFs;
+
+ /** \brief
+ * Defines the interioir boundaries of the domain that result from partitioning
+ * the whole mesh.
+ */
+ InteriorBoundary interiorBoundary;
+
+ std::map<const DegreeOfFreedom*, int> boundaryDofs;
};
class ParallelDomainProblemScal : public ParallelDomainProblemBase