From 0ccef766690c8c328d216455e744316651b71336 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 14 Oct 2009 15:54:39 +0000
Subject: [PATCH] Work on periodic boundary conditions in parallel computation.
---
AMDiS/libtool | 64 ++---
AMDiS/src/AMDiS.h | 6 +
AMDiS/src/BoundaryCondition.h | 6 +
AMDiS/src/BoundaryManager.cc | 94 +++----
AMDiS/src/BoundaryManager.h | 47 +++-
AMDiS/src/DOFMatrix.cc | 48 ++--
AMDiS/src/DOFVector.cc | 2 +
AMDiS/src/DataCollector.cc | 3 +-
AMDiS/src/ElInfo2d.cc | 7 +-
AMDiS/src/InteriorBoundary.cc | 2 +-
AMDiS/src/InteriorBoundary.h | 2 +-
AMDiS/src/LeafData.h | 3 +-
AMDiS/src/MacroReader.cc | 19 +-
AMDiS/src/ParallelDomainBase.cc | 470 +++++++++++++++++++++++++------
AMDiS/src/ParallelDomainBase.h | 38 ++-
AMDiS/src/ParallelDomainVec.cc | 32 ++-
AMDiS/src/ParallelDomainVec.h | 3 +
AMDiS/src/PeriodicBC.cc | 110 ++++----
AMDiS/src/PeriodicBC.h | 22 +-
AMDiS/src/ProblemVec.cc | 4 -
AMDiS/src/RCNeighbourList.cc | 317 ++++++++++-----------
AMDiS/src/RefinementManager.cc | 5 +-
AMDiS/src/RefinementManager2d.cc | 34 ++-
AMDiS/src/VertexVector.cc | 39 ++-
AMDiS/src/VertexVector.h | 8 +-
AMDiS/src/VtkWriter.cc | 2 +
26 files changed, 879 insertions(+), 508 deletions(-)
diff --git a/AMDiS/libtool b/AMDiS/libtool
index b24ee23d..4c12058b 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -30,10 +30,10 @@
# the same distribution terms that you use for the rest of that program.
# A sed program that does not truncate output.
-SED="/bin/sed"
+SED="/usr/bin/sed"
# Sed that helps us avoid accidentally triggering echo(1) options like -n.
-Xsed="/bin/sed -e 1s/^X//"
+Xsed="/usr/bin/sed -e 1s/^X//"
# The HP-UX ksh and POSIX shell print the target directory to stdout
# if CDPATH is set.
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host deimos101:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -66,12 +66,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -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
@@ -97,7 +97,7 @@ with_gcc=yes
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -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"
@@ -325,10 +325,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
+sys_lib_search_path_spec=" /usr/lib64/gcc/x86_64-suse-linux/4.1.2/ /usr/lib/gcc/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/ /lib/x86_64-suse-linux/4.1.2/ /lib/../lib64/ /usr/lib/x86_64-suse-linux/4.1.2/ /usr/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../ /lib/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -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 NWRW15:
+# Libtool was configured on host deimos101:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6782,12 +6782,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -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
@@ -6813,7 +6813,7 @@ with_gcc=yes
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -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"
@@ -6942,11 +6942,11 @@ striplib="strip --strip-unneeded"
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects="/usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crti.o /usr/lib/gcc/i386-redhat-linux/4.1.2/crtbeginS.o"
+predep_objects="/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/crti.o /usr/lib64/gcc/x86_64-suse-linux/4.1.2/crtbeginS.o"
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects="/usr/lib/gcc/i386-redhat-linux/4.1.2/crtendS.o /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crtn.o"
+postdep_objects="/usr/lib64/gcc/x86_64-suse-linux/4.1.2/crtendS.o /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/crtn.o"
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -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/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-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/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"
@@ -7038,10 +7038,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
+sys_lib_search_path_spec=" /usr/lib64/gcc/x86_64-suse-linux/4.1.2/ /usr/lib/gcc/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/ /lib/x86_64-suse-linux/4.1.2/ /lib/../lib64/ /usr/lib/x86_64-suse-linux/4.1.2/ /usr/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../ /lib/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host deimos101:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7087,12 +7087,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -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"
@@ -7112,13 +7112,13 @@ LTCFLAGS="-g -O2"
CC="g77"
# Is the compiler the GNU C compiler?
-with_gcc=yes
+with_gcc=
# An ERE matcher.
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -7346,10 +7346,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../ /lib/i386-redhat-linux/3.4.6/ /lib/ /usr/lib/i386-redhat-linux/3.4.6/ /usr/lib/"
+sys_lib_search_path_spec=" /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/ /usr/lib/gcc/x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../ /lib/x86_64-suse-linux/3.3.5/ /lib/ /usr/lib/x86_64-suse-linux/3.3.5/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
diff --git a/AMDiS/src/AMDiS.h b/AMDiS/src/AMDiS.h
index 3b43a746..f0d73e1a 100644
--- a/AMDiS/src/AMDiS.h
+++ b/AMDiS/src/AMDiS.h
@@ -92,4 +92,10 @@
#include "VtkWriter.h"
#include "parareal/ProblemBase.h"
#include "parareal/AdaptParaReal.h"
+
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+#include "ParallelDomainScal.h"
+#include "ParallelDomainVec.h"
+#endif
+
#endif
diff --git a/AMDiS/src/BoundaryCondition.h b/AMDiS/src/BoundaryCondition.h
index b48e407c..fb7f312b 100644
--- a/AMDiS/src/BoundaryCondition.h
+++ b/AMDiS/src/BoundaryCondition.h
@@ -118,6 +118,12 @@ namespace AMDiS {
return false;
}
+ /// Returns whether the boundary condition is a periodic condition or not.
+ virtual bool isPeriodic()
+ {
+ return false;
+ }
+
/** \brief
* In some situations it may be required to set Dirichlet boundary conditions,
* but not to apply them to the matrix. This is for example the case, if the
diff --git a/AMDiS/src/BoundaryManager.cc b/AMDiS/src/BoundaryManager.cc
index e29a76c7..0aa2e763 100644
--- a/AMDiS/src/BoundaryManager.cc
+++ b/AMDiS/src/BoundaryManager.cc
@@ -9,6 +9,9 @@
namespace AMDiS {
+ std::map<BoundaryType, std::vector<BoundaryCondition*> >
+ BoundaryManager::globalBoundaryMap;
+
BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
{
localBounds.resize(omp_get_overall_max_threads());
@@ -39,8 +42,7 @@ namespace AMDiS {
const DOFVectorBase<double> *dv)
{
double result = 0.0;
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
if ((*it).second)
result += (*it).second->boundResidual(elInfo, matrix, dv);
@@ -55,7 +57,6 @@ namespace AMDiS {
Vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
// get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()];
@@ -65,21 +66,20 @@ namespace AMDiS {
basisFcts->getLocalIndicesVec(elInfo->getElement(), feSpace->getAdmin(), &dofVec);
// apply non dirichlet boundary conditions
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
if ((*it).second && !(*it).second->isDirichlet())
(*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
localBound, nBasFcts);
// apply dirichlet boundary conditions
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
if ((*it).second && (*it).second->isDirichlet())
(*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0],
localBound, nBasFcts);
}
}
- void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo,
- DOFMatrix *mat)
+ void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, DOFMatrix *mat)
{
if (localBCs.size() <= 0)
return;
@@ -88,7 +88,6 @@ namespace AMDiS {
Vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
// get boundaries of all DOFs
BoundaryType *localBound = localBounds[omp_get_thread_num()];
@@ -99,75 +98,60 @@ namespace AMDiS {
feSpace->getAdmin(), &dofVec);
// apply non dirichlet boundary conditions
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if (!(*it).second->isDirichlet())
- (*it).second->fillBoundaryCondition(mat, elInfo, &dofVec[0],
- localBound, nBasFcts);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && !(*it).second->isDirichlet())
+ (*it).second->fillBoundaryCondition(mat, elInfo, &dofVec[0],
+ localBound, nBasFcts);
// apply dirichlet boundary conditions
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if ((*it).second->isDirichlet())
- (*it).second->fillBoundaryCondition(mat, elInfo, &dofVec[0],
- localBound, nBasFcts);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && (*it).second->isDirichlet())
+ (*it).second->fillBoundaryCondition(mat, elInfo, &dofVec[0],
+ localBound, nBasFcts);
}
void BoundaryManager::initMatrix(DOFMatrix *matrix)
{
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
-
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if (!(*it).second->isDirichlet())
- (*it).second->initMatrix(matrix);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && !(*it).second->isDirichlet())
+ (*it).second->initMatrix(matrix);
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if ((*it).second->isDirichlet())
- (*it).second->initMatrix(matrix);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && (*it).second->isDirichlet())
+ (*it).second->initMatrix(matrix);
}
void BoundaryManager::exitMatrix(DOFMatrix *matrix)
{
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if (!(*it).second->isDirichlet())
- (*it).second->exitMatrix(matrix);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && !(*it).second->isDirichlet())
+ (*it).second->exitMatrix(matrix);
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if ((*it).second->isDirichlet())
- (*it).second->exitMatrix(matrix);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && (*it).second->isDirichlet())
+ (*it).second->exitMatrix(matrix);
}
void BoundaryManager::initVector(DOFVectorBase<double> *vector)
{
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if (!(*it).second->isDirichlet())
- (*it).second->initVector(vector);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && !(*it).second->isDirichlet())
+ (*it).second->initVector(vector);
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if ((*it).second->isDirichlet())
- (*it).second->initVector(vector);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && (*it).second->isDirichlet())
+ (*it).second->initVector(vector);
}
void BoundaryManager::exitVector(DOFVectorBase<double> *vector)
{
- std::map<BoundaryType, BoundaryCondition*>::iterator it;
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if (!(*it).second->isDirichlet())
- (*it).second->exitVector(vector);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && !(*it).second->isDirichlet())
+ (*it).second->exitVector(vector);
- for (it = localBCs.begin(); it != localBCs.end(); ++it)
- if ((*it).second)
- if ((*it).second->isDirichlet())
- (*it).second->exitVector(vector);
+ for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+ if ((*it).second && (*it).second->isDirichlet())
+ (*it).second->exitVector(vector);
}
}
diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h
index a027d680..a68aa641 100644
--- a/AMDiS/src/BoundaryManager.h
+++ b/AMDiS/src/BoundaryManager.h
@@ -30,6 +30,8 @@
namespace AMDiS {
+ typedef std::map<BoundaryType, BoundaryCondition*> BoundaryIndexMap;
+
/**
* \ingroup Assembler
*
@@ -50,9 +52,15 @@ namespace AMDiS {
/// Adds a local boundary condition to the list of managed conditions.
void addBoundaryCondition(BoundaryCondition *localBC)
{
+ FUNCNAME("BoundaryManager::addBoundaryCondition()");
+
BoundaryType type = localBC->getBoundaryType();
- TEST_EXIT(localBCs[type] == NULL)("there is already a condition for this type\n");
+ TEST_EXIT(localBCs[type] == NULL)
+ ("There is already a condition for this type.\n");
localBCs[type] = localBC;
+
+ std::vector<BoundaryCondition*>& boundMap = globalBoundaryMap[type];
+ boundMap.push_back(localBC);
}
void initMatrix(DOFMatrix *matrix);
@@ -79,8 +87,7 @@ namespace AMDiS {
* Calls BoundaryCondition::boundResidual() for each boundary condition in
* \ref localBCs.
*/
- double boundResidual(ElInfo *elInfo,
- DOFMatrix *matrix,
+ double boundResidual(ElInfo *elInfo, DOFMatrix *matrix,
const DOFVectorBase<double> *dv);
inline BoundaryCondition *getBoundaryCondition(BoundaryType type)
@@ -88,19 +95,32 @@ namespace AMDiS {
return localBCs[type];
}
- const std::map<BoundaryType, BoundaryCondition*>& getBoundaryConditionMap()
+ const BoundaryIndexMap& getBoundaryConditionMap()
{
return localBCs;
}
- void setBoundaryConditionMap(const std::map<BoundaryType, BoundaryCondition*>& bcs)
+ void setBoundaryConditionMap(const BoundaryIndexMap& bcs)
{
localBCs = bcs;
}
+
+ /** \brief
+ * Returns true, if there is at least one boundary object with the given boundary
+ * id, which implements a periodic boundary.
+ */
+ static bool isBoundaryPeriodic(BoundaryType b)
+ {
+ for (int i = 0; i < static_cast<int>(globalBoundaryMap[b].size()); i++)
+ if (globalBoundaryMap[b][i]->isPeriodic())
+ return true;
+
+ return false;
+ }
protected:
/// Map of managed local boundary conditions.
- std::map<BoundaryType, BoundaryCondition*> localBCs;
+ BoundaryIndexMap localBCs;
/// Temporary thread-safe variable for functions fillBoundaryconditions.
std::vector<BoundaryType*> localBounds;
@@ -113,6 +133,21 @@ namespace AMDiS {
* each localBounds value. Is used to free the memory in the destructor.
*/
int allocatedMemoryLocalBounds;
+
+ /** \brief
+ * For every boundary id we store here all possible boundary object (although
+ * it's not clear if it is meaningful to have different boundary conditions on the
+ * same boundary id).
+ *
+ * We have to use this global variable, because the mesh traverse interface does
+ * not provide more information about traversed boundaries at elements than the
+ * boundary id.
+ *
+ * TODO: Change interface such that mesh traverse returns the boundary objects
+ * directly and we can remove this global variable. The biggest problem will be
+ * than serialization and deserialization of the mesh.
+ */
+ static std::map<BoundaryType, std::vector<BoundaryCondition*> > globalBoundaryMap;
};
}
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 7ea53621..d6606ee8 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -30,16 +30,18 @@ namespace AMDiS {
inserter(NULL)
{}
- DOFMatrix::DOFMatrix(const FiniteElemSpace* rowFESpace_,
- const FiniteElemSpace* colFESpace_,
- std::string name_)
- : rowFESpace(rowFESpace_),
- colFESpace(colFESpace_),
- name(name_),
+ DOFMatrix::DOFMatrix(const FiniteElemSpace* rowSpace,
+ const FiniteElemSpace* colSpace,
+ std::string n)
+ : rowFESpace(rowSpace),
+ colFESpace(colSpace),
+ name(n),
coupleMatrix(false),
inserter(NULL)
{
- TEST_EXIT(rowFESpace)("no rowFESpace\n");
+ FUNCNAME("DOFMatrix::DOFMatrix()");
+
+ TEST_EXIT(rowFESpace)("No fe space for row!\n");
if (!colFESpace)
colFESpace = rowFESpace;
@@ -47,7 +49,7 @@ namespace AMDiS {
if (rowFESpace && rowFESpace->getAdmin())
(const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
- boundaryManager = new BoundaryManager(rowFESpace_);
+ boundaryManager = new BoundaryManager(rowFESpace);
nRow = rowFESpace->getBasisFcts()->getNumber();
nCol = colFESpace->getBasisFcts()->getNumber();
@@ -61,12 +63,14 @@ namespace AMDiS {
DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
: name(rhs.name + "copy")
{
+ FUNCNAME("DOFMatrix::DOFMatrix()");
+
*this = rhs;
if (rowFESpace && rowFESpace->getAdmin())
(const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
- TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion");
- inserter= 0;
+ TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
+ inserter = 0;
}
DOFMatrix::~DOFMatrix()
@@ -87,26 +91,6 @@ namespace AMDiS {
if (inserter)
inserter->print();
-
- /* using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
- namespace traits= mtl::traits;
- typedef base_matrix_type Matrix;
-
- traits::row<Matrix>::type row(matrix);
- traits::col<Matrix>::type col(matrix);
- traits::const_value<Matrix>::type value(matrix);
-
- typedef traits::range_generator<major, Matrix>::type cursor_type;
- typedef traits::range_generator<nz, cursor_type>::type icursor_type;
-
- std::cout.precision(10);
- for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) {
- for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
- if (value(*icursor) != 0.0)
- std::cout << "(" << row(*icursor) << "," << col(*icursor) << "," << value(*icursor) << ") ";
-
- std::cout << "\n";
- }*/
}
bool DOFMatrix::symmetric()
@@ -141,9 +125,9 @@ namespace AMDiS {
int non_symmetric = !symmetric();
if (non_symmetric)
- MSG("matrix `%s' not symmetric.\n", name.data());
+ MSG("Matrix `%s' not symmetric.\n", name.data());
else
- MSG("matrix `%s' is symmetric.\n", name.data());
+ MSG("Matrix `%s' is symmetric.\n", name.data());
}
DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index c2eb996c..fb0740e0 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -714,11 +714,13 @@ namespace AMDiS {
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof)
{
+ /*
std::vector<DegreeOfFreedom>::iterator it;
std::vector<DegreeOfFreedom>::iterator end = vec.end();
DegreeOfFreedom pos = 0;
for (it = vec.begin(); it != end; ++it, ++pos)
if (*it == dof) *it = pos;
+ */
}
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
diff --git a/AMDiS/src/DataCollector.cc b/AMDiS/src/DataCollector.cc
index c158402f..4e023542 100644
--- a/AMDiS/src/DataCollector.cc
+++ b/AMDiS/src/DataCollector.cc
@@ -171,9 +171,8 @@ namespace AMDiS {
getElementData()->
getElementData(PERIODIC));
- if (ldp) {
+ if (ldp)
nConnection += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
- }
periodicConnections.push_back(DimVec<bool>(dim, DEFAULT_VALUE, false));
}
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index a486147b..8563f4e1 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -276,11 +276,10 @@ namespace AMDiS {
fillFlag_.isSet(Mesh::FILL_DET) ||
fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- if (elem->isNewCoordSet()) {
+ if (elem->isNewCoordSet())
coord_[2] = *(elem->getNewCoord());
- } else {
- coord_[2].setMidpoint(elInfoOld->coord_[0], elInfoOld->coord_[1]);
- }
+ else
+ coord_[2].setMidpoint(elInfoOld->coord_[0], elInfoOld->coord_[1]);
if (ichild == 0) {
coord_[0] = elInfoOld->coord_[2];
diff --git a/AMDiS/src/InteriorBoundary.cc b/AMDiS/src/InteriorBoundary.cc
index 3d219649..bf90070e 100644
--- a/AMDiS/src/InteriorBoundary.cc
+++ b/AMDiS/src/InteriorBoundary.cc
@@ -3,7 +3,7 @@
namespace AMDiS {
- AtomicBoundary& InteriorBoundary::getNewAtomicBoundary(int rank)
+ AtomicBoundary& InteriorBoundary::getNewAtomic(int rank)
{
boundary[rank].resize(boundary[rank].size() + 1);
return boundary[rank][boundary[rank].size() - 1];
diff --git a/AMDiS/src/InteriorBoundary.h b/AMDiS/src/InteriorBoundary.h
index 148ccb34..79a48fa4 100644
--- a/AMDiS/src/InteriorBoundary.h
+++ b/AMDiS/src/InteriorBoundary.h
@@ -80,7 +80,7 @@ namespace AMDiS {
public:
InteriorBoundary() {}
- AtomicBoundary& getNewAtomicBoundary(int rank);
+ AtomicBoundary& getNewAtomic(int rank);
/// Writes this object to a file.
void serialize(std::ostream &out);
diff --git a/AMDiS/src/LeafData.h b/AMDiS/src/LeafData.h
index 84127837..4966318f 100644
--- a/AMDiS/src/LeafData.h
+++ b/AMDiS/src/LeafData.h
@@ -407,8 +407,9 @@ namespace AMDiS {
{
if (type == PERIODIC)
return true;
+
return false;
- };
+ }
class PeriodicInfo : public Serializable
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index c5e64c2e..9e58ec7a 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -141,15 +141,8 @@ namespace AMDiS {
element2->setElementData(ldp2);
}
- ldp1->addPeriodicInfo(mode,
- boundaryType,
- sideEl1,
- &periodicCoordsEl1);
-
- ldp2->addPeriodicInfo(mode,
- boundaryType,
- sideEl2,
- &periodicCoordsEl2);
+ ldp1->addPeriodicInfo(mode, boundaryType, sideEl1, &periodicCoordsEl1);
+ ldp2->addPeriodicInfo(mode, boundaryType, sideEl2, &periodicCoordsEl2);
if (mode != 0) {
VertexVector *associated = mesh->periodicAssociations[boundaryType];
@@ -184,9 +177,7 @@ namespace AMDiS {
std::map<BoundaryType, VertexVector*>::iterator assocEnd =
mesh->periodicAssociations.end();
- for (assoc = mesh->periodicAssociations.begin();
- assoc != assocEnd;
- ++assoc) {
+ for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
DegreeOfFreedom a = (*(assoc->second))[i];
if (a != i) {
@@ -201,9 +192,7 @@ namespace AMDiS {
std::map<BoundaryType, VertexVector*>::iterator assoc;
std::map<BoundaryType, VertexVector*>::iterator assocEnd =
mesh->periodicAssociations.end();
- for (assoc = mesh->periodicAssociations.begin();
- assoc != assocEnd;
- ++assoc) {
+ for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
if (i != (*(assoc->second))[i])
diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc
index 4733416c..806b26df 100644
--- a/AMDiS/src/ParallelDomainBase.cc
+++ b/AMDiS/src/ParallelDomainBase.cc
@@ -16,6 +16,7 @@
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
#include "ElementFileWriter.h"
+#include "VertexVector.h"
#include "petscksp.h"
@@ -118,11 +119,13 @@ namespace AMDiS {
DbgTestInteriorBoundary();
#endif
-
// === Reset all DOFAdmins of the mesh. ===
updateDofAdmins();
+ // === ===
+
+ createPeriodicMap();
// === Global refinements. ===
@@ -148,11 +151,44 @@ namespace AMDiS {
#if (DEBUG != 0)
- DbgTestCommonDofs(true);
+ // DbgTestCommonDofs(true);
#endif
nRankRows = nRankDofs * nComponents;
nOverallRows = nOverallDOFs * nComponents;
+
+#if 0
+ DOFVector<double> *tmp = new DOFVector<double>(feSpace, "tmp");
+ tmp->set(0.0);
+ if (mpiRank == 0)
+ VtkWriter::writeFile(tmp, "test0.vtu");
+ else
+ VtkWriter::writeFile(tmp, "test1.vtu");
+
+ if (mpiRank == 1) {
+ for (PeriodicDofMap::iterator it = periodicDof.begin();
+ it != periodicDof.end(); ++it) {
+ std::cout << "DOF MAP " << it->first << ": ";
+ for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
+ dofit != it->second.end(); ++dofit)
+ std::cout << *dofit << " ";
+ std::cout << std::endl;
+
+ DegreeOfFreedom localdof;
+ for (DofMapping::iterator dofIt = mapLocalGlobalDOFs.begin();
+ dofIt != mapLocalGlobalDOFs.end(); ++dofIt)
+ if (dofIt->second == it->first)
+ localdof = dofIt->first;
+
+ WorldVector<double> coords;
+ mesh->getDofIndexCoords(localdof, feSpace, coords);
+ coords.print();
+ }
+ }
+#endif
+
+ if (mpiRank == 0)
+ writePartitioningMesh("part.vtu");
}
@@ -188,7 +224,7 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
TEST_EXIT(elInfo->getLevel() == 0)
- ("Mesh is already refined! This does nor work with parallelization!\n");
+ ("Mesh is already refined! This does not work with parallelization!\n");
nMacroElements++;
@@ -227,16 +263,57 @@ namespace AMDiS {
cols.clear();
values.clear();
- int r = mapLocalGlobalDOFs[*cursor] * dispMult + dispAddRow;
+ int globalRowDof = mapLocalGlobalDOFs[*cursor];
+ bool periodicRow = (periodicDof.count(globalRowDof) > 0);
- for (icursor_type icursor = begin<nz>(cursor),
- icend = end<nz>(cursor); icursor != icend; ++icursor)
+ for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
+ icursor != icend; ++icursor) {
if (value(*icursor) != 0.0) {
- cols.push_back(mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol);
- values.push_back(value(*icursor));
+ int globalColDof = mapLocalGlobalDOFs[col(*icursor)];
+ int colIndex = globalColDof * dispMult + dispAddCol;
+
+ if (periodicRow == false && periodicDof.count(globalColDof) > 0) {
+ cols.push_back(colIndex);
+ values.push_back(value(*icursor) * 0.5);
+
+ for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalColDof].begin();
+ it != periodicDof[globalColDof].end(); ++it) {
+ cols.push_back(*it * dispMult + dispAddCol);
+ values.push_back(value(*icursor) * 0.5);
+ }
+ } else {
+ cols.push_back(colIndex);
+ values.push_back(value(*icursor));
+ }
+ }
+ }
+
+ int rowIndex = globalRowDof * dispMult + dispAddRow;
+
+ if (periodicRow) {
+ int diagIndex = -1;
+
+ for (int i = 0; i < static_cast<int>(values.size()); i++) {
+ if (cols[i] != rowIndex)
+ values[i] *= 0.5;
+ else
+ diagIndex = i;
+ }
+
+ MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
+
+ if (diagIndex != -1)
+ values[diagIndex] = 0.0;
+
+ for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRowDof].begin();
+ it != periodicDof[globalRowDof].end(); ++it) {
+ int perRowIndex = *it * dispMult + dispAddRow;
+ MatSetValues(petscMatrix, 1, &perRowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
}
- MatSetValues(petscMatrix, 1, &r, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
+ } else {
+ MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
+ }
}
}
@@ -246,10 +323,23 @@ namespace AMDiS {
{
DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
- int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()] * dispMult + dispAdd;
- double value = *dofIt;
+ int globalRow = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
+ int index = globalRow * dispMult + dispAdd;
+
+ if (periodicDof.count(globalRow) > 0) {
+ double value = *dofIt * 0.5;
+ VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
+
+ for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRow].begin();
+ it != periodicDof[globalRow].end(); ++it) {
+ index = *it * dispMult + dispAdd;
+ VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
+ }
- VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
+ } else {
+ double value = *dofIt;
+ VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
+ }
}
}
@@ -476,8 +566,7 @@ namespace AMDiS {
TEST_EXIT_DBG(r >= 0 && r < nRankRows)("Should not happen!\n");
- if (c < rstart * nComponents ||
- c >= rstart * nComponents + nRankRows)
+ if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
o_nnz[r]++;
else
d_nnz[r]++;
@@ -788,9 +877,9 @@ namespace AMDiS {
(element->getElementData(PARTITION_ED));
if (partitionData && partitionData->getPartitionStatus() == IN) {
- if (partitionData->getLevel() == 0) {
+ if (partitionData->getLevel() == 0)
elNum = element->getIndex();
- }
+
TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
if (element->isLeaf()) {
elemWeights[elNum] += 1.0;
@@ -824,11 +913,17 @@ namespace AMDiS {
{
FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
+ // To be general, we have to use elInfo->getBoundary(EDGE/FACE, i) instead of
+ // elInfo->getBoundar(i).
+ TEST_EXIT(mesh->getDim() == 2)
+ ("getBoundary(i) returns always i-th edge, generalize for 3d!\n");
+
// === First, create all the information about the interior boundaries. ===
TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
+ ElInfo *elInfo =
+ stack.traverseFirst(mesh, -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
while (elInfo) {
Element *element = elInfo->getElement();
@@ -844,16 +939,8 @@ namespace AMDiS {
dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED));
if (neighbourPartitionData->getPartitionStatus() == OUT) {
- // We have found an element that is at an interior boundary.
-
- int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];
- bool ranksBoundary = (mpiRank > otherElementRank);
-
- AtomicBoundary& bound =
- (ranksBoundary ?
- myIntBoundary.getNewAtomicBoundary(otherElementRank) :
- otherIntBoundary.getNewAtomicBoundary(otherElementRank));
+ AtomicBoundary bound;
bound.rankObject.el = element;
bound.rankObject.elIndex = element->getIndex();
bound.rankObject.subObjAtBoundary = EDGE;
@@ -864,10 +951,24 @@ namespace AMDiS {
bound.neighbourObject.elIndex = elInfo->getNeighbour(i)->getIndex();
bound.neighbourObject.subObjAtBoundary = EDGE;
bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);
-
+
// i == 2 => getSideOfNeighbour(i) == 2
TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
("Should not happen!\n");
+
+ int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];
+
+ if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(i))) {
+ // We have found an element that is at an interior, periodic boundary.
+ AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
+ b = bound;
+ } else {
+ // We have found an element that is at an interior, non-periodic boundary.
+ AtomicBoundary& b = ((mpiRank > otherElementRank) ?
+ myIntBoundary.getNewAtomic(otherElementRank) :
+ otherIntBoundary.getNewAtomic(otherElementRank));
+ b = bound;
+ }
}
}
}
@@ -875,15 +976,15 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
-
// === Once we have this information, we must care about their order. Eventually ===
// === all the boundaries have to be in the same order on both ranks that share ===
// === the bounday. ===
std::vector<int*> sendBuffers, recvBuffers;
- MPI::Request request[myIntBoundary.boundary.size() +
- otherIntBoundary.boundary.size()];
+ MPI::Request request[max(myIntBoundary.boundary.size() +
+ otherIntBoundary.boundary.size(),
+ periodicBoundary.boundary.size())];
int requestCounter = 0;
for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
@@ -943,11 +1044,84 @@ namespace AMDiS {
}
}
- delete [] recvBuffers[i++];
+ delete [] recvBuffers[i];
+ i++;
}
for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
delete [] sendBuffers[i];
+
+ recvBuffers.clear();
+ sendBuffers.clear();
+
+
+ // === Deal with periodic boundaries. ===
+
+ requestCounter = 0;
+
+ for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
+ rankIt != periodicBoundary.boundary.end(); ++rankIt) {
+
+ TEST_EXIT_DBG(rankIt->first != mpiRank)
+ ("It is no possible to have an interior boundary within a rank partition!\n");
+
+ if (rankIt->first < mpiRank) {
+ int nSendInt = rankIt->second.size();
+ int* buffer = new int[nSendInt * 2];
+ for (int i = 0; i < nSendInt; i++) {
+ buffer[i * 2] = (rankIt->second)[i].rankObject.elIndex;
+ buffer[i * 2 + 1] = (rankIt->second)[i].rankObject.ithObjAtBoundary;
+ }
+ sendBuffers.push_back(buffer);
+
+ request[requestCounter++] =
+ mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
+ } else {
+ int nRecvInt = rankIt->second.size();
+ int *buffer = new int[nRecvInt * 2];
+ recvBuffers.push_back(buffer);
+
+ request[requestCounter++] =
+ mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
+ }
+ }
+
+
+ MPI::Request::Waitall(requestCounter, request);
+
+
+ i = 0;
+ for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
+ rankIt != periodicBoundary.boundary.end(); ++rankIt) {
+ if (rankIt->first > mpiRank) {
+
+ for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++)
+ if (periodicBoundary.boundary[rankIt->first][j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
+ periodicBoundary.boundary[rankIt->first][j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
+ int k = j + 1;
+
+ for (; k < static_cast<int>(rankIt->second.size()); k++)
+ if (periodicBoundary.boundary[rankIt->first][k].neighbourObject.elIndex == recvBuffers[i][j * 2] &&
+ periodicBoundary.boundary[rankIt->first][k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
+ break;
+
+ // The element must always be found, because the list is just in another order.
+ TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
+ ("Should never happen!\n");
+
+ // Swap the current with the found element.
+ AtomicBoundary tmpBound = (rankIt->second)[k];
+ (rankIt->second)[k] = (rankIt->second)[j];
+ (rankIt->second)[j] = tmpBound;
+ }
+
+ delete [] recvBuffers[i];
+ i++;
+ }
+ }
+
+ for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
+ delete [] sendBuffers[i];
}
@@ -983,6 +1157,30 @@ namespace AMDiS {
createDofMemberInfo(partitionDOFs, rankDofs, rankAllDofs, boundaryDofs, vertexDof);
+ // === BEGIN EXP
+
+ DofIndexToBool rankAllDofIndices;
+ for (DofContainer::iterator dofIt = rankAllDofs.begin();
+ dofIt != rankAllDofs.end(); ++dofIt) {
+ rankAllDofIndices[**dofIt] = true;
+ }
+
+ for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
+ it != mesh->getPeriodicAssociations().end(); ++it) {
+ VertexVector *tmp = it->second;
+ DOFIteratorBase it(const_cast<DOFAdmin*>(tmp->getAdmin()), USED_DOFS);
+ for (it.reset(); !it.end(); ++it) {
+ if (!it.isDOFFree()) {
+ if (!rankAllDofIndices[(*tmp)[it.getDOFIndex()]] &&
+ rankAllDofIndices[it.getDOFIndex()])
+ // (*tmp)[it.getDOFIndex()] = -1;
+ (*tmp)[it.getDOFIndex()] = it.getDOFIndex();
+ }
+ }
+ }
+
+ // === ENDE EXP
+
nRankDofs = rankDofs.size();
nOverallDOFs = partitionDOFs.size();
@@ -1095,7 +1293,7 @@ namespace AMDiS {
int requestCounter = 0;
i = 0;
- for (std::map<int, DofIndexMap >::iterator sendIt = sendNewDofs.begin();
+ for (std::map<int, DofIndexMap>::iterator sendIt = sendNewDofs.begin();
sendIt != sendNewDofs.end(); ++sendIt, i++) {
int nSendDofs = sendIt->second.size() * 2;
sendBuffers[i] = new int[nSendDofs];
@@ -1182,6 +1380,23 @@ namespace AMDiS {
delete [] recvBuffers[i];
}
+ // ====== BEGIN TEST
+
+ std::map<int, int> dofIndexMap;
+ for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
+ dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
+ dofIndexMap[*(dofIt->first)] = dofIt->second;
+ }
+
+ for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
+ it != mesh->getPeriodicAssociations().end(); ++it) {
+ it->second->changeDofIndices(dofIndexMap);
+ }
+
+
+ // ====== ENDE TEST
+
+
// === Create now the local to global index and local to dof index mappings. ===
createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
@@ -1239,29 +1454,28 @@ namespace AMDiS {
sendDofs.clear();
recvDofs.clear();
- for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
- it != myIntBoundary.boundary.end(); ++it) {
+ for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
+ it != myIntBoundary.boundary.end(); ++it) {
+
+ DofContainer &dofsToSend = sendDofs[it->first];
+
+ for (DofContainer::iterator iit = oldSendDofs[it->first].begin();
+ iit != oldSendDofs[it->first].end(); ++iit)
+ if (oldVertexDof[*iit])
+ dofsToSend.push_back(*iit);
for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
boundIt != it->second.end(); ++boundIt) {
-
- DofContainer &dofsToSend = sendDofs[it->first];
-
- for (DofContainer::iterator iit = oldSendDofs[it->first].begin();
- iit != oldSendDofs[it->first].end(); ++iit)
- if (oldVertexDof[*iit])
- dofsToSend.push_back(*iit);
-
DofContainer dofs;
- addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
- dofs);
- addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
- dofs);
+ addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ dofs);
+ addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ dofs);
for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
TEST_EXIT_DBG(find(dofsToSend.begin(), dofsToSend.end(), dofs[i]) == dofsToSend.end())
("Should not happen!\n");
- dofsToSend.push_back(dofs[i]);
+ dofsToSend.push_back(dofs[i]);
}
}
}
@@ -1270,27 +1484,26 @@ namespace AMDiS {
for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
it != otherIntBoundary.boundary.end(); ++it) {
- for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
- boundIt != it->second.end(); ++boundIt) {
-
- DofContainer &dofsToRecv = recvDofs[it->first];
-
- for (DofContainer::iterator iit = oldRecvDofs[it->first].begin();
- iit != oldRecvDofs[it->first].end(); ++iit)
- if (oldVertexDof[*iit]) {
- dofsToRecv.push_back(*iit);
+ DofContainer &dofsToRecv = recvDofs[it->first];
+
+ for (DofContainer::iterator iit = oldRecvDofs[it->first].begin();
+ iit != oldRecvDofs[it->first].end(); ++iit)
+ if (oldVertexDof[*iit]) {
+ dofsToRecv.push_back(*iit);
- DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), *iit);
- if (eraseIt != rankDofs.end())
- rankDofs.erase(eraseIt);
- }
+ DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), *iit);
+ if (eraseIt != rankDofs.end())
+ rankDofs.erase(eraseIt);
+ }
+ for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+ boundIt != it->second.end(); ++boundIt) {
DofContainer dofs;
- addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
- dofs);
- addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
- dofs);
+ addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ dofs);
+ addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ dofs);
for (int i = static_cast<int>(dofs.size()) - 1; i >= 0; i--) {
TEST_EXIT_DBG(find(dofsToRecv.begin(), dofsToRecv.end(), dofs[i]) == dofsToRecv.end())
@@ -1398,7 +1611,6 @@ namespace AMDiS {
delete [] recvBuffers[i++];
}
-
// === Create now the local to global index and local to dof index mappings. ===
createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
@@ -1426,31 +1638,52 @@ namespace AMDiS {
mapLocalToDofIndex[dofIt->second] = *(dofIt->first);
}
- void ParallelDomainBase::addAllVertexDOFs(Element *el, int ithEdge,
- DofContainer& dofs)
+ void ParallelDomainBase::addVertexDofs(Element *el, int ithEdge, DofContainer& dofs,
+ bool parentVertices)
{
- FUNCNAME("ParallelDomainBase::addAllVertexDOFs()");
+ FUNCNAME("ParallelDomainBase::addVertexDofs()");
+
+ TEST_EXIT_DBG(mesh->getDim() == 2)("Not yet implemented for 3D!\n");
+
+ if (parentVertices) {
+ switch (ithEdge) {
+ case 0:
+ dofs.push_back(el->getDOF(1));
+ dofs.push_back(el->getDOF(2));
+ break;
+ case 1:
+ dofs.push_back(el->getDOF(0));
+ dofs.push_back(el->getDOF(2));
+ break;
+ case 2:
+ dofs.push_back(el->getDOF(0));
+ dofs.push_back(el->getDOF(1));
+ break;
+ default:
+ ERROR_EXIT("Should never happen!\n");
+ }
+ }
switch (ithEdge) {
case 0:
if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
- addAllVertexDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
+ addVertexDofs(el->getSecondChild()->getFirstChild(), 0, dofs);
dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
- addAllVertexDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
+ addVertexDofs(el->getSecondChild()->getSecondChild(), 1, dofs);
}
break;
case 1:
if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
- addAllVertexDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
+ addVertexDofs(el->getFirstChild()->getFirstChild(), 0, dofs);
dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
- addAllVertexDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
+ addVertexDofs(el->getFirstChild()->getSecondChild(), 1, dofs);
}
break;
case 2:
if (el->getFirstChild()) {
- addAllVertexDOFs(el->getFirstChild(), 0, dofs);
+ addVertexDofs(el->getFirstChild(), 0, dofs);
dofs.push_back(el->getFirstChild()->getDOF(2));
- addAllVertexDOFs(el->getSecondChild(), 1, dofs);
+ addVertexDofs(el->getSecondChild(), 1, dofs);
}
break;
default:
@@ -1459,32 +1692,33 @@ namespace AMDiS {
}
- void ParallelDomainBase::addAllEdgeDOFs(Element *el, int ithEdge,
- DofContainer& dofs)
+ void ParallelDomainBase::addEdgeDofs(Element *el, int ithEdge, DofContainer& dofs)
{
- FUNCNAME("ParallelDomainBase::addAllEdgeDOFs()");
+ FUNCNAME("ParallelDomainBase::addEdgeDofs()");
+
+ TEST_EXIT_DBG(mesh->getDim() == 2)("Not yet implemented for 3D!\n");
bool addThisEdge = false;
switch (ithEdge) {
case 0:
if (el->getSecondChild())
- addAllEdgeDOFs(el->getSecondChild(), 2, dofs);
+ addEdgeDofs(el->getSecondChild(), 2, dofs);
else
addThisEdge = true;
break;
case 1:
if (el->getFirstChild())
- addAllEdgeDOFs(el->getFirstChild(), 2, dofs);
+ addEdgeDofs(el->getFirstChild(), 2, dofs);
else
addThisEdge = true;
break;
case 2:
if (el->getFirstChild()) {
- addAllEdgeDOFs(el->getFirstChild(), 0, dofs);
- addAllEdgeDOFs(el->getSecondChild(), 1, dofs);
+ addEdgeDofs(el->getFirstChild(), 0, dofs);
+ addEdgeDofs(el->getSecondChild(), 1, dofs);
} else {
addThisEdge = true;
}
@@ -1594,6 +1828,76 @@ namespace AMDiS {
}
+ void ParallelDomainBase::createPeriodicMap()
+ {
+ FUNCNAME("ParallelDomainBase::createPeriodicMap()");
+
+ MPI::Request request[periodicBoundary.boundary.size() * 2];
+ int requestCounter = 0;
+ std::vector<int*> sendBuffers, recvBuffers;
+
+ for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
+ it != periodicBoundary.boundary.end(); ++it) {
+
+ TEST_EXIT_DBG(it->first != mpiRank)
+ ("Periodic interior boundary within the rank itself is not possible!\n");
+
+ DofContainer dofs;
+
+ for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+ boundIt != it->second.end(); ++boundIt) {
+ addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ dofs, true);
+ addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ dofs);
+ }
+
+ int *sendbuf = new int[dofs.size()];
+ for (int i = 0; i < static_cast<int>(dofs.size()); i++)
+ sendbuf[i] = mapLocalGlobalDOFs[*(dofs[i])];
+
+ request[requestCounter++] =
+ mpiComm.Isend(sendbuf, dofs.size(), MPI_INT, it->first, 0);
+
+ sendBuffers.push_back(sendbuf);
+
+ int *recvbuf = new int[dofs.size()];
+ request[requestCounter++] =
+ mpiComm.Irecv(recvbuf, dofs.size(), MPI_INT, it->first, 0);
+
+ recvBuffers.push_back(recvbuf);
+ }
+
+
+ MPI::Request::Waitall(requestCounter, request);
+
+
+ int i = 0;
+ for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
+ it != periodicBoundary.boundary.end(); ++it) {
+ DofContainer dofs;
+
+ for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+ boundIt != it->second.end(); ++boundIt) {
+ DofContainer tmpdofs;
+ addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ tmpdofs);
+ addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+ tmpdofs, true);
+ for (int j = static_cast<int>(tmpdofs.size()) - 1; j >= 0; j--)
+ dofs.push_back(tmpdofs[j]);
+ }
+
+ for (int j = 0; j < static_cast<int>(dofs.size()); j++)
+ periodicDof[mapLocalGlobalDOFs[*(dofs[j])]].insert(recvBuffers[i][j]);
+
+ delete [] sendBuffers[i];
+ delete [] recvBuffers[i];
+ i++;
+ }
+ }
+
+
Flag ParallelDomainBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
{
FUNCNAME("ParallelDomainBase::oneIteration()");
@@ -1799,7 +2103,7 @@ namespace AMDiS {
{
FUNCNAME("ParallelDomainBase::DbgTestCommonDofs()");
- // Maps to each neighbour rank an array of WorldVectors. This array contain the
+ // 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
// neighbour rank. A rank sends the coordinates to another rank, if it owns the
// boundarys dof.
diff --git a/AMDiS/src/ParallelDomainBase.h b/AMDiS/src/ParallelDomainBase.h
index 6e9ff028..34fbca52 100644
--- a/AMDiS/src/ParallelDomainBase.h
+++ b/AMDiS/src/ParallelDomainBase.h
@@ -81,6 +81,8 @@ namespace AMDiS {
typedef std::map<const DegreeOfFreedom*, DegreeOfFreedom> DofIndexMap;
+ typedef std::map<DegreeOfFreedom, std::set<DegreeOfFreedom> > PeriodicDofMap;
+
public:
ParallelDomainBase(ProblemIterationInterface *iterationIF,
ProblemTimeInterface *timeIF,
@@ -93,7 +95,7 @@ namespace AMDiS {
void exitParallelization(AdaptInfo *adaptInfo);
- void updateDofAdmins();
+ void updateDofAdmins();
/** \brief
* Test, if the mesh consists of macro elements only. The mesh partitioning of
@@ -212,6 +214,8 @@ namespace AMDiS {
void updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs);
+ void createPeriodicMap();
+
/** \brief
* This function create new mappings from local to global indices,
* \ref mapLocalGlobalDOFs, and from local to dof indices, \ref mapLocalToDofIndex.
@@ -229,9 +233,30 @@ namespace AMDiS {
DofIndexMap &rankOwnedDofsNewLocalIndex,
DofIndexMap &rankDofsNewGlobalIndex);
- void addAllVertexDOFs(Element *el, int ithEdge, DofContainer& dofs);
+ /** \brief
+ * Traverses an edge on a given element and all its child elements. All vertex
+ * dofs alonge this edge are assembled and put together to a list.
+ *
+ * \param[in] el Element for the edge traverse.
+ * \param[in] ithEdge Defines the edge on which all the vertex dofs
+ * are assembled.
+ * \param[out] dofs All dofs are put to this dof list.
+ * \param[in] parentVertices If true, also the two vertices of the parent
+ * element el are put into the result list.
+ */
+ void addVertexDofs(Element *el, int ithEdge, DofContainer& dofs,
+ bool parentVertices = false);
- void addAllEdgeDOFs(Element *el, int ithEdge, DofContainer& dofs);
+ /** \brief
+ * Traverses an edge on a given element and all its child elements. All non vertex
+ * dofs alonge this edge are assembled and put together to a list.
+ *
+ * \param[in] el Element for the edge traverse.
+ * \param[in] ithEdge Defines the edge on which all the non vertex
+ * dofs are assembled.
+ * \param[out] dofs All dofs are put to this dof list.
+ */
+ void addEdgeDofs(Element *el, int ithEdge, DofContainer& dofs);
/** \brief
* This function traverses the whole mesh, i.e. before it is really partitioned,
@@ -456,6 +481,11 @@ namespace AMDiS {
*/
InteriorBoundary otherIntBoundary;
+ /** \brief
+ *
+ */
+ InteriorBoundary periodicBoundary;
+
/** \brief
* This map contains for each rank the list of dofs the current rank must send
* to exchange solution dofs at the interior boundaries.
@@ -487,6 +517,8 @@ namespace AMDiS {
*/
DofToBool vertexDof;
+ PeriodicDofMap periodicDof;
+
/// Is the index of the first row of the linear system, which is owned by the rank.
int rstart;
diff --git a/AMDiS/src/ParallelDomainVec.cc b/AMDiS/src/ParallelDomainVec.cc
index 12edf33b..6fb80f1a 100644
--- a/AMDiS/src/ParallelDomainVec.cc
+++ b/AMDiS/src/ParallelDomainVec.cc
@@ -4,6 +4,7 @@
#include "ProblemInstat.h"
#include "SystemVector.h"
#include "Serializer.h"
+#include "BoundaryManager.h"
namespace AMDiS {
@@ -58,7 +59,7 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++)
if (probVec->getSystemMatrix(i, j))
- probVec->getSystemMatrix(i, j)->setRankDofs(isRankDof);
+ probVec->getSystemMatrix(i, j)->setRankDofs(isRankDof);
TEST_EXIT_DBG(probVec->getRHS()->getDOFVector(i))("No rhs vector!\n");
TEST_EXIT_DBG(probVec->getSolution()->getDOFVector(i))("No solution vector!\n");
@@ -66,6 +67,23 @@ namespace AMDiS {
probVec->getRHS()->getDOFVector(i)->setRankDofs(isRankDof);
probVec->getSolution()->getDOFVector(i)->setRankDofs(isRankDof);
}
+
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++)
+ if (probVec->getSystemMatrix(i, j) && probVec->getSystemMatrix(i, j)->getBoundaryManager())
+ removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getSystemMatrix(i, j)->getBoundaryManager()->getBoundaryConditionMap()));
+
+ if (probVec->getSolution()->getDOFVector(i)->getBoundaryManager())
+ removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
+ }
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
+ while (elInfo) {
+ elInfo->getElement()->deleteElementData(PERIODIC);
+ elInfo = stack.traverseNext(elInfo);
+ }
+
}
@@ -91,4 +109,16 @@ namespace AMDiS {
#endif
}
+
+ void ParallelDomainVec::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
+ {
+ BoundaryIndexMap::iterator it = boundaryMap.begin();
+ while (it != boundaryMap.end()) {
+ if (it->second->isPeriodic())
+ boundaryMap.erase(it++);
+ else
+ ++it;
+ }
+ }
+
}
diff --git a/AMDiS/src/ParallelDomainVec.h b/AMDiS/src/ParallelDomainVec.h
index 4d6a7f4f..ef79ce30 100644
--- a/AMDiS/src/ParallelDomainVec.h
+++ b/AMDiS/src/ParallelDomainVec.h
@@ -24,6 +24,7 @@
#include "ParallelDomainBase.h"
#include "ProblemVec.h"
+#include "BoundaryManager.h"
namespace AMDiS {
@@ -60,6 +61,8 @@ namespace AMDiS {
/// Starts the solution of the linear system using Petsc.
void solve();
+ void removeBoundaryCondition(BoundaryIndexMap& boundaryMap);
+
protected:
/// Pointer to the stationary problem.
ProblemVec *probVec;
diff --git a/AMDiS/src/PeriodicBC.cc b/AMDiS/src/PeriodicBC.cc
index a252749d..6e178db4 100644
--- a/AMDiS/src/PeriodicBC.cc
+++ b/AMDiS/src/PeriodicBC.cc
@@ -9,40 +9,39 @@
namespace AMDiS {
- std::vector<PeriodicDOFMapping*> PeriodicDOFMapping::mappings_;
+ std::vector<PeriodicDOFMapping*> PeriodicDOFMapping::mappings;
PeriodicDOFMapping*
- PeriodicDOFMapping::providePeriodicDOFMapping(const BasisFunction *basFcts)
+ PeriodicDOFMapping::providePeriodicDOFMapping(const BasisFunction *fcts)
{
- std::vector<PeriodicDOFMapping*>::iterator it;
- std::vector<PeriodicDOFMapping*>::iterator end = mappings_.end();
- for (it = mappings_.begin(); it != end; ++it)
- if ((*it)->basFcts_ == basFcts)
+ for (std::vector<PeriodicDOFMapping*>::iterator it = mappings.begin();
+ it != mappings.end(); ++it)
+ if ((*it)->basFcts == fcts)
return *it;
- PeriodicDOFMapping *newMapping = new PeriodicDOFMapping(basFcts);
- mappings_.push_back(newMapping);
+ PeriodicDOFMapping *newMapping = new PeriodicDOFMapping(fcts);
+ mappings.push_back(newMapping);
return newMapping;
}
- PeriodicDOFMapping::PeriodicDOFMapping(const BasisFunction *basFcts)
- : basFcts_(basFcts)
+ PeriodicDOFMapping::PeriodicDOFMapping(const BasisFunction *fcts)
+ : basFcts(fcts)
{
FUNCNAME("PeriodicDOFMapping::PeriodicDOFMapping()");
- TEST_EXIT(basFcts_->getDim() > 1)("dim == 1\n");
- int num = basFcts_->getNumber();
+ TEST_EXIT(basFcts->getDim() > 1)("dim == 1\n");
+ int nBasFcts = basFcts->getNumber();
DimVec<double> *lambda;
- for (int i = 0; i < num; i++) {
- lambda = basFcts_->getCoords(i);
- indexOfCoords_[*lambda] = i;
+ for (int i = 0; i < nBasFcts; i++) {
+ lambda = basFcts->getCoords(i);
+ indexOfCoords[*lambda] = i;
}
}
PeriodicDOFMapping::~PeriodicDOFMapping()
{
std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> >::iterator it;
- for (it = dofPermutation_.begin(); it != dofPermutation_.end(); ++it)
+ for (it = dofPermutation.begin(); it != dofPermutation.end(); ++it)
if (it->second)
delete [] it->second;
}
@@ -51,11 +50,9 @@ namespace AMDiS {
{
FUNCNAME("PeriodicDOFMapping::getDOFPermutation()");
- if (dofPermutation_[vertexPermutation] == NULL) {
- //MSG("new dof permutation needed\n");
-
- int dim = basFcts_->getDim();
- int num = basFcts_->getNumber();
+ if (dofPermutation[vertexPermutation] == NULL) {
+ int dim = basFcts->getDim();
+ int nBasFcts = basFcts->getNumber();
int sum = 0;
for (int i = 0; i < dim + 1; i++) {
sum += i - vertexPermutation[i];
@@ -68,33 +65,32 @@ namespace AMDiS {
DimVec<double> *lambda;
DimVec<double> newLambda(dim, NO_INIT);
- DegreeOfFreedom *mapping = new DegreeOfFreedom[num];
+ DegreeOfFreedom *mapping = new DegreeOfFreedom[nBasFcts];
- for (int i = 0; i < num; i++) {
- lambda = basFcts_->getCoords(i);
+ for (int i = 0; i < nBasFcts; i++) {
+ lambda = basFcts->getCoords(i);
for (int j = 0; j < dim + 1; j++)
newLambda[vertexPermutation[j]] = (*lambda)[j];
- mapping[i] = indexOfCoords_[newLambda];
+ mapping[i] = indexOfCoords[newLambda];
}
- dofPermutation_[vertexPermutation] = mapping;
+ dofPermutation[vertexPermutation] = mapping;
}
- return dofPermutation_[vertexPermutation];
+ return dofPermutation[vertexPermutation];
}
- PeriodicBC::PeriodicBC(BoundaryType type, FiniteElemSpace *rowFESpace_)
- : BoundaryCondition(type, rowFESpace_, NULL),
- masterMatrix_(NULL)
+ PeriodicBC::PeriodicBC(BoundaryType type, FiniteElemSpace *rowSpace)
+ : BoundaryCondition(type, rowSpace, NULL),
+ masterMatrix(NULL)
{
- if (rowFESpace->getMesh()->getDim() > 1) {
- periodicDOFMapping_ =
- PeriodicDOFMapping::providePeriodicDOFMapping(rowFESpace_->getBasisFcts());
- } else {
- periodicDOFMapping_ = NULL;
- }
+ if (rowFESpace->getMesh()->getDim() > 1)
+ periodicDOFMapping =
+ PeriodicDOFMapping::providePeriodicDOFMapping(rowFESpace->getBasisFcts());
+ else
+ periodicDOFMapping = NULL;
}
PeriodicBC::~PeriodicBC()
@@ -104,18 +100,16 @@ namespace AMDiS {
{
FUNCNAME("PeriodicBC::initMatrix()");
- if (!masterMatrix_) {
- masterMatrix_ = matrix;
-
+ if (!masterMatrix) {
+ masterMatrix = matrix;
Mesh *mesh = matrix->getRowFESpace()->getMesh();
-
associated = mesh->getPeriodicAssociations()[boundaryType];
TEST_EXIT(associated)("No associations for periodic boundary condition %d!\n",
boundaryType);
const BasisFunction *basFcts = rowFESpace->getBasisFcts();
- neighIndices_ = new DegreeOfFreedom[basFcts->getNumber()];
+ neighIndices = new DegreeOfFreedom[basFcts->getNumber()];
}
}
@@ -125,40 +119,31 @@ namespace AMDiS {
const BoundaryType *localBound,
int nBasFcts)
{
- if (matrix == masterMatrix_) {
+ if (matrix == masterMatrix) {
int dim = rowFESpace->getMesh()->getDim();
if (dim > 1) {
DOFAdmin *admin = rowFESpace->getAdmin();
-
FixVec<int, WORLD> elFace(dim, NO_INIT);
FixVec<int, WORLD> neighFace(dim, NO_INIT);
DimVec<int> vertexPermutation(dim, NO_INIT);
-
const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int num = basFcts->getNumber();
-
Element *element = elInfo->getElement();
+ DimVec<DegreeOfFreedom> periodicDOFs(dim - 1, NO_INIT);
+ GeoIndex sideGeoIndex = INDEX_OF_DIM(dim - 1, dim);
- DimVec<DegreeOfFreedom> periodicDOFs(dim-1, NO_INIT);
-
- int vertex, index, side;
-
- GeoIndex sideGeoIndex = INDEX_OF_DIM(dim-1, dim);
-
- for (side = 0; side < dim + 1; side++) {
+ for (int side = 0; side < dim + 1; side++) {
if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
- for (vertex = 0; vertex < dim; vertex++) {
- index = element->getVertexOfPosition(sideGeoIndex,
- side,
- vertex);
+ for (int vertex = 0; vertex < dim; vertex++) {
+ int index = element->getVertexOfPosition(sideGeoIndex, side, vertex);
periodicDOFs[vertex] = (*associated)[dofIndices[index]];
}
Element *neigh = elInfo->getNeighbour(side);
- basFcts->getLocalIndices(neigh, admin, neighIndices_);
+ basFcts->getLocalIndices(neigh, admin, neighIndices);
int oppVertex = 0;
for (int i = 0; i < dim + 1; i++) {
@@ -182,12 +167,12 @@ namespace AMDiS {
// get DOF permutation
const DegreeOfFreedom *dofPermutation =
- periodicDOFMapping_->getDOFPermutation(vertexPermutation);
+ periodicDOFMapping->getDOFPermutation(vertexPermutation);
// set associated dofs
for (int i = 0; i < num; i++)
if ((*(basFcts->getCoords(i)))[side] == 0)
- (*associated)[dofIndices[i]] = neighIndices_[dofPermutation[i]];
+ (*associated)[dofIndices[i]] = neighIndices[dofPermutation[i]];
}
}
}
@@ -200,9 +185,9 @@ namespace AMDiS {
TEST_EXIT(matrix)("no matrix\n");
- if (matrix == masterMatrix_) {
- delete [] neighIndices_;
- masterMatrix_ = NULL;
+ if (matrix == masterMatrix) {
+ delete [] neighIndices;
+ masterMatrix = NULL;
}
using namespace mtl;
@@ -236,7 +221,6 @@ namespace AMDiS {
for (vecIt.reset(); !vecIt.end(); ++vecIt) {
dof = vecIt.getDOFIndex();
-
newDOF = (*associated)[dof];
if (dof < newDOF) {
diff --git a/AMDiS/src/PeriodicBC.h b/AMDiS/src/PeriodicBC.h
index f0444f16..c58cb8ba 100644
--- a/AMDiS/src/PeriodicBC.h
+++ b/AMDiS/src/PeriodicBC.h
@@ -60,13 +60,13 @@ namespace AMDiS {
const DegreeOfFreedom *getDOFPermutation(const DimVec<int> &vertexPermutation);
protected:
- const BasisFunction *basFcts_;
+ const BasisFunction *basFcts;
- std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> > dofPermutation_;
+ std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> > dofPermutation;
- std::map<DimVec<double>, int, DimVecLess<double> > indexOfCoords_;
+ std::map<DimVec<double>, int, DimVecLess<double> > indexOfCoords;
- static std::vector<PeriodicDOFMapping*> mappings_;
+ static std::vector<PeriodicDOFMapping*> mappings;
};
@@ -80,7 +80,7 @@ namespace AMDiS {
{
public:
/// Constructor.
- PeriodicBC(BoundaryType type, FiniteElemSpace *rowFESpace_);
+ PeriodicBC(BoundaryType type, FiniteElemSpace *rowFESpace);
~PeriodicBC();
@@ -96,14 +96,20 @@ namespace AMDiS {
void exitVector(DOFVectorBase<double>* vector);
+ /// We are defining periodic boundary conditions, so return always true here.
+ bool isPeriodic()
+ {
+ return true;
+ }
+
protected:
VertexVector *associated;
- PeriodicDOFMapping *periodicDOFMapping_;
+ PeriodicDOFMapping *periodicDOFMapping;
- DegreeOfFreedom *neighIndices_;
+ DegreeOfFreedom *neighIndices;
- DOFMatrix *masterMatrix_;
+ DOFMatrix *masterMatrix;
};
}
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 273629f5..7ea4f6ce 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -954,7 +954,6 @@ namespace AMDiS {
FUNCNAME("ProblemVec::addPeriodicBC()");
FiniteElemSpace *feSpace = componentSpaces[row];
-
PeriodicBC *periodic = new PeriodicBC(type, feSpace);
if (systemMatrix && (*systemMatrix)[row][col])
@@ -1027,7 +1026,6 @@ namespace AMDiS {
}
#endif
-
// == Traverse mesh (either sequentially or in parallel) and assemble. ==
@@ -1061,7 +1059,6 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
-
// == Finally, if we have assembled in parallel, we have to add the thread ==
// == private matrix and vector to the global one. ==
@@ -1110,7 +1107,6 @@ namespace AMDiS {
delete [] bound;
} // pragma omp parallel
-
}
void ProblemVec::assembleOnDifMeshes(FiniteElemSpace *rowFeSpace,
diff --git a/AMDiS/src/RCNeighbourList.cc b/AMDiS/src/RCNeighbourList.cc
index 9ec42078..862a6df1 100644
--- a/AMDiS/src/RCNeighbourList.cc
+++ b/AMDiS/src/RCNeighbourList.cc
@@ -14,17 +14,17 @@
namespace AMDiS {
- RCNeighbourList::RCNeighbourList(int maxEdgeNeigh) {
+ RCNeighbourList::RCNeighbourList(int maxEdgeNeigh)
+ {
rclist.resize(maxEdgeNeigh);
- for (int i = 0; i < maxEdgeNeigh; i++) {
+ for (int i = 0; i < maxEdgeNeigh; i++)
rclist[i] = new RCListElement;
- }
}
- RCNeighbourList::~RCNeighbourList() {
- for (int i = 0; i < static_cast<int>(rclist.size()); i++) {
+ RCNeighbourList::~RCNeighbourList()
+ {
+ for (int i = 0; i < static_cast<int>(rclist.size()); i++)
delete rclist[i];
- }
}
/****************************************************************************/
@@ -39,7 +39,7 @@ namespace AMDiS {
for (int i = 0; i < n_neigh; i++) {
Element *lel = rclist[i]->el;
- if (lel->getMark() >= 0 || lel->isLeaf()) {
+ if (lel->getMark() >= 0 || lel->isLeaf()) {
/****************************************************************************/
/* element must not be coarsend or element is a leaf element, reset the */
/* the coarsening flag on all those elements that have to be coarsend with */
@@ -51,7 +51,7 @@ namespace AMDiS {
rclist[j]->el->setMark(0);
return false;
- } else if (lel->getFirstChild()->getMark() >= 0 ||
+ } else if (lel->getFirstChild()->getMark() >= 0 ||
lel->getSecondChild()->getMark() >= 0) {
/****************************************************************************/
@@ -64,8 +64,8 @@ namespace AMDiS {
rclist[j]->el->setMark(0);
return false;
- } else if (!lel->getFirstChild()->isLeaf() ||
- !lel->getSecondChild()->isLeaf()) {
+ } else if (!lel->getFirstChild()->isLeaf() ||
+ !lel->getSecondChild()->isLeaf()) {
/****************************************************************************/
/* one of the element's children is not a leaf element; */
/* element may be coarsend after coarsening one of the children; try again */
@@ -78,18 +78,15 @@ namespace AMDiS {
/* either one element is a macro element or we can coarsen the patch */
/****************************************************************************/
if (rclist[i]->flag == 0) {
-
- std::deque<MacroElement*>::const_iterator mel;
-
Mesh* coarse_mesh = coarseningManager->getMesh();
+ std::deque<MacroElement*>::const_iterator mel;
// set in Mesh::coarsen()
for (mel = coarse_mesh->firstMacroElement();
- mel!=coarse_mesh->endOfMacroElements();
- mel++)
+ mel!=coarse_mesh->endOfMacroElements(); ++mel)
if ((*mel)->getElement() == lel)
break;
- TEST_EXIT(mel!=coarse_mesh->endOfMacroElements())
+ TEST_EXIT(mel != coarse_mesh->endOfMacroElements())
("incompatible coarsening patch found\n");
}
}
@@ -103,19 +100,19 @@ namespace AMDiS {
{
FUNCNAME("RCNeighbourList::getNeighOnPatch()");
- int j, k, dir;
- Element *el, *neigh;
-
for (int i = 0; i < n_neigh; i++) {
- el = rclist[i]->el;
+ Element *el = rclist[i]->el;
+ Element *neigh = NULL;
rclist[i]->no = i;
- for (dir = 0; dir < 2; dir++) {
- for (j = 0; j < n_neigh; j++) {
+ for (int dir = 0; dir < 2; dir++) {
+ int j = 0;
+ for (; j < n_neigh; j++) {
if ((neigh = rclist[j]->el) == el)
continue;
- for (k = 0; k < 2; k++) {
+ int k = 0;
+ for (; k < 2; k++) {
if (neigh->getDOF(2+k) == el->getDOF(3-dir)) {
rclist[i]->neigh[dir] = rclist[j];
rclist[i]->oppVertex[dir] = 3-k;
@@ -138,197 +135,166 @@ namespace AMDiS {
void RCNeighbourList::addDOFParent(int elIndex, DegreeOfFreedom* dof) // 3d
{
Element *el = rclist[elIndex]->el;
- RCListElement *neighbour;
- int node;
-
+ RCListElement *neighbour = NULL;
Mesh *coarse_mesh = coarseningManager->getMesh();
RCListElement *coarse_list = rclist[elIndex];
- if (coarse_mesh->getNumberOfDOFs(EDGE))
- {
- node = coarse_mesh->getNode(EDGE);
-
- /****************************************************************************/
- /* set the dof in the coarsening edge */
- /****************************************************************************/
+ if (coarse_mesh->getNumberOfDOFs(EDGE)) {
+ int node = coarse_mesh->getNode(EDGE);
- el->setDOF(node, dof);
-
- /****************************************************************************/
- /* and now those handed on by the children */
- /****************************************************************************/
-
- el->setDOF(node+1, const_cast<int*>( el->getFirstChild()->getDOF(node)));
- el->setDOF(node+2, const_cast<int*>( el->getFirstChild()->getDOF(node+1)));
- el->setDOF(node+5, const_cast<int*>( el->getFirstChild()->getDOF(node+3)));
-
- if (coarse_list->elType)
- {
- el->setDOF(node+3, const_cast<int*>( el->getSecondChild()->getDOF(node)));
- el->setDOF(node+4, const_cast<int*>( el->getSecondChild()->getDOF(node+1)));
- }
- else
- {
- el->setDOF(node+3, const_cast<int*>( el->getSecondChild()->getDOF(node+1)));
- el->setDOF(node+4, const_cast<int*>( el->getSecondChild()->getDOF(node)));
- }
+ /****************************************************************************/
+ /* set the dof in the coarsening edge */
+ /****************************************************************************/
+
+ el->setDOF(node, dof);
+
+ /****************************************************************************/
+ /* and now those handed on by the children */
+ /****************************************************************************/
+
+ el->setDOF(node + 1, const_cast<int*>( el->getFirstChild()->getDOF(node)));
+ el->setDOF(node + 2, const_cast<int*>( el->getFirstChild()->getDOF(node + 1)));
+ el->setDOF(node + 5, const_cast<int*>( el->getFirstChild()->getDOF(node + 3)));
+
+ if (coarse_list->elType) {
+ el->setDOF(node + 3, const_cast<int*>( el->getSecondChild()->getDOF(node)));
+ el->setDOF(node + 4, const_cast<int*>( el->getSecondChild()->getDOF(node + 1)));
+ } else {
+ el->setDOF(node + 3, const_cast<int*>( el->getSecondChild()->getDOF(node + 1)));
+ el->setDOF(node + 4, const_cast<int*>( el->getSecondChild()->getDOF(node)));
}
+ }
- if (coarse_mesh->getNumberOfDOFs(FACE))
- {
- node = coarse_mesh->getNode(FACE);
- /****************************************************************************/
- /* dof's at the faces within the patch: add new dof's if it is a boundary */
- /* face or neighbour is an element behind this element in the corsen list */
- /****************************************************************************/
-
- neighbour = coarse_list->neigh[0];
- if ((!neighbour) || (neighbour > coarse_list))
- {
- if (!el->getDOF(node+2))
- {
- el->setDOF(node+2, const_cast<int*>( coarse_mesh->getDOF(FACE))); /* face 2 */
- if (neighbour)
- neighbour->el->setDOF(node+coarse_list->oppVertex[0],
- const_cast<int*>( el->getDOF(node+2)));
- }
- }
- neighbour = coarse_list->neigh[1];
- if ((!neighbour) || (neighbour > coarse_list))
- {
- if (!el->getDOF(node+3))
- {
- el->setDOF(node+3, const_cast<int*>( coarse_mesh->getDOF(FACE))); /* face 3 */
- if (neighbour)
- neighbour->el->setDOF(node+coarse_list->oppVertex[1],
- const_cast<int*>( el->getDOF(node+3)));
- }
- }
- /****************************************************************************/
- /* and now those handed on by the children */
- /****************************************************************************/
-
- el->setDOF(node, const_cast<int*>( el->getSecondChild()->getDOF(node+3)));
- el->setDOF(node+1, const_cast<int*>( el->getFirstChild()->getDOF(node+3)));
+ if (coarse_mesh->getNumberOfDOFs(FACE)) {
+ int node = coarse_mesh->getNode(FACE);
+ /****************************************************************************/
+ /* dof's at the faces within the patch: add new dof's if it is a boundary */
+ /* face or neighbour is an element behind this element in the corsen list */
+ /****************************************************************************/
+
+ neighbour = coarse_list->neigh[0];
+ if (!neighbour || neighbour > coarse_list) {
+ if (!el->getDOF(node + 2)) {
+ // face 2
+ el->setDOF(node + 2, const_cast<int*>(coarse_mesh->getDOF(FACE)));
+ if (neighbour)
+ neighbour->el->setDOF(node + coarse_list->oppVertex[0],
+ const_cast<int*>(el->getDOF(node + 2)));
+ }
}
-
- if (coarse_mesh->getNumberOfDOFs(CENTER))
- {
- node = coarse_mesh->getNode(CENTER);
- if (!el->getDOF(node))
- el->setDOF(node, const_cast<int*>( coarse_mesh->getDOF(CENTER)));
+
+ neighbour = coarse_list->neigh[1];
+ if (!neighbour || neighbour > coarse_list) {
+ if (!el->getDOF(node + 3)) {
+ // face 3
+ el->setDOF(node + 3, const_cast<int*>(coarse_mesh->getDOF(FACE)));
+ if (neighbour)
+ neighbour->el->setDOF(node + coarse_list->oppVertex[1],
+ const_cast<int*>(el->getDOF(node + 3)));
+ }
}
- return;
+ /****************************************************************************/
+ /* and now those handed on by the children */
+ /****************************************************************************/
+
+ el->setDOF(node, const_cast<int*>(el->getSecondChild()->getDOF(node + 3)));
+ el->setDOF(node + 1, const_cast<int*>(el->getFirstChild()->getDOF(node + 3)));
+ }
+
+ if (coarse_mesh->getNumberOfDOFs(CENTER)) {
+ int node = coarse_mesh->getNode(CENTER);
+ if (!el->getDOF(node))
+ el->setDOF(node, const_cast<int*>(coarse_mesh->getDOF(CENTER)));
+ }
}
void RCNeighbourList::addDOFParents(int n_neigh) // 2d
{
- int i, node;
-
Mesh *coarse_mesh = coarseningManager->getMesh();
- if (coarse_mesh->getNumberOfDOFs(EDGE))
- {
- node = coarse_mesh->getNode(EDGE);
-
- /****************************************************************************/
- /* get dofs on the boundary of the coarsening patch from the children */
- /****************************************************************************/
- for (i = 0; i < n_neigh; i++)
- {
- rclist[i]->el->
- setDOF(node, const_cast<int*>( rclist[i]->el->getSecondChild()->getDOF(node+2)));
- rclist[i]->el->
- setDOF(node+1, const_cast<int*>( rclist[i]->el->getFirstChild()->getDOF(node+2)));
- }
- }
-
- if (coarse_mesh->getNumberOfDOFs(CENTER))
- {
- int node = coarse_mesh->getNode(CENTER);
-
- /****************************************************************************/
- /* get new dof on parents at the barycenter */
- /****************************************************************************/
- for (i = 0; i < n_neigh; i++)
- {
- if (!rclist[i]->el->getDOF(node))
- rclist[i]->el->setDOF(node, const_cast<int*>( coarse_mesh->getDOF(CENTER)));
- }
+ if (coarse_mesh->getNumberOfDOFs(EDGE)) {
+ int node = coarse_mesh->getNode(EDGE);
+
+ /****************************************************************************/
+ /* get dofs on the boundary of the coarsening patch from the children */
+ /****************************************************************************/
+ for (int i = 0; i < n_neigh; i++) {
+ rclist[i]->el->setDOF(node, const_cast<int*>(rclist[i]->el->getSecondChild()->getDOF(node + 2)));
+ rclist[i]->el->setDOF(node + 1, const_cast<int*>(rclist[i]->el->getFirstChild()->getDOF(node + 2)));
}
+ }
- return;
+ if (coarse_mesh->getNumberOfDOFs(CENTER)) {
+ int node = coarse_mesh->getNode(CENTER);
+
+ /****************************************************************************/
+ /* get new dof on parents at the barycenter */
+ /****************************************************************************/
+ for (int i = 0; i < n_neigh; i++)
+ if (!rclist[i]->el->getDOF(node))
+ rclist[i]->el->setDOF(node, const_cast<int*>(coarse_mesh->getDOF(CENTER)));
+ }
}
void RCNeighbourList::removeDOFParents(int n_neigh)
{
- int i, j, node;
Mesh *mesh = rclist[0]->el->getMesh();
int edges = mesh->getGeo(EDGE);
- if (mesh->getNumberOfDOFs(EDGE))
- {
- node = mesh->getNode(EDGE);
+ if (mesh->getNumberOfDOFs(EDGE)) {
+ int node = mesh->getNode(EDGE);
- for (i = 0; i < n_neigh; i++)
- {
- for (j = 0; j < edges; j++)
- rclist[i]->el->setDOF(node+j, NULL);
- }
- }
+ for (int i = 0; i < n_neigh; i++)
+ for (int j = 0; j < edges; j++)
+ rclist[i]->el->setDOF(node + j, NULL);
+ }
- if (mesh->getNumberOfDOFs(CENTER))
- {
- node = mesh->getNode(CENTER);
- for (i = 0; i < n_neigh; i++)
- {
- mesh->freeDOF(const_cast<int*>( rclist[i]->el->getDOF(node)), CENTER);
- rclist[i]->el->setDOF(node, NULL);
- }
+ if (mesh->getNumberOfDOFs(CENTER)) {
+ int node = mesh->getNode(CENTER);
+ for (int i = 0; i < n_neigh; i++) {
+ mesh->freeDOF(const_cast<int*>(rclist[i]->el->getDOF(node)), CENTER);
+ rclist[i]->el->setDOF(node, NULL);
}
+ }
}
void RCNeighbourList::removeDOFParent(int index)
{
Tetrahedron *el = dynamic_cast<Tetrahedron*>(rclist[index]->el);
- RCListElement *neigh;
- int j, node;
Mesh *mesh = el->getMesh();
int edges = mesh->getGeo(EDGE);
int faces = mesh->getGeo(FACE);
- if (mesh->getNumberOfDOFs(EDGE))
- {
- node = mesh->getNode(EDGE);
- for (j = 0; j < edges; j++)
- el->setDOF(node+j, NULL);
- }
-
- if (mesh->getNumberOfDOFs(FACE))
- {
- node = mesh->getNode(FACE);
-
- neigh = rclist[index]->neigh[0];
- if ((!neigh) || (neigh > rclist[index]))
- mesh->freeDOF(const_cast<int*>( el->getDOF(node+2)), FACE); /* face 2 */
+ if (mesh->getNumberOfDOFs(EDGE)) {
+ int node = mesh->getNode(EDGE);
+ for (int j = 0; j < edges; j++)
+ el->setDOF(node + j, NULL);
+ }
- neigh = rclist[index]->neigh[1];
- if ((!neigh) || (neigh > rclist[index]))
- mesh->freeDOF(const_cast<int*>( el->getDOF(node+3)), FACE); /* face 3 */
+ if (mesh->getNumberOfDOFs(FACE)) {
+ int node = mesh->getNode(FACE);
+ RCListElement *neigh = rclist[index]->neigh[0];
- for (j = 0; j < faces; j++)
- el->setDOF(node+j, NULL);
- }
+ // face 2
+ if (!neigh || neigh > rclist[index])
+ mesh->freeDOF(const_cast<int*>(el->getDOF(node + 2)), FACE);
+
+ neigh = rclist[index]->neigh[1];
+ // face 3
+ if (!neigh || neigh > rclist[index])
+ mesh->freeDOF(const_cast<int*>(el->getDOF(node + 3)), FACE);
+
+ for (int j = 0; j < faces; j++)
+ el->setDOF(node + j, NULL);
+ }
- if (mesh->getNumberOfDOFs(CENTER))
- {
- node = mesh->getNode(CENTER);
- mesh->freeDOF(const_cast<int*>( el->getDOF(node)), CENTER);
- el->setDOF(node, NULL);
- }
- return;
+ if (mesh->getNumberOfDOFs(CENTER)) {
+ int node = mesh->getNode(CENTER);
+ mesh->freeDOF(const_cast<int*>(el->getDOF(node)), CENTER);
+ el->setDOF(node, NULL);
+ }
}
RCNeighbourList *RCNeighbourList::periodicSplit(DegreeOfFreedom *edge[2],
@@ -388,9 +354,8 @@ namespace AMDiS {
count++;
}
- for (int i = 0; i < *n_neigh_periodic; i++) {
+ for (int i = 0; i < *n_neigh_periodic; i++)
periodicList->rclist[i]->no = i;
- }
return periodicList;
}
diff --git a/AMDiS/src/RefinementManager.cc b/AMDiS/src/RefinementManager.cc
index 0dc3649f..02c072cd 100644
--- a/AMDiS/src/RefinementManager.cc
+++ b/AMDiS/src/RefinementManager.cc
@@ -67,9 +67,8 @@ namespace AMDiS {
}
}
- if (newCoords) {
- setNewCoords(); // call of sub-class method
- }
+ if (newCoords)
+ setNewCoords(); // call of sub-class method
delete stack;
diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc
index 668151c0..3b9db073 100644
--- a/AMDiS/src/RefinementManager2d.cc
+++ b/AMDiS/src/RefinementManager2d.cc
@@ -19,11 +19,17 @@ namespace AMDiS {
{
FUNCNAME("RefinementManager::refineFunction()");
+// if (MPI::COMM_WORLD.Get_rank() == 0)
+// std::cout << "IN EL = " << el_info->getElement()->getIndex() << std::endl;
+
bool bound = false;
DegreeOfFreedom *edge[2];
RCNeighbourList* refineList = new RCNeighbourList(2);
if (el_info->getElement()->getMark() <= 0) {
+// if (MPI::COMM_WORLD.Get_rank() == 0)
+// std::cout << "RETURN!" << std::endl;
+
delete refineList;
// Element may not be refined.
@@ -84,16 +90,22 @@ namespace AMDiS {
firstNewDOF = newDOF;
if (lastNewDOF != -1) {
+ int i = 0;
for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
+// if (MPI::COMM_WORLD.Get_rank() == 0)
+// std::cout << "i = " << i++ << std::endl;
+
if (it->second) {
+// if (MPI::COMM_WORLD.Get_rank() == 0)
+// std::cout << "*: " << newDOF << " " << lastNewDOF << std::endl;
+
if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
(*(it->second))[edge[1][0]] == last_edge[1][0]) ||
((*(it->second))[edge[0][0]] == last_edge[1][0] &&
- (*(it->second))[edge[1][0]] == last_edge[0][0]))
- {
- (*(it->second))[lastNewDOF] = newDOF;
- (*(it->second))[newDOF] = lastNewDOF;
- }
+ (*(it->second))[edge[1][0]] == last_edge[0][0])) {
+ (*(it->second))[lastNewDOF] = newDOF;
+ (*(it->second))[newDOF] = lastNewDOF;
+ }
}
}
}
@@ -112,11 +124,13 @@ namespace AMDiS {
if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
(*(it->second))[first_edge[1][0]] == last_edge[1][0]) ||
((*(it->second))[first_edge[0][0]] == last_edge[1][0] &&
- (*(it->second))[first_edge[1][0]] == last_edge[0][0]))
- {
- (*(it->second))[lastNewDOF] = firstNewDOF;
- (*(it->second))[firstNewDOF] = lastNewDOF;
- }
+ (*(it->second))[first_edge[1][0]] == last_edge[0][0])) {
+// if (MPI::COMM_WORLD.Get_rank() == 0)
+// std::cout << "**: " << newDOF << " " << lastNewDOF << std::endl;
+
+ (*(it->second))[lastNewDOF] = firstNewDOF;
+ (*(it->second))[firstNewDOF] = lastNewDOF;
+ }
}
}
}
diff --git a/AMDiS/src/VertexVector.cc b/AMDiS/src/VertexVector.cc
index 919458be..4f17a8ec 100644
--- a/AMDiS/src/VertexVector.cc
+++ b/AMDiS/src/VertexVector.cc
@@ -4,12 +4,12 @@
namespace AMDiS {
- VertexVector::VertexVector(const DOFAdmin *admin_, std::string name_)
+ VertexVector::VertexVector(const DOFAdmin *a, std::string n)
: DOFVectorDOF()
{
- name = name_;
+ name = n;
feSpace = NULL;
- admin = admin_;
+ admin = a;
const_cast<DOFAdmin*>(admin)->addDOFIndexed(this);
const_cast<DOFAdmin*>(admin)->addDOFContainer(this);
}
@@ -23,9 +23,36 @@ namespace AMDiS {
void VertexVector::set(DegreeOfFreedom val)
{
DOFIteratorBase it(const_cast<DOFAdmin*>(admin), USED_DOFS);
- for(it.reset(); !it.end(); ++it) {
- if(!it.isDOFFree()) operator[](it.getDOFIndex()) = val;
- }
+ for (it.reset(); !it.end(); ++it)
+ if (!it.isDOFFree())
+ operator[](it.getDOFIndex()) = val;
}
+ void VertexVector::print()
+ {
+ DOFIteratorBase it(const_cast<DOFAdmin*>(admin), USED_DOFS);
+ for (it.reset(); !it.end(); ++it)
+ if (!it.isDOFFree())
+ std::cout << "DOF " << it.getDOFIndex() << " -> "
+ << operator[](it.getDOFIndex()) << "\n";
+ }
+
+ void VertexVector::changeDofIndices(std::map<int, int>& dofIndexMap)
+ {
+ FUNCNAME("VertexVector::changedofIndices()");
+
+ std::vector<DegreeOfFreedom> tmp(vec.size());
+ for (int i = 0; i < static_cast<int>(tmp.size()); i++)
+ tmp[i] = i;
+
+ for (std::map<int, int>::iterator it = dofIndexMap.begin();
+ it != dofIndexMap.end(); ++it)
+ if (vec[it->first] == -1)
+ tmp[it->second] = -1;
+ else
+ tmp[it->second] = dofIndexMap[vec[it->first]];
+
+ vec = tmp;
+ }
+
}
diff --git a/AMDiS/src/VertexVector.h b/AMDiS/src/VertexVector.h
index 6601f692..f44a587e 100644
--- a/AMDiS/src/VertexVector.h
+++ b/AMDiS/src/VertexVector.h
@@ -17,7 +17,7 @@ namespace AMDiS {
{}
};
- VertexVector(const DOFAdmin *admin_, std::string name_);
+ VertexVector(const DOFAdmin *admin, std::string name);
~VertexVector();
@@ -40,10 +40,14 @@ namespace AMDiS {
{
DOFContainer::compressDOFContainer(size, newDOF);
int totalSize = getAdmin()->getSize();
- for(int i = size; i < totalSize; i++)
+ for (int i = size; i < totalSize; i++)
vec[i] = i;
}
+ void changeDofIndices(std::map<int, int>& dofIndexMap);
+
+ void print();
+
protected:
const DOFAdmin *admin;
};
diff --git a/AMDiS/src/VtkWriter.cc b/AMDiS/src/VtkWriter.cc
index 0230281f..9937dae8 100644
--- a/AMDiS/src/VtkWriter.cc
+++ b/AMDiS/src/VtkWriter.cc
@@ -16,6 +16,8 @@ namespace AMDiS {
int VtkWriter::writeFile(std::string name)
{
+ FUNCNAME("VtkWriter::writeFile()");
+
#ifdef HAVE_BOOST
boost::iostreams::filtering_ostream file;
--
GitLab