diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index c642788e9c5a11e8723b9074361783291466096a..994e8ff88dda0ca5ec32dd1d7ace1e2b7f489e3f 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/AMDiS/src/AMDiS_fwd.h
@@ -99,6 +99,7 @@ namespace AMDiS {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
class FeSpaceDofMap;
+ class MeshLevelData;
#endif
struct BoundaryObject;
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 53a6d21e9be023a4792117870ca3e9bf4d648e01..6b32912b01c5102ae5f362edb1a61365d18856cc 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -238,7 +238,7 @@ namespace AMDiS {
if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- if (dofMap->isRankDof(rowIndices[i]))
+ if (dofMap->isRankDof(rowIndices[i], 0))
applyDBCs.insert(static_cast<int>(row));
#else
applyDBCs.insert(static_cast<int>(row));
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 94813642e05f4b2000987ac204f8f66795960638..8625f12cdfd6cd08759432e2f71778f4cdb2e474 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -248,7 +248,7 @@ namespace AMDiS {
{
TEST_EXIT_DBG(dofMap)("No rank dofs set!\n");
- return dofMap->isRankDof(dof);
+ return dofMap->isRankDof(dof, 0);
}
#endif
diff --git a/AMDiS/src/io/Spreadsheet.cc b/AMDiS/src/io/Spreadsheet.cc
index 5cb39df11ceebf3a60061c0b99107a179b80c1f7..96ef2bb3028d052623ce361619c21eb918b65282 100644
--- a/AMDiS/src/io/Spreadsheet.cc
+++ b/AMDiS/src/io/Spreadsheet.cc
@@ -54,6 +54,8 @@ namespace AMDiS {
void Spreadsheet::read(string filename)
{
+ FUNCNAME("Spreadsheet::read()");
+
data.clear();
string line;
diff --git a/AMDiS/src/parallel/DofComm.cc b/AMDiS/src/parallel/DofComm.cc
index aa973a91a4f6d94e35a8ee16e9f4063930a2d99d..28fa6378293f93eec4a0e53385afdaebd5b40713 100644
--- a/AMDiS/src/parallel/DofComm.cc
+++ b/AMDiS/src/parallel/DofComm.cc
@@ -41,7 +41,26 @@ namespace AMDiS {
}
}
-
+
+ int DofComm::getNumberDofs(int level, const FiniteElemSpace *feSpace)
+ {
+ FUNCNAME("DofComm::getNumberDofs()");
+
+ TEST_EXIT_DBG(level < data.size())("Should not happen!\n");
+
+ DofContainerSet dofs;
+
+ for (DataIter rankIt = data[level].begin();
+ rankIt != data[level].end(); ++rankIt)
+ for (FeMapIter feIt = rankIt->second.begin();
+ feIt != rankIt->second.end(); ++feIt)
+ if (feIt->first == feSpace)
+ dofs.insert(feIt->second.begin(), feIt->second.end());
+
+ return static_cast<int>(dofs.size());
+ }
+
+
bool DofComm::Iterator::setNextFeMap()
{
FUNCNAME("DofComm::Iterator::setNextFeMap()");
diff --git a/AMDiS/src/parallel/DofComm.h b/AMDiS/src/parallel/DofComm.h
index dd33200176e5a1df2d81ffe0232eb6e2e82958c3..6e5467a12434fb778f0dc32ed1bec3d0ccde6b39 100644
--- a/AMDiS/src/parallel/DofComm.h
+++ b/AMDiS/src/parallel/DofComm.h
@@ -65,6 +65,8 @@ namespace AMDiS {
return data[level];
}
+ int getNumberDofs(int level, const FiniteElemSpace *feSpace);
+
protected:
LevelDataType data;
diff --git a/AMDiS/src/parallel/ElementObjectDatabase.cc b/AMDiS/src/parallel/ElementObjectDatabase.cc
index 93308fb6937cc6621173c1f5f5507418076748d4..f169e8fe7a7ad7f95c434af5375c1968ae47fe00 100644
--- a/AMDiS/src/parallel/ElementObjectDatabase.cc
+++ b/AMDiS/src/parallel/ElementObjectDatabase.cc
@@ -376,9 +376,13 @@ namespace AMDiS {
}
- void ElementObjectDatabase::createRankData(map<int, int>& macroElementRankMap)
+ void ElementObjectDatabase::createRankData(map<int, int>& macroElementRankMap,
+ MeshLevelData& levelData)
{
FUNCNAME("ElementObjectDatabase::createRankData()");
+
+ int nLevel = levelData.getLevelNumber();
+ TEST_EXIT_DBG(nLevel > 0)("Should not happen!\n");
vertexOwner.clear();
vertexInRank.clear();
@@ -391,7 +395,11 @@ namespace AMDiS {
if (it2->elIndex > vertexInRank[it->first][elementInRank].elIndex)
vertexInRank[it->first][elementInRank] = *it2;
- vertexOwner[it->first] = std::max(vertexOwner[it->first], elementInRank);
+ for (int level = 0; level < nLevel; level++) {
+ int levelId = levelData.getLevelId(level, elementInRank);
+ vertexOwner[it->first][level] =
+ std::max(vertexOwner[it->first][level], levelId);
+ }
}
}
@@ -407,7 +415,11 @@ namespace AMDiS {
if (it2->elIndex > edgeInRank[it->first][elementInRank].elIndex)
edgeInRank[it->first][elementInRank] = *it2;
- edgeOwner[it->first] = std::max(edgeOwner[it->first], elementInRank);
+ for (int level = 0; level < nLevel; level++) {
+ int levelId = levelData.getLevelId(level, elementInRank);
+ edgeOwner[it->first][level] =
+ std::max(edgeOwner[it->first][level], levelId);
+ }
}
}
@@ -423,7 +435,11 @@ namespace AMDiS {
if (it2->elIndex > faceInRank[it->first][elementInRank].elIndex)
faceInRank[it->first][elementInRank] = *it2;
- faceOwner[it->first] = std::max(faceOwner[it->first], elementInRank);
+ for (int level = 0; level < nLevel; level++) {
+ int levelId = levelData.getLevelId(level, elementInRank);
+ faceOwner[it->first][level] =
+ std::max(faceOwner[it->first][level], levelId);
+ }
}
}
}
@@ -597,11 +613,12 @@ namespace AMDiS {
}
-
+ ERROR_EXIT("REWRITE SERIALIZATION!\n");
+ /*
SerUtil::serialize(out, vertexOwner);
SerUtil::serialize(out, edgeOwner);
SerUtil::serialize(out, faceOwner);
-
+ */
nSize = vertexInRank.size();
SerUtil::serialize(out, nSize);
@@ -737,6 +754,7 @@ namespace AMDiS {
}
+ ERROR_EXIT("REWRITE DESERIALIZATION!\n");
SerUtil::deserialize(in, vertexOwner);
SerUtil::deserialize(in, edgeOwner);
SerUtil::deserialize(in, faceOwner);
diff --git a/AMDiS/src/parallel/ElementObjectDatabase.h b/AMDiS/src/parallel/ElementObjectDatabase.h
index b68005ad799eac238b51ac7c7150dd239e7297d1..a7ecb5989c2a9c6d7f39b2955cde10dcd92dace7 100644
--- a/AMDiS/src/parallel/ElementObjectDatabase.h
+++ b/AMDiS/src/parallel/ElementObjectDatabase.h
@@ -28,6 +28,7 @@
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_comparison.hpp>
+#include "AMDiS_fwd.h"
#include "Containers.h"
#include "Global.h"
#include "Boundary.h"
@@ -144,7 +145,8 @@ namespace AMDiS {
* \param[in] macroElementRankMap Maps to each macro element of the mesh
* the rank that owns this macro element.
*/
- void createRankData(map<int, int>& macroElementRankMap);
+ void createRankData(map<int, int>& macroElementRankMap,
+ MeshLevelData& levelData);
/** \brief
@@ -257,17 +259,17 @@ namespace AMDiS {
/// Returns the rank owner of the current iterator position.
- int getIterateOwner()
+ int getIterateOwner(int level)
{
switch (iterGeoPos) {
case VERTEX:
- return vertexOwner[vertexIter->first];
+ return vertexOwner[vertexIter->first][level];
break;
case EDGE:
- return edgeOwner[edgeIter->first];
+ return edgeOwner[edgeIter->first][level];
break;
case FACE:
- return faceOwner[faceIter->first];
+ return faceOwner[faceIter->first][level];
break;
default:
ERROR_EXIT("Should not happen!\n");
@@ -279,21 +281,21 @@ namespace AMDiS {
/// Returns the rank owner of a vertex DOF.
- int getOwner(DegreeOfFreedom vertex)
+ int getOwner(DegreeOfFreedom vertex, int level)
{
- return vertexOwner[vertex];
+ return vertexOwner[vertex][level];
}
/// Returns the rank owner of an edge.
- int getOwner(DofEdge edge)
+ int getOwner(DofEdge edge, int level)
{
- return edgeOwner[edge];
+ return edgeOwner[edge][level];
}
/// Returns the rank owner of an face.
- int getOwner(DofFace face)
+ int getOwner(DofFace face, int level)
{
- return faceOwner[face];
+ return faceOwner[face][level];
}
@@ -540,15 +542,17 @@ namespace AMDiS {
/// Maps to an element object the corresponding face.
map<ElementObjectData, DofFace> faceLocalMap;
+ /// Maps from level to rank number
+ typedef map<int, int> LevelRank;
/// Defines for all vertex DOFs the rank that ownes this vertex DOF.
- map<DegreeOfFreedom, int> vertexOwner;
+ map<DegreeOfFreedom, LevelRank> vertexOwner;
/// Defines for all edges the rank that ownes this edge.
- map<DofEdge, int> edgeOwner;
+ map<DofEdge, LevelRank> edgeOwner;
/// Defines for all faces the rank that ownes this face.
- map<DofFace, int> faceOwner;
+ map<DofFace, LevelRank> faceOwner;
/// Defines to each vertex DOF a map that maps to each rank number the element
diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc
index f3e694a4435b195970036e42300b825aed2f506b..62b13aa6faf0ab12a77a0b5d27d7a214cbfda4fa 100644
--- a/AMDiS/src/parallel/InteriorBoundary.cc
+++ b/AMDiS/src/parallel/InteriorBoundary.cc
@@ -18,10 +18,11 @@
namespace AMDiS {
- AtomicBoundary& InteriorBoundary::getNewAtomic(int rank)
+ AtomicBoundary& InteriorBoundary::getNewAtomic(int level, int rank)
{
- boundary[rank].resize(boundary[rank].size() + 1);
- return boundary[rank][boundary[rank].size() - 1];
+ int size = boundary[level][rank].size();
+ boundary[level][rank].resize(size + 1);
+ return boundary[level][rank][size];
}
@@ -29,32 +30,51 @@ namespace AMDiS {
{
InteriorBoundary& other2 = const_cast<InteriorBoundary&>(other);
- for (RankToBoundMap::const_iterator it = boundary.begin();
- it != boundary.end(); ++it) {
- if (other2.boundary.count(it->first) == 0)
- return false;
-
- if (other2.boundary[it->first].size() != it->second.size())
- return false;
-
- for (unsigned int i = 0; i < it->second.size(); i++) {
- std::vector<AtomicBoundary>::iterator bIt =
- find(other2.boundary[it->first].begin(),
- other2.boundary[it->first].end(), it->second[i]);
+ if (boundary.size() != other2.boundary.size())
+ return false;
+
+ for (unsigned int level = 0; level < boundary.size(); level++) {
+ for (RankToBoundMap::const_iterator it = boundary[level].begin();
+ it != boundary[level].end(); ++it) {
+ if (other2.boundary[level].count(it->first) == 0)
+ return false;
+
+ if (other2.boundary[level][it->first].size() != it->second.size())
+ return false;
- if (bIt == other2.boundary[it->first].end())
- return false;
- }
+ for (unsigned int i = 0; i < it->second.size(); i++) {
+ std::vector<AtomicBoundary>::iterator bIt =
+ find(other2.boundary[level][it->first].begin(),
+ other2.boundary[level][it->first].end(), it->second[i]);
+
+ if (bIt == other2.boundary[level][it->first].end())
+ return false;
+ }
+ }
}
return true;
}
+ void InteriorBoundary::reset(int level)
+ {
+ FUNCNAME("InteriorBoundary::reset()");
+
+ nLevel = level;
+
+ boundary.clear();
+ boundary.resize(nLevel);
+ }
+
+
void InteriorBoundary::serialize(std::ostream &out)
{
FUNCNAME("InteriorBoundary::serialize()");
+ ERROR_EXIT("REWRITE TO MULTILEVEL STRUCTURE!\n");
+
+#if 0
int mSize = boundary.size();
SerUtil::serialize(out, mSize);
for (RankToBoundMap::iterator it = boundary.begin();
@@ -83,6 +103,7 @@ namespace AMDiS {
SerUtil::serialize(out, bound.type);
}
}
+#endif
}
@@ -91,6 +112,9 @@ namespace AMDiS {
{
FUNCNAME("InteriorBoundary::deserialize()");
+ ERROR_EXIT("REWRITE TO MULTILEVEL STRUCTURE!\n");
+
+#if 0
int mSize = 0;
SerUtil::deserialize(in, mSize);
for (int i = 0; i < mSize; i++) {
@@ -137,6 +161,7 @@ namespace AMDiS {
bound.neighObj.el = NULL;
}
}
+#endif
}
diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h
index 3102165e7eb4828ac6d3595b322e8ff73a62b86d..bdc02b773648b55be2e199ff5ac2d5075512752a 100644
--- a/AMDiS/src/parallel/InteriorBoundary.h
+++ b/AMDiS/src/parallel/InteriorBoundary.h
@@ -47,34 +47,35 @@ namespace AMDiS {
public:
iterator(InteriorBoundary &b)
: bound(b),
- levelData(NULL),
level(0)
{
reset();
}
- iterator(InteriorBoundary &b, MeshLevelData &levelData, int level)
+ iterator(InteriorBoundary &b, int level)
: bound(b),
- levelData(&levelData),
level(level)
{
+ TEST_EXIT_DBG(level < bound.boundary.size())
+ ("Should not happen!\n");
+
reset();
}
/// Set the iterator to the first position.
void reset()
{
- mapIt = bound.boundary.begin();
+ mapIt = bound.boundary[level].begin();
nextNonempty();
- if (mapIt != bound.boundary.end())
+ if (mapIt != bound.boundary[level].end())
vecIt = mapIt->second.begin();
}
/// Test if iterator is at the final position.
bool end() const
{
- return (mapIt == bound.boundary.end());
+ return (mapIt == bound.boundary[level].end());
}
/// Move iterator to the next position.
@@ -85,7 +86,7 @@ namespace AMDiS {
++mapIt;
nextNonempty();
- if (mapIt != bound.boundary.end())
+ if (mapIt != bound.boundary[level].end())
vecIt = mapIt->second.begin();
}
}
@@ -105,22 +106,12 @@ namespace AMDiS {
++mapIt;
nextNonempty();
- if (mapIt != bound.boundary.end())
+ if (mapIt != bound.boundary[level].end())
vecIt = mapIt->second.begin();
}
inline int getRank()
{
- if (level > 0) {
- int r = levelData->mapRank(mapIt->first, level - 1, level);
-
- TEST_EXIT_DBG(r >= 0)
- ("Mapping rank %d from level % to level %d does not work!\n",
- mapIt->first, level - 1, level);
-
- return r;
- }
-
return mapIt->first;
}
@@ -128,25 +119,13 @@ namespace AMDiS {
inline void nextNonempty()
{
- if (mapIt == bound.boundary.end())
+ if (mapIt == bound.boundary[level].end())
return;
- if (level > 0) {
- TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");
- TEST_EXIT_DBG(level == 1)("Only 2-level method supported!\n");
-
- while (!levelData->rankInSubdomain(mapIt->first, level) ||
- mapIt->second.size() == 0) {
- ++mapIt;
- if (mapIt == bound.boundary.end())
- return;
- }
- } else {
- while (mapIt->second.size() == 0) {
- ++mapIt;
- if (mapIt == bound.boundary.end())
- return;
- }
+ while (mapIt->second.size() == 0) {
+ ++mapIt;
+ if (mapIt == bound.boundary[level].end())
+ return;
}
}
@@ -157,20 +136,19 @@ namespace AMDiS {
InteriorBoundary &bound;
- MeshLevelData *levelData;
-
int level;
};
public:
- InteriorBoundary() {}
-
- void clear()
+ InteriorBoundary(int l = 1)
+ : nLevel(l)
{
- boundary.clear();
+ boundary.resize(nLevel);
}
- AtomicBoundary& getNewAtomic(int rank);
+ void reset(int nLevel);
+
+ AtomicBoundary& getNewAtomic(int level, int rank);
/// Writes this object to a file.
void serialize(ostream &out);
@@ -185,10 +163,13 @@ namespace AMDiS {
protected:
void serializeExcludeList(ostream &out, ExcludeList &list);
- void deserializeExcludeList(istream &in, ExcludeList &list);
+ void deserializeExcludeList(istream &in, ExcludeList &list);
public:
- RankToBoundMap boundary;
+ vector<RankToBoundMap> boundary;
+
+ protected:
+ int nLevel;
};
}
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 027993db02295f18f69375e420708415f6c6439f..ee86021131c18282f24e499f0bc0f8425b3eae86 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -232,9 +232,10 @@ namespace AMDiS {
}
}
+ // If required, create hierarchical mesh level structure.
+ createMeshLevelStructure();
- // === Create interior boundary information. ===
-
+ // Create interior boundary information.
createInteriorBoundaryInfo();
#if (DEBUG != 0)
@@ -294,11 +295,6 @@ namespace AMDiS {
it != mesh->getPeriodicAssociations().end(); ++it)
const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
-
- // If required, create hierarchical mesh level structure.
- createMeshLevelStructure();
-
-
updateLocalGlobalNumbering();
// === In 3D we have to make some test, if the resulting mesh is valid. ===
@@ -784,18 +780,73 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::createMeshLevelStructure()");
+ if (mpiSize != 16)
+ return;
+
std::set<int> neighbours;
- for (InteriorBoundary::iterator it(rankIntBoundary); !it.end(); ++it)
- neighbours.insert(it.getRank());
+ switch (mpiRank) {
+ case 0:
+ neighbours.insert(1); neighbours.insert(4); neighbours.insert(5);
+ break;
+ case 1:
+ neighbours.insert(0); neighbours.insert(2); neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);
+ break;
+ case 2:
+ neighbours.insert(1); neighbours.insert(3); neighbours.insert(5); neighbours.insert(6); neighbours.insert(7);
+ break;
+ case 3:
+ neighbours.insert(2); neighbours.insert(6); neighbours.insert(7);
+ break;
+ case 4:
+ neighbours.insert(0); neighbours.insert(1); neighbours.insert(5); neighbours.insert(8); neighbours.insert(9);
+ break;
+ case 5:
+ neighbours.insert(0); neighbours.insert(1); neighbours.insert(2);
+ neighbours.insert(4); neighbours.insert(6);
+ neighbours.insert(8); neighbours.insert(9); neighbours.insert(10);
+ break;
+ case 6:
+ neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
+ neighbours.insert(5); neighbours.insert(7);
+ neighbours.insert(9); neighbours.insert(10); neighbours.insert(11);
+ break;
+ case 7:
+ neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); neighbours.insert(10); neighbours.insert(11);
+ break;
+ case 8:
+ neighbours.insert(4); neighbours.insert(5); neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);
+ break;
+ case 9:
+ neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);
+ neighbours.insert(8); neighbours.insert(10);
+ neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
+ break;
+ case 10:
+ neighbours.insert(5); neighbours.insert(6); neighbours.insert(7);
+ neighbours.insert(9); neighbours.insert(11);
+ neighbours.insert(13); neighbours.insert(14); neighbours.insert(15);
+ break;
+ case 11:
+ neighbours.insert(6); neighbours.insert(7); neighbours.insert(10); neighbours.insert(14); neighbours.insert(15);
+ break;
+ case 12:
+ neighbours.insert(8); neighbours.insert(9); neighbours.insert(13);
+ break;
+ case 13:
+ neighbours.insert(8); neighbours.insert(9); neighbours.insert(10); neighbours.insert(12); neighbours.insert(14);
+ break;
+ case 14:
+ neighbours.insert(9); neighbours.insert(10); neighbours.insert(11); neighbours.insert(13); neighbours.insert(15);
+ break;
+ case 15:
+ neighbours.insert(10); neighbours.insert(11); neighbours.insert(14);
+ break;
+ }
- for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
- neighbours.insert(it.getRank());
+ TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
levelData.init(neighbours);
- if (mpiSize != 16)
- return;
-
bool multiLevelTest = false;
Parameters::get("parallel->multi level test", multiLevelTest);
if (multiLevelTest) {
@@ -1472,7 +1523,7 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::updateInteriorBoundaryInfo()");
- elObjDb.createRankData(partitionMap);
+ elObjDb.createRankData(partitionMap, levelData);
createBoundaryData();
#if (DEBUG != 0)
@@ -1515,7 +1566,7 @@ namespace AMDiS {
macroElIndexTypeMap);
// Create mesh element data for this rank.
- elObjDb.createRankData(partitionMap);
+ elObjDb.createRankData(partitionMap, levelData);
}
@@ -1525,10 +1576,19 @@ namespace AMDiS {
// === Clear all relevant data structures. ===
- rankIntBoundary.clear();
- otherIntBoundary.clear();
- periodicBoundary.clear();
+ int nLevel = levelData.getLevelNumber();
+ rankIntBoundary.reset(nLevel);
+ otherIntBoundary.reset(nLevel);
+ periodicBoundary.reset(nLevel);
+ for (int level = 0; level < nLevel; level++)
+ createBoundaryData(level);
+ }
+
+
+ void MeshDistributor::createBoundaryData(int level)
+ {
+ FUNCNAME("MeshDistributor::createBoundaryData()");
// === Create interior boundary data structure. ===
@@ -1540,7 +1600,7 @@ namespace AMDiS {
if (!(objData.count(mpiRank) && objData.size() > 1))
continue;
- int owner = elObjDb.getIterateOwner();
+ int owner = elObjDb.getIterateOwner(level);
ElementObjectData& rankBoundEl = objData[mpiRank];
AtomicBoundary bound;
@@ -1574,7 +1634,7 @@ namespace AMDiS {
bound.type = INTERIOR;
- AtomicBoundary& b = rankIntBoundary.getNewAtomic(it2->first);
+ AtomicBoundary& b = rankIntBoundary.getNewAtomic(level, it2->first);
b = bound;
if (geoIndex == EDGE)
b.neighObj.reverseMode =
@@ -1598,7 +1658,7 @@ namespace AMDiS {
bound.type = INTERIOR;
- AtomicBoundary& b = otherIntBoundary.getNewAtomic(owner);
+ AtomicBoundary& b = otherIntBoundary.getNewAtomic(level, owner);
b = bound;
if (geoIndex == EDGE)
b.rankObj.reverseMode =
@@ -1618,16 +1678,12 @@ namespace AMDiS {
if (elObjDb.isInRank(it->first.first, mpiRank) == false)
continue;
- // MSG("PERIODIC DOF: %d -> %d with ASOC %d\n", it->first.first, it->first.second, it->second);
-
ElementObjectData& perDofEl0 =
elObjDb.getElementsInRank(it->first.first)[mpiRank];
for (map<int, ElementObjectData>::iterator elIt = elObjDb.getElementsInRank(it->first.second).begin();
elIt != elObjDb.getElementsInRank(it->first.second).end(); ++elIt) {
- // MSG(" in el %d %d\n", perDofEl0.elIndex, perDofEl0.ithObject);
-
int otherElementRank = elIt->first;
ElementObjectData& perDofEl1 = elIt->second;
@@ -1651,11 +1707,11 @@ namespace AMDiS {
bound.rankObj.ithObj == 1 ||
bound.rankObj.elIndex == 78 &&
bound.rankObj.ithObj == 2) {
- AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
+ AtomicBoundary& b = periodicBoundary.getNewAtomic(level, otherElementRank);
b = bound;
}
} else {
- AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
+ AtomicBoundary& b = periodicBoundary.getNewAtomic(level, otherElementRank);
b = bound;
}
}
@@ -1690,7 +1746,7 @@ namespace AMDiS {
bound.type = it->second;
- AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
+ AtomicBoundary& b = periodicBoundary.getNewAtomic(level, otherElementRank);
b = bound;
if (mpiRank > otherElementRank)
@@ -1736,7 +1792,7 @@ namespace AMDiS {
bound.type = it->second;
- AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
+ AtomicBoundary& b = periodicBoundary.getNewAtomic(level, otherElementRank);
b = bound;
if (mpiRank > otherElementRank)
@@ -1755,8 +1811,8 @@ namespace AMDiS {
// === share the bounday. ===
StdMpi<vector<AtomicBoundary> > stdMpi(mpiComm);
- stdMpi.send(rankIntBoundary.boundary);
- stdMpi.recv(otherIntBoundary.boundary);
+ stdMpi.send(rankIntBoundary.boundary[level]);
+ stdMpi.recv(otherIntBoundary.boundary[level]);
stdMpi.startCommunication();
@@ -1765,8 +1821,8 @@ namespace AMDiS {
// === the same order. If not, the atomic boundaries are swaped to the ===
// === correct order. ===
- for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
- rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
+ for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary[level].begin();
+ rankIt != otherIntBoundary.boundary[level].end(); ++rankIt) {
// === We have received from rank "rankIt->first" the ordered list of ===
// === element indices. Now, we have to sort the corresponding list in ===
@@ -1800,28 +1856,29 @@ namespace AMDiS {
// === Do the same for the periodic boundaries. ===
- if (periodicBoundary.boundary.size() > 0) {
+ if (periodicBoundary.boundary[level].size() > 0) {
stdMpi.clear();
- InteriorBoundary sendBounds, recvBounds;
- for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
- rankIt != periodicBoundary.boundary.end(); ++rankIt) {
+ InteriorBoundary sendBounds(levelData.getLevelNumber());
+ InteriorBoundary recvBounds(levelData.getLevelNumber());
+ for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary[level].begin();
+ rankIt != periodicBoundary.boundary[level].end(); ++rankIt) {
if (rankIt->first == mpiRank)
continue;
if (rankIt->first < mpiRank)
- sendBounds.boundary[rankIt->first] = rankIt->second;
+ sendBounds.boundary[level][rankIt->first] = rankIt->second;
else
- recvBounds.boundary[rankIt->first] = rankIt->second;
+ recvBounds.boundary[level][rankIt->first] = rankIt->second;
}
- stdMpi.send(sendBounds.boundary);
- stdMpi.recv(recvBounds.boundary);
+ stdMpi.send(sendBounds.boundary[level]);
+ stdMpi.recv(recvBounds.boundary[level]);
stdMpi.startCommunication();
- for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
- rankIt != periodicBoundary.boundary.end(); ++rankIt) {
+ for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary[level].begin();
+ rankIt != periodicBoundary.boundary[level].end(); ++rankIt) {
if (rankIt->first <= mpiRank)
continue;
@@ -1832,12 +1889,12 @@ namespace AMDiS {
BoundaryObject &recvNeighObj =
stdMpi.getRecvData()[rankIt->first][j].neighObj;
- if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvRankObj ||
- periodicBoundary.boundary[rankIt->first][j].rankObj != recvNeighObj) {
+ if (periodicBoundary.boundary[level][rankIt->first][j].neighObj != recvRankObj ||
+ periodicBoundary.boundary[level][rankIt->first][j].rankObj != recvNeighObj) {
unsigned int k = j + 1;
for (; k < rankIt->second.size(); k++)
- if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvRankObj &&
- periodicBoundary.boundary[rankIt->first][k].rankObj == recvNeighObj)
+ if (periodicBoundary.boundary[level][rankIt->first][k].neighObj == recvRankObj &&
+ periodicBoundary.boundary[level][rankIt->first][k].rankObj == recvNeighObj)
break;
// The element must always be found, because the list is just in
@@ -1851,7 +1908,7 @@ namespace AMDiS {
}
}
}
- } // periodicBoundary.boundary.size() > 0
+ } // periodicBoundary.boundary[level].size() > 0
}
@@ -1888,7 +1945,7 @@ namespace AMDiS {
// === Create send DOFs. ===
for (int geo = FACE; geo >= VERTEX; geo--) {
- for (InteriorBoundary::iterator it(rankIntBoundary, levelData, level);
+ for (InteriorBoundary::iterator it(rankIntBoundary, level);
!it.end(); ++it) {
if (it->rankObj.subObj == geo) {
DofContainer dofs;
@@ -1906,7 +1963,7 @@ namespace AMDiS {
// === Create recv DOFs. ===
for (int geo = FACE; geo >= VERTEX; geo--) {
- for (InteriorBoundary::iterator it(otherIntBoundary, levelData, level);
+ for (InteriorBoundary::iterator it(otherIntBoundary, level);
!it.end(); ++it) {
if (it->rankObj.subObj == geo) {
DofContainer dofs;
@@ -1922,12 +1979,12 @@ namespace AMDiS {
}
}
} else {
- for (InteriorBoundary::iterator it(rankIntBoundary, levelData, level);
+ for (InteriorBoundary::iterator it(rankIntBoundary, level);
!it.end(); ++it)
it->rankObj.el->getAllDofs(feSpace, it->rankObj,
sendDofs.getDofContainer(it.getRank(), feSpace, level));
- for (InteriorBoundary::iterator it(otherIntBoundary, levelData, level);
+ for (InteriorBoundary::iterator it(otherIntBoundary, level);
!it.end(); ++it)
it->rankObj.el->getAllDofs(feSpace, it->rankObj,
recvDofs.getDofContainer(it.getRank(), feSpace, level));
@@ -1966,14 +2023,12 @@ namespace AMDiS {
int nLevels = levelData.getLevelNumber();
TEST_EXIT_DBG(nLevels >= 1)("Should not happen!\n");
- sendDofs.init(nLevels);
- recvDofs.init(nLevels);
- boundaryDofInfo.resize(nLevels);
-
dofMap.init(levelData, feSpaces, feSpaces, true, true);
dofMap.setDofComm(sendDofs, recvDofs);
dofMap.clear();
+ createBoundaryDofs();
+
for (unsigned int i = 0; i < feSpaces.size(); i++)
updateLocalGlobalNumbering(feSpaces[i]);
@@ -2027,22 +2082,10 @@ namespace AMDiS {
int nLevels = levelData.getLevelNumber();
for (int level = 0; level < nLevels; level++) {
- createBoundaryDofs(feSpace, level);
-
- // All DOFs that must be received are DOFs not owned by rank and have
- // therefore to be removed from the set 'rankDofs'.
DofContainerSet nonRankDofs;
- for (DofComm::Iterator it(recvDofs, level, feSpace); !it.end(); it.nextRank()) {
- for (; !it.endDofIter(); it.nextDof()) {
-#if 0
- DofContainer::iterator eraseIt =
- find(rankDofs.begin(), rankDofs.end(), it.getDof());
- if (eraseIt != rankDofs.end())
- rankDofs.erase(eraseIt);
-#endif
+ for (DofComm::Iterator it(recvDofs, level, feSpace); !it.end(); it.nextRank())
+ for (; !it.endDofIter(); it.nextDof())
nonRankDofs.insert(it.getDof());
- }
- }
for (unsigned int i = 0; i < rankDofs.size(); i++)
if (nonRankDofs.count(rankDofs[i]) == 0)
@@ -2069,6 +2112,8 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::createPeriodicMap()");
+ TEST_EXIT(levelData.getLevelNumber() == 1)("Not yet implemented for multilevel!\n");
+
// Clear all periodic DOF mappings calculated before. We do it from scratch.
periodicDofs.init(levelData.getLevelNumber());
periodicMap.clear();
@@ -2077,7 +2122,7 @@ namespace AMDiS {
// periodicMap must be still cleared before: if we do repartitioning and
// there were periodic boundaries in subdomain before and after repartitioning
// there are no more periodic boundaries.
- if (periodicBoundary.boundary.size() == 0)
+ if (periodicBoundary.boundary[0].size() == 0)
return;
for (unsigned int i = 0; i < feSpaces.size(); i++)
@@ -2089,6 +2134,8 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::createPeriodicMap()");
+ TEST_EXIT(levelData.getLevelNumber() == 1)("No yet implemented for multilevel stuff!\n");
+
StdMpi<vector<int> > stdMpi(mpiComm, false);
// === Each rank traverse its periodic boundaries and sends the DOF ===
@@ -2096,8 +2143,8 @@ namespace AMDiS {
map<int, vector<int> > rankToDofType;
- for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
- it != periodicBoundary.boundary.end(); ++it) {
+ for (RankToBoundMap::iterator it = periodicBoundary.boundary[0].begin();
+ it != periodicBoundary.boundary[0].end(); ++it) {
if (it->first == mpiRank) {
// Here we have a periodic boundary within rank's subdomain. So we can
@@ -2122,8 +2169,8 @@ namespace AMDiS {
BoundaryType type = bound.type;
for (unsigned int j = 0; j < dofs0.size(); j++) {
- DegreeOfFreedom globalDof0 = dofMap[feSpace][*(dofs0[j])].global;
- DegreeOfFreedom globalDof1 = dofMap[feSpace][*(dofs1[j])].global;
+ DegreeOfFreedom globalDof0 = dofMap[feSpace][0][*(dofs0[j])].global;
+ DegreeOfFreedom globalDof1 = dofMap[feSpace][0][*(dofs1[j])].global;
if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDof0))
periodicMap.add(feSpace, type, globalDof0, globalDof1);
@@ -2148,7 +2195,7 @@ namespace AMDiS {
// 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(dofMap[feSpace][*(dofs[i])].global);
+ stdMpi.getSendData(it->first).push_back(dofMap[feSpace][0][*(dofs[i])].global);
// Receive from this rank the same number of dofs.
stdMpi.recv(it->first, dofs.size());
@@ -2165,8 +2212,8 @@ namespace AMDiS {
// === DOFs from the other ranks. ===
- for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
- it != periodicBoundary.boundary.end(); ++it) {
+ for (RankToBoundMap::iterator it = periodicBoundary.boundary[0].begin();
+ it != periodicBoundary.boundary[0].end(); ++it) {
DofContainer& dofs = periodicDofs.getDofContainer(it->first, feSpace);
vector<int>& types = rankToDofType[it->first];
@@ -2174,7 +2221,7 @@ namespace AMDiS {
// Added the received DOFs to the mapping.
for (unsigned int i = 0; i < dofs.size(); i++) {
- int globalDofIndex = dofMap[feSpace][*(dofs[i])].global;
+ int globalDofIndex = dofMap[feSpace][0][*(dofs[i])].global;
int mapGlobalDofIndex = stdMpi.getRecvData(it->first)[i];
BoundaryType type = types[i];
@@ -2189,8 +2236,8 @@ namespace AMDiS {
StdMpi<PeriodicDofMap> stdMpi2(mpiComm);
- for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
- it != periodicBoundary.boundary.end(); ++it) {
+ for (RankToBoundMap::iterator it = periodicBoundary.boundary[0].begin();
+ it != periodicBoundary.boundary[0].end(); ++it) {
if (it->first == mpiRank)
continue;
@@ -2207,7 +2254,7 @@ namespace AMDiS {
boundIt->rankObj.el->getAllDofs(feSpace, boundIt->rankObj, dofs);
for (unsigned int i = 0; i < dofs.size(); i++) {
- DegreeOfFreedom globalDof = dofMap[feSpace][*dofs[i]].global;
+ DegreeOfFreedom globalDof = dofMap[feSpace][0][*dofs[i]].global;
std::set<BoundaryType>& assoc =
periodicMap.getAssociations(feSpace, globalDof);
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index f1b0c2bfab692f88af56ba84794fc0c16d0937d3..0f584268b8f3a2bcde9776bd7f4d4ab4ad25d140 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -255,8 +255,8 @@ namespace AMDiS {
createBoundaryDofFlag = flag;
}
- BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace,
- int level = 0)
+ BoundaryDofInfo& getBoundaryDofInfo(const FiniteElemSpace *feSpace,
+ int level)
{
FUNCNAME("MeshDistributor::getBoundaryDofInfo()");
@@ -297,6 +297,8 @@ namespace AMDiS {
void createBoundaryData();
+ void createBoundaryData(int level);
+
void createBoundaryDofs();
void createBoundaryDofs(const FiniteElemSpace *feSpace, int level);
diff --git a/AMDiS/src/parallel/MeshLevelData.h b/AMDiS/src/parallel/MeshLevelData.h
index 9151705db4b4e87a47fb3f1e978c16ee5198c056..0bb461756e455a665ffed8030c1be4b57d4e2fa0 100644
--- a/AMDiS/src/parallel/MeshLevelData.h
+++ b/AMDiS/src/parallel/MeshLevelData.h
@@ -86,6 +86,31 @@ namespace AMDiS {
return mpiGroups[level];
}
+ int getLevelId(int level, int rank)
+ {
+ if (level == 0)
+ return rank;
+
+ TEST_EXIT_DBG(level == 1 && rank <= 15)("Should not happen!\n");
+ TEST_EXIT_DBG(MPI::COMM_WORLD.Get_size() == 16)("Should not happen!\n");
+
+ if (rank == 0 || rank == 1 || rank == 4 || rank == 5)
+ return 0;
+
+ if (rank == 2 || rank == 3 || rank == 6 || rank == 7)
+ return 1;
+
+ if (rank == 8 || rank == 9 || rank == 12 || rank == 13)
+ return 2;
+
+ if (rank == 10 || rank == 11 || rank == 14 || rank == 15)
+ return 3;
+
+ ERROR_EXIT("Should not happen!\n");
+
+ return -1;
+ }
+
int mapRank(int fromRank, int fromLevel, int toLevel)
{
int toRank = -1;
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index d9a63eedadd85a7202940a0c82f13e27ac241553..fbcac563997bff65bdae71cf7a25ca07fcdd2bfb 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/AMDiS/src/parallel/ParallelDebug.cc
@@ -39,8 +39,8 @@ namespace AMDiS {
// === Send rank's boundary information. ===
- for (RankToBoundMap::iterator rankIt = pdb.rankIntBoundary.boundary.begin();
- rankIt != pdb.rankIntBoundary.boundary.end(); ++rankIt) {
+ for (RankToBoundMap::iterator rankIt = pdb.rankIntBoundary.boundary[0].begin();
+ rankIt != pdb.rankIntBoundary.boundary[0].end(); ++rankIt) {
int nSendInt = rankIt->second.size();
int* buffer = new int[nSendInt];
@@ -56,8 +56,8 @@ namespace AMDiS {
// === Receive information from other ranks about the interior boundaries. ====
- for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary.begin();
- rankIt != pdb.otherIntBoundary.boundary.end(); ++rankIt) {
+ for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary[0].begin();
+ rankIt != pdb.otherIntBoundary.boundary[0].end(); ++rankIt) {
int nRecvInt = rankIt->second.size();
int *buffer = new int[nRecvInt];
recvBuffers.push_back(buffer);
@@ -69,8 +69,8 @@ namespace AMDiS {
// === To the last, do the same of periodic boundaries. ===
- for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
- rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
+ for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary[0].begin();
+ rankIt != pdb.periodicBoundary.boundary[0].end(); ++rankIt) {
if (rankIt->first == pdb.mpiRank)
continue;
@@ -103,8 +103,8 @@ namespace AMDiS {
// === and after this the periodic ones. ===
int bufCounter = 0;
- for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary.begin();
- rankIt != pdb.otherIntBoundary.boundary.end(); ++rankIt) {
+ for (RankToBoundMap::iterator rankIt = pdb.otherIntBoundary.boundary[0].begin();
+ rankIt != pdb.otherIntBoundary.boundary[0].end(); ++rankIt) {
TEST_EXIT(rankIt->second.size() ==
pdb.otherIntBoundary.boundary[rankIt->first].size())
@@ -112,7 +112,7 @@ namespace AMDiS {
for (unsigned int i = 0; i < rankIt->second.size(); i++) {
int elIndex1 = recvBuffers[bufCounter][i];
- int elIndex2 = pdb.otherIntBoundary.boundary[rankIt->first][i].neighObj.elIndex;
+ int elIndex2 = pdb.otherIntBoundary.boundary[0][rankIt->first][i].neighObj.elIndex;
TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n");
}
@@ -121,18 +121,18 @@ namespace AMDiS {
}
- for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary.begin();
- rankIt != pdb.periodicBoundary.boundary.end(); ++rankIt) {
+ for (RankToBoundMap::iterator rankIt = pdb.periodicBoundary.boundary[0].begin();
+ rankIt != pdb.periodicBoundary.boundary[0].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;
+ int elIndex2 = pdb.periodicBoundary.boundary[0][rankIt->first][i].neighObj.elIndex;
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,
+ pdb.periodicBoundary.boundary[0][rankIt->first][i].rankObj.elIndex,
rankIt->first, elIndex1, elIndex2);
}
@@ -256,8 +256,8 @@ namespace AMDiS {
RankToCoords sendCoords;
map<int, vector<BoundaryType> > rankToDofType;
- for (InteriorBoundary::RankToBoundMap::iterator it = pdb.periodicBoundary.boundary.begin();
- it != pdb.periodicBoundary.boundary.end(); ++it) {
+ for (InteriorBoundary::RankToBoundMap::iterator it = pdb.periodicBoundary.boundary[0].begin();
+ it != pdb.periodicBoundary.boundary[0].end(); ++it) {
if (it->first == pdb.mpiRank)
continue;
@@ -479,7 +479,7 @@ namespace AMDiS {
DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
for (it.reset(); !it.end(); ++it) {
- coordsToIndex[(*it)] = pdb.dofMap[feSpace][it.getDOFIndex()].global;
+ coordsToIndex[(*it)] = pdb.dofMap[feSpace][0][it.getDOFIndex()].global;
// MSG(" CHECK FOR DOF %d AT COORDS %f %f %f\n",
// coordsToIndex[(*it)], (*it)[0], (*it)[1], (*it)[2]);
}
@@ -780,8 +780,8 @@ namespace AMDiS {
DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
for (it.reset(); !it.end(); ++it) {
file << it.getDOFIndex() << " "
- << pdb.dofMap[feSpace][it.getDOFIndex()].global << " "
- << pdb.dofMap[feSpace].isRankDof(it.getDOFIndex());
+ << pdb.dofMap[feSpace][0][it.getDOFIndex()].global << " "
+ << pdb.dofMap[feSpace].isRankDof(it.getDOFIndex(), 0);
for (int i = 0; i < pdb.mesh->getDim(); i++)
file << " " << (*it)[i];
file << "\n";
diff --git a/AMDiS/src/parallel/ParallelDofMapping.cc b/AMDiS/src/parallel/ParallelDofMapping.cc
index 927847cd2c3746e32b0abc29bf915b6f959195b6..d981fee3782399be4c94a02d898b9e9282e7c97e 100644
--- a/AMDiS/src/parallel/ParallelDofMapping.cc
+++ b/AMDiS/src/parallel/ParallelDofMapping.cc
@@ -86,7 +86,8 @@ namespace AMDiS {
// === Send all global indices of DOFs that are owned by the rank to all ===
// === other ranks that also include this DOF. ===
- StdMpi<vector<int> > stdMpi(levelData->getMpiComm(level));
+ // StdMpi<vector<int> > stdMpi(levelData->getMpiComm(level));
+ StdMpi<vector<int> > stdMpi(levelData->getMpiComm(0));
for (DofComm::Iterator it(*sendDofs, level, feSpace);
!it.end(); it.nextRank())
@@ -171,15 +172,15 @@ namespace AMDiS {
nLocalDofs.resize(nLevel);
nOverallDofs.resize(nLevel);
rStartDofs.resize(nLevel);
+ dofToMatIndex.resize(nLevel);
for (int i = 0; i < nLevel; i++) {
nRankDofs[i] = -1;
nLocalDofs[i] = -1;
nOverallDofs[i] = -1;
rStartDofs[i] = -1;
+ dofToMatIndex[i].clear();
}
-
- dofToMatIndex.clear();
}
@@ -310,9 +311,7 @@ namespace AMDiS {
{
FUNCNAME("ParallelDofMapping::computeMatIndex()");
- TEST_EXIT(level == 0)("Not yet implemented for non-zero level!\n");
-
- dofToMatIndex.clear();
+ dofToMatIndex[level].clear();
// The offset is always added to the local matrix index. The offset for the
// DOFs in the first FE spaces is the smalled global index of a DOF that is
@@ -328,12 +327,12 @@ namespace AMDiS {
// a matrix index.
DofMap& dofMap = data[feSpaces[i]].getMap(level);
for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
- if (data[feSpaces[i]].isRankDof(it->first)) {
+ if (data[feSpaces[i]].isRankDof(it->first, level)) {
int globalMatIndex = it->second.local + offset;
if (globalIndex)
- dofToMatIndex.add(i, it->second.global, globalMatIndex);
+ dofToMatIndex[level].add(i, it->second.global, globalMatIndex);
else
- dofToMatIndex.add(i, it->first, globalMatIndex);
+ dofToMatIndex[level].add(i, it->first, globalMatIndex);
}
}
@@ -351,7 +350,8 @@ namespace AMDiS {
// === Communicate the matrix indices for all DOFs that are on some ===
// === interior boundaries. ===
- StdMpi<vector<DegreeOfFreedom> > stdMpi(levelData->getMpiComm(level));
+ // StdMpi<vector<DegreeOfFreedom> > stdMpi(levelData->getMpiComm(level));
+ StdMpi<vector<DegreeOfFreedom> > stdMpi(levelData->getMpiComm(0));
for (DofComm::Iterator it(*sendDofs, level, feSpaces[i]);
!it.end(); it.nextRank()) {
vector<DegreeOfFreedom> sendGlobalDofs;
@@ -359,9 +359,9 @@ namespace AMDiS {
for (; !it.endDofIter(); it.nextDof())
if (dofMap.count(it.getDofIndex()))
if (globalIndex)
- sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
+ sendGlobalDofs.push_back(dofToMatIndex[level].get(i, dofMap[it.getDofIndex()].global));
else
- sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
+ sendGlobalDofs.push_back(dofToMatIndex[level].get(i, it.getDofIndex()));
stdMpi.send(it.getRank(), sendGlobalDofs);
}
@@ -380,9 +380,9 @@ namespace AMDiS {
if (dofMap.count(it.getDofIndex())) {
DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
if (globalIndex)
- dofToMatIndex.add(i, dofMap[it.getDofIndex()].global, d);
+ dofToMatIndex[level].add(i, dofMap[it.getDofIndex()].global, d);
else
- dofToMatIndex.add(i, it.getDofIndex(), d);
+ dofToMatIndex[level].add(i, it.getDofIndex(), d);
}
}
}
diff --git a/AMDiS/src/parallel/ParallelDofMapping.h b/AMDiS/src/parallel/ParallelDofMapping.h
index 1043fe7535f1f1a4404f7600e28a0a4c7fe40a7f..c4e945113e54cbf54f265362ad11971d6c085836 100644
--- a/AMDiS/src/parallel/ParallelDofMapping.h
+++ b/AMDiS/src/parallel/ParallelDofMapping.h
@@ -121,7 +121,8 @@ namespace AMDiS {
/// Clears all data of the mapping.
void clear();
-
+
+#if 0
/// Maps a DOF index to both, the local and global index of the mapping. The
/// global index must not be set.
MultiIndex& operator[](DegreeOfFreedom d)
@@ -130,6 +131,14 @@ namespace AMDiS {
return dofMap[0][d];
}
+#else
+ DofMap& operator[](int level)
+ {
+ TEST_EXIT_DBG(level < static_cast<int>(dofMap.size()))("Should not happen!\n");
+
+ return dofMap[level];
+ }
+#endif
/// Inserts a new DOF to rank's mapping. The DOF is assumed to be owend by
/// the rank.
@@ -159,7 +168,7 @@ namespace AMDiS {
}
/// Checks if a given DOF is in the DOF mapping.
- bool isSet(DegreeOfFreedom dof, int level = 0)
+ bool isSet(DegreeOfFreedom dof, int level)
{
return static_cast<bool>(dofMap[level].count(dof));
}
@@ -167,15 +176,15 @@ namespace AMDiS {
/// Checks if a given DOF is a rank owned DOF of the DOF mapping. The DOF must
/// a DOF of the mapping (this is not checked here), otherwise the result is
/// meaningsless.
- bool isRankDof(DegreeOfFreedom dof)
+ bool isRankDof(DegreeOfFreedom dof, int level)
{
- return !(static_cast<bool>(nonRankDofs[0].count(dof)));
+ return !(static_cast<bool>(nonRankDofs[level].count(dof)));
}
/// Returns number of DOFs in the mapping.
- unsigned int size()
+ unsigned int size(int level)
{
- return dofMap[0].size();
+ return dofMap[level].size();
}
/// Returns the raw data of the mapping.
@@ -362,21 +371,21 @@ namespace AMDiS {
/// Returns the global matrix index of a given DOF for a given
/// component number.
- inline int getMatIndex(int ithComponent, DegreeOfFreedom d)
+ inline int getMatIndex(int level, int ithComponent, DegreeOfFreedom d)
{
- return dofToMatIndex.get(ithComponent, d);
+ return dofToMatIndex[level].get(ithComponent, d);
}
/// Returns the local matrix index of a given DOF for a given
/// component number.
- inline int getLocalMatIndex(int ithComponent, DegreeOfFreedom d)
+ inline int getLocalMatIndex(int level, int ithComponent, DegreeOfFreedom d)
{
FUNCNAME("ParallelDofMapping::getLocalMatIndex()");
- TEST_EXIT_DBG(data[feSpaces[ithComponent]].isRankDof(d))
+ TEST_EXIT_DBG(data[feSpaces[ithComponent]].isRankDof(d, level))
("Should not happen!\n");
- return dofToMatIndex.get(ithComponent, d) - rStartDofs[0];
+ return dofToMatIndex[level].get(ithComponent, d) - rStartDofs[level];
}
// Writes all data of this object to an output stream.
@@ -455,7 +464,7 @@ namespace AMDiS {
/// Mapping from global DOF indices to global matrix indices under
/// consideration of possibly multiple components.
- DofToMatIndex dofToMatIndex;
+ vector<DofToMatIndex> dofToMatIndex;
};
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 4cfd74ae934f72e33a4c82ff2a9c1bf1c86dc057..9705970ce56d52a085a33033d59f0e00fb331b82 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -295,8 +295,10 @@ namespace AMDiS {
for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
- createPrimals(feSpace);
+ createPrimals(feSpace);
+
createDuals(feSpace);
+
createIndexB(feSpace);
}
@@ -306,6 +308,15 @@ namespace AMDiS {
if (fetiPreconditioner != FETI_NONE)
interiorDofMap.update();
+ for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
+ const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
+
+ MSG("TEST DUALS %d %d %d\n",
+ dualDofMap[feSpace].nLocalDofs[meshLevel],
+ dualDofMap[feSpace].nRankDofs[meshLevel],
+ dualDofMap[feSpace].nOverallDofs[meshLevel]);
+ }
+
for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
createLagrange(feSpace);
@@ -313,20 +324,24 @@ namespace AMDiS {
lagrangeMap.update();
+ MSG("FETI-DP data created on mesh level %d\n", meshLevel);
for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
MSG("FETI-DP data for %d-ith FE space:\n", i);
MSG(" nRankPrimals = %d nOverallPrimals = %d\n",
- primalDofMap[feSpace].nRankDofs[0], primalDofMap[feSpace].nOverallDofs[0]);
+ primalDofMap[feSpace].nRankDofs[meshLevel],
+ primalDofMap[feSpace].nOverallDofs[meshLevel]);
MSG(" nRankDuals = %d nOverallDuals = %d\n",
- dualDofMap[feSpace].nRankDofs[0], dualDofMap[feSpace].nOverallDofs[0]);
+ dualDofMap[feSpace].nRankDofs[meshLevel],
+ dualDofMap[feSpace].nOverallDofs[meshLevel]);
MSG(" nRankLagrange = %d nOverallLagrange = %d\n",
- lagrangeMap[feSpace].nRankDofs[0], lagrangeMap[feSpace].nOverallDofs[0]);
+ lagrangeMap[feSpace].nRankDofs[meshLevel],
+ lagrangeMap[feSpace].nOverallDofs[meshLevel]);
- TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() ==
+ TEST_EXIT_DBG(localDofMap[feSpace].size(meshLevel) + primalDofMap[feSpace].size(meshLevel) ==
static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
("Should not happen!\n");
}
@@ -358,10 +373,10 @@ namespace AMDiS {
// === create local indices of the primals starting at zero. ===
for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
- if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
- primalDofMap[feSpace].insertRankDof(0, *it);
+ if (meshDistributor->getDofMap()[feSpace].isRankDof(*it, meshLevel))
+ primalDofMap[feSpace].insertRankDof(meshLevel, *it);
else
- primalDofMap[feSpace].insertNonRankDof(0, *it);
+ primalDofMap[feSpace].insertNonRankDof(meshLevel, *it);
}
@@ -377,7 +392,7 @@ namespace AMDiS {
for (DofContainer::iterator it = allBoundaryDofs.begin();
it != allBoundaryDofs.end(); ++it)
if (!isPrimal(feSpace, **it))
- dualDofMap[feSpace].insertRankDof(0, **it);
+ dualDofMap[feSpace].insertRankDof(meshLevel, **it);
}
@@ -385,11 +400,14 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createLagrange()");
+ boundaryDofRanks[feSpace].clear();
+
+ if (dualDofMap[feSpace].nLocalDofs[meshLevel] == 0)
+ return;
+
// === Create for each dual node that is owned by the rank, the set ===
// === of ranks that contain this node (denoted by W(x_j)). ===
- boundaryDofRanks[feSpace].clear();
-
for (DofComm::Iterator it(meshDistributor->getSendDofs(), meshLevel, feSpace);
!it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof()) {
@@ -442,18 +460,18 @@ namespace AMDiS {
// === appropriate number of Lagrange constraints. ===
int nRankLagrange = 0;
- DofMap& dualMap = dualDofMap[feSpace].getMap(0);
+ DofMap& dualMap = dualDofMap[feSpace].getMap(meshLevel);
for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
- if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
- lagrangeMap[feSpace].insertRankDof(0, it->first, nRankLagrange);
+ if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first, meshLevel)) {
+ lagrangeMap[feSpace].insertRankDof(meshLevel, it->first, nRankLagrange);
int degree = boundaryDofRanks[feSpace][it->first].size();
nRankLagrange += (degree * (degree - 1)) / 2;
} else {
- lagrangeMap[feSpace].insertNonRankDof(0, it->first);
+ lagrangeMap[feSpace].insertNonRankDof(meshLevel, it->first);
}
}
- lagrangeMap[feSpace].nRankDofs[0] = nRankLagrange;
+ lagrangeMap[feSpace].nRankDofs[meshLevel] = nRankLagrange;
}
@@ -471,11 +489,11 @@ namespace AMDiS {
for (int i = 0; i < admin->getUsedSize(); i++) {
if (admin->isDofFree(i) == false &&
isPrimal(feSpace, i) == false &&
- dualDofMap[feSpace].isSet(i) == false) {
- localDofMap[feSpace].insertRankDof(0, i, nLocalInterior);
+ isDual(feSpace, i) == false) {
+ localDofMap[feSpace].insertRankDof(meshLevel, i, nLocalInterior);
if (fetiPreconditioner != FETI_NONE)
- interiorDofMap[feSpace].insertRankDof(0, i, nLocalInterior);
+ interiorDofMap[feSpace].insertRankDof(meshLevel, i, nLocalInterior);
nLocalInterior++;
}
@@ -483,9 +501,9 @@ namespace AMDiS {
// === And finally, add the global indicies of all dual nodes. ===
- for (DofMap::iterator it = dualDofMap[feSpace].getMap(0).begin();
- it != dualDofMap[feSpace].getMap(0).end(); ++it)
- localDofMap[feSpace].insertRankDof(0, it->first);
+ for (DofMap::iterator it = dualDofMap[feSpace].getMap(meshLevel).begin();
+ it != dualDofMap[feSpace].getMap(meshLevel).end(); ++it)
+ localDofMap[feSpace].insertRankDof(meshLevel, it->first);
}
@@ -496,8 +514,10 @@ namespace AMDiS {
// === Create distributed matrix for Lagrange constraints. ===
MatCreateMPIAIJ(mpiComm,
- lagrangeMap.getRankDofs(0), localDofMap.getRankDofs(0),
- lagrangeMap.getOverallDofs(0), localDofMap.getOverallDofs(0),
+ lagrangeMap.getRankDofs(meshLevel),
+ localDofMap.getRankDofs(meshLevel),
+ lagrangeMap.getOverallDofs(meshLevel),
+ localDofMap.getOverallDofs(meshLevel),
2, PETSC_NULL, 2, PETSC_NULL,
&mat_lagrange);
@@ -509,14 +529,14 @@ namespace AMDiS {
// === constraint. ===
for (unsigned int k = 0; k < feSpaces.size(); k++) {
- DofMap &dualMap = dualDofMap[feSpaces[k]].getMap(0);
+ DofMap &dualMap = dualDofMap[feSpaces[k]].getMap(meshLevel);
for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
("Should not happen!\n");
// Global index of the first Lagrange constriant for this node.
- int index = lagrangeMap.getMatIndex(k, it->first);
+ int index = lagrangeMap.getMatIndex(meshLevel, k, it->first);
// Copy set of all ranks that contain this dual node.
vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(),
@@ -527,7 +547,7 @@ namespace AMDiS {
for (int i = 0; i < degree; i++) {
for (int j = i + 1; j < degree; j++) {
if (W[i] == mpiRank || W[j] == mpiRank) {
- int colIndex = localDofMap.getMatIndex(k, it->first);
+ int colIndex = localDofMap.getMatIndex(meshLevel, k, it->first);
double value = (W[i] == mpiRank ? 1.0 : -1.0);
MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
}
@@ -544,7 +564,7 @@ namespace AMDiS {
void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
{
- FUNCNAME("PetscSolverFeti::createSchurPrimal()");
+ FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
if (schurPrimalSolver == 0) {
MSG("Create iterative schur primal solver!\n");
@@ -552,15 +572,19 @@ namespace AMDiS {
schurPrimalData.subSolver = subDomainSolver;
VecCreateMPI(mpiComm,
- localDofMap.getRankDofs(0), localDofMap.getOverallDofs(0),
+ localDofMap.getRankDofs(meshLevel),
+ localDofMap.getOverallDofs(meshLevel),
&(schurPrimalData.tmp_vec_b));
VecCreateMPI(mpiComm,
- primalDofMap.getRankDofs(0), primalDofMap.getOverallDofs(0),
+ primalDofMap.getRankDofs(meshLevel),
+ primalDofMap.getOverallDofs(meshLevel),
&(schurPrimalData.tmp_vec_primal));
MatCreateShell(mpiComm,
- primalDofMap.getRankDofs(0), primalDofMap.getRankDofs(0),
- primalDofMap.getOverallDofs(0), primalDofMap.getOverallDofs(0),
+ primalDofMap.getRankDofs(meshLevel),
+ primalDofMap.getRankDofs(meshLevel),
+ primalDofMap.getOverallDofs(meshLevel),
+ primalDofMap.getOverallDofs(meshLevel),
&schurPrimalData,
&mat_schur_primal);
MatShellSetOperation(mat_schur_primal, MATOP_MULT,
@@ -576,10 +600,10 @@ namespace AMDiS {
double wtime = MPI::Wtime();
- int nRowsRankPrimal = primalDofMap.getRankDofs(0);
- int nRowsOverallPrimal = primalDofMap.getOverallDofs(0);
- int nRowsRankB = localDofMap.getRankDofs(0);
- int nRowsOverallB = localDofMap.getOverallDofs(0);
+ int nRowsRankPrimal = primalDofMap.getRankDofs(meshLevel);
+ int nRowsOverallPrimal = primalDofMap.getOverallDofs(meshLevel);
+ int nRowsRankB = localDofMap.getRankDofs(meshLevel);
+ int nRowsOverallB = localDofMap.getOverallDofs(meshLevel);
Mat matBPi;
MatCreateMPIAIJ(mpiComm,
@@ -595,7 +619,7 @@ namespace AMDiS {
map<int, vector<pair<int, double> > > mat_b_primal_cols;
for (int i = 0; i < nRowsRankB; i++) {
- PetscInt row = localDofMap.getStartDofs(0) + i;
+ PetscInt row = localDofMap.getStartDofs(meshLevel) + i;
MatGetRow(subDomainSolver->getMatIntCoarse(), row, &nCols, &cols, &values);
for (int j = 0; j < nCols; j++)
@@ -606,7 +630,7 @@ namespace AMDiS {
}
TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) ==
- primalDofMap.getLocalDofs(0))
+ primalDofMap.getLocalDofs(meshLevel))
("Should not happen!\n");
for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
@@ -627,7 +651,7 @@ namespace AMDiS {
VecGetArray(tmpVec, &tmpValues);
for (int i = 0; i < nRowsRankB; i++)
MatSetValue(matBPi,
- localDofMap.getStartDofs(0) + i,
+ localDofMap.getStartDofs(meshLevel) + i,
it->first,
tmpValues[i],
ADD_VALUES);
@@ -694,18 +718,23 @@ namespace AMDiS {
fetiData.ksp_schur_primal = &ksp_schur_primal;
VecCreateMPI(mpiComm,
- localDofMap.getRankDofs(0), localDofMap.getOverallDofs(0),
+ localDofMap.getRankDofs(meshLevel),
+ localDofMap.getOverallDofs(meshLevel),
&(fetiData.tmp_vec_b));
VecCreateMPI(mpiComm,
- lagrangeMap.getRankDofs(0), lagrangeMap.getOverallDofs(0),
+ lagrangeMap.getRankDofs(meshLevel),
+ lagrangeMap.getOverallDofs(meshLevel),
&(fetiData.tmp_vec_lagrange));
VecCreateMPI(mpiComm,
- primalDofMap.getRankDofs(0), primalDofMap.getOverallDofs(0),
+ primalDofMap.getRankDofs(meshLevel),
+ primalDofMap.getOverallDofs(meshLevel),
&(fetiData.tmp_vec_primal));
MatCreateShell(mpiComm,
- lagrangeMap.getRankDofs(0), lagrangeMap.getRankDofs(0),
- lagrangeMap.getOverallDofs(0), lagrangeMap.getOverallDofs(0),
+ lagrangeMap.getRankDofs(meshLevel),
+ lagrangeMap.getRankDofs(meshLevel),
+ lagrangeMap.getOverallDofs(meshLevel),
+ lagrangeMap.getOverallDofs(meshLevel),
&fetiData, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);
@@ -746,7 +775,8 @@ namespace AMDiS {
fetiDirichletPreconData.ksp_interior = &ksp_interior;
VecCreateMPI(mpiComm,
- localDofMap.getRankDofs(0),localDofMap.getOverallDofs(0),
+ localDofMap.getRankDofs(meshLevel),
+ localDofMap.getOverallDofs(meshLevel),
&(fetiDirichletPreconData.tmp_vec_b));
MatGetVecs(mat_duals_duals, PETSC_NULL,
&(fetiDirichletPreconData.tmp_vec_duals0));
@@ -756,11 +786,11 @@ namespace AMDiS {
&(fetiDirichletPreconData.tmp_vec_interior));
for (unsigned int i = 0; i < feSpaces.size(); i++) {
- DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(0);
+ DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(meshLevel);
for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
DegreeOfFreedom d = it->first;
- int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
- int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
+ int matIndexLocal = localDofMap.getLocalMatIndex(meshLevel, i, d);
+ int matIndexDual = dualDofMap.getLocalMatIndex(meshLevel, i, d);
fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
}
}
@@ -777,18 +807,18 @@ namespace AMDiS {
fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;
for (unsigned int i = 0; i < feSpaces.size(); i++) {
- DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(0);
+ DofMap &dualMap = dualDofMap[feSpaces[i]].getMap(meshLevel);
for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
DegreeOfFreedom d = it->first;
- int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
- int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
+ int matIndexLocal = localDofMap.getLocalMatIndex(meshLevel, i, d);
+ int matIndexDual = dualDofMap.getLocalMatIndex(meshLevel, i, d);
fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
}
}
VecCreateMPI(mpiComm,
- localDofMap.getRankDofs(0),
- localDofMap.getOverallDofs(0),
+ localDofMap.getRankDofs(meshLevel),
+ localDofMap.getOverallDofs(meshLevel),
&(fetiLumpedPreconData.tmp_vec_b));
MatGetVecs(mat_duals_duals, PETSC_NULL,
&(fetiLumpedPreconData.tmp_vec_duals0));
@@ -876,20 +906,20 @@ namespace AMDiS {
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);
vector<PetscInt> globalIsIndex, localIsIndex;
- globalIsIndex.reserve(primalDofMap.getLocalDofs(0));
- localIsIndex.reserve(primalDofMap.getLocalDofs(0));
+ globalIsIndex.reserve(primalDofMap.getLocalDofs(meshLevel));
+ localIsIndex.reserve(primalDofMap.getLocalDofs(meshLevel));
{
int cnt = 0;
for (unsigned int i = 0; i < feSpaces.size(); i++) {
- DofMap& dofMap = primalDofMap[feSpaces[i]].getMap(0);
+ DofMap& dofMap = primalDofMap[feSpaces[i]].getMap(meshLevel);
for (DofMap::iterator it = dofMap.begin();it != dofMap.end(); ++it) {
- globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
+ globalIsIndex.push_back(primalDofMap.getMatIndex(meshLevel, i, it->first));
localIsIndex.push_back(cnt++);
}
}
- TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs(0))
+ TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs(meshLevel))
("Should not happen!\n");
}
@@ -929,14 +959,14 @@ namespace AMDiS {
for (int i = 0; i < vec.getSize(); i++) {
DOFVector<double>& dofVec = *(vec.getDOFVector(i));
- for (DofMap::iterator it = localDofMap[feSpaces[i]].getMap(0).begin();
- it != localDofMap[feSpaces[i]].getMap(0).end(); ++it) {
- int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
+ for (DofMap::iterator it = localDofMap[feSpaces[i]].getMap(meshLevel).begin();
+ it != localDofMap[feSpaces[i]].getMap(meshLevel).end(); ++it) {
+ int petscIndex = localDofMap.getLocalMatIndex(meshLevel, i, it->first);
dofVec[it->first] = localSolB[petscIndex];
}
- for (DofMap::iterator it = primalDofMap[feSpaces[i]].getMap(0).begin();
- it != primalDofMap[feSpaces[i]].getMap(0).end(); ++it)
+ for (DofMap::iterator it = primalDofMap[feSpaces[i]].getMap(meshLevel).begin();
+ it != primalDofMap[feSpaces[i]].getMap(meshLevel).end(); ++it)
dofVec[it->first] = localSolPrimal[cnt++];
}
@@ -960,25 +990,20 @@ namespace AMDiS {
// === Create matrices for the FETI-DP method. ===
- int nRowsRankB = localDofMap.getRankDofs(0);
- int nRowsOverallB = localDofMap.getOverallDofs(0);
- int nRowsRankPrimal = primalDofMap.getRankDofs(0);
- int nRowsOverallPrimal = primalDofMap.getOverallDofs(0);
-
subDomainSolver->fillPetscMatrix(mat);
// === Create matrices for FETI-DP preconditioner. ===
if (fetiPreconditioner != FETI_NONE) {
- int nRowsDual = dualDofMap.getRankDofs(0);
+ int nRowsDual = dualDofMap.getRankDofs(meshLevel);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsDual, nRowsDual, 60, PETSC_NULL,
&mat_duals_duals);
if (fetiPreconditioner == FETI_DIRICHLET) {
- int nRowsInterior = interiorDofMap.getRankDofs(0);
+ int nRowsInterior = interiorDofMap.getRankDofs(meshLevel);
MatCreateSeqAIJ(PETSC_COMM_SELF,
nRowsInterior, nRowsInterior, 60, PETSC_NULL,
@@ -992,10 +1017,6 @@ namespace AMDiS {
nRowsDual, nRowsInterior, 60, PETSC_NULL,
&mat_duals_interior);
}
- }
-
-
- if (fetiPreconditioner != FETI_NONE) {
// === Prepare traverse of sequentially created matrices. ===
@@ -1047,21 +1068,21 @@ namespace AMDiS {
if (!rowPrimal && !colPrimal) {
if (!isDual(feSpaces[i], *cursor)) {
if (!isDual(feSpaces[j], col(*icursor))) {
- int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor));
+ int colIndex = interiorDofMap.getLocalMatIndex(meshLevel, j, col(*icursor));
colsLocal.push_back(colIndex);
valuesLocal.push_back(value(*icursor));
} else {
- int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor));
+ int colIndex = dualDofMap.getLocalMatIndex(meshLevel, j, col(*icursor));
colsLocalOther.push_back(colIndex);
valuesLocalOther.push_back(value(*icursor));
}
} else {
if (!isDual(feSpaces[j], col(*icursor))) {
- int colIndex = interiorDofMap.getLocalMatIndex(j, col(*icursor));
+ int colIndex = interiorDofMap.getLocalMatIndex(meshLevel, j, col(*icursor));
colsLocalOther.push_back(colIndex);
valuesLocalOther.push_back(value(*icursor));
} else {
- int colIndex = dualDofMap.getLocalMatIndex(j, col(*icursor));
+ int colIndex = dualDofMap.getLocalMatIndex(meshLevel, j, col(*icursor));
colsLocal.push_back(colIndex);
valuesLocal.push_back(value(*icursor));
}
@@ -1076,14 +1097,14 @@ namespace AMDiS {
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
if (!isDual(feSpaces[i], *cursor)) {
- int rowIndex = interiorDofMap.getLocalMatIndex(i, *cursor);
+ int rowIndex = interiorDofMap.getLocalMatIndex(meshLevel, i, *cursor);
MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
- int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
+ int rowIndex = dualDofMap.getLocalMatIndex(meshLevel, i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
@@ -1094,7 +1115,7 @@ namespace AMDiS {
case FETI_LUMPED:
if (isDual(feSpaces[i], *cursor)) {
- int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
+ int rowIndex = dualDofMap.getLocalMatIndex(meshLevel, i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
}
@@ -1171,13 +1192,17 @@ namespace AMDiS {
Vec tmp_b0, tmp_b1, tmp_lagrange0, tmp_primal0, tmp_primal1;
VecCreateMPI(mpiComm,
- localDofMap.getRankDofs(0), localDofMap.getOverallDofs(0), &tmp_b0);
+ localDofMap.getRankDofs(meshLevel),
+ localDofMap.getOverallDofs(meshLevel), &tmp_b0);
VecCreateMPI(mpiComm,
- localDofMap.getRankDofs(0), localDofMap.getOverallDofs(0), &tmp_b1);
+ localDofMap.getRankDofs(meshLevel),
+ localDofMap.getOverallDofs(meshLevel), &tmp_b1);
VecCreateMPI(mpiComm,
- primalDofMap.getRankDofs(0), primalDofMap.getOverallDofs(0), &tmp_primal0);
+ primalDofMap.getRankDofs(meshLevel),
+ primalDofMap.getOverallDofs(meshLevel), &tmp_primal0);
VecCreateMPI(mpiComm,
- primalDofMap.getRankDofs(0), primalDofMap.getOverallDofs(0), &tmp_primal1);
+ primalDofMap.getRankDofs(meshLevel),
+ primalDofMap.getOverallDofs(meshLevel), &tmp_primal1);
MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0);
MatGetVecs(mat_lagrange, PETSC_NULL, &vec_rhs);
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index e640e7d699e145f04e3e624cfc168500d1fb15f1..e1dfc73121f82296fdee7560fe25846006ea2d02 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -171,15 +171,17 @@ namespace AMDiS {
void printStatistics();
/// Checks whether a given DOF in a given FE space is a primal DOF.
- inline bool isPrimal(const FiniteElemSpace *feSpace, DegreeOfFreedom dof)
+ inline bool isPrimal(const FiniteElemSpace *feSpace,
+ DegreeOfFreedom dof)
{
- return primalDofMap[feSpace].isSet(dof);
+ return primalDofMap[feSpace].isSet(dof, meshLevel);
}
/// Checks whether a given DOF in a give FE space is a dual DOF.
- inline bool isDual(const FiniteElemSpace *feSpace, DegreeOfFreedom dof)
+ inline bool isDual(const FiniteElemSpace *feSpace,
+ DegreeOfFreedom dof)
{
- return dualDofMap[feSpace].isSet(dof);
+ return dualDofMap[feSpace].isSet(dof, meshLevel);
}
protected:
diff --git a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
index 74026bb680381414deb4baa8f4fbea734fcb50b8..5a74eaae43d972d2f5cc4d9d1a1f58ac13c46618 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalBlockMatrix.cc
@@ -232,7 +232,7 @@ namespace AMDiS {
cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {
// Global index of the current row DOF.
- int rowIndex = (*dofMap)[feSpace][*cursor].global + dispRowIndex;
+ int rowIndex = (*dofMap)[feSpace][0][*cursor].global + dispRowIndex;
cols.clear();
values.clear();
@@ -240,7 +240,7 @@ namespace AMDiS {
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
// Global index of the current column index.
- int colIndex = (*dofMap)[feSpace][col(*icursor)].global + dispColIndex;
+ int colIndex = (*dofMap)[feSpace][0][col(*icursor)].global + dispColIndex;
// Ignore all zero entries, expect it is a diagonal entry.
if (value(*icursor) == 0.0 && rowIndex != colIndex)
@@ -267,7 +267,7 @@ namespace AMDiS {
// Traverse all used DOFs in the dof vector.
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
- int index = (*dofMap)[feSpace][dofIt.getDOFIndex()].global;
+ int index = (*dofMap)[feSpace][0][dofIt.getDOFIndex()].global;
double value = *dofIt;
VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 58eaaa653174c32ee183d057b110b7623aa48979..f1dafc82d61dc239522763eca8d25398eae08396 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -268,7 +268,7 @@ namespace AMDiS {
const FiniteElemSpace *colFe = mat->getColFeSpace();
// Global index of the current row DOF.
- int globalRowDof = (*dofMap)[rowFe][*cursor].global;
+ int globalRowDof = (*dofMap)[rowFe][0][*cursor].global;
// Test if the current row DOF is a periodic DOF.
bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof);
@@ -277,7 +277,7 @@ namespace AMDiS {
// === Row DOF index is not periodic. ===
// Get PETSc's mat row index.
- int rowIndex = dofMap->getMatIndex(nRowMat, globalRowDof);
+ int rowIndex = dofMap->getMatIndex(0, nRowMat, globalRowDof);
cols.clear();
values.clear();
@@ -286,11 +286,11 @@ namespace AMDiS {
icursor != icend; ++icursor) {
// Global index of the current column index.
- int globalColDof = (*dofMap)[colFe][col(*icursor)].global;
+ int globalColDof = (*dofMap)[colFe][0][col(*icursor)].global;
// Test if the current col dof is a periodic dof.
bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
// Get PETSc's mat col index.
- int colIndex = dofMap->getMatIndex(nColMat, globalColDof);
+ int colIndex = dofMap->getMatIndex(0, nColMat, globalColDof);
// Ignore all zero entries, expect it is a diagonal entry.
if (value(*icursor) == 0.0 && rowIndex != colIndex)
@@ -339,7 +339,7 @@ namespace AMDiS {
}
for (unsigned int i = 0; i < newCols.size(); i++) {
- cols.push_back(dofMap->getMatIndex(nColMat, newCols[i]));
+ cols.push_back(dofMap->getMatIndex(0, nColMat, newCols[i]));
values.push_back(scaledValue);
}
}
@@ -361,7 +361,7 @@ namespace AMDiS {
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
// Global index of the current column index.
- int globalColDof = (*dofMap)[colFe][col(*icursor)].global;
+ int globalColDof = (*dofMap)[colFe][0][col(*icursor)].global;
// Ignore all zero entries, expect it is a diagonal entry.
if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
@@ -426,8 +426,8 @@ namespace AMDiS {
// === Translate the matrix entries to PETSc's matrix.
for (unsigned int i = 0; i < entry.size(); i++) {
- int rowIdx = dofMap->getMatIndex(nRowMat, entry[i].first);
- int colIdx = dofMap->getMatIndex(nColMat, entry[i].second);
+ int rowIdx = dofMap->getMatIndex(0, nRowMat, entry[i].first);
+ int colIdx = dofMap->getMatIndex(0, nColMat, entry[i].second);
colsMap[rowIdx].push_back(colIdx);
valsMap[rowIdx].push_back(scaledValue);
@@ -465,14 +465,14 @@ namespace AMDiS {
// Traverse all used DOFs in the dof vector.
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
- if (rankOnly && !(*dofMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
+ if (rankOnly && !(*dofMap)[feSpace].isRankDof(dofIt.getDOFIndex(), 0))
continue;
// Calculate global row index of the DOF.
DegreeOfFreedom globalRowDof =
- (*dofMap)[feSpace][dofIt.getDOFIndex()].global;
+ (*dofMap)[feSpace][0][dofIt.getDOFIndex()].global;
// Get PETSc's mat index of the row DOF.
- int index = dofMap->getMatIndex(nRowVec, globalRowDof);
+ int index = dofMap->getMatIndex(0, nRowVec, globalRowDof);
if (perMap.isPeriodic(feSpace, globalRowDof)) {
std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
@@ -482,7 +482,7 @@ namespace AMDiS {
for (std::set<int>::iterator perIt = perAsc.begin();
perIt != perAsc.end(); ++perIt) {
int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
- int mappedIndex = dofMap->getMatIndex(nRowVec, mappedDof);
+ int mappedIndex = dofMap->getMatIndex(0, nRowVec, mappedDof);
VecSetValues(petscVec, 1, &mappedIndex, &value, ADD_VALUES);
}
} else {
@@ -572,11 +572,11 @@ namespace AMDiS {
for (cursor_type cursor = begin<row>(bmat),
cend = end<row>(bmat); cursor != cend; ++cursor) {
- int globalRowDof = (*dofMap)[feSpaces[i]][*cursor].global;
+ int globalRowDof = (*dofMap)[feSpaces[i]][0][*cursor].global;
// The corresponding global matrix row index of the current row DOF.
- int petscRowIdx = dofMap->getMatIndex(i, globalRowDof);
- if ((*dofMap)[feSpaces[i]].isRankDof(*cursor)) {
+ int petscRowIdx = dofMap->getMatIndex(0, i, globalRowDof);
+ if ((*dofMap)[feSpaces[i]].isRankDof(*cursor, 0)) {
// === The current row DOF is a rank DOF, so create the ===
// === corresponding nnz values directly on rank's nnz data. ===
@@ -587,7 +587,7 @@ namespace AMDiS {
TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
("Should not happen! \n Debug info: localRowIdx = %d globalRowIndx = %d petscRowIdx = %d localPetscRowIdx = %d rStart = %d nCompontens = %d nRankRows = %d\n",
*cursor,
- (*dofMap)[feSpaces[i]][*cursor].global,
+ (*dofMap)[feSpaces[i]][0][*cursor].global,
petscRowIdx,
localPetscRowIdx,
rankStartIndex,
@@ -598,8 +598,8 @@ namespace AMDiS {
// Traverse all non zero entries in this row.
for (icursor_type icursor = begin<nz>(cursor),
icend = end<nz>(cursor); icursor != icend; ++icursor) {
- int globalColDof = (*dofMap)[feSpaces[j]][col(*icursor)].global;
- int petscColIdx = dofMap->getMatIndex(j, globalColDof);
+ int globalColDof = (*dofMap)[feSpaces[j]][0][col(*icursor)].global;
+ int petscColIdx = dofMap->getMatIndex(0, j, globalColDof);
if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) {
// The row DOF is a rank DOF, if also the column is a rank DOF,
@@ -626,8 +626,8 @@ namespace AMDiS {
icend = end<nz>(cursor); icursor != icend; ++icursor) {
if (value(*icursor) != 0.0) {
int globalColDof =
- (*dofMap)[feSpaces[j]][col(*icursor)].global;
- int petscColIdx = dofMap->getMatIndex(j, globalColDof);
+ (*dofMap)[feSpaces[j]][0][col(*icursor)].global;
+ int petscColIdx = dofMap->getMatIndex(0, j, globalColDof);
sendMatrixEntry[sendToRank].
push_back(make_pair(petscRowIdx, petscColIdx));
diff --git a/AMDiS/src/parallel/PetscSolverSchur.cc b/AMDiS/src/parallel/PetscSolverSchur.cc
index 1ef550dd72929d38ad5e3848c25dbc47f403f91c..38813de583e526909a9eabb60fbfe753b1f3e493 100644
--- a/AMDiS/src/parallel/PetscSolverSchur.cc
+++ b/AMDiS/src/parallel/PetscSolverSchur.cc
@@ -35,7 +35,7 @@ namespace AMDiS {
!it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof()) {
boundaryLocalDofs.insert(it.getDofIndex());
- boundaryDofs.insert((*dofMap)[feSpace][it.getDofIndex()].global);
+ boundaryDofs.insert((*dofMap)[feSpace][0][it.getDofIndex()].global);
}
@@ -45,9 +45,9 @@ namespace AMDiS {
DofContainerSet& edgeDofs =
- meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[EDGE];
+ meshDistributor->getBoundaryDofInfo(feSpace, 0).geoDofs[EDGE];
DofContainerSet& vertexDofs =
- meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
+ meshDistributor->getBoundaryDofInfo(feSpace, 0).geoDofs[VERTEX];
int nEdgeDofs = edgeDofs.size();
int nVertexDofs = vertexDofs.size();
@@ -72,14 +72,14 @@ namespace AMDiS {
int counter = rStartEdgeDofs;
for (DofContainerSet::iterator it = edgeDofs.begin();
it != edgeDofs.end(); ++it)
- mapGlobalBoundaryDof[(*dofMap)[feSpace][**it].global] =
+ mapGlobalBoundaryDof[(*dofMap)[feSpace][0][**it].global] =
counter++;
}
{
int counter = nOverallEdgeDofs + rStartVertexDofs;
for (DofContainerSet::iterator it = vertexDofs.begin();
it != vertexDofs.end(); ++it)
- mapGlobalBoundaryDof[(*dofMap)[feSpace][**it].global] =
+ mapGlobalBoundaryDof[(*dofMap)[feSpace][0][**it].global] =
counter++;
}
#else
@@ -135,7 +135,7 @@ namespace AMDiS {
stdMpi.getSendData(it.getRank()).reserve(it.getDofs().size());
for (; !it.endDofIter(); it.nextDof()) {
- int globalSendDof = (*dofMap)[feSpace][it.getDofIndex()].global;
+ int globalSendDof = (*dofMap)[feSpace][0][it.getDofIndex()].global;
TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalSendDof))
("No mapping for boundary DOF %d!\n", globalSendDof);
@@ -155,7 +155,7 @@ namespace AMDiS {
for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
!it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof()) {
- int globalRecvDof = (*dofMap)[feSpace][it.getDofIndex()].global;
+ int globalRecvDof = (*dofMap)[feSpace][0][it.getDofIndex()].global;
mapGlobalBoundaryDof[globalRecvDof] =
stdMpi.getRecvData(it.getRank())[it.getDofCounter()];
boundaryDofs.insert(globalRecvDof);
@@ -302,7 +302,7 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
DOFVector<double>::Iterator dofIt(vec.getDOFVector(i), USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
- DegreeOfFreedom globalRowDof = (*dofMap)[feSpace][dofIt.getDOFIndex()].global;
+ DegreeOfFreedom globalRowDof = (*dofMap)[feSpace][0][dofIt.getDOFIndex()].global;
if (boundaryDofs.count(globalRowDof)) {
int index =
(mapGlobalBoundaryDof[globalRowDof] - rStartBoundaryDofs + nInteriorDofs) * (i + 1);
@@ -373,7 +373,7 @@ namespace AMDiS {
cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {
// Global index of the current row DOF.
- int globalRowDof = (*dofMap)[feSpace][*cursor].global;
+ int globalRowDof = (*dofMap)[feSpace][0][*cursor].global;
colsBoundary.clear();
colsInterior.clear();
@@ -382,7 +382,7 @@ namespace AMDiS {
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
- int globalColDof = (*dofMap)[feSpace][col(*icursor)].global;
+ int globalColDof = (*dofMap)[feSpace][0][col(*icursor)].global;
if (boundaryDofs.count(globalColDof)) {
TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalColDof))
@@ -441,12 +441,12 @@ namespace AMDiS {
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
- if (rankOnly && !(*dofMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
+ if (rankOnly && !(*dofMap)[feSpace].isRankDof(dofIt.getDOFIndex(), 0))
continue;
// Calculate global row index of the DOF.
DegreeOfFreedom globalRowDof =
- (*dofMap)[feSpace][dofIt.getDOFIndex()].global;
+ (*dofMap)[feSpace][0][dofIt.getDOFIndex()].global;
double value = *dofIt;
if (boundaryDofs.count(globalRowDof)) {
diff --git a/AMDiS/src/parallel/SubDomainSolver.cc b/AMDiS/src/parallel/SubDomainSolver.cc
index c9fa411a379773e801b8be2f3d2ca5dc5abfb04f..9d0f7a879b3a8a87f19ebbc518b3208d3c021382 100644
--- a/AMDiS/src/parallel/SubDomainSolver.cc
+++ b/AMDiS/src/parallel/SubDomainSolver.cc
@@ -29,6 +29,10 @@ namespace AMDiS {
int nRowsRankCoarse = coarseSpaceMap->getRankDofs(level);
int nRowsOverallCoarse = coarseSpaceMap->getOverallDofs(level);
+ MSG("FILL MATRIX IN MPI COMM WITH SIZES %d %d\n",
+ mpiCommInterior.Get_size(),
+ mpiCommCoarseSpace.Get_size());
+
if (mpiCommInterior.Get_size() == 1) {
MatCreateSeqAIJ(mpiCommInterior, nRowsRankInterior, nRowsRankInterior,
60, PETSC_NULL, &matIntInt);
@@ -122,33 +126,33 @@ namespace AMDiS {
// === Set matrix values. ===
if (rowPrimal) {
- int rowIndex = coarseSpaceMap->getMatIndex(i, *cursor);
+ int rowIndex = coarseSpaceMap->getMatIndex(level, i, *cursor);
for (unsigned int k = 0; k < cols.size(); k++)
- cols[k] = coarseSpaceMap->getMatIndex(j, cols[k]);
+ cols[k] = coarseSpaceMap->getMatIndex(level, j, cols[k]);
MatSetValues(matCoarseCoarse, 1, &rowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
if (colsOther.size()) {
for (unsigned int k = 0; k < colsOther.size(); k++)
- colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]);
+ colsOther[k] = interiorMap->getMatIndex(level, j, colsOther[k]);
MatSetValues(matCoarseInt, 1, &rowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
}
} else {
- int localRowIndex = interiorMap->getLocalMatIndex(i, *cursor);
+ int localRowIndex = interiorMap->getLocalMatIndex(level, i, *cursor);
for (unsigned int k = 0; k < cols.size(); k++)
- cols[k] = interiorMap->getLocalMatIndex(j, cols[k]);
+ cols[k] = interiorMap->getLocalMatIndex(level, j, cols[k]);
MatSetValues(matIntInt, 1, &localRowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
if (colsOther.size()) {
- int globalRowIndex = interiorMap->getMatIndex(i, *cursor);
+ int globalRowIndex = interiorMap->getMatIndex(level, i, *cursor);
for (unsigned int k = 0; k < colsOther.size(); k++)
- colsOther[k] = coarseSpaceMap->getMatIndex(j, colsOther[k]);
+ colsOther[k] = coarseSpaceMap->getMatIndex(level, j, colsOther[k]);
MatSetValues(matIntCoarse, 1, &globalRowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
@@ -207,10 +211,10 @@ namespace AMDiS {
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
int index = dofIt.getDOFIndex();
if (isCoarseSpace(feSpace, index)) {
- index = coarseSpaceMap->getMatIndex(i, index);
+ index = coarseSpaceMap->getMatIndex(level, i, index);
VecSetValue(rhsCoarseSpace, index, *dofIt, ADD_VALUES);
} else {
- index = interiorMap->getMatIndex(i, index);
+ index = interiorMap->getMatIndex(level, i, index);
VecSetValue(rhsInterior, index, *dofIt, INSERT_VALUES);
}
}
diff --git a/AMDiS/src/parallel/SubDomainSolver.h b/AMDiS/src/parallel/SubDomainSolver.h
index 52df9910a23d505b2349f687fa1450b17a3a4fa1..3e9eb8258ff7a10b8d01a381484e74f7c356307b 100644
--- a/AMDiS/src/parallel/SubDomainSolver.h
+++ b/AMDiS/src/parallel/SubDomainSolver.h
@@ -73,9 +73,10 @@ namespace AMDiS {
void solveGlobal(Vec &rhs, Vec &sol);
- inline bool isCoarseSpace(const FiniteElemSpace *feSpace, DegreeOfFreedom dof)
+ inline bool isCoarseSpace(const FiniteElemSpace *feSpace,
+ DegreeOfFreedom dof)
{
- return (*coarseSpaceMap)[feSpace].isSet(dof);
+ return (*coarseSpaceMap)[feSpace].isSet(dof, level);
}
inline Vec& getRhsCoarseSpace()