From 29fd364cfa6b6d9079372644319d96c2dacad67d Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 18 Jun 2010 09:32:08 +0000
Subject: [PATCH] Fixed a serious bug for mesh coarsening in parallel code.
---
AMDiS/libtool | 24 +--
AMDiS/src/BasisFunction.h | 3 +
AMDiS/src/DOFMatrix.cc | 5 +-
AMDiS/src/DOFVector.cc | 17 +-
AMDiS/src/Debug.cc | 89 +++++++-
AMDiS/src/Debug.h | 11 +
AMDiS/src/ElementFileWriter.cc | 21 +-
AMDiS/src/Lagrange.cc | 17 +-
AMDiS/src/Mesh.cc | 5 +-
AMDiS/src/Mesh.h | 3 +
AMDiS/src/MeshStructure.cc | 23 ++-
AMDiS/src/MeshStructure.h | 30 +--
AMDiS/src/ProblemInstat.cc | 10 +
AMDiS/src/ProblemVec.cc | 13 +-
AMDiS/src/RefinementManager.cc | 2 +-
AMDiS/src/Traverse.cc | 23 ++-
AMDiS/src/parallel/GlobalMatrixSolver.cc | 17 +-
AMDiS/src/parallel/MpiHelper.cc | 6 +
AMDiS/src/parallel/MpiHelper.h | 2 +
AMDiS/src/parallel/ParallelDomainBase.cc | 248 ++++++++++++++++-------
AMDiS/src/parallel/ParallelDomainBase.h | 14 +-
AMDiS/src/parallel/ParallelDomainDbg.cc | 43 +++-
AMDiS/src/parallel/ParallelDomainDbg.h | 3 +
23 files changed, 475 insertions(+), 154 deletions(-)
diff --git a/AMDiS/libtool b/AMDiS/libtool
index 076d8ba8..ae88f17c 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host p2q081:
+# Libtool was configured on host p2q084:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -82,13 +82,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="gcc"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -171,7 +171,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host p2q081:
+# Libtool was configured on host p2q084:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6798,13 +6798,13 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="g++"
+CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpiCC"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -6887,7 +6887,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6954,11 +6954,11 @@ predeps=""
# Dependencies to place after the objects being linked to create a
# shared library.
-postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
+postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -libverbs -lrt -lnuma -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path="-L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/fastfs/wir/local/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
+compiler_lib_search_path="-L/usr/lib64 -L/licsoft/libraries/openmpi/1.2.6/64bit/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/fastfs/wir/local/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host p2q081:
+# Libtool was configured on host p2q084:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7103,7 +7103,7 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="gcc"
+LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
# LTCC compiler flags.
LTCFLAGS="-g -O2"
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index d8038de6..e5612afb 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -41,6 +41,7 @@ namespace AMDiS {
virtual double operator()(const DimVec<double>&) const = 0;
};
+
/// Function interface for evaluating gradients of basis functions.
class GrdBasFctType
{
@@ -52,6 +53,7 @@ namespace AMDiS {
virtual void operator()(const DimVec<double>&, DimVec<double>&) const = 0;
};
+
/// Function interface for evaluating second derivative of basis functions.
class D2BasFctType
{
@@ -62,6 +64,7 @@ namespace AMDiS {
virtual void operator()(const DimVec<double>&, DimMat<double>&) const = 0;
};
+
typedef BasFctType *BFptr;
typedef GrdBasFctType *GBFptr;
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index b8227d97..64b339ff 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -219,10 +219,7 @@ namespace AMDiS {
} else {
for (int j = 0; j < nCol; j++) {
DegreeOfFreedom col = colIndices[j];
- double entry = elMat[i][j];
-
- // std::cout << "ADD at " << row << " " << col << std::endl;
- ins[row][col] += entry;
+ ins[row][col] += elMat[i][j];
}
}
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 593b8603..a1a6a2d9 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -9,6 +9,8 @@ namespace AMDiS {
template<>
void DOFVector<double>::coarseRestrict(RCNeighbourList& list, int n)
{
+ FUNCNAME("DOFVector<double>::coarseRestrict()");
+
switch (coarsenOperation) {
case NO_OPERATION:
return;
@@ -17,11 +19,8 @@ namespace AMDiS {
(const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
break;
case COARSE_INTERPOL:
- if (feSpace == NULL || !feSpace)
- ERROR_EXIT("ERR1\n");
-
- if (feSpace->getBasisFcts() == NULL || !(feSpace->getBasisFcts()))
- ERROR_EXIT("ERR2\n");
+ TEST_EXIT_DBG(feSpace)("Should not happen!\n");
+ TEST_EXIT_DBG(feSpace->getBasisFcts())("Shoud not happen!\n");
(const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
break;
@@ -31,12 +30,14 @@ namespace AMDiS {
}
}
+
template<>
void DOFVector<double>::refineInterpol(RCNeighbourList& list, int n)
{
(const_cast<BasisFunction*>(feSpace->getBasisFcts()))->refineInter(this, &list, n);
}
+
template<>
void DOFVector<WorldVector<double> >::refineInterpol(RCNeighbourList& list, int n)
{
@@ -53,6 +54,7 @@ namespace AMDiS {
vec[dof_new] *= 0.5;
}
+
template<>
DOFVector<WorldVector<double> >*
DOFVector<double>::getGradient(DOFVector<WorldVector<double> > *grad) const
@@ -138,6 +140,7 @@ namespace AMDiS {
return result;
}
+
template<>
DOFVector<WorldVector<double> >*
DOFVector<double>::getRecoveryGradient(DOFVector<WorldVector<double> > *grad) const
@@ -211,6 +214,7 @@ namespace AMDiS {
return result;
}
+
template<>
const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
@@ -296,6 +300,7 @@ namespace AMDiS {
return const_cast<const WorldVector<double>*>(result);
}
+
template<>
const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
@@ -373,6 +378,7 @@ namespace AMDiS {
return const_cast<const WorldVector<double>*>(result);
}
+
template<>
const WorldMatrix<double> *DOFVectorBase<double>::getD2AtQPs(const ElInfo *elInfo,
const Quadrature *quad,
@@ -469,6 +475,7 @@ namespace AMDiS {
return const_cast<const WorldMatrix<double>*>(result);
}
+
template<>
void DOFVector<double>::interpol(DOFVector<double> *source, double factor)
{
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
index 3727dfee..d7291ed4 100644
--- a/AMDiS/src/Debug.cc
+++ b/AMDiS/src/Debug.cc
@@ -6,6 +6,7 @@
#include "MacroElement.h"
#include "VtkWriter.h"
#include "ElementFileWriter.h"
+#include "ElementDofIterator.h"
namespace AMDiS {
@@ -165,6 +166,57 @@ namespace AMDiS {
return NULL;
}
+
+ Element* getParentElement(Mesh *mesh, Element *el)
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ if (elInfo->getElement()->getChild(0) == el ||
+ elInfo->getElement()->getChild(1) == el)
+ return elInfo->getElement();
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ return NULL;
+
+ }
+
+
+ Element* getParentElement(Mesh *mesh, int elIndex)
+ {
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ if (elInfo->getElement()->isLeaf() == false)
+ if (elInfo->getElement()->getChild(0)->getIndex() == elIndex ||
+ elInfo->getElement()->getChild(1)->getIndex() == elIndex)
+ return elInfo->getElement();
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ return NULL;
+
+ }
+
+
+ void printElementInfo(Element *el)
+ {
+ FUNCNAME("printElementInfo()");
+
+ if (el->isLeaf()) {
+ MSG("el %d is leaf!\n", el->getIndex());
+ } else {
+ MSG("el child0 %d child1 %d\n",
+ el->getChild(0)->getIndex(),
+ el->getChild(1)->getIndex());
+ printElementInfo(el->getChild(0));
+ printElementInfo(el->getChild(1));
+ }
+ }
+
void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof)
{
@@ -181,7 +233,7 @@ namespace AMDiS {
while (elInfo) {
if (elInfo->getElement()->getIndex() == parEl->getIndex())
std::cout << "[DBG] EL INFO TO " << parEl->getIndex() << ": "
- << " tzpe = " << elInfo->getType() << std::endl;
+ << " type = " << elInfo->getType() << std::endl;
elInfo = stack.traverseNext(elInfo);
}
@@ -230,6 +282,41 @@ namespace AMDiS {
}
+ void printAllDofCoords(FiniteElemSpace *feSpace)
+ {
+ FUNCNAME("printAllDofCoords()");
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+ while (elInfo) {
+ Element *el = elInfo->getElement();
+ for (int i = 0; i <= feSpace->getMesh()->getDim(); i++) {
+ MSG("Coords for DOF %d = %f %f\n",
+ *(el->getDOF(i)), (elInfo->getCoord(i))[0], (elInfo->getCoord(i))[1]);
+ }
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+
+ void getAllDofs(FiniteElemSpace *feSpace, std::set<const DegreeOfFreedom*>& dofs)
+ {
+ FUNCNAME("getAllDofs()");
+
+ ElementDofIterator elDofIter(feSpace);
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ elDofIter.reset(elInfo->getElement());
+ do {
+ dofs.insert(elDofIter.getDofPtr());
+ } while (elDofIter.next());
+ elInfo = stack.traverseNext(elInfo);
+ }
+ }
+
+
void writeMatlabMatrix(DOFMatrix &mat, std::string filename)
{
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h
index 8c502e30..0c724cd2 100644
--- a/AMDiS/src/Debug.h
+++ b/AMDiS/src/Debug.h
@@ -22,6 +22,7 @@
#ifndef AMDIS_DEBUG_H
#define AMDIS_DEBUG_H
+#include <set>
#include "AMDiS_fwd.h"
#include "Global.h"
@@ -85,11 +86,21 @@ namespace AMDiS {
Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
Element* getLevel0ParentElement(Mesh *mesh, Element *el);
+
+ Element* getParentElement(Mesh *mesh, Element *el);
+
+ Element* getParentElement(Mesh *mesh, int elIndex);
+
+ void printElementInfo(Element *el);
void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof);
void printMatValuesStatistics(Matrix<DOFMatrix*> *mat);
+ void printAllDofCoords(FiniteElemSpace *feSpace);
+
+ void getAllDofs(FiniteElemSpace *feSpace, std::set<const DegreeOfFreedom*>& dofs);
+
/** \brief
* Creates a text file storing the value of a sparse matrix. Each line of the file
* has three columns:
diff --git a/AMDiS/src/ElementFileWriter.cc b/AMDiS/src/ElementFileWriter.cc
index f312cc96..188840cc 100644
--- a/AMDiS/src/ElementFileWriter.cc
+++ b/AMDiS/src/ElementFileWriter.cc
@@ -280,6 +280,7 @@ namespace AMDiS {
fout.close();
}
+
void ElementFileWriter::writeVtkValues(std::string filename)
{
FUNCNAME("ElementFileWriter::writeVtkValues()");
@@ -291,7 +292,6 @@ namespace AMDiS {
int dimOfWorld = Global::getGeo(WORLD);
int vertices = mesh->getGeo(VERTEX);
int nElements = mesh->getNumberOfLeaves();
- double val;
fout << "<?xml version=\"1.0\"?>" << std::endl;
fout << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
@@ -304,10 +304,8 @@ namespace AMDiS {
// === Write vertex coordinates (every vertex for every element). ===
TraverseStack stack;
-
- ElInfo *elInfo = stack.traverseFirst(mesh,
- -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+ ElInfo *elInfo =
+ stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
// Write coordinates of all element vertices.
@@ -369,18 +367,19 @@ namespace AMDiS {
fout.setf(std::ios::scientific,std::ios::floatfield);
- elInfo = stack.traverseFirst(mesh,
- -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+ elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
int vc = 0;
while (elInfo) {
// Get element value.
- val = vec[elInfo->getElement()->getIndex()];
+ double val = vec[elInfo->getElement()->getIndex()];
+ if (vc == 1101) {
+ MSG("FOUND EL %d\n", elInfo->getElement()->getIndex());
+ }
// Write value for each vertex of each element.
- for (int i = 0; i <= dim; i++) {
+ for (int i = 0; i <= dim; i++)
fout << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
- }
+
vc++;
elInfo = stack.traverseNext(elInfo);
}
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 79415e2e..a57bed36 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -4217,15 +4217,15 @@ namespace AMDiS {
}
}
- // ====== coarseInter functions ======
void Lagrange::coarseInter0(DOFIndexed<double> *drv, RCNeighbourList* list,
int n, BasisFunction* basFct)
{
FUNCNAME("Lagrange::coarseInter0()");
- TEST_EXIT(drv->getFeSpace())("No fe_space in dof_real_vec!\n");
- TEST_EXIT(drv->getFeSpace()->getBasisFcts())("No basis functions in fe_space!\n");
+ TEST_EXIT_DBG(drv->getFeSpace())("No fe_space in dof_real_vec!\n");
+ TEST_EXIT_DBG(drv->getFeSpace()->getBasisFcts())
+ ("No basis functions in fe_space!\n");
if (n < 1)
return;
@@ -4241,13 +4241,15 @@ namespace AMDiS {
(*drv)[pdof] = (*drv)[cdof];
}
+
void Lagrange::coarseInter2_1d(DOFIndexed<double> *drv, RCNeighbourList *list,
int n, BasisFunction* basFct)
{
- FUNCNAME("real_coarseInter2()");
+ FUNCNAME("Lagrange::coarseInter2_1d()");
- TEST_EXIT(drv->getFeSpace())("No fe_space in dof_real_vec!\n");
- TEST_EXIT(drv->getFeSpace()->getBasisFcts())("No basis functions in fe_space!\n");
+ TEST_EXIT_DBG(drv->getFeSpace())("No fe_space in dof_real_vec!\n");
+ TEST_EXIT_DBG(drv->getFeSpace()->getBasisFcts())
+ ("No basis functions in fe_space!\n");
if (n < 1)
return;
@@ -4257,12 +4259,13 @@ namespace AMDiS {
Mesh *mesh = const_cast<Mesh*>(drv->getFeSpace()->getMesh());
// values on child[0]
- DegreeOfFreedom cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX)+1,
+ DegreeOfFreedom cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX) + 1,
admin->getNumberOfPreDOFs(VERTEX));
DegreeOfFreedom pdof = el->getDOF(mesh->getNode(CENTER), admin->getNumberOfPreDOFs(CENTER));
(*drv)[pdof] = (*drv)[cdof];
}
+
void Lagrange::coarseInter2_2d(DOFIndexed<double> *drv, RCNeighbourList* list,
int n, BasisFunction* basFct)
{
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 35612050..0d4433a1 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -58,11 +58,10 @@ namespace AMDiS {
const Flag Mesh::CALL_LEAF_EL = 0X0800L;
const Flag Mesh::CALL_LEAF_EL_LEVEL = 0X1000L;
const Flag Mesh::CALL_EL_LEVEL = 0X2000L;
- const Flag Mesh::CALL_MG_LEVEL = 0X4000L ; // used in mg methods
+ const Flag Mesh::CALL_MG_LEVEL = 0X4000L;
+ const Flag Mesh::CALL_REVERSE_MODE = 0X8000L;
- // const Flag Mesh::USE_PARAMETRIC = 0X8000L ; // used in mg methods
-
std::vector<DegreeOfFreedom> Mesh::dof_used;
const int Mesh::MAX_DOF = 100;
std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index 53578c78..b4cca8cf 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -644,6 +644,9 @@ namespace AMDiS {
///
static const Flag CALL_MG_LEVEL;
+ /// If set, left and right children are swapped in traverse.
+ static const Flag CALL_REVERSE_MODE;
+
protected:
///
bool findElementAtPointRecursive(ElInfo *elinfo,
diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc
index ecce7dea..94637989 100644
--- a/AMDiS/src/MeshStructure.cc
+++ b/AMDiS/src/MeshStructure.cc
@@ -69,6 +69,13 @@ namespace AMDiS {
int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
+
+ if (debugMode) {
+ MSG("addAlondSide(%d, %d, %d, %d)\n",
+ bound.elIndex, bound.ithObj, bound.elType, bound.reverseMode);
+ MSG("Element is leaf: %d\n", el->isLeaf());
+ MSG("s1 = %d s2 = %d\n", s1, s2);
+ }
if (!el->isLeaf()) {
if (s1 == -1)
@@ -102,16 +109,24 @@ namespace AMDiS {
int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
+ if (debugMode) {
+ MSG("Child index %d %d\n",
+ el->getFirstChild()->getIndex(),
+ el->getSecondChild()->getIndex());
+ MSG("s1 = %d s2 = %d\n", s1, s2);
+ MSG(" \n");
+ }
+
if (!reverseOrder) {
if (s1 != -1)
- addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), false);
+ addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
if (s2 != -1)
- addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), false);
+ addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
} else {
if (s2 != -1)
- addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), false);
+ addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
if (s1 != -1)
- addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), false);
+ addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
}
}
}
diff --git a/AMDiS/src/MeshStructure.h b/AMDiS/src/MeshStructure.h
index 3594b7a5..2fac234a 100644
--- a/AMDiS/src/MeshStructure.h
+++ b/AMDiS/src/MeshStructure.h
@@ -112,22 +112,24 @@ namespace AMDiS {
{
FUNCNAME("MeshStructure::print()");
+ std::stringstream oss;
+
if (empty()) {
- std::cout << "-" << std::endl;
- return;
- }
-
- reset();
- bool cont = true;
- while (cont) {
- if (isLeafElement())
- std::cout << "0";
- else
- std::cout << "1";
-
- cont = nextElement();
+ oss << "-" << std::endl;
+ } else {
+ reset();
+ bool cont = true;
+ while (cont) {
+ if (isLeafElement())
+ oss << "0";
+ else
+ oss << "1";
+
+ cont = nextElement();
+ }
}
- std::cout << std::endl;
+
+ MSG("Mesh structure code: %s\n", oss.str().c_str());
}
/// Returns the mesh structure code.
diff --git a/AMDiS/src/ProblemInstat.cc b/AMDiS/src/ProblemInstat.cc
index 1ed2ffe0..e7d6ffeb 100644
--- a/AMDiS/src/ProblemInstat.cc
+++ b/AMDiS/src/ProblemInstat.cc
@@ -94,6 +94,7 @@ namespace AMDiS {
WARNING("no oldSolution created\n");
}
+
void ProblemInstatScal::createUhOld()
{
if (oldSolution) {
@@ -115,12 +116,14 @@ namespace AMDiS {
problemStat->writeFiles(adaptInfo, force);
}
+
void ProblemInstatVec::closeTimestep(AdaptInfo *adaptInfo)
{
bool force = (adaptInfo->getTime() >= adaptInfo->getEndTime());
problemStat->writeFiles(adaptInfo, force);
}
+
ProblemInstatVec::ProblemInstatVec(std::string sname,
ProblemVec *prob, ProblemStatBase *initialProb)
: ProblemInstat(sname, initialProb),
@@ -128,12 +131,14 @@ namespace AMDiS {
oldSolution(NULL)
{}
+
ProblemInstatVec::ProblemInstatVec(std::string sname, ProblemVec &prob)
: ProblemInstat(sname, NULL),
problemStat(&prob),
oldSolution(NULL)
{}
+
ProblemInstatVec::ProblemInstatVec(std::string sname,
ProblemVec &prob, ProblemStatBase &initialProb)
: ProblemInstat(sname, &initialProb),
@@ -141,11 +146,13 @@ namespace AMDiS {
oldSolution(NULL)
{}
+
ProblemInstatVec::~ProblemInstatVec()
{
delete oldSolution;
}
+
void ProblemInstatVec::initialize(Flag initFlag, ProblemInstat *adoptProblem,
Flag adoptFlag)
{
@@ -173,6 +180,7 @@ namespace AMDiS {
WARNING("no oldSolution created\n");
}
+
void ProblemInstatVec::createUhOld()
{
if (oldSolution) {
@@ -193,11 +201,13 @@ namespace AMDiS {
}
}
+
void ProblemInstatScal::initTimestep(AdaptInfo *adaptInfo)
{
oldSolution->copy(*(problemStat->getSolution()));
}
+
void ProblemInstatVec::initTimestep(AdaptInfo *adaptInfo)
{
oldSolution->copy(*(problemStat->getSolution()));
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 03baff91..a8f12021 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -226,6 +226,7 @@ namespace AMDiS {
doOtherStuff();
}
+
void ProblemVec::createMesh()
{
FUNCNAME("ProblemVec::createMesh()");
@@ -273,6 +274,7 @@ namespace AMDiS {
}
}
+
void ProblemVec::createFeSpace(DOFAdmin *admin)
{
FUNCNAME("ProblemVec::createFeSpace()");
@@ -322,6 +324,7 @@ namespace AMDiS {
}
}
+
void ProblemVec::createMatricesAndVectors()
{
FUNCNAME("ProblemVec::createMatricesAndVectors()");
@@ -349,6 +352,7 @@ namespace AMDiS {
}
}
+
void ProblemVec::createSolver()
{
FUNCNAME("ProblemVec::createSolver()");
@@ -424,6 +428,7 @@ namespace AMDiS {
}
}
+
void ProblemVec::createFileWriter()
{
FUNCNAME("ProblemVec::createFileWriter()");
@@ -479,10 +484,12 @@ namespace AMDiS {
#endif
}
+
void ProblemVec::doOtherStuff()
{
}
+
void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
{
FUNCNAME("Problem::solve()");
@@ -670,7 +677,7 @@ namespace AMDiS {
for (int j = 0; j < nComponents; j++) {
if (writeAsmInfo)
- std::cout << "-------" << i << " " << j << "----------------" << std::endl;
+ MSG("--------- %d %d -------------\n", i, j);
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
@@ -750,10 +757,10 @@ namespace AMDiS {
solution->getDOFVector(i),
componentMeshes[i],
assembleFlag);
- }
+ }
if (asmMatrix) {
- solverMatrix.setMatrix(*systemMatrix);
+ solverMatrix.setMatrix(*systemMatrix);
createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
diff --git a/AMDiS/src/RefinementManager.cc b/AMDiS/src/RefinementManager.cc
index c6dd1e13..9386be79 100644
--- a/AMDiS/src/RefinementManager.cc
+++ b/AMDiS/src/RefinementManager.cc
@@ -68,7 +68,7 @@ namespace AMDiS {
delete stack;
nElements -= mesh->getNumberOfLeaves();
-
+
if (nElements != 0) {
aMesh->incChangeIndex();
return MESH_REFINED;
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 6a48d8dd..07fb4b5a 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -16,6 +16,10 @@ namespace AMDiS {
ElInfo* TraverseStack::traverseFirst(Mesh *mesh, int level, Flag fill_flag)
{
FUNCNAME("TraverseStack::traverseFirst()");
+
+ TEST_EXIT_DBG(fill_flag.isSet(Mesh::CALL_REVERSE_MODE) == false ||
+ fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER))
+ ("Not yet implemented!\n");
traverse_mesh = mesh;
traverse_level = level;
@@ -343,13 +347,20 @@ namespace AMDiS {
if (stack_used >= stack_size - 1)
enlargeTraverseStack();
- int i = info_stack[stack_used];
+ int fillIthChild = info_stack[stack_used];
info_stack[stack_used]++;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) {
+ if (fillIthChild == 0)
+ fillIthChild = 1;
+ else
+ fillIthChild = 0;
+ }
+
+ elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]);
stack_used++;
TEST_EXIT_DBG(stack_used < stack_size)
- ("stack_size=%d too small, level=%d\n",
+ ("stack_size = %d too small, level = %d\n",
stack_size, elinfo_stack[stack_used]->getLevel());
info_stack[stack_used] = 0;
@@ -763,7 +774,7 @@ namespace AMDiS {
info_stack[stack_used] = i + 1;
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
- nb = 1-i;
+ nb = 1 - i;
}
/****************************************************************************/
@@ -780,8 +791,8 @@ namespace AMDiS {
if (nb < 2) { /* go down one level in hierarchy */
if (stack_used >= stack_size-1)
enlargeTraverseStack();
- elinfo_stack[stack_used+1]->fillElInfo(1-nb, elinfo_stack[stack_used]);
- info_stack[stack_used] = 2-nb;
+ elinfo_stack[stack_used + 1]->fillElInfo(1 - nb, elinfo_stack[stack_used]);
+ info_stack[stack_used] = 2 - nb;
stack_used++;
nb = 2;
}
diff --git a/AMDiS/src/parallel/GlobalMatrixSolver.cc b/AMDiS/src/parallel/GlobalMatrixSolver.cc
index 19640da4..5616bf83 100644
--- a/AMDiS/src/parallel/GlobalMatrixSolver.cc
+++ b/AMDiS/src/parallel/GlobalMatrixSolver.cc
@@ -3,6 +3,7 @@
#include "Debug.h"
#include "SystemVector.h"
#include "parallel/StdMpi.h"
+#include "parallel/ParallelDomainDbg.h"
#include "petscksp.h"
@@ -190,6 +191,8 @@ namespace AMDiS {
void GlobalMatrixSolver::setDofVector(Vec& petscVec, DOFVector<double>* vec,
int dispMult, int dispAdd)
{
+ FUNCNAME("GlobalMatrixSolver::setDofVector()");
+
// Traverse all used dofs in the dof vector.
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
@@ -215,7 +218,7 @@ namespace AMDiS {
double value = *dofIt;
VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
}
- }
+ }
}
@@ -402,7 +405,9 @@ namespace AMDiS {
if (!d_nnz || meshDistributor->getLastMeshChangeIndex() != lastMeshNnz) {
if (d_nnz) {
delete [] d_nnz;
+ d_nnz = NULL;
delete [] o_nnz;
+ o_nnz = NULL;
}
createPetscNnzStructure(mat);
@@ -430,7 +435,7 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++)
for (int j = 0; j < nComponents; j++)
if ((*mat)[i][j])
- setDofMatrix((*mat)[i][j], nComponents, i, j);
+ setDofMatrix((*mat)[i][j], nComponents, i, j);
INFO(info, 8)("Fill petsc matrix 2 needed %.5f seconds\n", TIME_USED(first, clock()));
@@ -474,12 +479,18 @@ namespace AMDiS {
// Do not delete the solution vector, use it for the initial guess.
// KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
+#if (DEBUG != 0)
+ // ParallelDomainDbg::writeCoordsFile(*meshDistributor, "mpi-coords", "dat");
+#endif
+
// === Run Petsc. ===
KSPSolve(solver, petscRhsVec, petscSolVec);
+
// === Transfere values from Petsc's solution vectors to the dof vectors.
+
PetscScalar *vecPointer;
VecGetArray(petscSolVec, &vecPointer);
@@ -488,7 +499,7 @@ namespace AMDiS {
DOFVector<double> *dofvec = vec.getDOFVector(i);
for (int j = 0; j < nRankDofs; j++)
(*dofvec)[meshDistributor->mapLocalToDofIndex(j)] =
- vecPointer[j * nComponents + i];
+ vecPointer[j * nComponents + i];
}
VecRestoreArray(petscSolVec, &vecPointer);
diff --git a/AMDiS/src/parallel/MpiHelper.cc b/AMDiS/src/parallel/MpiHelper.cc
index b4663288..ff0d193c 100644
--- a/AMDiS/src/parallel/MpiHelper.cc
+++ b/AMDiS/src/parallel/MpiHelper.cc
@@ -11,6 +11,12 @@ namespace AMDiS {
double valCopy = value;
MPI::COMM_WORLD.Allreduce(&valCopy, &value, 1, MPI_DOUBLE, MPI_SUM);
}
+
+ void globalAdd(int &value)
+ {
+ int valCopy = value;
+ MPI::COMM_WORLD.Allreduce(&valCopy, &value, 1, MPI_INT, MPI_SUM);
+ }
void globalMin(double &value)
{
diff --git a/AMDiS/src/parallel/MpiHelper.h b/AMDiS/src/parallel/MpiHelper.h
index b7e3f229..bf8557e4 100644
--- a/AMDiS/src/parallel/MpiHelper.h
+++ b/AMDiS/src/parallel/MpiHelper.h
@@ -31,6 +31,8 @@ namespace AMDiS {
namespace mpi {
void globalAdd(double &value);
+ void globalAdd(int &value);
+
void globalMin(double &value);
void globalMax(double &value);
diff --git a/AMDiS/src/parallel/ParallelDomainBase.cc b/AMDiS/src/parallel/ParallelDomainBase.cc
index 4d5a872f..a7ddc962 100644
--- a/AMDiS/src/parallel/ParallelDomainBase.cc
+++ b/AMDiS/src/parallel/ParallelDomainBase.cc
@@ -49,7 +49,8 @@ namespace AMDiS {
nOverallDofs(0),
rstart(0),
deserialized(false),
- lastMeshChangeIndex(0)
+ lastMeshChangeIndex(0),
+ macroElementStructureConsisten(false)
{
FUNCNAME("MeshDistributor::ParalleDomainBase()");
@@ -117,14 +118,18 @@ namespace AMDiS {
createLocalGlobalNumbering();
+
// === Remove all macro elements that are not part of the rank partition. ===
removeMacroElements();
+ macroElementStructureConsisten = true;
+
// === Reset all DOFAdmins of the mesh. ===
updateDofAdmins();
+
// === If in debug mode, make some tests. ===
#if (DEBUG != 0)
@@ -137,10 +142,12 @@ namespace AMDiS {
debug::writeMesh(feSpace, -1, "macro_mesh");
#endif
+
// === Create periodic dof mapping, if there are periodic boundaries. ===
createPeriodicMap();
+
// === Global refinements. ===
int globalRefinement = 0;
@@ -160,6 +167,7 @@ namespace AMDiS {
createPeriodicMap();
}
+
/// === Set DOF rank information to all matrices and vectors. ===
for (unsigned int i = 0; i < probStat.size(); i++) {
@@ -177,6 +185,7 @@ namespace AMDiS {
}
}
+
// === Remove periodic boundary conditions in sequential problem definition. ===
// Remove periodic boundaries in boundary manager on matrices and vectors.
@@ -403,6 +412,10 @@ namespace AMDiS {
{
FUNCNAME("MeshDistributor::checkMeshChange()");
+#if (DEBUG != 0)
+ debug::writeMesh(feSpace, -1, "before_check_mesh");
+#endif
+
// === If mesh has not been changed on all ranks, return. ===
int recvAllValues = 0;
@@ -460,6 +473,7 @@ namespace AMDiS {
// === Because the mesh has been changed, update the DOF numbering and mappings. ===
+ mesh->dofCompress();
updateLocalGlobalNumbering();
}
@@ -508,7 +522,8 @@ namespace AMDiS {
boundIt->rankObj.el,
boundIt->rankObj.subObj,
boundIt->rankObj.ithObj,
- boundIt->rankObj.elType);
+ boundIt->rankObj.elType,
+ boundIt->rankObj.reverseMode);
if (b)
meshFitTogether = false;
@@ -523,10 +538,11 @@ namespace AMDiS {
bool MeshDistributor::fitElementToMeshCode(MeshStructure &code,
- Element *el,
- GeoIndex subObj,
- int ithObj,
- int elType)
+ Element *el,
+ GeoIndex subObj,
+ int ithObj,
+ int elType,
+ bool reverseMode)
{
FUNCNAME("MeshDistributor::fitElementToMeshCode()");
@@ -541,16 +557,22 @@ namespace AMDiS {
TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");
bool meshChanged = false;
+ Flag traverseFlag =
+ Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;
+
+ if (reverseMode) {
+ std::swap(s1, s2);
+ traverseFlag |= Mesh::CALL_REVERSE_MODE;
+ }
if (s1 != -1 && s2 != -1) {
TraverseStack stack;
- ElInfo *elInfo =
- stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
+ ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
while (elInfo && elInfo->getElement() != el) {
elInfo = stack.traverseNext(elInfo);
}
- meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType);
+ meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
return meshChanged;
}
@@ -561,8 +583,7 @@ namespace AMDiS {
meshChanged = true;
TraverseStack stack;
- ElInfo *elInfo =
- stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
+ ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
while (elInfo && elInfo->getElement() != el) {
elInfo = stack.traverseNext(elInfo);
@@ -576,28 +597,31 @@ namespace AMDiS {
refineManager->refineFunction(elInfo);
}
+ Element *child0 = el->getFirstChild();
+ Element *child1 = el->getSecondChild();
+ if (reverseMode)
+ std::swap(child0, child1);
+
if (s1 != -1) {
TraverseStack stack;
- ElInfo *elInfo =
- stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
+ ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
- while (elInfo && elInfo->getElement() != el->getFirstChild()) {
+ while (elInfo && elInfo->getElement() != child0) {
elInfo = stack.traverseNext(elInfo);
}
meshChanged |=
- fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType));
+ fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
} else {
TraverseStack stack;
- ElInfo *elInfo =
- stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
+ ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
- while (elInfo && elInfo->getElement() != el->getSecondChild()) {
+ while (elInfo && elInfo->getElement() != child1) {
elInfo = stack.traverseNext(elInfo);
}
meshChanged |=
- fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType));
+ fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
}
return meshChanged;
@@ -605,10 +629,11 @@ namespace AMDiS {
bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code,
- TraverseStack &stack,
- GeoIndex subObj,
- int ithObj,
- int elType)
+ TraverseStack &stack,
+ GeoIndex subObj,
+ int ithObj,
+ int elType,
+ bool reverseMode)
{
FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
@@ -643,24 +668,29 @@ namespace AMDiS {
int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
+ Element *child0 = el->getFirstChild();
+ Element *child1 = el->getSecondChild();
+ if (reverseMode) {
+ std::swap(s1, s2);
+ std::swap(child0, child1);
+ }
if (s1 != -1) {
stack.traverseNext(elInfo);
code.nextElement();
- value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType));
+ value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
elInfo = stack.getElInfo();
} else {
do {
elInfo = stack.traverseNext(elInfo);
- } while (elInfo && elInfo->getElement() != el->getSecondChild());
+ } while (elInfo && elInfo->getElement() != child1);
}
- TEST_EXIT_DBG(elInfo->getElement() == el->getSecondChild())
- ("This should not happen!\n");
+ TEST_EXIT_DBG(elInfo->getElement() == child1)("This should not happen!\n");
if (s2 != -1) {
code.nextElement();
- value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType));
+ value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
} else {
int level = elInfo->getLevel();
@@ -685,7 +715,7 @@ namespace AMDiS {
void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
- std::map<int, const DegreeOfFreedom*> &dofMap)
+ std::map<int, const DegreeOfFreedom*> &dofMap)
{
FUNCNAME("MeshDistributor::deserialize()");
@@ -717,7 +747,7 @@ namespace AMDiS {
void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
- std::map<int, const DegreeOfFreedom*> &dofMap)
+ std::map<int, const DegreeOfFreedom*> &dofMap)
{
int mapSize = 0;
SerUtil::deserialize(in, mapSize);
@@ -1223,7 +1253,7 @@ namespace AMDiS {
// === the necessary boundary objects. ===
for (std::map<int, std::vector<BoundaryObject> >::iterator it = recvObjs.begin();
- it != recvObjs.end(); ++it) {
+ it != recvObjs.end(); ++it) {
if (it->second.size() > 0) {
TEST_EXIT_DBG(it->second.size() == stdMpiObjs.getRecvData(it->first).size())
("This should not happen!\n");
@@ -1297,7 +1327,7 @@ namespace AMDiS {
for (std::map<int, std::vector<BoundaryObject> >::iterator it = recvObjs.begin();
- it != recvObjs.end(); ++it) {
+ it != recvObjs.end(); ++it) {
if (it->second.size() > 0) {
TEST_EXIT_DBG(it->second.size() == stdMpiObjs.getRecvData(it->first).size())
("Should not happen!\n");
@@ -1389,7 +1419,7 @@ namespace AMDiS {
// === Create information which dof indices must be send and which must ===
- // === be received. ===
+ // === be received to/from other ranks. ===
// Stores to each rank a map from DOF pointers (which are all owned by the rank
// and lie on an interior boundary) to new global DOF indices.
@@ -1447,7 +1477,7 @@ namespace AMDiS {
std::map<int, DofMapVec> &dofMap = stdMpi.getRecvData();
- // === Change dof indices at boundary from other ranks. ===
+ // === Change global dof indices at boundary from other ranks. ===
// Within this small data structure we track which dof index was already changed.
// This is used to avoid the following situation: Assume, there are two dof indices
@@ -1463,10 +1493,10 @@ namespace AMDiS {
++dofIt)
dofChanged[dofIt->first] = false;
- for (std::map<int, std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > >::iterator rankIt = dofMap.begin();
+ for (std::map<int, DofMapVec>::iterator rankIt = dofMap.begin();
rankIt != dofMap.end(); ++rankIt) {
- for (std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> >::iterator recvDofIt = rankIt->second.begin();
+ for (DofMapVec::iterator recvDofIt = rankIt->second.begin();
recvDofIt != rankIt->second.end(); ++recvDofIt) {
DegreeOfFreedom oldDof = recvDofIt->first;
@@ -1511,7 +1541,7 @@ namespace AMDiS {
#if (DEBUG != 0)
debug::ElementIdxToDofs elMap;
- debug::createSortedDofs(mesh, elMap);
+ debug::createSortedDofs(mesh, elMap);
#endif
typedef std::set<const DegreeOfFreedom*> DofSet;
@@ -1530,17 +1560,25 @@ namespace AMDiS {
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
- Element *element = elInfo->getElement();
- elDofIt.reset(element);
+ elDofIt.reset(elInfo->getElement());
do {
rankDofSet.insert(elDofIt.getDofPtr());
+ vertexDof[elDofIt.getDofPtr()] = false;
+ } while(elDofIt.next());
- if (elDofIt.getCurrentPos() == 0)
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
+ while (elInfo) {
+ Element *element = elInfo->getElement();
+ elDofIt.reset(element);
+ elDofIt.reset(elInfo->getElement());
+ do {
+ if (elDofIt.getCurrentPos() == 0)
vertexDof[elDofIt.getDofPtr()] = true;
- else
- vertexDof[elDofIt.getDofPtr()] = false;
} while(elDofIt.next());
-
+
elInfo = stack.traverseNext(elInfo);
}
@@ -1566,7 +1604,7 @@ namespace AMDiS {
dofIt != it->second.end(); ++dofIt)
if (oldVertexDof[*dofIt] && vertexDof[*dofIt])
sendDofs[it->first].push_back(*dofIt);
-
+
for (RankToDofContainer::iterator it = oldRecvDofs.begin();
it != oldRecvDofs.end(); ++it) {
@@ -1578,22 +1616,27 @@ namespace AMDiS {
DofContainer::iterator eraseIt =
find(rankDofs.begin(), rankDofs.end(), *dofIt);
if (eraseIt != rankDofs.end())
- rankDofs.erase(eraseIt);
+ rankDofs.erase(eraseIt);
}
}
}
+
for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
DofContainer dofs;
it->rankObj.el->getVertexDofs(feSpace, it->rankObj, dofs);
it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs);
-
+
for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
- /* FOR DEBUGGING
+ // FOR DEBUGGING
+ /*
WorldVector<double> cs;
mesh->getDofIndexCoords(dofs[i], feSpace, cs);
MSG("SEND EL %d DOF %d TO %d\n", it->rankObj.elIndex, *(dofs[i]), it.getRank());
- MSG("COORDS: %f %f %f\n", cs[0], cs[1], cs[2]);
+ if (cs.getSize() == 2)
+ MSG("COORDS-s2: %f %f\n", cs[0], cs[1]);
+ else
+ MSG("COORDS-s2: %f %f %f\n", cs[0], cs[1], cs[2]);
*/
sendDofs[it.getRank()].push_back(dofs[i]);
@@ -1612,17 +1655,22 @@ namespace AMDiS {
if (eraseIt != rankDofs.end())
rankDofs.erase(eraseIt);
- /* FOR DEBUGGING
- WorldVector<double> cs;
- mesh->getDofIndexCoords(dofs[i], feSpace, cs);
- MSG("RECV EL %d DOF %d FROM %d\n", it->rankObj.elIndex, *(dofs[i]), it.getRank());
- MSG("COORDS: %f %f %f\n", cs[0], cs[1], cs[2]);
+ // FOR DEBUGGING
+ /*
+ WorldVector<double> cs;
+ mesh->getDofIndexCoords(dofs[i], feSpace, cs);
+ MSG("RECV EL %d DOF %d FROM %d\n", it->rankObj.elIndex, *(dofs[i]), it.getRank());
+ if (cs.getSize() == 2)
+ MSG("COORDS-r2: %f %f\n", cs[0], cs[1]);
+ else
+ MSG("COORDS-r2: %f %f %f\n", cs[0], cs[1], cs[2]);
*/
recvDofs[it.getRank()].push_back(dofs[i]);
}
}
+
nRankDofs = rankDofs.size();
// === Get starting position for global rank dof ordering. ====
@@ -1643,7 +1691,6 @@ namespace AMDiS {
int i = 0;
for (DofContainer::iterator dofIt = rankAllDofs.begin();
dofIt != rankAllDofs.end(); ++dofIt) {
-
rankDofsNewLocalIndex[*dofIt] = i;
// First, we set all dofs in ranks partition to be owend by the rank. Later,
// the dofs in ranks partition that are owned by other rank are set to false.
@@ -1709,45 +1756,104 @@ namespace AMDiS {
#if (DEBUG != 0)
debug::testSortedDofs(mesh, elMap);
- debug::writeMesh(feSpace, -1, "last_update_mesh");
ParallelDomainDbg::testCommonDofs(*this, true);
+
+#if 0
+ MSG("------------- Debug information -------------\n");
+ MSG("nRankDofs = %d\n", nRankDofs);
+ MSG("nOverallDofs = %d\n", nOverallDofs);
+ MSG("rstart %d\n", rstart);
+ for (RankToDofContainer::iterator it = sendDofs.begin();
+ it != sendDofs.end(); ++it) {
+ MSG("send %d DOFs to rank %d\n", it->second.size(), it->first);
+
+ for (unsigned int i = 0; i < it->second.size(); i++)
+ MSG("%d send DOF: %d\n", i, *(it->second[i]));
+ }
+ for (RankToDofContainer::iterator it = recvDofs.begin();
+ it != recvDofs.end(); ++it) {
+ MSG("recv %d DOFs from rank %d\n", it->second.size(), it->first);
+
+ for (unsigned int i = 0; i < it->second.size(); i++)
+ MSG("%d recv DOF: %d\n", i, *(it->second[i]));
+ }
+
+ std::set<const DegreeOfFreedom*> testDofs;
+ debug::getAllDofs(feSpace, testDofs);
+ for (std::set<const DegreeOfFreedom*>::iterator it = testDofs.begin();
+ it != testDofs.end(); ++it) {
+ MSG("DOF %d: mapLocalGlobalDofs = %d vertexDof = %d isRankDof = %d\n", **it,
+ mapLocalGlobalDofs[**it],
+ vertexDof[*it],
+ isRankDof[**it]);
+ }
+ MSG("\n");
+ for (DofMapping::iterator it = mapLocalDofIndex.begin();
+ it != mapLocalDofIndex.end(); ++it) {
+ MSG("mapLocalDofIndex[%d] = %d\n", it->first, it->second);
+ }
+#endif
#endif
}
void MeshDistributor::createLocalMappings(DofIndexMap &rankDofsNewLocalIndex,
- DofIndexMap &rankOwnedDofsNewLocalIndex,
- DofIndexMap &rankDofsNewGlobalIndex)
+ DofIndexMap &rankOwnedDofsNewLocalIndex,
+ DofIndexMap &rankDofsNewGlobalIndex)
{
+ FUNCNAME("MeshDistributor::createLocalMappings()");
+
+#if (DEBUG != 0)
+ if (macroElementStructureConsisten) {
+ std::set<const DegreeOfFreedom*> allDofs;
+ debug::getAllDofs(feSpace, allDofs);
+
+ for (std::set<const DegreeOfFreedom*>::iterator it = allDofs.begin();
+ it != allDofs.end(); ++it) {
+ TEST_EXIT(rankDofsNewLocalIndex.count(*it) == 1)
+ ("Missing DOF %d in rankDofsNewLocalIndex!\n", **it);
+
+ TEST_EXIT(rankDofsNewGlobalIndex.count(*it) == 1)
+ ("Missing DOF %d in rankDofsNewGlobalIndex!\n", **it);
+
+ if (isRankDof[**it]) {
+ TEST_EXIT(rankOwnedDofsNewLocalIndex.count(*it) == 1)
+ ("Missing DOF %d in rankOwnedDofsNewLocalIndex!\n", **it);
+ }
+ }
+ }
+#endif
+
mapLocalGlobalDofs.clear();
mapLocalDofIndex.clear();
+ // Iterate over all DOFs in ranks partition.
for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
- DegreeOfFreedom localDof = dofIt->second;
- DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];
+ DegreeOfFreedom newLocalIndex = dofIt->second;
+ DegreeOfFreedom newGlobalIndex = rankDofsNewGlobalIndex[dofIt->first];
- *const_cast<DegreeOfFreedom*>(dofIt->first) = localDof;
- mapLocalGlobalDofs[localDof] = globalDof;
+ *const_cast<DegreeOfFreedom*>(dofIt->first) = newLocalIndex;
+ mapLocalGlobalDofs[newLocalIndex] = newGlobalIndex;
}
for (DofIndexMap::iterator dofIt = rankOwnedDofsNewLocalIndex.begin();
dofIt != rankOwnedDofsNewLocalIndex.end(); ++dofIt)
- mapLocalDofIndex[dofIt->second] = *(dofIt->first);
+ mapLocalDofIndex[dofIt->second] = *(dofIt->first);
}
void MeshDistributor::createDofMemberInfo(DofToPartitions& partitionDofs,
- DofContainer& rankOwnedDofs,
- DofContainer& rankAllDofs,
- DofToRank& boundaryDofs,
- DofToBool& vertexDof)
+ DofContainer& rankOwnedDofs,
+ DofContainer& rankAllDofs,
+ DofToRank& boundaryDofs,
+ DofToBool& isVertexDof)
{
partitionDofs.clear();
rankOwnedDofs.clear();
rankAllDofs.clear();
boundaryDofs.clear();
- vertexDof.clear();
+ isVertexDof.clear();
// === Determine to each dof the set of partitions the dof belongs to. ===
@@ -1763,9 +1869,9 @@ namespace AMDiS {
partitionDofs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
if (elDofIt.getCurrentPos() == 0)
- vertexDof[elDofIt.getDofPtr()] = true;
+ isVertexDof[elDofIt.getDofPtr()] = true;
else
- vertexDof[elDofIt.getDofPtr()] = false;
+ isVertexDof[elDofIt.getDofPtr()] = false;
} while (elDofIt.next());
elInfo = stack.traverseNext(elInfo);
@@ -1817,7 +1923,7 @@ namespace AMDiS {
}
if (!isInRank)
- vertexDof.erase(it->first);
+ isVertexDof.erase(it->first);
}
sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
@@ -1997,7 +2103,7 @@ namespace AMDiS {
elDofIter.reset(el);
do {
dofMap[elDofIter.getDof()] = elDofIter.getDofPtr();
- } while(elDofIter.next());
+ } while (elDofIter.next());
}
myIntBoundary.deserialize(in, elIndexMap);
diff --git a/AMDiS/src/parallel/ParallelDomainBase.h b/AMDiS/src/parallel/ParallelDomainBase.h
index 3b291e5d..c287bdb9 100644
--- a/AMDiS/src/parallel/ParallelDomainBase.h
+++ b/AMDiS/src/parallel/ParallelDomainBase.h
@@ -203,7 +203,6 @@ namespace AMDiS {
*/
void synchVector(SystemVector &vec);
-
protected:
/** \brief
* Determines the interior boundaries, i.e. boundaries between ranks, and stores
@@ -276,7 +275,7 @@ namespace AMDiS {
DofContainer& rankOwnedDofs,
DofContainer& rankAllDofs,
DofToRank& boundaryDofs,
- DofToBool& vertexDof);
+ DofToBool& isVertexDof);
/** \brief
* Checks for all given interior boundaries if the elements fit together on both
@@ -306,13 +305,15 @@ namespace AMDiS {
Element *el,
GeoIndex subObj,
int ithObj,
- int elType);
+ int elType,
+ bool reverseMode);
bool fitElementToMeshCode2(MeshStructure &code,
TraverseStack &stack,
GeoIndex subObj,
int ithObj,
- int elType);
+ int elType,
+ bool reverseMode);
/// Writes a vector of dof pointers to an output stream.
void serialize(std::ostream &out, DofContainer &data);
@@ -516,6 +517,11 @@ namespace AMDiS {
*/
long lastMeshChangeIndex;
+ /// This variable is true, if the macro elements are consistent with all other
+ /// data structures. Within the initialization and during mesh redistribution this
+ /// may not be the case.
+ bool macroElementStructureConsisten;
+
friend class ParallelDomainDbg;
};
}
diff --git a/AMDiS/src/parallel/ParallelDomainDbg.cc b/AMDiS/src/parallel/ParallelDomainDbg.cc
index 4d0026fe..eeed93ba 100644
--- a/AMDiS/src/parallel/ParallelDomainDbg.cc
+++ b/AMDiS/src/parallel/ParallelDomainDbg.cc
@@ -6,6 +6,7 @@
#include "FixVec.h"
#include "StdMpi.h"
#include "Debug.h"
+#include "MpiHelper.h"
namespace AMDiS {
@@ -255,12 +256,19 @@ namespace AMDiS {
for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it)
stdMpi.recv(it->first);
stdMpi.startCommunication<int>(MPI_INT);
-
+
+ int foundError = 0;
for (std::map<int, int>::iterator it = stdMpi.getRecvData().begin();
- it != stdMpi.getRecvData().end(); ++it)
- TEST_EXIT(it->second == static_cast<int>(recvDofs[it->first].size()))
- ("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n",
- recvDofs[it->first].size(), it->first, it->second);
+ it != stdMpi.getRecvData().end(); ++it) {
+ if (it->second != static_cast<int>(recvDofs[it->first].size())) {
+ ERROR("Rank expectes %d DOFs to receive from rank %d, but %d DOFs are received!\n",
+ recvDofs[it->first].size(), it->first, it->second);
+ foundError = 1;
+ }
+ }
+
+ mpi::globalAdd(foundError);
+ TEST_EXIT(foundError == 0)("Error found on at least one rank!\n");
}
@@ -391,4 +399,29 @@ namespace AMDiS {
}
}
+
+ void ParallelDomainDbg::writeCoordsFile(MeshDistributor &pdb,
+ std::string prefix, std::string postfix)
+ {
+ FUNCNAME("ParallelDomainDbg::writeCoordsFile()");
+
+ std::stringstream filename;
+ filename << prefix << "-" << pdb.mpiRank << "." << postfix;
+
+ DOFVector<WorldVector<double> > coords(pdb.feSpace, "tmp");
+ pdb.mesh->getDofIndexCoords(pdb.feSpace, coords);
+
+ std::ofstream file;
+ file.open(filename.str().c_str());
+ file << coords.getUsedSize() << "\n";
+ DOFIterator<WorldVector<double> > it(&coords, USED_DOFS);
+ for (it.reset(); !it.end(); ++it) {
+ file << pdb.mapLocalGlobalDofs[it.getDOFIndex()];
+ for (int i = 0; i < pdb.mesh->getDim(); i++)
+ file << " " << (*it)[i];
+ file << "\n";
+ }
+ file.close();
+ }
+
}
diff --git a/AMDiS/src/parallel/ParallelDomainDbg.h b/AMDiS/src/parallel/ParallelDomainDbg.h
index 1102dd98..d5523b91 100644
--- a/AMDiS/src/parallel/ParallelDomainDbg.h
+++ b/AMDiS/src/parallel/ParallelDomainDbg.h
@@ -112,6 +112,9 @@ namespace AMDiS {
* \param[in] pdb Parallel problem definition used for debugging.
*/
static void printBoundaryInfo(MeshDistributor &pdb);
+
+ static void writeCoordsFile(MeshDistributor &pdb,
+ std::string prefix, std::string postfix);
};
} // namespace AMDiS
--
GitLab