diff --git a/AMDiS/Makefile.in b/AMDiS/Makefile.in
index 4389e9635f6e4807a61d63e5c7cbf531dbea8ee1..55848b904d8a0272df0b4dee6362639d5ae825cc 100644
--- a/AMDiS/Makefile.in
+++ b/AMDiS/Makefile.in
@@ -142,6 +142,7 @@ PACKAGE_VERSION = @PACKAGE_VERSION@
PATH_SEPARATOR = @PATH_SEPARATOR@
PETSC_DIR = @PETSC_DIR@
RANLIB = @RANLIB@
+SED = @SED@
SET_MAKE = @SET_MAKE@
SHELL = @SHELL@
STRIP = @STRIP@
@@ -418,7 +419,7 @@ distdir: $(DISTFILES)
|| exit 1; \
fi; \
done
- -find $(distdir) -type d ! -perm -777 -exec chmod a+rwx {} \; -o \
+ -find $(distdir) -type d ! -perm -755 -exec chmod a+rwx,go+rx {} \; -o \
! -type d ! -perm -444 -links 1 -exec chmod a+r {} \; -o \
! -type d ! -perm -400 -exec chmod a+r {} \; -o \
! -type d ! -perm -444 -exec $(SHELL) $(install_sh) -c -m a+r {} {} \; \
diff --git a/AMDiS/aclocal.m4 b/AMDiS/aclocal.m4
index c6b83dae620f453eff109d72b180b3a18307599f..602417c17f403a97d4d9054388be743c913ef513 100644
--- a/AMDiS/aclocal.m4
+++ b/AMDiS/aclocal.m4
@@ -1578,10 +1578,27 @@ 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)); 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"
+ 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"
fi
# We used to test for /lib/ld.so.1 and disable shared libraries on
@@ -4288,6 +4305,9 @@ 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
@@ -4421,11 +4441,11 @@ striplib=$lt_striplib
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=$lt_[]_LT_AC_TAGVAR(predep_objects, $1)
+predep_objects=\`echo $lt_[]_LT_AC_TAGVAR(predep_objects, $1) | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=$lt_[]_LT_AC_TAGVAR(postdep_objects, $1)
+postdep_objects=\`echo $lt_[]_LT_AC_TAGVAR(postdep_objects, $1) | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -4437,7 +4457,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=$lt_[]_LT_AC_TAGVAR(compiler_lib_search_path, $1)
+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"\`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method=$lt_deplibs_check_method
@@ -4517,7 +4537,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=$lt_sys_lib_search_path_spec
+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"\`
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec
@@ -6319,7 +6339,7 @@ ifelse([AC_DISABLE_FAST_INSTALL])
AC_DEFUN([LT_AC_PROG_GCJ],
[AC_CHECK_TOOL(GCJ, gcj, no)
- test "x${GCJFLAGS+set}" = xset || GCJFLAGS="-g -O2"
+ test "x${GCJFLAGS+set}" = xset || GCJFLAGS="-g -O3"
AC_SUBST(GCJFLAGS)
])
@@ -6353,6 +6373,7 @@ 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
@@ -6385,6 +6406,7 @@ 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.am b/AMDiS/bin/Makefile.am
index 7f4612bd309aba321e2659f57df7a4a6329aab2b..b8e0ee11736d7657cb590fc4aab1d70620675e30 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -59,7 +59,7 @@ INCLUDES = $(AMDIS_INCLUDES) $(PARALLEL_INCLUDES) $(TEMPLATE_INCLUDES)
if AMDIS_DEBUG
libamdis_la_CXXFLAGS += -g -O0 -Wall -DDEBUG=1 $(OPENMP_FLAG) $(INCLUDES) #-pedantic
else
- libamdis_la_CXXFLAGS += -O2 -Wall -DDEBUG=0 -DNDEBUG $(OPENMP_FLAG) -ftemplate-depth-100 $(INCLUDES) #-pedantic
+ libamdis_la_CXXFLAGS += -O3 -Wall -DDEBUG=0 -DNDEBUG $(OPENMP_FLAG) -ftemplate-depth-100 $(INCLUDES) #-pedantic
endif
diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in
index 93c6c2648605a95c38147f32ee97966cf89d570a..76a8dd5fbc02e2de35acf182c66bf127296a4c66 100644
--- a/AMDiS/bin/Makefile.in
+++ b/AMDiS/bin/Makefile.in
@@ -64,7 +64,7 @@ host_triplet = @host@
@ENABLE_DUNE_TRUE@am__append_8 = -I$(DUNE_DIR)
@ENABLE_BOOST_TRUE@am__append_9 = -DHAVE_BOOST=1
@AMDIS_DEBUG_TRUE@am__append_10 = -g -O0 -Wall -DDEBUG=1 $(OPENMP_FLAG) $(INCLUDES) #-pedantic
-@AMDIS_DEBUG_FALSE@am__append_11 = -O2 -Wall -DDEBUG=0 -DNDEBUG $(OPENMP_FLAG) -ftemplate-depth-100 $(INCLUDES) #-pedantic
+@AMDIS_DEBUG_FALSE@am__append_11 = -O3 -Wall -DDEBUG=0 -DNDEBUG $(OPENMP_FLAG) -ftemplate-depth-100 $(INCLUDES) #-pedantic
subdir = bin
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
@@ -438,6 +438,7 @@ 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/compositeFEM/src/CompositeFEMOperator.cc b/AMDiS/compositeFEM/src/CompositeFEMOperator.cc
index ce4b67c5a36b381499ee81f62c189cf786fe8dec..114d321c9c400ff4ddaaa6a2283f086131277f81 100644
--- a/AMDiS/compositeFEM/src/CompositeFEMOperator.cc
+++ b/AMDiS/compositeFEM/src/CompositeFEMOperator.cc
@@ -123,10 +123,8 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
subElementAssembler->getNCol());
set_to_zero(subPolMat2);
- int myRank = omp_get_thread_num();
-
- if (!assembler[myRank]) {
- assembler[myRank] =
+ if (!assembler) {
+ assembler =
new StandardAssembler(this, NULL, NULL, NULL, NULL, rowFeSpace, colFeSpace);
}
@@ -138,7 +136,7 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
}
- assembler[myRank]->calculateElementMatrix(elInfo, elMat, 1.0);
+ assembler->calculateElementMatrix(elInfo, elMat, 1.0);
subElementAssembler->getSubPolytopeMatrix(subPolytope,
subElementAssembler,
elInfo,
@@ -258,12 +256,9 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
ElementVector subPolVec2(subElementAssembler->getNRow());
set_to_zero(subPolVec2);
- int myRank = omp_get_thread_num();
-
- if (!assembler[myRank]) {
- assembler[myRank] =
+ if (!assembler)
+ assembler =
new StandardAssembler(this, NULL, NULL, NULL, NULL, rowFeSpace, colFeSpace);
- }
if (elLS->getLevelSetDomain() ==
ElementLevelSet::LEVEL_SET_INTERIOR) {
@@ -272,7 +267,7 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
}
- assembler[myRank]->calculateElementVector(elInfo, elVec, 1.0);
+ assembler->calculateElementVector(elInfo, elVec, 1.0);
subElementAssembler->getSubPolytopeVector(subPolytope,
subElementAssembler,
elInfo,
diff --git a/AMDiS/configure b/AMDiS/configure
index e41f34bb6c5b2d5091e27b256a85000e395b2028..3f5fcb510ebf1bb1a5af4ac2536ea02252fc8f6b 100755
--- a/AMDiS/configure
+++ b/AMDiS/configure
@@ -462,7 +462,7 @@ ac_includes_default="\
# include <unistd.h>
#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_ZOLTAN_TRUE ENABLE_ZOLTAN_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_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_ZOLTAN_TRUE ENABLE_ZOLTAN_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_files=''
# Initialize some variables set by options.
@@ -2918,7 +2918,7 @@ if test "$ac_test_CFLAGS" = set; then
CFLAGS=$ac_save_CFLAGS
elif test $ac_cv_prog_cc_g = yes; then
if test "$GCC" = yes; then
- CFLAGS="-g -O2"
+ CFLAGS="-g -O3"
else
CFLAGS="-g"
fi
@@ -3599,7 +3599,7 @@ if test "$ac_test_CXXFLAGS" = set; then
CXXFLAGS=$ac_save_CXXFLAGS
elif test $ac_cv_prog_cxx_g = yes; then
if test "$GXX" = yes; then
- CXXFLAGS="-g -O2"
+ CXXFLAGS="-g -O3"
else
CXXFLAGS="-g"
fi
@@ -3977,6 +3977,7 @@ 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
@@ -4011,6 +4012,7 @@ done
fi
SED=$lt_cv_path_SED
+
echo "$as_me:$LINENO: result: $SED" >&5
echo "${ECHO_T}$SED" >&6
@@ -4451,7 +4453,7 @@ ia64-*-hpux*)
;;
*-*-irix6*)
# Find out which ABI we are using.
- echo '#line 4454 "configure"' > conftest.$ac_ext
+ echo '#line 4456 "configure"' > conftest.$ac_ext
if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
(eval $ac_compile) 2>&5
ac_status=$?
@@ -5586,7 +5588,7 @@ fi
# Provide some information about the compiler.
-echo "$as_me:5589:" \
+echo "$as_me:5591:" \
"checking for Fortran 77 compiler version" >&5
ac_compiler=`set X $ac_compile; echo $2`
{ (eval echo "$as_me:$LINENO: \"$ac_compiler --version </dev/null >&5\"") >&5
@@ -5711,7 +5713,7 @@ if test "$ac_test_FFLAGS" = set; then
FFLAGS=$ac_save_FFLAGS
elif test $ac_cv_prog_f77_g = yes; then
if test "x$ac_cv_f77_compiler_gnu" = xyes; then
- FFLAGS="-g -O2"
+ FFLAGS="-g -O3"
else
FFLAGS="-g"
fi
@@ -6649,11 +6651,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:6652: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:6654: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
- echo "$as_me:6656: \$? = $ac_status" >&5
+ echo "$as_me:6658: \$? = $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.
@@ -6917,11 +6919,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:6920: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:6922: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
- echo "$as_me:6924: \$? = $ac_status" >&5
+ echo "$as_me:6926: \$? = $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.
@@ -7021,11 +7023,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:7024: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:7026: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
- echo "$as_me:7028: \$? = $ac_status" >&5
+ echo "$as_me:7030: \$? = $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
@@ -8486,10 +8488,31 @@ 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 8495 "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)); 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"
+ 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"
fi
# We used to test for /lib/ld.so.1 and disable shared libraries on
@@ -9366,7 +9389,7 @@ else
lt_dlunknown=0; lt_dlno_uscore=1; lt_dlneed_uscore=2
lt_status=$lt_dlunknown
cat > conftest.$ac_ext <<EOF
-#line 9369 "configure"
+#line 9392 "configure"
#include "confdefs.h"
#if HAVE_DLFCN_H
@@ -9466,7 +9489,7 @@ else
lt_dlunknown=0; lt_dlno_uscore=1; lt_dlneed_uscore=2
lt_status=$lt_dlunknown
cat > conftest.$ac_ext <<EOF
-#line 9469 "configure"
+#line 9492 "configure"
#include "confdefs.h"
#if HAVE_DLFCN_H
@@ -9797,6 +9820,9 @@ CC=$lt_compiler
# Is the compiler the GNU C compiler?
with_gcc=$GCC
+gcc_dir=\`gcc -print-file-name=. | $SED 's,/\.$,,'\`
+gcc_ver=\`gcc -dumpversion\`
+
# An ERE matcher.
EGREP=$lt_EGREP
@@ -9930,11 +9956,11 @@ striplib=$lt_striplib
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=$lt_predep_objects
+predep_objects=\`echo $lt_predep_objects | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=$lt_postdep_objects
+postdep_objects=\`echo $lt_postdep_objects | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -9946,7 +9972,7 @@ postdeps=$lt_postdeps
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path=$lt_compiler_lib_search_path
+compiler_lib_search_path=\`echo $lt_compiler_lib_search_path | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method=$lt_deplibs_check_method
@@ -10026,7 +10052,7 @@ variables_saved_for_relink="$variables_saved_for_relink"
link_all_deplibs=$link_all_deplibs
# Compile-time system search path for libraries
-sys_lib_search_path_spec=$lt_sys_lib_search_path_spec
+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"\`
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec
@@ -11806,11 +11832,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:11809: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:11835: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
- echo "$as_me:11813: \$? = $ac_status" >&5
+ echo "$as_me:11839: \$? = $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.
@@ -11910,11 +11936,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:11913: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:11939: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
- echo "$as_me:11917: \$? = $ac_status" >&5
+ echo "$as_me:11943: \$? = $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
@@ -12442,10 +12468,31 @@ 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 12475 "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)); 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"
+ 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"
fi
# We used to test for /lib/ld.so.1 and disable shared libraries on
@@ -12829,6 +12876,9 @@ 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
@@ -12962,11 +13012,11 @@ striplib=$lt_striplib
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=$lt_predep_objects_CXX
+predep_objects=\`echo $lt_predep_objects_CXX | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=$lt_postdep_objects_CXX
+postdep_objects=\`echo $lt_postdep_objects_CXX | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -12978,7 +13028,7 @@ postdeps=$lt_postdeps_CXX
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path=$lt_compiler_lib_search_path_CXX
+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"\`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method=$lt_deplibs_check_method
@@ -13058,7 +13108,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=$lt_sys_lib_search_path_spec
+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"\`
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec
@@ -13480,11 +13530,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:13483: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:13533: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
- echo "$as_me:13487: \$? = $ac_status" >&5
+ echo "$as_me:13537: \$? = $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.
@@ -13584,11 +13634,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:13587: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:13637: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
- echo "$as_me:13591: \$? = $ac_status" >&5
+ echo "$as_me:13641: \$? = $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
@@ -15029,10 +15079,31 @@ 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 15086 "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)); 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"
+ 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"
fi
# We used to test for /lib/ld.so.1 and disable shared libraries on
@@ -15416,6 +15487,9 @@ 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
@@ -15549,11 +15623,11 @@ striplib=$lt_striplib
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=$lt_predep_objects_F77
+predep_objects=\`echo $lt_predep_objects_F77 | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=$lt_postdep_objects_F77
+postdep_objects=\`echo $lt_postdep_objects_F77 | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -15565,7 +15639,7 @@ postdeps=$lt_postdeps_F77
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path=$lt_compiler_lib_search_path_F77
+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"\`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method=$lt_deplibs_check_method
@@ -15645,7 +15719,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=$lt_sys_lib_search_path_spec
+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"\`
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec
@@ -15787,11 +15861,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:15790: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:15864: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
- echo "$as_me:15794: \$? = $ac_status" >&5
+ echo "$as_me:15868: \$? = $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.
@@ -16055,11 +16129,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:16058: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:16132: $lt_compile\"" >&5)
(eval "$lt_compile" 2>conftest.err)
ac_status=$?
cat conftest.err >&5
- echo "$as_me:16062: \$? = $ac_status" >&5
+ echo "$as_me:16136: \$? = $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.
@@ -16159,11 +16233,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:16162: $lt_compile\"" >&5)
+ (eval echo "\"\$as_me:16236: $lt_compile\"" >&5)
(eval "$lt_compile" 2>out/conftest.err)
ac_status=$?
cat out/conftest.err >&5
- echo "$as_me:16166: \$? = $ac_status" >&5
+ echo "$as_me:16240: \$? = $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
@@ -17624,10 +17698,31 @@ 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 17705 "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)); 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"
+ 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"
fi
# We used to test for /lib/ld.so.1 and disable shared libraries on
@@ -18011,6 +18106,9 @@ 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
@@ -18144,11 +18242,11 @@ striplib=$lt_striplib
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=$lt_predep_objects_GCJ
+predep_objects=\`echo $lt_predep_objects_GCJ | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=$lt_postdep_objects_GCJ
+postdep_objects=\`echo $lt_postdep_objects_GCJ | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -18160,7 +18258,7 @@ postdeps=$lt_postdeps_GCJ
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path=$lt_compiler_lib_search_path_GCJ
+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"\`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method=$lt_deplibs_check_method
@@ -18240,7 +18338,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=$lt_sys_lib_search_path_spec
+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"\`
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec
@@ -18492,6 +18590,9 @@ 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
@@ -18625,11 +18726,11 @@ striplib=$lt_striplib
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=$lt_predep_objects_RC
+predep_objects=\`echo $lt_predep_objects_RC | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=$lt_postdep_objects_RC
+postdep_objects=\`echo $lt_postdep_objects_RC | \$SED -e "s@\${gcc_dir}@\\\${gcc_dir}@g;s@\${gcc_ver}@\\\${gcc_ver}@g"\`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -18641,7 +18742,7 @@ postdeps=$lt_postdeps_RC
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path=$lt_compiler_lib_search_path_RC
+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"\`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method=$lt_deplibs_check_method
@@ -18721,7 +18822,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=$lt_sys_lib_search_path_spec
+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"\`
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec=$lt_sys_lib_dlsearch_path_spec
@@ -19654,6 +19755,7 @@ 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/libtool b/AMDiS/libtool
index 0215125986d6e6ec1968564bf7d9883394d45bd2..c5ffc30f6b0c85b7f91e0e4dbd12ec005c494a6c 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -30,10 +30,10 @@
# the same distribution terms that you use for the rest of that program.
# A sed program that does not truncate output.
-SED="/usr/bin/sed"
+SED="/bin/sed"
# Sed that helps us avoid accidentally triggering echo(1) options like -n.
-Xsed="/usr/bin/sed -e 1s/^X//"
+Xsed="/bin/sed -e 1s/^X//"
# The HP-UX ksh and POSIX shell print the target directory to stdout
# if CDPATH is set.
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host deimos103:
+# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -66,12 +66,12 @@ fast_install=yes
# The host system.
host_alias=
-host=x86_64-unknown-linux-gnu
+host=i686-pc-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=x86_64-unknown-linux-gnu
+build=i686-pc-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -82,22 +82,25 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
+LTCC="gcc"
# LTCC compiler flags.
-LTCFLAGS="-g -O2"
+LTCFLAGS="-g -O3"
# A language-specific compiler.
-CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
+CC="gcc"
# Is the compiler the GNU C compiler?
with_gcc=yes
+gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_ver=`gcc -dumpversion`
+
# An ERE matcher.
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
+LD="/usr/bin/ld"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -171,7 +174,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag=""
+link_static_flag="-static"
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -229,11 +232,11 @@ striplib="strip --strip-unneeded"
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=""
+predep_objects=`echo "" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=""
+postdep_objects=`echo "" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -245,7 +248,7 @@ postdeps=""
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path=""
+compiler_lib_search_path=`echo "" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -325,10 +328,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /fastfs/wir/local/lib/x86_64-suse-linux/4.1.2/ /fastfs/wir/local/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/ /usr/lib/gcc/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/ /lib/x86_64-suse-linux/4.1.2/ /lib/../lib64/ /usr/lib/x86_64-suse-linux/4.1.2/ /usr/lib/../lib64/ /fastfs/wir/local/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../ /lib/ /usr/lib/"
+sys_lib_search_path_spec=`echo " /u/witkowski/local/lib/i386-redhat-linux/4.1.2/ /u/witkowski/local/lib/ /u/witkowski/local/intel/mkl/10.0.1.014/lib/32/i386-redhat-linux/4.1.2/ /u/witkowski/local/intel/mkl/10.0.1.014/lib/32/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/mysql /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib /usr/lib/xulrunner-1.9.2 "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -6760,7 +6763,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host deimos103:
+# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6782,12 +6785,12 @@ fast_install=yes
# The host system.
host_alias=
-host=x86_64-unknown-linux-gnu
+host=i686-pc-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=x86_64-unknown-linux-gnu
+build=i686-pc-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -6798,22 +6801,25 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
+LTCC="gcc"
# LTCC compiler flags.
-LTCFLAGS="-g -O2"
+LTCFLAGS="-g -O3"
# A language-specific compiler.
-CC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicxx"
+CC="g++"
# Is the compiler the GNU C compiler?
with_gcc=yes
+gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_ver=`gcc -dumpversion`
+
# An ERE matcher.
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
+LD="/usr/bin/ld"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -6887,7 +6893,7 @@ dlopen_self=unknown
dlopen_self_static=unknown
# Compiler flag to prevent dynamic linking.
-link_static_flag=""
+link_static_flag="-static"
# Compiler flag to turn off builtin functions.
no_builtin_flag=" -fno-builtin"
@@ -6942,11 +6948,11 @@ striplib="strip --strip-unneeded"
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects="/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/crti.o /usr/lib64/gcc/x86_64-suse-linux/4.1.2/crtbeginS.o"
+predep_objects=`echo "/usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crti.o /usr/lib/gcc/i386-redhat-linux/4.1.2/crtbeginS.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects="/usr/lib64/gcc/x86_64-suse-linux/4.1.2/crtendS.o /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/crtn.o"
+postdep_objects=`echo "/usr/lib/gcc/i386-redhat-linux/4.1.2/crtendS.o /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crtn.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -6954,11 +6960,11 @@ predeps=""
# Dependencies to place after the objects being linked to create a
# shared library.
-postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -libverbs -lrt -lnuma -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
+postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path="-L/usr/lib64 -L/licsoft/libraries/openmpi/1.2.6/64bit/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2 -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/fastfs/wir/local/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../.."
+compiler_lib_search_path=`echo "-L/u/witkowski/local/lib -L/u/witkowski/local/intel/mkl/10.0.1.014/lib/32 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7038,10 +7044,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /fastfs/wir/local/lib/x86_64-suse-linux/4.1.2/ /fastfs/wir/local/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/ /usr/lib/gcc/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/../lib64/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../x86_64-suse-linux/4.1.2/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../lib64/ /lib/x86_64-suse-linux/4.1.2/ /lib/../lib64/ /usr/lib/x86_64-suse-linux/4.1.2/ /usr/lib/../lib64/ /fastfs/wir/local/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc/x86_64-suse-linux/4.1.2/../../../ /lib/ /usr/lib/"
+sys_lib_search_path_spec=`echo " /u/witkowski/local/lib/i386-redhat-linux/4.1.2/ /u/witkowski/local/lib/ /u/witkowski/local/intel/mkl/10.0.1.014/lib/32/i386-redhat-linux/4.1.2/ /u/witkowski/local/intel/mkl/10.0.1.014/lib/32/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/mysql /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib /usr/lib/xulrunner-1.9.2 "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -7065,7 +7071,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host deimos103:
+# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7087,12 +7093,12 @@ fast_install=yes
# The host system.
host_alias=
-host=x86_64-unknown-linux-gnu
+host=i686-pc-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=x86_64-unknown-linux-gnu
+build=i686-pc-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -7103,22 +7109,25 @@ AR="ar"
AR_FLAGS="cru"
# A C compiler.
-LTCC="/licsoft/libraries/openmpi/1.2.6/64bit/bin/mpicc"
+LTCC="gcc"
# LTCC compiler flags.
-LTCFLAGS="-g -O2"
+LTCFLAGS="-g -O3"
# A language-specific compiler.
CC="g77"
# Is the compiler the GNU C compiler?
-with_gcc=
+with_gcc=yes
+
+gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_ver=`gcc -dumpversion`
# An ERE matcher.
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
+LD="/usr/bin/ld"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -7250,11 +7259,11 @@ striplib="strip --strip-unneeded"
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=""
+predep_objects=`echo "" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=""
+postdep_objects=`echo "" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -7266,7 +7275,7 @@ postdeps=""
# The library search path used internally by the compiler when linking
# a shared library.
-compiler_lib_search_path=""
+compiler_lib_search_path=`echo "" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7346,10 +7355,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /fastfs/wir/local/lib/x86_64-suse-linux/3.3.5/ /fastfs/wir/local/lib/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/ /usr/lib/gcc/x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../../x86_64-suse-linux/lib/x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../../x86_64-suse-linux/lib/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../x86_64-suse-linux/3.3.5/ /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.5/../../../ /lib/x86_64-suse-linux/3.3.5/ /lib/ /usr/lib/x86_64-suse-linux/3.3.5/ /usr/lib/"
+sys_lib_search_path_spec=`echo " /u/witkowski/local/lib/i386-redhat-linux/3.4.6/ /u/witkowski/local/lib/ /u/witkowski/local/intel/mkl/10.0.1.014/lib/32/i386-redhat-linux/3.4.6/ /u/witkowski/local/intel/mkl/10.0.1.014/lib/32/ /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../ /lib/i386-redhat-linux/3.4.6/ /lib/ /usr/lib/i386-redhat-linux/3.4.6/ /usr/lib/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/x86_64-suse-linux/lib /usr/local/lib64 /usr/local/lib /opt/kde3/lib64 /opt/kde3/lib /opt/gnome/lib64 /opt/gnome/lib /lib64 /lib /usr/lib64 /usr/lib /opt/cluster/intel/cce/9.1.042/lib /opt/cluster/intel/fce/9.1.036/lib /opt/cluster/Pathscale3.0/lib/2.9.99 /opt/cluster/Pathscale3.0/lib/2.9.99/32 /work/licsoft/compilers/pgi/linux86-64/6.2/lib /work/licsoft/compilers/pgi/linux86-64/6.2/libso "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/mysql /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib /usr/lib/xulrunner-1.9.2 "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index ba8b2182e8097f619198c6d6e3e9fb31f1991cc5..ce5e77ce9638d97357d7689820b67f8ebc61fbd0 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -379,7 +379,7 @@ namespace AMDiS {
if (secondOrderAssembler)
secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
if (firstOrderAssemblerGrdPsi)
- firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
+ firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
if (firstOrderAssemblerGrdPhi)
firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
if (zeroOrderAssembler)
diff --git a/AMDiS/src/BasisFunction.cc b/AMDiS/src/BasisFunction.cc
index 4cc38a80ca206a28e4922767d41dfc1fe6566ff6..2c5c6ab85965a829310d8171e6a7202cad9ea7fc 100644
--- a/AMDiS/src/BasisFunction.cc
+++ b/AMDiS/src/BasisFunction.cc
@@ -36,24 +36,11 @@ namespace AMDiS {
nDOF = new DimVec<int>(dim, DEFAULT_VALUE, -1);
dow = Global::getGeo(WORLD);
-
- grdTmpVec1.resize(omp_get_overall_max_threads());
- grdTmpVec2.resize(omp_get_overall_max_threads());
-
- for (int i = 0; i < omp_get_overall_max_threads(); i++) {
- grdTmpVec1[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
- grdTmpVec2[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
- }
}
BasisFunction::~BasisFunction()
{
delete nDOF;
-
- for (int i = 0; i < static_cast<int>(grdTmpVec1.size()); i++) {
- delete grdTmpVec1[i];
- delete grdTmpVec2[i];
- }
}
/****************************************************************************/
@@ -99,27 +86,22 @@ namespace AMDiS {
const ElementVector& uh_loc,
WorldVector<double>* val) const
{
- TEST_EXIT_DBG(val)("return value is NULL\n");
-
- int myRank = omp_get_thread_num();
-
- DimVec<double> *grdTmp1 = grdTmpVec1[myRank];
- DimVec<double> *grdTmp2 = grdTmpVec2[myRank];
+ TEST_EXIT_DBG(val)("Return value is NULL\n");
- for (int j = 0; j < dim + 1; j++)
- (*grdTmp2)[j] = 0.0;
+ mtl::dense_vector<double> grdTmp1(dim + 1);
+ mtl::dense_vector<double> grdTmp2(dim + 1, 0.0);
for (int i = 0; i < nBasFcts; i++) {
- (*(*grdPhi)[i])(lambda, *grdTmp1);
+ (*(*grdPhi)[i])(lambda, grdTmp1);
for (int j = 0; j < dim + 1; j++)
- (*grdTmp2)[j] += uh_loc[i] * (*grdTmp1)[j];
+ grdTmp2[j] += uh_loc[i] * grdTmp1[j];
}
for (int i = 0; i < dow; i++) {
(*val)[i] = 0.0;
for (int j = 0; j < dim + 1; j++)
- (*val)[i] += grd_lambda[j][i] * (*grdTmp2)[j];
+ (*val)[i] += grd_lambda[j][i] * grdTmp2[j];
}
return ((*val));
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index 6274740bcacac21dbe66bd643c9292447a058474..55ad9ee80ed9ae8d200b3284f6eb24a606841c34 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -51,7 +51,8 @@ namespace AMDiS {
virtual ~GrdBasFctType() {}
- virtual void operator()(const DimVec<double>&, DimVec<double>&) const = 0;
+ virtual void operator()(const DimVec<double>&,
+ mtl::dense_vector<double>&) const = 0;
};
@@ -371,18 +372,6 @@ namespace AMDiS {
/// Vector of second derivatives
std::vector<D2BasFctType*> *d2Phi;
-
- /** \brief
- * Is used by function evalGrdUh. To make it thread safe, we need a
- * temporary DimVec for every thread.
- */
- std::vector<DimVec<double>* > grdTmpVec1;
-
- /** \brief
- * Is used by function evalGrdUh. To make it thread safe, we need a
- * temporary DimVec for every thread.
- */
- std::vector<DimVec<double>* > grdTmpVec2;
};
}
diff --git a/AMDiS/src/BoundaryManager.cc b/AMDiS/src/BoundaryManager.cc
index 72d9f5a66eff8729ae3b79c8d952384a549e7221..8248f2fcf289cb9ff4bf6fa968ae727cace7863d 100644
--- a/AMDiS/src/BoundaryManager.cc
+++ b/AMDiS/src/BoundaryManager.cc
@@ -17,7 +17,6 @@
#include "Traverse.h"
#include "BasisFunction.h"
#include "ElInfo.h"
-#include "OpenMP.h"
namespace AMDiS {
@@ -27,11 +26,8 @@ namespace AMDiS {
BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
{
- localBounds.resize(omp_get_overall_max_threads());
- dofIndices.resize(omp_get_overall_max_threads());
allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber();
- for (unsigned int i = 0; i < localBounds.size(); i++)
- localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
+ localBound = new BoundaryType[allocatedMemoryLocalBounds];
}
@@ -39,17 +35,13 @@ namespace AMDiS {
{
localBCs = bm.localBCs;
allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds;
- localBounds.resize(bm.localBounds.size());
- dofIndices.resize(bm.localBounds.size());
- for (unsigned int i = 0; i < localBounds.size(); i++)
- localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
+ localBound = new BoundaryType[allocatedMemoryLocalBounds];
}
BoundaryManager::~BoundaryManager()
{
- for (unsigned int i = 0; i < localBounds.size(); i++)
- delete [] localBounds[i];
+ delete [] localBound;
}
@@ -71,13 +63,11 @@ namespace AMDiS {
{
if (localBCs.size() > 0) {
const FiniteElemSpace *feSpace = vec->getFeSpace();
- std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
dofVec.resize(nBasFcts);
// get boundaries of all DOFs
- BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound);
// get dof indices
@@ -104,13 +94,11 @@ namespace AMDiS {
return;
const FiniteElemSpace *feSpace = mat->getRowFeSpace();
- std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
dofVec.resize(nBasFcts);
// get boundaries of all DOFs
- BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound);
// get dof indices
diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h
index 685c9a33ea9c004f6a88a95da2f2fe1284ad9a53..85fffc2aeaa304cba6f3ebf734965a13f74c3a6c 100644
--- a/AMDiS/src/BoundaryManager.h
+++ b/AMDiS/src/BoundaryManager.h
@@ -123,11 +123,11 @@ namespace AMDiS {
/// Map of managed local boundary conditions.
BoundaryIndexMap localBCs;
- /// Temporary thread-safe variable for functions fillBoundaryconditions.
- std::vector<BoundaryType*> localBounds;
+ /// Temporary variable for functions fillBoundaryconditions.
+ BoundaryType* localBound;
- /// Temporary thread-safe variable for functions fillBoundaryconditions.
- std::vector<std::vector<DegreeOfFreedom> > dofIndices;
+ /// Temporary variable for functions fillBoundaryconditions.
+ std::vector<DegreeOfFreedom> dofVec;
/** \brief
* Stores the number of byte that were allocated in the constructor for
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index bac9d45b6defb4f804d3885dbeaeefad26474d9e..11faf18cb92f56a399f7f06fa67b484ec0073502 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -265,8 +265,8 @@ namespace AMDiS {
std::vector<Operator*>::iterator it = operators.begin();
std::vector<double*>::iterator factorIt = operatorFactor.begin();
for (; it != operators.end(); ++it, ++factorIt)
- if ((*it)->getNeedDualTraverse() == false)
- (*it)->getElementMatrix(elInfo, elementMatrix, *factorIt ? **factorIt : 1.0);
+ if ((*it)->getNeedDualTraverse() == false)
+ (*it)->getElementMatrix(elInfo, elementMatrix, *factorIt ? **factorIt : 1.0);
if (factor != 1.0)
elementMatrix *= factor;
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 7e90223dd626be9a26e4b74603d4b7d60f8bf506..d79dbf32d09fa6bebd6363c10e7a1db6ff8ee1c3 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -264,14 +264,13 @@ namespace AMDiS {
const ElementVector& localVec = localVectors[myRank];
getLocalVector(elInfo->getElement(), const_cast<ElementVector&>(localVec));
- DimVec<double> &grd1 = *grdTmp[myRank];
+ mtl::dense_vector<double> grd1(dim + 1);
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (quadFast) {
for (int i = 0; i < nPoints; i++) {
- for (int j = 0; j < dim + 1; j++)
- grd1[j] = 0.0;
+ grd1 = 0.0;
for (int j = 0; j < nBasFcts; j++)
for (int k = 0; k < parts; k++)
@@ -286,16 +285,15 @@ namespace AMDiS {
} else {
const BasisFunction *basFcts = feSpace->getBasisFcts();
- DimVec<double>* grdPhi = grdPhis[myRank];
+ mtl::dense_vector<double> grdPhi(dim + 1);
for (int i = 0; i < nPoints; i++) {
- for (int j = 0; j < dim + 1; j++)
- grd1[j] = 0.0;
+ grd1 = 0.0;
for (int j = 0; j < nBasFcts; j++) {
- (*(basFcts->getGrdPhi(j)))(quad->getLambda(i), *grdPhi);
+ (*(basFcts->getGrdPhi(j)))(quad->getLambda(i), grdPhi);
for (int k = 0; k < parts; k++)
- grd1[k] += (*grdPhi)[k] * localVec[j];
+ grd1[k] += grdPhi[k] * localVec[j];
}
for (int l = 0; l < dow; l++) {
@@ -353,28 +351,27 @@ namespace AMDiS {
const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
- DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();
- DimVec<double>* grdPhi = grdPhis[myRank];
- DimVec<double> tmp(dim, DEFAULT_VALUE, 0.0);
+ mtl::dense_vector<double> grd1(dim + 1);
+ mtl::dense_vector<double> grdPhi(dim + 1);
+ mtl::dense_vector<double> tmp(dim + 1);
for (int i = 0; i < nPoints; i++) {
- for (int j = 0; j < dim + 1; j++)
- grd1[j] = 0.0;
+ grd1 = 0.0;
for (int j = 0; j < nBasFcts; j++) {
+ grdPhi = 0.0;
- grdPhi->set(0.0);
for (int k = 0; k < nBasFcts; k++) {
(*(basFcts->getGrdPhi(k)))(quad->getLambda(i), tmp);
tmp *= m[j][k];
- *grdPhi += tmp;
+ grdPhi += tmp;
}
for (int k = 0; k < parts; k++)
- grd1[k] += (*grdPhi)[k] * localVec[j];
+ grd1[k] += grdPhi[k] * localVec[j];
}
for (int l = 0; l < dow; l++) {
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index f0ec85ffe62b695714272a7b962f3784e99e4245..e93d0bc10f0e118254bf09e77379c05e2df0df07 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -269,12 +269,6 @@ namespace AMDiS {
/// Are used to store temporary local values of an element. Thread safe.
std::vector<mtl::dense_vector<T> > localVectors;
- /// Temporarly used in \ref getGrdAtQPs. Thread safe.
- std::vector<DimVec<double>*> grdTmp;
-
- /// Temporarly used in \ref getGrdAtQPs. Thread safe.
- std::vector<DimVec<double>*> grdPhis;
-
/// Temporarly used in \ref getD2AtQPs. Thread safe.
std::vector<DimMat<double>*> D2Phis;
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 95ad682cccb1f90513e18828bc9c9d0375899534..57d7eb46039f10e14b952597098e5a4db8730604 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -98,8 +98,6 @@ namespace AMDiS {
{
for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
delete [] localIndices[i];
- delete grdPhis[i];
- delete grdTmp[i];
delete D2Phis[i];
}
}
@@ -110,15 +108,11 @@ namespace AMDiS {
{
localIndices.resize(omp_get_overall_max_threads());
localVectors.resize(omp_get_overall_max_threads());
- grdPhis.resize(omp_get_overall_max_threads());
- grdTmp.resize(omp_get_overall_max_threads());
D2Phis.resize(omp_get_overall_max_threads());
for (int i = 0; i < omp_get_overall_max_threads(); i++) {
localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
localVectors[i].change_dim(this->nBasFcts);
- grdPhis[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
- grdTmp[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
D2Phis[i] = new DimMat<double>(dim, NO_INIT);
}
}
@@ -1010,7 +1004,6 @@ namespace AMDiS {
d[i] = (*this)[localIndices[i]];
}
-
template<typename T>
void DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
@@ -1018,35 +1011,35 @@ namespace AMDiS {
mtl::dense_vector<T>& vecAtQPs) const
{
FUNCNAME("DOFVector<T>::getVecAtQPs()");
-
- TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
-
- if (quad && quadFast) {
- TEST_EXIT_DBG(quad == quadFast->getQuadrature())
- ("quad != quadFast->quadrature\n");
- }
-
+
+ TEST_EXIT_DBG(quad || quadFast)
+ ("Neither quad nor quadFast defined!\n");
+ TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
+ ("quad != quadFast->quadrature\n");
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
- ("invalid basis functions");
-
- Element *el = elInfo->getElement();
- const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
- const BasisFunction *basFcts = feSpace->getBasisFcts();
- int nPoints = quadrature->getNumPoints();
- int nBasFcts = basFcts->getNumber();
-
- vecAtQPs.change_dim(nPoints);
+ ("Invalid basis functions!");
const mtl::dense_vector<T>& localVec = localVectors[omp_get_thread_num()];
- getLocalVector(el, const_cast<mtl::dense_vector<T>& >(localVec));
-
- for (int i = 0; i < nPoints; i++) {
- vecAtQPs[i] = 0.0;
- for (int j = 0; j < nBasFcts; j++)
- vecAtQPs[i] += localVec[j] *
- (quadFast ?
- (quadFast->getPhi(i, j)) :
- ((*(basFcts->getPhi(j)))(quad->getLambda(i))));
+ getLocalVector(elInfo->getElement(),
+ const_cast<mtl::dense_vector<T>& >(localVec));
+
+ if (quadFast) {
+ const mtl::dense2D<double>& phi = quadFast->getPhi();
+ vecAtQPs.change_dim(num_rows(phi));
+ vecAtQPs = phi * localVec;
+ } else {
+ const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
+ const BasisFunction *basFcts = feSpace->getBasisFcts();
+ int nPoints = quadrature->getNumPoints();
+ int nBasFcts = basFcts->getNumber();
+ vecAtQPs.change_dim(nPoints);
+
+ for (int i = 0; i < nPoints; i++) {
+ vecAtQPs[i] = 0.0;
+ for (int j = 0; j < nBasFcts; j++)
+ vecAtQPs[i] +=
+ localVec[j] * (*(basFcts->getPhi(j)))(quad->getLambda(i));
+ }
}
}
diff --git a/AMDiS/src/DirichletBC.cc b/AMDiS/src/DirichletBC.cc
index 010dd2d3fa4a720da2ef9b91eade4b67ce10e207..a739647d71268304ad26490c7ac3ebdb016efb67 100644
--- a/AMDiS/src/DirichletBC.cc
+++ b/AMDiS/src/DirichletBC.cc
@@ -27,9 +27,7 @@ namespace AMDiS {
f(fct),
dofVec(NULL),
applyBC(apply)
- {
- worldCoords.resize(omp_get_overall_max_threads());
- }
+ {}
DirichletBC::DirichletBC(BoundaryType type,
@@ -64,7 +62,6 @@ namespace AMDiS {
TEST_EXIT_DBG(vector->getFeSpace() == rowFeSpace)("invalid row fe space\n");
const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- int myRank = omp_get_thread_num();
for (int i = 0; i < nBasFcts; i++) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -72,8 +69,8 @@ namespace AMDiS {
#endif
if (localBound[i] == boundaryType) {
if (f) {
- elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords[myRank]);
- (*vector)[dofIndices[i]] = (*f)(worldCoords[myRank]);
+ elInfo->coordToWorld(*(basFcts->getCoords(i)), worldCoords);
+ (*vector)[dofIndices[i]] = (*f)(worldCoords);
}
if (dofVec)
(*vector)[dofIndices[i]] = (*dofVec)[dofIndices[i]];
diff --git a/AMDiS/src/DirichletBC.h b/AMDiS/src/DirichletBC.h
index 956cdad4bf983abe4beb3d7fb1e6582f3fb0efff..312bde445999fc784a0e414f3c2317d091c6e554 100644
--- a/AMDiS/src/DirichletBC.h
+++ b/AMDiS/src/DirichletBC.h
@@ -23,9 +23,10 @@
#ifndef AMDIS_DIRICHLETBC_H
#define AMDIS_DIRICHLETBC_H
+#include "AMDiS_fwd.h"
#include "BoundaryCondition.h"
#include "AbstractFunction.h"
-#include "AMDiS_fwd.h"
+#include "FixVec.h"
namespace AMDiS {
@@ -102,7 +103,7 @@ namespace AMDiS {
AbstractFunction<double, WorldVector<double> > *f;
///
- std::vector<WorldVector<double> > worldCoords;
+ WorldVector<double> worldCoords;
/// DOFVector containing the boundary values
DOFVectorBase<double> *dofVec;
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 2b543bf22a07e59b9f1be51bd3d0949defe525c2..40708e05a75896cdeda9ffcb5772baf2e976dbc3 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -37,9 +37,12 @@ namespace AMDiS {
FirstOrderType type)
: SubAssembler(op, assembler, quad, 1, optimized, type)
{
- VectorOfFixVecs<DimVec<double> >
- Lb(assembler->getRowFeSpace()->getMesh()->getDim(), 1, NO_INIT);
- tmpLb.resize(omp_get_overall_max_threads(), Lb);
+ FUNCNAME("FirstOrderAssmebler::FirstOrderAssembler()");
+
+ TEST_EXIT_DBG(dim > 0)("Should not happen!\n");
+
+ Lb.resize(1);
+ Lb[0].change_dim(dim + 1);
}
@@ -58,9 +61,8 @@ namespace AMDiS {
&standardSubAssemblersGrdPsi :
&standardSubAssemblersGrdPhi);
- int myRank = omp_get_thread_num();
std::vector<OperatorTerm*> opTerms =
- (type == GRD_PSI) ? op->firstOrderGrdPsi[myRank] : op->firstOrderGrdPhi[myRank];
+ (type == GRD_PSI) ? op->firstOrderGrdPsi : op->firstOrderGrdPhi;
// check if a assembler is needed at all
if (opTerms.size() == 0)
@@ -126,18 +128,17 @@ namespace AMDiS {
void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
- DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
-
+ mtl::dense_vector<double> grdPsi(dim + 1, 0.0);
int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
std::vector<double> phival(nCol);
- for (int iq = 0; iq < nPoints; iq++)
- Lb[iq].set(0.0);
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].change_dim(dim + 1);
+ Lb[iq] = 0.0;
+ }
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
@@ -148,7 +149,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
for (int j = 0; j < nCol; j++)
- mat[i][j] += quadrature->getWeight(iq) * phival[j] * (Lb[iq] * grdPsi);
+ mat[i][j] += quadrature->getWeight(iq) * phival[j] * dot(Lb[iq], grdPsi);
}
}
}
@@ -156,23 +157,24 @@ namespace AMDiS {
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
- DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
+ mtl::dense_vector<double> grdPsi(dim + 1, 0.0);
int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
- for (int iq = 0; iq < nPoints; iq++)
- Lb[iq].set(0.0);
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].change_dim(dim + 1);
+ Lb[iq] = 0.0;
+ }
+
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
- vec[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
+ vec[i] += quadrature->getWeight(iq) * dot(Lb[iq], grdPsi);
}
}
}
@@ -187,80 +189,68 @@ namespace AMDiS {
void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
- VectorOfFixVecs<DimVec<double> > *grdPsi;
-
if (firstCall) {
-#ifdef _OPENMP
-#pragma omp critical
-#endif
- {
- const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
- basFcts = colFeSpace->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
- }
+ const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
+ basFcts = colFeSpace->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
}
int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].change_dim(dim + 1);
+ Lb[iq] = 0.0;
+ }
- for (int iq = 0; iq < nPoints; iq++)
- Lb[iq].set(0.0);
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
+
+ const mtl::dense2D<double>& phi(phiFast->getPhi());
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
-
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
- const double *phi = phiFast->getPhi(iq);
- grdPsi = psiFast->getGradient(iq);
+ const vector<mtl::dense_vector<double> >& grdPsi = psiFast->getGradient(iq);
double factor = quadrature->getWeight(iq);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++)
- mat[i][j] += factor * (Lb[iq] * (*grdPsi)[i]) * phi[j];
+ mat[i][j] += factor * phi[iq][j] * dot(Lb[iq], grdPsi[i]);
}
}
+
+ // TODO: Do the same performance update es in Quad01::calculateElementMatrix
}
void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
- VectorOfFixVecs<DimVec<double> > *grdPsi;
-
if (firstCall) {
-#ifdef _OPENMP
-#pragma omp critical
-#endif
- {
- const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
- psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
- basFcts = colFeSpace->getBasisFcts();
- phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
- firstCall = false;
- }
+ const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
+ psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
+ basFcts = colFeSpace->getBasisFcts();
+ phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
+ firstCall = false;
}
int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].change_dim(dim + 1);
+ Lb[iq] = 0.0;
+ }
- for (int iq = 0; iq < nPoints; iq++)
- Lb[iq].set(0.0);
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
- grdPsi = psiFast->getGradient(iq);
+ const vector<mtl::dense_vector<double> >& grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++)
- vec[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
+ vec[i] += quadrature->getWeight(iq) * dot(Lb[iq], grdPsi[i]);
}
}
@@ -291,13 +281,11 @@ namespace AMDiS {
}
const int **nEntries = q10->getNumberEntries();
- int myRank = omp_get_thread_num();
// Do not need do resize Lb, because it's size is always at least one.
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
- Lb[0].set(0.0);
+ Lb[0] = 0.0;
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
Lb[0] *= elInfo->getDet();
@@ -313,30 +301,34 @@ namespace AMDiS {
}
}
+
Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
{
- VectorOfFixVecs<DimVec<double> >
- grdPhi(assembler->getRowFeSpace()->getMesh()->getDim(), nCol, NO_INIT);
- tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi);
+ FUNCNAME("Stand01::Stand01()");
+
+ TEST_EXIT_DBG(dim > 0)("Should not happen!\n");
+
+ grdPhi.resize(nCol);
+ for (int i = 0; i < nCol; i++)
+ grdPhi[i].change_dim(dim + 1);
psi = rowFeSpace->getBasisFcts();
phi = colFeSpace->getBasisFcts();
}
+
void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
- VectorOfFixVecs<DimVec<double> >& Lb = tmpLb[myRank];
- VectorOfFixVecs<DimVec<double> >& grdPhi = tmpGrdPhi[myRank];
Lb.resize(nPoints);
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].change_dim(dim + 1);
+ Lb[iq] = 0.0;
+ }
- for (int iq = 0; iq < nPoints; iq++)
- Lb[iq].set(0.0);
-
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
+ for (int i = 0; i < static_cast<int>(terms.size()); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
@@ -347,7 +339,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
for (int j = 0; j < nCol; j++)
- mat[i][j] += quadrature->getWeight(iq) * psival * (Lb[iq] * grdPhi[j]);
+ mat[i][j] += quadrature->getWeight(iq) * psival * dot(Lb[iq], grdPhi[j]);
}
}
}
@@ -376,27 +368,30 @@ namespace AMDiS {
}
int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints);
+ for (int iq = 0; iq < nPoints; iq++) {
+ Lb[iq].change_dim(dim + 1);
+ Lb[iq] = 0.0;
+ }
- for (int iq = 0; iq < nPoints; iq++)
- Lb[iq].set(0.0);
+ for (int i = 0; i < static_cast<int>(terms.size()); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
+ const mtl::dense2D<double>& psi = psiFast->getPhi();
+ mtl::dense2D<double> facMat(nPoints, nCol);
for (int iq = 0; iq < nPoints; iq++) {
- Lb[iq] *= elInfo->getDet();
- const double *psi = psiFast->getPhi(iq);
+ double weight = quadrature->getWeight(iq);
+ mtl::dense_vector<double>& Lb_iq = Lb[iq];
+ const vector<mtl::dense_vector<double> >& grdPsi = phiFast->getGradient(iq);
- for (int i = 0; i < nCol; i++) {
- double factor =
- quadrature->getWeight(iq) * (Lb[iq] * phiFast->getGradient(iq, i));
- for (int j = 0; j < nRow; j++)
- mat[j][i] += factor * psi[j];
- }
+ Lb_iq *= elInfo->getDet();
+
+ for (int i = 0; i < nCol; i++)
+ facMat[iq][i] = weight * dot(Lb_iq, grdPsi[i]);
}
+
+ mat += trans(psi) * facMat;
}
@@ -426,14 +421,11 @@ namespace AMDiS {
}
const int **nEntries = q01->getNumberEntries();
- int myRank = omp_get_thread_num();
// Do not need to resize Lb, because it's size is always at least one!
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
- Lb[0].set(0.0);
+ Lb[0] = 0.0;
- for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
- (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
- }
+ for (int i = 0; i < static_cast<int>( terms.size()); i++)
+ (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, Lb);
Lb[0] *= elInfo->getDet();
@@ -470,13 +462,11 @@ namespace AMDiS {
}
const int *nEntries = q1->getNumberEntries();
- int myRank = omp_get_thread_num();
// Do not need to resize Lb, because it's size is always at least one!
- VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
- Lb[0].set(0.0);
+ Lb[0] = 0.0;
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
- (static_cast<FirstOrderTerm*>(terms[myRank][i]))->getLb(elInfo, 1, Lb);
+ for (int i = 0; i < static_cast<int>(terms.size()); i++)
+ (static_cast<FirstOrderTerm*>(terms[i]))->getLb(elInfo, Lb);
Lb[0] *= elInfo->getDet();
diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h
index aa398fd6dba51fbf56e67d9a9e4a9d2e625d8025..9f9a4347706798ca3ccfbbb3fe00b2b7da5e71d3 100644
--- a/AMDiS/src/FirstOrderAssembler.h
+++ b/AMDiS/src/FirstOrderAssembler.h
@@ -63,11 +63,8 @@ namespace AMDiS {
protected:
- /** \brief
- * Thread safe temporary vector of DimMats for calculation in
- * function calculateElementMatrix().
- */
- std::vector<VectorOfFixVecs<DimVec<double> > > tmpLb;
+ /// Vector of DimMats for calculation in function calculateElementMatrix().
+ std::vector<mtl::dense_vector<double> > Lb;
/// List of all yet created optimized zero order assemblers for grdPsi.
static std::vector<SubAssembler*> optimizedSubAssemblersGrdPsi;
@@ -128,7 +125,7 @@ namespace AMDiS {
}
protected:
- std::vector<VectorOfFixVecs<DimVec<double> > > tmpGrdPhi;
+ vector<mtl::dense_vector<double> > grdPhi;
const BasisFunction *psi, *phi;
};
diff --git a/AMDiS/src/FirstOrderTerm.cc b/AMDiS/src/FirstOrderTerm.cc
index ec34f630c40895c1bea0889300a08cb2770cee01..206959571630a43abf702ddff2154542f83d5d49 100644
--- a/AMDiS/src/FirstOrderTerm.cc
+++ b/AMDiS/src/FirstOrderTerm.cc
@@ -53,20 +53,21 @@ namespace AMDiS {
getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs);
}
- void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void VecAtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
- lb_one(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
+ lb_one(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
- l1(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
+ l1(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
}
}
@@ -101,10 +102,12 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void CoordsAtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], (*g)(coordsAtQPs[iq]));
}
@@ -140,10 +143,12 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], (*g)(coordsAtQPs[iq]));
}
@@ -188,17 +193,19 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
- void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void VectorGradient_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- if (f) {
+ const int nPoints = static_cast<int>(Lb.size());
+
+ if (f)
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
- } else {
+ else
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, gradAtQPs[iq], Lb[iq], 1.0);
- }
+
}
void VectorGradient_FOT::eval(int nPoints,
@@ -240,10 +247,12 @@ namespace AMDiS {
getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
}
- void VectorFct_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void VectorFct_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
}
@@ -274,10 +283,12 @@ namespace AMDiS {
}
- void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void VecFctAtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
}
@@ -333,9 +344,12 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad);
}
- void VecGrad_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const {
+ void VecGrad_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
+ {
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
}
@@ -366,10 +380,12 @@ namespace AMDiS {
grad2AtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
}
- void FctGrad2_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void FctGrad2_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
}
@@ -414,8 +430,12 @@ namespace AMDiS {
gradAtQPs1 = getGradientsAtQPs(vec1, elInfo, subAssembler, quad);
}
- void Vec2AndGradAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ void Vec2AndGradAtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
+ {
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++) {
if (b)
lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq],gradAtQPs1[iq]));
@@ -458,10 +478,12 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void FctVecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void FctVecAtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++) {
if (b)
lb(grdLambda, *b, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
@@ -519,10 +541,11 @@ namespace AMDiS {
subAssembler->getVectorAtQPs(vec2, elInfo, quad, vec2AtQPs);
}
- void Vec2AtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void Vec2AtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
@@ -594,10 +617,12 @@ namespace AMDiS {
subAssembler->getVectorAtQPs(vec3, elInfo, quad, vec3AtQPs);
}
- void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
@@ -627,12 +652,12 @@ namespace AMDiS {
// =========== General_FOT ==========
- General_FOT::General_FOT(std::vector<DOFVectorBase<double>*> vecs,
- std::vector<DOFVectorBase<double>*> grads,
+ General_FOT::General_FOT(vector<DOFVectorBase<double>*> vecs,
+ vector<DOFVectorBase<double>*> grads,
TertiaryAbstractFunction<WorldVector<double>,
WorldVector<double>,
- std::vector<double>,
- std::vector<WorldVector<double> > > *af)
+ vector<double>,
+ vector<WorldVector<double> > > *af)
: FirstOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
{
vecsAtQPs.resize(vecs_.size());
@@ -651,6 +676,7 @@ namespace AMDiS {
}
}
+
void General_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -666,17 +692,18 @@ namespace AMDiS {
gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
+
void General_FOT::getLb(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >& result) const
+ vector<mtl::dense_vector<double> >& Lb) const
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
- std::vector<double> vecsArg(nVecs);
- std::vector<WorldVector<double> > gradsArg(nGrads);
+ vector<double> vecsArg(nVecs);
+ vector<WorldVector<double> > gradsArg(nGrads);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
for (int iq = 0; iq < nPoints; iq++) {
for (int i = 0; i < nVecs; i++)
@@ -684,11 +711,11 @@ namespace AMDiS {
for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
- lb(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
- result[iq], 1.0);
+ lb(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), Lb[iq], 1.0);
}
}
+
void General_FOT::eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -700,8 +727,8 @@ namespace AMDiS {
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
- std::vector<double> vecsArg(nVecs);
- std::vector<WorldVector<double> > gradsArg(nGrads);
+ vector<double> vecsArg(nVecs);
+ vector<WorldVector<double> > gradsArg(nGrads);
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
@@ -725,13 +752,13 @@ namespace AMDiS {
// =========== GeneralParametric_FOT ==========
- GeneralParametric_FOT::GeneralParametric_FOT(std::vector<DOFVectorBase<double>*> vecs,
- std::vector<DOFVectorBase<double>*> grads,
+ GeneralParametric_FOT::GeneralParametric_FOT(vector<DOFVectorBase<double>*> vecs,
+ vector<DOFVectorBase<double>*> grads,
QuartAbstractFunction<WorldVector<double>,
WorldVector<double>,
WorldVector<double>,
- std::vector<double>,
- std::vector<WorldVector<double> > > *af)
+ vector<double>,
+ vector<WorldVector<double> > > *af)
: FirstOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
{
vecsAtQPs.resize(vecs_.size());
@@ -769,16 +796,17 @@ namespace AMDiS {
void GeneralParametric_FOT::getLb(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >& result) const
+ vector<mtl::dense_vector<double> >& Lb) const
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
- std::vector<double> vecsArg(nVecs);
- std::vector<WorldVector<double> > gradsArg(nGrads);
+ vector<double> vecsArg(nVecs);
+ vector<WorldVector<double> > gradsArg(nGrads);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
+
for (int iq = 0; iq < nPoints; iq++) {
for (int i = 0; i < nVecs; i++)
vecsArg[i] = vecsAtQPs[i][iq];
@@ -786,7 +814,7 @@ namespace AMDiS {
gradsArg[i] = gradsAtQPs_[i][iq];
lb(grdLambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg),
- result[iq], 1.0);
+ Lb[iq], 1.0);
}
}
@@ -802,8 +830,8 @@ namespace AMDiS {
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
- std::vector<double> vecsArg(nVecs);
- std::vector<WorldVector<double> > gradsArg(nGrads);
+ vector<double> vecsArg(nVecs);
+ vector<WorldVector<double> > gradsArg(nGrads);
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
diff --git a/AMDiS/src/FirstOrderTerm.h b/AMDiS/src/FirstOrderTerm.h
index ace037e00314c9563cad9d14c33b0b83ab52be98..1eb3ce622b208a916daa0de43a77e74fffb85b05 100644
--- a/AMDiS/src/FirstOrderTerm.h
+++ b/AMDiS/src/FirstOrderTerm.h
@@ -23,6 +23,7 @@
#ifndef AMDIS_FIRST_ORDER_TERM_H
#define AMDIS_FIRST_ORDER_TERM_H
+#include <vector>
#include "AMDiS_fwd.h"
#include "OperatorTerm.h"
#include "AbstractFunction.h"
@@ -30,6 +31,8 @@
namespace AMDiS {
+ using namespace std;
+
/**
* \ingroup Assembler
*
@@ -49,8 +52,7 @@ namespace AMDiS {
/// Evaluation of \f$ \Lambda b \f$.
virtual void getLb(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >& result) const = 0;
+ vector<mtl::dense_vector<double> >& result) const = 0;
/// Implenetation of FirstOrderTerm::eval().
void eval(int nPoints,
@@ -79,7 +81,7 @@ namespace AMDiS {
* each component.
*/
inline void l1(const DimVec<WorldVector<double> >& Lambda,
- DimVec<double>& Lb,
+ mtl::dense_vector<double>& Lb,
double factor) const
{
const int dim = Lambda.getSize();
@@ -97,7 +99,7 @@ namespace AMDiS {
/// Evaluation of \f$ \Lambda \cdot b\f$.
inline void lb(const DimVec<WorldVector<double> >& Lambda,
const WorldVector<double>& b,
- DimVec<double>& Lb,
+ mtl::dense_vector<double>& Lb,
double factor) const
{
const int dim = Lambda.getSize();
@@ -114,7 +116,7 @@ namespace AMDiS {
/// Evaluation of \f$ \Lambda \cdot b\f$.
inline void lb_one(const DimVec<WorldVector<double> >& Lambda,
- DimVec<double>& Lb,
+ mtl::dense_vector<double>& Lb,
double factor) const
{
const int dim = Lambda.getSize();
@@ -142,10 +144,10 @@ namespace AMDiS {
/// Implements FirstOrderTerm::getLb().
inline void getLb(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], 1.0);
@@ -176,9 +178,11 @@ namespace AMDiS {
{}
/// Implements FirstOrderTerm::getLb().
- inline void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
+ inline void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], (*factor));
@@ -211,11 +215,11 @@ namespace AMDiS {
}
/// Implements FirstOrderTerm::getLb().
- inline void getLb(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >&Lb) const
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(Lb.size());
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
@@ -274,8 +278,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -318,7 +322,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FistOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints,VectorOfFixVecs<DimVec<double> >&Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -356,8 +361,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FistOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -396,8 +401,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -435,8 +440,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int qPoint,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -476,8 +481,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FistOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >&Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -520,8 +525,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int qPoint,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -566,8 +571,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int qPoint,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -607,8 +612,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
void getLb(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ vector<mtl::dense_vector<double> >& Lb) const;
void eval(int nPoints,
const ElementVector& uhAtQP,
@@ -640,8 +644,8 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
void eval(int nPoints,
const ElementVector& uhAtQP,
@@ -673,8 +677,8 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
void eval(int nPoints,
const ElementVector& uhAtQP,
@@ -708,8 +712,8 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& Lb) const;
void eval(int nPoints,
const ElementVector& uhAtQP,
@@ -733,12 +737,12 @@ namespace AMDiS {
{
public:
/// Constructor
- General_FOT(std::vector<DOFVectorBase<double>*> vecs,
- std::vector<DOFVectorBase<double>*> grads,
+ General_FOT(vector<DOFVectorBase<double>*> vecs,
+ vector<DOFVectorBase<double>*> grads,
TertiaryAbstractFunction<WorldVector<double>,
WorldVector<double>,
- std::vector<double>,
- std::vector<WorldVector<double> > > *f);
+ vector<double>,
+ vector<WorldVector<double> > > *f);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo*,
@@ -746,8 +750,8 @@ namespace AMDiS {
Quadrature *quad= NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& result) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& result) const;
/// Implenetation of FirstOrderTerm::eval().
void eval(int nPoints,
@@ -758,20 +762,20 @@ namespace AMDiS {
double factor);
protected:
- std::vector<DOFVectorBase<double>*> vecs_;
+ vector<DOFVectorBase<double>*> vecs_;
- std::vector<DOFVectorBase<double>*> grads_;
+ vector<DOFVectorBase<double>*> grads_;
TertiaryAbstractFunction<WorldVector<double>,
WorldVector<double>,
- std::vector<double>,
- std::vector<WorldVector<double> > > *f_;
+ vector<double>,
+ vector<WorldVector<double> > > *f_;
WorldVector<double> *coordsAtQPs_;
- std::vector<mtl::dense_vector<double> > vecsAtQPs;
+ vector<mtl::dense_vector<double> > vecsAtQPs;
- std::vector<WorldVector<double>*> gradsAtQPs_;
+ vector<WorldVector<double>*> gradsAtQPs_;
};
@@ -779,20 +783,20 @@ namespace AMDiS {
{
public:
/// Constructor
- GeneralParametric_FOT(std::vector<DOFVectorBase<double>*> vecs,
- std::vector<DOFVectorBase<double>*> grads,
+ GeneralParametric_FOT(vector<DOFVectorBase<double>*> vecs,
+ vector<DOFVectorBase<double>*> grads,
QuartAbstractFunction<WorldVector<double>,
WorldVector<double>,
WorldVector<double>,
- std::vector<double>,
- std::vector<WorldVector<double> > > *f);
+ vector<double>,
+ vector<WorldVector<double> > > *f);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints,
- VectorOfFixVecs<DimVec<double> >& result) const;
+ void getLb(const ElInfo *elInfo,
+ vector<mtl::dense_vector<double> >& result) const;
/// Implenetation of FirstOrderTerm::eval().
void eval(int nPoints,
@@ -803,22 +807,22 @@ namespace AMDiS {
double factor);
protected:
- std::vector<DOFVectorBase<double>*> vecs_;
+ vector<DOFVectorBase<double>*> vecs_;
- std::vector<DOFVectorBase<double>*> grads_;
+ vector<DOFVectorBase<double>*> grads_;
QuartAbstractFunction<WorldVector<double>,
WorldVector<double>,
WorldVector<double>,
- std::vector<double>,
- std::vector<WorldVector<double> > > *f_;
+ vector<double>,
+ vector<WorldVector<double> > > *f_;
WorldVector<double> *coordsAtQPs_;
WorldVector<double> elementNormal_;
- std::vector<mtl::dense_vector<double> > vecsAtQPs;
+ vector<mtl::dense_vector<double> > vecsAtQPs;
- std::vector<WorldVector<double>*> gradsAtQPs_;
+ vector<WorldVector<double>*> gradsAtQPs_;
};
}
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index 16c2e98aaa38360b6eea74dc65b9818b22f0c170..79c5e38b4679b375087117243aa9981fc21f41fe 100644
--- a/AMDiS/src/Lagrange.h
+++ b/AMDiS/src/Lagrange.h
@@ -23,10 +23,11 @@
#ifndef AMDIS_LAGRANGE_H
#define AMDIS_LAGRANGE_H
+#include <list>
+#include <boost/numeric/mtl/mtl.hpp>
+#include "AbstractFunction.h"
#include "BasisFunction.h"
#include "FixVec.h"
-#include "AbstractFunction.h"
-#include <list>
namespace AMDiS {
@@ -416,10 +417,12 @@ namespace AMDiS {
private:
int* vertices;
- void (*func)(const DimVec<double>& lambda, int* vertices_,
- DimVec<double>& result);
+ void (*func)(const DimVec<double>& lambda,
+ int* vertices_,
+ mtl::dense_vector<double>& result);
- inline void operator()(const DimVec<double>& lambda, DimVec<double>& result) const
+ inline void operator()(const DimVec<double>& lambda,
+ mtl::dense_vector<double>& result) const
{
func(lambda, vertices, result);
}
@@ -428,18 +431,18 @@ namespace AMDiS {
// center
inline static void grdPhi0c(const DimVec<double>&,
int*,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
}
// ====== Lagrange1 ================================================
// vertex
inline static void grdPhi1v(const DimVec<double>&,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = 1.0;
}
@@ -447,18 +450,18 @@ namespace AMDiS {
// vertex
inline static void grdPhi2v(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = 4.0 * lambda[vertices[0]] - 1.0;
}
// edge
inline static void grdPhi2e(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = 4.0 * lambda[vertices[1]];
result[vertices[1]] = 4.0 * lambda[vertices[0]];
}
@@ -467,9 +470,9 @@ namespace AMDiS {
// vertex
inline static void grdPhi3v(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = (13.5 * lambda[vertices[0]] - 9.0) *
lambda[vertices[0]] + 1.0;
}
@@ -477,9 +480,9 @@ namespace AMDiS {
// edge
inline static void grdPhi3e(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = (27.0 * lambda[vertices[0]] - 4.5) *
lambda[vertices[1]];
result[vertices[1]] = (13.5 * lambda[vertices[0]] - 4.5) *
@@ -489,9 +492,9 @@ namespace AMDiS {
// face
inline static void grdPhi3f(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = 27.0 * lambda[vertices[1]] * lambda[vertices[2]];
result[vertices[1]] = 27.0 * lambda[vertices[0]] * lambda[vertices[2]];
result[vertices[2]] = 27.0 * lambda[vertices[0]] * lambda[vertices[1]];
@@ -502,9 +505,9 @@ namespace AMDiS {
// vertex
inline static void grdPhi4v(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] =
((128.0 * lambda[vertices[0]] - 144.0) * lambda[vertices[0]] + 44.0) *
lambda[vertices[0]] / 3.0 - 1.0;
@@ -513,9 +516,9 @@ namespace AMDiS {
// edge
inline static void grdPhi4e0(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = ((128.0 * lambda[vertices[0]] - 64.0) *
lambda[vertices[0]] + 16.0 / 3.0) * lambda[vertices[1]];
result[vertices[1]] = ((128.0 * lambda[vertices[0]] - 96.0) *
@@ -524,9 +527,9 @@ namespace AMDiS {
inline static void grdPhi4e1(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = 4.0 * (8.0 * lambda[vertices[0]] - 1.0) *
lambda[vertices[1]] * (4.0 * lambda[vertices[1]] - 1.0);
result[vertices[1]] = 4.0 * lambda[vertices[0]] *
@@ -536,9 +539,9 @@ namespace AMDiS {
// face
inline static void grdPhi4f(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
+ mtl::dense_vector<double>& result)
{
- result.set(0.0);
+ result = 0.0;
result[vertices[0]] = 32.0 * (8.0 * lambda[vertices[0]] - 1.0) *
lambda[vertices[1]] * lambda[vertices[2]];
result[vertices[1]] = 32.0 * (4.0 * lambda[vertices[0]] - 1.0) *
@@ -550,13 +553,17 @@ namespace AMDiS {
// center
inline static void grdPhi4c(const DimVec<double>& lambda,
int* vertices,
- DimVec<double>& result)
- {
- result.set(0.0);
- result[0] = 256.0 * lambda[vertices[1]] * lambda[vertices[2]] * lambda[vertices[3]];
- result[1] = 256.0 * lambda[vertices[0]] * lambda[vertices[2]] * lambda[vertices[3]];
- result[2] = 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[3]];
- result[3] = 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[2]];
+ mtl::dense_vector<double>& result)
+ {
+ result = 0.0;
+ result[0] =
+ 256.0 * lambda[vertices[1]] * lambda[vertices[2]] * lambda[vertices[3]];
+ result[1] =
+ 256.0 * lambda[vertices[0]] * lambda[vertices[2]] * lambda[vertices[3]];
+ result[2] =
+ 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[3]];
+ result[3] =
+ 256.0 * lambda[vertices[0]] * lambda[vertices[1]] * lambda[vertices[2]];
}
};
diff --git a/AMDiS/src/MatrixVector.h b/AMDiS/src/MatrixVector.h
index e7e1525658d0a1e520c88149e6ec989cc24675ae..26d3eb02675ed4bfc30db5a04d9869895afcc220 100644
--- a/AMDiS/src/MatrixVector.h
+++ b/AMDiS/src/MatrixVector.h
@@ -352,49 +352,6 @@ namespace AMDiS {
return result;
}
- /// Scaled matrix vector multiplication: result = a * (m * v)
- template<typename T>
- inline const Vector<T>& amv(double a, const Matrix<T>& m, const Vector<T>& v,
- Vector<T>& result)
- {
- TEST_EXIT_DBG(m.getNumCols() == v.getSize())("m and v not compatible\n");
- TEST_EXIT_DBG(v.getSize() == result.getSize())("wrong result size\n");
-
- T *resultIt, *mIt, *vIt;
-
- for (resultIt = result.begin(), mIt = m.begin();
- resultIt != result.end();
- ++resultIt) {
- *resultIt = 0;
- for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt)
- *resultIt += *mIt * *vIt;
- *resultIt *= a;
- }
-
- return result;
- }
-
- /// Computes x(Ay)^T, with x and y vectors, and A a matrix.
- template<typename T>
- inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y)
- {
- TEST_EXIT_DBG(a.getNumCols() == x.getSize())("A and x not compatible\n");
- TEST_EXIT_DBG(a.getNumCols() == y.getSize())("A and y not compatible\n");
-
- T result = 0;
- T tmp = 0;
- T *aIt, *xIt, *yIt;
-
- for (aIt = a.begin(), xIt = x.begin(); aIt != a.end(); ++xIt) {
- tmp = 0;
- for (yIt = y.begin(); yIt != y.end(); ++yIt, ++aIt)
- tmp += *aIt * *yIt;
- result += *xIt * tmp;
- }
-
- return result;
- }
-
/// Matrix vector multiplication.
template<typename T>
inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m)
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 57b8d76d6ea3429d23e858cab9464ff5fa95e9e3..08fdcc2285e1aa3bc4fadae978aff84af1f38257 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -16,7 +16,6 @@
#include "FixVec.h"
#include "DOFVector.h"
#include "Quadrature.h"
-#include "OpenMP.h"
namespace AMDiS {
@@ -27,20 +26,19 @@ namespace AMDiS {
int Operator::getQuadratureDegree(int order, FirstOrderType firstOrderType)
{
std::vector<OperatorTerm*>* terms = NULL;
- int myRank = omp_get_thread_num();
switch(order) {
case 0:
- terms = &zeroOrder[myRank];
+ terms = &zeroOrder;
break;
case 1:
if (firstOrderType == GRD_PHI)
- terms = &firstOrderGrdPhi[myRank];
+ terms = &firstOrderGrdPhi;
else
- terms = &firstOrderGrdPsi[myRank];
+ terms = &firstOrderGrdPsi;
break;
case 2:
- terms = &secondOrder[myRank];
+ terms = &secondOrder;
break;
}
@@ -84,21 +82,11 @@ namespace AMDiS {
uhOld(NULL),
optimized(true)
{
- int maxThreads = omp_get_overall_max_threads();
-
- assembler.resize(maxThreads);
- secondOrder.resize(maxThreads);
- firstOrderGrdPsi.resize(maxThreads);
- firstOrderGrdPhi.resize(maxThreads);
- zeroOrder.resize(maxThreads);
-
- for (int i = 0; i < maxThreads; i++) {
- assembler[i] = NULL;
- secondOrder[i].resize(0);
- firstOrderGrdPsi[i].resize(0);
- firstOrderGrdPhi[i].resize(0);
- zeroOrder[i].resize(0);
- }
+ assembler = NULL;
+ secondOrder.resize(0);
+ firstOrderGrdPsi.resize(0);
+ firstOrderGrdPhi.resize(0);
+ zeroOrder.resize(0);
nRow = rowFeSpace->getBasisFcts()->getNumber();
nCol = colFeSpace->getBasisFcts()->getNumber();
@@ -116,12 +104,10 @@ namespace AMDiS {
ElementMatrix& userMat,
double factor)
{
- int myRank = omp_get_thread_num();
-
- if (!assembler[myRank])
- initAssembler(myRank, NULL, NULL, NULL, NULL);
+ if (!assembler)
+ initAssembler(NULL, NULL, NULL, NULL);
- assembler[myRank]->calculateElementMatrix(elInfo, userMat, factor);
+ assembler->calculateElementMatrix(elInfo, userMat, factor);
}
@@ -131,15 +117,13 @@ namespace AMDiS {
ElementMatrix& userMat,
double factor)
{
- int myRank = omp_get_thread_num();
-
- if (!assembler[myRank])
- initAssembler(myRank, NULL, NULL, NULL, NULL);
+ if (!assembler)
+ initAssembler(NULL, NULL, NULL, NULL);
- assembler[myRank]->calculateElementMatrix(rowElInfo, colElInfo,
- smallElInfo, largeElInfo,
- rowColFeSpaceEqual,
- userMat, factor);
+ assembler->calculateElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo,
+ rowColFeSpaceEqual,
+ userMat, factor);
}
@@ -147,12 +131,10 @@ namespace AMDiS {
ElementVector& userVec,
double factor)
{
- int myRank = omp_get_thread_num();
+ if (!assembler)
+ initAssembler(NULL, NULL, NULL, NULL);
- if (!assembler[myRank])
- initAssembler(myRank, NULL, NULL, NULL, NULL);
-
- assembler[myRank]->calculateElementVector(elInfo, userVec, factor);
+ assembler->calculateElementVector(elInfo, userVec, factor);
}
@@ -161,68 +143,55 @@ namespace AMDiS {
ElementVector& userVec,
double factor)
{
- int myRank = omp_get_thread_num();
-
- if (!assembler[myRank])
- initAssembler(myRank, NULL, NULL, NULL, NULL);
+ if (!assembler)
+ initAssembler(NULL, NULL, NULL, NULL);
- assembler[myRank]->calculateElementVector(mainElInfo, auxElInfo,
- smallElInfo, largeElInfo,
- userVec, factor);
+ assembler->calculateElementVector(mainElInfo, auxElInfo,
+ smallElInfo, largeElInfo,
+ userVec, factor);
}
- void Operator::initAssembler(int rank,
- Quadrature *quad2,
+ void Operator::initAssembler(Quadrature *quad2,
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0)
{
-#ifdef _OPENMP
-#pragma omp critical (initAssembler)
-#endif
- {
- if (optimized)
- assembler[rank] =
- new OptimizedAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFeSpace, colFeSpace);
- else
- assembler[rank] =
- new StandardAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFeSpace, colFeSpace);
- }
+ if (optimized)
+ assembler = new OptimizedAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
+ rowFeSpace, colFeSpace);
+ else
+ assembler = new StandardAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
+ rowFeSpace, colFeSpace);
}
void Operator::finishAssembling()
{
- assembler[omp_get_thread_num()]->finishAssembling();
+ assembler->finishAssembling();
}
- Assembler *Operator::getAssembler(int rank)
+ Assembler *Operator::getAssembler()
{
- if (!assembler[rank])
- initAssembler(rank, NULL, NULL, NULL, NULL);
+ if (!assembler)
+ initAssembler(NULL, NULL, NULL, NULL);
- return assembler[rank];
+ return assembler;
}
- void Operator::setAssembler(int rank, Assembler *a)
+ void Operator::setAssembler(Assembler *a)
{
- assembler[rank] = a;
+ assembler = a;
}
void Operator::weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
- int myRank = omp_get_thread_num();
-
- for (std::vector<OperatorTerm*>::const_iterator termIt = secondOrder[myRank].begin();
- termIt != secondOrder[myRank].end();
- ++termIt)
+ for (std::vector<OperatorTerm*>::const_iterator termIt = secondOrder.begin();
+ termIt != secondOrder.end(); ++termIt)
static_cast<SecondOrderTerm*>(*termIt)->weakEval(grdUhAtQP, result);
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 368e2b04adb25dd6cae2ce705ea00716426af291..cff634d535c2e25571d1d528e66a22a806f86207 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -39,6 +39,8 @@
namespace AMDiS {
+ using namespace std;
+
/** \brief
* An Operator holds all information needed to assemble the system matrix
* and the right hand side. It consists of four OperatorTerm lists each storing
@@ -152,10 +154,10 @@ namespace AMDiS {
}
/// Returns \ref assembler
- Assembler *getAssembler(int rank = 0);
+ Assembler *getAssembler();
/// Sets \ref assembler
- void setAssembler(int rank, Assembler *ass);
+ void setAssembler(Assembler *ass);
/// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
inline void setFillFlag(Flag f)
@@ -182,8 +184,7 @@ namespace AMDiS {
}
/// Initialization of the needed SubAssemblers using the given quadratures.
- void initAssembler(int rank,
- Quadrature *quad2,
+ void initAssembler(Quadrature *quad2,
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0);
@@ -200,12 +201,8 @@ namespace AMDiS {
ElementVector& result,
double factor) const
{
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = zeroOrder[myRank].begin();
- termIt != zeroOrder[myRank].end();
- ++termIt)
+ for (vector<OperatorTerm*>::const_iterator termIt = zeroOrder.begin();
+ termIt != zeroOrder.end(); ++termIt)
(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
@@ -216,13 +213,9 @@ namespace AMDiS {
const WorldMatrix<double> *D2UhAtQP,
ElementVector& result,
double factor) const
- {
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = firstOrderGrdPsi[myRank].begin();
- termIt != firstOrderGrdPsi[myRank].end();
- ++termIt)
+ {
+ for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPsi.begin();
+ termIt != firstOrderGrdPsi.end(); ++termIt)
(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
@@ -234,12 +227,8 @@ namespace AMDiS {
ElementVector& result,
double factor) const
{
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = firstOrderGrdPhi[myRank].begin();
- termIt != firstOrderGrdPhi[myRank].end();
- ++termIt)
+ for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPhi.begin();
+ termIt != firstOrderGrdPhi.end(); ++termIt)
(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
@@ -250,94 +239,72 @@ namespace AMDiS {
const WorldMatrix<double> *D2UhAtQP,
ElementVector& result,
double factor) const
- {
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = secondOrder[myRank].begin();
- termIt != secondOrder[myRank].end();
- ++termIt)
+ {
+ for (vector<OperatorTerm*>::const_iterator termIt = secondOrder.begin();
+ termIt != secondOrder.end(); ++termIt)
(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
}
/// Weak evaluation of all terms in \ref secondOrder.
- void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
- std::vector<WorldVector<double> > &result) const;
+ void weakEvalSecondOrder(const vector<WorldVector<double> > &grdUhAtQP,
+ vector<WorldVector<double> > &result) const;
/// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
- {
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = secondOrder[myRank].begin();
- termIt != secondOrder[myRank].end();
- ++termIt)
- static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
+ void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
+ {
+ for (vector<OperatorTerm*>::const_iterator termIt = secondOrder.begin();
+ termIt != secondOrder.end(); ++termIt)
+ static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, LALt);
}
/// Calls getLb() for each term in \ref firstOrderGrdPsi and adds the results to Lb.
void getLbGrdPsi(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
- {
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = firstOrderGrdPsi[myRank].begin();
- termIt != firstOrderGrdPsi[myRank].end();
- ++termIt)
- static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
+ vector<mtl::dense_vector<double> >& Lb) const
+ {
+ for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPsi.begin();
+ termIt != firstOrderGrdPsi.end(); ++termIt)
+ static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, Lb);
}
/// Calls getLb() for each term in \ref firstOrderGrdPhi and adds the results to Lb.
void getLbGrdPhi(const ElInfo *elInfo,
- int nPoints,
- VectorOfFixVecs<DimVec<double> >& Lb) const
- {
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = firstOrderGrdPhi[myRank].begin();
- termIt != firstOrderGrdPhi[myRank].end();
- ++termIt)
- static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
+ vector<mtl::dense_vector<double> >& Lb) const
+ {
+ for (vector<OperatorTerm*>::const_iterator termIt = firstOrderGrdPhi.begin();
+ termIt != firstOrderGrdPhi.end(); ++termIt)
+ static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, Lb);
}
/// Calls getC() for each term in \ref zeroOrder and adds the results to c.
void getC(const ElInfo *elInfo, int nPoints, ElementVector& c)
- {
- int myRank = omp_get_thread_num();
-
- std::vector<OperatorTerm*>::const_iterator termIt;
- for (termIt = zeroOrder[myRank].begin();
- termIt != zeroOrder[myRank].end();
- ++termIt)
+ {
+ for (vector<OperatorTerm*>::const_iterator termIt = zeroOrder.begin();
+ termIt != zeroOrder.end(); ++termIt)
static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
}
/// Returns true, if there are second order terms. Returns false otherwise.
inline bool secondOrderTerms()
{
- return secondOrder[omp_get_thread_num()].size() != 0;
+ return secondOrder.size() != 0;
}
/// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
inline bool firstOrderTermsGrdPsi()
{
- return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
+ return firstOrderGrdPsi.size() != 0;
}
/// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
inline bool firstOrderTermsGrdPhi()
{
- return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
+ return firstOrderGrdPhi.size() != 0;
}
/// Returns true, if there are zero order terms. Returns false otherwise.
inline bool zeroOrderTerms()
{
- return zeroOrder[omp_get_thread_num()].size() != 0;
+ return zeroOrder.size() != 0;
}
public:
@@ -374,19 +341,19 @@ namespace AMDiS {
* created especially for this Operator, when \ref getElementMatrix()
* or \ref getElementVector is called for the first time.
*/
- std::vector<Assembler*> assembler;
+ Assembler* assembler;
/// List of all second order terms
- std::vector< std::vector<OperatorTerm*> > secondOrder;
+ vector<OperatorTerm*> secondOrder;
/// List of all first order terms derived to psi
- std::vector< std::vector<OperatorTerm*> > firstOrderGrdPsi;
+ vector<OperatorTerm*> firstOrderGrdPsi;
/// List of all first order terms derived to phi
- std::vector< std::vector<OperatorTerm*> > firstOrderGrdPhi;
+ vector<OperatorTerm*> firstOrderGrdPhi;
/// List of all zero order terms
- std::vector< std::vector<OperatorTerm*> > zeroOrder;
+ vector<OperatorTerm*> zeroOrder;
/** \brief
* Pointer to the solution of the last timestep. Can be used for a more efficient
diff --git a/AMDiS/src/Operator.hh b/AMDiS/src/Operator.hh
index 0399065160bc8fe026250f212fbf4de45b4c8492..61e92dfe855c4f03730865263e88040dfda9d5a8 100644
--- a/AMDiS/src/Operator.hh
+++ b/AMDiS/src/Operator.hh
@@ -15,16 +15,11 @@ namespace AMDiS {
template <typename T>
void Operator::addSecondOrderTerm(T *term)
{
- secondOrder[0].push_back(term);
+ secondOrder.push_back(term);
term->operat = this;
for (std::set<const FiniteElemSpace*>::iterator it = term->getAuxFeSpaces().begin();
it != term->getAuxFeSpaces().end(); ++it)
auxFeSpaces.insert(*it);
-
- for (int i = 1; i < omp_get_overall_max_threads(); i++) {
- T *newTerm = new T(static_cast<const T>(*term));
- secondOrder[i].push_back(newTerm);
- }
}
template <typename T>
@@ -32,38 +27,24 @@ namespace AMDiS {
FirstOrderType type)
{
if (type == GRD_PSI) {
- firstOrderGrdPsi[0].push_back(term);
+ firstOrderGrdPsi.push_back(term);
} else {
- firstOrderGrdPhi[0].push_back(term);
+ firstOrderGrdPhi.push_back(term);
}
term->operat = this;
for (std::set<const FiniteElemSpace*>::iterator it = term->getAuxFeSpaces().begin();
it != term->getAuxFeSpaces().end(); ++it)
auxFeSpaces.insert(*it);
-
- for (int i = 1; i < omp_get_overall_max_threads(); i++) {
- T *newTerm = new T(static_cast<const T>(*term));
- if (type == GRD_PSI)
- firstOrderGrdPsi[i].push_back(newTerm);
- else
- firstOrderGrdPhi[i].push_back(newTerm);
- }
}
template <typename T>
void Operator::addZeroOrderTerm(T *term)
{
- zeroOrder[0].push_back(term);
+ zeroOrder.push_back(term);
term->operat = this;
for (std::set<const FiniteElemSpace*>::iterator it = term->getAuxFeSpaces().begin();
it != term->getAuxFeSpaces().end(); ++it)
auxFeSpaces.insert(*it);
-
- for (int i = 1; i < omp_get_overall_max_threads(); i++) {
- T *newTerm = new T(static_cast<const T>(*term));
- zeroOrder[i].push_back(newTerm);
- }
}
-
}
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 63710b1108518a1d9f2051a71a5d23dac78769b1..5e7ce79102dc879dae2e8d297738bad7709a7258 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -651,6 +651,9 @@ namespace AMDiS {
{
FUNCNAME("ProblemVec::buildAfterCoarsen()");
+// buildAfterCoarsen_sebastianMode(adaptInfo, flag);
+// return;
+
if (dualMeshTraverseRequired()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
ERROR_EXIT("Dual mesh assemble does not work in parallel code!\n");
@@ -777,13 +780,12 @@ namespace AMDiS {
assembledMatrix[i][j] = true;
-
if (assembleMatrix)
matrix->finishInsertion();
if (assembleMatrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->exitMatrix(matrix);
-
+
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
}
@@ -817,6 +819,196 @@ namespace AMDiS {
#endif
}
+ void ProblemVec::buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag)
+ {
+ FUNCNAME("ProblemVec::buildAfterCoarsen()");
+
+ clock_t first = clock();
+#ifdef _OPENMP
+ double wtime = omp_get_wtime();
+#endif
+
+ for (int i = 0; i < static_cast<int>(meshes.size()); i++)
+ meshes[i]->dofCompress();
+
+ Flag assembleFlag =
+ flag |
+ (*systemMatrix)[0][0]->getAssembleFlag() |
+ rhs->getDOFVector(0)->getAssembleFlag() |
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_NEIGH;
+
+ if (useGetBound)
+ assembleFlag |= Mesh::FILL_BOUND;
+
+ traverseInfo.updateStatus();
+
+ // Used to calculate the overall number of non zero entries.
+ int nnz = 0;
+
+
+ /// === INITIALIZE ===
+
+ for (int i = 0; i < nComponents; i++) {
+ MSG("%d DOFs for %s\n",
+ componentSpaces[i]->getAdmin()->getUsedSize(),
+ componentSpaces[i]->getName().c_str());
+
+ rhs->getDOFVector(i)->set(0.0);
+
+ for (int j = 0; j < nComponents; j++) {
+ // Only if this variable is true, the current matrix will be assembled.
+ bool assembleMatrix = true;
+ // The DOFMatrix which should be assembled (or not, if assembleMatrix
+ // will be set to false).
+ DOFMatrix *matrix = (*systemMatrix)[i][j];
+
+ if (matrix)
+ matrix->calculateNnz();
+
+ // If the matrix was assembled before and it is marked to be assembled
+ // only once, it will not be assembled.
+ if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
+ assembleMatrix = false;
+ } else if (matrix) {
+ matrix->getBaseMatrix().
+ change_dim(componentSpaces[i]->getAdmin()->getUsedSize(),
+ componentSpaces[j]->getAdmin()->getUsedSize());
+
+ set_to_zero(matrix->getBaseMatrix());
+ }
+
+ // If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
+ if (!matrix || !assembleMatrix)
+ assembleMatrix = false;
+
+ // If the matrix should not be assembled, the rhs vector has to be considered.
+ // This will be only done, if i == j. So, if both is not true, we can jump
+ // to the next matrix.
+ if (!assembleMatrix && i != j) {
+ if (matrix)
+ nnz += matrix->getBaseMatrix().nnz();
+
+ continue;
+ }
+
+ if (assembleMatrix && matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->initMatrix(matrix);
+
+ if (matrix && assembleMatrix)
+ matrix->startInsertion(matrix->getNnz());
+ }
+ }
+
+ // === TRAVERSE ===
+
+ Mesh *mesh = componentMeshes[0];
+ const FiniteElemSpace *feSpace = componentSpaces[0];
+ const BasisFunction *basisFcts = feSpace->getBasisFcts();
+ ElementMatrix elMat(basisFcts->getNumber(), basisFcts->getNumber());
+ ElementMatrix tmpElMat(elMat);
+ ElementVector elVec(basisFcts->getNumber());
+ ElementVector tmpElVec(elVec);
+ TraverseStack stack;
+ BoundaryType *bound =
+ useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
+ while (elInfo) {
+ if (useGetBound)
+ basisFcts->getBound(elInfo, bound);
+
+ for (std::map<Operator*, std::vector<OperatorPos> >::iterator opIt = operators.begin();
+ opIt != operators.end(); ++opIt) {
+ if (opIt->first->getNeedDualTraverse() == true)
+ continue;
+
+ if (opFlags[opIt->first].isSet(Operator::MATRIX_OPERATOR)) {
+ set_to_zero(elMat);
+ opIt->first->getElementMatrix(elInfo, elMat, 1.0);
+ }
+ if (opFlags[opIt->first].isSet(Operator::VECTOR_OPERATOR)) {
+ set_to_zero(elVec);
+ opIt->first->getElementVector(elInfo, elVec, 1.0);
+ }
+
+ for (std::vector<OperatorPos>::iterator posIt = opIt->second.begin();
+ posIt != opIt->second.end(); ++posIt) {
+
+ if (posIt->operatorType.isSet(Operator::MATRIX_OPERATOR)) {
+ if (!posIt->factor || *(posIt->factor) == 1.0) {
+ (*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(elMat, bound, elInfo, NULL);
+ } else {
+ tmpElMat = *(posIt->factor) * elMat;
+ (*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(tmpElMat, bound, elInfo, NULL);
+ }
+ }
+
+ if (posIt->operatorType.isSet(Operator::VECTOR_OPERATOR)) {
+ if (!posIt->factor || *(posIt->factor) == 1.0) {
+ rhs->getDOFVector(posIt->row)->addElementVector(1.0, elVec, bound, elInfo);
+ } else {
+ tmpElVec = *(posIt->factor) * elVec;
+ rhs->getDOFVector(posIt->row)->addElementVector(1.0, tmpElVec, bound, elInfo);
+ }
+ }
+ }
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ if (useGetBound)
+ delete [] bound;
+
+ // === FINELIZE ===
+
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++) {
+ bool assembleMatrix = true;
+ DOFMatrix *matrix = (*systemMatrix)[i][j];
+
+ if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j])
+ assembleMatrix = false;
+ if (!matrix || !assembleMatrix)
+ assembleMatrix = false;
+ if (!assembleMatrix && i != j)
+ continue;
+
+ assembledMatrix[i][j] = true;
+
+ if (assembleMatrix) {
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+ matrix->finishInsertion();
+ }
+ if (assembleMatrix && matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->exitMatrix(matrix);
+ if (matrix)
+ nnz += matrix->getBaseMatrix().nnz();
+ }
+
+ assembleBoundaryConditions(rhs->getDOFVector(i),
+ solution->getDOFVector(i),
+ componentMeshes[i],
+ assembleFlag);
+ }
+
+ solverMatrix.setMatrix(*systemMatrix);
+
+ createPrecon();
+
+ INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
+
+#ifdef _OPENMP
+ INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
+ TIME_USED(first, clock()), omp_get_wtime() - wtime);
+#else
+ INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
+ TIME_USED(first, clock()));
+#endif
+ }
bool ProblemVec::dualMeshTraverseRequired()
{
@@ -1192,6 +1384,10 @@ namespace AMDiS {
break;
}
}
+
+ OperatorPos opPos = {i, j, factor, estFactor, Operator::MATRIX_OPERATOR};
+ operators[op].push_back(opPos);
+ opFlags[op].setFlag(Operator::MATRIX_OPERATOR);
}
@@ -1221,6 +1417,10 @@ namespace AMDiS {
break;
}
}
+
+ OperatorPos opPos = {i, -1, factor, estFactor, Operator::VECTOR_OPERATOR};
+ operators[op].push_back(opPos);
+ opFlags[op].setFlag(Operator::VECTOR_OPERATOR);
}
@@ -1391,6 +1591,10 @@ namespace AMDiS {
Flag assembleFlag,
DOFMatrix *matrix, DOFVector<double> *vector)
{
+ FUNCNAME("ProblemVec::assembleOnOnMesh()");
+
+ // clock_t first = clock();
+
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basisFcts = feSpace->getBasisFcts();
@@ -1536,6 +1740,11 @@ namespace AMDiS {
delete [] bound;
} // pragma omp parallel
+
+// INFO(info, 8)("Assemble matrix with %d mat ops and %d vec ops needed %.5f seconds\n",
+// matrix ? matrix->getOperators().size() : 0,
+// vector ? vector->getOperators().size() : 0,
+// TIME_USED(first, clock()));
}
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index a4bf93ff8c27c79c702c1f6c3ce15a1c426415de..57b058f4724e08a49c1bb5b69f4c264c88f78310 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -39,6 +39,13 @@
namespace AMDiS {
+ struct OperatorPos
+ {
+ int row, col;
+ double *factor, *estFactor;
+ Flag operatorType;
+ };
+
class ProblemVec : public ProblemStatBase,
public StandardProblemIteration
{
@@ -175,6 +182,8 @@ namespace AMDiS {
bool assembleMatrix = true,
bool assembleVector = true);
+ void buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag);
+
bool dualMeshTraverseRequired();
void dualAssemble(AdaptInfo *adaptInfo, Flag flag,
@@ -644,6 +653,10 @@ namespace AMDiS {
/// If true, AMDiS prints information about the assembling procedure to the screen.
bool writeAsmInfo;
+
+ std::map<Operator*, std::vector<OperatorPos> > operators;
+
+ std::map<Operator*, Flag> opFlags;
};
}
diff --git a/AMDiS/src/QPsiPhi.h b/AMDiS/src/QPsiPhi.h
index c285a3d88a51cc79a799a855393321ab8bbe62ea..455af55a607025914bbf03bdc38abfe43576f46e 100644
--- a/AMDiS/src/QPsiPhi.h
+++ b/AMDiS/src/QPsiPhi.h
@@ -39,35 +39,29 @@ namespace AMDiS {
class compareQPsiPhi : public std::unary_function<bool, T*>
{
public:
- /** \brief
- * Constructor.
- */
+ /// Constructor.
compareQPsiPhi(const BasisFunction* psi_,
const BasisFunction* phi_,
const Quadrature* quad_)
- : psi(psi_), phi(phi_), quadrature(quad_)
+ : psi(psi_),
+ phi(phi_),
+ quadrature(quad_)
{}
- /** \brief
- * Returns true, if *q is equivalent to *this.
- */
- bool operator()(T* q) {
+ /// Returns true, if *q is equivalent to *this.
+ bool operator()(T* q)
+ {
return (q->psi == psi) && (q->phi == phi) && (q->quadrature == quadrature);
}
private:
- /** \brief
- * Basis functions of QPsiPhi.
- */
+ /// Basis functions of QPsiPhi.
const BasisFunction *psi, *phi;
- /** \brief
- * Quadrature of QPsiPhi.
- */
+ /// Quadrature of QPsiPhi.
const Quadrature *quadrature;
};
- // ============================================================================
/**
* \ingroup Integration
@@ -79,36 +73,27 @@ namespace AMDiS {
class compareQPsi : public std::unary_function<bool, T*>
{
public:
- /** \brief
- * Constructor.
- */
+ /// Constructor.
compareQPsi(const BasisFunction* psi_,
const Quadrature* quad_)
- : psi(psi_), quadrature(quad_)
+ : psi(psi_),
+ quadrature(quad_)
{}
- /** \brief
- * Returns true, if *q is equivalent to *this.
- */
- bool operator()(T* q) {
+ /// Returns true, if *q is equivalent to *this.
+ bool operator()(T* q)
+ {
return (q->psi == psi) && (q->quadrature == quadrature);
}
private:
- /** \brief
- * Basis functions of the QPsi object.
- */
+ /// Basis functions of the QPsi object.
const BasisFunction *psi;
- /** \brief
- * Quadrature of the QPsi object.
- */
+ /// Quadrature of the QPsi object.
const Quadrature *quadrature;
};
- // ===========================================================================
- // ===== class Q11PsiPhi =====================================================
- // ===========================================================================
/** \ingroup Integration
* \brief
@@ -119,130 +104,114 @@ namespace AMDiS {
class Q11PsiPhi
{
protected:
- /** \brief
- * Constructor
- */
+ /// Constructor
Q11PsiPhi(const BasisFunction *psi,
const BasisFunction *phi,
- const Quadrature *q);
-
+ const Quadrature *q);
public:
- /** \brief
- * Destructor
- */
+ /// Destructor
~Q11PsiPhi();
- /** \brief
- * Returns a Q11PsiPhi object.
- */
+ /// Returns a Q11PsiPhi object.
static const Q11PsiPhi* provideQ11PsiPhi(const BasisFunction *,
const BasisFunction *,
const Quadrature*);
- /** \brief
- * Compares two Q11PsiPhi objects.
- */
+ /// Compares two Q11PsiPhi objects.
const bool operator==(const Q11PsiPhi&) const;
- /** \brief
- * Compares two Q11PsiPhi objects.
- */
- const bool operator!=(const Q11PsiPhi& q11pp) const {
+ /// Compares two Q11PsiPhi objects.
+ const bool operator!=(const Q11PsiPhi& q11pp) const
+ {
return !(operator==(q11pp));
}
- /** \brief
- * Returns \ref values[i][j][k]
- */
+ /// Returns \ref values[i][j][k]
inline const double getValue(unsigned int i,
unsigned int j,
unsigned int v) const
{
- if ((values)&&(values[i])&&(values[i][j])&&((static_cast<int>(v))<nrEntries[i][j]))
+ if (values && values[i] && values[i][j] &&
+ (static_cast<int>(v) < nrEntries[i][j]))
return values[i][j][v];
- return 0.;
+
+ return 0.0;
}
- /** \brief
- * Returns \ref nrEntries[i][j]
- */
+ /// Returns \ref nrEntries[i][j]
inline const int getNumberEntries(unsigned int i, unsigned int j) const
{
- if ((nrEntries)&&(nrEntries[i])) return nrEntries[i][j];
+ if (nrEntries && nrEntries[i])
+ return nrEntries[i][j];
+
return 0;
}
- /** \brief
- * Returns \ref nrEntries
- */
- inline const int** getNumberEntries() const {
- return const_cast<const int**>( nrEntries);
+ /// Returns \ref nrEntries
+ inline const int** getNumberEntries() const
+ {
+ return const_cast<const int**>(nrEntries);
}
- /** \brief
- * Returns \ref k[i1][i2][i3]
- */
- inline const int
- getK(unsigned int i1, unsigned int i2, unsigned int i3) const
+ /// Returns \ref k[i1][i2][i3]
+ inline const int getK(unsigned int i1,
+ unsigned int i2,
+ unsigned int i3) const
{
- if ((k)&&(k[i1])&&(k[i1][i2])&&((static_cast<int>(i3))<nrEntries[i1][i2]))
+ if (k && k[i1] && k[i1][i2] && (static_cast<int>(i3) < nrEntries[i1][i2]))
return k[i1][i2][i3];
+
return 0;
- };
+ }
- /** \brief
- * Returns \ref l[i][j][v]
- */
+ /// Returns \ref l[i][j][v]
inline const int getL(unsigned int i, unsigned int j, unsigned int v) const
{
- if ((l)&&(l[i])&&(l[i][j])&&((static_cast<int>(v))<nrEntries[i][j])) return l[i][j][v];
+ if (l && l[i] && l[i][j] && (static_cast<int>(v) < nrEntries[i][j]))
+ return l[i][j][v];
+
return 0;
}
- /** \brief
- * Returns \values[i][j]
- */
- inline const double *getValVec(unsigned int i, unsigned int j) const {
- if ((values)&&(values[i])&&(values[i][j])) return values[i][j];
+ /// Returns \values[i][j]
+ inline const double *getValVec(unsigned int i, unsigned int j) const
+ {
+ if (values && values[i] && values[i][j])
+ return values[i][j];
+
return NULL;
}
- /** \brief
- * Returns \ref k[i][j]
- */
- inline const int *getKVec(unsigned int i, unsigned int j) const {
- if ((k)&&(k[i])&&(k[i][j])) return k[i][j];
+ /// Returns \ref k[i][j]
+ inline const int *getKVec(unsigned int i, unsigned int j) const
+ {
+ if (k && k[i] && k[i][j])
+ return k[i][j];
+
return NULL;
}
- /** \brief
- * Returns \ref l[i][j]
- */
- inline const int *getLVec(unsigned int i,unsigned int j) const {
- if ((l)&&(l[i])&&(l[i][j])) return l[i][j];
+ /// Returns \ref l[i][j]
+ inline const int *getLVec(unsigned int i, unsigned int j) const
+ {
+ if (l && l[i] && l[i][j])
+ return l[i][j];
+
return NULL;
}
protected:
- /** \brief
- * List of pointers to all Q11PsiPhi objects
- */
+ /// List of pointers to all Q11PsiPhi objects.
static std::list<Q11PsiPhi*> preList;
- /** \brief
- * Pointer to the first set of basis functions
- */
+ /// Pointer to the first set of basis functions.
const BasisFunction *psi;
- /** \brief
- * pointer to the second set of basis functions
- */
+ /// Pointer to the second set of basis functions.
const BasisFunction *phi;
- /** \brief
- * Pointer to the Quadrature which is used for the integration
- */
+ /// Pointer to the Quadrature which is used for the integration.
const Quadrature *quadrature;
/** \brief
@@ -261,45 +230,28 @@ namespace AMDiS {
*/
double ***values;
- /** \brief
- * tensor storing the indices k of the non zero integrals;
- */
+ /// Tensor storing the indices k of the non zero integrals.
int ***k;
- /** \brief
- * tensor storing the indices l of the non zero integrals;
- */
+ /// Tensor storing the indices l of the non zero integrals.
int ***l;
- /** \brief
- * Number of non zero entries.
- */
+ /// Number of non zero entries.
int allEntries;
- /** \brief
- * Pointer to an array in values.
- */
+ /// Pointer to an array in values.
double *val_vec;
- /** \brief
- * Pointer to an array in k.
- */
+ /// Pointer to an array in k.
int *k_vec;
- /** \brief
- * Pointer to an array in l.
- */
+ /// Pointer to an array in l.
int *l_vec;
friend class compareQPsiPhi<Q11PsiPhi>;
};
-
- // ===========================================================================
- // ===== class Q10PsiPhi =====================================================
- // ===========================================================================
-
/** \ingroup Integration
* \brief
* Calculates element stiffness matrices by preevaluated integrals over the
@@ -309,111 +261,96 @@ namespace AMDiS {
class Q10PsiPhi
{
protected:
- /** \brief
- * Constructor
- */
+ /// Constructor
Q10PsiPhi(const BasisFunction *psi,
const BasisFunction *phi,
- const Quadrature *q);
+ const Quadrature *q);
public:
- /** \brief
- * Destructor
- */
+ /// Destructor
~Q10PsiPhi();
- /** \brief
- * Returns a Q10PsiPhi object.
- */
+ /// Returns a Q10PsiPhi object.
static const Q10PsiPhi* provideQ10PsiPhi(const BasisFunction *,
const BasisFunction *,
const Quadrature*);
- /** \brief
- * Compares two Q10PsiPhi objects.
- */
+ /// Compares two Q10PsiPhi objects.
const bool operator==(const Q10PsiPhi&) const;
- /** \brief
- * Compares two Q10PsiPhi objects.
- */
- const bool operator!=(const Q10PsiPhi& q10pp) const {
+ /// Compares two Q10PsiPhi objects.
+ const bool operator!=(const Q10PsiPhi& q10pp) const
+ {
return !(operator==(q10pp));
}
-
- /** \brief
- * Returns \ref values[i][j][k]
- */
+ /// Returns \ref values[i][j][k]
inline const double getValue(unsigned int i,
unsigned int j,
unsigned int v) const
{
- if ((values)&&(values[i])&&(values[i][j])&&((static_cast<int>(v))<nrEntries[i][j])) return values[i][j][v];
- return 0.;
+ if (values && values[i] && values[i][j] &&
+ (static_cast<int>(v) < nrEntries[i][j]))
+ return values[i][j][v];
+
+ return 0.0;
}
- /** \brief
- * Returns \ref nrEntries[i][j]
- */
- inline const int getNumberEntries(unsigned int i, unsigned int j) const {
- if ((nrEntries)&&(nrEntries[i])) return nrEntries[i][j];
+ /// Returns \ref nrEntries[i][j]
+ inline const int getNumberEntries(unsigned int i, unsigned int j) const
+ {
+ if (nrEntries && nrEntries[i])
+ return nrEntries[i][j];
+
return 0;
}
- /** \brief
- * Returns \ref nrEntries
- */
- inline const int** getNumberEntries() const {
+ /// Returns \ref nrEntries
+ inline const int** getNumberEntries() const
+ {
return const_cast<const int**>(nrEntries);
}
- /** \brief
- * Returns \ref k[i1][i2][i3]
- */
- inline const int
- getK(unsigned int i1, unsigned int i2, unsigned int i3) const
+ /// Returns \ref k[i1][i2][i3]
+ inline const int getK(unsigned int i1,
+ unsigned int i2,
+ unsigned int i3) const
{
- if ((k)&&(k[i1])&&(k[i1][i2])&&((static_cast<int>(i3))<nrEntries[i1][i2]))
+ if (k && k[i1] && k[i1][i2] && (static_cast<int>(i3) < nrEntries[i1][i2]))
return k[i1][i2][i3];
+
return 0;
}
- /** \brief
- * Returns \values[i][j]
- */
- inline const double *getValVec(unsigned int i, unsigned int j) const {
- if ((values)&&(values[i])&&(values[i][j])) return values[i][j];
+ /// Returns \values[i][j]
+ inline const double *getValVec(unsigned int i, unsigned int j) const
+ {
+ if (values && values[i] && values[i][j])
+ return values[i][j];
+
return NULL;
}
- /** \brief
- * Returns \ref k[i][j]
- */
- inline const int *getKVec(unsigned int i, unsigned int j) const {
- if ((k)&&(k[i])&&(k[i][j])) return k[i][j];
+ /// Returns \ref k[i][j]
+ inline const int *getKVec(unsigned int i, unsigned int j) const
+ {
+ if (k && k[i] && k[i][j])
+ return k[i][j];
+
return NULL;
}
protected:
- /** \brief
- * List of pointers to all Q11PsiPhi objects
- */
+ /// List of pointers to all Q11PsiPhi objects
static std::list<Q10PsiPhi*> preList;
- /** \brief
- * Pointer to the first set of basis functions
- */
+ /// Pointer to the first set of basis functions
const BasisFunction *psi;
- /** \brief
- * pointer to the second set of basis functions
- */
+ /// Pointer to the second set of basis functions
const BasisFunction *phi;
- /** \brief
- * Pointer to the Quadrature which is used for the integration
- */
+ /// Pointer to the Quadrature which is used for the integration
const Quadrature *quadrature;
/** \brief
@@ -432,32 +369,21 @@ namespace AMDiS {
*/
double ***values;
- /** \brief
- * tensor storing the indices k of the non zero integrals;
- */
+ /// Tensor storing the indices k of the non zero integrals;
int ***k;
- /** \brief
- * Number of all non zero entries.
- */
+ /// Number of all non zero entries.
int allEntries;
- /** \brief
- * Pointer to an array in values.
- */
+ /// Pointer to an array in values.
double *val_vec;
- /** \brief
- * Pointer to an array in k.
- */
+ /// Pointer to an array in k.
int *k_vec;
friend class compareQPsiPhi<Q10PsiPhi>;
};
- // ===========================================================================
- // ===== class Q01PsiPhi =====================================================
- // ===========================================================================
/** \ingroup Integration
* \brief
@@ -468,119 +394,101 @@ namespace AMDiS {
class Q01PsiPhi
{
protected:
- /** \brief
- * Constructor
- */
+ /// Constructor
Q01PsiPhi(const BasisFunction *psi,
const BasisFunction *phi,
- const Quadrature *q);
+ const Quadrature *q);
public:
- /** \brief
- * Destructor
- */
+ /// Destructor
~Q01PsiPhi();
- /** \brief
- * Returns a Q01PsiPhi object.
- */
+ /// Returns a Q01PsiPhi object.
static const Q01PsiPhi* provideQ01PsiPhi(const BasisFunction *,
const BasisFunction *,
const Quadrature*);
- /** \brief
- * Compares two Q01PsiPhi objects.
- */
+ /// Compares two Q01PsiPhi objects.
const bool operator==(const Q01PsiPhi&) const;
- /** \brief
- * Compares two Q01PsiPhi objects.
- */
- const bool operator!=(const Q01PsiPhi& q01pp) const {
+ /// Compares two Q01PsiPhi objects.
+ const bool operator!=(const Q01PsiPhi& q01pp) const
+ {
return !(operator==(q01pp));
}
-
- /** \brief
- * Returns \ref values[i][j][k]
- */
+ /// Returns \ref values[i][j][k]
inline const double getValue(unsigned int i,
unsigned int j,
unsigned int v) const
{
- if ((values)&&(values[i])&&(values[i][j])&&((static_cast<int>(v))<nrEntries[i][j]))
+ if (values && values[i] && values[i][j] &&
+ (static_cast<int>(v) < nrEntries[i][j]))
return values[i][j][v];
- return 0.;
+
+ return 0.0;
}
- /** \brief
- * Returns \ref nrEntries[i][j]
- */
+ /// Returns \ref nrEntries[i][j]
inline const int getNumberEntries(unsigned int i, unsigned int j) const
{
- if ((nrEntries)&&(nrEntries[i])) return nrEntries[i][j];
+ if (nrEntries && nrEntries[i])
+ return nrEntries[i][j];
+
return 0;
}
- /** \brief
- * Returns \ref nrEntries
- */
- inline const int** getNumberEntries() const {
+ /// Returns \ref nrEntries
+ inline const int** getNumberEntries() const
+ {
return const_cast<const int**>(nrEntries);
}
- /** \brief
- * Returns \ref k[i1][i2][i3]
- */
- inline const int
- getK(unsigned int i1, unsigned int i2, unsigned int i3) const;
+ /// Returns \ref k[i1][i2][i3]
+ inline const int getK(unsigned int i1,
+ unsigned int i2,
+ unsigned int i3) const;
+
+ /// Returns \values[i][j]
+ inline const double *getValVec(unsigned int i, unsigned int j) const
+ {
+ if (values && values[i] && values[i][j])
+ return values[i][j];
- /** \brief
- * Returns \values[i][j]
- */
- inline const double *getValVec(unsigned int i, unsigned int j) const {
- if ((values)&&(values[i])&&(values[i][j])) return values[i][j];
return NULL;
}
- /** \brief
- * Returns \ref k[i][j]
- */
- inline const int *getLVec(unsigned int i, unsigned int j) const {
- if ((l)&&(l[i])&&(l[i][j])) return l[i][j];
+ /// Returns \ref k[i][j]
+ inline const int *getLVec(unsigned int i, unsigned int j) const
+ {
+ if (l && l[i] && l[i][j])
+ return l[i][j];
+
return NULL;
}
- /** \brief
- * Returns \ref k[i][j][v]
- */
+ /// Returns \ref k[i][j][v]
inline const int getL(unsigned int i,
unsigned int j,
unsigned int v) const
{
- if ((l)&&(l[i])&&(l[i][j])&&((static_cast<int>(v))<nrEntries[i][j])) return l[i][j][v];
+ if (l && l[i] && l[i][j] && (static_cast<int>(v) < nrEntries[i][j]))
+ return l[i][j][v];
+
return 0;
}
protected:
- /** \brief
- * List of pointers to all Q11PsiPhi objects
- */
+ /// List of pointers to all Q11PsiPhi objects
static std::list<Q01PsiPhi*> preList;
- /** \brief
- * Pointer to the first set of basis functions
- */
+ /// Pointer to the first set of basis functions
const BasisFunction *psi;
- /** \brief
- * pointer to the second set of basis functions
- */
+ /// Pointer to the second set of basis functions
const BasisFunction *phi;
- /** \brief
- * Pointer to the Quadrature which is used for the integration
- */
+ /// Pointer to the Quadrature which is used for the integration
const Quadrature *quadrature;
/** \brief
@@ -599,32 +507,21 @@ namespace AMDiS {
*/
double ***values;
- /** \brief
- * tensor storing the indices l of the non zero integrals;
- */
+ /// Tensor storing the indices l of the non zero integrals;
int ***l;
- /** \brief
- * Number of all non zero entries.
- */
+ /// Number of all non zero entries.
int allEntries;
- /** \brief
- * Pointer to an array in values
- */
+ /// Pointer to an array in values
double *val_vec;
- /** \brief
- * Pointer to an array in l.
- */
+ /// Pointer to an array in l.
int *l_vec;
friend class compareQPsiPhi<Q01PsiPhi>;
};
- // ===========================================================================
- // ===== class Q00PsiPhi =====================================================
- // ===========================================================================
/** \ingroup Integration
* \brief
@@ -633,24 +530,16 @@ namespace AMDiS {
*/
class Q00PsiPhi
{
- /** \brief
- * List of pointers to all Q11PsiPhi objects
- */
+ /// List of pointers to all Q11PsiPhi objects
static std::list<Q00PsiPhi*> preList;
- /** \brief
- * Pointer to the first set of basis functions
- */
+ /// Pointer to the first set of basis functions
const BasisFunction *psi;
- /** \brief
- * pointer to the second set of basis functions
- */
+ /// Pointer to the second set of basis functions
const BasisFunction *phi;
- /** \brief
- * Pointer to the Quadrature which is used for the integration
- */
+ /// Pointer to the Quadrature which is used for the integration
const Quadrature *quadrature;
/** \brief
@@ -660,60 +549,40 @@ namespace AMDiS {
* for the pair (psi[i], phi[j]), 0 <= i <= psi->getNumber(),
* 0 <= j <= phi->getNumber()
*/
- double **values;
+ double **values;
protected:
- /** \brief
- * Constructor
- */
- Q00PsiPhi(const BasisFunction *,
- const BasisFunction *,
- const Quadrature*);
+ /// Constructor
+ Q00PsiPhi(const BasisFunction *, const BasisFunction *, const Quadrature*);
public:
- /** \brief
- * Destructor
- */
+ /// Destructor
~Q00PsiPhi();
- /** \brief
- * Returns a Q00PsiPhi object.
- */
+ /// Returns a Q00PsiPhi object.
static Q00PsiPhi* provideQ00PsiPhi(const BasisFunction *,
const BasisFunction *,
const Quadrature*);
-
- /** \brief
- * Compares two Q00PsiPhi objects.
- */
+ /// Compares two Q00PsiPhi objects.
const bool operator==(const Q00PsiPhi&) const;
- /** \brief
- * Compares two Q00PsiPhi objects.
- */
- const bool operator!=(const Q00PsiPhi& q00pp) const {
+ /// Compares two Q00PsiPhi objects.
+ const bool operator!=(const Q00PsiPhi& q00pp) const
+ {
return !(operator==(q00pp));
}
- /** \brief
- * Returns \ref values[i][j]
- */
- const double getValue(unsigned int i,unsigned int j) const;
+ /// Returns \ref values[i][j]
+ const double getValue(unsigned int i, unsigned int j) const;
- /** \brief
- * Returns \ref values[i]
- */
+ /// Returns \ref values[i]
const double *getValVec(unsigned int i) const;
friend class compareQPsiPhi<Q00PsiPhi>;
};
- // ===========================================================================
- // ===== class Q0Psi =========================================================
- // ===========================================================================
-
/** \ingroup Integration
* \brief
* Integral of psi.
@@ -721,76 +590,53 @@ namespace AMDiS {
class Q0Psi
{
protected:
- /** \brief
- * List of pointers to all Q0Psi objects
- */
+ /// List of pointers to all Q0Psi objects
static std::list<Q0Psi*> preList;
- /** \brief
- * Pointer to the first set of basis functions
- */
+ /// Pointer to the first set of basis functions
const BasisFunction *psi;
- /** \brief
- * Pointer to the Quadrature which is used for the integration
- */
+ /// Pointer to the Quadrature which is used for the integration
const Quadrature *quadrature;
- /** \brief
- * Vector storing the integrals
- */
+ /// Vector storing the integrals
double* values;
protected:
- /** \brief
- * Constructor
- */
- Q0Psi(const BasisFunction *,
- const Quadrature*);
+ /// Constructor
+ Q0Psi(const BasisFunction *, const Quadrature*);
public:
- /** \brief
- * Destructor
- */
+ /// Destructor
~Q0Psi();
- /** \brief
- * Returns a Q0Psi object.
- */
+ /// Returns a Q0Psi object.
static Q0Psi* provideQ0Psi(const BasisFunction *, const Quadrature*);
- /** \brief
- * Compares two Q0Psi objects.
- */
+ /// Compares two Q0Psi objects.
const bool operator==(const Q0Psi&) const;
- /** \brief
- * Compares two Q0Psi objects.
- */
- const bool operator!=(const Q0Psi& q0p) const {
+ /// Compares two Q0Psi objects.
+ const bool operator!=(const Q0Psi& q0p) const
+ {
return !(operator==(q0p));
}
- /** \brief
- * Returns \ref value
- */
- inline const double getValue(int i) const {
+ /// Returns \ref value
+ inline const double getValue(int i) const
+ {
return values[i];
}
- /** \brief
- * Returns \ref values[i]
- */
- inline const double *getValVec() const {
+ /// Returns \ref values[i]
+ inline const double *getValVec() const
+ {
return values;
}
friend class compareQPsi<Q0Psi>;
};
- // ===========================================================================
- // ===== class Q1Psi =========================================================
- // ===========================================================================
/** \ingroup Integration
* \brief
@@ -799,105 +645,85 @@ namespace AMDiS {
class Q1Psi
{
protected:
- /** \brief
- * Constructor
- */
- Q1Psi(const BasisFunction *psi,
- const Quadrature *q);
+ /// Constructor
+ Q1Psi(const BasisFunction *psi, const Quadrature *q);
public:
- /** \brief
- * Destructor
- */
+ /// Destructor
~Q1Psi();
- /** \brief
- * Returns a Q1Psi object.
- */
+ /// Returns a Q1Psi object.
static const Q1Psi* provideQ1Psi(const BasisFunction *,
const Quadrature*);
- /** \brief
- * Compares two Q1Psi objects.
- */
+ /// Compares two Q1Psi objects.
const bool operator==(const Q1Psi&) const;
- /** \brief
- * Compares two Q1Psi objects.
- */
- const bool operator!=(const Q1Psi& q1p) const {
+ /// Compares two Q1Psi objects.
+ const bool operator!=(const Q1Psi& q1p) const
+ {
return !(operator==(q1p));
}
- /** \brief
- * Returns \ref values[i][j]
- */
+ /// Returns \ref values[i][j]
inline const double getValue(unsigned int i,
unsigned int j) const
{
- if ((values)&&(values[i])&&((static_cast<int>(j))<nrEntries[i]))
+ if (values && values[i] && (static_cast<int>(j) < nrEntries[i]))
return values[i][j];
- return 0.;
+
+ return 0.0;
}
- /** \brief
- * Returns \ref nrEntries[i]
- */
+ /// Returns \ref nrEntries[i]
inline const int getNumberEntries(unsigned int i) const
{
- if (nrEntries) return nrEntries[i];
- return 0;
+ return (nrEntries ? nrEntries[i] : 0);
}
- /** \brief
- * Returns \ref nrEntries
- */
- inline const int* getNumberEntries() const {
+ /// Returns \ref nrEntries
+ inline const int* getNumberEntries() const
+ {
return const_cast<const int*>(nrEntries);
}
- /** \brief
- * Returns \ref k[i1][i2]
- */
+ /// Returns \ref k[i1][i2]
inline const int
getK(unsigned int i1, unsigned int i2) const
{
- if ((k)&&(k[i1])&&((static_cast<int>(i2))<nrEntries[i1]))
+ if (k && k[i1] && (static_cast<int>(i2) < nrEntries[i1]))
return k[i1][i2];
+
return 0;
- };
+ }
+
+ /// Returns \ref k[i]
+ inline const int *getKVec(unsigned int i) const
+ {
+ if (k && k[i])
+ return k[i];
- /** \brief
- * Returns \ref k[i]
- */
- inline const int *getKVec(unsigned int i) const {
- if ((k)&&(k[i])) return k[i];
return NULL;
}
- /** \brief
- * Returns \values[i]
- */
- inline const double *getValVec(unsigned int i) const {
- if ((values)&&(values[i])) return values[i];
+ /// Returns \values[i]
+ inline const double *getValVec(unsigned int i) const
+ {
+ if (values && values[i])
+ return values[i];
+
return NULL;
}
protected:
- /** \brief
- * List of pointers to all Q1Psi objects
- */
+ /// List of pointers to all Q1Psi objects
static std::list<Q1Psi*> preList;
- /** \brief
- * Pointer to the first set of basis functions
- */
+ /// Pointer to the first set of basis functions
const BasisFunction *psi;
- /** \brief
- * Pointer to the Quadrature which is used for the integration
- */
+ /// Pointer to the Quadrature which is used for the integration
const Quadrature *quadrature;
/** \brief
@@ -908,30 +734,20 @@ namespace AMDiS {
*/
int *nrEntries;
- /** \brief
- * Number of all non zero entries.
- */
+ /// Number of all non zero entries.
int allEntries;
- /** \brief
- * tensor storing the non zero integrals; values[i] is a vector of length
- * \ref nrEntries[i] storing the non zero values for psi[i].
- */
+ /// Tensor storing the non zero integrals; values[i] is a vector of length
+ /// \ref nrEntries[i] storing the non zero values for psi[i].
double **values;
- /** \brief
- * Matrix storing the indices k of the non zero integrals;
- */
+ /// Matrix storing the indices k of the non zero integrals;
int **k;
- /** \brief
- * Pointer to an array in values.
- */
+ /// Pointer to an array in values.
double *val_vec;
- /** \brief
- * Pointer to an array in k.
- */
+ /// Pointer to an array in k.
int *k_vec;
friend class compareQPsi<Q1Psi>;
diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc
index 19e9a27cfdddb6f602cc9bf830c7bc3b3ef0efae..73f35c21fb5311c3496148dfbfdb3d3e09649abf 100644
--- a/AMDiS/src/Quadrature.cc
+++ b/AMDiS/src/Quadrature.cc
@@ -19,7 +19,7 @@ namespace AMDiS {
const int Quadrature::maxNQuadPoints[4] = {0, 10, 61, 64};
- std::list<FastQuadrature*> FastQuadrature::fastQuadList;
+ list<FastQuadrature*> FastQuadrature::fastQuadList;
int FastQuadrature::max_points = 0;
Quadrature::Quadrature(const Quadrature& q)
@@ -1490,7 +1490,7 @@ namespace AMDiS {
#pragma omp critical (provideFastQuad)
#endif
{
- std::list<FastQuadrature*>::iterator fast = fastQuadList.begin();
+ list<FastQuadrature*>::iterator fast = fastQuadList.begin();
for (fast = fastQuadList.begin(); fast != fastQuadList.end(); fast++)
@@ -1530,12 +1530,8 @@ namespace AMDiS {
// ----- initialize phi ---------------------------------------------
- if (!phi && init_flag.isSet(INIT_PHI)) { // check flag
-
- // allocate memory
- phi = new double*[nPoints];
- for (int i = 0; i < nPoints; i++)
- phi[i] = new double[nBasFcts];
+ if (num_rows(phi) == 0 && init_flag.isSet(INIT_PHI)) {
+ phi.change_dim(nPoints, nBasFcts);
// fill memory
for (int i = 0; i< nPoints; i++) {
@@ -1549,17 +1545,20 @@ namespace AMDiS {
}
// initialize grd_phi
-
- if (!grdPhi && init_flag.isSet(INIT_GRD_PHI)) { // check flag
+ if (grdPhi.empty() && init_flag.isSet(INIT_GRD_PHI)) {
// allocate memory
- grdPhi = new MatrixOfFixVecs<DimVec<double> >(dim, nPoints, nBasFcts, NO_INIT);
+ grdPhi.resize(nPoints);
// fill memory
for (int i = 0; i< nPoints; i++) {
+ grdPhi[i].resize(nBasFcts);
lambda = quadrature->getLambda(i);
- for (int j = 0; j < nBasFcts; j++)
- (*(basisFunctions->getGrdPhi(j)))(lambda, (*(grdPhi))[i][j]);
+
+ for (int j = 0; j < nBasFcts; j++) {
+ grdPhi[i][j].change_dim(dim + 1);
+ (*(basisFunctions->getGrdPhi(j)))(lambda, grdPhi[i][j]);
+ }
}
// update flag
@@ -1568,7 +1567,7 @@ namespace AMDiS {
// initialize D2_phi
- if (!D2Phi && init_flag.isSet(INIT_D2_PHI)) { // check flag
+ if (!D2Phi && init_flag.isSet(INIT_D2_PHI)) {
// allocate memory
D2Phi = new MatrixOfFixVecs<DimMat<double> >(dim, nPoints, nBasFcts, NO_INIT);
@@ -1604,20 +1603,18 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
int nBasFcts = basisFunctions->getNumber();
- if (fastQuad.phi) {
- phi = new double*[nPoints];
- for (int i = 0; i < nPoints; i++) {
- phi[i] = new double[nBasFcts];
- for (int j = 0; j < nBasFcts; j++)
- phi[i][j] = fastQuad.phi[i][j];
- }
+ if (num_rows(fastQuad.phi) > 0) {
+ phi.change_dim(nPoints, nBasFcts);
+ phi = fastQuad.phi;
}
- if (fastQuad.grdPhi) {
- grdPhi = new MatrixOfFixVecs<DimVec<double> >(dim, nPoints, nBasFcts, NO_INIT);
- for (int i = 0; i < nPoints; i++)
+ if (!fastQuad.grdPhi.empty()) {
+ grdPhi.resize(nPoints);
+ for (int i = 0; i < nPoints; i++) {
+ grdPhi[i].resize(nBasFcts);
for (int j = 0; j < nBasFcts; j++)
- (*grdPhi)[i][j] = (*(fastQuad.grdPhi))[i][j];
+ grdPhi[i][j] = fastQuad.grdPhi[i][j];
+ }
}
if (fastQuad.D2Phi) {
@@ -1631,13 +1628,6 @@ namespace AMDiS {
FastQuadrature::~FastQuadrature()
{
- int nPoints = quadrature->getNumPoints();
-
- for (int i = 0; i < nPoints; i++)
- delete [] phi[i];
-
- delete [] phi;
- delete grdPhi;
delete D2Phi;
}
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index fd3fc7c1d37e547660f5b1856bde86454f81aee6..e2a254ca39143a4fd5539ec5f1f8ffb39b218feb 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -24,6 +24,8 @@
#define AMDIS_QUADRATURE_H
#include <list>
+#include <vector>
+#include <boost/numeric/mtl/mtl.hpp>
#include "AMDiS_fwd.h"
#include "BasisFunction.h"
#include "Flag.h"
@@ -31,6 +33,8 @@
namespace AMDiS {
+ using namespace std;
+
/**
* \ingroup Assembler
*
@@ -93,7 +97,7 @@ namespace AMDiS {
double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f);
/// Returns \ref name
- inline std::string getName()
+ inline string getName()
{
return name;
}
@@ -184,7 +188,7 @@ namespace AMDiS {
protected:
/// Name of this Quadrature
- std::string name;
+ string name;
/// Quadrature is exact of this degree
int degree;
@@ -326,9 +330,8 @@ namespace AMDiS {
* flag.
*/
FastQuadrature(BasisFunction* basFcts, Quadrature* quad, Flag flag)
- : init_flag(flag),
- phi(NULL),
- grdPhi(NULL),
+ : init_flag(flag),
+ phi(0, 0),
D2Phi(NULL),
quadrature(quad),
basisFunctions(basFcts)
@@ -355,10 +358,10 @@ namespace AMDiS {
inline bool initialized(Flag flag)
{
if (flag == INIT_PHI)
- return (phi != NULL);
+ return (num_rows(phi) > 0);
if (flag == INIT_GRD_PHI)
- return (grdPhi != NULL);
+ return (!grdPhi.empty());
if (flag == INIT_D2_PHI)
return (D2Phi != NULL);
@@ -388,30 +391,29 @@ namespace AMDiS {
/// Returns (*\ref grdPhi)[q][i][j]
inline const double getGradient(int q, int i ,int j) const
{
- return (grdPhi) ? (*grdPhi)[q][i][j] : 0.0;
+ return (!grdPhi.empty()) ? grdPhi[q][i][j] : 0.0;
}
/// Returns (*\ref grdPhi)[q]
- inline VectorOfFixVecs<DimVec<double> >* getGradient(int q) const
+ inline const vector<mtl::dense_vector<double> >& getGradient(int q) const
{
- return (grdPhi) ? &((*grdPhi)[q]) : NULL;
+ return grdPhi[q];
}
- inline DimVec<double>& getGradient(int q, int i) const
+ inline const mtl::dense_vector<double>& getGradient(int q, int i) const
{
- return (*grdPhi)[q][i];
+ return grdPhi[q][i];
}
- /// Returns \ref phi[q][i]
- inline const double getPhi(int q, int i) const
+ inline const mtl::dense2D<double>& getPhi() const
{
- return (phi) ? phi[q][i] : 0.0;
+ return phi;
}
- /// Returns \ref phi[q]
- inline const double *getPhi(int q) const
+ /// Returns \ref phi[q][i]
+ inline const double getPhi(int q, int i) const
{
- return (phi) ? phi[q] : NULL;
+ return phi[q][i];
}
/// Returns \ref quadrature ->integrateStdSimplex(f)
@@ -491,7 +493,7 @@ namespace AMDiS {
* (quadrature->lambda[i]), 0 <= j < basisFunctions->getNumber() and
* 0 <= i < n_points
*/
- double **phi;
+ mtl::dense2D<double> phi;
/** \brief
* Matrix storing all gradients (with respect to the barycentric coordinates)
@@ -500,7 +502,7 @@ namespace AMDiS {
* for 0 <= j < basisFunctions->getNumber(),
* 0 <= i < . . . , n_points, and 0 <= k < DIM
*/
- MatrixOfFixVecs<DimVec<double> > *grdPhi;
+ vector<vector<mtl::dense_vector<double> > > grdPhi;
/** \brief
* Matrix storing all second derivatives (with respect to the barycentric
@@ -514,7 +516,7 @@ namespace AMDiS {
/** \brief
* List of all used FastQuadratures
*/
- static std::list<FastQuadrature*> fastQuadList;
+ static list<FastQuadrature*> fastQuadList;
/** \brief
* Maximal number of quadrature points for all yet initialised FastQuadrature
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index 2b9e8fbc1afe331c747171777bc3626eb08b1ab5..01ac52f9e1affc0392fd86075b986ab7f621d66a 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -225,21 +225,21 @@ namespace AMDiS {
continue;
if (dualElInfo)
- (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo,
- dualElInfo->largeElInfo,
- quad);
+ (*it)->getAssembler()->initElement(dualElInfo->smallElInfo,
+ dualElInfo->largeElInfo,
+ quad);
else
- (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
+ (*it)->getAssembler()->initElement(elInfo, NULL, quad);
}
if (C0 > 0.0)
for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it) {
if (dualElInfo)
- (*it)->getAssembler(omp_get_thread_num())->initElement(dualElInfo->smallElInfo,
- dualElInfo->largeElInfo,
- quad);
+ (*it)->getAssembler()->initElement(dualElInfo->smallElInfo,
+ dualElInfo->largeElInfo,
+ quad);
else
- (*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
+ (*it)->getAssembler()->initElement(elInfo, NULL, quad);
}
}
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index 1935bf8035a6b09fe537d67c904b9cfa56b142f6..c03c216515fb440e5cb3e14bdcfdc8f0f998b7cf 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -15,7 +15,6 @@
#include "Assembler.h"
#include "DOFVector.h"
#include "DOFMatrix.h"
-#include "OpenMP.h"
#include "SurfaceOperator.h"
#include <math.h>
@@ -225,9 +224,9 @@ namespace AMDiS {
neumannQuad = true;
} else {
if (neumannOperators) {
- if ((*neumannOperators)[0]->getAssembler(omp_get_thread_num())->
+ if ((*neumannOperators)[0]->getAssembler()->
getZeroOrderAssembler()->getQuadrature()->getNumPoints() >
- (*robinOperators)[0]->getAssembler(omp_get_thread_num())->
+ (*robinOperators)[0]->getAssembler()->
getZeroOrderAssembler()->getQuadrature()->getNumPoints())
neumannQuad = true;
}
@@ -235,20 +234,20 @@ namespace AMDiS {
std::vector<Operator*>::iterator op;
for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
- (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
+ (*op)->getAssembler()->initElement(elInfo);
for (int face = 0; face < dim + 1; face++) {
elInfo->getNormal(face, normal);
Quadrature *quadrature =
neumannQuad ?
- (*neumannOperators)[face]->getAssembler(omp_get_thread_num())->
+ (*neumannOperators)[face]->getAssembler()->
getZeroOrderAssembler()->getQuadrature() :
- (*robinOperators)[face]->getAssembler(omp_get_thread_num())->
+ (*robinOperators)[face]->getAssembler()->
getZeroOrderAssembler()->getQuadrature();
if (elInfo->getBoundary(face) == boundaryType) {
- (*neumannOperators)[face]->getAssembler(omp_get_thread_num())->
+ (*neumannOperators)[face]->getAssembler()->
getZeroOrderAssembler()->initElement(elInfo);
int nPoints = quadrature->getNumPoints();
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 683fc2e528d63a5a533fe9d29218585e7c113d70..08a29b9d4bb095c3273ab2aa3e20db793e75dc87 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -39,10 +39,8 @@ namespace AMDiS {
Quadrature *quad,
bool optimized)
{
- int myRank = omp_get_thread_num();
-
// check if a assembler is needed at all
- if (op->secondOrder[myRank].size() == 0)
+ if (op->secondOrder.size() == 0)
return NULL;
SecondOrderAssembler *newAssembler;
@@ -52,7 +50,7 @@ namespace AMDiS {
&optimizedSubAssemblers :
&standardSubAssemblers;
- std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
+ std::vector<OperatorTerm*> opTerms = op->zeroOrder;
sort(opTerms.begin(), opTerms.end());
@@ -68,8 +66,8 @@ namespace AMDiS {
// check if all terms are pw_const
bool pwConst = true;
- for (unsigned int i = 0; i < op->secondOrder[myRank].size(); i++) {
- if (!op->secondOrder[myRank][i]->isPWConst()) {
+ for (unsigned int i = 0; i < op->secondOrder.size(); i++) {
+ if (!op->secondOrder[i]->isPWConst()) {
pwConst = false;
break;
}
@@ -98,21 +96,9 @@ namespace AMDiS {
q11 = Q11PsiPhi::provideQ11PsiPhi(rowFeSpace->getBasisFcts(),
colFeSpace->getBasisFcts(),
- quadrature);
- tmpLALt.resize(omp_get_overall_max_threads());
- for (int i = 0; i < omp_get_overall_max_threads(); i++) {
- tmpLALt[i] = new DimMat<double>*;
- *(tmpLALt[i]) = new DimMat<double>(dim, NO_INIT);
- }
- }
-
-
- Pre2::~Pre2()
- {
- for (unsigned int i = 0; i < tmpLALt.size(); i++) {
- delete *(tmpLALt[i]);
- delete tmpLALt[i];
- }
+ quadrature);
+ LALt.resize(1);
+ LALt[0].change_dim(dim + 1, dim + 1);
}
@@ -122,15 +108,19 @@ namespace AMDiS {
const int *k, *l;
const double *values;
- int myRank = omp_get_thread_num();
- DimMat<double> **LALt = tmpLALt[0];
- DimMat<double> &tmpMat = *LALt[0];
- tmpMat.set(0.0);
+ mtl::dense2D<double> &tmpMat = LALt[0];
+ tmpMat = 0.0;
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, LALt);
+
+ // Compute: LALt[0] *= elInfo->getDet();
+ {
+ for (int i = 0; i <= dim; i++)
+ for (int j = 0; j <= dim; j++)
+ tmpMat[i][j] *= elInfo->getDet();
+ }
- tmpMat *= elInfo->getDet();
nEntries = q11->getNumberEntries();
if (symmetric) {
@@ -141,9 +131,8 @@ namespace AMDiS {
l = q11->getLVec(i, i);
values = q11->getValVec(i, i);
double val = 0.0;
- for (int m = 0; m < nEntries[i][i]; m++) {
- val += values[m] * tmpMat[k[m]][l[m]];
- }
+ for (int m = 0; m < nEntries[i][i]; m++)
+ val += values[m] * tmpMat[k[m]][l[m]];
mat[i][i] += val;
for (int j = i + 1; j < nCol; j++) {
@@ -151,9 +140,8 @@ namespace AMDiS {
l = q11->getLVec(i, j);
values = q11->getValVec(i, j);
val = 0.0;
- for (int m = 0; m < nEntries[i][j]; m++) {
- val += values[m] * tmpMat[k[m]][l[m]];
- }
+ for (int m = 0; m < nEntries[i][j]; m++)
+ val += values[m] * tmpMat[k[m]][l[m]];
mat[i][j] += val;
mat[j][i] += val;
}
@@ -165,9 +153,8 @@ namespace AMDiS {
l = q11->getLVec(i, j);
values = q11->getValVec(i, j);
double val = 0.0;
- for (int m = 0; m < nEntries[i][j]; m++) {
+ for (int m = 0; m < nEntries[i][j]; m++)
val += values[m] * tmpMat[k[m]][l[m]];
- }
mat[i][j] += val;
}
}
@@ -179,91 +166,102 @@ namespace AMDiS {
: SecondOrderAssembler(op, assembler, quad, true)
{
name = "fast quadrature second order assembler";
-
- tmpVec.resize(omp_get_overall_max_threads());
- tmpLALt.resize(omp_get_overall_max_threads());
- }
-
-
- Quad2::~Quad2()
- {
- if (!firstCall) {
- int nPoints = quadrature->getNumPoints();
- for (unsigned int i = 0; i < tmpLALt.size(); i++) {
- for (int j = 0; j < nPoints; j++)
- delete tmpLALt[i][j];
-
- delete [] tmpLALt[i];
- }
- }
}
void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
- int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
+ const int nPoints = quadrature->getNumPoints();
if (firstCall) {
- tmpVec[myRank] = new DimVec<double>(dim, NO_INIT);
- tmpLALt[myRank] = new DimMat<double>*[nPoints];
+ dimVec.change_dim(dim + 1);
+ LALt.resize(nPoints);
for (int j = 0; j < nPoints; j++)
- tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);
-
-#ifdef _OPENMP
-#pragma omp critical
-#endif
- {
- psiFast =
- updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
- phiFast =
- updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
- }
+ LALt[j].change_dim(dim + 1, dim + 1);
+
+ psiFast =
+ updateFastQuadrature(psiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
+ phiFast =
+ updateFastQuadrature(phiFast, rowFeSpace->getBasisFcts(), INIT_GRD_PHI);
+
firstCall = false;
}
- DimMat<double> **LALt = tmpLALt[myRank];
- DimVec<double> &dimVec = *(tmpVec[myRank]);
-
for (int i = 0; i < nPoints; i++)
- LALt[i]->set(0.0);
+ LALt[i] = 0.0;
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
-
- VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, LALt);
if (symmetric) {
// === Symmetric assembling. ===
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
for (int iq = 0; iq < nPoints; iq++) {
- (*LALt[iq]) *= elInfo->getDet();
- grdPsi = psiFast->getGradient(iq);
- grdPhi = phiFast->getGradient(iq);
-
- for (int i = 0; i < nCol; i++) {
- amv(quadrature->getWeight(iq), (*LALt[iq]), (*grdPhi)[i], dimVec);
-
- mat[i][i] += (*grdPsi)[i] * dimVec;
- for (int j = i + 1; j < nRow; j++) {
- double tmp = (*grdPhi)[j] * dimVec;
- mat[i][j] += tmp;
- mat[j][i] += tmp;
- }
+ // Compute: LALt[iq] *= elInfo->getDet();
+ for (int i = 0; i <= dim; i++)
+ for (int j = 0; j <= dim; j++)
+ LALt[iq][i][j] *= elInfo->getDet();
+
+ const vector<mtl::dense_vector<double> >& grdPsi = psiFast->getGradient(iq);
+ const vector<mtl::dense_vector<double> >& grdPhi = phiFast->getGradient(iq);
+ double weight = quadrature->getWeight(iq);
+
+ for (int i = 0; i < nCol; i++) {
+ // Compute: dimVec = quadrate->getWeight(iq) * (LALt[iq] * grdPhi[i])
+ // dimVec = LALt[iq] * grdPhi[i];
+ // dimVec *= quadrature->getWeight(iq);
+ double v = 0.0;
+ for (int j = 0; j <= dim; j++) {
+ v = 0.0;
+ for (int k = 0; k <= dim; k++)
+ v += LALt[iq][j][k] * grdPhi[i][k];
+ dimVec[j] = weight * v;
+ }
+
+ mat[i][i] += dot(dimVec, grdPsi[i]);
+ for (int j = i + 1; j < nRow; j++) {
+ double tmp = dot(dimVec, grdPhi[j]);
+ mat[i][j] += tmp;
+ mat[j][i] += tmp;
+ }
}
}
} else {
// === Non symmetric assembling. ===
+
for (int iq = 0; iq < nPoints; iq++) {
- (*LALt[iq]) *= elInfo->getDet();
- grdPsi = psiFast->getGradient(iq);
- grdPhi = phiFast->getGradient(iq);
-
- for (int i = 0; i < nRow; i++)
- for (int j = 0; j < nCol; j++)
- mat[i][j] += quadrature->getWeight(iq) *
- xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
+ // Compute: LALt[iq] *= elInfo->getDet();
+ {
+ for (int i = 0; i <= dim; i++)
+ for (int j = 0; j <= dim; j++)
+ LALt[iq][i][j] *= elInfo->getDet();
+ }
+
+ const vector<mtl::dense_vector<double> >& grdPsi = psiFast->getGradient(iq);
+ const vector<mtl::dense_vector<double> >& grdPhi = phiFast->getGradient(iq);
+
+ for (int i = 0; i < nRow; i++) {
+ const mtl::dense_vector<double>& grdPsi_i = grdPsi[i];
+
+ for (int j = 0; j < nCol; j++) {
+ const mtl::dense_vector<double>& grdPhi_j = grdPhi[j];
+
+ // Compute: mat[i][j] += quadrature->getWeight(iq) * (grdPsi[i] * (LALt[iq] * grdPhi[j]))
+ // dimVec = LALt[iq] * grdPhi[j];
+ // mat[i][j] += quadrature->getWeight(iq) * dot(grdPsi[i], dimVec);
+ {
+ double v = 0.0;
+ for (int k = 0; k <= dim; k++) {
+ double w = 0.0;
+ for (int l = 0; l <= dim; l++)
+ w += LALt[iq][k][l] * grdPhi_j[l];
+ v += grdPsi_i[k] * w;
+ }
+ mat[i][j] += quadrature->getWeight(iq) * v;
+ }
+ }
+ }
}
}
}
@@ -278,30 +276,32 @@ namespace AMDiS {
void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
- DimVec<double> grdPsi(dim, NO_INIT);
- VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
+ mtl::dense_vector<double> grdPsi(dim + 1);
+ vector<mtl::dense_vector<double> > grdPhi(nCol);
+ for (int i = 0; i < nCol; i++)
+ grdPhi[i].change_dim(dim + 1);
const BasisFunction *psi = rowFeSpace->getBasisFcts();
const BasisFunction *phi = colFeSpace->getBasisFcts();
int nPoints = quadrature->getNumPoints();
- DimMat<double> **LALt = new DimMat<double>*[nPoints];
+ std::vector<mtl::dense2D<double> > LALt(nPoints);
for (int iq = 0; iq < nPoints; iq++) {
- LALt[iq] = new DimMat<double>(dim, NO_INIT);
- LALt[iq]->set(0.0);
+ LALt[iq].change_dim(dim + 1, dim + 1);
+ LALt[iq] = 0.0;
}
- int myRank = omp_get_thread_num();
+ mtl::dense_vector<double> tmpVec;
- for (unsigned int i = 0; i < terms[myRank].size(); i++)
- (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
+ for (unsigned int i = 0; i < terms.size(); i++)
+ (static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, LALt);
if (symmetric) {
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
for (int iq = 0; iq < nPoints; iq++) {
- (*LALt[iq]) *= elInfo->getDet();
+ LALt[iq] *= elInfo->getDet();
for (int i = 0; i < nCol; i++)
(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
@@ -309,12 +309,12 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
- mat[i][i] += quadrature->getWeight(iq) *
- (grdPsi * ((*LALt[iq]) * grdPhi[i]));
+ tmpVec = LALt[iq] * grdPhi[i];
+ mat[i][i] += quadrature->getWeight(iq) * dot(grdPsi, tmpVec);
for (int j = i + 1; j < nCol; j++) {
- double val = quadrature->getWeight(iq) *
- (grdPsi * ((*LALt[iq]) * grdPhi[j]));
+ tmpVec = (LALt[iq] * grdPhi[j]);
+ double val = quadrature->getWeight(iq) * dot(grdPsi, tmpVec);
mat[i][j] += val;
mat[j][i] += val;
}
@@ -322,7 +322,7 @@ namespace AMDiS {
}
} else { /* non symmetric assembling */
for (int iq = 0; iq < nPoints; iq++) {
- (*LALt[iq]) *= elInfo->getDet();
+ LALt[iq] *= elInfo->getDet();
for (int i = 0; i < nCol; i++)
(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
@@ -330,17 +330,12 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
for (int j = 0; j < nCol; j++) {
- mat[i][j] += quadrature->getWeight(iq) *
- (grdPsi * ((*LALt[iq]) * grdPhi[j]));
+ tmpVec = (LALt[iq] * grdPhi[j]);
+ mat[i][j] += quadrature->getWeight(iq) * dot(grdPsi, tmpVec);
}
}
}
}
-
- for (int iq = 0; iq < nPoints; iq++)
- delete LALt[iq];
-
- delete [] LALt;
}
}
diff --git a/AMDiS/src/SecondOrderAssembler.h b/AMDiS/src/SecondOrderAssembler.h
index e3da7f5618ddf9e0366cf0de5db0c0a0f6561c81..8fc6bc6ed26c9eb9b43284b914ace191d74eb61e 100644
--- a/AMDiS/src/SecondOrderAssembler.h
+++ b/AMDiS/src/SecondOrderAssembler.h
@@ -23,11 +23,14 @@
#ifndef AMDIS_SECOND_ORDER_ASSEMBLER_H
#define AMDIS_SECOND_ORDER_ASSEMBLER_H
+#include <vector>
+#include "AMDiS_fwd.h"
+#include "QPsiPhi.h"
#include "SubAssembler.h"
namespace AMDiS {
- class Q11PsiPhi;
+ using namespace std;
/**
* \ingroup Assembler
@@ -62,10 +65,10 @@ namespace AMDiS {
protected:
/// List of all yet created optimized second order assemblers.
- static std::vector<SubAssembler*> optimizedSubAssemblers;
+ static vector<SubAssembler*> optimizedSubAssemblers;
/// List of all yet created standard second order assemblers.
- static std::vector<SubAssembler*> standardSubAssemblers;
+ static vector<SubAssembler*> standardSubAssemblers;
};
@@ -103,8 +106,6 @@ namespace AMDiS {
/// Constructor.
Quad2(Operator *op, Assembler *assembler, Quadrature *quad);
- ~Quad2();
-
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
@@ -115,13 +116,9 @@ namespace AMDiS {
}
protected:
- std::vector<DimVec<double>*> tmpVec;
+ mtl::dense_vector<double> dimVec;
- /** \brief
- * Thread safe temporary vector of DimMats for calculation in
- * function calculateElementMatrix().
- */
- std::vector<DimMat<double>**> tmpLALt;
+ vector<mtl::dense2D<double> > LALt;
};
/**
@@ -136,9 +133,6 @@ namespace AMDiS {
/// Constructor.
Pre2(Operator *op, Assembler *assembler, Quadrature *quad);
- /// Destructor.
- ~Pre2();
-
/// Implements SubAssembler::calculateElementMatrix().
void calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat);
@@ -152,11 +146,7 @@ namespace AMDiS {
/// Integral of the product of the derivative of psi and the derivative of phi.
const Q11PsiPhi *q11;
- /** \brief
- * Thread safe temporary vector of DimMats for calculation in
- * function calculateElementMatrix().
- */
- std::vector< DimMat<double>** > tmpLALt;
+ vector<mtl::dense2D<double> > LALt;
friend class SecondOrderAssembler;
};
diff --git a/AMDiS/src/SecondOrderTerm.cc b/AMDiS/src/SecondOrderTerm.cc
index 57af28831f7a497b2716139330a34227578981da..e681db44c4f64e17b96475f41ead151c1b0bd0a8 100644
--- a/AMDiS/src/SecondOrderTerm.cc
+++ b/AMDiS/src/SecondOrderTerm.cc
@@ -17,11 +17,11 @@ namespace AMDiS {
void SecondOrderTerm::lalt(const DimVec<WorldVector<double> >& Lambda,
const WorldMatrix<double>& matrix,
- DimMat<double>& LALt,
+ mtl::dense2D<double>& LALt,
bool symm,
double factor) const
{
- int dim = LALt.getNumRows();
+ int dim = num_rows(LALt);
if (symm) {
for (int i = 0; i < dim; i++) {
@@ -58,10 +58,10 @@ namespace AMDiS {
void SecondOrderTerm::lalt_kl(const DimVec<WorldVector<double> >& Lambda,
int k, int l,
- DimMat<double>& LALt,
+ mtl::dense2D<double>& LALt,
double factor)
{
- int dim = LALt.getNumRows();
+ int dim = num_rows(LALt);
for (int i = 0; i < dim; i++)
for (int j = 0; j < dim; j++)
@@ -95,12 +95,14 @@ namespace AMDiS {
getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
}
- void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void MatrixFct_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt(grdLambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*matrixFct)(vecAtQPs[iq]), LALt[iq], symmetric, 1.0);
}
void MatrixFct_SOT::eval(int nPoints,
@@ -201,12 +203,14 @@ namespace AMDiS {
getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs);
}
- void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void VecAtQP_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
+ l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq]));
}
void VecAtQP_SOT::eval(int nPoints,
@@ -220,11 +224,10 @@ namespace AMDiS {
if (D2UhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
- double factor = (*f)(vecAtQPs[iq]);
double resultQP = 0.0;
for (int i = 0; i < dow; i++)
resultQP += D2UhAtQP[iq][i][i];
- result[iq] += fac * factor * resultQP;
+ result[iq] += fac * (*f)(vecAtQPs[iq]) * resultQP;
}
}
}
@@ -263,12 +266,14 @@ namespace AMDiS {
getVectorAtQPs(vec2, elInfo, subAssembler, quad, vecAtQPs2);
}
- void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
+ l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
}
void Vec2AtQP_SOT::eval(int nPoints,
@@ -311,12 +316,14 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
+ l1lt(grdLambda, LALt[iq], (*g)(coordsAtQPs[iq]));
}
void CoordsAtQP_SOT::eval(int nPoints,
@@ -385,12 +392,14 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
}
- void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void MatrixGradient_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt(grdLambda, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*f)(gradAtQPs[iq]), LALt[iq], symmetric, 1.0);
}
void MatrixGradient_SOT::eval(int nPoints,
@@ -454,12 +463,14 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
- void FctGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void FctGradient_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
+ l1lt(grdLambda, LALt[iq], (*f)(gradAtQPs[iq]));
}
void FctGradient_SOT::eval(int nPoints,
@@ -517,12 +528,14 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
- void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt(grdLambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), LALt[iq], symmetric, 1.0);
}
void VecMatrixGradientAtQP_SOT::eval(int nPoints,
@@ -584,12 +597,14 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
+ l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
}
void VecGradCoordsAtQP_SOT::eval(int nPoints,
@@ -644,12 +659,14 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
+ l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
}
void VecAndCoordsAtQP_SOT::eval(int nPoints,
@@ -708,12 +725,14 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt(grdLambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), LALt[iq], symmetric, 1.0);
}
void MatrixGradientAndCoords_SOT::eval(int nPoints,
@@ -786,12 +805,14 @@ namespace AMDiS {
getVectorAtQPs(vec2, elInfo, subAssembler, quad, vec2AtQPs);
}
- void MatrixVec2_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void MatrixVec2_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt(grdLambda, A, *(LALt[iq]), symmetric, (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]));
+ lalt(grdLambda, A, LALt[iq], symmetric, (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]));
}
void MatrixVec2_SOT::eval(int nPoints,
@@ -874,14 +895,14 @@ namespace AMDiS {
}
void General_SOT::getLALt(const ElInfo *elInfo,
- int nPoints,
- DimMat<double> **LALt) const
+ vector<mtl::dense2D<double> > &LALt) const
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
+ const int nPoints = static_cast<int>(LALt.size());
- std::vector<double> vecsArg(nVecs);
- std::vector<WorldVector<double> > gradsArg(nGrads);
+ vector<double> vecsArg(nVecs);
+ vector<WorldVector<double> > gradsArg(nGrads);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
@@ -892,7 +913,7 @@ namespace AMDiS {
gradsArg[i] = gradsAtQPs_[i][iq];
lalt(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
- *(LALt[iq]), symmetric_, 1.0);
+ LALt[iq], symmetric_, 1.0);
}
}
@@ -1007,14 +1028,15 @@ namespace AMDiS {
gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
- void GeneralParametric_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void GeneralParametric_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
+ const int nPoints = static_cast<int>(LALt.size());
- std::vector<double> vecsArg(nVecs);
- std::vector<WorldVector<double> > gradsArg(nGrads);
+ vector<double> vecsArg(nVecs);
+ vector<WorldVector<double> > gradsArg(nGrads);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
@@ -1024,7 +1046,7 @@ namespace AMDiS {
for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
lalt(grdLambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg),
- *(LALt[iq]), symmetric_, 1.0);
+ LALt[iq], symmetric_, 1.0);
}
}
@@ -1108,12 +1130,14 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
- void VecAndGradAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void VecAndGradAtQP_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
+ l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq], gradAtQPs[iq]));
}
void VecAndGradAtQP_SOT::eval(int nPoints,
@@ -1158,12 +1182,14 @@ namespace AMDiS {
gradAtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
}
- void VecGrad_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const
+ void VecGrad_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
+ l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq], gradAtQPs[iq]));
}
void VecGrad_SOT::eval(int nPoints,
@@ -1207,12 +1233,13 @@ namespace AMDiS {
}
void CoordsAtQP_IJ_SOT::getLALt(const ElInfo *elInfo,
- int nPoints,
- DimMat<double> **LALt) const
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > Lambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*g)(coordsAtQPs[iq]));
+ lalt_kl(Lambda, xi, xj, LALt[iq], (*g)(coordsAtQPs[iq]));
}
void CoordsAtQP_IJ_SOT::eval(int nPoints,
@@ -1263,13 +1290,13 @@ namespace AMDiS {
}
void VecAtQP_IJ_SOT::getLALt(const ElInfo *elInfo,
- int nPoints,
- DimMat<double> **LALt) const
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
- lalt_kl(grdLambda, xi, xj, *(LALt[iq]), (*f)(vecAtQPs[iq]));
- }
+ const int nPoints = static_cast<int>(LALt.size());
+
+ for (int iq = 0; iq < nPoints; iq++)
+ lalt_kl(grdLambda, xi, xj, LALt[iq], (*f)(vecAtQPs[iq]));
}
void VecAtQP_IJ_SOT::eval(int nPoints,
diff --git a/AMDiS/src/SecondOrderTerm.h b/AMDiS/src/SecondOrderTerm.h
index d2ad39918d259d9ca41767dddf4c8454a7e5841b..bf70218b0f48473215045912a2e2673740e869c4 100644
--- a/AMDiS/src/SecondOrderTerm.h
+++ b/AMDiS/src/SecondOrderTerm.h
@@ -30,6 +30,8 @@
namespace AMDiS {
+ using namespace std;
+
/**
* \ingroup Assembler
*
@@ -50,8 +52,7 @@ namespace AMDiS {
/// Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
virtual void getLALt(const ElInfo *elInfo,
- int nPoints,
- DimMat<double> **result) const = 0;
+ vector<mtl::dense2D<double> > &result) const = 0;
/// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
virtual void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -61,7 +62,7 @@ namespace AMDiS {
/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$.
void lalt(const DimVec<WorldVector<double> >& Lambda,
const WorldMatrix<double>& matrix,
- DimMat<double>& LALt,
+ mtl::dense2D<double>& LALt,
bool symm,
double factor) const;
@@ -72,15 +73,15 @@ namespace AMDiS {
*/
static void lalt_kl(const DimVec<WorldVector<double> >& Lambda,
int k, int l,
- DimMat<double>& LALt,
+ mtl::dense2D<double>& LALt,
double factor);
/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
inline void l1lt(const DimVec<WorldVector<double> >& Lambda,
- DimMat<double>& LALt,
+ mtl::dense2D<double>& LALt,
double factor) const
{
- const int dim = LALt.getNumRows();
+ const int dim = num_rows(LALt);
for (int i = 0; i < dim; i++) {
double val = 0.0;
@@ -119,12 +120,13 @@ namespace AMDiS {
}
/// Implenetation of SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
+ inline void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), 1.0);
+ l1lt(grdLambda, LALt[iq], 1.0);
}
/// Implementation of SecondOrderTerm::eval().
@@ -187,11 +189,13 @@ namespace AMDiS {
}
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
+ inline void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- l1lt(grdLambda, *(LALt[iq]), (*factor));
+ l1lt(grdLambda, LALt[iq], (*factor));
}
/// Implenetation of SecondOrderTerm::eval().
@@ -252,7 +256,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -302,11 +306,13 @@ namespace AMDiS {
}
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
+ inline void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> >& Lambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0);
+ lalt(Lambda, matrix, LALt[iq], symmetric, 1.0);
}
/// Implenetation of SecondOrderTerm::eval().
@@ -359,11 +365,14 @@ namespace AMDiS {
}
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
+ inline void getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
{
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
for (int iq = 0; iq < nPoints; iq++)
- lalt_kl(grdLambda, xi, xj, *(LALt[iq]), (*factor));
+ lalt_kl(grdLambda, xi, xj, LALt[iq], (*factor));
}
/// Implenetation of SecondOrderTerm::eval().
@@ -405,7 +414,8 @@ namespace AMDiS {
* points of a given DOFVector:
* \f$ f(v(\vec{x})) \Delta u(\vec{x}) \f$
*/
- class VecAtQP_SOT : public SecondOrderTerm {
+ class VecAtQP_SOT : public SecondOrderTerm
+ {
public:
/// Constructor.
VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af);
@@ -421,7 +431,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -454,7 +464,8 @@ namespace AMDiS {
* points of a given DOFVector:
* \f$ f(v(\vec{x}), w(\vec{x})) \Delta u(\vec{x}) \f$
*/
- class Vec2AtQP_SOT : public SecondOrderTerm {
+ class Vec2AtQP_SOT : public SecondOrderTerm
+ {
public:
/// Constructor.
Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
@@ -465,7 +476,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -498,7 +509,8 @@ namespace AMDiS {
* Laplace multiplied with a function evaluated at each quadrature point:
* \f$ f(\vec{x}) \Delta u(\vec{x}) \f$
*/
- class CoordsAtQP_SOT : public SecondOrderTerm {
+ class CoordsAtQP_SOT : public SecondOrderTerm
+ {
public:
/// Constructor.
CoordsAtQP_SOT(AbstractFunction<double, WorldVector<double> > *af)
@@ -512,7 +524,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -564,7 +576,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -621,7 +633,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -671,7 +683,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -727,7 +739,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -779,7 +791,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::eval().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
void eval(int nPoints,
const ElementVector& uhAtQP,
@@ -830,7 +842,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -880,7 +892,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -930,7 +942,7 @@ namespace AMDiS {
void initElement(const ElInfo*, SubAssembler*, Quadrature *quad= NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -988,7 +1000,7 @@ namespace AMDiS {
Quadrature *quad= NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
@@ -1048,7 +1060,7 @@ namespace AMDiS {
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implements SecondOrderTerm::eval().
void eval(int nPoints,
@@ -1098,7 +1110,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implements SecondOrderTerm::eval().
void eval(int nPoints,
@@ -1152,7 +1164,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implements SecondOrderTerm::eval().
void eval(int nPoints,
@@ -1198,7 +1210,7 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getLALt().
- void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
/// Implements SecondOrderTerm::eval().
void eval(int nPoints,
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index ed6725a7fb0c4bdb7ea5836f5381682e6c4eca4c..b6c06b3e297c95ee740e393b808d7889542b1a95 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -53,9 +53,6 @@ namespace AMDiS {
nRow = psi->getNumber();
nCol = phi->getNumber();
- int maxThreads = omp_get_overall_max_threads();
- terms.resize(maxThreads);
-
switch (order) {
case 0:
terms = op->zeroOrder;
@@ -77,8 +74,8 @@ namespace AMDiS {
// quadratic, and therefore cannot be symmetric.
if (nRow == nCol) {
symmetric = true;
- for (int i = 0; i < static_cast<int>(terms[0].size()); i++) {
- if (!(terms[0][i])->isSymmetric()) {
+ for (unsigned int i = 0; i < terms.size(); i++) {
+ if (!(terms[i])->isSymmetric()) {
symmetric = false;
break;
}
@@ -125,10 +122,9 @@ namespace AMDiS {
for (it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2)
((*it2).second)->valid = false;
- int myRank = omp_get_thread_num();
// calls initElement of each term
- std::vector<OperatorTerm*>::iterator it;
- for (it = terms[myRank].begin(); it != terms[myRank].end(); ++it) {
+ for (std::vector<OperatorTerm*>::iterator it = terms.begin();
+ it != terms.end(); ++it) {
if (largeElInfo == NULL)
(*it)->initElement(smallElInfo, this, quad);
else
diff --git a/AMDiS/src/SubAssembler.h b/AMDiS/src/SubAssembler.h
index b69a4bf9be2a984d4c3e786956b13b23d63ae499..94da76c82b705b57ca483e46e4c96891fd635ca0 100644
--- a/AMDiS/src/SubAssembler.h
+++ b/AMDiS/src/SubAssembler.h
@@ -28,7 +28,6 @@
#include "AMDiS_fwd.h"
#include "FixVec.h"
#include "Flag.h"
-#include "OpenMP.h"
namespace AMDiS {
@@ -81,7 +80,7 @@ namespace AMDiS {
/// Returns \ref terms
inline std::vector<OperatorTerm*> *getTerms()
{
- return &terms[omp_get_thread_num()];
+ return &terms;
}
/// Returns \ref quadrature.
@@ -246,7 +245,7 @@ namespace AMDiS {
bool symmetric;
/// List of all terms with a contribution to this SubAssembler
- std::vector< std::vector<OperatorTerm*> > terms;
+ std::vector<OperatorTerm*> terms;
///
bool opt;
diff --git a/AMDiS/src/SurfaceOperator.h b/AMDiS/src/SurfaceOperator.h
index c907eaedda9085cea803a631b91e9cfb37d97054..469c8bc6d4592aabee5ec0aba2be2962bdf3e158 100644
--- a/AMDiS/src/SurfaceOperator.h
+++ b/AMDiS/src/SurfaceOperator.h
@@ -27,7 +27,6 @@
#include "Mesh.h"
#include "Operator.h"
#include "SurfaceQuadrature.h"
-#include "OpenMP.h"
namespace AMDiS {
@@ -56,32 +55,31 @@ namespace AMDiS {
quad1GrdPhi(NULL),
quad0(NULL)
{
- assembler[omp_get_thread_num()] = NULL;
+ assembler = NULL;
int dim = rowFeSpace->getMesh()->getDim();
int degree;
- int myRank = omp_get_thread_num();
// create surface quadratures
- if (secondOrder[myRank].size() > 0) {
+ if (secondOrder.size() > 0) {
degree = getQuadratureDegree(2);
quad2 =
new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree), coords);
}
- if (firstOrderGrdPsi[myRank].size() > 0) {
+ if (firstOrderGrdPsi.size() > 0) {
degree = getQuadratureDegree(1, GRD_PSI);
quad1GrdPsi =
new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree), coords);
}
- if (firstOrderGrdPhi[myRank].size() > 0) {
+ if (firstOrderGrdPhi.size() > 0) {
degree = getQuadratureDegree(1, GRD_PHI);
quad1GrdPhi =
new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree), coords);
}
- if (zeroOrder[myRank].size() > 0) {
+ if (zeroOrder.size() > 0) {
degree = getQuadratureDegree(0);
quad0 =
new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree), coords);
@@ -89,7 +87,7 @@ namespace AMDiS {
// initialize assembler with surface quadratures
optimized = false;
- initAssembler(myRank, quad2, quad1GrdPsi, quad1GrdPhi, quad0);
+ initAssembler(quad2, quad1GrdPsi, quad1GrdPhi, quad0);
}
/// Adapt surface quadratures to \ref coords.
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index a0cba42562425c91625c7a5f3f433522baeb1780..51e68be0968ad2eb5f6e7feba439cecfb63ae913 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -39,10 +39,8 @@ namespace AMDiS {
Quadrature *quad,
bool optimized)
{
- int myRank = omp_get_thread_num();
-
// check if an assembler is needed at all
- if (op->zeroOrder[myRank].size() == 0)
+ if (op->zeroOrder.size() == 0)
return NULL;
ZeroOrderAssembler *newAssembler;
@@ -50,7 +48,7 @@ namespace AMDiS {
std::vector<SubAssembler*> *subAssemblers =
optimized ? &optimizedSubAssemblers : &standardSubAssemblers;
- std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
+ std::vector<OperatorTerm*> opTerms = op->zeroOrder;
sort(opTerms.begin(), opTerms.end());
@@ -68,8 +66,8 @@ namespace AMDiS {
// check if all terms are pw_const
bool pwConst = true;
- for (int i = 0; i < static_cast<int>(op->zeroOrder[myRank].size()); i++) {
- if (!op->zeroOrder[myRank][i]->isPWConst()) {
+ for (int i = 0; i < static_cast<int>(op->zeroOrder.size()); i++) {
+ if (!op->zeroOrder[i]->isPWConst()) {
pwConst = false;
break;
}
@@ -106,9 +104,8 @@ namespace AMDiS {
ElementVector c(nPoints, 0.0);
std::vector<double> phival(nCol);
- int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
if (symmetric) {
@@ -153,10 +150,9 @@ namespace AMDiS {
{
int nPoints = quadrature->getNumPoints();
ElementVector c(nPoints, 0.0);
- int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
for (int iq = 0; iq < nPoints; iq++) {
@@ -175,14 +171,12 @@ namespace AMDiS {
: ZeroOrderAssembler(op, assembler, quad, true)
{
name = "fast quadrature zero order assembler";
- tmpC.resize(omp_get_overall_max_threads());
}
void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
int nPoints = quadrature->getNumPoints();
- int myRank = omp_get_thread_num();
if (firstCall) {
#ifdef _OPENMP
@@ -197,49 +191,49 @@ namespace AMDiS {
}
}
- ElementVector &c = tmpC[myRank];
if (num_rows(c) < static_cast<unsigned int>(nPoints))
c.change_dim(nPoints);
c = 0.0;
- for (std::vector<OperatorTerm*>::iterator termIt = terms[myRank].begin();
- termIt != terms[myRank].end(); ++termIt)
+ for (std::vector<OperatorTerm*>::iterator termIt = terms.begin();
+ termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
+ const mtl::dense2D<double>& psi = psiFast->getPhi();
+ const mtl::dense2D<double>& phi = phiFast->getPhi();
+
if (symmetric) {
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);
- const double *psi = psiFast->getPhi(iq);
- const double *phi = phiFast->getPhi(iq);
for (int i = 0; i < nRow; i++) {
- mat[i][i] += c[iq] * psi[i] * phi[i];
+ mat[i][i] += c[iq] * psi[iq][i] * phi[iq][i];
for (int j = i + 1; j < nCol; j++) {
- double val = c[iq] * psi[i] * phi[j];
+ double val = c[iq] * psi[iq][i] * phi[iq][j];
mat[i][j] += val;
mat[j][i] += val;
}
}
}
} else { /* non symmetric assembling */
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);
- const double *psi = psiFast->getPhi(iq);
- const double *phi = phiFast->getPhi(iq);
- for (int i = 0; i < nRow; i++)
- for (int j = 0; j < nCol; j++)
- mat[i][j] += c[iq] * psi[i] * phi[j];
- }
+ for (int i = 0; i < nRow; i++)
+ for (int j = 0; j < nCol; j++) {
+ double v = 0.0;
+ for (int iq = 0; iq < nPoints; iq++)
+ v += c[iq] * psi[iq][i] * phi[iq][j];
+ mat[i][j] += v;
+ }
}
}
void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
- int myRank = omp_get_thread_num();
int nPoints = quadrature->getNumPoints();
if (firstCall) {
@@ -257,15 +251,16 @@ namespace AMDiS {
ElementVector c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
+ const mtl::dense2D<double>& psi = psiFast->getPhi();
+
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
- const double *psi = psiFast->getPhi(iq);
for (int i = 0; i < nRow; i++)
- vec[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
+ vec[i] += quadrature->getWeight(iq) * c[iq] * psi[iq][i];
}
}
@@ -293,11 +288,10 @@ namespace AMDiS {
}
ElementVector c(1, 0.0);
- int myRank = omp_get_thread_num();
- int size = static_cast<int>(terms[myRank].size());
+ int size = static_cast<int>(terms.size());
for (int i = 0; i < size; i++)
- (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, c);
+ (static_cast<ZeroOrderTerm*>((terms[i])))->getC(elInfo, 1, c);
c[0] *= elInfo->getDet();
@@ -337,9 +331,8 @@ namespace AMDiS {
std::vector<OperatorTerm*>::iterator termIt;
- int myRank = omp_get_thread_num();
ElementVector c(1, 0.0);
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>(*termIt))->getC(elInfo, 1, c);
c[0] *= elInfo->getDet();
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index c14e129d407a1a8959295beedad75ecd832d0c0c..703479fadaa88d2e7392615e3de2fb5629d5de74 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -109,7 +109,7 @@ namespace AMDiS {
void calculateElementVector(const ElInfo *elInfo, ElementVector& vec);
protected:
- std::vector<ElementVector> tmpC;
+ ElementVector c;
};
diff --git a/AMDiS/src/io/FileWriter.cc b/AMDiS/src/io/FileWriter.cc
index b7f1cdf7f2d9fd8c0f5df70a8664a0bd13a7b61f..44a43d083fbf81d8402b9dd573b8544b18efc78c 100644
--- a/AMDiS/src/io/FileWriter.cc
+++ b/AMDiS/src/io/FileWriter.cc
@@ -120,6 +120,7 @@ namespace AMDiS {
writePeriodicFormat = 0;
writePngFormat = 0;
writePovrayFormat = 0;
+ writeDofFormat = 0;
pngType = 0;
appendIndex = 0;
indexLength = 5;