From 7908ff204c0e2dd4651893c753a80b12d300ec4c Mon Sep 17 00:00:00 2001 From: Thomas Witkowski Date: Tue, 19 Jan 2010 11:09:48 +0000 Subject: [PATCH] 3d adaptivity for parallelization. --- AMDiS/Makefile.in | 1 - AMDiS/aclocal.m4 | 34 +- AMDiS/bin/Makefile.in | 1 - AMDiS/configure | 208 +++-------- AMDiS/other/include/Makefile_AMDiS.mk | 5 +- AMDiS/src/AMDiS_fwd.h | 4 + AMDiS/src/BasisFunction.h | 12 +- AMDiS/src/BoundaryManager.cc | 9 +- AMDiS/src/BoundaryManager.h | 2 +- AMDiS/src/DOFAdmin.cc | 2 +- AMDiS/src/DOFAdmin.h | 2 +- AMDiS/src/DOFVector.hh | 14 +- AMDiS/src/Element.cc | 59 ++- AMDiS/src/Element.h | 12 +- AMDiS/src/Global.cc | 11 +- AMDiS/src/Global.h | 8 +- AMDiS/src/InteriorBoundary.cc | 43 +++ AMDiS/src/InteriorBoundary.h | 7 + AMDiS/src/Lagrange.cc | 14 +- AMDiS/src/Lagrange.h | 5 - AMDiS/src/Line.h | 4 +- AMDiS/src/MacroElement.cc | 2 +- AMDiS/src/Mesh.cc | 138 ++++--- AMDiS/src/Mesh.h | 23 +- AMDiS/src/MeshStructure.cc | 68 +++- AMDiS/src/MeshStructure.h | 26 +- AMDiS/src/ParMetisPartitioner.cc | 69 ++-- AMDiS/src/ParMetisPartitioner.h | 6 +- AMDiS/src/ParallelDomainBase.cc | 512 ++++++++++++++++---------- AMDiS/src/ParallelDomainBase.h | 38 +- AMDiS/src/ParallelDomainVec.cc | 4 - AMDiS/src/PngWriter.cc | 8 +- AMDiS/src/ProblemVec.cc | 12 +- AMDiS/src/RefinementManager.cc | 41 ++- AMDiS/src/RefinementManager.h | 2 + AMDiS/src/RefinementManager3d.cc | 135 +++---- AMDiS/src/StdMpi.h | 54 ++- AMDiS/src/Tetrahedron.cc | 71 ++-- AMDiS/src/Tetrahedron.h | 4 +- AMDiS/src/Triangle.cc | 16 +- AMDiS/src/Triangle.h | 4 +- AMDiS/src/VtkWriter.h | 6 + 42 files changed, 993 insertions(+), 703 deletions(-) diff --git a/AMDiS/Makefile.in b/AMDiS/Makefile.in index e7566807..e9c5bf1a 100644 --- a/AMDiS/Makefile.in +++ b/AMDiS/Makefile.in @@ -140,7 +140,6 @@ PACKAGE_VERSION = @PACKAGE_VERSION@ PATH_SEPARATOR = @PATH_SEPARATOR@ PETSC_DIR = @PETSC_DIR@ RANLIB = @RANLIB@ -SED = @SED@ SET_MAKE = @SET_MAKE@ SHELL = @SHELL@ STRIP = @STRIP@ diff --git a/AMDiS/aclocal.m4 b/AMDiS/aclocal.m4 index 6cec4be8..c6b83dae 100644 --- a/AMDiS/aclocal.m4 +++ b/AMDiS/aclocal.m4 @@ -1578,27 +1578,10 @@ linux*) # before this can be enabled. hardcode_into_libs=yes - # find out which ABI we are using - libsuff= - case "$host_cpu" in - x86_64*|s390x*|powerpc64*) - echo '[#]line __oline__ "configure"' > conftest.$ac_ext - if AC_TRY_EVAL(ac_compile); then - case `/usr/bin/file conftest.$ac_objext` in - *64-bit*) - libsuff=64 - sys_lib_search_path_spec="/lib${libsuff} /usr/lib${libsuff} /usr/local/lib${libsuff}" - ;; - esac - fi - rm -rf conftest* - ;; - esac - # Append ld.so.conf contents to the search path if test -f /etc/ld.so.conf; then - lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s 2>/dev/null", \[$]2)); skip = 1; } { if (!skip) print \[$]0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` - sys_lib_dlsearch_path_spec="/lib${libsuff} /usr/lib${libsuff} $lt_ld_extra" + lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s", \[$]2)); skip = 1; } { if (!skip) print \[$]0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` + sys_lib_dlsearch_path_spec="/lib /usr/lib $lt_ld_extra" fi # We used to test for /lib/ld.so.1 and disable shared libraries on @@ -4305,9 +4288,6 @@ CC=$lt_[]_LT_AC_TAGVAR(compiler, $1) # Is the compiler the GNU C compiler? with_gcc=$_LT_AC_TAGVAR(GCC, $1) -gcc_dir=\`gcc -print-file-name=. | $SED 's,/\.$,,'\` -gcc_ver=\`gcc -dumpversion\` - # An ERE matcher. EGREP=$lt_EGREP @@ -4441,11 +4421,11 @@ striplib=$lt_striplib # Dependencies to place before the objects being linked to create a # shared library. -predep_objects=\`echo $lt_[]_LT_AC_TAGVAR(predep_objects, $1) | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +predep_objects=$lt_[]_LT_AC_TAGVAR(predep_objects, $1) # Dependencies to place after the objects being linked to create a # shared library. -postdep_objects=\`echo $lt_[]_LT_AC_TAGVAR(postdep_objects, $1) | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +postdep_objects=$lt_[]_LT_AC_TAGVAR(postdep_objects, $1) # Dependencies to place before the objects being linked to create a # shared library. @@ -4457,7 +4437,7 @@ postdeps=$lt_[]_LT_AC_TAGVAR(postdeps, $1) # The library search path used internally by the compiler when linking # a shared library. -compiler_lib_search_path=\`echo $lt_[]_LT_AC_TAGVAR(compiler_lib_search_path, $1) | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +compiler_lib_search_path=$lt_[]_LT_AC_TAGVAR(compiler_lib_search_path, $1) # Method to check whether dependent libraries are shared objects. deplibs_check_method=$lt_deplibs_check_method @@ -4537,7 +4517,7 @@ variables_saved_for_relink="$variables_saved_for_relink" link_all_deplibs=$_LT_AC_TAGVAR(link_all_deplibs, $1) # Compile-time system search path for libraries -sys_lib_search_path_spec=\`echo $lt_sys_lib_search_path_spec | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +sys_lib_search_path_spec=$lt_sys_lib_search_path_spec # Run-time system search path for libraries sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec @@ -6373,7 +6353,6 @@ do done done done -IFS=$as_save_IFS lt_ac_max=0 lt_ac_count=0 # Add /usr/xpg4/bin/sed as it is typically found on Solaris @@ -6406,7 +6385,6 @@ for lt_ac_sed in $lt_ac_sed_list /usr/xpg4/bin/sed; do done ]) SED=$lt_cv_path_SED -AC_SUBST([SED]) AC_MSG_RESULT([$SED]) ]) diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in index 019b4cca..55098118 100644 --- a/AMDiS/bin/Makefile.in +++ b/AMDiS/bin/Makefile.in @@ -398,7 +398,6 @@ PACKAGE_VERSION = @PACKAGE_VERSION@ PATH_SEPARATOR = @PATH_SEPARATOR@ PETSC_DIR = @PETSC_DIR@ RANLIB = @RANLIB@ -SED = @SED@ SET_MAKE = @SET_MAKE@ SHELL = @SHELL@ STRIP = @STRIP@ diff --git a/AMDiS/configure b/AMDiS/configure index 5d40cff8..b4b86c5b 100755 --- a/AMDiS/configure +++ b/AMDiS/configure @@ -462,7 +462,7 @@ ac_includes_default="\ # include #endif" -ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CYGPATH_W PACKAGE VERSION ACLOCAL AUTOCONF AUTOMAKE AUTOHEADER MAKEINFO install_sh STRIP ac_ct_STRIP INSTALL_STRIP_PROGRAM mkdir_p AWK SET_MAKE am__leading_dot AMTAR am__tar am__untar MAINTAINER_MODE_TRUE MAINTAINER_MODE_FALSE MAINT AMDIS_DEBUG_TRUE AMDIS_DEBUG_FALSE CXX CC AMDIS_INTEL_TRUE AMDIS_INTEL_FALSE AMDIS_OPENMP_TRUE AMDIS_OPENMP_FALSE OPENMP_FLAG MPI_DIR PETSC_DIR USE_PARALLEL_AMDIS_TRUE USE_PARALLEL_AMDIS_FALSE USE_PARALLEL_DOMAIN_AMDIS_TRUE USE_PARALLEL_DOMAIN_AMDIS_FALSE ENABLE_UMFPACK_TRUE ENABLE_UMFPACK_FALSE ENABLE_MKL_TRUE ENABLE_MKL_FALSE DUNE_DIR ENABLE_DUNE_TRUE ENABLE_DUNE_FALSE ENABLE_BOOST_TRUE ENABLE_BOOST_FALSE CFLAGS LDFLAGS CPPFLAGS ac_ct_CC EXEEXT OBJEXT DEPDIR am__include am__quote AMDEP_TRUE AMDEP_FALSE AMDEPBACKSLASH CCDEPMODE am__fastdepCC_TRUE am__fastdepCC_FALSE CXXFLAGS ac_ct_CXX CXXDEPMODE am__fastdepCXX_TRUE am__fastdepCXX_FALSE build build_cpu build_vendor build_os host host_cpu host_vendor host_os SED EGREP LN_S ECHO AR ac_ct_AR RANLIB ac_ct_RANLIB CPP CXXCPP F77 FFLAGS ac_ct_F77 LIBTOOL LIBOBJS LTLIBOBJS' +ac_subst_vars='SHELL PATH_SEPARATOR PACKAGE_NAME PACKAGE_TARNAME PACKAGE_VERSION PACKAGE_STRING PACKAGE_BUGREPORT exec_prefix prefix program_transform_name bindir sbindir libexecdir datadir sysconfdir sharedstatedir localstatedir libdir includedir oldincludedir infodir mandir build_alias host_alias target_alias DEFS ECHO_C ECHO_N ECHO_T LIBS INSTALL_PROGRAM INSTALL_SCRIPT INSTALL_DATA CYGPATH_W PACKAGE VERSION ACLOCAL AUTOCONF AUTOMAKE AUTOHEADER MAKEINFO install_sh STRIP ac_ct_STRIP INSTALL_STRIP_PROGRAM mkdir_p AWK SET_MAKE am__leading_dot AMTAR am__tar am__untar MAINTAINER_MODE_TRUE MAINTAINER_MODE_FALSE MAINT AMDIS_DEBUG_TRUE AMDIS_DEBUG_FALSE CXX CC AMDIS_INTEL_TRUE AMDIS_INTEL_FALSE AMDIS_OPENMP_TRUE AMDIS_OPENMP_FALSE OPENMP_FLAG MPI_DIR PETSC_DIR USE_PARALLEL_AMDIS_TRUE USE_PARALLEL_AMDIS_FALSE USE_PARALLEL_DOMAIN_AMDIS_TRUE USE_PARALLEL_DOMAIN_AMDIS_FALSE ENABLE_UMFPACK_TRUE ENABLE_UMFPACK_FALSE ENABLE_MKL_TRUE ENABLE_MKL_FALSE DUNE_DIR ENABLE_DUNE_TRUE ENABLE_DUNE_FALSE ENABLE_BOOST_TRUE ENABLE_BOOST_FALSE CFLAGS LDFLAGS CPPFLAGS ac_ct_CC EXEEXT OBJEXT DEPDIR am__include am__quote AMDEP_TRUE AMDEP_FALSE AMDEPBACKSLASH CCDEPMODE am__fastdepCC_TRUE am__fastdepCC_FALSE CXXFLAGS ac_ct_CXX CXXDEPMODE am__fastdepCXX_TRUE am__fastdepCXX_FALSE build build_cpu build_vendor build_os host host_cpu host_vendor host_os EGREP LN_S ECHO AR ac_ct_AR RANLIB ac_ct_RANLIB CPP CXXCPP F77 FFLAGS ac_ct_F77 LIBTOOL LIBOBJS LTLIBOBJS' ac_subst_files='' # Initialize some variables set by options. @@ -3967,7 +3967,6 @@ do done done done -IFS=$as_save_IFS lt_ac_max=0 lt_ac_count=0 # Add /usr/xpg4/bin/sed as it is typically found on Solaris @@ -4002,7 +4001,6 @@ done fi SED=$lt_cv_path_SED - echo "$as_me:$LINENO: result: $SED" >&5 echo "${ECHO_T}$SED" >&6 @@ -4443,7 +4441,7 @@ ia64-*-hpux*) ;; *-*-irix6*) # Find out which ABI we are using. - echo '#line 4446 "configure"' > conftest.$ac_ext + echo '#line 4444 "configure"' > conftest.$ac_ext if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 (eval $ac_compile) 2>&5 ac_status=$? @@ -5578,7 +5576,7 @@ fi # Provide some information about the compiler. -echo "$as_me:5581:" \ +echo "$as_me:5579:" \ "checking for Fortran 77 compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` { (eval echo "$as_me:$LINENO: \"$ac_compiler --version &5\"") >&5 @@ -6641,11 +6639,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:6644: $lt_compile\"" >&5) + (eval echo "\"\$as_me:6642: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:6648: \$? = $ac_status" >&5 + echo "$as_me:6646: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings other than the usual output. @@ -6909,11 +6907,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:6912: $lt_compile\"" >&5) + (eval echo "\"\$as_me:6910: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:6916: \$? = $ac_status" >&5 + echo "$as_me:6914: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings other than the usual output. @@ -7013,11 +7011,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:7016: $lt_compile\"" >&5) + (eval echo "\"\$as_me:7014: $lt_compile\"" >&5) (eval "$lt_compile" 2>out/conftest.err) ac_status=$? cat out/conftest.err >&5 - echo "$as_me:7020: \$? = $ac_status" >&5 + echo "$as_me:7018: \$? = $ac_status" >&5 if (exit $ac_status) && test -s out/conftest2.$ac_objext then # The compiler can only warn and ignore the option if not recognized @@ -8478,31 +8476,10 @@ linux*) # before this can be enabled. hardcode_into_libs=yes - # find out which ABI we are using - libsuff= - case "$host_cpu" in - x86_64*|s390x*|powerpc64*) - echo '#line 8485 "configure"' > conftest.$ac_ext - if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 - (eval $ac_compile) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; then - case `/usr/bin/file conftest.$ac_objext` in - *64-bit*) - libsuff=64 - sys_lib_search_path_spec="/lib${libsuff} /usr/lib${libsuff} /usr/local/lib${libsuff}" - ;; - esac - fi - rm -rf conftest* - ;; - esac - # Append ld.so.conf contents to the search path if test -f /etc/ld.so.conf; then - lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s 2>/dev/null", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` - sys_lib_dlsearch_path_spec="/lib${libsuff} /usr/lib${libsuff} $lt_ld_extra" + lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` + sys_lib_dlsearch_path_spec="/lib /usr/lib $lt_ld_extra" fi # We used to test for /lib/ld.so.1 and disable shared libraries on @@ -9379,7 +9356,7 @@ else lt_dlunknown=0; lt_dlno_uscore=1; lt_dlneed_uscore=2 lt_status=$lt_dlunknown cat > conftest.$ac_ext < conftest.$ac_ext <&5) + (eval echo "\"\$as_me:11799: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:11829: \$? = $ac_status" >&5 + echo "$as_me:11803: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings other than the usual output. @@ -11926,11 +11900,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:11929: $lt_compile\"" >&5) + (eval echo "\"\$as_me:11903: $lt_compile\"" >&5) (eval "$lt_compile" 2>out/conftest.err) ac_status=$? cat out/conftest.err >&5 - echo "$as_me:11933: \$? = $ac_status" >&5 + echo "$as_me:11907: \$? = $ac_status" >&5 if (exit $ac_status) && test -s out/conftest2.$ac_objext then # The compiler can only warn and ignore the option if not recognized @@ -12458,31 +12432,10 @@ linux*) # before this can be enabled. hardcode_into_libs=yes - # find out which ABI we are using - libsuff= - case "$host_cpu" in - x86_64*|s390x*|powerpc64*) - echo '#line 12465 "configure"' > conftest.$ac_ext - if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 - (eval $ac_compile) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; then - case `/usr/bin/file conftest.$ac_objext` in - *64-bit*) - libsuff=64 - sys_lib_search_path_spec="/lib${libsuff} /usr/lib${libsuff} /usr/local/lib${libsuff}" - ;; - esac - fi - rm -rf conftest* - ;; - esac - # Append ld.so.conf contents to the search path if test -f /etc/ld.so.conf; then - lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s 2>/dev/null", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` - sys_lib_dlsearch_path_spec="/lib${libsuff} /usr/lib${libsuff} $lt_ld_extra" + lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` + sys_lib_dlsearch_path_spec="/lib /usr/lib $lt_ld_extra" fi # We used to test for /lib/ld.so.1 and disable shared libraries on @@ -12866,9 +12819,6 @@ CC=$lt_compiler_CXX # Is the compiler the GNU C compiler? with_gcc=$GCC_CXX -gcc_dir=\`gcc -print-file-name=. | $SED 's,/\.$,,'\` -gcc_ver=\`gcc -dumpversion\` - # An ERE matcher. EGREP=$lt_EGREP @@ -13002,11 +12952,11 @@ striplib=$lt_striplib # Dependencies to place before the objects being linked to create a # shared library. -predep_objects=\`echo $lt_predep_objects_CXX | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +predep_objects=$lt_predep_objects_CXX # Dependencies to place after the objects being linked to create a # shared library. -postdep_objects=\`echo $lt_postdep_objects_CXX | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +postdep_objects=$lt_postdep_objects_CXX # Dependencies to place before the objects being linked to create a # shared library. @@ -13018,7 +12968,7 @@ postdeps=$lt_postdeps_CXX # The library search path used internally by the compiler when linking # a shared library. -compiler_lib_search_path=\`echo $lt_compiler_lib_search_path_CXX | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +compiler_lib_search_path=$lt_compiler_lib_search_path_CXX # Method to check whether dependent libraries are shared objects. deplibs_check_method=$lt_deplibs_check_method @@ -13098,7 +13048,7 @@ variables_saved_for_relink="$variables_saved_for_relink" link_all_deplibs=$link_all_deplibs_CXX # Compile-time system search path for libraries -sys_lib_search_path_spec=\`echo $lt_sys_lib_search_path_spec | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +sys_lib_search_path_spec=$lt_sys_lib_search_path_spec # Run-time system search path for libraries sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec @@ -13520,11 +13470,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:13523: $lt_compile\"" >&5) + (eval echo "\"\$as_me:13473: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:13527: \$? = $ac_status" >&5 + echo "$as_me:13477: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings other than the usual output. @@ -13624,11 +13574,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:13627: $lt_compile\"" >&5) + (eval echo "\"\$as_me:13577: $lt_compile\"" >&5) (eval "$lt_compile" 2>out/conftest.err) ac_status=$? cat out/conftest.err >&5 - echo "$as_me:13631: \$? = $ac_status" >&5 + echo "$as_me:13581: \$? = $ac_status" >&5 if (exit $ac_status) && test -s out/conftest2.$ac_objext then # The compiler can only warn and ignore the option if not recognized @@ -15069,31 +15019,10 @@ linux*) # before this can be enabled. hardcode_into_libs=yes - # find out which ABI we are using - libsuff= - case "$host_cpu" in - x86_64*|s390x*|powerpc64*) - echo '#line 15076 "configure"' > conftest.$ac_ext - if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 - (eval $ac_compile) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; then - case `/usr/bin/file conftest.$ac_objext` in - *64-bit*) - libsuff=64 - sys_lib_search_path_spec="/lib${libsuff} /usr/lib${libsuff} /usr/local/lib${libsuff}" - ;; - esac - fi - rm -rf conftest* - ;; - esac - # Append ld.so.conf contents to the search path if test -f /etc/ld.so.conf; then - lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s 2>/dev/null", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` - sys_lib_dlsearch_path_spec="/lib${libsuff} /usr/lib${libsuff} $lt_ld_extra" + lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` + sys_lib_dlsearch_path_spec="/lib /usr/lib $lt_ld_extra" fi # We used to test for /lib/ld.so.1 and disable shared libraries on @@ -15477,9 +15406,6 @@ CC=$lt_compiler_F77 # Is the compiler the GNU C compiler? with_gcc=$GCC_F77 -gcc_dir=\`gcc -print-file-name=. | $SED 's,/\.$,,'\` -gcc_ver=\`gcc -dumpversion\` - # An ERE matcher. EGREP=$lt_EGREP @@ -15613,11 +15539,11 @@ striplib=$lt_striplib # Dependencies to place before the objects being linked to create a # shared library. -predep_objects=\`echo $lt_predep_objects_F77 | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +predep_objects=$lt_predep_objects_F77 # Dependencies to place after the objects being linked to create a # shared library. -postdep_objects=\`echo $lt_postdep_objects_F77 | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +postdep_objects=$lt_postdep_objects_F77 # Dependencies to place before the objects being linked to create a # shared library. @@ -15629,7 +15555,7 @@ postdeps=$lt_postdeps_F77 # The library search path used internally by the compiler when linking # a shared library. -compiler_lib_search_path=\`echo $lt_compiler_lib_search_path_F77 | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +compiler_lib_search_path=$lt_compiler_lib_search_path_F77 # Method to check whether dependent libraries are shared objects. deplibs_check_method=$lt_deplibs_check_method @@ -15709,7 +15635,7 @@ variables_saved_for_relink="$variables_saved_for_relink" link_all_deplibs=$link_all_deplibs_F77 # Compile-time system search path for libraries -sys_lib_search_path_spec=\`echo $lt_sys_lib_search_path_spec | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +sys_lib_search_path_spec=$lt_sys_lib_search_path_spec # Run-time system search path for libraries sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec @@ -15851,11 +15777,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:15854: $lt_compile\"" >&5) + (eval echo "\"\$as_me:15780: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:15858: \$? = $ac_status" >&5 + echo "$as_me:15784: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings other than the usual output. @@ -16119,11 +16045,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:16122: $lt_compile\"" >&5) + (eval echo "\"\$as_me:16048: $lt_compile\"" >&5) (eval "$lt_compile" 2>conftest.err) ac_status=$? cat conftest.err >&5 - echo "$as_me:16126: \$? = $ac_status" >&5 + echo "$as_me:16052: \$? = $ac_status" >&5 if (exit $ac_status) && test -s "$ac_outfile"; then # The compiler can only warn and ignore the option if not recognized # So say no if there are warnings other than the usual output. @@ -16223,11 +16149,11 @@ else -e 's:.*FLAGS}\{0,1\} :&$lt_compiler_flag :; t' \ -e 's: [^ ]*conftest\.: $lt_compiler_flag&:; t' \ -e 's:$: $lt_compiler_flag:'` - (eval echo "\"\$as_me:16226: $lt_compile\"" >&5) + (eval echo "\"\$as_me:16152: $lt_compile\"" >&5) (eval "$lt_compile" 2>out/conftest.err) ac_status=$? cat out/conftest.err >&5 - echo "$as_me:16230: \$? = $ac_status" >&5 + echo "$as_me:16156: \$? = $ac_status" >&5 if (exit $ac_status) && test -s out/conftest2.$ac_objext then # The compiler can only warn and ignore the option if not recognized @@ -17688,31 +17614,10 @@ linux*) # before this can be enabled. hardcode_into_libs=yes - # find out which ABI we are using - libsuff= - case "$host_cpu" in - x86_64*|s390x*|powerpc64*) - echo '#line 17695 "configure"' > conftest.$ac_ext - if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5 - (eval $ac_compile) 2>&5 - ac_status=$? - echo "$as_me:$LINENO: \$? = $ac_status" >&5 - (exit $ac_status); }; then - case `/usr/bin/file conftest.$ac_objext` in - *64-bit*) - libsuff=64 - sys_lib_search_path_spec="/lib${libsuff} /usr/lib${libsuff} /usr/local/lib${libsuff}" - ;; - esac - fi - rm -rf conftest* - ;; - esac - # Append ld.so.conf contents to the search path if test -f /etc/ld.so.conf; then - lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s 2>/dev/null", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` - sys_lib_dlsearch_path_spec="/lib${libsuff} /usr/lib${libsuff} $lt_ld_extra" + lt_ld_extra=`awk '/^include / { system(sprintf("cd /etc; cat %s", \$2)); skip = 1; } { if (!skip) print \$0; skip = 0; }' < /etc/ld.so.conf | $SED -e 's/#.*//;s/[:, ]/ /g;s/=[^=]*$//;s/=[^= ]* / /g;/^$/d' | tr '\n' ' '` + sys_lib_dlsearch_path_spec="/lib /usr/lib $lt_ld_extra" fi # We used to test for /lib/ld.so.1 and disable shared libraries on @@ -18096,9 +18001,6 @@ CC=$lt_compiler_GCJ # Is the compiler the GNU C compiler? with_gcc=$GCC_GCJ -gcc_dir=\`gcc -print-file-name=. | $SED 's,/\.$,,'\` -gcc_ver=\`gcc -dumpversion\` - # An ERE matcher. EGREP=$lt_EGREP @@ -18232,11 +18134,11 @@ striplib=$lt_striplib # Dependencies to place before the objects being linked to create a # shared library. -predep_objects=\`echo $lt_predep_objects_GCJ | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +predep_objects=$lt_predep_objects_GCJ # Dependencies to place after the objects being linked to create a # shared library. -postdep_objects=\`echo $lt_postdep_objects_GCJ | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +postdep_objects=$lt_postdep_objects_GCJ # Dependencies to place before the objects being linked to create a # shared library. @@ -18248,7 +18150,7 @@ postdeps=$lt_postdeps_GCJ # The library search path used internally by the compiler when linking # a shared library. -compiler_lib_search_path=\`echo $lt_compiler_lib_search_path_GCJ | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +compiler_lib_search_path=$lt_compiler_lib_search_path_GCJ # Method to check whether dependent libraries are shared objects. deplibs_check_method=$lt_deplibs_check_method @@ -18328,7 +18230,7 @@ variables_saved_for_relink="$variables_saved_for_relink" link_all_deplibs=$link_all_deplibs_GCJ # Compile-time system search path for libraries -sys_lib_search_path_spec=\`echo $lt_sys_lib_search_path_spec | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +sys_lib_search_path_spec=$lt_sys_lib_search_path_spec # Run-time system search path for libraries sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec @@ -18580,9 +18482,6 @@ CC=$lt_compiler_RC # Is the compiler the GNU C compiler? with_gcc=$GCC_RC -gcc_dir=\`gcc -print-file-name=. | $SED 's,/\.$,,'\` -gcc_ver=\`gcc -dumpversion\` - # An ERE matcher. EGREP=$lt_EGREP @@ -18716,11 +18615,11 @@ striplib=$lt_striplib # Dependencies to place before the objects being linked to create a # shared library. -predep_objects=\`echo $lt_predep_objects_RC | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +predep_objects=$lt_predep_objects_RC # Dependencies to place after the objects being linked to create a # shared library. -postdep_objects=\`echo $lt_postdep_objects_RC | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +postdep_objects=$lt_postdep_objects_RC # Dependencies to place before the objects being linked to create a # shared library. @@ -18732,7 +18631,7 @@ postdeps=$lt_postdeps_RC # The library search path used internally by the compiler when linking # a shared library. -compiler_lib_search_path=\`echo $lt_compiler_lib_search_path_RC | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +compiler_lib_search_path=$lt_compiler_lib_search_path_RC # Method to check whether dependent libraries are shared objects. deplibs_check_method=$lt_deplibs_check_method @@ -18812,7 +18711,7 @@ variables_saved_for_relink="$variables_saved_for_relink" link_all_deplibs=$link_all_deplibs_RC # Compile-time system search path for libraries -sys_lib_search_path_spec=\`echo $lt_sys_lib_search_path_spec | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\` +sys_lib_search_path_spec=$lt_sys_lib_search_path_spec # Run-time system search path for libraries sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec @@ -19736,7 +19635,6 @@ s,@host@,$host,;t t s,@host_cpu@,$host_cpu,;t t s,@host_vendor@,$host_vendor,;t t s,@host_os@,$host_os,;t t -s,@SED@,$SED,;t t s,@EGREP@,$EGREP,;t t s,@LN_S@,$LN_S,;t t s,@ECHO@,$ECHO,;t t diff --git a/AMDiS/other/include/Makefile_AMDiS.mk b/AMDiS/other/include/Makefile_AMDiS.mk index 5df5370e..59908071 100644 --- a/AMDiS/other/include/Makefile_AMDiS.mk +++ b/AMDiS/other/include/Makefile_AMDiS.mk @@ -79,6 +79,7 @@ ifeq ($(strip $(USE_PARALLEL_AMDIS)), 1) endif LIBS += $(PARMETIS_LIB) -lmpi + CPPFLAGS += -DHAVE_PARALLEL_DOMAIN_AMDIS else ifeq ($(strip $(USE_COMPILER)), gcc) COMPILE = g++ @@ -96,9 +97,9 @@ endif # ============================================================================ ifeq ($(strip $(DEBUG)), 0) - CPPFLAGS = -O2 + CPPFLAGS += -O2 else - CPPFLAGS = -g -O0 + CPPFLAGS += -g -O0 endif ifeq ($(strip $(USE_OPENMP)), 1) diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h index 64e103d2..6c55e9fb 100644 --- a/AMDiS/src/AMDiS_fwd.h +++ b/AMDiS/src/AMDiS_fwd.h @@ -51,6 +51,7 @@ namespace AMDiS { class FiniteElemSpace; class Flag; class IdentityPreconditioner; + class InteriorBoundary; class ITL_BasePreconditioner; class LeafDataPeriodic; class LevelAdmin; @@ -94,6 +95,9 @@ namespace AMDiS { class VertexInfo; class VertexVector; + struct BoundaryObject; + struct AtomicBoundary; + template class AbstractFunction; template class DOFIndexed; template class DOFVectorBase; diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h index 3eb7b05e..6203bb6c 100644 --- a/AMDiS/src/BasisFunction.h +++ b/AMDiS/src/BasisFunction.h @@ -49,8 +49,7 @@ namespace AMDiS { virtual ~GrdBasFctType() {} - virtual void operator()(const DimVec&, - DimVec&) const = 0; + virtual void operator()(const DimVec&, DimVec&) const = 0; }; /// Function interface for evaluating second derivative of basis functions. @@ -61,8 +60,7 @@ namespace AMDiS { virtual ~D2BasFctType() {} - virtual void operator()(const DimVec&, - DimMat&) const = 0; + virtual void operator()(const DimVec&, DimMat&) const = 0; }; typedef BasFctType *BFptr; @@ -297,12 +295,6 @@ namespace AMDiS { getLocalIndices(el, admin, &(indices[0])); } - /// Returns local dof indices of the element for the given fe space. - virtual void getLocalIndicesVec(const Element *el, - const DOFAdmin *admin, - Vector *ve) const - {} - virtual void getLocalDofPtrVec(const Element *el, const DOFAdmin *admin, std::vector& vec) const diff --git a/AMDiS/src/BoundaryManager.cc b/AMDiS/src/BoundaryManager.cc index 0aa2e763..c36f6c31 100644 --- a/AMDiS/src/BoundaryManager.cc +++ b/AMDiS/src/BoundaryManager.cc @@ -54,7 +54,7 @@ namespace AMDiS { { if (localBCs.size() > 0) { const FiniteElemSpace *feSpace = vec->getFESpace(); - Vector &dofVec = dofIndices[omp_get_thread_num()]; + std::vector &dofVec = dofIndices[omp_get_thread_num()]; const BasisFunction *basisFcts = feSpace->getBasisFcts(); int nBasFcts = basisFcts->getNumber(); @@ -63,7 +63,7 @@ namespace AMDiS { basisFcts->getBound(elInfo, localBound); // get dof indices - basisFcts->getLocalIndicesVec(elInfo->getElement(), feSpace->getAdmin(), &dofVec); + basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec); // apply non dirichlet boundary conditions for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it) @@ -85,7 +85,7 @@ namespace AMDiS { return; const FiniteElemSpace *feSpace = mat->getRowFESpace(); - Vector &dofVec = dofIndices[omp_get_thread_num()]; + std::vector &dofVec = dofIndices[omp_get_thread_num()]; const BasisFunction *basisFcts = feSpace->getBasisFcts(); int nBasFcts = basisFcts->getNumber(); @@ -94,8 +94,7 @@ namespace AMDiS { basisFcts->getBound(elInfo, localBound); // get dof indices - basisFcts->getLocalIndicesVec(elInfo->getElement(), - feSpace->getAdmin(), &dofVec); + basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec); // apply non dirichlet boundary conditions for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it) diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h index a68aa641..c9f391a0 100644 --- a/AMDiS/src/BoundaryManager.h +++ b/AMDiS/src/BoundaryManager.h @@ -126,7 +126,7 @@ namespace AMDiS { std::vector localBounds; /// Temporary thread-safe variable for functions fillBoundaryconditions. - std::vector > dofIndices; + std::vector > dofIndices; /** \brief * Stores the number of byte that were allocated in the constructor for diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc index d3007521..9e8ddd69 100755 --- a/AMDiS/src/DOFAdmin.cc +++ b/AMDiS/src/DOFAdmin.cc @@ -139,7 +139,7 @@ namespace AMDiS { firstHole = i; } else { // if there is no hole // enlarge dof-list - enlargeDOFLists(0); + enlargeDOFLists(); TEST_EXIT_DBG(firstHole < static_cast(dofFree.size())) ("no free entry after enlargeDOFLists\n"); diff --git a/AMDiS/src/DOFAdmin.h b/AMDiS/src/DOFAdmin.h index 098251aa..8214f3bd 100644 --- a/AMDiS/src/DOFAdmin.h +++ b/AMDiS/src/DOFAdmin.h @@ -66,7 +66,7 @@ namespace AMDiS { * Enlarges the number of DOFs that can be managed at least to minsize by * a step size of \ref sizeIncrement. */ - void enlargeDOFLists(int minsize); + void enlargeDOFLists(int minsize = 0); /// assignment operator DOFAdmin& operator=(const DOFAdmin&); diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh index 929946ea..141a3dc6 100644 --- a/AMDiS/src/DOFVector.hh +++ b/AMDiS/src/DOFVector.hh @@ -147,10 +147,9 @@ namespace AMDiS { { FUNCNAME("DOFVector::addElementVector()"); - Vector indices(nBasFcts); - feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), - feSpace->getAdmin(), - &indices); + std::vector indices(nBasFcts); + feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), + indices); for (DegreeOfFreedom i = 0; i < nBasFcts; i++) { BoundaryCondition *condition = @@ -159,11 +158,10 @@ namespace AMDiS { if (!(condition && condition->isDirichlet())) { DegreeOfFreedom irow = indices[i]; - if (add) { + if (add) (*this)[irow] += factor * elVec[i]; - } else { - (*this)[irow] = factor * elVec[i]; - } + else + (*this)[irow] = factor * elVec[i]; } } } diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc index c47f82d4..e251c974 100644 --- a/AMDiS/src/Element.cc +++ b/AMDiS/src/Element.cc @@ -566,30 +566,71 @@ namespace AMDiS { return result; } - void fitElementToMeshCode(RefinementManager *refineManager, MeshStructure &code, - Element *el, int ithSide, int elType) + bool fitElementToMeshCode(RefinementManager *refineManager, + MeshStructure &code, + Element *el, + int ithSide, + int elType) + { + FUNCNAME("fitElementToMeshCode()"); + + if (code.empty()) + return false; + + int s1 = el->getSideOfChild(0, ithSide, elType); + int s2 = el->getSideOfChild(1, ithSide, elType); + + TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n"); + + if (s1 != -1 && s2 != -1) + return fitElementToMeshCode2(refineManager, code, el, ithSide, elType); + + if (el->isLeaf()) { + if (code.getNumElements() == 1 && code.isLeafElement()) + return false; + + el->setMark(1); + refineManager->refineElement(el->getMesh(), el); + } + + if (s1 != -1) + return fitElementToMeshCode2(refineManager, code, + el->getFirstChild(), s1, el->getChildType(elType)); + else + return fitElementToMeshCode2(refineManager, code, + el->getSecondChild(), s2, el->getChildType(elType)); + } + + bool fitElementToMeshCode2(RefinementManager *refineManager, MeshStructure &code, + Element *el, int ithSide, int elType) { if (code.isLeafElement()) - return; + return false; + + bool value = false; if (el->isLeaf()) { el->setMark(1); - refineManager->refineMesh(el->getMesh()); + refineManager->refineElement(el->getMesh(), el); + value = true; } int s1 = el->getSideOfChild(0, ithSide, elType); int s2 = el->getSideOfChild(1, ithSide, elType); - + if (s1 != -1) { code.nextElement(); - fitElementToMeshCode(refineManager, code, - el->getFirstChild(), s1, el->getChildType(elType)); + value |= fitElementToMeshCode2(refineManager, code, el->getFirstChild(), + s1, el->getChildType(elType)); } + if (s2 != -1) { code.nextElement(); - fitElementToMeshCode(refineManager, code, - el->getSecondChild(), s2, el->getChildType(elType)); + value |= fitElementToMeshCode2(refineManager, code, el->getSecondChild(), + s2, el->getChildType(elType)); } + + return value; } } diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h index 9f9d60d4..6de14fc4 100644 --- a/AMDiS/src/Element.h +++ b/AMDiS/src/Element.h @@ -393,7 +393,9 @@ namespace AMDiS { */ virtual void getVertexDofs(FiniteElemSpace* feSpace, int ith, int elType, - DofContainer& dofs, bool parentVertices = 0) const = 0; + DofContainer& dofs, + bool reverseMode, + bool parentVertices = false) const = 0; /** \brief * Traverses an edge/face of a given element (this includes also all children of @@ -574,11 +576,17 @@ namespace AMDiS { friend class Mesh; }; - void fitElementToMeshCode(RefinementManager *refineManager, + bool fitElementToMeshCode(RefinementManager *refineManager, MeshStructure &code, Element *el, int ithSide, int elType); + + bool fitElementToMeshCode2(RefinementManager *refineManager, + MeshStructure &code, + Element *el, + int ithSide, + int elType); } #endif // AMDIS_ELEMENT_H diff --git a/AMDiS/src/Global.cc b/AMDiS/src/Global.cc index b23e86db..e5c02d9c 100644 --- a/AMDiS/src/Global.cc +++ b/AMDiS/src/Global.cc @@ -188,8 +188,9 @@ namespace AMDiS { oss << "WARNING in " << file << ", line " << line << std::endl; } - PRINT_LINE((*out), oss.str()); - + if (oss.str() != "") + PRINT_LINE((*out), oss.str()); + old_line = line; } @@ -283,6 +284,12 @@ namespace AMDiS { return i * fac(i - 1); } + void waitSec(int seconds) + { + clock_t endwait = clock () + seconds * CLOCKS_PER_SEC; + while (clock() < endwait) {} + } + std::string memSizeStr(int size) { std::string result; diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h index 2c4f8aa8..21a9ad77 100644 --- a/AMDiS/src/Global.h +++ b/AMDiS/src/Global.h @@ -46,6 +46,7 @@ #include #include #include +#include #if HAVE_PARALLEL_DOMAIN_AMDIS #include "mpi.h" @@ -92,6 +93,8 @@ namespace AMDiS { /// Calculates factorial of i int fac(int i); + void waitSec(int seconds); + /// Content comparision of two pointers. Used e.g. for std::find_if template struct comparePtrContents : public std::binary_function @@ -314,7 +317,6 @@ namespace AMDiS { */ #define WAIT_REALLY Msg::wait(true) -#include #define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC) /** \brief @@ -443,12 +445,8 @@ namespace AMDiS { GRD_PHI }; -#define MEMORY_MANAGED(className) ; #define NEW new #define DELETE delete -#define GET_MEMORY(typename, number) new typename[number] -#define FREE_MEMORY(ptr, typename, number) delete [] ptr - } #endif // AMDIS_GLOBAL_H diff --git a/AMDiS/src/InteriorBoundary.cc b/AMDiS/src/InteriorBoundary.cc index d751a1d0..aa0a1b2f 100644 --- a/AMDiS/src/InteriorBoundary.cc +++ b/AMDiS/src/InteriorBoundary.cc @@ -1,14 +1,48 @@ #include "InteriorBoundary.h" +#include "FiniteElemSpace.h" +#include "BasisFunction.h" #include "Serializer.h" namespace AMDiS { + void BoundaryObject::setReverseMode(BoundaryObject &otherBound, + FiniteElemSpace *feSpace) + { + FUNCNAME("BoundaryObject::setReverseMode()"); + + bool otherMode = false; + + const BasisFunction *basFcts = feSpace->getBasisFcts(); + int nBasFcts = basFcts->getNumber(); + std::vector localDofs0(nBasFcts), localDofs1(nBasFcts); + + switch (feSpace->getMesh()->getDim()) { + case 2: + ERROR_EXIT("Not yet implemented!\n"); + break; + case 3: + if (ithObj == 2 || ithObj == 3) { + basFcts->getLocalIndices(el, feSpace->getAdmin(), localDofs0); + basFcts->getLocalIndices(otherBound.el, feSpace->getAdmin(), localDofs1); + otherMode = (localDofs0[0] != localDofs1[0]); + } + + break; + default: + ERROR_EXIT("This should not happen!\n"); + } + + otherBound.reverseMode = otherMode; + } + + AtomicBoundary& InteriorBoundary::getNewAtomic(int rank) { boundary[rank].resize(boundary[rank].size() + 1); return boundary[rank][boundary[rank].size() - 1]; } + void InteriorBoundary::serialize(std::ostream &out) { int mSize = boundary.size(); @@ -22,16 +56,21 @@ namespace AMDiS { AtomicBoundary &bound = (it->second)[i]; SerUtil::serialize(out, bound.rankObj.elIndex); + SerUtil::serialize(out, bound.rankObj.elType); SerUtil::serialize(out, bound.rankObj.subObj); SerUtil::serialize(out, bound.rankObj.ithObj); + SerUtil::serialize(out, bound.rankObj.reverseMode); SerUtil::serialize(out, bound.neighObj.elIndex); + SerUtil::serialize(out, bound.neighObj.elType); SerUtil::serialize(out, bound.neighObj.subObj); SerUtil::serialize(out, bound.neighObj.ithObj); + SerUtil::serialize(out, bound.neighObj.reverseMode); } } } + void InteriorBoundary::deserialize(std::istream &in, std::map &elIndexMap) { @@ -48,12 +87,16 @@ namespace AMDiS { AtomicBoundary &bound = boundary[rank][i]; SerUtil::deserialize(in, bound.rankObj.elIndex); + SerUtil::deserialize(in, bound.rankObj.elType); SerUtil::deserialize(in, bound.rankObj.subObj); SerUtil::deserialize(in, bound.rankObj.ithObj); + SerUtil::deserialize(in, bound.rankObj.reverseMode); SerUtil::deserialize(in, bound.neighObj.elIndex); + SerUtil::deserialize(in, bound.neighObj.elType); SerUtil::deserialize(in, bound.neighObj.subObj); SerUtil::deserialize(in, bound.neighObj.ithObj); + SerUtil::deserialize(in, bound.neighObj.reverseMode); bound.rankObj.el = elIndexMap[bound.rankObj.elIndex]; bound.neighObj.el = NULL; diff --git a/AMDiS/src/InteriorBoundary.h b/AMDiS/src/InteriorBoundary.h index cb12f255..b155b473 100644 --- a/AMDiS/src/InteriorBoundary.h +++ b/AMDiS/src/InteriorBoundary.h @@ -32,6 +32,9 @@ namespace AMDiS { /// Defines the geometrical objects that forms the boundary; struct BoundaryObject { + BoundaryObject() + : reverseMode(false) + {} /// The macro element to which the boundary element corresponds to. Element* el; @@ -58,6 +61,10 @@ namespace AMDiS { * boundary. */ int ithObj; + + bool reverseMode; + + void setReverseMode(BoundaryObject &otherBound, FiniteElemSpace *feSpace); }; /** \brief diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc index f3fe9d50..847b8cd3 100644 --- a/AMDiS/src/Lagrange.cc +++ b/AMDiS/src/Lagrange.cc @@ -712,7 +712,7 @@ namespace AMDiS { return &sortedVertex; // edge - if ((dimOfPosition == 1) && (degree == 2)) + if (dimOfPosition == 1 && degree == 2) return &sortedEdgeDeg2; int vertex[3]; @@ -762,7 +762,7 @@ namespace AMDiS { } // center - if ((dimOfPosition == 3) && (degree == 4)) + if (dimOfPosition == 3 && degree == 4) return &sortedCenterDeg4; ERROR_EXIT("should not be reached\n"); @@ -953,16 +953,6 @@ namespace AMDiS { return result; } - void Lagrange::getLocalIndicesVec(const Element* el, - const DOFAdmin *admin, - Vector *indices) const - { - if (indices->getSize() < nBasFcts) - indices->resize(nBasFcts); - - getLocalIndices(el, admin, &((*indices)[0])); - } - void Lagrange::getLocalDofPtrVec(const Element *el, const DOFAdmin *admin, std::vector& vec) const diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h index abb7bfe2..f4a7bddb 100644 --- a/AMDiS/src/Lagrange.h +++ b/AMDiS/src/Lagrange.h @@ -123,11 +123,6 @@ namespace AMDiS { const DOFAdmin *admin, DegreeOfFreedom *dofs) const; - /// - void getLocalIndicesVec(const Element *el, - const DOFAdmin *admin, - Vector *vec) const; - void getLocalDofPtrVec(const Element *el, const DOFAdmin *admin, std::vector& vec) const; diff --git a/AMDiS/src/Line.h b/AMDiS/src/Line.h index a32573f1..9eed856a 100644 --- a/AMDiS/src/Line.h +++ b/AMDiS/src/Line.h @@ -154,7 +154,9 @@ namespace AMDiS { } void getVertexDofs(FiniteElemSpace* feSpace, int ithEdge, int elType, - DofContainer& dofs, bool parentVertices = 0) const + DofContainer& dofs, + bool reverseMode, + bool parentVertices = false) const { FUNCNAME("Line::getVertexDofs()"); diff --git a/AMDiS/src/MacroElement.cc b/AMDiS/src/MacroElement.cc index 507002ad..49cef625 100644 --- a/AMDiS/src/MacroElement.cc +++ b/AMDiS/src/MacroElement.cc @@ -28,7 +28,7 @@ namespace AMDiS { MacroElement::~MacroElement() { if (element) - delete element; + delete element; } MacroElement& MacroElement::operator=(const MacroElement &el) diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index c6ad8174..059afb14 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -262,7 +262,7 @@ namespace AMDiS { me->setIndex(macroElements.size()); } - void Mesh::removeMacroElements(std::vector& macros, + void Mesh::removeMacroElements(std::set& macros, const FiniteElemSpace *feSpace) { FUNCNAME("Mesh::removeMacroElement()"); @@ -276,11 +276,13 @@ namespace AMDiS { ElementDofIterator elDofIter(feSpace); // Map that stores for each dof pointer (which may have a list of dofs) - // all macro element indices that own the dof. + // all macro element indices that own this dof. DofElMap dofsOwner; DofPosMap dofsPosIndex; - - // Determine all dof owner macro elements. + + + // === Determine to all dofs the macro elements where the dof is part of. === + for (std::deque::iterator macroIt = macroElements.begin(); macroIt != macroElements.end(); ++macroIt) { elDofIter.reset((*macroIt)->getElement()); @@ -290,16 +292,31 @@ namespace AMDiS { } while (elDofIter.nextStrict()); } - // Remove all the given macro elements. - for (std::vector::iterator macroIt = macros.begin(); - macroIt != macros.end(); - ++macroIt) { + // === Remove macro elements from mesh macro element list. === - std::deque::iterator mEl = - find(macroElements.begin(), macroElements.end(), *macroIt); - TEST_EXIT(mEl != macroElements.end()) - ("Cannot find MacroElement that should be removed!\n"); - macroElements.erase(mEl); + // Removing arbitrary elements from an std::deque is very slow. Therefore, we + // create a new deque with all macro elements that should not be deleted. The + // macro element deque is than replaced by the new created one. + + std::deque newMacroElements; + for (std::deque::iterator elIter = macroElements.begin(); + elIter != macroElements.end(); ++elIter) { + // If the current mesh macro element should not be deleted, i.e., it is not a + // member of the list of macro elements to be deleted, is is inserted to the new + // macro element list. + if (macros.find(*elIter) == macros.end()) + newMacroElements.push_back(*elIter); + } + // And replace the macro element list with the new one. + macroElements.clear(); + macroElements = newMacroElements; + + + // === For all macro elements to be deleted, delete them also to be neighbours === + // === of some other macro elements. === + + for (std::set::iterator macroIt = macros.begin(); + macroIt != macros.end(); ++macroIt) { // Go through all neighbours of the macro element and remove this macro element // to be neighbour of some other macro element. @@ -315,33 +332,44 @@ namespace AMDiS { } } + // Decrease also the number of elements of the mesh. nLeaves--; nElements--; + } - // Remove this macro element from the dof owner list. - for (DofElMap::iterator dofsIt = dofsOwner.begin(); - dofsIt != dofsOwner.end(); ++dofsIt) { - std::set::iterator mIt = dofsIt->second.find(*macroIt); - if (mIt != dofsIt->second.end()) - dofsIt->second.erase(mIt); - } - // And remove the macro element from memory - delete *macroIt; - } + // === Check now all the dofs, that have no owner anymore and therefore have === + // === to be removed. === int nRemainDofs = 0; - // Check now all the dofs, that have no owner anymore and therefore have to - // be removed. for (DofElMap::iterator dofsIt = dofsOwner.begin(); - dofsIt != dofsOwner.end(); ++dofsIt) { - if (dofsIt->second.size() == 0) + dofsIt != dofsOwner.end(); ++dofsIt) { + + bool deleteDof = true; + + for (std::set::iterator elIter = dofsIt->second.begin(); + elIter != dofsIt->second.end(); ++elIter) { + std::set::iterator mIt = macros.find(*elIter); + if (mIt == macros.end()) { + deleteDof = false; + break; + } + } + + if (deleteDof) freeDOF(const_cast(dofsIt->first), dofsPosIndex[dofsIt->first]); else nRemainDofs++; } + + // === Finally, remove the macro elements from memory. === + + for (std::set::iterator macroIt = macros.begin(); + macroIt != macros.end(); ++macroIt) + delete *macroIt; + nVertices = nRemainDofs; } @@ -784,9 +812,9 @@ namespace AMDiS { int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2); MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n", - ichild, c_outside, c_lambda[0],c_lambda[1],c_lambda[2],c_lambda[3]); + ichild, c_outside, c_lambda[0], c_lambda[1], c_lambda[2], c_lambda[3]); MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n", - 1-ichild, c_outside2, c_lambda2[0],c_lambda2[1],c_lambda2[2], + 1 - ichild, c_outside2, c_lambda2[0], c_lambda2[1], c_lambda2[2], c_lambda2[3]); if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) { @@ -830,23 +858,24 @@ namespace AMDiS { return inside; } - bool Mesh::getDofIndexCoords(const DegreeOfFreedom* dof, + + bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, const FiniteElemSpace* feSpace, WorldVector& coords) { DimVec* baryCoords; bool found = false; TraverseStack stack; - Vector dofVec(feSpace->getBasisFcts()->getNumber()); + std::vector dofVec(feSpace->getBasisFcts()->getNumber()); ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while (elInfo) { - feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), - feSpace->getAdmin(), - &dofVec); + feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), + feSpace->getAdmin(), + dofVec); for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) { - if (dofVec[i] == *dof) { + if (dofVec[i] == dof) { baryCoords = feSpace->getBasisFcts()->getCoords(i); elInfo->coordToWorld(*baryCoords, coords); found = true; @@ -863,39 +892,30 @@ namespace AMDiS { return found; } - bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, - const FiniteElemSpace* feSpace, - WorldVector& coords) + + void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace, + DOFVector >& coords) { - DimVec* baryCoords; - bool found = false; - TraverseStack stack; - Vector dofVec(feSpace->getBasisFcts()->getNumber()); + const BasisFunction* basFcts = feSpace->getBasisFcts(); + int nBasFcts = basFcts->getNumber(); + std::vector dofVec(nBasFcts); - ElInfo *elInfo = stack.traverseFirst(this, -1, - Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); + TraverseStack stack; + ElInfo *elInfo = + stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while (elInfo) { - feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), - feSpace->getAdmin(), - &dofVec); - for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) { - if (dofVec[i] == dof) { - baryCoords = feSpace->getBasisFcts()->getCoords(i); - elInfo->coordToWorld(*baryCoords, coords); - found = true; - break; - } + basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec); + for (int i = 0; i < nBasFcts; i++) { + DimVec *baryCoords = basFcts->getCoords(i); + elInfo->coordToWorld(*baryCoords, coords[dofVec[i]]); } - - if (found) - break; - + elInfo = stack.traverseNext(elInfo); } - return found; } + void Mesh::setDiameter(const WorldVector& w) { diam = w; diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index 80e050f4..fa5c6953 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -414,7 +414,7 @@ namespace AMDiS { * that there are no global or local refinements, i.e., all macro elements have * no children. */ - void removeMacroElements(std::vector& macros, + void removeMacroElements(std::set& macros, const FiniteElemSpace* feSpace); /// Frees the array of DOF pointers (see \ref createDOFPtrs) @@ -485,12 +485,31 @@ namespace AMDiS { */ bool getDofIndexCoords(const DegreeOfFreedom* dof, const FiniteElemSpace* feSpace, - WorldVector& coords); + WorldVector& coords) + { + return getDofIndexCoords(*dof, feSpace, coords); + } + + /** \brief + * This function is equal to \ref getDofIndexCoords as defined above, but takes + * a DOF index instead of a DOF pointer. + */ bool getDofIndexCoords(DegreeOfFreedom dof, const FiniteElemSpace* feSpace, WorldVector& coords); + /** \brief + * Traverse the whole mesh and stores to each DOF of the given finite element space + * the coordinates in a given DOFVector. Works in the same way as the function + * \ref getDofIndexCoords defined above. + * + * @param[in] feSpace The fe soace to be used for the search. + * @param[out] coords DOF vector that stores the coordinates to each dof. + */ + void getDofIndexCoords(const FiniteElemSpace* feSpace, + DOFVector >& coords); + /// Returns FILL_ANY_?D inline static const Flag& getFillAnyFlag(int dim) { diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc index 59489b48..3fedf06d 100644 --- a/AMDiS/src/MeshStructure.cc +++ b/AMDiS/src/MeshStructure.cc @@ -6,6 +6,7 @@ #include "Traverse.h" #include "ElInfo.h" #include "RefinementManager.h" +#include "InteriorBoundary.h" namespace AMDiS { @@ -59,7 +60,46 @@ namespace AMDiS { clear(); - addAlongSide(el, ithSide, elType, reverseOrder); + int s1 = el->getSideOfChild(0, ithSide, elType); + int s2 = el->getSideOfChild(1, ithSide, elType); + + TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n"); + + if (!el->isLeaf()) { + if (s1 == -1) + addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder); + else if (s2 == -1) + addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder); + else + addAlongSide(el, ithSide, elType, reverseOrder); + } + + commit(); + } + + void MeshStructure::init(BoundaryObject &bound) + { + FUNCNAME("MeshStructure::init()"); + + Element *el = bound.el; + + clear(); + + int s1 = el->getSideOfChild(0, bound.ithObj, bound.elType); + int s2 = el->getSideOfChild(1, bound.ithObj, bound.elType); + + TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n"); + + if (!el->isLeaf()) { + if (s1 == -1) + addAlongSide(el->getSecondChild(), s2, + el->getChildType(bound.elType), bound.reverseMode); + else if (s2 == -1) + addAlongSide(el->getFirstChild(), s1, + el->getChildType(bound.elType), bound.reverseMode); + else + addAlongSide(el, bound.ithObj, bound.elType, bound.reverseMode); + } commit(); } @@ -67,6 +107,14 @@ namespace AMDiS { void MeshStructure::addAlongSide(Element *el, int ithSide, int elType, bool reverseOrder) { + FUNCNAME("MeshStructure::addAlongSide()"); + + if (debugMode) { + MSG("addAlondSide(%d, %d, %d, %d)\n", + el->getIndex(), ithSide, elType, reverseOrder); + MSG("Element is leaf: %d\n", el->isLeaf()); + } + insertElement(el->isLeaf()); if (!el->isLeaf()) { @@ -75,14 +123,14 @@ namespace AMDiS { if (!reverseOrder) { if (s1 != -1) - addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder); + addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), false); if (s2 != -1) - addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder); + addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), false); } else { if (s2 != -1) - addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder); + addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), false); if (s1 != -1) - addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder); + addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), false); } } } @@ -90,9 +138,13 @@ namespace AMDiS { void MeshStructure::reset() { currentIndex = 0; - currentCode = code[0]; pos = 0; currentElement = 0; + + if (code.size() > 0) + currentCode = code[0]; + else + currentCode = 0; } bool MeshStructure::nextElement(MeshStructure *insert) @@ -238,4 +290,8 @@ namespace AMDiS { } while (!finished); } + bool MeshStructure::compare(MeshStructure &other) + { + return (other.getCode() == code); + } } diff --git a/AMDiS/src/MeshStructure.h b/AMDiS/src/MeshStructure.h index 7e79326d..165d2948 100644 --- a/AMDiS/src/MeshStructure.h +++ b/AMDiS/src/MeshStructure.h @@ -36,7 +36,8 @@ namespace AMDiS { currentCode(0), pos(0), currentElement(0), - nElements(0) + nElements(0), + debugMode(false) {} void clear(); @@ -46,6 +47,8 @@ namespace AMDiS { void init(Element *el, int ithSide, int elType, bool reverseOrder); + void init(BoundaryObject &bound); + void init(const std::vector& initCode, int n) { code = initCode; @@ -59,6 +62,12 @@ namespace AMDiS { */ void reset(); + /// Returns whether the code is empty or not. + inline bool empty() + { + return (nElements == 0); + } + inline void commit() { if (pos > 0) @@ -104,6 +113,11 @@ namespace AMDiS { { FUNCNAME("MeshStructure::print()"); + if (empty()) { + std::cout << "-" << std::endl; + return; + } + reset(); bool cont = true; while (cont) { @@ -133,6 +147,14 @@ namespace AMDiS { return currentElement; } + void setDebugMode(bool b) + { + debugMode = b; + } + + /// Returns true, if the given mesh structure code is equal to this one. + bool compare(MeshStructure &other); + protected: /// Insert a new element to the structure code. Is used by the init function. void insertElement(bool isLeaf); @@ -157,6 +179,8 @@ namespace AMDiS { int nElements; + bool debugMode; + static const int unsignedLongSize; }; diff --git a/AMDiS/src/ParMetisPartitioner.cc b/AMDiS/src/ParMetisPartitioner.cc index f06a4c82..57b98cb0 100644 --- a/AMDiS/src/ParMetisPartitioner.cc +++ b/AMDiS/src/ParMetisPartitioner.cc @@ -24,8 +24,7 @@ namespace AMDiS { int dow = Global::getGeo(WORLD); TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh, -1, - Mesh::CALL_EVERY_EL_PREORDER); + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { // get partition data PartitionElementData *partitionData = dynamic_cast @@ -33,9 +32,8 @@ namespace AMDiS { if (partitionData && partitionData->getPartitionStatus() == IN && - partitionData->getLevel() == 0) { - elementCounter++; - } + partitionData->getLevel() == 0) + elementCounter++; elInfo = stack.traverseNext(elInfo); } @@ -75,8 +73,7 @@ namespace AMDiS { elementCounter = 0; elInfo = stack.traverseFirst(mesh, -1, - Mesh::CALL_EVERY_EL_PREORDER | - Mesh::FILL_COORDS); + Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_COORDS); while (elInfo) { Element *element = elInfo->getElement(); int index = element->getIndex(); @@ -169,8 +166,7 @@ namespace AMDiS { void ParMetisPartitioner::deletePartitionData() { TraverseStack stack; - ElInfo *elInfo; - elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER); + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { Element *element = elInfo->getElement(); element->deleteElementData(PARTITION_ED); @@ -178,16 +174,17 @@ namespace AMDiS { } } - void ParMetisPartitioner::createPartitionData() { + void ParMetisPartitioner::createPartitionData() + { int mpiRank = mpiComm->Get_rank(); int mpiSize = mpiComm->Get_size(); + int nLeaves = mesh->getNumberOfLeaves(); + int elPerRank = nLeaves / mpiSize; - TraverseStack stack; - ElInfo *elInfo; + // === Create initial partitioning of the AMDiS mesh. === - // === create initial partitioning on AMDiS mesh === - int totalElementCounter = 0; - elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_LEAF_EL); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *element = elInfo->getElement(); @@ -198,13 +195,12 @@ namespace AMDiS { new PartitionElementData(element->getElementData()); element->setElementData(elData); - if (totalElementCounter % mpiSize == mpiRank) + if (element->getIndex() >= mpiRank * elPerRank && + element->getIndex() < (mpiRank + 1) * elPerRank) elData->setPartitionStatus(IN); else elData->setPartitionStatus(UNDEFINED); - totalElementCounter++; - elInfo = stack.traverseNext(elInfo); } } @@ -213,15 +209,15 @@ namespace AMDiS { PartitionMode mode, float itr) { - int mpiSize = mpiComm->Get_size(); + FUNCNAME("ParMetisPartitioner::partition()"); - TraverseStack stack; - ElInfo *elInfo; + int mpiSize = mpiComm->Get_size(); // === create parmetis mesh === if (parMetisMesh) delete parMetisMesh; - parMetisMesh = new ParMetisMesh(mesh_, mpiComm); + + parMetisMesh = new ParMetisMesh(mesh, mpiComm); int nElements = parMetisMesh->getNumElements(); @@ -231,7 +227,8 @@ namespace AMDiS { float maxWgt = 0.0; float *ptr_floatWgts = floatWgts; - elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER); + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { Element *element = elInfo->getElement(); @@ -242,6 +239,7 @@ namespace AMDiS { if (partitionData && partitionData->getPartitionStatus() == IN && partitionData->getLevel() == 0) { + int index = element->getIndex(); // get weight @@ -267,6 +265,7 @@ namespace AMDiS { int numflag = 0; // c numbering style! int ncon = elemWeights ? 1 : 0; // one weight at each vertex! int nparts = mpiSize; // number of partitions + float *tpwgts = elemWeights ? new float[mpiSize] : NULL; float ubvec = 1.05; int options[3] = {0, 0, 0}; // default options @@ -278,7 +277,7 @@ namespace AMDiS { for (int i = 0; i < mpiSize; i++) tpwgts[i] = 1.0 / nparts; - float scale = 10000 / maxWgt; + float scale = 10000.0 / maxWgt; // scale wgts for (int i = 0; i < nElements; i++) wgts[i] = static_cast(floatWgts[i] * scale); @@ -288,12 +287,14 @@ namespace AMDiS { switch(mode) { case INITIAL: + // TODO: Removed weight because it does not work correctly for very large + // macro meshes. ParMETIS_V3_PartKway(parMetisMesh->getElementDist(), parMetisGraph.getXAdj(), parMetisGraph.getAdjncy(), - wgts, NULL, - &wgtflag, + NULL, + NULL, &numflag, &ncon, &nparts, @@ -370,8 +371,8 @@ namespace AMDiS { partitionVec->clear(); // update ParMETIS mesh to new partitioning - if (!parMetisMesh) - parMetisMesh = new ParMetisMesh(mesh_, mpiComm); + if (!parMetisMesh) + parMetisMesh = new ParMetisMesh(mesh, mpiComm); int mpiRank = mpiComm->Get_rank(); int mpiSize = mpiComm->Get_size(); @@ -408,6 +409,8 @@ namespace AMDiS { void ParMetisPartitioner::distributePartitioning(int *part) { + FUNCNAME("ParMetisPartitioner::distributePartitioning()"); + int mpiSize = mpiComm->Get_size(); int mpiRank = mpiComm->Get_rank(); int nElements = parMetisMesh->getNumElements(); @@ -421,14 +424,12 @@ namespace AMDiS { // collect number of partition elements from all ranks for this rank int *nRankElements = new int[mpiSize]; - mpiComm->Alltoall(nPartitionElements, 1, MPI_INT, - nRankElements, 1, MPI_INT); + mpiComm->Alltoall(nPartitionElements, 1, MPI_INT, nRankElements, 1, MPI_INT); // sum up partition elements over all ranks int *sumPartitionElements = new int[mpiSize]; mpiComm->Allreduce(nPartitionElements, sumPartitionElements, mpiSize, MPI_INT, MPI_SUM); - // prepare distribution (fill partitionElements with AMDiS indices) int *bufferOffset = new int[mpiSize]; @@ -477,7 +478,7 @@ namespace AMDiS { } TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER); + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { Element *element = elInfo->getElement(); @@ -513,6 +514,8 @@ namespace AMDiS { void ParMetisPartitioner::descendPartitionData(Element *element) { + FUNCNAME("ParMetisPartitioner::descendPartitionData()"); + if (!element->isLeaf()) { Element *child0 = element->getChild(0); Element *child1 = element->getChild(1); @@ -541,7 +544,7 @@ namespace AMDiS { { int partition = -1; TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh_, -1, Mesh::CALL_EVERY_EL_PREORDER); + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { Element *element = elInfo->getElement(); PartitionElementData *partitionData = dynamic_cast diff --git a/AMDiS/src/ParMetisPartitioner.h b/AMDiS/src/ParMetisPartitioner.h index afd17740..5ca73d47 100644 --- a/AMDiS/src/ParMetisPartitioner.h +++ b/AMDiS/src/ParMetisPartitioner.h @@ -167,8 +167,8 @@ namespace AMDiS { class ParMetisPartitioner { public: - ParMetisPartitioner(Mesh *mesh, MPI::Intracomm *comm) - : mesh_(mesh), + ParMetisPartitioner(Mesh *m, MPI::Intracomm *comm) + : mesh(m), mpiComm(comm), parMetisMesh(NULL) {} @@ -200,7 +200,7 @@ namespace AMDiS { void descendPartitionData(Element *element); protected: - Mesh *mesh_; + Mesh *mesh; MPI::Intracomm *mpiComm; diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc index 31e21ec4..f0e4a0d9 100644 --- a/AMDiS/src/ParallelDomainBase.cc +++ b/AMDiS/src/ParallelDomainBase.cc @@ -1,4 +1,5 @@ #include +#include #include "ParallelDomainBase.h" #include "ParMetisPartitioner.h" @@ -83,6 +84,8 @@ namespace AMDiS { // in some way. testForMacroMesh(); + MSG("START PARTITIONING!\n"); + // create an initial partitioning of the mesh partitioner->createPartitionData(); // set the element weights, which are 1 at the very first begin @@ -93,16 +96,22 @@ namespace AMDiS { #if (DEBUG != 0) ElementIdxToDofs elMap; dbgCreateElementMap(elMap); + if (mpiRank == 0) { + int writePartMesh = 1; + GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh); - if (mpiRank == 0) - writePartitioningMesh("part.vtu"); + if (writePartMesh > 0) + writePartitioningMesh("part.vtu"); + else + MSG("Skip write part mesh!\n"); + } #endif // === Create new global and local DOF numbering. === createLocalGlobalNumbering(); - // === Create interior boundary information === + // === Create interior boundary information. === createInteriorBoundaryInfo(); @@ -110,10 +119,14 @@ namespace AMDiS { removeMacroElements(); + // === If in debug mode, make some tests. === + #if (DEBUG != 0) + MSG("AMDiS runs in debug mode, so make some test ...\n"); dbgTestElementMap(elMap); dbgTestInteriorBoundary(); dbgTestCommonDofs(true); + MSG("Debug mode tests finished!\n"); #endif // === Reset all DOFAdmins of the mesh. === @@ -147,13 +160,19 @@ namespace AMDiS { void ParallelDomainBase::updateDofAdmins() { - int nAdmins = mesh->getNumberOfDOFAdmin(); - for (int i = 0; i < nAdmins; i++) { + FUNCNAME("ParallelDomainBase::updateDofAdmins()"); + + for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) { DOFAdmin& admin = const_cast(mesh->getDOFAdmin(i)); + + // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise, + // it is not possible to define a first DOF hole. + if (static_cast(admin.getSize()) == mapLocalGlobalDOFs.size()) + admin.enlargeDOFLists(); for (int j = 0; j < admin.getSize(); j++) admin.setDOFFree(j, true); - for (int j = 0; j < static_cast(mapLocalGlobalDOFs.size()); j++) + for (unsigned int j = 0; j < mapLocalGlobalDOFs.size(); j++) admin.setDOFFree(j, false); admin.setUsedSize(mapLocalGlobalDOFs.size()); @@ -224,7 +243,7 @@ namespace AMDiS { // === Traverse all non zero entries of the row and produce vector cols === // === with the column indices of all row entries and vector values === - // === with the corresponding values. + // === with the corresponding values. === for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { @@ -735,27 +754,72 @@ namespace AMDiS { { FUNCNAME("ParallelDomainBase::checkMeshChange()"); - // === If mesh has not been changed, return. === + // === If mesh has not been changed on all ranks, return. === + + int recvAllValues = 0; + int sendValue = static_cast(mesh->getChangeIndex() != lastMeshChangeIndex); + mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM); - if (mesh->getChangeIndex() == lastMeshChangeIndex) + if (recvAllValues == 0) return; - MSG("MESH CHANGED!\n"); + // === At least one rank mesh has been changed, so the boundaries must be === + // === adapted to the new mesh structure. === - // === Create mesh structure codes for all ranks boundary elements. === + clock_t first = clock(); + + do { + // To check the interior boundaries, the ownership of the boundaries is not + // important. Therefore, we add all boundaries to one boundary container. + RankToBoundMap allBound = myIntBoundary.boundary; + for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); + it != otherIntBoundary.boundary.end(); ++it) + allBound[it->first] = it->second; - typedef std::vector MeshCodeVec; - std::map sendCodes; + // === Check the boundaries and adapt mesh if necessary. === - for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); - it != myIntBoundary.boundary.end(); ++it) { + bool meshChanged = checkAndAdaptBoundary(allBound); + + // === Check on all ranks if at least one rank's mesh has changed. === + + int sendValue = static_cast(!meshChanged); + recvAllValues = 0; + mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM); + } while (recvAllValues != 0); + +#if (DEBUG != 0) + writeMesh(feSpace); +#endif + + INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n", + TIME_USED(first, clock())); + + // === Because the mesh has been changed, update the DOF numbering and mappings. === + + updateLocalGlobalNumbering(); + + exit(0); + } + + + bool ParallelDomainBase::checkAndAdaptBoundary(RankToBoundMap &allBound) + { + FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()"); + + // === Create mesh structure codes for all ranks boundary elements. === + + std::map sendCodes; + + for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) { for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { - MeshStructure elCode; - elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, - boundIt->rankObj.elType, false); + elCode.init(boundIt->rankObj); +// elCode.init(boundIt->rankObj.el, +// boundIt->rankObj.ithObj, +// boundIt->rankObj.elType, +// boundIt->rankObj.reverseMode); sendCodes[it->first].push_back(elCode); } } @@ -763,85 +827,47 @@ namespace AMDiS { StdMpi stdMpi(mpiComm); stdMpi.prepareCommunication(true); stdMpi.send(sendCodes); - stdMpi.recv(otherIntBoundary.boundary); - stdMpi.startCommunication(); - - + stdMpi.recv(allBound); + stdMpi.startCommunication(MPI_UNSIGNED_LONG); + // === Compare received mesh structure codes. === - + bool meshFitTogether = true; - if (mpiRank == 0) { - DOFVector ddd(feSpace, "tmp"); - ddd.set(0.0); - VtkWriter::writeFile(&ddd, "testold0.vtu"); - } - - if (mpiRank == 1) { - DOFVector ddd(feSpace, "tmp"); - ddd.set(0.0); - VtkWriter::writeFile(&ddd, "testold1.vtu"); - } - - - for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); - it != otherIntBoundary.boundary.end(); ++it) { - + for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) { + MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first]; int i = 0; - + for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { - - std::cout << "[" << mpiRank << "] GET CODE: "; - recvCodes[i].print(); - - MeshStructure elCode; - elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, - boundIt->rankObj.elType, true); - - std::cout << "[" << mpiRank << "] CALC CODE: "; - elCode.print(); - + + MeshStructure elCode; + elCode.init(boundIt->rankObj); +// elCode.init(boundIt->rankObj.el, +// boundIt->rankObj.ithObj, +// boundIt->rankObj.elType, +// boundIt->rankObj.reverseMode); + if (elCode.getCode() != recvCodes[i].getCode()) { TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n"); - recvCodes[i].reset(); - fitElementToMeshCode(refineManager, recvCodes[i], boundIt->rankObj.el, - boundIt->rankObj.ithObj, boundIt->rankObj.elType); - - meshFitTogether = false; + + bool b = fitElementToMeshCode(refineManager, + recvCodes[i], + boundIt->rankObj.el, + boundIt->rankObj.ithObj, + boundIt->rankObj.elType); + if (b) + meshFitTogether = false; } - + i++; } } - MSG("MESH FIT ALGO FINISHED!\n"); - - if (mpiRank == 0) { - DOFVector ddd(feSpace, "tmp"); - ddd.set(0.0); - VtkWriter::writeFile(&ddd, "test0.vtu"); - } - - if (mpiRank == 1) { - DOFVector ddd(feSpace, "tmp"); - ddd.set(0.0); - VtkWriter::writeFile(&ddd, "test1.vtu"); - } - - - if (!meshFitTogether) { - MSG("MESH STRUCTURE CHANGED!\n"); - - std::cout << "MESH HAS BEEN CHANGED!" << std::endl; - - exit(0); - } - - updateLocalGlobalNumbering(); + return meshFitTogether; } - + void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data) { @@ -907,8 +933,7 @@ namespace AMDiS { elemWeights.clear(); TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh, -1, - Mesh::CALL_EVERY_EL_PREORDER); + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { Element *element = elInfo->getElement(); @@ -936,6 +961,8 @@ namespace AMDiS { void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo) { + FUNCNAME("ParallelDomainBase::partitionMesh()"); + if (initialPartitionMesh) { initialPartitionMesh = false; partitioner->fillCoarsePartitionVec(&oldPartitionVec); @@ -1003,9 +1030,8 @@ namespace AMDiS { bound.rankObj.elType = elInfo->getType(); bound.rankObj.subObj = subObj; bound.rankObj.ithObj = i; - // Do not set a pointer to the element, because if will be deleted from - // mesh after partitioning and the pointer would become invalid. - bound.neighObj.el = NULL; + + bound.neighObj.el = elInfo->getNeighbour(i); bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex(); bound.neighObj.elType = -1; bound.neighObj.subObj = subObj; @@ -1035,10 +1061,15 @@ namespace AMDiS { 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; + if (mpiRank > otherElementRank) { + AtomicBoundary& b = myIntBoundary.getNewAtomic(otherElementRank); + b = bound; + b.rankObj.setReverseMode(b.neighObj, feSpace); + } else { + AtomicBoundary& b = otherIntBoundary.getNewAtomic(otherElementRank); + b = bound; + b.neighObj.setReverseMode(b.rankObj, feSpace); + } } } } @@ -1047,7 +1078,6 @@ namespace AMDiS { elInfo = stack.traverseNext(elInfo); } - // === Once we have this information, we must care about the order of the atomic === // === bounds in the three boundary handling object. Eventually all the bound- === // === aries have to be in the same order on both ranks that share the bounday. === @@ -1206,11 +1236,12 @@ namespace AMDiS { } // periodicBoundary.boundary.size() > 0 } - void ParallelDomainBase::removeMacroElements() { - std::vector macrosToRemove; + FUNCNAME("ParallelDomainBase::removeMacroElements()"); + + std::set macrosToRemove; for (std::deque::iterator it = mesh->firstMacroElement(); it != mesh->endOfMacroElements(); ++it) { @@ -1218,7 +1249,7 @@ namespace AMDiS { dynamic_cast ((*it)->getElement()->getElementData(PARTITION_ED)); if (partitionData->getPartitionStatus() != IN) - macrosToRemove.push_back(*it); + macrosToRemove.insert(*it); } mesh->removeMacroElements(macrosToRemove, feSpace); @@ -1248,12 +1279,10 @@ namespace AMDiS { mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM); rstart -= nRankDofs; - // === Create for all dofs in rank new indices. The new index must start at === // === index 0, must be continues and have the same order as the indices === // === had before. === - // Do not change the indices now, but create a new indexing and store it here. DofIndexMap rankDofsNewLocalIndex; isRankDof.clear(); @@ -1267,6 +1296,7 @@ namespace AMDiS { i++; } + // === Create for all rank owned dofs a new global indexing. === // Stores for dofs in rank a new global index. @@ -1282,6 +1312,7 @@ namespace AMDiS { i++; } + // === Create information which dof indices must be send and which must === // === be received. === @@ -1330,16 +1361,16 @@ namespace AMDiS { dofIt != sendIt->second.end(); ++dofIt) sendDofs[sendIt->first].push_back(dofIt->first); - StdMpi > > stdMpi(mpiComm); stdMpi.prepareCommunication(false); stdMpi.send(sendNewDofs); for (std::map::iterator recvIt = recvNewDofs.begin(); recvIt != recvNewDofs.end(); ++recvIt, i++) stdMpi.recv(recvIt->first, recvIt->second * 2); - stdMpi.startCommunication(); + stdMpi.startCommunication(MPI_INT); std::map > > &dofMap = stdMpi.getRecvData(); + // === Change dof indices at boundary from other ranks. === // Within this small data structure we track which dof index was already changed. @@ -1388,6 +1419,7 @@ namespace AMDiS { } } + // === Create now the local to global index and local to dof index mappings. === createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex, @@ -1399,6 +1431,7 @@ namespace AMDiS { lastMeshChangeIndex = mesh->getChangeIndex(); } + void ParallelDomainBase::updateLocalGlobalNumbering() { FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()"); @@ -1470,7 +1503,8 @@ namespace AMDiS { if (oldVertexDof[*dofIt]) { recvDofs[it->first].push_back(*dofIt); - DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), *dofIt); + DofContainer::iterator eraseIt = + find(rankDofs.begin(), rankDofs.end(), *dofIt); if (eraseIt != rankDofs.end()) rankDofs.erase(eraseIt); } @@ -1490,11 +1524,12 @@ namespace AMDiS { boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, - dofs); - boundIt->rankObj.el->getNonVertexDofs(feSpace, - boundIt->rankObj.ithObj, - boundIt->rankObj.elType, - dofs); + dofs, + boundIt->rankObj.reverseMode); +// boundIt->rankObj.el->getNonVertexDofs(feSpace, +// boundIt->rankObj.ithObj, +// boundIt->rankObj.elType, +// dofs); for (int i = 0; i < static_cast(dofs.size()); i++) dofsToSend.push_back(dofs[i]); @@ -1510,33 +1545,23 @@ namespace AMDiS { boundIt != it->second.end(); ++boundIt) { DofContainer dofs; - boundIt->rankObj.el->getNonVertexDofs(feSpace, - boundIt->rankObj.ithObj, - boundIt->rankObj.elType, - dofs); +// boundIt->rankObj.el->getNonVertexDofs(feSpace, +// boundIt->rankObj.ithObj, +// boundIt->rankObj.elType, +// dofs); boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, - dofs); - - if (mesh->getDim() == 2) { - for (int i = static_cast(dofs.size()) - 1; i >= 0; i--) { - DofContainer::iterator eraseIt = - find(rankDofs.begin(), rankDofs.end(), dofs[i]); - if (eraseIt != rankDofs.end()) - rankDofs.erase(eraseIt); - - dofsToRecv.push_back(dofs[i]); - } - } else { - for (int i = 0; i < static_cast(dofs.size()); i++) { - DofContainer::iterator eraseIt = - find(rankDofs.begin(), rankDofs.end(), dofs[i]); - if (eraseIt != rankDofs.end()) - rankDofs.erase(eraseIt); - - dofsToRecv.push_back(dofs[i]); - } + dofs, + boundIt->rankObj.reverseMode); + + for (int i = 0; i < static_cast(dofs.size()); i++) { + DofContainer::iterator eraseIt = + find(rankDofs.begin(), rankDofs.end(), dofs[i]); + if (eraseIt != rankDofs.end()) + rankDofs.erase(eraseIt); + + dofsToRecv.push_back(dofs[i]); } } } @@ -2211,53 +2236,40 @@ namespace AMDiS { { FUNCNAME("ParallelDomainBase::dbgTestCommonDofs()"); + clock_t first = clock(); + + int testCommonDofs = 1; + GET_PARAMETER(0, "dbg->test common dofs", "%d", &testCommonDofs); + if (testCommonDofs == 0) { + MSG("Skip test common dofs!\n"); + return; + } + // 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. RankToCoords sendCoords; + // A rank receives all boundary dofs that are at its interior boundaries but are - // not owned by the rank. This map stores for each rank the the coordinates of dofs + // not owned by the rank. This map stores for each rank the coordinates of dofs // this rank expectes to receive from. RankToCoords recvCoords; - WorldVector coords; - + DOFVector > coords(feSpace, "dofCorrds"); + mesh->getDofIndexCoords(feSpace, coords); + for (RankToDofContainer::iterator it = sendDofs.begin(); - it != sendDofs.end(); ++it) { - for (std::vector::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) { - bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords); - - if (!b) { - std::cout << "ERROR (send): Cannot find DOF in mesh!" << std::endl; - std::cout << " mpi rank = " << mpiRank - << " send to rank = " << it->first - << " dof = " << *dofIt << " = " << **dofIt << std::endl; - ERROR_EXIT("Should not happen!\n"); - } - - sendCoords[it->first].push_back(coords); - } - } + it != sendDofs.end(); ++it) + for (DofContainer::iterator dofIt = it->second.begin(); + dofIt != it->second.end(); ++dofIt) + sendCoords[it->first].push_back(coords[**dofIt]); for (RankToDofContainer::iterator it = recvDofs.begin(); - it != recvDofs.end(); ++it) { - for (std::vector::iterator dofIt = it->second.begin(); - dofIt != it->second.end(); ++dofIt) { - bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords); - - if (!b) { - std::cout << "ERROR (recv): Cannot find DOF in mesh!" << std::endl; - std::cout << " mpi rank = " << mpiRank - << " send to rank = " << it->first - << " dof = " << *dofIt << " = " << **dofIt << std::endl; - ERROR_EXIT("Should not happen!\n"); - } - - recvCoords[it->first].push_back(coords); - } - } + it != recvDofs.end(); ++it) + for (DofContainer::iterator dofIt = it->second.begin(); + dofIt != it->second.end(); ++dofIt) + recvCoords[it->first].push_back(coords[**dofIt]); std::vector sendSize(mpiSize, 0); std::vector recvSize(mpiSize, 0); @@ -2301,59 +2313,43 @@ namespace AMDiS { } } - std::map sendCoordsBuffer, recvCoordsBuffer; - requestCounter = 0; - - int dimOfWorld = Global::getGeo(WORLD); - - for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) { - sendCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld]; - - for (int i = 0; i < static_cast(it->second.size()); i++) - for (int j = 0; j < dimOfWorld; j++) - sendCoordsBuffer[it->first][i * dimOfWorld + j] = (it->second)[i][j]; - - request[requestCounter++] = mpiComm.Isend(sendCoordsBuffer[it->first], - it->second.size() * dimOfWorld, - MPI_DOUBLE, it->first, 0); - } - - for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) { - recvCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld]; - - request[requestCounter++] = mpiComm.Irecv(recvCoordsBuffer[it->first], - it->second.size() * dimOfWorld, - MPI_DOUBLE, it->first, 0); - } - - MPI::Request::Waitall(requestCounter, request); + // === Now we know that the number of send and received DOFs fits together. === + // === So we can check if also the coordinates of the communicated DOFs are === + // === the same on both corresponding ranks. === - for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) - delete [] sendCoordsBuffer[it->first]; + typedef std::vector > CoordsVec; + StdMpi stdMpi(mpiComm); + stdMpi.prepareCommunication(true); + stdMpi.send(sendCoords); + stdMpi.recv(recvCoords); + stdMpi.startCommunication(MPI_DOUBLE); double eps = 1e-13; + int dimOfWorld = Global::getGeo(WORLD); - for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) { + // === Compare the received with the expected coordinates. === + + for (RankToCoords::iterator it = stdMpi.getRecvData().begin(); + it != stdMpi.getRecvData().end(); ++it) { for (int i = 0; i < static_cast(it->second.size()); i++) { for (int j = 0; j < dimOfWorld; j++) { - bool c = fabs((it->second)[i][j] - - recvCoordsBuffer[it->first][i * dimOfWorld + j]) < eps; + bool c = fabs((it->second)[i][j] - recvCoords[it->first][i][j]) < eps; + // === Print error message if the coordinates are not the same. === + if (printCoords && !c) { - std::cout << "[DBG] i = " << i << std::endl; - std::cout.precision(5); - std::cout << "[DBG] " - << "Rank " << mpiRank << " from rank " << it->first + std::cout << "[DBG] i = " << i << std::endl; + std::cout << "[DBG] Rank " << mpiRank << " from rank " << it->first << " expect coords ("; for (int k = 0; k < dimOfWorld; k++) { - std::cout << (it->second)[i][k]; + std::cout << recvCoords[it->first][i][k]; if (k + 1 < dimOfWorld) std::cout << " / "; } std::cout << ") received coords ("; for (int k = 0; k < dimOfWorld; k++) { - std::cout << recvCoordsBuffer[it->first][i * dimOfWorld + k]; + std::cout << (it->second)[i][k]; if (k + 1 < dimOfWorld) std::cout << " / "; } @@ -2363,9 +2359,9 @@ namespace AMDiS { TEST_EXIT(c)("Wrong DOFs in rank %d!\n", mpiRank); } } - - delete [] recvCoordsBuffer[it->first]; } + + INFO(info, 8)("Test common dofs needed %.5f seconds\n", TIME_USED(first, clock())); } @@ -2482,4 +2478,114 @@ namespace AMDiS { ElementFileWriter::writeFile(vec, feSpace, filename); } + void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace) + { + using boost::lexical_cast; + + if (MPI::COMM_WORLD.Get_rank() == rank) { + DOFVector tmp(feSpace, "tmp"); + colorDofVectorByLocalElementDofs(tmp, feSpace->getMesh(), elIdx); + VtkWriter::writeFile(tmp, "tmp" + lexical_cast(elIdx) + ".vtu"); + } + } + + void writeMesh(FiniteElemSpace *feSpace, int rank) + { + using boost::lexical_cast; + + int myRank = MPI::COMM_WORLD.Get_rank(); + if (rank == -1 || myRank == rank) { + DOFVector tmp(feSpace, "tmp"); + VtkWriter::writeFile(tmp, "mesh" + lexical_cast(myRank) + ".vtu"); + } + } + + void colorDofVectorByLocalElementDofs(DOFVector& vec, Element *el) + { + // === Get local indices of the given element. === + + const BasisFunction *basisFcts = vec.getFESpace()->getBasisFcts(); + int nBasisFcts = basisFcts->getNumber(); + std::vector localDofs(nBasisFcts); + basisFcts->getLocalIndices(el, vec.getFESpace()->getAdmin(), localDofs); + + // === Set the values of the dof vector. === + + vec.set(0.0); + for (int i = 0; i < nBasisFcts; i++) + vec[localDofs[i]] = static_cast(i); + } + + bool colorDofVectorByLocalElementDofs(DOFVector& vec, Mesh *mesh, int elIndex) + { + FUNCNAME("colorDofVectorByLocalElementDofs()"); + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + if (elInfo->getElement()->getIndex() == elIndex) { + colorDofVectorByLocalElementDofs(vec, elInfo->getElement()); + return true; + } + elInfo = stack.traverseNext(elInfo); + } + + return false; + } + + Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof) + { + const BasisFunction* basFcts = feSpace->getBasisFcts(); + int nBasFcts = basFcts->getNumber(); + std::vector dofVec(nBasFcts); + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, + Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec); + for (int i = 0; i < nBasFcts; i++) + if (dofVec[i] == dof) + return elInfo->getElement(); + + elInfo = stack.traverseNext(elInfo); + } + + return NULL; + } + + Element* getLevel0ParentElement(Mesh *mesh, Element *el) + { + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + if (elInfo->getElement() == el) + return elInfo->getMacroElement()->getElement(); + + elInfo = stack.traverseNext(elInfo); + } + + return NULL; + } + + void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof) + { + Element *el = getDofIndexElement(feSpace, dof); + Element *parEl = getLevel0ParentElement(feSpace->getMesh(), el); + + std::cout << "DOF-INFO: dof = " << dof + << " elidx = " << el->getIndex() + << " pelidx = " << parEl->getIndex() << std::endl; + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, + Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + if (elInfo->getElement()->getIndex() == parEl->getIndex()) + std::cout << "EL INFO TO " << parEl->getIndex() << ": " + << elInfo->getType() << std::endl; + + elInfo = stack.traverseNext(elInfo); + } + } } diff --git a/AMDiS/src/ParallelDomainBase.h b/AMDiS/src/ParallelDomainBase.h index 28533eca..9a45327b 100644 --- a/AMDiS/src/ParallelDomainBase.h +++ b/AMDiS/src/ParallelDomainBase.h @@ -86,6 +86,8 @@ namespace AMDiS { typedef std::map DofIndexMap; typedef std::map > PeriodicDofMap; + + typedef std::vector MeshCodeVec; public: ParallelDomainBase(ProblemIterationInterface *iterationIF, @@ -257,8 +259,29 @@ namespace AMDiS { void setDofVector(Vec& petscVec, DOFVector* vec, int disMult = 1, int dispAdd = 0); + /** \brief + * This function checks if the mesh has changed on at least on rank. In this case, + * the interior boundaries are adapted on all ranks such that they fit together on + * all ranks. Furthermore the function \ref updateLocalGlobalNumbering() is called + * to update the dof numberings and mappings on all rank due to the new mesh + * structure. + */ void checkMeshChange(); + /** \brief + * Checks for all given interior boundaries if the elements fit together on both + * sides of the boundaries. If this is not the case, the mesh is adapted. Because + * refinement of a certain element may forces the refinement of other elements, + * it is not guaranteed that all rank's meshes fit together after this function + * terminates. Hence, it must be called until a stable mesh refinement is reached. + * If the mesh has been changed by this function, it returns true. Otherwise, it + * returns false, i.e., the given interior boundaries fit together on both sides. + * + * \param[in] allBound Defines a map from rank to interior boundaries which + * should be checked. + */ + bool checkAndAdaptBoundary(RankToBoundMap &allBound); + void dbgCreateElementMap(ElementIdxToDofs &elMap); void dbgTestElementMap(ElementIdxToDofs &elMap); @@ -268,7 +291,7 @@ namespace AMDiS { /** \brief * This function is used for debugging only. It traverses all interior boundaries * and compares the dof indices on them with the dof indices of the boundarys - * neighbours. The function fails, when dof indices on an interior boundary does + * neighbours. The function fails, when dof indices on an interior boundary do * not fit together. * * \param printCoords If true, the coords of all common dofs are printed to @@ -575,6 +598,19 @@ namespace AMDiS { long lastMeshChangeIndex; }; + void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace); + + void writeMesh(FiniteElemSpace *feSpace, int rank = -1); + + void colorDofVectorByLocalElementDofs(DOFVector& vec, Element *el); + + bool colorDofVectorByLocalElementDofs(DOFVector& vec, Mesh *mesh, int elIndex); + + Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof); + + Element* getLevel0ParentElement(Mesh *mesh, Element *el); + + void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof); } #endif // AMDIS_PARALLELDOMAINBASE_H diff --git a/AMDiS/src/ParallelDomainVec.cc b/AMDiS/src/ParallelDomainVec.cc index 47a89498..54fed899 100644 --- a/AMDiS/src/ParallelDomainVec.cc +++ b/AMDiS/src/ParallelDomainVec.cc @@ -54,12 +54,8 @@ namespace AMDiS { { FUNCNAME("ParallelDomainVec::initParallelization()"); - std::cout << "TEST RHS 1 " << probVec->getRHS()->getDOFVector(0) << ": " << probVec->getRHS()->getDOFVector(0)->getCoarsenOperation() << std::endl; - ParallelDomainBase::initParallelization(adaptInfo); - std::cout << "TEST RHS 2 " << probVec->getRHS()->getDOFVector(0) << ": " << probVec->getRHS()->getDOFVector(0)->getCoarsenOperation() << std::endl; - for (int i = 0; i < nComponents; i++) { for (int j = 0; j < nComponents; j++) if (probVec->getSystemMatrix(i, j)) diff --git a/AMDiS/src/PngWriter.cc b/AMDiS/src/PngWriter.cc index 63659bed..3e6f83f8 100644 --- a/AMDiS/src/PngWriter.cc +++ b/AMDiS/src/PngWriter.cc @@ -53,15 +53,15 @@ namespace AMDiS { } const BasisFunction* basisFcts = dataCollector->getFeSpace()->getBasisFcts(); - Vector localDofs(3); + std::vector localDofs(3); DOFVector* dofvalues = dataCollector->getValues(); elInfo = stack.traverseFirst(dataCollector->getMesh(), -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while (elInfo) { - basisFcts->getLocalIndicesVec(elInfo->getElement(), - dataCollector->getFeSpace()->getAdmin(), - &localDofs); + basisFcts->getLocalIndices(elInfo->getElement(), + dataCollector->getFeSpace()->getAdmin(), + localDofs); for (int i = 0; i < 3; i++) { if (imageType == 0) { diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index 2c5140a1..8a65c066 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/AMDiS/src/ProblemVec.cc @@ -534,7 +534,7 @@ namespace AMDiS { adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i); adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i); } else { - WARNING("no estimator for component %d\n" , i); + WARNING("No estimator for component %d\n" , i); } } } @@ -562,7 +562,7 @@ namespace AMDiS { if (marker[i]) markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]); else - WARNING("no marker for component %d\n", i); + WARNING("No marker for component %d\n", i); return markFlag; } @@ -1381,13 +1381,13 @@ namespace AMDiS { } // Compute estimate for every mesh element - Vector locInd(componentSpaces[i]->getBasisFcts()->getNumber()); + std::vector locInd(componentSpaces[i]->getBasisFcts()->getNumber()); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, Mesh::CALL_LEAF_EL); while (elInfo) { - componentSpaces[i]->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), - componentSpaces[i]->getAdmin(), - &locInd); + componentSpaces[i]->getBasisFcts()->getLocalIndices(elInfo->getElement(), + componentSpaces[i]->getAdmin(), + locInd); double estimate = 0.0; for (int j = 0; j < componentSpaces[i]->getBasisFcts()->getNumber(); j++) estimate += (*tmp)[locInd[j]]; diff --git a/AMDiS/src/RefinementManager.cc b/AMDiS/src/RefinementManager.cc index 757c0dd3..fb383c17 100644 --- a/AMDiS/src/RefinementManager.cc +++ b/AMDiS/src/RefinementManager.cc @@ -42,14 +42,13 @@ namespace AMDiS { mesh = aMesh; int nElements = mesh->getNumberOfLeaves(); - ElInfo *elInfo; newCoords = false; stack = new TraverseStack; doMoreRecursiveRefine = true; while (doMoreRecursiveRefine) { doMoreRecursiveRefine = false; - elInfo = + ElInfo *elInfo = stack->traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND); while (elInfo) { @@ -77,4 +76,42 @@ namespace AMDiS { } } + + void RefinementManager::refineElement(Mesh *aMesh, Element *el) + { + FUNCNAME("RefineManager::refineElement()"); + + mesh = aMesh; + int nElements = mesh->getNumberOfLeaves(); + newCoords = false; + doMoreRecursiveRefine = false; + stack = new TraverseStack; + + ElInfo *elInfo = + stack->traverseFirst(mesh, -1, + Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND); + while (elInfo) { + if (elInfo->getElement() == el) { + refineFunction(elInfo); + + TEST_EXIT(!doMoreRecursiveRefine)("Okay, this will not work in that way!\n"); + + break; + } + elInfo = stack->traverseNext(elInfo); + } + +#if 0 + if (newCoords) + setNewCoords(); // call of sub-class method +#endif + + delete stack; + + nElements -= mesh->getNumberOfLeaves(); + + if (nElements != 0) + aMesh->incChangeIndex(); + } + } diff --git a/AMDiS/src/RefinementManager.h b/AMDiS/src/RefinementManager.h index 4015eba2..5376545b 100644 --- a/AMDiS/src/RefinementManager.h +++ b/AMDiS/src/RefinementManager.h @@ -66,6 +66,8 @@ namespace AMDiS { */ virtual Flag refineMesh(Mesh *aMesh); + void refineElement(Mesh *amesh, Element *el); + /// All elements of the mesh will be refined. Flag globalRefine(Mesh *aMesh, int mark); diff --git a/AMDiS/src/RefinementManager3d.cc b/AMDiS/src/RefinementManager3d.cc index 487ce595..e8966900 100644 --- a/AMDiS/src/RefinementManager3d.cc +++ b/AMDiS/src/RefinementManager3d.cc @@ -24,13 +24,12 @@ namespace AMDiS { Tetrahedron *el = dynamic_cast(const_cast(ref_list->getElement(index))); Tetrahedron *child[2]; - int i, node; int el_type = ref_list->getType(index); child[0] = dynamic_cast(mesh->createNewElement(el)); child[1] = dynamic_cast(mesh->createNewElement(el)); - int mark = max(0, el->getMark()-1); + int mark = max(0, el->getMark() - 1); child[0]->setMark(mark); child[1]->setMark(mark); el->setMark(0); @@ -44,13 +43,14 @@ namespace AMDiS { el->setFirstChild(child[0]); el->setSecondChild(child[1]); - if (child[0]->getMark() > 0) doMoreRecursiveRefine = true; + if (child[0]->getMark() > 0) + doMoreRecursiveRefine = true; int n_vertices = mesh->getGeo(VERTEX); - child[0]->setDOF(n_vertices-1, dof[0]); - child[1]->setDOF(n_vertices-1, dof[0]); + child[0]->setDOF(n_vertices - 1, dof[0]); + child[1]->setDOF(n_vertices - 1, dof[0]); - for (i = 0; i < n_vertices-1; i++) { + for (int i = 0; i < n_vertices - 1; i++) { child[0]-> setDOF(i, const_cast(el->getDOF(Tetrahedron::childVertex[el_type][0][i]))); child[1]-> @@ -69,7 +69,7 @@ namespace AMDiS { /****************************************************************************/ if (mesh->getNumberOfDOFs(EDGE)) { - node = mesh->getNode(EDGE); + int node = mesh->getNode(EDGE); /****************************************************************************/ /* set pointers to those dof's that are handed on from the parant */ @@ -108,31 +108,31 @@ namespace AMDiS { } if (mesh->getNumberOfDOFs(FACE)) { - node = mesh->getNode(FACE); + int node = mesh->getNode(FACE); /****************************************************************************/ /* set pointers to those dof's that are handed on from the parant */ /****************************************************************************/ - child[0]->setDOF(node+3, const_cast( el->getDOF(node+1))); - child[1]->setDOF(node+3, const_cast( el->getDOF(node+0))); + child[0]->setDOF(node + 3, const_cast(el->getDOF(node + 1))); + child[1]->setDOF(node + 3, const_cast(el->getDOF(node + 0))); /****************************************************************************/ /* get new dof for the common face of child0 and child1 */ /****************************************************************************/ DegreeOfFreedom *newDOF = mesh->getDOF(FACE); - child[0]->setDOF(node, static_cast( newDOF)); - child[1]->setDOF(node, static_cast( newDOF)); + child[0]->setDOF(node, static_cast(newDOF)); + child[1]->setDOF(node, static_cast(newDOF)); } if (mesh->getNumberOfDOFs(CENTER)) { - node = mesh->getNode(CENTER); - child[0]->setDOF(node, const_cast( mesh->getDOF(CENTER))); - child[1]->setDOF(node, const_cast( mesh->getDOF(CENTER))); + int node = mesh->getNode(CENTER); + child[0]->setDOF(node, const_cast(mesh->getDOF(CENTER))); + child[1]->setDOF(node, const_cast(mesh->getDOF(CENTER))); } - if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(FACE)) + if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(FACE)) fillPatchConnectivity(ref_list, index); } @@ -197,7 +197,7 @@ namespace AMDiS { j = Tetrahedron::adjacentChild[adjc][i]; i_neigh = Tetrahedron::nChildFace[el_type][i][dir]; - j_neigh = Tetrahedron::nChildFace[n_type][j][opp_v-2]; + j_neigh = Tetrahedron::nChildFace[n_type][j][opp_v - 2]; /****************************************************************************/ /* adjust dof pointer in the edge in the common face of el and neigh and */ @@ -205,16 +205,14 @@ namespace AMDiS { /****************************************************************************/ if (mesh->getNumberOfDOFs(EDGE)) { - node0 = - mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir]; - node1 = - mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][opp_v-2]; + node0 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][i][dir]; + node1 = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[n_type][j][opp_v - 2]; TEST_EXIT_DBG(neigh->getChild(j)->getDOF(node1)) ("no dof on neighbour %d at node %d\n", neigh->getChild(j)->getIndex(), node1); - (const_cast( el->getChild(i)))-> + (const_cast(el->getChild(i)))-> setDOF(node0, const_cast( neigh->getChild(j)->getDOF(node1))); } if (mesh->getNumberOfDOFs(FACE)) { @@ -225,8 +223,8 @@ namespace AMDiS { ("no dof on neighbour %d at node %d\n", neigh->getChild(j)->getIndex(), node1); - (const_cast( el->getChild(i)))-> - setDOF(node0, const_cast( neigh->getChild(j)->getDOF(node1))); + (const_cast(el->getChild(i)))-> + setDOF(node0, const_cast(neigh->getChild(j)->getDOF(node1))); } } /* for (i = 0; i < 2; i++) */ @@ -285,24 +283,21 @@ namespace AMDiS { return 0; } + void RefinementManager3d::setNewCoords() { - ElInfo *el_info; - - Flag fillFlag = Mesh::CALL_EVERY_EL_PREORDER| - Mesh::FILL_BOUND| - Mesh::FILL_COORDS; - if (refList) delete refList; refList = new RCNeighbourList(mesh->getMaxEdgeNeigh()); - fillFlag |= Mesh::FILL_NEIGH; - el_info = stack->traverseFirst(mesh, -1, fillFlag); - while (el_info) { - newCoordsFct(el_info); - el_info = stack->traverseNext(el_info); + ElInfo *elInfo = stack->traverseFirst(mesh, -1, + Mesh::CALL_EVERY_EL_PREORDER | + Mesh::FILL_BOUND | Mesh::FILL_COORDS | + Mesh::FILL_NEIGH); + while (elInfo) { + newCoordsFct(elInfo); + elInfo = stack->traverseNext(elInfo); } } @@ -311,8 +306,8 @@ namespace AMDiS { RCNeighbourList* refineList, int n_neigh, bool bound) { - int i; - Tetrahedron *el = dynamic_cast(const_cast( refineList->getElement(0))); + Tetrahedron *el = + dynamic_cast(const_cast(refineList->getElement(0))); /* first element in the list */ DegreeOfFreedom *dof[3] = {NULL, NULL, NULL}; @@ -328,7 +323,7 @@ namespace AMDiS { dof[2] = mesh->getDOF(EDGE); } - for (i = 0; i < n_neigh; i++) + for (int i = 0; i < n_neigh; i++) bisectTetrahedron(refineList, i, dof, edge); /****************************************************************************/ @@ -354,14 +349,14 @@ namespace AMDiS { /* remove dof of the midpoint of the common refinement edge */ /****************************************************************************/ - el = dynamic_cast(const_cast( refineList->getElement(0))); - mesh->freeDOF(const_cast( el->getDOF(mesh->getNode(EDGE))), EDGE); + el = dynamic_cast(const_cast(refineList->getElement(0))); + mesh->freeDOF(const_cast(el->getDOF(mesh->getNode(EDGE))), EDGE); } - if (mesh->getNumberOfDOFs(EDGE) || - mesh->getNumberOfDOFs(FACE) || + if (mesh->getNumberOfDOFs(EDGE) || + mesh->getNumberOfDOFs(FACE) || mesh->getNumberOfDOFs(CENTER)) { - for (i = 0; i < n_neigh; i++) + for (int i = 0; i < n_neigh; i++) refineList->removeDOFParent(i); } } @@ -374,11 +369,10 @@ namespace AMDiS { if (bound) { mesh->incrementNumberOfEdges(n_neigh + 2); - mesh->incrementNumberOfFaces(2*n_neigh + 1); - newCoords = true; + mesh->incrementNumberOfFaces(2 * n_neigh + 1); } else { mesh->incrementNumberOfEdges(n_neigh + 1); - mesh->incrementNumberOfFaces(2*n_neigh); + mesh->incrementNumberOfFaces(2 * n_neigh); } return dof[0][0]; @@ -437,36 +431,6 @@ namespace AMDiS { } } - // LeafDataPeriodic *pd = - // dynamic_cast(el->getElementData(PERIODIC)); - - // if (pd) { - // std::list::iterator it; - // std::list::iterator end = - // pd->getInfoList().end(); - - // for (it = pd->getInfoList().begin(); it != end; ++it) { - // if (it->periodicMode != 0) { - // PeriodicBC *periodicCondition = mesh->getPeriodicBCMap()[it->type]; - // if (periodicCondition) { - // VertexVector *associated = mesh->getPeriodicAssociations()[it->type]; - // for (j = 0; j < vertices; j++) - // if (neigh->getDOF(j, 0) == (*associated)[edge[0][0]]) break; - // for (k = 0; k < vertices; k++) - // if (neigh->getDOF(k, 0) == (*associated)[edge[1][0]]) break; - - // TEST_EXIT_DBG(j < vertices && k < vertices) - // ("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", - // edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDOF(0,0), - // neigh->getDOF(1,0), neigh->getDOF(2,0), neigh->getDOF(3,0)); - - // } else { - // ERROR_EXIT("???\n"); - // } - // } - // } - // } - TEST_EXIT_DBG(j < mesh->getGeo(VERTEX) && k < mesh->getGeo(VERTEX)) ("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", @@ -584,9 +548,8 @@ namespace AMDiS { ElInfo* RefinementManager3d::refineFunction(ElInfo* el_info) { FUNCNAME("RefinementManager3d::refineFunction()"); - int n_neigh, bound = false; + int bound = false; DegreeOfFreedom *edge[2]; - RCNeighbourList *ref_list; if (el_info->getElement()->getMark() <= 0) return el_info; /* element may not be refined */ @@ -594,21 +557,25 @@ namespace AMDiS { /****************************************************************************/ /* get memory for a list of all elements at the refinement edge */ /****************************************************************************/ - ref_list = new RCNeighbourList(mesh->getMaxEdgeNeigh()); + RCNeighbourList *ref_list = new RCNeighbourList(mesh->getMaxEdgeNeigh()); ref_list->setElType(0, el_info->getType()); ref_list->setElement(0, el_info->getElement()); - n_neigh = 1; + int n_neigh = 1; + + if (el_info->getProjection(0) && + el_info->getProjection(0)->getType() == VOLUME_PROJECTION) + newCoords = true; /****************************************************************************/ /* give the refinement edge the right orientation */ /****************************************************************************/ if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1,0)) { - edge[0] = const_cast( el_info->getElement()->getDOF(0)); - edge[1] = const_cast( el_info->getElement()->getDOF(1)); + edge[0] = const_cast(el_info->getElement()->getDOF(0)); + edge[1] = const_cast(el_info->getElement()->getDOF(1)); } else { - edge[1] = const_cast( el_info->getElement()->getDOF(0)); - edge[0] = const_cast( el_info->getElement()->getDOF(1)); + edge[1] = const_cast(el_info->getElement()->getDOF(0)); + edge[0] = const_cast(el_info->getElement()->getDOF(1)); } /****************************************************************************/ diff --git a/AMDiS/src/StdMpi.h b/AMDiS/src/StdMpi.h index bb56836e..c44fb5d0 100644 --- a/AMDiS/src/StdMpi.h +++ b/AMDiS/src/StdMpi.h @@ -42,6 +42,11 @@ namespace AMDiS { return s; } + int intSizeOf(std::vector > &data) + { + return data.size() * Global::getGeo(WORLD); + } + void makeBuf(std::map &data, int* buf) { int i = 0; @@ -69,20 +74,22 @@ namespace AMDiS { } void makeBuf(std::vector &data, unsigned long int *buf) - { + { int pos = 0; for (unsigned int i = 0; i < data.size(); i++) { buf[pos++] = data[i].getCode().size(); buf[pos++] = data[i].getNumElements(); + for (unsigned int j = 0; j < data[i].getCode().size(); j++) buf[pos++] = data[i].getCode()[j]; } } - void makeFromBuf(std::vector &data, unsigned long int *buf, int bufSize) + void makeFromBuf(std::vector &data, + unsigned long int *buf, int bufSize) { int pos = 0; - + while (pos < bufSize) { int codeSize = buf[pos++]; int nElements = buf[pos++]; @@ -98,6 +105,27 @@ namespace AMDiS { } } + void makeBuf(std::vector > &data, double *buf) + { + int dimOfWorld = Global::getGeo(WORLD); + int pos = 0; + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + buf[pos++] = data[i][j]; + } + + void makeFromBuf(std::vector > &data, double *buf, int bufSize) + { + int dimOfWorld = Global::getGeo(WORLD); + TEST_EXIT(bufSize % Global::getGeo(WORLD) == 0)("This should not happen!\n"); + + int pos = 0; + data.resize(bufSize / Global::getGeo(WORLD)); + for (unsigned int i = 0; i < data.size(); i++) + for (int j = 0; j < dimOfWorld; j++) + data[i][j] = buf[pos++]; + } + template class StdMpi { @@ -190,7 +218,7 @@ namespace AMDiS { } template - void startCommunication() + void startCommunication(MPI_Datatype m) { FUNCNAME("StdMpi::startCommunication()"); @@ -210,7 +238,7 @@ namespace AMDiS { makeBuf(sendIt->second, buf); request[requestCounter++] = - mpiComm.Isend(buf, bufferSize, MPI_INT, sendIt->first, 0); + mpiComm.Isend(buf, bufferSize, m, sendIt->first, 0); sendBuffers.push_back(buf); } @@ -220,8 +248,8 @@ namespace AMDiS { T* buf = new T[recvIt->second]; request[requestCounter++] = - mpiComm.Irecv(buf, recvIt->second, MPI_INT, recvIt->first, 0); - + mpiComm.Irecv(buf, recvIt->second, m, recvIt->first, 0); + recvBuffers.push_back(buf); } @@ -238,6 +266,18 @@ namespace AMDiS { delete [] recvBuffers[i]; i++; } + + commPrepared = false; + } + + void clear() + { + sendData.clear(); + recvData.clear(); + sendDataSize.clear(); + recvDataSize.clear(); + + commPrepared = true; } protected: diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index fa8151f6..11c6e175 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -174,23 +174,32 @@ namespace AMDiS { void Tetrahedron::getVertexDofs(FiniteElemSpace* feSpace, int ith, int elType, - DofContainer& dofs, bool parentVertices) const + DofContainer& dofs, + bool reverseMode, + bool parentVertices) const { FUNCNAME("Tetrahedron::getVertexDofs()"); + bool p = false; +// if (MPI::COMM_WORLD.Get_rank() == 0) +// p = true; + + if (p) + MSG("GET-VD: idx=%d ith=%d type=%d\n", this->getIndex(), ith, elType); + if (parentVertices) { switch (ith) { case 0: - dofs.push_back(dof[1]); dofs.push_back(dof[2]), dofs.push_back(dof[3]); + dofs.push_back(dof[1]); dofs.push_back(dof[2]); dofs.push_back(dof[3]); break; case 1: - dofs.push_back(dof[0]); dofs.push_back(dof[2]), dofs.push_back(dof[3]); + dofs.push_back(dof[0]); dofs.push_back(dof[2]); dofs.push_back(dof[3]); break; case 2: - dofs.push_back(dof[0]); dofs.push_back(dof[1]), dofs.push_back(dof[3]); + dofs.push_back(dof[0]); dofs.push_back(dof[1]); dofs.push_back(dof[3]); break; case 3: - dofs.push_back(dof[0]); dofs.push_back(dof[1]), dofs.push_back(dof[2]); + dofs.push_back(dof[0]); dofs.push_back(dof[1]); dofs.push_back(dof[2]); break; default: ERROR_EXIT("Should never happen!\n"); @@ -205,16 +214,12 @@ namespace AMDiS { exit(0); } - Element *child1 = child[1]; - if (child1 && child1->getFirstChild()) { - int c1 = sideOfChild[elType][1][ith]; - int c2 = sideOfChild[(elType + 1) % 3][0][c1]; - int c3 = sideOfChild[(elType + 1) % 3][1][c1]; - - child1->getFirstChild()->getVertexDofs(feSpace, c2, (elType + 2) % 3, dofs); - dofs.push_back(child1->getFirstChild()->getDOF(3)); - child1->getSecondChild()->getVertexDofs(feSpace, c3, (elType + 2) % 3, dofs); - } + if (child[1]) + child[1]->getVertexDofs(feSpace, + sideOfChild[elType][1][ith], + (elType + 1) % 3, + dofs, + reverseMode); } break; case 1: @@ -224,16 +229,12 @@ namespace AMDiS { exit(0); } - Element *child0 = child[0]; - if (child0 && child0->getFirstChild()) { - int c1 = sideOfChild[elType][0][ith]; - int c2 = sideOfChild[(elType + 1) % 3][0][c1]; - int c3 = sideOfChild[(elType + 1) % 3][1][c1]; - - child0->getFirstChild()->getVertexDofs(feSpace, c2, (elType + 2) % 3, dofs); - dofs.push_back(child0->getFirstChild()->getDOF(3)); - child0->getSecondChild()->getVertexDofs(feSpace, c3, (elType + 2) % 3, dofs); - } + if (child[0]) + child[0]->getVertexDofs(feSpace, + sideOfChild[elType][0][ith], + (elType + 1) % 3, + dofs, + reverseMode); } break; @@ -244,14 +245,26 @@ namespace AMDiS { int c1 = sideOfChild[elType][0][ith]; int c2 = sideOfChild[elType][1][ith]; + if (p) + MSG("GO TO CHILD WITH SIDE %d %d\n", c1, c2); + if (c1 == -1 || c2 == -1) { std::cout << "ERROR 3 WITH elType = " << elType << std::endl; exit(0); } - child[0]->getVertexDofs(feSpace, c1, (elType + 1) % 3, dofs); - dofs.push_back(child[0]->getDOF(3)); - child[1]->getVertexDofs(feSpace, c2, (elType + 1) % 3, dofs); + if (reverseMode) { + child[1]->getVertexDofs(feSpace, c2, (elType + 1) % 3, dofs, false); + dofs.push_back(child[0]->getDOF(3)); + child[0]->getVertexDofs(feSpace, c1, (elType + 1) % 3, dofs, false); + } else { + child[0]->getVertexDofs(feSpace, c1, (elType + 1) % 3, dofs, false); + dofs.push_back(child[0]->getDOF(3)); + child[1]->getVertexDofs(feSpace, c2, (elType + 1) % 3, dofs, false); + } + } else { + if (p) + MSG("NO CHILD\n"); } } break; @@ -268,7 +281,7 @@ namespace AMDiS { if (child[0]) { int childFace0 = sideOfChild[elType][0][ith]; - int childFace1 = sideOfChild[elType][0][ith]; + int childFace1 = sideOfChild[elType][1][ith]; TEST_EXIT(childFace0 != -1 || childFace1 != -1) ("No new face for child elements!\n"); diff --git a/AMDiS/src/Tetrahedron.h b/AMDiS/src/Tetrahedron.h index 4de965e2..930291df 100644 --- a/AMDiS/src/Tetrahedron.h +++ b/AMDiS/src/Tetrahedron.h @@ -136,7 +136,9 @@ namespace AMDiS { } void getVertexDofs(FiniteElemSpace* feSpace, int ith, int elType, - DofContainer& dofs, bool parentVertices = 0) const; + DofContainer& dofs, + bool reverseMode, + bool parentVertices = false) const; void getNonVertexDofs(FiniteElemSpace* feSpace, int ith, int elType, DofContainer& dofs) const; diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc index 8b20f104..b9b476ae 100644 --- a/AMDiS/src/Triangle.cc +++ b/AMDiS/src/Triangle.cc @@ -77,7 +77,9 @@ namespace AMDiS { void Triangle::getVertexDofs(FiniteElemSpace* feSpace, int ithEdge, int elType, - DofContainer& dofs, bool parentVertices) const + DofContainer& dofs, + bool reverseMode, + bool parentVertices) const { FUNCNAME("Triangle::getVertexDofs()"); @@ -102,9 +104,9 @@ namespace AMDiS { { Element *child1 = child[1]; if (child1 && child1->getFirstChild()) { - child1->getFirstChild()->getVertexDofs(feSpace, 0, elType, dofs); + child1->getFirstChild()->getVertexDofs(feSpace, 0, elType, dofs, reverseMode); dofs.push_back(child1->getFirstChild()->getDOF(2)); - child1->getSecondChild()->getVertexDofs(feSpace, 1, elType, dofs); + child1->getSecondChild()->getVertexDofs(feSpace, 1, elType, dofs, reverseMode); } } break; @@ -112,17 +114,17 @@ namespace AMDiS { { Element *child0 = child[0]; if (child0 && child0->getFirstChild()) { - child0->getFirstChild()->getVertexDofs(feSpace, 0, elType, dofs); + child0->getFirstChild()->getVertexDofs(feSpace, 0, elType, dofs, reverseMode); dofs.push_back(child[0]->getFirstChild()->getDOF(2)); - child0->getSecondChild()->getVertexDofs(feSpace, 1, elType, dofs); + child0->getSecondChild()->getVertexDofs(feSpace, 1, elType, dofs, reverseMode); } } break; case 2: if (child[0]) { - child[0]->getVertexDofs(feSpace, 0, elType, dofs); + child[0]->getVertexDofs(feSpace, 0, elType, dofs, reverseMode); dofs.push_back(child[0]->getDOF(2)); - child[1]->getVertexDofs(feSpace, 1, elType, dofs); + child[1]->getVertexDofs(feSpace, 1, elType, dofs, reverseMode); } break; default: diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h index fa105389..423bfa1c 100644 --- a/AMDiS/src/Triangle.h +++ b/AMDiS/src/Triangle.h @@ -158,7 +158,9 @@ namespace AMDiS { void getVertexDofs(FiniteElemSpace* feSpace, int ithEdge, int elType, - DofContainer& dofs, bool parentVertices = 0) const; + DofContainer& dofs, + bool reverseMode, + bool parentVertices = false) const; void getNonVertexDofs(FiniteElemSpace* feSpace, int ithEdge, int elType, diff --git a/AMDiS/src/VtkWriter.h b/AMDiS/src/VtkWriter.h index 12243da0..6f5910d6 100644 --- a/AMDiS/src/VtkWriter.h +++ b/AMDiS/src/VtkWriter.h @@ -57,6 +57,12 @@ namespace AMDiS { /// May be used to simply write ParaView files. static void writeFile(DOFVector *values, std::string filename); + /// May be used to simply write ParaView files. + static void writeFile(DOFVector &values, std::string filename) + { + writeFile(&values, filename); + } + /// Set a compressing method for file output. void setCompression(FileCompression c) { -- GitLab