diff --git a/AMDiS/libtool b/AMDiS/libtool
index cce1135148bfc12b7832dc2f0ba3eb5c551d83d2..7d1660395dd36a2f3331e850fc6f2f4255dd82de 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host deimos101:
+# Libtool was configured on host deimos103:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7263,7 +7263,7 @@ disable_libs=static
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host deimos101:
+# Libtool was configured on host deimos103:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7568,7 +7568,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host deimos101:
+# Libtool was configured on host deimos103:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
diff --git a/AMDiS/other/include/Makefile_AMDiS.mk b/AMDiS/other/include/Makefile_AMDiS.mk
index 11ffa5997d150a32d525c1853b97d9d476f0fead..6931fd50a09c55ebaac0d4dc50f0843ef2fb5723 100644
--- a/AMDiS/other/include/Makefile_AMDiS.mk
+++ b/AMDiS/other/include/Makefile_AMDiS.mk
@@ -39,7 +39,7 @@ MPCCI_LIB = -L$(MPCCI_DIR)/lib/linux-x86-glibc22 -lmpcci
PARMETIS_LIB = -L$(PARMETIS_DIR) -lparmetis -lmetis
LIBS = $(AMDIS_LIB) $(PNG_LIB)
-LIBS += -lboost_iostreams -lboost_filesystem -lboost_system
+LIBS += -lboost_iostreams -lboost_filesystem -lboost_system -lboost_date_time
ifeq ($(strip $(USE_UMFPACK)), 1)
LIBS += $(UMFPACK_LIB)
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index bdb97fd8721ad2fc8246f82c7daccac671c2f13e..25cceb930c026e017b7abba91c0515f306d08ba3 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -1,3 +1,5 @@
+#include <typeinfo>
+
#include "ElInfo3d.h"
#include "BasisFunction.h"
#include "Element.h"
@@ -478,6 +480,8 @@ namespace AMDiS {
{
FUNCNAME("ElInfo3d::fillElInfo()");
+ TEST_EXIT_DBG(elInfoOld)("Missing old elInfo!\n");
+
int ochild = 0; /* index of other child = 1-ichild */
int *cv = NULL; /* cv = child_vertex[el_type][ichild] */
const int (*cvg)[4] = NULL; /* cvg = child_vertex[el_type] */
@@ -497,7 +501,13 @@ namespace AMDiS {
parent = el_old;
level = elInfoOld->level + 1;
iChild = ichild;
- int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
+ int el_type_local = 0;
+ try {
+ el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
+ } catch (const std::bad_cast& e) {
+ ERROR_EXIT("ElInfo is not of type ElInfo3D but %s!\n", typeid(*elInfoOld).name());
+ }
+
elType = (el_type_local + 1) % 3;
TEST_EXIT_DBG(element)("missing child %d?\n", ichild);
diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc
index b8aa8abd130cdd45078ce4f13440d6208c33d468..bff5e1f2aca3026e72bcf0a140a24860aa920539 100644
--- a/AMDiS/src/Element.cc
+++ b/AMDiS/src/Element.cc
@@ -11,6 +11,7 @@ namespace AMDiS {
std::map<DegreeOfFreedom*, bool> Element::deletedDOFs;
+
int Element::getRegion() const
{
if (!elementData)
@@ -25,6 +26,7 @@ namespace AMDiS {
return -1;
}
+
void Element::setDOFPtrs()
{
FUNCNAME("Element::setDOFPtrs()");
@@ -34,6 +36,7 @@ namespace AMDiS {
dof = mesh->createDOFPtrs();
}
+
Element::Element(Mesh *aMesh)
{
mesh = aMesh;
@@ -51,6 +54,7 @@ namespace AMDiS {
}
}
+
// call destructor through Mesh::freeElement !!!
Element::~Element()
{
@@ -69,6 +73,7 @@ namespace AMDiS {
}
}
+
bool Element::deleteElementData(int typeID)
{
FUNCNAME("Element::deleteElementData()");
@@ -87,6 +92,7 @@ namespace AMDiS {
return false;
}
+
void Element::deleteElementDOFs()
{
int dim = mesh->getDim();
@@ -310,6 +316,7 @@ namespace AMDiS {
#undef CHANGE_DOFS_1
#undef CHANGE_DOFS_2
+
/****************************************************************************/
/* opp_vertex checks whether the face with vertices dof[0],..,dof[DIM-1] is */
/* part of mel's boundary. returns the opposite vertex if true, -1 else */
@@ -362,13 +369,16 @@ namespace AMDiS {
}
}
- void Element::eraseNewCoord() {
+
+ void Element::eraseNewCoord()
+ {
if (newCoord != NULL) {
delete newCoord;
newCoord = NULL;
}
}
+
void Element::serialize(std::ostream &out)
{
// write children
@@ -445,6 +455,7 @@ namespace AMDiS {
}
}
+
void Element::deserialize(std::istream &in)
{
FUNCNAME("Element::deserialize()");
@@ -555,6 +566,7 @@ namespace AMDiS {
}
}
+
int Element::calcMemoryUsage()
{
int result = 0;
@@ -568,4 +580,55 @@ namespace AMDiS {
return result;
}
+
+ void writeDotFile(Element *el, std::string filename, int maxLevels)
+ {
+ std::vector<int> graphNodes;
+ std::vector<std::pair<int, int> > graphEdges;
+
+ std::vector<Element*> traverseElements, nextTraverseElements;
+ traverseElements.push_back(el);
+
+ int nLevels = 0;
+
+ while (!traverseElements.empty() && (maxLevels == -1 || nLevels < maxLevels)) {
+ nextTraverseElements.clear();
+
+ for (unsigned int i = 0; i < traverseElements.size(); i++) {
+ graphNodes.push_back(traverseElements[i]->getIndex());
+
+ if (!traverseElements[i]->isLeaf() &&
+ (maxLevels == -1 || nLevels + 1 < maxLevels)) {
+ graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
+ traverseElements[i]->getChild(0)->getIndex()));
+ graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
+ traverseElements[i]->getChild(1)->getIndex()));
+ nextTraverseElements.push_back(traverseElements[i]->getChild(0));
+ nextTraverseElements.push_back(traverseElements[i]->getChild(1));
+ }
+ }
+
+ traverseElements = nextTraverseElements;
+ nLevels++;
+ }
+
+ std::ofstream file;
+ file.open(filename.c_str());
+
+ file << "digraph G\n";
+ file << "{\n";
+
+ for (unsigned int i = 0; i < graphNodes.size(); i++)
+ file << " node" << graphNodes[i]
+ << " [ label = \"" << graphNodes[i] << "\"];\n";
+
+ for (unsigned int i = 0; i < graphEdges.size(); i++)
+ file << " \"node" << graphEdges[i].first << "\" -> \"node"
+ << graphEdges[i].second << "\";\n";
+
+ file << "}\n";
+
+ file.close();
+ }
+
}
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index 2d585f545b9d3fb3f25e9b1b28d8cb315a48c5e8..9a733d1651b708e48c820da1d766f0cfa2da2557 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -583,6 +583,11 @@ namespace AMDiS {
friend class Mesh;
};
+
+
+ /// Writes the element hierarchie to a Graphviz dot-file. Using the dot-tool from
+ /// Graphvis, this dot-file can be converted to a ps-file. Useful for debugging!
+ void writeDotFile(Element *el, std::string filename, int maxLevels = -1);
}
#endif // AMDIS_ELEMENT_H
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index 5c33274875d4fcd446881370573c51e8c6114932..8d8e7c8d0d37f18d2fcf9bfb42a4c067589d06c1 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -504,7 +504,7 @@ namespace AMDiS {
/// returns the euclidian distance of a and b
template<typename T, GeoIndex d>
- double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b);
+ double absteukl(const FixVec<T,d>& a, const FixVec<T,d>& b);
/// FixVec operator for stream output
template<typename T, GeoIndex d>
@@ -587,6 +587,15 @@ namespace AMDiS {
template<typename T>
const WorldMatrix<T>& operator+=(WorldMatrix<T>& m1, const WorldMatrix<T>& m2);
+
+
+ inline double norm(const WorldVector<double>& v)
+ {
+ double val = 0.0;
+ for (int i = 0; i < Global::getGeo(WORLD); i++)
+ val += v[i] * v[i];
+ return sqrt(val);
+ }
}
#include "FixVec.hh"
diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc
index 859fb2d4c6855c63727580627a09bf9528ff0ed0..472364483a8548e206b87ba73b6921e996a964c2 100644
--- a/AMDiS/src/MeshStructure.cc
+++ b/AMDiS/src/MeshStructure.cc
@@ -330,8 +330,16 @@ namespace AMDiS {
cont = nextElement();
}
}
-
- MSG("Mesh structure code: %s\n", oss.str().c_str());
+
+ if (oss.str().length() < 255) {
+ MSG("Mesh structure code: %s\n", oss.str().c_str());
+ } else {
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "] Mesh structure code: " << oss.str() << "\n";
+#else
+ std::cout << " Mesh structure code: " << oss.str() << "\n";
+#endif
+ }
}
diff --git a/AMDiS/src/parallel/InteriorBoundary.cc b/AMDiS/src/parallel/InteriorBoundary.cc
index b29ccfd9ba3ac15b0546ad6a32ac03e96a320c60..abf05375fb23e4da6741a542d2260899a98c7314 100644
--- a/AMDiS/src/parallel/InteriorBoundary.cc
+++ b/AMDiS/src/parallel/InteriorBoundary.cc
@@ -18,6 +18,9 @@ namespace AMDiS {
break;
case 3:
+ TEST_EXIT_DBG(otherBound.elType == 0)
+ ("Only 3D macro elements with level 0 are supported!\n");
+
if (subObj == EDGE) {
int el0_v0 = el->getVertexOfEdge(ithObj, 0);
int el0_v1 = el->getVertexOfEdge(ithObj, 1);
@@ -52,7 +55,7 @@ namespace AMDiS {
otherMode = (localDofs0[0] != localDofs1[0]);
if (ithObj == 0)
- otherMode = localDofs0[1] != localDofs1[1];
+ otherMode = (localDofs0[1] != localDofs1[1]);
}
break;
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 9413d263d8c7a22443ab8e1d60e44388c9450245..a63ac2fec6f3496f28b61e652e0da7b5acb9c96c 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -105,10 +105,12 @@ namespace AMDiS {
int writePartMesh = 1;
GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
- if (writePartMesh > 0)
+ if (writePartMesh > 0) {
+ debug::writeElementIndexMesh(mesh, "elementIndex.vtu");
writePartitioningMesh("part.vtu");
- else
+ } else {
MSG("Skip write part mesh!\n");
+ }
}
ParallelDebug::testAllElements(*this);
@@ -157,7 +159,6 @@ namespace AMDiS {
createPeriodicMap();
-
// === Global refinements. ===
int globalRefinement = 0;
@@ -172,7 +173,8 @@ namespace AMDiS {
mesh->dofCompress();
updateLocalGlobalNumbering();
-
+
+
// === Update periodic mapping, if there are periodic boundaries. ===
createPeriodicMap();
@@ -492,7 +494,7 @@ namespace AMDiS {
RankToBoundMap allBound;
for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
- if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+ if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
allBound[it.getRank()].push_back(*it);
for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
@@ -500,7 +502,7 @@ namespace AMDiS {
allBound[it.getRank()].push_back(*it);
for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
- if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+ if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
allBound[it.getRank()].push_back(*it);
@@ -579,7 +581,7 @@ namespace AMDiS {
MeshStructure elCode;
elCode.init(boundIt->rankObj);
-
+
if (elCode.getCode() != recvCodes[i].getCode()) {
TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
@@ -1083,7 +1085,7 @@ namespace AMDiS {
bound.neighObj.el = elIndexMap[it2->second.elIndex];
bound.neighObj.elIndex = it2->second.elIndex;
- bound.neighObj.elType = -1;
+ bound.neighObj.elType = elIndexTypeMap[it2->second.elIndex];
bound.neighObj.subObj = geoIndex;
bound.neighObj.ithObj = it2->second.ithObject;
@@ -1435,7 +1437,7 @@ namespace AMDiS {
it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
for (unsigned int i = 0; i < dofs.size(); i++)
- sendDofs[it.getRank()].push_back(dofs[i]);
+ sendDofs[it.getRank()].push_back(dofs[i]);
}
@@ -1448,7 +1450,7 @@ namespace AMDiS {
DofContainer::iterator eraseIt =
find(rankDofs.begin(), rankDofs.end(), dofs[i]);
if (eraseIt != rankDofs.end())
- rankDofs.erase(eraseIt);
+ rankDofs.erase(eraseIt);
recvDofs[it.getRank()].push_back(dofs[i]);
}
@@ -1796,7 +1798,7 @@ namespace AMDiS {
int *sendbuf = new int[2];
sendbuf[0] = periodicDof[type0][it->first];
sendbuf[1] = periodicDof[type1][it->first];
-
+
request[requestCounter++] =
mpiComm.Isend(sendbuf, 2, MPI_INT, *(it->second.begin()), 0);
request[requestCounter++] =
diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc
index 094428459d894ffc42fd1a0e4f3463b32fddf59f..d058d0252989ac77377ec5faf241a4568f86ecf9 100644
--- a/AMDiS/src/parallel/ParallelDebug.cc
+++ b/AMDiS/src/parallel/ParallelDebug.cc
@@ -139,13 +139,13 @@ namespace AMDiS {
typedef std::map<int, DofContainer> RankToDofContainer;
// Maps to each neighbour rank an array of WorldVectors. This array contains the
- // coordinates of all dofs this rank shares on the interior boundary with the
+ // coordinates of all DOFs this rank shares on the interior boundary with the
// neighbour rank. A rank sends the coordinates to another rank, if it owns the
- // boundarys dof.
+ // boundarys DOFs.
RankToCoords sendCoords;
- // A rank receives all boundary dofs that are at its interior boundaries but are
- // not owned by the rank. This map stores for each rank the coordinates of dofs
+ // A rank receives all boundary DOFs that are at its interior boundaries but are
+ // not owned by the rank. This map stores for each rank the coordinates of DOFs
// this rank expectes to receive from.
RankToCoords recvCoords;
@@ -180,14 +180,16 @@ namespace AMDiS {
if (i == pdb.mpiRank)
continue;
- request[requestCounter++] = pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
+ request[requestCounter++] =
+ pdb.mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
}
for (int i = 0; i < pdb.mpiSize; i++) {
if (i == pdb.mpiRank)
continue;
- request[requestCounter++] = pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
+ request[requestCounter++] =
+ pdb.mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
}
MPI::Request::Waitall(requestCounter, request);
@@ -217,26 +219,24 @@ namespace AMDiS {
stdMpi.recv(recvCoords);
stdMpi.startCommunication<double>(MPI_DOUBLE);
- double eps = 1e-13;
int dimOfWorld = Global::getGeo(WORLD);
// === Compare the received with the expected coordinates. ===
for (RankToCoords::iterator it = stdMpi.getRecvData().begin();
it != stdMpi.getRecvData().end(); ++it) {
- for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
- for (int j = 0; j < dimOfWorld; j++) {
- bool c = fabs((it->second)[i][j] - recvCoords[it->first][i][j]) < eps;
-
- // === Print error message if the coordinates are not the same. ===
-
- if (printCoords && !c) {
- MSG("[DBG] i = %d\n", i);
+ for (unsigned int i = 0; i < it->second.size(); i++) {
+ WorldVector<double> tmp = (it->second)[i];
+ tmp -= recvCoords[it->first][i];
+ if (norm(tmp) > 1e-13) {
+ // === Print error message if the coordinates are not the same. ===
+ if (printCoords) {
+ MSG("[DBG] i = %d\n", i);
std::stringstream oss;
oss.precision(5);
oss << "[DBG] Rank " << pdb.mpiRank << " from rank " << it->first
- << " expect coords (";
+ << " expect coords (";
for (int k = 0; k < dimOfWorld; k++) {
oss << recvCoords[it->first][i][k];
if (k + 1 < dimOfWorld)
@@ -250,16 +250,13 @@ namespace AMDiS {
}
oss << ")";
MSG("%s\n", oss.str().c_str());
-
+
debug::printInfoByDof(pdb.feSpace, *(pdb.recvDofs[it->first][i]));
}
-
- if (!c) {
- ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
- foundError = 1;
- }
- }
- }
+ ERROR("Wrong DOFs in rank %d!\n", pdb.mpiRank);
+ foundError = 1;
+ }
+ }
}
mpi::globalAdd(foundError);
TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
@@ -459,6 +456,14 @@ namespace AMDiS {
MSG(" neigh obj-ind: %d sub-obj: %d ith-obj: %d\n",
it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
}
+
+ for (InteriorBoundary::iterator it(pdb.periodicBoundary); !it.end(); ++it) {
+ MSG("Periodic boundary with rank %d: \n", it.getRank());
+ MSG(" ranks obj-ind: %d sub-obj: %d ith-obj: %d\n",
+ it->rankObj.elIndex, it->rankObj.subObj, it->rankObj.ithObj);
+ MSG(" neigh obj-ind: %d sub-obj: %d ith-obj: %d\n",
+ it->neighObj.elIndex, it->neighObj.subObj, it->neighObj.ithObj);
+ }
}