diff --git a/AMDiS/libtool b/AMDiS/libtool
index e96f67065e881c5fbc2155783211a810082e9175..36f4947b7437916e018387a8df916c68568f5eaf 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -30,21 +30,21 @@
# the same distribution terms that you use for the rest of that program.
# A sed program that does not truncate output.
-SED="/bin/sed"
+SED="/usr/bin/sed"
# Sed that helps us avoid accidentally triggering echo(1) options like -n.
-Xsed="/bin/sed -e 1s/^X//"
+Xsed="/usr/bin/sed -e 1s/^X//"
# The HP-UX ksh and POSIX shell print the target directory to stdout
# if CDPATH is set.
(unset CDPATH) >/dev/null 2>&1 && unset CDPATH
# The names of the tagged configurations supported by this script.
-available_tags=" CXX F77"
+available_tags=" CXX"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host linux-5k82:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -66,12 +66,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -93,14 +93,14 @@ CC="gcc"
# Is the compiler the GNU C compiler?
with_gcc=yes
-gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_dir=`gcc -print-file-name=. | /usr/bin/sed 's,/\.$,,'`
gcc_ver=`gcc -dumpversion`
# An ERE matcher.
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -328,10 +328,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=`echo " /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+sys_lib_search_path_spec=`echo "/lib64 /usr/lib64 /usr/local/lib64" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib64 /usr/lib64 /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/lib64/Xaw3d /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/lib/Xaw3d /usr/x86_64-suse-linux/lib /usr/local/lib /opt/kde3/lib /lib64 /lib /usr/lib64 /usr/lib /usr/local/lib64 /opt/kde3/lib64 /usr/lib64/graphviz /usr/lib64/graphviz/sharp /usr/lib64/graphviz/java /usr/lib64/graphviz/perl /usr/lib64/graphviz/php /usr/lib64/graphviz/ocaml /usr/lib64/graphviz/python /usr/lib64/graphviz/lua /usr/lib64/graphviz/tcl /usr/lib64/graphviz/guile /usr/lib64/graphviz/ruby "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -6763,7 +6763,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host NWRW15:
+# Libtool was configured on host linux-5k82:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -6785,12 +6785,12 @@ fast_install=yes
# The host system.
host_alias=
-host=i686-pc-linux-gnu
+host=x86_64-unknown-linux-gnu
host_os=linux-gnu
# The build system.
build_alias=
-build=i686-pc-linux-gnu
+build=x86_64-unknown-linux-gnu
build_os=linux-gnu
# An echo program that does not interpret backslashes.
@@ -6812,14 +6812,14 @@ CC="g++"
# Is the compiler the GNU C compiler?
with_gcc=yes
-gcc_dir=`gcc -print-file-name=. | /bin/sed 's,/\.$,,'`
+gcc_dir=`gcc -print-file-name=. | /usr/bin/sed 's,/\.$,,'`
gcc_ver=`gcc -dumpversion`
# An ERE matcher.
EGREP="grep -E"
# The linker used to build libraries.
-LD="/usr/bin/ld"
+LD="/usr/x86_64-suse-linux/bin/ld -m elf_x86_64"
# Whether we need hard or soft links.
LN_S="ln -s"
@@ -6948,11 +6948,11 @@ striplib="strip --strip-unneeded"
# Dependencies to place before the objects being linked to create a
# shared library.
-predep_objects=`echo "/usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crti.o /usr/lib/gcc/i386-redhat-linux/4.1.2/crtbeginS.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+predep_objects=`echo "/usr/lib64/gcc/x86_64-suse-linux/4.4/../../../../lib64/crti.o /usr/lib64/gcc/x86_64-suse-linux/4.4/crtbeginS.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place after the objects being linked to create a
# shared library.
-postdep_objects=`echo "/usr/lib/gcc/i386-redhat-linux/4.1.2/crtendS.o /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../crtn.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+postdep_objects=`echo "/usr/lib64/gcc/x86_64-suse-linux/4.4/crtendS.o /usr/lib64/gcc/x86_64-suse-linux/4.4/../../../../lib64/crtn.o" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Dependencies to place before the objects being linked to create a
# shared library.
@@ -6964,7 +6964,7 @@ 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=`echo "-L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+compiler_lib_search_path=`echo "-L/usr/lib64/gcc/x86_64-suse-linux/4.4 -L/usr/lib64/gcc/x86_64-suse-linux/4.4/../../../../lib64 -L/lib/../lib64 -L/usr/lib/../lib64 -L/usr/lib64/gcc/x86_64-suse-linux/4.4/../../../../x86_64-suse-linux/lib -L/usr/lib64/gcc/x86_64-suse-linux/4.4/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Method to check whether dependent libraries are shared objects.
deplibs_check_method="pass_all"
@@ -7044,10 +7044,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=`echo " /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+sys_lib_search_path_spec=`echo "/lib64 /usr/lib64 /usr/local/lib64" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
+sys_lib_dlsearch_path_spec="/lib64 /usr/lib64 /usr/X11R6/lib64/Xaw3d /usr/X11R6/lib64 /usr/lib64/Xaw3d /usr/X11R6/lib/Xaw3d /usr/X11R6/lib /usr/lib/Xaw3d /usr/x86_64-suse-linux/lib /usr/local/lib /opt/kde3/lib /lib64 /lib /usr/lib64 /usr/lib /usr/local/lib64 /opt/kde3/lib64 /usr/lib64/graphviz /usr/lib64/graphviz/sharp /usr/lib64/graphviz/java /usr/lib64/graphviz/perl /usr/lib64/graphviz/php /usr/lib64/graphviz/ocaml /usr/lib64/graphviz/python /usr/lib64/graphviz/lua /usr/lib64/graphviz/tcl /usr/lib64/graphviz/guile /usr/lib64/graphviz/ruby "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -7069,314 +7069,3 @@ include_expsyms=""
# ### END LIBTOOL TAG CONFIG: CXX
-# ### BEGIN LIBTOOL TAG CONFIG: F77
-
-# Libtool was configured on host NWRW15:
-
-# Shell to use when invoking shell scripts.
-SHELL="/bin/sh"
-
-# Whether or not to build shared libraries.
-build_libtool_libs=yes
-
-# Whether or not to build static libraries.
-build_old_libs=yes
-
-# Whether or not to add -lc for building shared libraries.
-build_libtool_need_lc=no
-
-# Whether or not to disallow shared libs when runtime libs are static
-allow_libtool_libs_with_static_runtimes=no
-
-# Whether or not to optimize for fast installation.
-fast_install=yes
-
-# The host system.
-host_alias=
-host=i686-pc-linux-gnu
-host_os=linux-gnu
-
-# The build system.
-build_alias=
-build=i686-pc-linux-gnu
-build_os=linux-gnu
-
-# An echo program that does not interpret backslashes.
-echo="echo"
-
-# The archiver.
-AR="ar"
-AR_FLAGS="cru"
-
-# A C compiler.
-LTCC="gcc"
-
-# LTCC compiler flags.
-LTCFLAGS="-g -O2"
-
-# A language-specific compiler.
-CC="g77"
-
-# 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/bin/ld"
-
-# Whether we need hard or soft links.
-LN_S="ln -s"
-
-# A BSD-compatible nm program.
-NM="/usr/bin/nm -B"
-
-# A symbol stripping program
-STRIP="strip"
-
-# Used to examine libraries when file_magic_cmd begins "file"
-MAGIC_CMD=file
-
-# Used on cygwin: DLL creation program.
-DLLTOOL="dlltool"
-
-# Used on cygwin: object dumper.
-OBJDUMP="objdump"
-
-# Used on cygwin: assembler.
-AS="as"
-
-# The name of the directory that contains temporary libtool files.
-objdir=.libs
-
-# How to create reloadable object files.
-reload_flag=" -r"
-reload_cmds="\$LD\$reload_flag -o \$output\$reload_objs"
-
-# How to pass a linker flag through the compiler.
-wl="-Wl,"
-
-# Object file suffix (normally "o").
-objext="o"
-
-# Old archive suffix (normally "a").
-libext="a"
-
-# Shared library suffix (normally ".so").
-shrext_cmds='.so'
-
-# Executable file suffix (normally "").
-exeext=""
-
-# Additional compiler flags for building library objects.
-pic_flag=" -fPIC"
-pic_mode=default
-
-# What is the maximum length of a command?
-max_cmd_len=32768
-
-# Does compiler simultaneously support -c and -o options?
-compiler_c_o="yes"
-
-# Must we lock files when doing compilation?
-need_locks="no"
-
-# Do we need the lib prefix for modules?
-need_lib_prefix=no
-
-# Do we need a version for libraries?
-need_version=no
-
-# Whether dlopen is supported.
-dlopen_support=unknown
-
-# Whether dlopen of programs is supported.
-dlopen_self=unknown
-
-# Whether dlopen of statically linked programs is supported.
-dlopen_self_static=unknown
-
-# Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
-
-# Compiler flag to turn off builtin functions.
-no_builtin_flag=""
-
-# Compiler flag to allow reflexive dlopens.
-export_dynamic_flag_spec="\${wl}--export-dynamic"
-
-# Compiler flag to generate shared objects directly from archives.
-whole_archive_flag_spec="\${wl}--whole-archive\$convenience \${wl}--no-whole-archive"
-
-# Compiler flag to generate thread-safe objects.
-thread_safe_flag_spec=""
-
-# Library versioning type.
-version_type=linux
-
-# Format of library name prefix.
-libname_spec="lib\$name"
-
-# List of archive names. First name is the real one, the rest are links.
-# The last name is the one that the linker finds with -lNAME.
-library_names_spec="\${libname}\${release}\${shared_ext}\$versuffix \${libname}\${release}\${shared_ext}\$major \$libname\${shared_ext}"
-
-# The coded name of the library, if different from the real name.
-soname_spec="\${libname}\${release}\${shared_ext}\$major"
-
-# Commands used to build and install an old-style archive.
-RANLIB="ranlib"
-old_archive_cmds="\$AR \$AR_FLAGS \$oldlib\$oldobjs\$old_deplibs~\$RANLIB \$oldlib"
-old_postinstall_cmds="chmod 644 \$oldlib~\$RANLIB \$oldlib"
-old_postuninstall_cmds=""
-
-# Create an old-style archive from a shared archive.
-old_archive_from_new_cmds=""
-
-# Create a temporary old-style archive to link instead of a shared archive.
-old_archive_from_expsyms_cmds=""
-
-# Commands used to build and install a shared archive.
-archive_cmds="\$CC -shared \$libobjs \$deplibs \$compiler_flags \${wl}-soname \$wl\$soname -o \$lib"
-archive_expsym_cmds="\$echo \\\"{ global:\\\" > \$output_objdir/\$libname.ver~
- cat \$export_symbols | sed -e \\\"s/\\\\(.*\\\\)/\\\\1;/\\\" >> \$output_objdir/\$libname.ver~
- \$echo \\\"local: *; };\\\" >> \$output_objdir/\$libname.ver~
- \$CC -shared \$libobjs \$deplibs \$compiler_flags \${wl}-soname \$wl\$soname \${wl}-version-script \${wl}\$output_objdir/\$libname.ver -o \$lib"
-postinstall_cmds=""
-postuninstall_cmds=""
-
-# Commands used to build a loadable module (assumed same as above if empty)
-module_cmds=""
-module_expsym_cmds=""
-
-# Commands to strip libraries.
-old_striplib="strip --strip-debug"
-striplib="strip --strip-unneeded"
-
-# Dependencies to place before the objects being linked to create a
-# shared library.
-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=`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.
-predeps=""
-
-# Dependencies to place after the objects being linked to create a
-# shared library.
-postdeps=""
-
-# The library search path used internally by the compiler when linking
-# a shared library.
-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"
-
-# Command to use when deplibs_check_method == file_magic.
-file_magic_cmd="\$MAGIC_CMD"
-
-# Flag that allows shared libraries with undefined symbols to be built.
-allow_undefined_flag=""
-
-# Flag that forces no undefined symbols.
-no_undefined_flag=""
-
-# Commands used to finish a libtool library installation in a directory.
-finish_cmds="PATH=\\\"\\\$PATH:/sbin\\\" ldconfig -n \$libdir"
-
-# Same as above, but a single script fragment to be evaled but not shown.
-finish_eval=""
-
-# Take the output of nm and produce a listing of raw symbols and C names.
-global_symbol_pipe="sed -n -e 's/^.*[ ]\\([ABCDGIRSTW][ABCDGIRSTW]*\\)[ ][ ]*\\([_A-Za-z][_A-Za-z0-9]*\\)\$/\\1 \\2 \\2/p'"
-
-# Transform the output of nm in a proper C declaration
-global_symbol_to_cdecl="sed -n -e 's/^. .* \\(.*\\)\$/extern int \\1;/p'"
-
-# Transform the output of nm in a C name address pair
-global_symbol_to_c_name_address="sed -n -e 's/^: \\([^ ]*\\) \$/ {\\\"\\1\\\", (lt_ptr) 0},/p' -e 's/^[BCDEGRST] \\([^ ]*\\) \\([^ ]*\\)\$/ {\"\\2\", (lt_ptr) \\&\\2},/p'"
-
-# This is the shared library runtime path variable.
-runpath_var=LD_RUN_PATH
-
-# This is the shared library path variable.
-shlibpath_var=LD_LIBRARY_PATH
-
-# Is shlibpath searched before the hard-coded library search path?
-shlibpath_overrides_runpath=no
-
-# How to hardcode a shared library path into an executable.
-hardcode_action=immediate
-
-# Whether we should hardcode library paths into libraries.
-hardcode_into_libs=yes
-
-# Flag to hardcode $libdir into a binary during linking.
-# This must work even if $libdir does not exist.
-hardcode_libdir_flag_spec="\${wl}--rpath \${wl}\$libdir"
-
-# If ld is used when linking, flag to hardcode $libdir into
-# a binary during linking. This must work even if $libdir does
-# not exist.
-hardcode_libdir_flag_spec_ld=""
-
-# Whether we need a single -rpath flag with a separated argument.
-hardcode_libdir_separator=""
-
-# Set to yes if using DIR/libNAME during linking hardcodes DIR into the
-# resulting binary.
-hardcode_direct=no
-
-# Set to yes if using the -LDIR flag during linking hardcodes DIR into the
-# resulting binary.
-hardcode_minus_L=no
-
-# Set to yes if using SHLIBPATH_VAR=DIR during linking hardcodes DIR into
-# the resulting binary.
-hardcode_shlibpath_var=unsupported
-
-# Set to yes if building a shared library automatically hardcodes DIR into the library
-# and all subsequent libraries and executables linked against it.
-hardcode_automatic=no
-
-# Variables whose values should be saved in libtool wrapper scripts and
-# restored at relink time.
-variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COMPILER_PATH LIBRARY_PATH"
-
-# Whether libtool must link a program against all its dependency libraries.
-link_all_deplibs=unknown
-
-# Compile-time system search path for libraries
-sys_lib_search_path_spec=`echo " /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../ /lib/i386-redhat-linux/3.4.6/ /lib/ /usr/lib/i386-redhat-linux/3.4.6/ /usr/lib/" | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
-
-# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib /usr/lib/qt4/lib "
-
-# Fix the shell variable $srcfile for the compiler.
-fix_srcfile_path=""
-
-# Set to yes if exported symbols are required.
-always_export_symbols=no
-
-# The commands to list exported symbols.
-export_symbols_cmds="\$NM \$libobjs \$convenience | \$global_symbol_pipe | \$SED 's/.* //' | sort | uniq > \$export_symbols"
-
-# The commands to extract the exported symbol list from a shared archive.
-extract_expsyms_cmds=""
-
-# Symbols that should not be listed in the preloaded symbols.
-exclude_expsyms="_GLOBAL_OFFSET_TABLE_"
-
-# Symbols that must always be exported.
-include_expsyms=""
-
-# ### END LIBTOOL TAG CONFIG: F77
-
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index fc2c89e38d8851127b3e2b5f60037510bc8acf59..297ba661b12e244e4905173ee5a389848fa42ff0 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -73,6 +73,7 @@ namespace AMDiS {
userMat += factor * elementMatrix;
}
+
void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
@@ -110,9 +111,10 @@ namespace AMDiS {
ElementMatrix m(nRow, nRow);
smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
+
ElementMatrix tmpMat(nRow, nRow);
tmpMat = m * mat;
- mat = tmpMat;
+ mat = tmpMat;
}
if (firstOrderAssemblerGrdPsi) {
@@ -130,8 +132,17 @@ namespace AMDiS {
ElementMatrix m(nRow, nRow);
smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
+
+#if 0
+ std::cout << "ASM MAT:" << std::endl;
+ std::cout << mat << std::endl;
+ std::cout << "MULT MAT:" << std::endl;
+ std::cout << m << std::endl;
+#endif
+
ElementMatrix tmpMat(nRow, nRow);
- tmpMat = m * mat;
+ //tmpMat = m * mat;
+ tmpMat = mat * m;
mat = tmpMat;
}
@@ -139,6 +150,7 @@ namespace AMDiS {
userMat += factor * elementMatrix;
}
+
void Assembler::calculateElementVector(const ElInfo *elInfo,
ElementVector& userVec,
double factor)
diff --git a/AMDiS/src/ComponentTraverseInfo.h b/AMDiS/src/ComponentTraverseInfo.h
index 4eb511cd61e3e2656f08b025326f9ae3706592cb..4174599244e50029a90013d13d166c619872cb4d 100644
--- a/AMDiS/src/ComponentTraverseInfo.h
+++ b/AMDiS/src/ComponentTraverseInfo.h
@@ -185,22 +185,6 @@ namespace AMDiS {
return vectorComponents[row].getAuxFESpace(0);
}
- void fakeFESpace(const FiniteElemSpace *feOld,
- const FiniteElemSpace *feNew)
- {
- for (int i = 0; i < nComponents; i++) {
- for (int j = 0; j < nComponents; j++) {
- if (matrixComponents[i][j].getAuxFESpace(0) == feOld) {
- matrixComponents[i][j].setAuxFESpace(feNew, 0);
- }
- }
-
- if (vectorComponents[i].getAuxFESpace(0) == feOld) {
- vectorComponents[i].setAuxFESpace(feNew, 0);
- }
- }
- }
-
protected:
int nComponents;
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index fcc84633746208c484027612b498c13ad333fe34..0ee1b5a45a0be2e6d81c0f8e41b82ed221a56f87 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -186,6 +186,21 @@ namespace AMDiS {
}
}
+ using namespace mtl;
+
+#if 0
+ std::cout << "----- PRINT MAT --------" << std::endl;
+ std::cout << elMat << std::endl;
+ std::cout << "rows: ";
+ for (int i = 0; i < rowIndices.size(); i++)
+ std::cout << rowIndices[i] << " ";
+ std::cout << std::endl;
+ std::cout << "cols: ";
+ for (int i = 0; i < colIndices.size(); i++)
+ std::cout << colIndices[i] << " ";
+ std::cout << std::endl;
+#endif
+
for (int i = 0; i < nRow; i++) {
DegreeOfFreedom row = rowIndices[i];
@@ -212,14 +227,17 @@ namespace AMDiS {
}
}
+
double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
{
return matrix[a][b];
}
+
void DOFMatrix::freeDOFContent(int index)
{}
+
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
{
FUNCNAME("DOFMatrix::assemble()");
@@ -235,9 +253,10 @@ namespace AMDiS {
if (factor != 1.0)
elementMatrix *= factor;
- addElementMatrix(elementMatrix, bound, elInfo, NULL);
+ addElementMatrix(elementMatrix, bound, elInfo, NULL);
}
+
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
Operator *op)
{
@@ -254,6 +273,7 @@ namespace AMDiS {
addElementMatrix(elementMatrix, bound, elInfo, NULL);
}
+
void DOFMatrix::assemble(double factor,
ElInfo *rowElInfo, ElInfo *colElInfo,
ElInfo *smallElInfo, ElInfo *largeElInfo,
@@ -272,12 +292,11 @@ namespace AMDiS {
} else {
std::vector<Operator*>::iterator it = operators.begin();
std::vector<double*>::iterator factorIt = operatorFactor.begin();
- for (; it != operators.end(); ++it, ++factorIt) {
+ for (; it != operators.end(); ++it, ++factorIt)
(*it)->getElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo,
elementMatrix,
- *factorIt ? **factorIt : 1.0);
- }
+ *factorIt ? **factorIt : 1.0);
}
if (factor != 1.0)
@@ -286,10 +305,11 @@ namespace AMDiS {
addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);
}
+
void DOFMatrix::assemble2(double factor,
ElInfo *mainElInfo, ElInfo *auxElInfo,
- ElInfo *smallElInfo, ElInfo *largeElInfo,
- const BoundaryType *bound, Operator *op)
+ ElInfo *smallElInfo, ElInfo *largeElInfo,
+ const BoundaryType *bound, Operator *op)
{
FUNCNAME("DOFMatrix::assemble2()");
@@ -419,27 +439,32 @@ namespace AMDiS {
}
}
+
void DOFMatrix::copy(const DOFMatrix& rhs)
{
matrix = rhs.matrix;
}
+
void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
{
- inserter_type &ins= *inserter;
-
+ FUNCNAME("DOFMatrix::removeRowsWithDBC()");
+
+ inserter_type &ins = *inserter;
for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
- ins[*it][*it] = 1.0;
+ ins[*it][*it] = 1.0;
rows->clear();
}
+
int DOFMatrix::memsize()
{
return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
+ matrix.nnz() * sizeof(base_matrix_type::value_type);
}
+
void DOFMatrix::startInsertion(int nnz_per_row)
{
if (inserter) {
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 243ed39f341d15f1b565741cd47c63e96b120141..44a7a05a91e25fe3515bee781e88d925a20c73cd 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -307,8 +307,9 @@ namespace AMDiS {
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
- if (quad && quadFast)
+ if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
+ }
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
@@ -708,16 +709,10 @@ namespace AMDiS {
DOFVector<DegreeOfFreedom>()
{}
+
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof)
- {
- /*
- std::vector<DegreeOfFreedom>::iterator it;
- std::vector<DegreeOfFreedom>::iterator end = vec.end();
- DegreeOfFreedom pos = 0;
- for (it = vec.begin(); it != end; ++it, ++pos)
- if (*it == dof) *it = pos;
- */
- }
+ {}
+
WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
WorldVector<DOFVector<double>*> *res)
@@ -745,6 +740,7 @@ namespace AMDiS {
return r;
}
+
template<>
const double *DOFVectorBase<double>::getVecAtQPs(const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
@@ -765,7 +761,8 @@ namespace AMDiS {
("invalid basis functions");
const BasisFunction *basFcts = feSpace->getBasisFcts();
- int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
+ int nPoints =
+ quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static double *localvec = NULL;
double *result;
@@ -800,6 +797,7 @@ namespace AMDiS {
return const_cast<const double*>(result);
}
+
template<>
void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
const BoundaryType *bound,
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 21c9357f56dab7f556ac8d30caec6360b6dd19e1..ef41f136f907885c1d9abb1825f21ab6011637cc 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -987,9 +987,10 @@ namespace AMDiS {
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
- if (quad && quadFast)
+ if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())
("quad != quadFast->quadrature\n");
+ }
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
index b71621be537a607a81006b76f56540ecc0c9a4ac..ec187a044cbc2b07422a97886ad881504c799b99 100644
--- a/AMDiS/src/Debug.cc
+++ b/AMDiS/src/Debug.cc
@@ -44,12 +44,39 @@ namespace AMDiS {
int myRank = MPI::COMM_WORLD.Get_rank();
if (rank == -1 || myRank == rank) {
DOFVector<double> tmp(feSpace, "tmp");
- VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu", false);
+ VtkWriter::writeFile(tmp, filename + lexical_cast<std::string>(myRank) + ".vtu",
+ false);
}
}
#endif
+ void writeDofIndexMesh(FiniteElemSpace *feSpace)
+ {
+ DOFVector<double> tmp(feSpace, "tmp");
+ DOFIterator<double> it(&tmp, USED_DOFS);
+ for (it.reset(); !it.end(); ++it)
+ *it = it.getDOFIndex();
+ VtkWriter::writeFile(tmp, "dofindex.vtu");
+ }
+
+
+ void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename)
+ {
+ std::map<int, double> vec;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
+
+ while (elInfo) {
+ int index = elInfo->getElement()->getIndex();
+ vec[index] = index;
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ ElementFileWriter::writeFile(vec, feSpace, filename);
+ }
+
+
void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el)
{
// === Get local indices of the given element. ===
@@ -187,22 +214,6 @@ namespace AMDiS {
}
- void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename)
- {
- std::map<int, double> vec;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL);
-
- while (elInfo) {
- int index = elInfo->getElement()->getIndex();
- vec[index] = index;
- elInfo = stack.traverseNext(elInfo);
- }
-
- ElementFileWriter::writeFile(vec, feSpace, filename);
- }
-
-
void writeMatlabMatrix(DOFMatrix &mat, std::string filename)
{
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h
index 14a8cfe78cda16d7c26e0710f64ff8b3b49c902a..2c7927a28a54d6d6eda1e40db0fea52bdbd1c6ad 100644
--- a/AMDiS/src/Debug.h
+++ b/AMDiS/src/Debug.h
@@ -43,9 +43,36 @@ namespace AMDiS {
void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename);
+ /** \brief
+ * Writes a vtu file with the mesh, where all DOFs are set to zero, and only
+ * one given DOF is set to one. This can be used to easily identify DOFs in
+ * a mesh.
+ *
+ * \param[in] rank If set to -1, the vtu files are written on all ranks.
+ * Otherwise, only on the given rank the mesh is written.
+ * \param[in] dof Defines the DOF, which value is set to one in the mesh file.
+ * \param[in] feSpace The FE space to be used.
+ */
void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace);
#endif
-
+
+ /** \brief
+ * Create a vtu file with name 'dofindex.vtu'. All nodes in the mesh are colored
+ * by the global dof index.
+ *
+ * \param[in] feSpace The FE space to be used.
+ */
+ void writeDofIndexMesh(FiniteElemSpace *feSpace);
+
+ /** \brief
+ * Creates a vtu file where all elements in the mesh are colored by the global
+ * element indices.
+ *
+ * \param[in] feSpace The FE space to be used.
+ * \param[in] filename Name of the file.
+ */
+ void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename);
+
void colorDofVectorByLocalElementDofs(DOFVector<double>& vec, Element *el);
bool colorDofVectorByLocalElementDofs(DOFVector<double>& vec,
@@ -60,8 +87,6 @@ namespace AMDiS {
void printMatValuesStatistics(Matrix<DOFMatrix*> *mat);
- void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename);
-
/** \brief
* Creates a text file storing the value of a sparse matrix. Each line of the file
* has three columns:
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index c8f9212ff9b4f869950e52a121cd7911bc64edeb..4038ff857016df73c087325b39d5c649b299d4b0 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -294,10 +294,12 @@ namespace AMDiS {
}
}
- void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat,
- int basisFctsDegree) const
+
+ void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
{
- switch (basisFctsDegree) {
+ FUNCNAME("ElInfo1d::getSubElemCoordsMat()");
+
+ switch (degree) {
case 1:
mat = mat_d1;
@@ -323,14 +325,16 @@ namespace AMDiS {
break;
default:
- ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
+ ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
}
- void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
- int basisFctsDegree) const
+
+ void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
{
- switch (basisFctsDegree) {
+ FUNCNAME("ElInfo1d::getSubElemCoordsMat_so()");
+
+ switch (degree) {
case 1:
mat = test2_val;
@@ -339,7 +343,7 @@ namespace AMDiS {
break;
default:
- ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
+ ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
}
}
diff --git a/AMDiS/src/ElInfo1d.h b/AMDiS/src/ElInfo1d.h
index e3327af1cfe50f31fe53d04f23538092c2cf4366..35c9ff28420d466549e54a289a7392a2a262a712 100644
--- a/AMDiS/src/ElInfo1d.h
+++ b/AMDiS/src/ElInfo1d.h
@@ -62,11 +62,9 @@ namespace AMDiS {
return (i + 1) % 2;
}
- void getSubElemCoordsMat(mtl::dense2D<double>& mat,
- int basisFctsDegree) const;
+ void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
- void getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
- int basisFctsDegree) const;
+ void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
protected:
static double mat_d1_val[2][2];
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index f7f541ff7fbf4750adbc6ca1c206fc965d628333..9a837edc5a99bfc25b706fa09e30f4915fbf1878 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -19,16 +19,31 @@ namespace AMDiS {
{0.0, 0.0, 1.0}};
mtl::dense2D<double> ElInfo2d::mat_d1(mat_d1_val);
+ /*
double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5},
- {0.0, 0.0, 0.5},
- {1.0, 0.0, 0.0}};
+ {0.0, 0.0, 0.5},
+ {1.0, 0.0, 0.0}};
+ */
+
+ double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 0.0, 1.0},
+ {1.0, 0.0, 0.0},
+ {0.5, 0.5, 0.0}};
+
mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);
+ /*
double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5},
- {1.0, 0.0, 0.5},
- {0.0, 1.0, 0.0}};
+ {1.0, 0.0, 0.5},
+ {0.0, 1.0, 0.0}};
+ */
+
+ double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 1.0, 0.0},
+ {0.0, 0.0, 1.0},
+ {0.5, 0.5, 0.0}};
+
mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);
+
ElInfo2d::ElInfo2d(Mesh *aMesh)
: ElInfo(aMesh)
{
@@ -37,6 +52,7 @@ namespace AMDiS {
normal = new WorldVector<double>;
}
+
ElInfo2d::~ElInfo2d()
{
delete e1;
@@ -44,6 +60,7 @@ namespace AMDiS {
delete normal;
}
+
void ElInfo2d::fillMacroInfo(const MacroElement * mel)
{
FUNCNAME("ElInfo::fillMacroInfo()");
@@ -571,7 +588,8 @@ namespace AMDiS {
return adet;
}
- const int ElInfo2d::worldToCoord(const WorldVector<double>& xy,
+
+ const int ElInfo2d::worldToCoord(const WorldVector<double>& xy,
DimVec<double>* lambda) const
{
FUNCNAME("ElInfo::worldToCoord()");
@@ -654,6 +672,7 @@ namespace AMDiS {
return det;
}
+
/****************************************************************************/
/* calculate the normal of the element for dim of world = dim + 1 */
/* return the absulute value of the determinant from the */
@@ -680,4 +699,43 @@ namespace AMDiS {
return det;
}
+
+ void ElInfo2d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
+ {
+ FUNCNAME("ElInfo2d::getSubElemCoordsMat()");
+
+ using namespace mtl;
+
+ switch (degree) {
+ case 1:
+ {
+ mat = mat_d1;
+ dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+
+ for (int i = 0; i < refinementPathLength; i++) {
+ if (refinementPath & (1 << i)) {
+ tmpMat = mat_d1_right * mat;
+ mat = tmpMat;
+ // mat *= mat_d1_right;
+ } else {
+ tmpMat = mat_d1_left * mat;
+ mat = tmpMat;
+
+ // mat *= mat_d1_left;
+ }
+ }
+ }
+ break;
+ default:
+ ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
+ }
+ }
+
+
+ void ElInfo2d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
+ {
+ FUNCNAME("ElInfo2d::getSubElemCoordsMat_so()");
+
+ ERROR_EXIT("Not yet implemented!\n");
+ }
}
diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h
index e0141092979794e90a4160a1921ba8a5f9ab6dc3..4d530ab4b5e5e25617424000624212a439bdf955 100644
--- a/AMDiS/src/ElInfo2d.h
+++ b/AMDiS/src/ElInfo2d.h
@@ -58,6 +58,10 @@ namespace AMDiS {
/// 2-dimensional realisation of ElInfo's getElementNormal method.
double getElementNormal(WorldVector<double> &normal) const;
+ void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
+
+ void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
+
protected:
/// Temp vectors for function \ref calcGrdLambda.
WorldVector<double> *e1, *e2, *normal;
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 44d495d23547bf3e8d8c63f597612381ed4885e6..9d72a1056a2c7177a317d21452809e38b7a8d3f4 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -590,67 +590,20 @@ namespace AMDiS {
}
}
-#if 0
- void ElInfo3d::getRefSimplexCoords(const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords) const
+
+ void ElInfo3d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
{
-// (*coords)[0][0] = 1.0;
-// (*coords)[1][0] = 0.0;
-// (*coords)[2][0] = 0.0;
-// (*coords)[3][0] = 0.0;
-
-// (*coords)[0][1] = 0.0;
-// (*coords)[1][1] = 1.0;
-// (*coords)[2][1] = 0.0;
-// (*coords)[3][1] = 0.0;
-
-// (*coords)[0][2] = 0.0;
-// (*coords)[1][2] = 0.0;
-// (*coords)[2][2] = 1.0;
-// (*coords)[3][2] = 0.0;
-
-// (*coords)[0][3] = 0.0;
-// (*coords)[1][3] = 0.0;
-// (*coords)[2][3] = 0.0;
-// (*coords)[3][3] = 1.0;
+ FUNCNAME("ElInfo3d::getSubElemCoordsMat()");
+
+ ERROR_EXIT("Not yet implemented!\n");
}
- void ElInfo3d::getSubElementCoords(const BasisFunction *basisFcts,
- int iChild,
- mtl::dense2D<double>& coords) const
+
+ void ElInfo3d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
{
- FUNCNAME("ElInfo3d::getSubElementCoords()");
-
- ERROR_EXIT("Not yet\n");
-
-// double c0 = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
-// double c1 = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
-// double c2 = ((*coords)[2][0] + (*coords)[2][1]) * 0.5;
-// double c3 = ((*coords)[3][0] + (*coords)[3][1]) * 0.5;
-
-// if (iChild == 0) {
-// (*coords)[0][1] = (*coords)[0][0];
-// (*coords)[1][1] = (*coords)[1][0];
-// (*coords)[2][1] = (*coords)[2][0];
-
-// (*coords)[0][0] = (*coords)[0][2];
-// (*coords)[1][0] = (*coords)[1][2];
-// (*coords)[2][0] = (*coords)[2][2];
-// } else {
-// (*coords)[0][1] = (*coords)[0][2];
-// (*coords)[1][1] = (*coords)[1][2];
-// (*coords)[2][1] = (*coords)[2][2];
-
-// (*coords)[0][0] = (*coords)[0][1];
-// (*coords)[1][0] = (*coords)[1][1];
-// (*coords)[2][0] = (*coords)[2][1];
-// }
-
-// (*coords)[0][3] = c0;
-// (*coords)[1][3] = c1;
-// (*coords)[2][3] = c2;
-// (*coords)[3][3] = c3;
+ FUNCNAME("ElInfo3d::getSubElemCoordsMat_so()");
+
+ ERROR_EXIT("Not yet implemented!\n");
}
-#endif
}
diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h
index 50ae715420ed2c6b00a72dadbac174aec9fabdca..f7364b39c27983f679d7f6e1a0e4a6896a578fa5 100644
--- a/AMDiS/src/ElInfo3d.h
+++ b/AMDiS/src/ElInfo3d.h
@@ -80,6 +80,10 @@ namespace AMDiS {
orientation = o;
}
+ void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
+
+ void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
+
protected:
/** \brief
* +/- 1: sign of the determinant of the transformation to the reference
diff --git a/AMDiS/src/ElementFileWriter.cc b/AMDiS/src/ElementFileWriter.cc
index 0ba5dd5654b6555402a9627a5efd7331338b6926..f312cc96cc99c6affa79d017d1537dba8387bba7 100644
--- a/AMDiS/src/ElementFileWriter.cc
+++ b/AMDiS/src/ElementFileWriter.cc
@@ -91,6 +91,7 @@ namespace AMDiS {
}
}
+
void ElementFileWriter::writeFile(std::map<int, double> &vec,
const FiniteElemSpace *feSpace,
std::string filename)
@@ -179,6 +180,7 @@ namespace AMDiS {
fout.close();
}
+
void ElementFileWriter::writeMeshDatValues(std::string filename, double time)
{
FUNCNAME("ElementFileWriter::writeMeshDatValues()");
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 847b8cd3bda1daa1507e7f0217fa8abae455a81c..9d10f8cc07aa738d39d08b49d196cfef15ca49aa 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -2759,23 +2759,20 @@ namespace AMDiS {
TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n");
TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n");
- Element *el;
int node, n0;
- const DegreeOfFreedom *cdof;
- DegreeOfFreedom pdof[4], dof9;
- const DOFAdmin *admin;
+ DegreeOfFreedom pdof[10], dof9;
if (n < 1)
return;
- el = list->getElement(0);
- admin = drv->getFESpace()->getAdmin();
+ Element *el = list->getElement(0);
+ const DOFAdmin *admin = drv->getFESpace()->getAdmin();
basFct->getLocalIndices(el, admin, pdof);
/****************************************************************************/
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[pdof[0]] +=
(0.0625*(-(*drv)[cdof[2]] + (*drv)[cdof[6]] - (*drv)[cdof[9]])
@@ -3222,15 +3219,12 @@ namespace AMDiS {
TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n");
TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n");
- const Element *el;
- DegreeOfFreedom pdof[5];
- const DegreeOfFreedom *cdof;
- const DOFAdmin *admin;
-
- if (n < 1) return;
- el = list->getElement(0);
+ if (n < 1)
+ return;
- admin = drv->getFESpace()->getAdmin();
+ DegreeOfFreedom pdof[15];
+ const Element *el = list->getElement(0);
+ const DOFAdmin *admin = drv->getFESpace()->getAdmin();
basFct->getLocalIndices(el, admin, pdof);
@@ -3238,7 +3232,7 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[pdof[0]] +=
(0.2734375*(*drv)[cdof[3]]
@@ -3552,17 +3546,13 @@ namespace AMDiS {
TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n");
TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n");
- DegreeOfFreedom pd[32];
- const DegreeOfFreedom *cd;
- const Element *el;
- int i, typ, lr_set;
- const DOFAdmin *admin;
-
- if (n < 1) return;
- el = list->getElement(0);
- typ = list->getType(0);
+ if (n < 1)
+ return;
- admin = drv->getFESpace()->getAdmin();
+ DegreeOfFreedom pd[35];
+ const Element *el = list->getElement(0);
+ int typ = list->getType(0);
+ const DOFAdmin *admin = drv->getFESpace()->getAdmin();
basFct->getLocalIndices(el, admin, pd);
@@ -3570,7 +3560,7 @@ namespace AMDiS {
/* values on child[0] */
/****************************************************************************/
- cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
+ const DegreeOfFreedom *cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[pd[0]] +=
(0.2734375*(*drv)[cd[10]]
@@ -3801,12 +3791,12 @@ namespace AMDiS {
/* adjust neighbour values */
/****************************************************************************/
- for (i = 1; i < n; i++) {
+ for (int i = 1; i < n; i++) {
el = list->getElement(i);
typ = list->getType(i);
basFct->getLocalIndices(el, admin, pd);
- lr_set = 0;
+ int lr_set = 0;
if (list->getNeighbourElement(i,0) && list->getNeighbourNr(i,0) < i)
lr_set = 1;
@@ -4556,7 +4546,7 @@ namespace AMDiS {
Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFESpace()->getAdmin();
- DegreeOfFreedom pdof[5];
+ DegreeOfFreedom pdof[15];
basFct->getLocalIndices(el, admin, pdof);
// values on child[0]
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index 1a4cb610ab4bbedf99191f3951074b2db00dc391..7233ee557f21c4fe04befe945172a14dfa387d20 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -1337,12 +1337,13 @@ namespace AMDiS {
break;
}
- if (mesh->getNumberOfDOFs(EDGE))
+ if (mesh->getNumberOfDOFs(EDGE)) {
TEST_EXIT(nei->index > mel_index)
("neighbour index %d < element index %d\n", nei->getIndex(), mel_index);
+ }
- if (!mesh->getNumberOfDOFs(EDGE) && nei->getIndex() < mel_index) return false;
-
+ if (!mesh->getNumberOfDOFs(EDGE) && nei->getIndex() < mel_index)
+ return false;
edge_no = Tetrahedron::edgeOfDOFs[j][k];
@@ -1417,10 +1418,11 @@ namespace AMDiS {
if (j == 4 || k == 4)
return false;
- if (mesh->getNumberOfDOFs(EDGE))
+ if (mesh->getNumberOfDOFs(EDGE)) {
TEST_EXIT(nei->getIndex() > mel_index)
("neighbour index %d < element index %d\n", nei->getIndex(),
mel_index);
+ }
if (nei->getIndex() < mel_index)
return false;
diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc
index 032418a77c9dc472ee1f6dbad4fb719e2b90d5fc..c2643ec3fd38755027214aeb034f987c59f86a1c 100644
--- a/AMDiS/src/MeshStructure.cc
+++ b/AMDiS/src/MeshStructure.cc
@@ -237,8 +237,9 @@ namespace AMDiS {
Element *element = elInfo->getElement();
- if (isLeafElement())
+ if (isLeafElement()) {
TEST_EXIT_DBG(element->isLeaf())("mesh finer than code\n");
+ }
if (element->isLeaf() && !isLeafElement()) {
MeshStructure *structure = new MeshStructure();
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 239ce2ee5d2218dd719a415669ca384f7a788afd..7d8136d79c1135522b13c07b5b6cc89df006ccf8 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -251,6 +251,7 @@ namespace AMDiS {
uhOld = uhOld_;
}
+
void Operator::getElementMatrix(const ElInfo *elInfo,
ElementMatrix& userMat,
double factor)
@@ -263,6 +264,7 @@ namespace AMDiS {
assembler[myRank]->calculateElementMatrix(elInfo, userMat, factor);
}
+
void Operator::getElementMatrix(const ElInfo *rowElInfo, const ElInfo *colElInfo,
const ElInfo *smallElInfo, const ElInfo *largeElInfo,
ElementMatrix& userMat,
@@ -278,6 +280,7 @@ namespace AMDiS {
userMat, factor);
}
+
void Operator::getElementVector(const ElInfo *elInfo,
ElementVector& userVec,
double factor)
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 821767500eaa5d625986a19af3db19affdf7401a..f64dfa51fde6774004bf7db3d701b4151d37b59b 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -760,20 +760,14 @@ namespace AMDiS {
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
bool symm = false);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -781,9 +775,7 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -791,9 +783,7 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
- */
+ /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
AbstractFunction<WorldMatrix<double>, WorldVector<double> > *f;
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;
@@ -804,9 +794,7 @@ namespace AMDiS {
*/
WorldVector<double>* gradAtQPs;
- /** \brief
- * True, if \ref f provides always a symmetric WorldMatrix<double>.
- */
+ /// True, if \ref f provides always a symmetric WorldMatrix<double>.
bool symmetric;
};
@@ -825,20 +813,14 @@ namespace AMDiS {
FctGradient_SOT(DOFVectorBase<double> *dv,
AbstractFunction<double, WorldVector<double> > *af);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -846,9 +828,7 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -856,9 +836,7 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function wich maps \ref gradAtQPs to a double.
- */
+ /// Function wich maps \ref gradAtQPs to a double.
AbstractFunction<double, WorldVector<double> > *f;
/** \brief
@@ -884,20 +862,14 @@ namespace AMDiS {
VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<double, double, WorldVector<double> > *af);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -905,9 +877,7 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -915,9 +885,7 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function wich maps \ref gradAtQPs to a double.
- */
+ /// Function wich maps \ref gradAtQPs to a double.
BinaryAbstractFunction<double, double, WorldVector<double> > *f;
/** \brief
@@ -945,20 +913,14 @@ namespace AMDiS {
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
bool symm = false);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -966,9 +928,7 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -976,9 +936,7 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function wich maps \ref gradAtQPs to a double.
- */
+ /// Function wich maps \ref gradAtQPs to a double.
BinaryAbstractFunction<WorldMatrix<double>,
double, WorldVector<double> > *f;
@@ -1010,20 +968,14 @@ namespace AMDiS {
TertiaryAbstractFunction<double, double,
WorldVector<double>, WorldVector<double> > *af);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1031,9 +983,7 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -1041,9 +991,7 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function wich maps \ref gradAtQPs to a double.
- */
+ /// Function wich maps \ref gradAtQPs to a double.
TertiaryAbstractFunction<double, double, WorldVector<double>,
WorldVector<double> > *f;
@@ -1072,15 +1020,11 @@ namespace AMDiS {
VecAndCoordsAtQP_SOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<double, double, WorldVector<double> > *af);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::eval().
- */
+ /// Implements SecondOrderTerm::eval().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
void eval(int nPoints,
@@ -1095,22 +1039,15 @@ namespace AMDiS {
WorldVector<double> *result) const;
protected:
- /** \brief
- * DOFVector to be evaluated at quadrature points.
- */
+ /// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec;
- /** \brief
- * Pointer to an array containing the DOFVector evaluated at quadrature
- * points.
- */
+ /// Pointer to an array containing the DOFVector evaluated at quadrature points.
double* vecAtQPs;
WorldVector<double>* coordsAtQPs;
- /** \brief
- * Function evaluated at \ref vecAtQPs.
- */
+ /// Function evaluated at \ref vecAtQPs.
BinaryAbstractFunction<double, double, WorldVector<double> > *f;
};
@@ -1134,20 +1071,14 @@ namespace AMDiS {
WorldMatrix<double> > *divAf,
bool symm = false);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1155,9 +1086,7 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -1165,9 +1094,7 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
- */
+ /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
BinaryAbstractFunction<WorldMatrix<double>,
WorldVector<double>, WorldVector<double> > *f;
@@ -1181,9 +1108,7 @@ namespace AMDiS {
WorldVector<double>* coordsAtQPs;
- /** \brief
- * True, if \ref f provides always a symmetric WorldMatrix<double>.
- */
+ /// True, if \ref f provides always a symmetric WorldMatrix<double>.
bool symmetric;
};
@@ -1306,21 +1231,15 @@ namespace AMDiS {
WorldMatrix<double> > *divFct,
bool symmetric);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo*,
SubAssembler* ,
Quadrature *quad= NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1328,9 +1247,7 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
@@ -1459,9 +1376,7 @@ namespace AMDiS {
: FirstOrderTerm(0), factor(fptr)
{}
- /** \brief
- * Implements FirstOrderTerm::getLb().
- */
+ /// Implements FirstOrderTerm::getLb().
inline void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
@@ -1471,9 +1386,7 @@ namespace AMDiS {
}
private:
- /** \brief
- * Pointer to the factor.
- */
+ /// Pointer to the factor.
double *factor;
};
@@ -1595,20 +1508,14 @@ namespace AMDiS {
: FirstOrderTerm(af->getDegree()), g(af)
{}
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements FistOrderTerm::getLb().
- */
+ /// Implements FistOrderTerm::getLb().
void getLb(const ElInfo *elInfo, int nPoints,VectorOfFixVecs<DimVec<double> >&Lb) const;
- /** \brief
- * Implements FirstOrderTerm::eval().
- */
+ /// Implements FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1617,14 +1524,10 @@ namespace AMDiS {
double factor) const;
protected:
- /** \brief
- * Stores coordinates at quadrature points. Set in \ref initElement().
- */
+ /// Stores coordinates at quadrature points. Set in \ref initElement().
WorldVector<double>* coordsAtQPs;
- /** \brief
- * Function avaluated at world coordinates.
- */
+ /// Function avaluated at world coordinates.
AbstractFunction<double, WorldVector<double> > *g;
};
@@ -1643,21 +1546,15 @@ namespace AMDiS {
: FirstOrderTerm(af->getDegree()), g(af), b(wv)
{}
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements FistOrderTerm::getLb().
- */
+ /// Implements FistOrderTerm::getLb().
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const;
- /** \brief
- * Implements FirstOrderTerm::eval().
- */
+ /// Implements FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1666,19 +1563,13 @@ namespace AMDiS {
double factor) const;
protected:
- /** \brief
- * Stores coordinates at quadrature points. Set in \ref initElement().
- */
+ /// Stores coordinates at quadrature points. Set in \ref initElement().
WorldVector<double>* coordsAtQPs;
- /** \brief
- * Function evaluated at world coordinates.
- */
+ /// Function evaluated at world coordinates.
AbstractFunction<double, WorldVector<double> > *g;
- /** \brief
- * Coefficient vector.
- */
+ /// Coefficient vector.
WorldVector<double> b;
};
@@ -1695,21 +1586,15 @@ namespace AMDiS {
VectorGradient_FOT(DOFVectorBase<double> *dv,
AbstractFunction<WorldVector<double>, WorldVector<double> > *af);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements FirstOrderTerm::getLb().
- */
+ /// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const;
- /** \brief
- * Implements FirstOrderTerm::eval().
- */
+ /// Implements FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1720,14 +1605,10 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function for b.
- */
+ /// Function for b.
AbstractFunction<WorldVector<double>, WorldVector<double> > *f;
- /** \brief
- * Gradient of v at quadrature points.
- */
+ /// Gradient of v at quadrature points.
WorldVector<double> *gradAtQPs;
};
@@ -1744,20 +1625,15 @@ namespace AMDiS {
VectorFct_FOT(DOFVectorBase<double> *dv,
AbstractFunction<WorldVector<double>, double> *fct);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements FirstOrderTerm::getLb().
- */
- void getLb(const ElInfo *elInfo, int qPoint, VectorOfFixVecs<DimVec<double> >& Lb) const;
+ /// Implements FirstOrderTerm::getLb().
+ void getLb(const ElInfo *elInfo, int qPoint,
+ VectorOfFixVecs<DimVec<double> >& Lb) const;
- /** \brief
- * Implements FirstOrderTerm::eval().
- */
+ /// Implements FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1766,19 +1642,13 @@ namespace AMDiS {
double factor) const;
protected:
- /** \brief
- * DOFVector to be evaluated at quadrature points.
- */
+ /// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec;
- /** \brief
- * Vector v at quadrature points.
- */
+ /// Vector v at quadrature points.
double *vecAtQPs;
- /** \brief
- * Function for b.
- */
+ /// Function for b.
AbstractFunction<WorldVector<double>, double> *vecFct;
};
@@ -1796,20 +1666,15 @@ namespace AMDiS {
: FirstOrderTerm(af->getDegree()), g(af)
{}
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements FistOrderTerm::getLb().
- */
- void getLb(const ElInfo *elInfo, int nPoints,VectorOfFixVecs<DimVec<double> >&Lb) const;
+ /// Implements FistOrderTerm::getLb().
+ void getLb(const ElInfo *elInfo, int nPoints,
+ VectorOfFixVecs<DimVec<double> >&Lb) const;
- /** \brief
- * Implements FirstOrderTerm::eval().
- */
+ /// Implements FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1818,14 +1683,10 @@ namespace AMDiS {
double factor) const;
protected:
- /** \brief
- * Stores coordinates at quadrature points. Set in \ref initElement().
- */
+ /// Stores coordinates at quadrature points. Set in \ref initElement().
WorldVector<double>* coordsAtQPs;
- /** \brief
- * Function avaluated at world coordinates.
- */
+ /// Function avaluated at world coordinates.
AbstractFunction<WorldVector<double>, WorldVector<double> > *g;
};
@@ -1866,22 +1727,16 @@ namespace AMDiS {
double factor) const;
protected:
- /** \brief
- * DOFVector to be evaluated at quadrature points.
- */
+ /// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
- /** \brief
- * Vector v at quadrature points.
- */
+ /// Vector v at quadrature points.
double *vecAtQPs;
WorldVector<double> *gradAtQPs;
- /** \brief
- * Function for b.
- */
+ /// Function for b.
BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *vecFct;
};
@@ -2121,21 +1976,14 @@ namespace AMDiS {
std::vector<double>,
std::vector<WorldVector<double> > > *f);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
- void initElement(const ElInfo*, SubAssembler*,
- Quadrature *quad= NULL);
+ /// Implementation of \ref OperatorTerm::initElement().
+ void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL);
- /** \brief
- * Implements FirstOrderTerm::getLb().
- */
+ /// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& result) const;
- /** \brief
- * Implenetation of FirstOrderTerm::eval().
- */
+ /// Implenetation of FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -2179,16 +2027,11 @@ namespace AMDiS {
/// Constructor.
ZeroOrderTerm(int deg) : OperatorTerm(deg) {}
- /** \brief
- * Destructor.
- */
+ /// Destructor.
virtual ~ZeroOrderTerm() {}
- /** \brief
- * Evaluates \f$ c \f$
- */
- virtual void getC(const ElInfo *elInfo,
- int nPoints,
+ /// Evaluates \f$ c \f$
+ virtual void getC(const ElInfo *elInfo, int nPoints,
std::vector<double> &C) const = 0;
};
@@ -2294,14 +2137,10 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements ZeroOrderTerm::getC().
- */
+ /// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
- /** \brief
- * Implements ZeroOrderTerm::eval().
- */
+ /// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -2310,21 +2149,15 @@ namespace AMDiS {
double fac) const;
protected:
- /** \brief
- * DOFVectorBase to be evaluated at quadrature points.
- */
+ /// DOFVectorBase to be evaluated at quadrature points.
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
- /** \brief
- * Vector v at quadrature points.
- */
+ /// Vector v at quadrature points.
double *vecAtQPs1;
double *vecAtQPs2;
- /** \brief
- * Function for c.
- */
+ /// Function for c.
AbstractFunction<double, double> *f1;
AbstractFunction<double, double> *f2;
};
@@ -2395,20 +2228,14 @@ namespace AMDiS {
DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *f);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements ZeroOrderTerm::getC().
- */
+ /// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
- /** \brief
- * Implements ZeroOrderTerm::eval().
- */
+ /// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -2417,23 +2244,13 @@ namespace AMDiS {
double fac) const;
protected:
- /** \brief
- * DOFVector to be evaluated at quadrature points.
- */
- DOFVectorBase<double>* vec1;
- DOFVectorBase<double>* vec2;
- DOFVectorBase<double>* vec3;
+ /// DOFVectors to be evaluated at quadrature points.
+ DOFVectorBase<double> *vec1, *vec2, *vec3;
- /** \brief
- * Vector v at quadrature points.
- */
- double *vecAtQPs1;
- double *vecAtQPs2;
- double *vecAtQPs3;
+ /// Vectors at quadrature points.
+ double *vecAtQPs1, *vecAtQPs2, *vecAtQPs3;
- /** \brief
- * Function for c.
- */
+ /// Function for c.
TertiaryAbstractFunction<double, double, double, double> *f;
};
@@ -2451,20 +2268,14 @@ namespace AMDiS {
FctGradientCoords_ZOT(DOFVectorBase<double> *dv,
BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> > *f);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getC().
- */
+ /// Implements SecondOrderTerm::getC().
void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C) const;
- /** \brief
- * Implements ZeroOrderTerm::eval().
- */
+ /// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -2475,9 +2286,7 @@ namespace AMDiS {
protected:
DOFVectorBase<double>* vec;
- /** \brief
- * Function wich maps \ref gradAtQPs to a double.
- */
+ /// Function wich maps \ref gradAtQPs to a double.
BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> > *f;
/** \brief
@@ -3191,14 +3000,11 @@ namespace AMDiS {
double fac) const;
protected:
- DOFVectorBase<double>* vec1;
- DOFVectorBase<double>* vec2;
+ DOFVectorBase<double> *vec1, *vec2;
- double *vecAtQPs1;
- double *vecAtQPs2;
+ double *vecAtQPs1, *vecAtQPs2;
- WorldVector<double> *gradAtQPs1;
- WorldVector<double> *gradAtQPs2;
+ WorldVector<double> *gradAtQPs1, *gradAtQPs2;
QuartAbstractFunction<double, double, double, WorldVector<double>,WorldVector<double> > *f;
};
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 227495d830d725a3702175fff88c74c5297e3eec..4699b657aaf5d205fdcff76a51a4e9913e27c14a 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -363,6 +363,7 @@ namespace AMDiS {
solver->initParameters();
}
+
void ProblemVec::createEstimator()
{
FUNCNAME("ProblemVec::createEstimator()");
@@ -385,23 +386,29 @@ namespace AMDiS {
if (estimatorCreator) {
estimatorCreator->setName(estName);
estimatorCreator->setRow(i);
- if (estimatorType == "recovery") {
+ if (estimatorType == "recovery")
dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->
- setSolution(solution->getDOFVector(i));
- }
+ setSolution(solution->getDOFVector(i));
estimator[i] = estimatorCreator->create();
}
if (estimator[i]) {
- for (int j = 0; j < nComponents; j++)
- estimator[i]->addSystem((*systemMatrix)[i][j],
- solution->getDOFVector(j),
- rhs->getDOFVector(j));
+ MSG("Estimator fake to test different meshes ...\n");
+ WAIT_REALLY;
+ estimator[i]->addSystem((*systemMatrix)[i][i],
+ solution->getDOFVector(i),
+ rhs->getDOFVector(i));
+
+ // for (int j = 0; j < nComponents; j++)
+ // estimator[i]->addSystem((*systemMatrix)[i][j],
+ // solution->getDOFVector(j),
+ // rhs->getDOFVector(j));
}
}
}
+
void ProblemVec::createMarker()
{
FUNCNAME("ProblemVec::createMarker()");
@@ -513,6 +520,7 @@ namespace AMDiS {
adaptInfo->setSolverResidual(solver->getResidual());
}
+
void ProblemVec::estimate(AdaptInfo *adaptInfo)
{
FUNCNAME("ProblemVec::estimate()");
@@ -547,10 +555,10 @@ namespace AMDiS {
#else
INFO(info, 8)("estimation of the error needed %.5f seconds\n",
TIME_USED(first, clock()));
-
#endif
}
+
Flag ProblemVec::markElements(AdaptInfo *adaptInfo)
{
FUNCNAME("ProblemVec::markElements()");
@@ -658,6 +666,8 @@ namespace AMDiS {
for (int j = 0; j < nComponents; j++) {
+ std::cout << "-------" << i << " " << j << "----------------" << std::endl;
+
// 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
@@ -707,6 +717,7 @@ namespace AMDiS {
// The simplest case: either the right hand side has no operaters, no aux
// fe spaces, or all aux fe spaces are equal to the row and col fe space.
+ // std::cout << "ASM 1" << std::endl;
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
@@ -717,28 +728,31 @@ namespace AMDiS {
// Row fe space and col fe space are both equal, but right hand side has at
// least one another aux fe space.
- assembleOnOneMesh(componentSpaces[i],
- assembleFlag,
- assembleMatrix ? matrix : NULL,
- ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
-
- assembleOnDifMeshes2(componentSpaces[i],
- traverseInfo.getAuxFESpace(i, j),
- assembleFlag,
- NULL,
- ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
-
+ // std::cout << "ASM 2" << std::endl;
+
+ if (assembleMatrix)
+ assembleOnOneMesh(componentSpaces[i],
+ assembleFlag,
+ assembleMatrix ? matrix : NULL,
+ ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
+
+ if (i == j && asmVector)
+ assembleOnDifMeshes2(componentSpaces[i],
+ traverseInfo.getAuxFESpace(i),
+ assembleFlag,
+ NULL,
+ rhs->getDOFVector(i));
} else {
ERROR_EXIT("Possible? If yes, not yet implemented!\n");
}
} else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {
-
+ // std::cout << "ASM 4" << std::endl;
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
-
+ // std::cout << "ASM 5" << std::endl;
assembleOnDifMeshes2(componentSpaces[i],
traverseInfo.getAuxFESpace(i, j),
assembleFlag,
@@ -747,7 +761,7 @@ namespace AMDiS {
} else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_NO_AUX ||
traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_WITH_AUX) {
-
+ // std::cout << "ASM 6" << std::endl;
assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
assembleFlag,
assembleMatrix ? matrix : NULL,
@@ -838,11 +852,13 @@ namespace AMDiS {
#endif
}
+
void ProblemVec::writeFiles(AdaptInfo &adaptInfo, bool force)
{
writeFiles(&adaptInfo, force);
}
+
void ProblemVec::interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct)
{
FUNCNAME("ProblemVec::interpolInitialSolution()");
@@ -850,10 +866,14 @@ namespace AMDiS {
solution->interpol(fct);
}
+
void ProblemVec::addMatrixOperator(Operator *op, int i, int j,
double *factor, double *estFactor)
{
FUNCNAME("ProblemVec::addMatrixOperator()");
+
+ TEST_EXIT(!boundaryConditionSet)
+ ("Do not add operators after boundary conditions were set!\n");
if (!(*systemMatrix)[i][j]) {
TEST_EXIT(i != j)("should have been created already\n");
@@ -884,17 +904,22 @@ namespace AMDiS {
opFlags[op].setFlag(Operator::MATRIX_OPERATOR);
}
+
void ProblemVec::addMatrixOperator(Operator &op, int i, int j,
double *factor, double *estFactor)
{
addMatrixOperator(&op, i, j, factor, estFactor);
}
+
void ProblemVec::addVectorOperator(Operator *op, int i,
double *factor, double *estFactor)
{
FUNCNAME("ProblemVec::addVectorOperator()");
+ TEST_EXIT(!boundaryConditionSet)
+ ("Do not add operators after boundary conditions were set!\n");
+
rhs->getDOFVector(i)->addOperator(op, factor, estFactor);
traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces());
@@ -911,28 +936,34 @@ namespace AMDiS {
opFlags[op].setFlag(Operator::VECTOR_OPERATOR);
}
+
void ProblemVec::addVectorOperator(Operator &op, int i,
double *factor, double *estFactor)
{
addVectorOperator(&op, i, factor, estFactor);
}
+
void ProblemVec::addDirichletBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> >* b)
{
FUNCNAME("ProblemVec::addDirichletBC()");
+ boundaryConditionSet = true;
+
DirichletBC *dirichletApply =
new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], true);
DirichletBC *dirichletNotApply =
new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], false);
- for (int i = 0; i < nComponents; i++)
- if (systemMatrix && (*systemMatrix)[row][i])
+ for (int i = 0; i < nComponents; i++) {
+ if (systemMatrix && (*systemMatrix)[row][i]) {
if (i == col)
(*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletApply);
else
- (*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletNotApply);
+ (*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletNotApply);
+ }
+ }
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
@@ -940,11 +971,14 @@ namespace AMDiS {
solution->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
}
+
void ProblemVec::addNeumannBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n)
{
FUNCNAME("ProblemVec::addNeumannBC()");
+ boundaryConditionSet = true;
+
NeumannBC *neumann =
new NeumannBC(type, n, componentSpaces[row], componentSpaces[col]);
@@ -952,12 +986,15 @@ namespace AMDiS {
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
}
+
void ProblemVec::addRobinBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n,
AbstractFunction<double, WorldVector<double> > *r)
{
FUNCNAME("ProblemVec::addRobinBC()");
+ boundaryConditionSet = true;
+
RobinBC *robin =
new RobinBC(type, n, r, componentSpaces[row], componentSpaces[col]);
@@ -967,10 +1004,13 @@ namespace AMDiS {
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
}
+
void ProblemVec::addPeriodicBC(BoundaryType type, int row, int col)
{
FUNCNAME("ProblemVec::addPeriodicBC()");
+ boundaryConditionSet = true;
+
FiniteElemSpace *feSpace = componentSpaces[row];
PeriodicBC *periodic = new PeriodicBC(type, feSpace);
@@ -981,6 +1021,7 @@ namespace AMDiS {
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(periodic);
}
+
void ProblemVec::assembleOnOneMesh(FiniteElemSpace *feSpace,
Flag assembleFlag,
DOFMatrix *matrix, DOFVector<double> *vector)
@@ -1122,7 +1163,7 @@ namespace AMDiS {
#else
if (matrix)
- matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
#endif
if (useGetBound)
@@ -1131,12 +1172,15 @@ namespace AMDiS {
} // pragma omp parallel
}
+
void ProblemVec::assembleOnDifMeshes(FiniteElemSpace *rowFeSpace,
FiniteElemSpace *colFeSpace,
Flag assembleFlag,
DOFMatrix *matrix,
DOFVector<double> *vector)
{
+ FUNCNAME("ProblemVec::assembleOnDifMeshes()");
+
const BasisFunction *basisFcts = rowFeSpace->getBasisFcts();
BoundaryType *bound = NULL;
if (useGetBound)
@@ -1157,6 +1201,12 @@ namespace AMDiS {
&rowElInfo, &colElInfo,
&smallElInfo, &largeElInfo);
while (cont) {
+#if 0
+ std::cout << "ROW = " << rowElInfo->getElement()->getIndex() << " "
+ << "COL = " << colElInfo->getElement()->getIndex() << " "
+ << "SMA = " << smallElInfo->getElement()->getIndex() << " "
+ << "LAR = " << largeElInfo->getElement()->getIndex() << std::endl;
+#endif
if (useGetBound)
basisFcts->getBound(rowElInfo, bound);
@@ -1174,6 +1224,9 @@ namespace AMDiS {
&smallElInfo, &largeElInfo);
}
+ if (matrix)
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+
if (useGetBound)
delete [] bound;
}
@@ -1185,6 +1238,11 @@ namespace AMDiS {
DOFMatrix *matrix,
DOFVector<double> *vector)
{
+ FUNCNAME("ProblemVec::assembleOnDifMeshes2()");
+
+ TEST_EXIT(mainFeSpace)("No main FE space!\n");
+ TEST_EXIT(auxFeSpace)("No aux FE space!\n");
+
Mesh *mainMesh = mainFeSpace->getMesh();
Mesh *auxMesh = auxFeSpace->getMesh();
@@ -1219,15 +1277,21 @@ namespace AMDiS {
&smallElInfo, &largeElInfo);
}
+ if (matrix)
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+
if (useGetBound)
delete [] bound;
}
+
void ProblemVec::assembleBoundaryConditions(DOFVector<double> *rhs,
DOFVector<double> *solution,
Mesh *mesh,
Flag assembleFlag)
{
+ FUNCNAME("ProblemVec::assembleBoundaryConditions()");
+
// === Initialization of vectors ===
if (rhs->getBoundaryManager())
@@ -1321,10 +1385,12 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- ElementFileWriter fw(name, this->getFESpace(comp), vec);
- fw.writeFiles(adaptInfo, true);
+ ElementFileWriter::writeFile(vec, this->getFESpace(comp), name);
+// ElementFileWriter fw(name, this->getFESpace(comp), vec);
+// fw.writeFiles(adaptInfo, true);
}
+
void ProblemVec::serialize(std::ostream &out)
{
FUNCNAME("ProblemVec::serialize()");
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index b65c4e6974c5ef21ca9f9d980affcf6de2d934cb..6199f9ace6ce3b1e189eceac0a58951e5788fc7a 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -70,7 +70,8 @@ namespace AMDiS {
useGetBound(true),
info(10),
allowFirstRef(false),
- computeExactError(false)
+ computeExactError(false),
+ boundaryConditionSet(false)
{
GET_PARAMETER(0, name + "->components", "%d", &nComponents);
TEST_EXIT(nComponents > 0)("components not set!\n");
@@ -607,6 +608,13 @@ namespace AMDiS {
* defined by \ref exactSolutionFcts.
*/
bool computeExactError;
+
+ /** \brief
+ * If at least on boundary condition is set, this variable is true. It is used
+ * to ensure that no operators are added after boundary condition were set. If
+ * this would happen, boundary conditions could set wrong on off diagonal matrices.
+ */
+ bool boundaryConditionSet;
std::map<Operator*, std::vector<OperatorPos> > operators;
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index 35a9cee0a1d7c8af2f824e7dbbba8efe50128899..8244dfd71471da85467cdbaa2f914c1d0a91b103 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -34,15 +34,17 @@ namespace AMDiS {
TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n");
}
+
void ResidualEstimator::init(double ts)
{
FUNCNAME("ResidualEstimator::init()");
-
+
timestep = ts;
mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();
nSystems = static_cast<int>(uh.size());
+
TEST_EXIT_DBG(nSystems > 0)("no system set\n");
dim = mesh->getDim();
@@ -85,11 +87,9 @@ namespace AMDiS {
grdUh_qp = NULL;
D2uhqp = NULL;
- TraverseStack stack;
- ElInfo *elInfo = NULL;
-
// clear error indicators and mark elements for jumpRes
- elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
elInfo->getElement()->setEstimation(0.0, row);
elInfo->getElement()->setMark(1);
@@ -109,7 +109,6 @@ namespace AMDiS {
Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_DET |
Mesh::CALL_LEAF_EL;
-
neighInfo = mesh->createNewElInfo();
// prepare date for computing jump residual
@@ -126,6 +125,7 @@ namespace AMDiS {
}
}
+
void ResidualEstimator::exit(bool output)
{
FUNCNAME("ResidualEstimator::exit()");
@@ -187,6 +187,7 @@ namespace AMDiS {
delete neighInfo;
}
+
void ResidualEstimator::estimateElement(ElInfo *elInfo)
{
FUNCNAME("ResidualEstimator::estimateElement()");
@@ -206,23 +207,21 @@ namespace AMDiS {
riq[iq] = 0.0;
for (int system = 0; system < nSystems; system++) {
-
if (matrix[system] == NULL)
continue;
+ DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[system]);
+ DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[system]);
+
// === init assemblers ===
- for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
- itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
- it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
- ++it, ++itfac)
+ for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
+ it != dofMat->getOperatorsEnd(); ++it, ++itfac)
if (*itfac == NULL || **itfac != 0.0)
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
if (C0 > 0.0)
- for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
- it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd();
- ++it)
+ for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it)
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
if (timestep && uhOld[system]) {
@@ -249,12 +248,9 @@ namespace AMDiS {
}
// === Compute element residual. ===
-
if (C0 > 0.0) {
- for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
- itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
- it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
- ++it, ++itfac) {
+ for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
+ it != dofMat->getOperatorsEnd(); ++it, ++itfac) {
if (*itfac == NULL || **itfac != 0.0) {
if (uhQP == NULL && (*it)->zeroOrderTerms()) {
uhQP = new double[nPoints];
@@ -304,12 +300,12 @@ namespace AMDiS {
// === Compute jump residuals. ===
- if (C1 && (dim > 1)) {
+ if (C1 && dim > 1) {
int dow = Global::getGeo(WORLD);
for (int face = 0; face < nNeighbours; face++) {
Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
- if (neigh && neigh->getMark()) {
+ if (neigh && neigh->getMark()) {
int oppV = elInfo->getOppVertex(face);
el->sortFaceIndices(face, &faceIndEl);
@@ -356,7 +352,7 @@ namespace AMDiS {
}
}
}
-
+
if (!periodicCoords) {
for (int i = 0; i < dim; i++) {
int i1 = faceIndEl[i];
@@ -365,7 +361,7 @@ namespace AMDiS {
neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
}
}
-
+
Parametric *parametric = mesh->getParametric();
if (parametric)
neighInfo = parametric->addParametricInfo(neighInfo);
@@ -428,7 +424,7 @@ namespace AMDiS {
}
}
}
-
+
val = 0.0;
for (int iq = 0; iq < nPointsSurface; iq++)
val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
diff --git a/AMDiS/src/SubAssembler.cc b/AMDiS/src/SubAssembler.cc
index 7fe786d9766f72669b870d07ebf3503be02652c6..e780d6e2a6f404f943a09c72b47443f9b2b33341 100644
--- a/AMDiS/src/SubAssembler.cc
+++ b/AMDiS/src/SubAssembler.cc
@@ -71,6 +71,7 @@ namespace AMDiS {
dim = assembler->rowFESpace->getMesh()->getDim();
}
+
FastQuadrature *SubAssembler::updateFastQuadrature(FastQuadrature *quadFast,
const BasisFunction *psi,
Flag updateFlag)
@@ -85,6 +86,7 @@ namespace AMDiS {
return quadFast;
}
+
void SubAssembler::initElement(const ElInfo* smallElInfo,
const ElInfo *largeElInfo,
Quadrature *quad)
@@ -106,14 +108,14 @@ namespace AMDiS {
// calls initElement of each term
std::vector<OperatorTerm*>::iterator it;
for (it = terms[myRank].begin(); it != terms[myRank].end(); ++it) {
- if (largeElInfo == NULL) {
+ if (largeElInfo == NULL)
(*it)->initElement(smallElInfo, this, quad);
- } else {
- (*it)->initElement(smallElInfo, largeElInfo, this, quad);
- }
+ else
+ (*it)->initElement(smallElInfo, largeElInfo, this, quad);
}
}
+
WorldVector<double>* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo,
Quadrature *quad)
{
@@ -147,6 +149,7 @@ namespace AMDiS {
return coordsAtQPs;
}
+
double* SubAssembler::getVectorAtQPs(DOFVectorBase<double>* dv,
const ElInfo* elInfo,
Quadrature *quad)
@@ -166,7 +169,6 @@ namespace AMDiS {
if (!valuesAtQPs[vec])
valuesAtQPs[vec] = new ValuesAtQPs;
-
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
double *values = valuesAtQPs[vec]->values.getValArray();
@@ -204,6 +206,7 @@ namespace AMDiS {
return values;
}
+
double* SubAssembler::getVectorAtQPs(DOFVectorBase<double>* dv,
const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
@@ -221,16 +224,14 @@ namespace AMDiS {
valuesAtQPs[vec] = new ValuesAtQPs;
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
-
double *values = valuesAtQPs[vec]->values.getValArray();
-
valuesAtQPs[vec]->valid = true;
-
vec->getVecAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values);
return values;
}
+
WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* dv,
const ElInfo* elInfo,
Quadrature *quad)
@@ -261,22 +262,20 @@ namespace AMDiS {
(vec->getFESpace() == owner->colFESpace);
if (opt && !quad && sameFESpaces) {
- if (vec->getFESpace()->getBasisFcts() == psi) {
+ if (vec->getFESpace()->getBasisFcts() == psi)
psiFast = updateFastQuadrature(psiFast, psi, INIT_GRD_PHI);
- } else if(vec->getFESpace()->getBasisFcts() == phi) {
+ else if(vec->getFESpace()->getBasisFcts() == phi)
phiFast = updateFastQuadrature(phiFast, phi, INIT_GRD_PHI);
- }
}
// calculate new values
const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();
if (opt && !quad && sameFESpaces) {
- if (psiFast->getBasisFunctions() == basFcts) {
+ if (psiFast->getBasisFunctions() == basFcts)
vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
- } else {
- vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
- }
+ else
+ vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
} else {
vec->getGrdAtQPs(elInfo, localQuad, NULL, values);
}
@@ -286,6 +285,7 @@ namespace AMDiS {
return values;
}
+
WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* dv,
const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
@@ -302,9 +302,8 @@ namespace AMDiS {
Quadrature *localQuad = quad ? quad : quadrature;
- if (!gradientsAtQPs[vec]) {
- gradientsAtQPs[vec] = new GradientsAtQPs;
- }
+ if (!gradientsAtQPs[vec])
+ gradientsAtQPs[vec] = new GradientsAtQPs;
gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());
WorldVector<double> *values = gradientsAtQPs[vec]->values.getValArray();
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index f2d48199749ba4ada2342714e75f7a7b91d00609..8c9cb62f78c84fabeaf38075dabff0252fe86d86 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -32,8 +32,9 @@ namespace AMDiS {
elinfo_stack[0]->setMesh(mesh);
elinfo_stack[1]->setMesh(mesh);
- if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL))
+ if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
TEST_EXIT_DBG(level >= 0)("invalid level: %d\n", level);
+ }
traverse_mel = NULL;
stack_used = 0;
@@ -41,6 +42,7 @@ namespace AMDiS {
return traverseNext(NULL);
}
+
ElInfo* TraverseStack::traverseNext(ElInfo* elinfo_old)
{
FUNCNAME("TraverseStack::traverseNext()");
@@ -86,6 +88,7 @@ namespace AMDiS {
return elinfo;
}
+
void TraverseStack::enlargeTraverseStack()
{
FUNCNAME("TraverseStack::enlargeTraverseStack()");
@@ -197,14 +200,17 @@ namespace AMDiS {
return elinfo_stack[stack_used];
}
+
ElInfo* TraverseStack::traverseLeafElementLevel()
{
FUNCNAME("TraverseStack::traverseLeafElementLevel()");
+
ERROR_EXIT("not yet implemented\n");
return NULL;
}
+
ElInfo* TraverseStack::traverseElementLevel()
{
FUNCNAME("TraverseStack::traverseElementLevel()");
@@ -217,6 +223,7 @@ namespace AMDiS {
return elInfo;
}
+
ElInfo* TraverseStack::traverseMultiGridLevel()
{
FUNCNAME("TraverseStack::traverseMultiGridLevel()");
@@ -289,12 +296,11 @@ namespace AMDiS {
return traverseMultiGridLevel();
}
+
ElInfo* TraverseStack::traverseEveryElementPreorder()
{
FUNCNAME("TraverseStack::traverseEveryElementPreorder()");
- Element *el;
-
if (stack_used == 0) { /* first call */
currentMacro = traverse_mesh->firstMacroElement();
traverse_mel = *currentMacro;
@@ -308,7 +314,7 @@ namespace AMDiS {
return(elinfo_stack[stack_used]);
}
- el = elinfo_stack[stack_used]->getElement();
+ Element *el = elinfo_stack[stack_used]->getElement();
/* go up in tree until we can go down again */
while (stack_used > 0 &&
@@ -353,6 +359,7 @@ namespace AMDiS {
return(elinfo_stack[stack_used]);
}
+
ElInfo* TraverseStack::traverseEveryElementInorder()
{
FUNCNAME("TraverseStack::traverseEveryElementInorder");
@@ -360,12 +367,10 @@ namespace AMDiS {
return NULL;
}
+
ElInfo* TraverseStack::traverseEveryElementPostorder()
{
- FUNCNAME("TraverseStack::traverseEveryElementPostorder");
-
- Element *el;
- int i;
+ FUNCNAME("TraverseStack::traverseEveryElementPostorder()");
if (stack_used == 0) { /* first call */
currentMacro = traverse_mesh->firstMacroElement();
@@ -379,12 +384,11 @@ namespace AMDiS {
//return(elinfo_stack[stack_used]);
} else { /* don't go up on first call */
-
- el = elinfo_stack[stack_used]->getElement();
+ Element *el = elinfo_stack[stack_used]->getElement();
/* go up in tree until we can go down again */ /* postorder!!! */
- while ((stack_used > 0) &&
- ((info_stack[stack_used] >= 3) || (el->getFirstChild()==NULL))) {
+ while (stack_used > 0 &&
+ (info_stack[stack_used] >= 3 || el->getFirstChild() == NULL)) {
stack_used--;
el = elinfo_stack[stack_used]->getElement();
}
@@ -393,7 +397,8 @@ namespace AMDiS {
/* goto next macro element */
if (stack_used < 1) {
currentMacro++;
- if(currentMacro == traverse_mesh->endOfMacroElements()) return NULL;
+ if (currentMacro == traverse_mesh->endOfMacroElements())
+ return NULL;
traverse_mel = *currentMacro;
stack_used = 1;
@@ -405,12 +410,12 @@ namespace AMDiS {
}
/* go down tree */
- while (elinfo_stack[stack_used]->getElement()->getFirstChild()
- && (info_stack[stack_used] < 2)) {
+ while (elinfo_stack[stack_used]->getElement()->getFirstChild() &&
+ info_stack[stack_used] < 2) {
if (stack_used >= stack_size-1)
enlargeTraverseStack();
- i = info_stack[stack_used];
+ int i = info_stack[stack_used];
info_stack[stack_used]++;
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
stack_used++;
@@ -422,6 +427,7 @@ namespace AMDiS {
return(elinfo_stack[stack_used]);
}
+
ElInfo* TraverseStack::traverseNeighbour(ElInfo* elinfo_old, int neighbour)
{
int dim = elinfo_old->getMesh()->getDim();
@@ -440,6 +446,7 @@ namespace AMDiS {
return NULL;
}
+
ElInfo* TraverseStack::traverseNeighbour3d(ElInfo* elinfo_old, int neighbour)
{
FUNCNAME("TraverseStackLLtraverseNeighbour3d()");
@@ -483,7 +490,6 @@ namespace AMDiS {
}
}
-
/* save information about current element and its position in the tree */
save_traverse_mel = traverse_mel;
save_stack_used = stack_used;
@@ -512,7 +518,8 @@ namespace AMDiS {
if (nb >= 0) { /* go to macro element neighbour */
i = traverse_mel->getOppVertex(nb);
traverse_mel = traverse_mel->getNeighbour(nb);
- if (traverse_mel == NULL) return(NULL);
+ if (traverse_mel == NULL)
+ return(NULL);
if (nb < 2 && save_stack_used > 1)
stack2_used = 2; /* go down one level in OLD hierarchy */
@@ -527,9 +534,9 @@ namespace AMDiS {
nb = i;
} else { /* goto other child */
stack2_used = stack_used + 1;
- if (save_stack_used > stack2_used) {
+ if (save_stack_used > stack2_used)
stack2_used++; /* go down one level in OLD hierarchy */
- }
+
elinfo2 = save_elinfo_stack[stack2_used];
el2 = elinfo2->getElement();
@@ -541,10 +548,8 @@ namespace AMDiS {
stack_used++;
info_stack[stack_used] = 0;
nb = 0;
-
}
-
elinfo = elinfo_stack[stack_used];
el = elinfo->getElement();
@@ -649,6 +654,7 @@ namespace AMDiS {
return(elinfo);
}
+
ElInfo* TraverseStack::traverseNeighbour2d(ElInfo* elinfo_old, int neighbour)
{
FUNCNAME("TraverseStack::traverseNeighbour2d()");
@@ -810,7 +816,6 @@ namespace AMDiS {
}
}
-
if (elinfo->getNeighbour(opp_vertex) != old_elinfo->getElement()) {
MSG(" looking for neighbour %d of element %d at %8X\n",
neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
@@ -845,6 +850,7 @@ namespace AMDiS {
return(elinfo);
}
+
void TraverseStack::update()
{
FUNCNAME("TraverseStack::update()");
@@ -855,6 +861,7 @@ namespace AMDiS {
dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
}
+
void TraverseStack::fillRefinementPath(ElInfo& elInfo, const ElInfo& upperElInfo)
{
int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo.getLevel();