diff --git a/AMDiS/libtool b/AMDiS/libtool
index 3b6673371f6b9a9c8f45abefeac61e4c5b70b66e..5d13712732963994e0b5b46edc18f556108c2c8d 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 p2q024:
# 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,25 +82,25 @@ 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
-gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_dir=`gcc -print-file-name=. | /usr/bin/sed 's,/\.$,,'`
gcc_ver=`gcc -dumpversion`
# 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"
@@ -174,7 +174,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"
@@ -328,10 +328,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=`echo " /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/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+sys_lib_search_path_spec=`echo "/lib64 /usr/lib64 /usr/local/lib64" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# 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 "
+sys_lib_dlsearch_path_spec="/lib64 /usr/lib64 /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=""
@@ -6763,7 +6763,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 p2q024:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6785,12 +6785,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.
@@ -6801,25 +6801,25 @@ 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
-gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_dir=`gcc -print-file-name=. | /usr/bin/sed 's,/\.$,,'`
gcc_ver=`gcc -dumpversion`
# 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"
@@ -6893,7 +6893,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"
@@ -6948,11 +6948,11 @@ striplib="strip --strip-unneeded"
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=`echo "/usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crti.o /usr/lib/gcc/i386-redhat-linux/4.1.2/crtbeginS.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+predep_objects=`echo "/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" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=`echo "/usr/lib/gcc/i386-redhat-linux/4.1.2/crtendS.o /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crtn.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+postdep_objects=`echo "/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" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -6960,11 +6960,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=`echo "-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/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+compiler_lib_search_path=`echo "-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/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7044,10 +7044,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=`echo " /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/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+sys_lib_search_path_spec=`echo "/lib64 /usr/lib64 /usr/local/lib64" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# 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 "
+sys_lib_dlsearch_path_spec="/lib64 /usr/lib64 /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=""
@@ -7071,7 +7071,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host p2q024:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7093,12 +7093,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.
@@ -7109,7 +7109,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"
@@ -7118,16 +7118,16 @@ LTCFLAGS="-g -O2"
CC="g77"
# Is the compiler the GNU C compiler?
-with_gcc=yes
+with_gcc=
-gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_dir=`gcc -print-file-name=. | /usr/bin/sed 's,/\.$,,'`
gcc_ver=`gcc -dumpversion`
# 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"
@@ -7355,10 +7355,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=`echo " /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/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+sys_lib_search_path_spec=`echo "/lib64 /usr/lib64 /usr/local/lib64" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# 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 "
+sys_lib_dlsearch_path_spec="/lib64 /usr/lib64 /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/BasisFunction.h b/AMDiS/src/BasisFunction.h
index b48adf2ad49c5198821576faca1a44a10c3f5372..9b55d279ebc60f826368b8d8f4226ec5298c4122 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -133,7 +133,7 @@ namespace AMDiS {
* all routines using this function on the elements need the FILL_BOUND
* flag during mesh traversal;
*/
- virtual void getBound(const ElInfo*, BoundaryType *) const {};
+ virtual void getBound(const ElInfo*, BoundaryType *) const {}
/// Returns \ref degree of BasisFunction
inline const int getDegree() const
@@ -197,11 +197,11 @@ namespace AMDiS {
double *coeff) = 0;
/// WorldVector<double> valued interpol function.
- virtual const WorldVector<double>*
- interpol(const ElInfo *el_info, int no,
- const int *b_no,
- AbstractFunction<WorldVector<double>,WorldVector<double> > *f,
- WorldVector<double> *vec) = 0;
+ virtual const WorldVector<double>* interpol(const ElInfo *el_info, int no,
+ const int *b_no,
+ AbstractFunction<WorldVector<double>,
+ WorldVector<double> > *f,
+ WorldVector<double> *vec) = 0;
/// Returns the i-th local basis function
inline BasFctType *getPhi(int i) const
@@ -305,7 +305,8 @@ namespace AMDiS {
* will be overwritten after the next call.
*/
const WorldVector<double>& evalUh(const DimVec<double>& lambda,
- const WorldVector<double>* uh, WorldVector<double>* val) const;
+ const WorldVector<double>* uh,
+ WorldVector<double>* val) const;
/** \brief
* Evaluates the gradient at barycentric coordinates lambda. Lambda is the
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index bfc2430cc86162c1fd418d8bad22c7dbc1a52492..85e2e9dd9b3e7ad89877f6895e48eec2c0d0bb20 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -47,6 +47,8 @@ namespace AMDiS {
WorldVector<double>& w) const
{
+ testFlag(Mesh::FILL_COORDS);
+
double c = l[0];
for (int j = 0; j < dimOfWorld; j++)
diff --git a/AMDiS/src/ElementDofIterator.cc b/AMDiS/src/ElementDofIterator.cc
index b48ef8a6f0811ed78ab4ca276c52bbd387953273..338e37cd957822f34981164bbe788f9a69eb3f16 100644
--- a/AMDiS/src/ElementDofIterator.cc
+++ b/AMDiS/src/ElementDofIterator.cc
@@ -47,6 +47,10 @@ namespace AMDiS {
// We are finished with all element.
elementPos = 0;
+ // If we have iterated over all positions, we can finish the iteration.
+ if (pos >= dim)
+ return false;
+
// Increase position, i.e., go from vertices to edges to faces and search
// for the next position with dofs.
do {
@@ -55,16 +59,18 @@ namespace AMDiS {
posIndex = INDEX_OF_DIM(pos, dim);
// Get number of dofs in this position.
nDofs = admin->getNumberOfDOFs(posIndex);
- } while (nDofs == 0 && pos <= dim);
+ } while (nDofs == 0 && pos < dim);
- if (pos <= dim) {
+ if (nDofs > 0 && pos <= dim) {
// We have found on more position with dofs.
// Get number of elements in this position, i.e, the number of vertices,.
// edges and faces in the given dimension.
nElements = Global::getGeo(posIndex, dim);
+
// Calculate displacement. Is used if there is more than one dof admin on the mesh.
n0 = admin->getNumberOfPreDOFs(posIndex);
+
// Get first dof index position for the geo index position.
node0 = mesh->getNode(posIndex);
} else {
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index 41815863702e377cf39b26531e30024ff2ee9daa..c5e64c2ed068c665b6c229e2ecbc7bf6aa0036c3 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -292,10 +292,9 @@ namespace AMDiS {
macroInfo->fillBoundaryInfo(mesh);
if (mesh->getNumberOfDOFs(CENTER)) {
- for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
+ for (int i = 0; i < mesh->getNumberOfMacros(); i++)
const_cast<Element*>(mel[i]->getElement())->
setDOF(mesh->getNode(CENTER), mesh->getDOF(CENTER));
- }
}
/****************************************************************************/
@@ -1137,11 +1136,11 @@ namespace AMDiS {
mesh->incrementNumberOfFaces(1);
if (mesh->getNumberOfDOFs(FACE)) {
- TEST_EXIT(!(*(mel+i))->getElement()->getDOF(lnode+k))
+ TEST_EXIT(!(*(mel+i))->getElement()->getDOF(lnode + k))
("dof %d on element %d already set\n",
lnode+k, (*(mel+i))->getIndex());
- const_cast<Element*>((*(mel+i))->getElement())->setDOF(lnode+k,
+ const_cast<Element*>((*(mel+i))->getElement())->setDOF(lnode + k,
mesh->getDOF(FACE));
if (neigh) {
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 98b3535ef5473c0fbfe58bb2bc174166c1503058..b3093503906fbc9e62d9bf0599a67d9812579804 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -824,9 +824,8 @@ namespace AMDiS {
c_lambda2[3]);
if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
- for (int i = 0; i <= dim; i++) {
- c_lambda[i] = c_lambda2[i];
- }
+ for (int i = 0; i <= dim; i++)
+ c_lambda[i] = c_lambda2[i];
c_outside = c_outside2;
*c_el_info = *c_el_info2;
ichild = 1 - ichild;
@@ -862,6 +861,33 @@ namespace AMDiS {
return inside;
}
+ bool Mesh::getDofIndexCoords(const DegreeOfFreedom* dof,
+ const FiniteElemSpace* feSpace,
+ WorldVector<double>& coords)
+ {
+ DimVec<double>* baryCoords;
+ bool found = false;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(this, -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+ while (elInfo) {
+ for (int i = 0; i < nDOFEl; i++) {
+ if (elInfo->getElement()->getDOF(i) == dof) {
+ baryCoords = feSpace->getBasisFcts()->getCoords(i);
+ elInfo->coordToWorld(*baryCoords, coords);
+ found = true;
+ break;
+ }
+ }
+
+ if (found)
+ break;
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ return found;
+ }
void Mesh::setDiameter(const WorldVector<double>& w)
{
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index c6020f0c04b0c6920e612acbf50ecbf18a1a919b..4fc0f3e71311ccf67fd91618f59151ee6db924d7 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -484,6 +484,23 @@ namespace AMDiS {
const WorldVector<double> *xy0,
double *sp);
+ /** \brief
+ * Returns for a given dof its world coordinates in this mesh. Because we do
+ * not have any direct connection between dofs and coordinates, this function
+ * has to search for the element in this mesh, that contains the dof. Than the
+ * coordinates can be computed. Therefore, this function is very costly and
+ * should be used for debugging purpose only.
+ *
+ * @param[in] dof A pointer to the dof we have to search for.
+ * @param[in] feSpace The fe space to be used for the search.
+ * @param[out] coords World vector that stores the coordinates of the dof.
+ *
+ * The function returns true, if the dof was found, otherwise false.
+ */
+ bool getDofIndexCoords(const DegreeOfFreedom* dof,
+ const FiniteElemSpace* feSpace,
+ WorldVector<double>& coords);
+
/// Returns FILL_ANY_?D
inline static const Flag& getFillAnyFlag(int dim)
diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc
index ff54a07bfccba62732ea248c572e4d67d2e03581..e457e1007a13e3a2fb97e89ff4a4127198fa57f5 100644
--- a/AMDiS/src/ParallelDomainProblem.cc
+++ b/AMDiS/src/ParallelDomainProblem.cc
@@ -11,16 +11,25 @@
#include "DOFMatrix.h"
#include "DOFVector.h"
#include "VtkWriter.h"
+#include "ElementDofIterator.h"
#include "petscksp.h"
namespace AMDiS {
- ParallelDomainProblemBase::ParallelDomainProblemBase(const std::string& name,
- ProblemIterationInterface *iIF,
- ProblemTimeInterface *tIF,
- FiniteElemSpace *fe,
- RefinementManager *refineManager)
+ PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
+ {
+ if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
+ std::cout << " Iteration " << iter << ": " << rnorm << "\n";
+
+ return 0;
+ }
+
+ ParallelDomainBase::ParallelDomainBase(const std::string& name,
+ ProblemIterationInterface *iIF,
+ ProblemTimeInterface *tIF,
+ FiniteElemSpace *fe,
+ RefinementManager *refineManager)
: iterationIF(iIF),
timeIF(tIF),
feSpace(fe),
@@ -35,11 +44,16 @@ namespace AMDiS {
partitioner = new ParMetisPartitioner(mesh, &mpiComm);
}
- void ParallelDomainProblemBase::initParallelization(AdaptInfo *adaptInfo)
+ void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
{
if (mpiSize <= 1)
return;
+ // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
+ // only for macro meshes, so it will not work yet if the mesh is already refined
+ // in some way.
+ testForMacroMesh();
+
// create an initial partitioning of the mesh
partitioner->createPartitionData();
// set the element weights, which are 1 at the very first begin
@@ -86,6 +100,8 @@ namespace AMDiS {
admin.setFirstHole(mapLocalGlobalDOFs.size());
}
+ exit(0);
+
// === Global refinements. ===
int globalRefinement = 0;
@@ -98,9 +114,9 @@ namespace AMDiS {
}
-#if (DEBUG != 0)
- testInteriorBoundary();
-#endif
+// #if (DEBUG != 0)
+// testInteriorBoundary();
+// #endif
// === Create petsc matrix. ===
@@ -118,12 +134,36 @@ namespace AMDiS {
ierr = VecSetType(petscSolVec, VECMPI);
}
- void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
- {}
+
+ void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
+ {
+ }
- void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat,
- DOFVector<double> *vec)
+ void ParallelDomainBase::testForMacroMesh()
+ {
+ FUNCNAME("ParallelDomainBase::testForMacroMesh()");
+
+ int nMacroElements = 0;
+
+ TraverseStack stack;
+ 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");
+
+ nMacroElements++;
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ TEST_EXIT(nMacroElements >= mpiSize)
+ ("The mesh has less macro elements than number of mpi processes!\n");
+ }
+
+
+ void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat,
+ DOFVector<double> *vec)
{
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits= mtl::traits;
@@ -136,8 +176,10 @@ namespace AMDiS {
typedef traits::range_generator<major, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
- for (cursor_type cursor = begin<major>(mat->getBaseMatrix()), cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor)
- for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
+ for (cursor_type cursor = begin<major>(mat->getBaseMatrix()),
+ cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor)
+ for (icursor_type icursor = begin<nz>(cursor),
+ icend = end<nz>(cursor); icursor != icend; ++icursor)
if (value(*icursor) != 0.0) {
int r = mapLocalGlobalDOFs[row(*icursor)];
int c = mapLocalGlobalDOFs[col(*icursor)];
@@ -161,10 +203,8 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
+ void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
{
- clock_t t = clock();
-
KSP ksp;
PC pc;
@@ -174,7 +214,7 @@ namespace AMDiS {
PCSetType(pc, PCJACOBI);
KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSetType(ksp, KSPBCGS);
- KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, 0);
+ KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
KSPSolve(ksp, petscRhsVec, petscSolVec);
PetscScalar *vecPointer;
@@ -237,12 +277,10 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
delete [] sendBuffers[i];
-
- std::cout << "SOLUTION = " << TIME_USED(t,clock()) << std::endl;
}
- double ParallelDomainProblemBase::setElemWeights(AdaptInfo *adaptInfo)
+ double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo)
{
double localWeightSum = 0.0;
int elNum = -1;
@@ -277,7 +315,7 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::partitionMesh(AdaptInfo *adaptInfo)
+ void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
{
if (initialPartitionMesh) {
initialPartitionMesh = false;
@@ -292,15 +330,16 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
- std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
+ void ParallelDomainBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
+ std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
{
- FUNCNAME("ParallelDomainProblemBase::createInteriorBoundaryInfo()");
+ FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
// === 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);
while (elInfo) {
Element *element = elInfo->getElement();
@@ -313,11 +352,12 @@ namespace AMDiS {
PartitionElementData *neighbourPartitionData =
dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED));
+
if (neighbourPartitionData->getPartitionStatus() == OUT) {
// We have found an element that is at an interior boundary.
- // === Find out, if the boundary part of the element corresponds to the
- // rank or to the rank "on the other side" of the interoir boundary. ===
+ // === Find out, if the boundary part of the element corresponds to the ===
+ // === rank or to the rank "on the other side" of the interoir boundary. ===
const DegreeOfFreedom* boundDOF1 = NULL;
const DegreeOfFreedom* boundDOF2 = NULL;
@@ -442,7 +482,7 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::removeMacroElements()
+ void ParallelDomainBase::removeMacroElements()
{
std::vector<MacroElement*> macrosToRemove;
for (std::deque<MacroElement*>::iterator it = mesh->firstMacroElement();
@@ -459,13 +499,13 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::createLocalGlobalNumbering(
+ void ParallelDomainBase::createLocalGlobalNumbering(
std::vector<const DegreeOfFreedom*>& rankDOFs,
std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
int& nRankDOFs,
int& nOverallDOFs)
{
- FUNCNAME("ParallelDomainProblemBase::createLocalGlobalNumbering()");
+ FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
// === Get rank information about DOFs. ===
@@ -549,6 +589,7 @@ namespace AMDiS {
++sendIt, i++) {
int nSendDofs = sendIt->second.size() * 2;
sendBuffers[i] = new int[nSendDofs];
+
int c = 0;
for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator
dofIt = sendIt->second.begin();
@@ -641,10 +682,10 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::updateLocalGlobalNumbering(int& nRankDOFs,
- int& nOverallDOFs)
+ void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs,
+ int& nOverallDOFs)
{
- FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");
+ FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
std::set<const DegreeOfFreedom*> rankDOFs;
std::map<const DegreeOfFreedom*, int> newBoundaryDOFs;
@@ -664,7 +705,7 @@ namespace AMDiS {
// === Traverse on interior boundaries and move all not ranked owned DOFs from ===
- // === rankDOFs to boundaryDOFs ===
+ // === rankDOFs to boundaryDOFs. ===
std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;
@@ -871,10 +912,10 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::addAllDOFs(Element *el, int ithEdge,
- std::vector<const DegreeOfFreedom*>& dofs)
+ void ParallelDomainBase::addAllDOFs(Element *el, int ithEdge,
+ std::vector<const DegreeOfFreedom*>& dofs)
{
- FUNCNAME("ParallelDomainProblemBase::addAllDOFs()");
+ FUNCNAME("ParallelDomainBase::addAllDOFs()");
switch (ithEdge) {
case 0:
@@ -904,25 +945,30 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::createDOFMemberInfo(
+ void ParallelDomainBase::createDOFMemberInfo(
std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
std::vector<const DegreeOfFreedom*>& rankDOFs,
std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
{
// === Determine to each dof the set of partitions the dof belongs to. ===
+ ElementDofIterator elDofIt(feSpace->getAdmin());
+
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
Element *element = elInfo->getElement();
- // Determine to each dof the partition(s) it corresponds to.
- for (int i = 0; i < 3; i++)
- partitionDOFs[element->getDOF(i)].insert(partitionVec[element->getIndex()]);
-
+ elDofIt.reset(element);
+ do {
+ // Determine to each dof the partition(s) it corresponds to.
+ partitionDOFs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
+ } while(elDofIt.next());
+
elInfo = stack.traverseNext(elInfo);
}
+
// === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
// iterate over all DOFs
@@ -968,37 +1014,73 @@ namespace AMDiS {
}
- void ParallelDomainProblemBase::testInteriorBoundary()
+ void ParallelDomainBase::testInteriorBoundary()
{
+ FUNCNAME("ParallelDomainBase::testInteriorBoundary()");
+
// Maps to each neighbour rank an array of WorldVectors. This array contain the
// coordinates of all dofs this rank shares on the interior boundary with the
- // neighbour rank.
+ // neighbour rank. A rank sends the coordinates to another rank, if it owns the
+ // boundarys dof.
std::map<int, std::vector<WorldVector<double> > > sendCoords;
- std::map<int, int> nRecvCoords;
-
- for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDOFs.begin();
- it != boundaryDOFs.end();
- ++it) {
-
+ // A rank receives all boundary dofs that are at its interior boundaries but are
+ // not owned by the rank. This map stores for each rank the the coordinates of dofs
+ // this rank expectes to receive from.
+ std::map<int, std::vector<WorldVector<double> > > recvCoords;
+
+ WorldVector<double> coords;
+
+ for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator it =
+ sendDofs.begin();
+ it != sendDofs.end(); ++it) {
+ for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
+ dofIt != it->second.end(); ++dofIt) {
+ bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
+ TEST_EXIT(b)("Cannot find DOF in mesh!\n");
+ sendCoords[it->first].push_back(coords);
+ }
+ }
+
+ for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator it =
+ recvDofs.begin();
+ it != recvDofs.end(); ++it) {
+ for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
+ dofIt != it->second.end(); ++dofIt) {
+ bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
+ TEST_EXIT(b)("Cannot find DOF in mesh!\n");
+ recvCoords[it->first].push_back(coords);
+ }
+ }
+
+ for (std::map<int, std::vector<WorldVector<double> > >::iterator it
+ = sendCoords.begin();
+ it != sendCoords.end(); ++it) {
+ }
+
+ for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin();
+ it != recvCoords.end(); ++it) {
}
}
- ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name,
- ProblemScal *problem,
- ProblemInstatScal *problemInstat)
- : ParallelDomainProblemBase(name,
- problem,
- problemInstat,
- problem->getFESpace(),
- problem->getRefinementManager()),
+ ParallelDomainScal::ParallelDomainScal(const std::string& name,
+ ProblemScal *problem,
+ ProblemInstatScal *problemInstat)
+ : ParallelDomainBase(name,
+ problem,
+ problemInstat,
+ problem->getFESpace(),
+ problem->getRefinementManager()),
probScal(problem)
{
+ info = problem->getInfo();
}
- void ParallelDomainProblemScal::initParallelization(AdaptInfo *adaptInfo)
+ void ParallelDomainScal::initParallelization(AdaptInfo *adaptInfo)
{
- ParallelDomainProblemBase::initParallelization(adaptInfo);
+ FUNCNAME("ParallelDomainScal::initParallelization()");
+
+ ParallelDomainBase::initParallelization(adaptInfo);
DOFMatrix* m = probScal->getSystemMatrix();
@@ -1007,21 +1089,40 @@ namespace AMDiS {
m->setIsRankDOF(isRankDOF);
}
- Flag ParallelDomainProblemScal::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
+ void ParallelDomainScal::solve()
+ {
+ FUNCNAME("ParallelDomainScal::solve()");
+
+#ifdef _OPENMP
+ double wtime = omp_get_wtime();
+#endif
+ clock_t first = clock();
+
+ fillPetscMatrix(probScal->getSystemMatrix(), probScal->getRHS());
+ solvePetscMatrix(probScal->getSolution());
+
+#ifdef _OPENMP
+ INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
+ TIME_USED(first, clock()),
+ omp_get_wtime() - wtime);
+#else
+ INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
+ TIME_USED(first, clock()));
+#endif
+ }
+
+ Flag ParallelDomainScal::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
{
- FUNCNAME("ParallelDomainProblemScal::oneIteration()");
+ FUNCNAME("ParallelDomainScal::oneIteration()");
Flag flag = dynamic_cast<StandardProblemIteration*>(iterationIF)->
buildAndAdapt(adaptInfo, toDo);
- if (toDo.isSet(SOLVE)) {
- fillPetscMatrix(probScal->getSystemMatrix(), probScal->getRHS());
- solvePetscMatrix(probScal->getSolution());
- }
+ if (toDo.isSet(SOLVE))
+ solve();
- if (toDo.isSet(SOLVE_RHS)) {
+ if (toDo.isSet(SOLVE_RHS))
ERROR_EXIT("Not yet implemented!\n");
- }
if (toDo.isSet(ESTIMATE))
iterationIF->getProblem()->estimate(adaptInfo);
diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h
index 333316a188aef541210a83c706b0d2db89ee06e8..56d3d030e8d1428a40e21165a953685498d1ba89 100644
--- a/AMDiS/src/ParallelDomainProblem.h
+++ b/AMDiS/src/ParallelDomainProblem.h
@@ -17,10 +17,10 @@
// == ==
// ============================================================================
-/** \file ParallelDomainProblem.h */
+/** \file ParallelDomain.h */
-#ifndef AMDIS_PARALLELDOMAINPROBLEM_H
-#define AMDIS_PARALLELDOMAINPROBLEM_H
+#ifndef AMDIS_PARALLELDOMAIN_H
+#define AMDIS_PARALLELDOMAIN_H
#include <map>
#include <set>
@@ -42,21 +42,29 @@ namespace AMDiS {
class ParMetisPartitioner;
- class ParallelDomainProblemBase : public ProblemIterationInterface,
+ class ParallelDomainBase : public ProblemIterationInterface,
public ProblemTimeInterface
{
public:
- ParallelDomainProblemBase(const std::string& name,
+ ParallelDomainBase(const std::string& name,
ProblemIterationInterface *iterationIF,
ProblemTimeInterface *timeIF,
FiniteElemSpace *feSpace,
RefinementManager *refineManager);
- virtual ~ParallelDomainProblemBase() {}
+ virtual ~ParallelDomainBase() {}
virtual void initParallelization(AdaptInfo *adaptInfo);
- void exitParallelization(AdaptInfo *adaptInfo);
+ void exitParallelization(AdaptInfo *adaptInfo);
+
+ /** \brief
+ * Test, if the mesh consists of macro elements only. The mesh partitioning of
+ * the parallelization works for macro meshes only and would fail, if the mesh
+ * is already refined in some way. Therefore, this function will exit the program
+ * if it finds a non macro element in the mesh.
+ */
+ void testForMacroMesh();
/// Set for each element on the partitioning level the number of leaf elements.
double setElemWeights(AdaptInfo *adaptInfo);
@@ -206,6 +214,9 @@ namespace AMDiS {
/// Mesh of the problem.
Mesh *mesh;
+ /// Info level.
+ int info;
+
/// Refinement manager for the mesh.
RefinementManager *refinementManager;
@@ -287,10 +298,10 @@ namespace AMDiS {
std::map<DegreeOfFreedom, bool> isRankDOF;
};
- class ParallelDomainProblemScal : public ParallelDomainProblemBase
+ class ParallelDomainScal : public ParallelDomainBase
{
public:
- ParallelDomainProblemScal(const std::string& name,
+ ParallelDomainScal(const std::string& name,
ProblemScal *problem,
ProblemInstatScal *problemInstat);
@@ -299,9 +310,13 @@ namespace AMDiS {
virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION);
protected:
- ProblemScal *probScal;
+ /// Starts the solution of the linear system using Petsc.
+ void solve();
+ protected:
+ /// Pointer to the stationary problem.
+ ProblemScal *probScal;
};
}
-#endif // AMDIS_PARALLELDOMAINPROBLEM_H
+#endif // AMDIS_PARALLELDOMAIN_H
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index c71a268e4e4e7b4fb9b4d701095f8de4ee7c7ab4..0a4b6df5eefaa007d9e0e3cec6b4bdcdf8139cf3 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -214,6 +214,7 @@ namespace AMDiS {
void ProblemScal::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
{
FUNCNAME("Problem::solve()");
+
if (!solver) {
WARNING("no solver\n");
return;
@@ -222,7 +223,6 @@ namespace AMDiS {
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
-
clock_t first = clock();
SolverMatrix<DOFMatrix> solverMatrix;
diff --git a/AMDiS/src/ProblemScal.h b/AMDiS/src/ProblemScal.h
index b3ecf94066a4f6efba77501c222ecdd1fbe9e2be..429c4d2610debea303331c1a77de67310de894e2 100644
--- a/AMDiS/src/ProblemScal.h
+++ b/AMDiS/src/ProblemScal.h
@@ -310,17 +310,24 @@ namespace AMDiS {
void writeResidualMesh(AdaptInfo *adaptInfo, const std::string name);
- /// serialization
+ /// Serialization of this class to a stream.
virtual void serialize(std::ostream &out);
- /// deserialization
+ /// Deserialization of this class from a stream.
virtual void deserialize(std::istream &in);
+ /// Returns \ref fileWriters.
std::vector<FileWriterInterface*> getFileWriterList()
{
return fileWriters;
}
+ /// Returns \ref info.
+ int getInfo()
+ {
+ return info;
+ }
+
protected:
/// Name of this problem.
std::string name;
@@ -371,7 +378,6 @@ namespace AMDiS {
/// Info level.
int info;
-
};
}