From 32b711485251f9b4076aee6dda1f7f836c9879a4 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 5 Jan 2011 11:13:08 +0000
Subject: [PATCH] Fixed some bugs for periodic BCs in 3D in parrallel
computing.
---
AMDiS/src/ElInfo.cc | 1 +
AMDiS/src/parallel/InteriorBoundary.cc | 32 +-
AMDiS/src/parallel/InteriorBoundary.h | 4 +-
AMDiS/src/parallel/MeshDistributor.cc | 442 +++++++++++++-----
AMDiS/src/parallel/MeshDistributor.h | 4 -
AMDiS/src/parallel/ParMetisPartitioner.cc | 2 +
AMDiS/src/parallel/ParallelDebug.cc | 13 +-
AMDiS/src/parallel/ParallelProblemStatBase.cc | 10 +-
AMDiS/src/parallel/PetscSolver.cc | 7 +-
AMDiS/src/parallel/StdMpi.cc | 26 +-
10 files changed, 378 insertions(+), 163 deletions(-)
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index fb213ef4..439dd568 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -172,6 +172,7 @@ namespace AMDiS {
}
}
+
BoundaryType ElInfo::getBoundary(GeoIndex pos, int i)
{
static int indexOffset[3][3] = {
diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc
index e483ba7e..612d82d7 100644
--- a/AMDiS/src/parallel/InteriorBoundary.cc
+++ b/AMDiS/src/parallel/InteriorBoundary.cc
@@ -14,11 +14,13 @@
#include "FiniteElemSpace.h"
#include "BasisFunction.h"
#include "Serializer.h"
+#include "VertexVector.h"
namespace AMDiS {
void BoundaryObject::setReverseMode(BoundaryObject &otherBound,
- FiniteElemSpace *feSpace)
+ FiniteElemSpace *feSpace,
+ BoundaryType boundary)
{
FUNCNAME("BoundaryObject::setReverseMode()");
@@ -31,7 +33,8 @@ namespace AMDiS {
case 3:
TEST_EXIT_DBG(otherBound.elType == 0)
- ("Only 3D macro elements with level 0 are supported!\n");
+ ("Only 3D macro elements with level 0 are supported. This element has level %d!\n", otherBound.elType);
+
if (subObj == EDGE) {
int el0_v0 = el->getVertexOfEdge(ithObj, 0);
@@ -45,15 +48,22 @@ namespace AMDiS {
basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0);
basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1);
- TEST_EXIT_DBG(localDofs0[el0_v0] == localDofs1[el1_v0] ||
- localDofs0[el0_v0] == localDofs1[el1_v1])
- ("This should not happen!\n");
- TEST_EXIT_DBG(localDofs0[el0_v1] == localDofs1[el1_v0] ||
- localDofs0[el0_v1] == localDofs1[el1_v1])
- ("This should not happen!\n");
-
- if (localDofs0[el0_v0] != localDofs1[el1_v0])
- otherMode = true;
+ Mesh *mesh = feSpace->getMesh();
+
+ if (mesh->isPeriodicAssociation(boundary) == false) {
+ TEST_EXIT_DBG(localDofs0[el0_v0] == localDofs1[el1_v0] ||
+ localDofs0[el0_v0] == localDofs1[el1_v1])
+ ("This should not happen!\n");
+ TEST_EXIT_DBG(localDofs0[el0_v1] == localDofs1[el1_v0] ||
+ localDofs0[el0_v1] == localDofs1[el1_v1])
+ ("This should not happen!\n");
+
+ if (localDofs0[el0_v0] != localDofs1[el1_v0])
+ otherMode = true;
+ } else {
+ if (mesh->associated(localDofs0[el0_v0], localDofs1[el1_v0]) == false)
+ otherMode = true;
+ }
}
if (subObj == FACE && ithObj != 1) {
diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h
index ea12b82f..81280de0 100644
--- a/AMDiS/src/parallel/InteriorBoundary.h
+++ b/AMDiS/src/parallel/InteriorBoundary.h
@@ -54,7 +54,9 @@ namespace AMDiS {
excludedSubstructures(0)
{}
- void setReverseMode(BoundaryObject &otherBound, FiniteElemSpace *feSpace);
+ void setReverseMode(BoundaryObject &otherBound,
+ FiniteElemSpace *feSpace,
+ BoundaryType boundary);
bool operator==(const BoundaryObject& other) const;
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 6235c168..a1762416 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -221,6 +221,7 @@ namespace AMDiS {
createPeriodicMap();
+
// === Global refinements. ===
int globalRefinement = 0;
@@ -544,10 +545,15 @@ namespace AMDiS {
(mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
- for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
- if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
- (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
- allBound[it.getRank()].push_back(*it);
+ for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
+ if (it.getRank() == mpiRank) {
+ // ERROR_EXIT("Na, du weisst schon!\n");
+ } else {
+ if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
+ (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
+ allBound[it.getRank()].push_back(*it);
+ }
+ }
// === Check the boundaries and adapt mesh if necessary. ===
@@ -579,8 +585,7 @@ namespace AMDiS {
// === Update periodic mapping, if there are periodic boundaries. ===
- createPeriodicMap();
-
+ createPeriodicMap();
INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n",
@@ -1256,6 +1261,8 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::createInteriorBoundaryInfo()");
+ MSG("CREATE BOUNDARY INFO!\n");
+
createMeshElementData();
createBoundaryData();
@@ -1267,6 +1274,8 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::updateInteriorBoundaryInfo()");
+ MSG("UPDATE BOUNDARY INFO!\n");
+
elObjects.createRankData(partitionVec);
createBoundaryData();
@@ -1277,6 +1286,13 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::createMeshElementData()");
+ // Stores to each vertex all its periodic associations.
+ std::map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
+
+ // Stores to each edge all its periodic associations.
+ std::map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc;
+
+
// === Fills macro element data structures. ===
TraverseStack stack;
@@ -1301,43 +1317,68 @@ namespace AMDiS {
switch (mesh->getDim()) {
case 2:
for (int i = 0; i < el->getGeo(EDGE); i++) {
- if (mesh->isPeriodicAssociation(elInfo->getBoundary(i))) {
+ if (mesh->isPeriodicAssociation(elInfo->getBoundary(EDGE, i))) {
+ Element *neigh = elInfo->getNeighbour(i);
+
DofEdge edge1 = el->getEdge(i);
- DofEdge edge2 = elInfo->getNeighbour(i)->getEdge(elInfo->getOppVertex(i));
- BoundaryType boundaryType = elInfo->getBoundary(i);
+ DofEdge edge2 = neigh->getEdge(elInfo->getOppVertex(i));
+ BoundaryType boundaryType = elInfo->getBoundary(EDGE, i);
periodicEdges[std::make_pair(edge1, edge2)] = boundaryType;
+ periodicEdgeAssoc[edge1].insert(edge2);
- for (int j = 0; j < 2; j++) {
- int ithDofRankObj = el->getVertexOfEdge(i, j);
- DegreeOfFreedom dof = el->getDof(ithDofRankObj, 0);
-
- DegreeOfFreedom neighDof = -1;
- int ithDofNeighObj = -1;
- for (int k = 0; k < 3; k++) {
- if (elInfo->getNeighbour(i)->getDof(k, 0) ==
- mesh->getPeriodicAssociations(boundaryType)[dof]) {
- neighDof = elInfo->getNeighbour(i)->getDof(k, 0);
- ithDofNeighObj = k;
- break;
- }
- }
- TEST_EXIT_DBG(neighDof > -1)("Should not happen!\n");
-
- periodicDofs[std::make_pair(dof, neighDof)] = boundaryType;
- periodicDofAssoc[dof].insert(boundaryType);
- }
+ periodicDofs[std::make_pair(edge1.first, edge2.first)] = boundaryType;
+ periodicDofs[std::make_pair(edge1.second, edge2.second)] = boundaryType;
+
+ periodicDofAssoc[edge1.first].insert(boundaryType);
+ periodicDofAssoc[edge1.second].insert(boundaryType);
+
+ TEST_EXIT_DBG(edge1.first == mesh->getPeriodicAssociations(boundaryType)[edge2.first] &&
+ edge1.second == mesh->getPeriodicAssociations(boundaryType)[edge2.second])
+ ("Should not happen!\n");
}
}
break;
case 3:
for (int i = 0; i < el->getGeo(FACE); i++) {
- if (mesh->isPeriodicAssociation(elInfo->getBoundary(i))) {
+ if (mesh->isPeriodicAssociation(elInfo->getBoundary(FACE, i))) {
+ Element *neigh = elInfo->getNeighbour(i);
+
DofFace face1 = el->getFace(i);
- DofFace face2 = elInfo->getNeighbour(i)->getFace(elInfo->getOppVertex(i));
-
+ DofFace face2 = neigh->getFace(elInfo->getOppVertex(i));
+ BoundaryType boundaryType = elInfo->getBoundary(FACE, i);
+
periodicFaces[std::make_pair(face1, face2)] = elInfo->getBoundary(i);
+
+ periodicDofs[std::make_pair(face1.get<0>(), face2.get<0>())] = boundaryType;
+ periodicDofs[std::make_pair(face1.get<1>(), face2.get<1>())] = boundaryType;
+ periodicDofs[std::make_pair(face1.get<2>(), face2.get<2>())] = boundaryType;
+
+ periodicDofAssoc[face1.get<0>()].insert(boundaryType);
+ periodicDofAssoc[face1.get<1>()].insert(boundaryType);
+ periodicDofAssoc[face1.get<2>()].insert(boundaryType);
+
+ TEST_EXIT_DBG(face1.get<0>() == mesh->getPeriodicAssociations(boundaryType)[face2.get<0>()] &&
+ face1.get<1>() == mesh->getPeriodicAssociations(boundaryType)[face2.get<1>()] &&
+ face1.get<2>() == mesh->getPeriodicAssociations(boundaryType)[face2.get<2>()])
+ ("Should not happen!\n");
+
+
+ DofEdge elEdge1 = std::make_pair(face1.get<0>(), face1.get<1>());
+ DofEdge elEdge2 = std::make_pair(face1.get<0>(), face1.get<2>());
+ DofEdge elEdge3 = std::make_pair(face1.get<1>(), face1.get<2>());
+ DofEdge neighEdge1 = std::make_pair(face2.get<0>(), face2.get<1>());
+ DofEdge neighEdge2 = std::make_pair(face2.get<0>(), face2.get<2>());
+ DofEdge neighEdge3 = std::make_pair(face2.get<1>(), face2.get<2>());
+
+ periodicEdges[std::make_pair(elEdge1, neighEdge1)] = boundaryType;
+ periodicEdges[std::make_pair(elEdge2, neighEdge2)] = boundaryType;
+ periodicEdges[std::make_pair(elEdge3, neighEdge3)] = boundaryType;
+
+ periodicEdgeAssoc[elEdge1].insert(neighEdge1);
+ periodicEdgeAssoc[elEdge2].insert(neighEdge2);
+ periodicEdgeAssoc[elEdge3].insert(neighEdge3);
}
}
break;
@@ -1355,35 +1396,90 @@ namespace AMDiS {
// === Search for interectly connected vertices in periodic boundaries. ===
if (periodicDofs.size() > 0) {
- TEST_EXIT(mesh->getDim() == 2)("Parallel boundaries not yet supported in 3D!\n");
- std::vector<DegreeOfFreedom> twoPeriodicDof;
+
+ // === Search for an unsed boundary index. ===
+
+ BoundaryType newPeriodicBoundaryType = 0;
+ for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
+ it != mesh->getPeriodicAssociations().end(); ++it)
+ newPeriodicBoundaryType = min(newPeriodicBoundaryType, it->first);
+
+ TEST_EXIT_DBG(newPeriodicBoundaryType < 0)("Should not happen!\n");
+ newPeriodicBoundaryType--;
+ mesh->getPeriodicAssociations()[newPeriodicBoundaryType] =
+ new VertexVector(feSpace->getAdmin(), "");
+
+
+ // === Get all vertex DOFs that have multiple periodic associations. ===
+
+ std::vector<DegreeOfFreedom> multiplePeriodicDof;
for (std::map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
it != periodicDofAssoc.end(); ++it) {
- TEST_EXIT_DBG(it->second.size() < 3)("Should not happen!\n");
+ if (mesh->getDim() == 2) {
+ TEST_EXIT_DBG(it->second.size() <= 2)("Should not happen!\n");
+ }
+ if (mesh->getDim() == 3) {
+ TEST_EXIT_DBG(it->second.size() <= 3)("Should not happen!\n");
+ }
+
- if (it->second.size() == 2)
- twoPeriodicDof.push_back(it->first);
+ if ((mesh->getDim() == 2 && it->second.size() == 2) ||
+ (mesh->getDim() == 3 && it->second.size() == 3))
+ multiplePeriodicDof.push_back(it->first);
}
- TEST_EXIT_DBG(twoPeriodicDof.size() == 0 || twoPeriodicDof.size() == 4)
- ("Should not happen (%d)!\n", twoPeriodicDof.size());
- if (twoPeriodicDof.size() == 4) {
- for (int i = 0; i < 4; i++) {
- for (int j = 0; j < 4; j++) {
- if (i == j)
- continue;
-
- std::pair<DegreeOfFreedom, DegreeOfFreedom> perDofs =
- std::make_pair(twoPeriodicDof[i], twoPeriodicDof[j]);
+ if (mesh->getDim() == 2) {
+ TEST_EXIT_DBG(multiplePeriodicDof.size() == 0 ||
+ multiplePeriodicDof.size() == 4)
+ ("Should not happen (%d)!\n", multiplePeriodicDof.size());
+ }
- if (periodicDofs.count(perDofs) == 0)
- periodicDofs[perDofs] = 3;
- }
+ if (mesh->getDim() == 3) {
+ TEST_EXIT_DBG(multiplePeriodicDof.size() == 0 ||
+ multiplePeriodicDof.size() == 8)
+ ("Should not happen (%d)!\n", multiplePeriodicDof.size());
+ }
+
+
+ int nMultiplePeriodicDofs = multiplePeriodicDof.size();
+ for (int i = 0; i < nMultiplePeriodicDofs; i++) {
+ for (int j = 0; j < nMultiplePeriodicDofs; j++) {
+ if (i == j)
+ continue;
+
+ std::pair<DegreeOfFreedom, DegreeOfFreedom> perDofs =
+ std::make_pair(multiplePeriodicDof[i], multiplePeriodicDof[j]);
+
+ if (periodicDofs.count(perDofs) == 0)
+ periodicDofs[perDofs] = newPeriodicBoundaryType;
+ }
+ }
+
+
+
+ // === Get all edges that have multiple periodic associations (3D only!). ===
+
+ for (std::map<DofEdge, std::set<DofEdge> >::iterator it = periodicEdgeAssoc.begin();
+ it != periodicEdgeAssoc.end(); ++it) {
+ if (it->second.size() > 1) {
+ TEST_EXIT_DBG(mesh->getDim() == 3)("Should not happen!\n");
+ TEST_EXIT_DBG(it->second.size() == 2)("Should not happen!\n");
+
+ std::pair<DofEdge, DofEdge> perEdge0 =
+ std::make_pair(*(it->second.begin()), *(++(it->second.begin())));
+ std::pair<DofEdge, DofEdge> perEdge1 =
+ std::make_pair(perEdge0.second, perEdge0.first);
+
+ periodicEdges[perEdge0] = newPeriodicBoundaryType;
+ periodicEdges[perEdge1] = newPeriodicBoundaryType;
}
}
+
+
}
+
}
@@ -1450,7 +1546,7 @@ namespace AMDiS {
AtomicBoundary& b = myIntBoundary.getNewAtomic(it2->first);
b = bound;
- b.rankObj.setReverseMode(b.neighObj, feSpace);
+ b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
}
} else {
@@ -1473,7 +1569,7 @@ namespace AMDiS {
AtomicBoundary& b = otherIntBoundary.getNewAtomic(owner);
b = bound;
- b.neighObj.setReverseMode(b.rankObj, feSpace);
+ b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
}
}
}
@@ -1487,11 +1583,12 @@ namespace AMDiS {
if (elObjects.isInRank(it->first.first, mpiRank) == false)
continue;
+ ElementObjectData& perDofEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
+
for (std::map<int, ElementObjectData>::iterator elIt = elObjects.getElementsInRank(it->first.second).begin();
elIt != elObjects.getElementsInRank(it->first.second).end(); ++elIt) {
int otherElementRank = elIt->first;
- ElementObjectData& perDofEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
ElementObjectData& perDofEl1 = elIt->second;
AtomicBoundary bound;
@@ -1503,7 +1600,7 @@ namespace AMDiS {
bound.neighObj.el = macroElIndexMap[perDofEl1.elIndex];
bound.neighObj.elIndex = perDofEl1.elIndex;
- bound.neighObj.elType = -1;
+ bound.neighObj.elType = macroElIndexTypeMap[perDofEl1.elIndex];
bound.neighObj.subObj = VERTEX;
bound.neighObj.ithObj = perDofEl1.ithObject;
@@ -1515,50 +1612,88 @@ namespace AMDiS {
}
-
for (std::map<std::pair<DofEdge, DofEdge>, BoundaryType>::iterator it = periodicEdges.begin();
it != periodicEdges.end(); ++it) {
- int perEdgeOwner0 = elObjects.getOwner(it->first.first);
- int perEdgeOwner1 = elObjects.getOwner(it->first.second);
+ if (elObjects.isInRank(it->first.first, mpiRank) == false)
+ continue;
+
+ ElementObjectData& perEdgeEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
+
+ for (std::map<int, ElementObjectData>::iterator elIt = elObjects.getElementsInRank(it->first.second).begin();
+ elIt != elObjects.getElementsInRank(it->first.second).end(); ++elIt) {
+
+ int otherElementRank = elIt->first;
+ ElementObjectData& perEdgeEl1 = elIt->second;
+
+ AtomicBoundary bound;
+ bound.rankObj.el = macroElIndexMap[perEdgeEl0.elIndex];
+ bound.rankObj.elIndex = perEdgeEl0.elIndex;
+ bound.rankObj.elType = macroElIndexTypeMap[perEdgeEl0.elIndex];
+ bound.rankObj.subObj = EDGE;
+ bound.rankObj.ithObj = perEdgeEl0.ithObject;
+
+ bound.neighObj.el = macroElIndexMap[perEdgeEl1.elIndex];
+ bound.neighObj.elIndex = perEdgeEl1.elIndex;
+ bound.neighObj.elType = macroElIndexTypeMap[perEdgeEl1.elIndex];
+ bound.neighObj.subObj = EDGE;
+ bound.neighObj.ithObj = perEdgeEl1.ithObject;
+
+ bound.type = it->second;
+
+ AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
+ b = bound;
+
+ if (mpiRank > otherElementRank)
+ b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
+ else
+ b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
+ }
+ }
- if (perEdgeOwner0 != mpiRank)
+
+ for (std::map<std::pair<DofFace, DofFace>, BoundaryType>::iterator it = periodicFaces.begin();
+ it != periodicFaces.end(); ++it) {
+ if (elObjects.isInRank(it->first.first, mpiRank) == false)
continue;
- TEST_EXIT_DBG(perEdgeOwner0 != perEdgeOwner1)("Should not happen!\n");
TEST_EXIT_DBG(elObjects.getElements(it->first.first).size() == 1)
- ("Should not happen!\n");
+ ("Should not happen!\n");
TEST_EXIT_DBG(elObjects.getElements(it->first.second).size() == 1)
- ("Should not happen!\n");
-
- int otherElementRank = perEdgeOwner1;
- ElementObjectData& perEdgeEl0 = elObjects.getElements(it->first.first)[0];
- ElementObjectData& perEdgeEl1 = elObjects.getElements(it->first.second)[0];
-
- AtomicBoundary bound;
- bound.rankObj.el = macroElIndexMap[perEdgeEl0.elIndex];
- bound.rankObj.elIndex = perEdgeEl0.elIndex;
- bound.rankObj.elType = macroElIndexTypeMap[perEdgeEl0.elIndex];
- bound.rankObj.subObj = EDGE;
- bound.rankObj.ithObj = perEdgeEl0.ithObject;
-
- bound.neighObj.el = macroElIndexMap[perEdgeEl1.elIndex];
- bound.neighObj.elIndex = perEdgeEl1.elIndex;
- bound.neighObj.elType = -1;
- bound.neighObj.subObj = EDGE;
- bound.neighObj.ithObj = perEdgeEl1.ithObject;
-
- bound.type = it->second;
-
- AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
- b = bound;
-
- if (mpiRank > otherElementRank)
- b.rankObj.setReverseMode(b.neighObj, feSpace);
- else
- b.neighObj.setReverseMode(b.rankObj, feSpace);
- }
+ ("Should not happen!\n");
+ ElementObjectData& perEdgeEl0 = elObjects.getElementsInRank(it->first.first)[mpiRank];
+ for (std::map<int, ElementObjectData>::iterator elIt = elObjects.getElementsInRank(it->first.second).begin();
+ elIt != elObjects.getElementsInRank(it->first.second).end(); ++elIt) {
+
+ int otherElementRank = elIt->first;
+ ElementObjectData& perEdgeEl1 = elIt->second;
+
+ AtomicBoundary bound;
+ bound.rankObj.el = macroElIndexMap[perEdgeEl0.elIndex];
+ bound.rankObj.elIndex = perEdgeEl0.elIndex;
+ bound.rankObj.elType = macroElIndexTypeMap[perEdgeEl0.elIndex];
+ bound.rankObj.subObj = FACE;
+ bound.rankObj.ithObj = perEdgeEl0.ithObject;
+
+ bound.neighObj.el = macroElIndexMap[perEdgeEl1.elIndex];
+ bound.neighObj.elIndex = perEdgeEl1.elIndex;
+ bound.neighObj.elType = macroElIndexTypeMap[perEdgeEl1.elIndex];
+ bound.neighObj.subObj = FACE;
+ bound.neighObj.ithObj = perEdgeEl1.ithObject;
+
+ bound.type = it->second;
+
+ AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
+ b = bound;
+
+ if (mpiRank > otherElementRank)
+ b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
+ else
+ b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
+ }
+ }
+
// === Once we have this information, we must care about the order of the atomic ===
// === bounds in the three boundary handling object. Eventually all the bound- ===
@@ -1569,6 +1704,7 @@ namespace AMDiS {
stdMpi.recv(otherIntBoundary.boundary);
stdMpi.startCommunication<int>(MPI_INT);
+
// === The information about all neighbouring boundaries has been received. So ===
// === the rank tests if its own atomic boundaries are in the same order. If ===
// === not, the atomic boundaries are swaped to the correct order. ===
@@ -1614,43 +1750,66 @@ namespace AMDiS {
for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
rankIt != periodicBoundary.boundary.end(); ++rankIt) {
- TEST_EXIT_DBG(rankIt->first != mpiRank)
- ("It is not allowed to have an interior boundary within a rank partition!\n");
+ if (rankIt->first == mpiRank)
+ continue;
if (rankIt->first < mpiRank)
sendBounds.boundary[rankIt->first] = rankIt->second;
else
- recvBounds.boundary[rankIt->first] = rankIt->second;
+ recvBounds.boundary[rankIt->first] = rankIt->second;
}
stdMpi.send(sendBounds.boundary);
stdMpi.recv(recvBounds.boundary);
stdMpi.startCommunication<int>(MPI_INT);
+ /*
for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
rankIt != periodicBoundary.boundary.end(); ++rankIt) {
- if (rankIt->first <= mpiRank)
- continue;
-
+
for (unsigned int j = 0; j < rankIt->second.size(); j++) {
-
- BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj;
-
- if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvedBound) {
+
+ MSG("%-ith boundary with rank %d\n", j, rankIt->first);
+ MSG(" bound: el = %d subobj = %d ithobj = %d\n",
+ periodicBoundary.boundary[rankIt->first][j].rankObj.elIndex,
+ periodicBoundary.boundary[rankIt->first][j].rankObj.subObj,
+ periodicBoundary.boundary[rankIt->first][j].rankObj.ithObj);
+ MSG(" neigh: el = %d subobj = %d ithobj = %d\n",
+ periodicBoundary.boundary[rankIt->first][j].neighObj.elIndex,
+ periodicBoundary.boundary[rankIt->first][j].neighObj.subObj,
+ periodicBoundary.boundary[rankIt->first][j].neighObj.ithObj);
+ }
+ }
+ */
+
+
+ for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
+ rankIt != periodicBoundary.boundary.end(); ++rankIt) {
+
+ if (rankIt->first <= mpiRank)
+ continue;
+
+ for (unsigned int j = 0; j < rankIt->second.size(); j++) {
+ BoundaryObject &recvRankObj = stdMpi.getRecvData()[rankIt->first][j].rankObj;
+ BoundaryObject &recvNeighObj = stdMpi.getRecvData()[rankIt->first][j].neighObj;
+
+ if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvRankObj ||
+ periodicBoundary.boundary[rankIt->first][j].rankObj != recvNeighObj) {
unsigned int k = j + 1;
for (; k < rankIt->second.size(); k++)
- if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvedBound)
+ if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvRankObj &&
+ periodicBoundary.boundary[rankIt->first][k].rankObj == recvNeighObj)
break;
// The element must always be found, because the list is just in
// another order.
TEST_EXIT_DBG(k < rankIt->second.size())("Should never happen!\n");
-
+
// Swap the current with the found element.
AtomicBoundary tmpBound = (rankIt->second)[k];
(rankIt->second)[k] = (rankIt->second)[j];
(rankIt->second)[j] = tmpBound;
- }
+ }
}
}
} // periodicBoundary.boundary.size() > 0
@@ -1873,29 +2032,56 @@ namespace AMDiS {
for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
it != periodicBoundary.boundary.end(); ++it) {
- TEST_EXIT_DBG(it->first != mpiRank)
- ("Periodic interior boundary within the rank itself is not possible!\n");
- // Create dof indices on the boundary.
- DofContainer& dofs = rankPeriodicDofs[it->first];
- for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
- boundIt != it->second.end(); ++boundIt) {
- int nDofs = dofs.size();
+ if (it->first == mpiRank) {
+ TEST_EXIT_DBG(it->second.size() % 2 == 0)("Should not happen!\n");
- boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs);
- boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs);
+ for (unsigned int i = 0; i < it->second.size(); i++) {
+ AtomicBoundary &bound = it->second[i];
+ DofContainer dofs0, dofs1;
+ bound.rankObj.el->getVertexDofs(feSpace, bound.rankObj, dofs0);
+ bound.rankObj.el->getNonVertexDofs(feSpace, bound.rankObj, dofs0);
+ bound.neighObj.el->getVertexDofs(feSpace, bound.neighObj, dofs1);
+ bound.neighObj.el->getNonVertexDofs(feSpace, bound.neighObj, dofs1);
- for (unsigned int i = 0; i < (dofs.size() - nDofs); i++)
- rankToDofType[it->first].push_back(boundIt->type);
- }
+ TEST_EXIT_DBG(dofs0.size() == dofs1.size())("Should not happen!\n");
- // Send the global indices to the rank on the other side.
- stdMpi.getSendData(it->first).reserve(dofs.size());
- for (unsigned int i = 0; i < dofs.size(); i++)
- stdMpi.getSendData(it->first).push_back(mapLocalGlobalDofs[*(dofs[i])]);
+ BoundaryType type = bound.type;
+
+ for (unsigned int j = 0; j < dofs0.size(); j++) {
+ DegreeOfFreedom globalDof0 = mapLocalGlobalDofs[*(dofs0[j])];
+ DegreeOfFreedom globalDof1 = mapLocalGlobalDofs[*(dofs1[j])];
+
+ if (periodicDofAssociations[globalDof0].count(type) == 0) {
+ periodicDof[type][globalDof0] = globalDof1;
+ periodicDofAssociations[globalDof0].insert(type);
+ }
+ }
+ }
+
+ } else {
+ // Create DOF indices on the boundary.
+ DofContainer& dofs = rankPeriodicDofs[it->first];
+ for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+ boundIt != it->second.end(); ++boundIt) {
- // Receive from this rank the same number of dofs.
- stdMpi.recv(it->first, dofs.size());
+ int nDofs = dofs.size();
+
+ boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs);
+ boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs);
+
+ for (unsigned int i = 0; i < (dofs.size() - nDofs); i++)
+ rankToDofType[it->first].push_back(boundIt->type);
+ }
+
+ // Send the global indices to the rank on the other side.
+ stdMpi.getSendData(it->first).reserve(dofs.size());
+ for (unsigned int i = 0; i < dofs.size(); i++)
+ stdMpi.getSendData(it->first).push_back(mapLocalGlobalDofs[*(dofs[i])]);
+
+ // Receive from this rank the same number of dofs.
+ stdMpi.recv(it->first, dofs.size());
+ }
}
stdMpi.updateSendDataSize();
@@ -1908,8 +2094,6 @@ namespace AMDiS {
// === dofs from the other ranks. ===
- std::map<DegreeOfFreedom, std::set<int> > dofFromRank;
-
for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
it != periodicBoundary.boundary.end(); ++it) {
DofContainer& dofs = rankPeriodicDofs[it->first];
@@ -1928,19 +2112,21 @@ namespace AMDiS {
if (periodicDofAssociations[globalDofIndex].count(type) == 0) {
periodicDof[type][globalDofIndex] = mapGlobalDofIndex;
periodicDofAssociations[globalDofIndex].insert(type);
- dofFromRank[globalDofIndex].insert(it->first);
}
}
}
#if (DEBUG != 0)
- for (std::map<int, std::set<BoundaryType> >::iterator it = periodicDofAssociations.begin();
- it != periodicDofAssociations.end(); ++it) {
- int nAssoc = it->second.size();
- TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3)
- ("Should not happen! DOF %d has %d periodic associations!\n",
- it->first, nAssoc);
+ // In 2D, a periodic DOF can have either 1 or 3 associations. Check this!
+ if (mesh->getDim() == 2) {
+ for (std::map<int, std::set<BoundaryType> >::iterator it = periodicDofAssociations.begin();
+ it != periodicDofAssociations.end(); ++it) {
+ int nAssoc = it->second.size();
+ TEST_EXIT_DBG(nAssoc == 1 || nAssoc == 3)
+ ("Should not happen! DOF %d has %d periodic associations!\n",
+ it->first, nAssoc);
+ }
}
#endif
}
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index d4597e8c..7ad15f37 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -483,10 +483,6 @@ namespace AMDiS {
std::map<std::pair<DofEdge, DofEdge>, BoundaryType> periodicEdges;
std::map<std::pair<DofFace, DofFace>, BoundaryType> periodicFaces;
- // Stores to each DOF all its periodic associations.
- std::map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
-
-
/** \brief
* Defines the interior boundaries of the domain that result from partitioning
* the whole mesh. Contains only the boundaries, which are owned by the rank, i.e.,
diff --git a/AMDiS/src/parallel/ParMetisPartitioner.cc b/AMDiS/src/parallel/ParMetisPartitioner.cc
index b1f88900..78d4d3f5 100644
--- a/AMDiS/src/parallel/ParMetisPartitioner.cc
+++ b/AMDiS/src/parallel/ParMetisPartitioner.cc
@@ -364,6 +364,8 @@ namespace AMDiS {
for (int i = 0; i < nElements; i++) {
wgts[i] = min(wgts[i], kpartMax);
+ wgts[i] = max(wgts[i], 1);
+
smin = std::min(smin, wgts[i]);
smax = std::max(smax, wgts[i]);
ssum += wgts[i];
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index f4d65426..8fbf82ab 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/AMDiS/src/parallel/ParallelDebug.cc
@@ -67,8 +67,11 @@ namespace AMDiS {
// === To the last, do the same of periodic boundaries. ===
- for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
+ for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
+ if (rankIt->first == pdb.mpiRank)
+ continue;
+
int nValues = rankIt->second.size();
int* sBuffer = new int[nValues];
for (int i = 0; i < nValues; i++)
@@ -118,11 +121,17 @@ namespace AMDiS {
for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
+ if (rankIt->first == pdb.mpiRank)
+ continue;
+
for (unsigned int i = 0; i < rankIt->second.size(); i++) {
int elIndex1 = recvBuffers[bufCounter][i];
int elIndex2 = pdb.periodicBoundary.boundary[rankIt->first][i].neighObj.elIndex;
- TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at periodic boundary!\n");
+ TEST_EXIT(elIndex1 == elIndex2)
+ ("Wrong element index at periodic boundary el %d with rank %d: %d %d\n",
+ pdb.periodicBoundary.boundary[rankIt->first][i].rankObj.elIndex,
+ rankIt->first, elIndex1, elIndex2);
}
delete [] recvBuffers[bufCounter++];
diff --git a/AMDiS/src/parallel/ParallelProblemStatBase.cc b/AMDiS/src/parallel/ParallelProblemStatBase.cc
index 015d1727..d7d95dc2 100644
--- a/AMDiS/src/parallel/ParallelProblemStatBase.cc
+++ b/AMDiS/src/parallel/ParallelProblemStatBase.cc
@@ -26,14 +26,16 @@ namespace AMDiS {
ProblemVec::buildAfterCoarsen(adaptInfo, flag, assembleMatrix, assembleVector);
double vm, rss;
- processMemUsage(vm, rss);
+ processMemUsage(vm, rss);
+ vm /= 1024.0;
+ rss /= 1024.0;
- MSG("My memory usage is VM = %f RSS = %f\n", vm, rss);
+ MSG("My memory usage is VM = %f MB RSS = %f MB\n", vm, rss);
mpi::globalAdd(vm);
mpi::globalAdd(rss);
-
- MSG("Overall memory usage is VM = %f RSS = %f\n", vm, rss);
+
+ MSG("Overall memory usage is VM = %f MB RSS = %f MB\n", vm, rss);
}
diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc
index e8e25b1e..202b77aa 100644
--- a/AMDiS/src/parallel/PetscSolver.cc
+++ b/AMDiS/src/parallel/PetscSolver.cc
@@ -128,8 +128,8 @@ namespace AMDiS {
}
- // === Up to now we have assembled on row. Now, the row must be send to the ===
- // === corresponding rows to the petsc matrix. ===
+ // === Up to now we have assembled one row. Now, the row must be send to the ===
+ // === corresponding rows in the petsc matrix. ===
if (periodicRow) {
// The row dof is periodic, so send dof to all the corresponding rows.
@@ -144,7 +144,8 @@ namespace AMDiS {
MatSetValues(petscMatrix, 1, &rowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
- for (std::set<int>::iterator perIt = perAsc.begin(); perIt != perAsc.end(); ++perIt) {
+ for (std::set<int>::iterator perIt = perAsc.begin();
+ perIt != perAsc.end(); ++perIt) {
std::vector<int> perCols;
perCols.reserve(300);
std::vector<double> perValues;
diff --git a/AMDiS/src/parallel/StdMpi.cc b/AMDiS/src/parallel/StdMpi.cc
index a95c2d5c..aab70545 100644
--- a/AMDiS/src/parallel/StdMpi.cc
+++ b/AMDiS/src/parallel/StdMpi.cc
@@ -60,7 +60,7 @@ namespace AMDiS {
int intSizeOf(std::vector<AtomicBoundary> &data)
{
- return data.size() * 3;
+ return data.size() * 6;
}
int intSizeOf(std::vector<BoundaryObject> &data)
@@ -242,9 +242,12 @@ namespace AMDiS {
void makeBuf(std::vector<AtomicBoundary> &data, int *buf)
{
for (unsigned int i = 0; i < data.size(); i++) {
- buf[i * 3] = data[i].rankObj.elIndex;
- buf[i * 3 + 1] = data[i].rankObj.subObj;
- buf[i * 3 + 2] = data[i].rankObj.ithObj;
+ buf[i * 6] = data[i].rankObj.elIndex;
+ buf[i * 6 + 1] = data[i].rankObj.subObj;
+ buf[i * 6 + 2] = data[i].rankObj.ithObj;
+ buf[i * 6 + 3] = data[i].neighObj.elIndex;
+ buf[i * 6 + 4] = data[i].neighObj.subObj;
+ buf[i * 6 + 5] = data[i].neighObj.ithObj;
}
}
@@ -253,13 +256,16 @@ namespace AMDiS {
if (bufSize == 0)
return;
- TEST_EXIT(bufSize % 3 == 0)("This should not happen!\n");
+ TEST_EXIT(bufSize % 6 == 0)("This should not happen!\n");
- data.resize(bufSize / 3);
- for (int i = 0; i < bufSize / 3; i++) {
- data[i].rankObj.elIndex = buf[i * 3];
- data[i].rankObj.subObj = static_cast<GeoIndex>(buf[i * 3 + 1]);
- data[i].rankObj.ithObj = buf[i * 3 + 2];
+ data.resize(bufSize / 6);
+ for (int i = 0; i < bufSize / 6; i++) {
+ data[i].rankObj.elIndex = buf[i * 6];
+ data[i].rankObj.subObj = static_cast<GeoIndex>(buf[i * 6 + 1]);
+ data[i].rankObj.ithObj = buf[i * 6 + 2];
+ data[i].neighObj.elIndex = buf[i * 6 + 3];
+ data[i].neighObj.subObj = static_cast<GeoIndex>(buf[i * 6 + 4]);
+ data[i].neighObj.ithObj = buf[i * 6 + 5];
}
}
--
GitLab