Skip to content

GitLab

Commit 29fd364c authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Fixed a serious bug for mesh coarsening in parallel code.

parent 0a9d458f
Hide whitespace changes
Inline Side-by-side
AMDiS/libtool
View file @ 29fd364c
......@@ -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"
......
AMDiS/src/BasisFunction.h
View file @ 29fd364c
......@@ -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;
......
AMDiS/src/DOFMatrix.cc
View file @ 29fd364c
......@@ -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];
}
}
}
......
AMDiS/src/DOFVector.cc
View file @ 29fd364c
......@@ -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)
{
......
AMDiS/src/Debug.cc
View file @ 29fd364c
......@@ -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;
......
AMDiS/src/Debug.h
View file @ 29fd364c
......@@ -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:
......
AMDiS/src/ElementFileWriter.cc
View file @ 29fd364c
......@@ -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);
}
......
AMDiS/src/Lagrange.cc
View file @ 29fd364c
......@@ -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)
{
......
AMDiS/src/Mesh.cc
View file @ 29fd364c
......@@ -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;
......
AMDiS/src/Mesh.h
View file @ 29fd364c
......@@ -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,
......
AMDiS/src/MeshStructure.cc
View file @ 29fd364c
......@@ -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);
}
}
}
......
AMDiS/src/MeshStructure.h
View file @ 29fd364c
......@@ -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.
......
AMDiS/src/ProblemInstat.cc
View file @ 29fd364c
......@@ -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;
}