diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index f140ca49fa53d2e0c329850e933ac15399709d29..5d3e7f499f2691b49158feb80e55dc0e0fed999c 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -333,9 +333,9 @@ namespace AMDiS {
// Go through all neighbours of the macro element and remove this macro element
// to be neighbour of some other macro element.
- for (int i = 0; i <= dim; i++)
+ for (int i = 0; i < getGeo(NEIGH); i++)
if ((*macroIt)->getNeighbour(i))
- for (int j = 0; j <= dim; j++)
+ for (int j = 0; j < getGeo(NEIGH); j++)
if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
(*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
@@ -352,13 +352,13 @@ namespace AMDiS {
// We will delete at least some element DOFs in the next but will keep the
// element object still in memory. Therefore the DOF pointer must be set to be
- // invalide.
+ // invalid.
(*macroIt)->getElement()->setDofValid(false);
}
- // === Check now all the dofs, that have no owner anymore and therefore have ===
- // === to be removed. ===
+ // === Check now all the DOFs that have no owner anymore and therefore have ===
+ // === to be removed. ===
for (DofElMap::iterator dofsIt = dofsOwner.begin();
dofsIt != dofsOwner.end(); ++dofsIt) {
diff --git a/AMDiS/src/io/DataCollector.cc b/AMDiS/src/io/DataCollector.cc
index 6c20185634e8dc775450ea5d4a54701dd674a99a..6125ec61af63a64018257b4a47dd38e6398d0bf2 100644
--- a/AMDiS/src/io/DataCollector.cc
+++ b/AMDiS/src/io/DataCollector.cc
@@ -217,11 +217,11 @@ namespace AMDiS {
FUNCNAME("DataCollector::addElementData()");
const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
- DegreeOfFreedom vertexDOF;
- WorldVector<double> vertexCoords;
- // MSG("ELEMENT: %d\n", elInfo->getElement()->getIndex());
- // MSG("DOFs: %d %d %d\n", dof[0][0], dof[1][0], dof[2][0]);
+// MSG("ELEMENT: %d %d\n",
+// elInfo->getElement()->getIndex(),
+// elInfo->getMacroElement()->getIndex());
+// MSG("DOFs: %d %d %d %d\n", dof[0][0], dof[1][0], dof[2][0], dof[3][0]);
// create ElementInfo
ElementInfo elementInfo(dim);
@@ -244,28 +244,28 @@ namespace AMDiS {
}
// for all vertices
- for (int i = 0; i < dim + 1; i++) {
+ for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
// get coords of this vertex
- vertexCoords = elInfo->getCoord(i);
+ WorldVector<double> &vertexCoords = elInfo->getCoord(i);
// get dof index of this vertex
- vertexDOF = dof[i][nPreDofs];
+ DegreeOfFreedom vertexDof = dof[i][nPreDofs];
// search for coords at this dof
std::list<VertexInfo>::iterator it =
- find((*vertexInfos)[vertexDOF].begin(), (*vertexInfos)[vertexDOF].end(),
+ find((*vertexInfos)[vertexDof].begin(), (*vertexInfos)[vertexDof].end(),
vertexCoords);
// coords not yet in list?
- if (it == (*vertexInfos)[vertexDOF].end()) {
+ if (it == (*vertexInfos)[vertexDof].end()) {
// create new vertex info and increment number of vertices
VertexInfo newVertexInfo = {vertexCoords, nVertices};
// add new vertex info to list
- (*vertexInfos)[vertexDOF].push_front(newVertexInfo);
+ (*vertexInfos)[vertexDof].push_front(newVertexInfo);
// set iterator to new vertex info
- it = (*vertexInfos)[vertexDOF].begin();
+ it = (*vertexInfos)[vertexDof].begin();
nVertices++;
}
diff --git a/AMDiS/src/io/VtkWriter.cc b/AMDiS/src/io/VtkWriter.cc
index 41db97c3cf895a51a82a5c16fac4f9180caa9cb1..1d0503610b588b7a22d0550fd8db1a02c027e4b3 100644
--- a/AMDiS/src/io/VtkWriter.cc
+++ b/AMDiS/src/io/VtkWriter.cc
@@ -56,16 +56,6 @@ namespace AMDiS {
file.push(boost::iostreams::file_descriptor_sink(name, std::ios::trunc));
writeFileToStream(file);
-#if 0
-
- std::ofstream file;
- file.open(name.c_str());
- TEST_EXIT(file.is_open())("Cannot open file %s for writing\n", name.c_str());
- writeFileToStream(file);
- file.close();
-
-#endif
-
return 0;
}
@@ -159,7 +149,7 @@ namespace AMDiS {
DataCollector dc(values->getFeSpace(), values);
std::vector<DataCollector*> dcList(0);
dcList.push_back(&dc);
- writeFile(dcList,filename,writeParallel);
+ writeFile(dcList, filename, writeParallel);
}
diff --git a/AMDiS/src/io/VtkWriter.hh b/AMDiS/src/io/VtkWriter.hh
index 473c028ef4d3753494ef26fb52f3bd4be5a90eb2..214e1e83d8a9aad0eec917267ae50f1c62cf942b 100644
--- a/AMDiS/src/io/VtkWriter.hh
+++ b/AMDiS/src/io/VtkWriter.hh
@@ -19,6 +19,8 @@ namespace AMDiS {
template<typename T>
void VtkWriter::writeFileToStream(T &file)
{
+ FUNCNAME("VtkWriter::writeFileToStream()");
+
int nVertices = (*dataCollector)[0]->getNumberVertices();
int nElements = (*dataCollector)[0]->getNumberElements();
int vertices = (*dataCollector)[0]->getMesh()->getGeo(VERTEX);
diff --git a/AMDiS/src/parallel/ElementObjectData.cc b/AMDiS/src/parallel/ElementObjectData.cc
index a97b17e5c38dc514850c2af2ce2c801ac7670f0f..708574506697a05d82d2a9775b382f3bf308049e 100644
--- a/AMDiS/src/parallel/ElementObjectData.cc
+++ b/AMDiS/src/parallel/ElementObjectData.cc
@@ -336,7 +336,7 @@ namespace AMDiS {
FUNCNAME("ElementObjects::provideConnectedPeriodicBoundary()");
std::pair<BoundaryType, BoundaryType> bConn =
- (b0 <= b1 ? std::make_pair(b0, b1) : std::make_pair(b1, b0));
+ (b0 <= b1 ? make_pair(b0, b1) : make_pair(b1, b0));
if (bConnMap.count(bConn) == 0) {
BoundaryType newPeriodicBoundaryType = getNewBoundaryType();
@@ -417,6 +417,65 @@ namespace AMDiS {
}
+ void ElementObjects::createReverseModeData(map<int, Element*> &elIndexMap,
+ map<int, int> &elIndexTypeMap)
+ {
+ FUNCNAME("ElementObjects::createReverseModeData()");
+
+ // In 2D, all reverse modes are always true!
+ if (mesh->getDim() == 2)
+ return;
+
+ for (map<DofEdge, vector<ElementObjectData> >::iterator edgeIt = edgeElements.begin();
+ edgeIt != edgeElements.end(); ++edgeIt) {
+
+ vector<ElementObjectData>& els = edgeIt->second;
+
+ for (unsigned int i = 0; i < els.size(); i++) {
+ BoundaryObject obj0(elIndexMap[els[i].elIndex],
+ elIndexTypeMap[els[i].elIndex],
+ EDGE, els[i].ithObject);
+
+ for (unsigned int j = i + 1; j < els.size(); j++) {
+ BoundaryObject obj1(elIndexMap[els[j].elIndex],
+ elIndexTypeMap[els[j].elIndex],
+ EDGE, els[j].ithObject);
+
+ bool reverseMode =
+ BoundaryObject::computeReverseMode(obj0, obj1, feSpace, INTERIOR);
+
+ edgeReverseMode[make_pair(els[i], els[j])] = reverseMode;
+ edgeReverseMode[make_pair(els[j], els[i])] = reverseMode;
+ }
+ }
+ }
+
+ for (map<DofFace, vector<ElementObjectData> >::iterator faceIt = faceElements.begin();
+ faceIt != faceElements.end(); ++faceIt) {
+
+ vector<ElementObjectData>& els = faceIt->second;
+
+ for (unsigned int i = 0; i < els.size(); i++) {
+ BoundaryObject obj0(elIndexMap[els[i].elIndex],
+ elIndexTypeMap[els[i].elIndex],
+ FACE, els[i].ithObject);
+
+ for (unsigned int j = i + 1; j < els.size(); j++) {
+ BoundaryObject obj1(elIndexMap[els[j].elIndex],
+ elIndexTypeMap[els[j].elIndex],
+ FACE, els[j].ithObject);
+
+ bool reverseMode =
+ BoundaryObject::computeReverseMode(obj0, obj1, feSpace, INTERIOR);
+
+ faceReverseMode[make_pair(els[i], els[j])] = reverseMode;
+ faceReverseMode[make_pair(els[j], els[i])] = reverseMode;
+ }
+ }
+ }
+ }
+
+
void ElementObjects::serialize(ostream &out)
{
FUNCNAME("ElementObjects::serialize()");
@@ -505,6 +564,42 @@ namespace AMDiS {
SerUtil::serialize(out, periodicVertices);
SerUtil::serialize(out, periodicEdges);
SerUtil::serialize(out, periodicFaces);
+
+
+ nSize = periodicDofAssoc.size();
+ SerUtil::serialize(out, nSize);
+ for (map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
+ it != periodicDofAssoc.end(); ++it) {
+ SerUtil::serialize(out, it->first);
+ SerUtil::serialize(out, it->second);
+ }
+
+ nSize = periodicEdgeAssoc.size();
+ SerUtil::serialize(out, nSize);
+ for (map<DofEdge, std::set<DofEdge> >::iterator it = periodicEdgeAssoc.begin();
+ it != periodicEdgeAssoc.end(); ++it) {
+ SerUtil::serialize(out, it->first);
+ SerUtil::serialize(out, it->second);
+ }
+
+
+ nSize = edgeReverseMode.size();
+ SerUtil::serialize(out, nSize);
+ for (map<pair<ElementObjectData, ElementObjectData>, bool>::iterator it = edgeReverseMode.begin();
+ it != edgeReverseMode.end(); ++it) {
+ it->first.first.serialize(out);
+ it->first.second.serialize(out);
+ SerUtil::serialize(out, it->second);
+ }
+
+ nSize = faceReverseMode.size();
+ SerUtil::serialize(out, nSize);
+ for (map<pair<ElementObjectData, ElementObjectData>, bool>::iterator it = faceReverseMode.begin();
+ it != faceReverseMode.end(); ++it) {
+ it->first.first.serialize(out);
+ it->first.second.serialize(out);
+ SerUtil::serialize(out, it->second);
+ }
}
@@ -613,7 +708,55 @@ namespace AMDiS {
SerUtil::deserialize(in, periodicVertices);
SerUtil::deserialize(in, periodicEdges);
- SerUtil::deserialize(in, periodicFaces);
+ SerUtil::deserialize(in, periodicFaces);
+
+
+
+ SerUtil::deserialize(in, nSize);
+ periodicDofAssoc.clear();
+ for (int i = 0; i < nSize; i++) {
+ DegreeOfFreedom dof;
+ std::set<DegreeOfFreedom> dofs;
+ SerUtil::deserialize(in, dof);
+ SerUtil::deserialize(in, dofs);
+ periodicDofAssoc[dof] = dofs;
+ }
+
+ SerUtil::deserialize(in, nSize);
+ periodicEdgeAssoc.clear();
+ for (int i = 0; i < nSize; i++) {
+ DofEdge edge;
+ std::set<DofEdge> edges;
+ SerUtil::deserialize(in, edge);
+ SerUtil::deserialize(in, edges);
+ periodicEdgeAssoc[edge] = edges;
+ }
+
+
+
+ SerUtil::deserialize(in, nSize);
+ edgeReverseMode.clear();
+ for (int i = 0; i < nSize; i++) {
+ ElementObjectData obj0, obj1;
+ bool reverseMode;
+ obj0.deserialize(in);
+ obj1.deserialize(in);
+ SerUtil::deserialize(in, reverseMode);
+
+ edgeReverseMode[make_pair(obj0, obj1)] = reverseMode;
+ }
+
+ SerUtil::deserialize(in, nSize);
+ faceReverseMode.clear();
+ for (int i = 0; i < nSize; i++) {
+ ElementObjectData obj0, obj1;
+ bool reverseMode;
+ obj0.deserialize(in);
+ obj1.deserialize(in);
+ SerUtil::deserialize(in, reverseMode);
+
+ faceReverseMode[make_pair(obj0, obj1)] = reverseMode;
+ }
}
diff --git a/AMDiS/src/parallel/ElementObjectData.h b/AMDiS/src/parallel/ElementObjectData.h
index d859595dd28f509af18b7b8691682cdccff770ef..2bf63bd2c48e82a724955b7a62d70a3b76dba43c 100644
--- a/AMDiS/src/parallel/ElementObjectData.h
+++ b/AMDiS/src/parallel/ElementObjectData.h
@@ -148,6 +148,9 @@ namespace AMDiS {
*/
void createRankData(map<int, int>& macroElementRankMap);
+ void createReverseModeData(map<int, Element*> &elIndexMap,
+ map<int, int> &elIndexTypeMap);
+
/** \brief
* Iterates over all elements for one geometrical index, i.e., over all vertices,
@@ -320,6 +323,32 @@ namespace AMDiS {
return faceElements[face];
}
+
+ /// Returns a vector with all macro elements that have a given vertex DOF in common.
+ vector<ElementObjectData>& getElementsVertex(int elIndex, int ithVertex)
+ {
+ ElementObjectData elObj(elIndex, ithVertex);
+ DegreeOfFreedom vertex = vertexLocalMap[elObj];
+ return vertexElements[vertex];
+ }
+
+ /// Returns a vector with all macro elements that have a given edge in common.
+ vector<ElementObjectData>& getElementsEdge(int elIndex, int ithEdge)
+ {
+ ElementObjectData elObj(elIndex, ithEdge);
+ DofEdge edge = edgeLocalMap[elObj];
+ return edgeElements[edge];
+ }
+
+ /// Returns a vector with all macro elements that have a given face in common.
+ vector<ElementObjectData>& getElementsFace(int elIndex, int ithFace)
+ {
+ ElementObjectData elObj(elIndex, ithFace);
+ DofFace face = faceLocalMap[elObj];
+ return faceElements[face];
+ }
+
+
/// Returns a map that maps to each rank all macro elements in this rank that
/// have a given vertex DOF in common.
@@ -381,6 +410,22 @@ namespace AMDiS {
return periodicFaces;
}
+ bool getEdgeReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+ {
+ TEST_EXIT_DBG(edgeReverseMode.count(make_pair(obj0, obj1)))
+ ("Should not happen!\n");
+
+ return edgeReverseMode[make_pair(obj0, obj1)];
+ }
+
+ bool getFaceReverseMode(ElementObjectData &obj0, ElementObjectData &obj1)
+ {
+ TEST_EXIT_DBG(faceReverseMode.count(make_pair(obj0, obj1)))
+ ("Should not happen!\n");
+
+ return faceReverseMode[make_pair(obj0, obj1)];
+ }
+
/// Write the element database to disk.
void serialize(ostream &out);
@@ -505,7 +550,7 @@ namespace AMDiS {
/// iterators, i.e., if no iteration is running.
GeoIndex iterGeoPos;
- std::map<std::pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap;
+ map<pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap;
// The following three data structures store periodic DOFs, edges and faces.
PerBoundMap<DegreeOfFreedom>::type periodicVertices;
@@ -513,10 +558,14 @@ namespace AMDiS {
PerBoundMap<DofFace>::type periodicFaces;
// Stores to each vertex all its periodic associations.
- std::map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
+ map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc;
// Stores to each edge all its periodic associations.
- std::map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc;
+ map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc;
+
+ map<pair<ElementObjectData, ElementObjectData>, bool> edgeReverseMode;
+
+ map<pair<ElementObjectData, ElementObjectData>, bool> faceReverseMode;
};
}
diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc
index 4c8e92469f7f88d62f5e3e3ffba86642d6f47ef2..9b4366245b3c8f0dcc02b0b83966bdb1b01bfa9f 100644
--- a/AMDiS/src/parallel/InteriorBoundary.cc
+++ b/AMDiS/src/parallel/InteriorBoundary.cc
@@ -18,35 +18,38 @@
namespace AMDiS {
- void BoundaryObject::setReverseMode(BoundaryObject &otherBound,
- FiniteElemSpace *feSpace,
- BoundaryType boundary)
+
+ bool BoundaryObject::computeReverseMode(BoundaryObject &obj0,
+ BoundaryObject &obj1,
+ FiniteElemSpace *feSpace,
+ BoundaryType boundary)
{
- FUNCNAME("BoundaryObject::setReverseMode()");
+ FUNCNAME("BoundaryObject::computeReverseMode()");
- bool otherMode = false;
+ bool reverseMode = false;
switch (feSpace->getMesh()->getDim()) {
case 2:
- otherMode = true;
+ reverseMode = true;
break;
case 3:
- TEST_EXIT_DBG(otherBound.elType == 0)
- ("Only 3D macro elements with level 0 are supported. This element has level %d!\n", otherBound.elType);
+ TEST_EXIT_DBG(obj1.elType == 0)
+ ("Only 3D macro elements with level 0 are supported. This element has level %d!\n",
+ obj1.elType);
- if (subObj == EDGE) {
- int el0_v0 = el->getVertexOfEdge(ithObj, 0);
- int el0_v1 = el->getVertexOfEdge(ithObj, 1);
- int el1_v0 = el->getVertexOfEdge(otherBound.ithObj, 0);
- int el1_v1 = el->getVertexOfEdge(otherBound.ithObj, 1);
+ if (obj0.subObj == EDGE) {
+ int el0_v0 = obj0.el->getVertexOfEdge(obj0.ithObj, 0);
+ int el0_v1 = obj0.el->getVertexOfEdge(obj0.ithObj, 1);
+ int el1_v0 = obj0.el->getVertexOfEdge(obj1.ithObj, 0);
+ int el1_v1 = obj0.el->getVertexOfEdge(obj1.ithObj, 1);
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> localDofs0(nBasFcts), localDofs1(nBasFcts);
- basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0);
- basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1);
+ basFcts->getLocalIndices(obj0.el, feSpace->getAdmin(), localDofs0);
+ basFcts->getLocalIndices(obj1.el, feSpace->getAdmin(), localDofs1);
Mesh *mesh = feSpace->getMesh();
@@ -59,30 +62,25 @@ namespace AMDiS {
("This should not happen!\n");
if (localDofs0[el0_v0] != localDofs1[el1_v0])
- otherMode = true;
+ reverseMode = true;
} else {
if (mesh->associated(localDofs0[el0_v0], localDofs1[el1_v0]) == false)
- otherMode = true;
-
- if ((elIndex == 28 && otherBound.elIndex == 41) ||
- (elIndex == 41 && otherBound.elIndex == 28)) {
- MSG("HERE TEST B (%d %d): %d\n", el0_v0, el1_v0, otherMode);
- }
+ reverseMode = true;
}
}
- if (subObj == FACE && ithObj != 1) {
+ if (obj0.subObj == FACE && obj0.ithObj != 1) {
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> localDofs0(nBasFcts), localDofs1(nBasFcts);
- basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0);
- basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1);
+ basFcts->getLocalIndices(obj0.el, feSpace->getAdmin(), localDofs0);
+ basFcts->getLocalIndices(obj1.el, feSpace->getAdmin(), localDofs1);
- if (ithObj == 2 || ithObj == 3)
- otherMode = (localDofs0[0] != localDofs1[0]);
+ if (obj0.ithObj == 2 || obj0.ithObj == 3)
+ reverseMode = (localDofs0[0] != localDofs1[0]);
- if (ithObj == 0)
- otherMode = (localDofs0[1] != localDofs1[1]);
+ if (obj0.ithObj == 0)
+ reverseMode = (localDofs0[1] != localDofs1[1]);
}
break;
@@ -90,14 +88,23 @@ namespace AMDiS {
ERROR_EXIT("This should not happen!\n");
}
- otherBound.reverseMode = otherMode;
+ return reverseMode;
+ }
+
+
+ void BoundaryObject::setReverseMode(BoundaryObject &otherBound,
+ FiniteElemSpace *feSpace,
+ BoundaryType boundary)
+ {
+ FUNCNAME("BoundaryObject::setReverseMode()");
+
+ otherBound.reverseMode = computeReverseMode(*this, otherBound, feSpace, boundary);
}
bool BoundaryObject::operator==(const BoundaryObject& other) const
{
return (other.elIndex == elIndex &&
- // other.elType == elType &&
other.subObj == subObj &&
other.ithObj == ithObj);
}
@@ -106,7 +113,6 @@ namespace AMDiS {
bool BoundaryObject::operator!=(const BoundaryObject& other) const
{
return (other.elIndex != elIndex ||
- // other.elType != elType ||
other.subObj != subObj ||
other.ithObj != ithObj);
}
diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h
index 81280de026a99c88033b6498747b0f20fa2f06fc..b68bff77c1722bb157827059fc7587a15aafbf78 100644
--- a/AMDiS/src/parallel/InteriorBoundary.h
+++ b/AMDiS/src/parallel/InteriorBoundary.h
@@ -54,6 +54,11 @@ namespace AMDiS {
excludedSubstructures(0)
{}
+ static bool computeReverseMode(BoundaryObject &obj0,
+ BoundaryObject &obj1,
+ FiniteElemSpace *feSpace,
+ BoundaryType boundary);
+
void setReverseMode(BoundaryObject &otherBound,
FiniteElemSpace *feSpace,
BoundaryType boundary);
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index f1d9fe6166bcb3181986100ca622c6c1db8ff555..9c64ffde0e9df42fffa8dfaf838a82b42c17b860 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -92,6 +92,8 @@ namespace AMDiS {
tmp = 0;
GET_PARAMETER(0, name + "->log main rank", "%d", &tmp);
Msg::outputMainRank = (tmp > 0);
+
+ partitioner = new ParMetisPartitioner(&mpiComm);
}
@@ -120,18 +122,12 @@ namespace AMDiS {
macroElIndexMap.clear();
macroElIndexTypeMap.clear();
- map<int, bool>& elementInRank = partitioner->getElementInRank();
for (vector<MacroElement*>::iterator it = allMacroElements.begin();
it != allMacroElements.end(); ++it) {
- elementInRank[(*it)->getIndex()] = false;
macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
}
- for (deque<MacroElement*>::iterator it = mesh->getMacroElements().begin();
- it != mesh->getMacroElements().end(); ++it)
- elementInRank[(*it)->getIndex()] = true;
-
return;
}
@@ -146,7 +142,7 @@ namespace AMDiS {
createMacroElementInfo();
// create an initial partitioning of the mesh
- partitioner->createPartitionData();
+ partitioner->createInitialPartitioning();
// set the element weights, which are 1 at the very first begin
setInitialElementWeights();
@@ -316,7 +312,7 @@ namespace AMDiS {
ERROR_EXIT("This should not happen for dim = %d!\n", mesh->getDim());
}
- partitioner = new ParMetisPartitioner(mesh, &mpiComm);
+ partitioner->setMesh(mesh);
}
// Create parallel serialization file writer, if needed.
@@ -699,9 +695,7 @@ namespace AMDiS {
updateLocalGlobalNumbering();
-
// === Update periodic mapping, if there are periodic boundaries. ===
-
createPeriodicMap();
INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n",
@@ -981,6 +975,8 @@ namespace AMDiS {
for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
it != newMacroEl.end(); ++it) {
+ MSG("NEW MACRO EL: %d\n", (*it)->getIndex());
+
MacroElement *mel = *it;
// First, reset all neighbour relations. The correct neighbours will be set later.
@@ -1100,6 +1096,8 @@ namespace AMDiS {
("Could not find macro element %d\n", *elIt);
deleteMacroElements.insert(elIndexMap[*elIt]);
+
+ MSG("REMOVE MACRO EL: %d\n", *elIt);
}
}
@@ -1114,12 +1112,93 @@ namespace AMDiS {
MeshManipulation meshManipulation(feSpace);
meshManipulation.deleteDoubleDofs(newMacroEl, elObjects);
- mesh->dofCompress();
+#if 0
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
+ while (elInfo) {
+ MSG("HAVE MACRO EL %d\n", elInfo->getMacroElement()->getIndex());
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ MSG("A-CHECK DOFs ON %d: %d %d %d %d\n",
+ elInfo->getElement()->getIndex(),
+ elInfo->getElement()->getDof(0, 0),
+ elInfo->getElement()->getDof(1, 0),
+ elInfo->getElement()->getDof(2, 0),
+ elInfo->getElement()->getDof(3, 0));
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+#endif
+ mesh->dofCompress();
partitioner->fillCoarsePartitionVec(&partitionVec);
+
+#if 0
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ MSG("B-CHECK DOFs ON %d: %d %d %d %d\n",
+ elInfo->getElement()->getIndex(),
+ elInfo->getElement()->getDof(0, 0),
+ elInfo->getElement()->getDof(1, 0),
+ elInfo->getElement()->getDof(2, 0),
+ elInfo->getElement()->getDof(3, 0));
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+#endif
+
updateInteriorBoundaryInfo();
+
+#if 0
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ MSG("C-CHECK DOFs ON %d: %d %d %d %d\n",
+ elInfo->getElement()->getIndex(),
+ elInfo->getElement()->getDof(0, 0),
+ elInfo->getElement()->getDof(1, 0),
+ elInfo->getElement()->getDof(2, 0),
+ elInfo->getElement()->getDof(3, 0));
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+#endif
+
updateLocalGlobalNumbering();
+#if 0
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ MSG("D-CHECK DOFs ON %d: %d %d %d %d\n",
+ elInfo->getElement()->getIndex(),
+ elInfo->getElement()->getDof(0, 0),
+ elInfo->getElement()->getDof(1, 0),
+ elInfo->getElement()->getDof(2, 0),
+ elInfo->getElement()->getDof(3, 0));
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+#endif
+
+
+
#if (DEBUG != 0)
MSG("AMDiS runs in debug mode, so make some test ...\n");
@@ -1183,6 +1262,9 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
+
+ elObjects.createReverseModeData(macroElIndexMap, macroElIndexTypeMap);
+
// === Create periodic data, if there are periodic boundary conditions. ===
elObjects.createPeriodicData();
@@ -1196,7 +1278,6 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::createBoundaryData()");
-
// === Clear all relevant data structures. ===
myIntBoundary.clear();
@@ -1250,7 +1331,11 @@ namespace AMDiS {
AtomicBoundary& b = myIntBoundary.getNewAtomic(it2->first);
b = bound;
- b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
+ // b.rankObj.setReverseMode(b.neighObj, feSpace, bound.type);
+ if (geoIndex == EDGE)
+ b.neighObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, it2->second);
+ if (geoIndex == FACE)
+ b.neighObj.reverseMode = elObjects.getFaceReverseMode(rankBoundEl, it2->second);
}
} else {
@@ -1269,13 +1354,19 @@ namespace AMDiS {
AtomicBoundary& b = otherIntBoundary.getNewAtomic(owner);
b = bound;
- b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
+
+ // b.neighObj.setReverseMode(b.rankObj, feSpace, bound.type);
+
+ if (geoIndex == EDGE)
+ b.rankObj.reverseMode = elObjects.getEdgeReverseMode(rankBoundEl, ownerBoundEl);
+ if (geoIndex == FACE)
+ b.rankObj.reverseMode = elObjects.getFaceReverseMode(rankBoundEl, ownerBoundEl);
}
}
}
}
-
+
// === Create periodic boundary data structure. ===
for (PerBoundMap<DegreeOfFreedom>::iterator it = elObjects.getPeriodicVertices().begin();
@@ -1902,6 +1993,8 @@ namespace AMDiS {
checkMeshChange();
+ partitioner->serialize(out);
+
SerUtil::serialize(out, elemWeights);
SerUtil::serialize(out, partitionVec);
SerUtil::serialize(out, oldPartitionVec);
@@ -1941,6 +2034,8 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::deserialize()");
+ partitioner->deserialize(in);
+
SerUtil::deserialize(in, elemWeights);
SerUtil::deserialize(in, partitionVec);
SerUtil::deserialize(in, oldPartitionVec);
diff --git a/AMDiS/src/parallel/MeshManipulation.cc b/AMDiS/src/parallel/MeshManipulation.cc
index e72c4cd96d0b4eaf986dde5f7ed5679ce1bb0139..e17f001194de3a3b01f6d01e5ea24d34f66b4934 100644
--- a/AMDiS/src/parallel/MeshManipulation.cc
+++ b/AMDiS/src/parallel/MeshManipulation.cc
@@ -47,22 +47,19 @@ namespace AMDiS {
{
FUNCNAME("MeshManipulation::deleteDoubleDofs()");
- TEST_EXIT(mesh->getDim() == 2)("Not yet supported for dim != 2!\n");
-
-
// Create data structure that maps macro element indices to the
// corresponding pointers.
- std::map<int, MacroElement*> macroIndexMap;
+ map<int, MacroElement*> macroIndexMap;
for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
it != newMacroEl.end(); ++it)
macroIndexMap[(*it)->getIndex()] = *it;
std::set<int> macrosProcessed;
- std::map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
+ map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;
// === Traverse mesh and put all "old" macro element to macrosProcessed ===
- // === that stores all macro elements which are really conected to the ===
+ // === that stores all macro elements which are really connected to the ===
// === overall mesh structure. ===
TraverseStack stack;
@@ -86,12 +83,13 @@ namespace AMDiS {
for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
it != newMacroEl.end(); ++it) {
+ // === Traverse all vertices of the new element. ===
+
for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
- ElementObjectData elObj((*it)->getIndex(), i);
- DegreeOfFreedom vertex = objects.getVertexLocalMap(elObj);
- std::vector<ElementObjectData> &vertexEl = objects.getElements(vertex);
+ vector<ElementObjectData> &vertexEl =
+ objects.getElementsVertex((*it)->getIndex(), i);
- for (std::vector<ElementObjectData>::iterator elIt = vertexEl.begin();
+ for (vector<ElementObjectData>::iterator elIt = vertexEl.begin();
elIt != vertexEl.end(); ++elIt) {
if (elIt->elIndex == (*it)->getIndex())
continue;
@@ -115,33 +113,32 @@ namespace AMDiS {
for (int i = 0; i < mesh->getGeo(EDGE); i++) {
ElementObjectData elObj((*it)->getIndex(), i);
- DofEdge edge = objects.getEdgeLocalMap(elObj);
- std::vector<ElementObjectData> &edgeEl = objects.getElements(edge);
- for (std::vector<ElementObjectData>::iterator elIt = edgeEl.begin();
+ vector<ElementObjectData> &edgeEl =
+ objects.getElementsEdge((*it)->getIndex(), i);
+
+ for (vector<ElementObjectData>::iterator elIt = edgeEl.begin();
elIt != edgeEl.end(); ++elIt) {
if (elIt->elIndex == (*it)->getIndex())
continue;
if (macrosProcessed.count(elIt->elIndex) == 1) {
- TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex) == 1)
- ("Should not happen!\n");
-
- TEST_EXIT(mesh->getDim() == 2)("In 3D set correct reverse mode!\n");
-
- TEST_EXIT(feSpace->getBasisFcts()->getDegree() == 1)
- ("Not yet implemented!\n");
+ TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))("Should not happen!\n");
Element *el0 = (*it)->getElement();
Element *el1 = macroIndexMap[elIt->elIndex]->getElement();
- BoundaryObject b0(el0, 0, EDGE, i, true);
+ bool reverseMode = objects.getEdgeReverseMode(elObj, *elIt);
+
+ BoundaryObject b0(el0, 0, EDGE, i, reverseMode);
BoundaryObject b1(el1, 0, EDGE, elIt->ithObject, false);
DofContainer dofs0, dofs1;
el0->getVertexDofs(feSpace, b0, dofs0);
+ el0->getNonVertexDofs(feSpace, b0, dofs0);
el1->getVertexDofs(feSpace, b1, dofs1);
+ el1->getNonVertexDofs(feSpace, b1, dofs1);
#if (DEBUG != 0)
debug::testDofsByCoords(feSpace, dofs0, dofs1);
@@ -155,8 +152,47 @@ namespace AMDiS {
}
}
- TEST_EXIT(mesh->getDim() == 2)("Add face traverse here!\n");
+ for (int i = 0; i < mesh->getGeo(FACE); i++) {
+ ElementObjectData elObj((*it)->getIndex(), i);
+
+ vector<ElementObjectData> &faceEl =
+ objects.getElementsFace((*it)->getIndex(), i);
+
+ for (vector<ElementObjectData>::iterator elIt = faceEl.begin();
+ elIt != faceEl.end(); ++elIt) {
+ if (elIt->elIndex == (*it)->getIndex())
+ continue;
+
+ if (macrosProcessed.count(elIt->elIndex) == 1) {
+ TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))("Should not happen!\n");
+
+ Element *el0 = (*it)->getElement();
+ Element *el1 = macroIndexMap[elIt->elIndex]->getElement();
+
+ bool reverseMode = objects.getFaceReverseMode(elObj, *elIt);
+
+ BoundaryObject b0(el0, 0, FACE, i, reverseMode);
+ BoundaryObject b1(el1, 0, FACE, elIt->ithObject, false);
+
+ DofContainer dofs0, dofs1;
+
+ el0->getVertexDofs(feSpace, b0, dofs0);
+ el0->getNonVertexDofs(feSpace, b0, dofs0);
+ el1->getVertexDofs(feSpace, b1, dofs1);
+ el1->getNonVertexDofs(feSpace, b1, dofs1);
+
+#if (DEBUG != 0)
+ debug::testDofsByCoords(feSpace, dofs0, dofs1);
+#endif
+
+ for (unsigned int i = 0; i < dofs0.size(); i++)
+ mapDelDofs[dofs0[i]] = dofs1[i];
+
+ break;
+ }
+ }
+ }
macrosProcessed.insert((*it)->getIndex());
}
@@ -167,9 +203,9 @@ namespace AMDiS {
bool changed = false;
do {
changed = false;
- for (std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
+ for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
it != mapDelDofs.end(); ++it) {
- std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second);
+ map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second);
if (findIt != mapDelDofs.end()) {
TEST_EXIT_DBG(it->first != findIt->second)
("Cycle %d -> %d and %d -> %d! Should not happen!\n",
@@ -195,7 +231,7 @@ namespace AMDiS {
// === And delete all double DOFs. ===
- for (std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
+ for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
it != mapDelDofs.end(); ++it)
mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), VERTEX);
}
diff --git a/AMDiS/src/parallel/ParMetisPartitioner.cc b/AMDiS/src/parallel/ParMetisPartitioner.cc
index 0c0181e672b461c1a04e1653ae6f05aae522fd63..f8d4dfd36389fe87520b9be0ad8ca86c03a9e8dd 100644
--- a/AMDiS/src/parallel/ParMetisPartitioner.cc
+++ b/AMDiS/src/parallel/ParMetisPartitioner.cc
@@ -13,6 +13,7 @@
#include <queue>
#include "parallel/ParMetisPartitioner.h"
#include "parallel/MpiHelper.h"
+#include "Serializer.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
@@ -221,9 +222,9 @@ namespace AMDiS {
}
- void ParMetisPartitioner::createPartitionData()
+ void ParMetisPartitioner::createInitialPartitioning()
{
- FUNCNAME("ParMetrisPartitioner::createPartitionData()");
+ FUNCNAME("ParMetrisPartitioner::createInitialPartitioning()");
int mpiRank = mpiComm->Get_rank();
int mpiSize = mpiComm->Get_size();
@@ -660,4 +661,15 @@ namespace AMDiS {
return true;
}
+
+ void ParMetisPartitioner::serialize(std::ostream &out)
+ {
+ SerUtil::serialize(out, elementInRank);
+ }
+
+
+ void ParMetisPartitioner::deserialize(std::istream &in)
+ {
+ SerUtil::deserialize(in, elementInRank);
+ }
}
diff --git a/AMDiS/src/parallel/ParMetisPartitioner.h b/AMDiS/src/parallel/ParMetisPartitioner.h
index 11c713d4307a2eed9f158c4280299f72395a8547..a96974fd2090710667fdc17487923af86c7806bd 100644
--- a/AMDiS/src/parallel/ParMetisPartitioner.h
+++ b/AMDiS/src/parallel/ParMetisPartitioner.h
@@ -172,15 +172,20 @@ namespace AMDiS {
class ParMetisPartitioner
{
public:
- ParMetisPartitioner(Mesh *m, MPI::Intracomm *comm)
- : mesh(m),
- mpiComm(comm),
+ ParMetisPartitioner(MPI::Intracomm *comm)
+ : mpiComm(comm),
+ mesh(NULL),
parMetisMesh(NULL),
mapLocalGlobal(NULL)
{}
~ParMetisPartitioner();
+ void setMesh(Mesh *m)
+ {
+ mesh = m;
+ }
+
bool partition(std::map<int, double> &elemWeights,
PartitionMode mode = INITIAL,
float itr = 1000000.0);
@@ -192,7 +197,7 @@ namespace AMDiS {
void fillCoarsePartitionVec(std::map<int, int> *partitionVec);
/// Creates an initial paritioning of the AMDiS mesh.
- void createPartitionData();
+ void createInitialPartitioning();
void useLocalGlobalDofMap(std::map<DegreeOfFreedom, DegreeOfFreedom> *m)
{
@@ -214,20 +219,27 @@ namespace AMDiS {
return sendElements;
}
+ /// Write partitioner state to disk.
+ void serialize(std::ostream &out);
+
+ /// Read partitioner state from disk.
+ void deserialize(std::istream &in);
+
protected:
// Returns true, if the mesh partitioning could be distributed. If this is not the
// case, i.e., because there are empty partitions, the function returns false.
bool distributePartitioning(int *part);
protected:
- Mesh *mesh;
-
MPI::Intracomm *mpiComm;
+ Mesh *mesh;
+
ParMetisMesh *parMetisMesh;
std::map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal;
+ /// Maps to each macro element index if it is in rank's partition or not.
std::map<int, bool> elementInRank;
std::map<int, std::vector<int> > recvElements, sendElements;