COMPONENT library
PATTERN ".svn" EXCLUDE
)
+ FILE(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/deal_ii_source_includes
+ "${BOOST_FOLDER}/include\n"
+ )
ADD_SUBDIRECTORY(${BOOST_FOLDER}/libs/serialization/src)
FILES_MATCHING PATTERN "*.h"
PATTERN ".svn" EXCLUDE
)
+ FILE(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/deal_ii_source_includes
+ "${TBB_FOLDER}/include\n"
+ )
ADD_SUBDIRECTORY(${TBB_FOLDER}/src)
ENDIF()
##
## ---------------------------------------------------------------------
+#
+# This file sets up the project configuration consisting of
+#
+# Make.global_options (legacy)
+# deal.IIConfig.cmake
+# deal.IIVersionConfig.cmake
+#
+# and copies it (a) to the build directory and (b) prepares it for later
+# installation (the copy in CMAKE_CURRENT_BINARY_DIR).
+#
+
+#
+# Read in auxiliary include directories for the build directory
+# configuration:
+#
+FILE(STRINGS
+ ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/deal_ii_source_includes
+ deal_ii_source_includes
+ )
#
# Configure the template-arguments file
#
-CONFIGURE_FILE(
+CONFIGURE_FILE( # for binary dir:
${CMAKE_CURRENT_SOURCE_DIR}/template-arguments.in
- ${CMAKE_CURRENT_BINARY_DIR}/template-arguments
+ ${CMAKE_BINARY_DIR}/${DEAL_II_COMMON_RELDIR}/template-arguments
)
# #
########################################################################
-INSTALL(FILES
+FILE(COPY # for binary dir:
+ ${CMAKE_SOURCE_DIR}/cmake/macros/macro_deal_ii_initialize_cached_variables.cmake
+ ${CMAKE_SOURCE_DIR}/cmake/macros/macro_deal_ii_invoke_autopilot.cmake
+ ${CMAKE_SOURCE_DIR}/cmake/macros/macro_deal_ii_setup_target.cmake
+ DESTINATION ${CMAKE_BINARY_DIR}/${DEAL_II_CMAKE_MACROS_RELDIR}
+ )
+INSTALL(FILES # for installation:
${CMAKE_SOURCE_DIR}/cmake/macros/macro_deal_ii_initialize_cached_variables.cmake
${CMAKE_SOURCE_DIR}/cmake/macros/macro_deal_ii_invoke_autopilot.cmake
${CMAKE_SOURCE_DIR}/cmake/macros/macro_deal_ii_setup_target.cmake
)
ENDFOREACH()
-SET(CONFIG_INCLUDE_DIRS
+#
+# For binary dir:
+#
+SET(CONFIG_INCLUDE_DIRS_BINARY
+ \${DEAL_II_PATH}/include
+ \${DEAL_II_PATH}/include/deal.II
+ ${CMAKE_SOURCE_DIR}/include/
+ ${CMAKE_SOURCE_DIR}/include/deal.II
+ ${deal_ii_source_includes}
+ ${DEAL_II_USER_INCLUDE_DIRS}
+ )
+#
+# For installation:
+#
+SET(CONFIG_INCLUDE_DIRS_INSTALLATION
\${DEAL_II_PATH}/\${DEAL_II_INCLUDE_RELDIR}
\${DEAL_II_PATH}/\${DEAL_II_INCLUDE_RELDIR}/deal.II
\${DEAL_II_PATH}/\${DEAL_II_INCLUDE_RELDIR}/deal.II/bundled
${DEAL_II_USER_INCLUDE_DIRS}
)
-
CONFIGURE_FILE(
${CMAKE_CURRENT_SOURCE_DIR}/Config.cmake.in
- ${CMAKE_CURRENT_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake
@ONLY
)
-
#
# Append feature configuration:
#
-
-SET(_file ${CMAKE_CURRENT_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake)
+SET(_file
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake
+ )
FILE(APPEND ${_file}
"\n\n#\n# Feature configuration:\n#\n\n"
)
ENDIF()
ENDFOREACH()
-
CONFIGURE_FILE(
${CMAKE_CURRENT_SOURCE_DIR}/ConfigVersion.cmake.in
- ${CMAKE_CURRENT_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_NAME}ConfigVersion.cmake
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}ConfigVersion.cmake
@ONLY
)
INSTALL(FILES
- ${CMAKE_CURRENT_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake
- ${CMAKE_CURRENT_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_NAME}ConfigVersion.cmake
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}ConfigVersion.cmake
DESTINATION ${DEAL_II_PROJECT_CONFIG_RELDIR}
COMPONENT library
)
# Transform some cmake lists into a string that the old Makefile
# mechanism actually understands:
#
- TO_STRING_AND_ADD_PREFIX(MAKEFILE_INCLUDE_DIRS "${CMAKE_INCLUDE_FLAG_CXX}"
- $(D)/${DEAL_II_INCLUDE_RELDIR}
- $(D)/${DEAL_II_INCLUDE_RELDIR}/deal.II
- $(D)/${DEAL_II_INCLUDE_RELDIR}/deal.II/bundled
- ${DEAL_II_USER_INCLUDE_DIRS}
- ${DEAL_II_INCLUDE_DIRS}
- )
FOREACH(build ${DEAL_II_BUILD_TYPES})
TO_STRING_AND_ADD_PREFIX(MAKEFILE_DEFINITIONS_${build}
COND_SET_TO_YES(DEAL_II_WITH_P4EST MAKEFILE_P4EST)
COND_SET_TO_YES(DEAL_II_WITH_MPI MAKEFILE_MPI)
+ #
+ # For binary dir:
+ #
+ SET(MAKEFILE_D_BINARY ${CMAKE_BINARY_DIR})
+ TO_STRING_AND_ADD_PREFIX(MAKEFILE_INCLUDE_DIRS_BINARY "${CMAKE_INCLUDE_FLAG_CXX}"
+ $(D)/install
+ $(D)/install/deal.II
+ ${CMAKE_SOURCE_DIR}/include/
+ ${CMAKE_SOURCE_DIR}/include/deal.II
+ ${deal_ii_source_includes}
+ ${DEAL_II_USER_INCLUDE_DIRS}
+ )
+ #
+ # For installation:
+ #
+ SET(MAKEFILE_D_INSTALLATION ${CMAKE_INSTALL_PREFIX})
+ TO_STRING_AND_ADD_PREFIX(MAKEFILE_INCLUDE_DIRS_INSTALLATION "${CMAKE_INCLUDE_FLAG_CXX}"
+ $(D)/${DEAL_II_INCLUDE_RELDIR}
+ $(D)/${DEAL_II_INCLUDE_RELDIR}/deal.II
+ $(D)/${DEAL_II_INCLUDE_RELDIR}/deal.II/bundled
+ ${DEAL_II_USER_INCLUDE_DIRS}
+ )
CONFIGURE_FILE(
${CMAKE_CURRENT_SOURCE_DIR}/Make.global_options.in
- ${CMAKE_CURRENT_BINARY_DIR}/Make.global_options
+ ${CMAKE_BINARY_DIR}/${DEAL_II_COMMON_RELDIR}/Make.global_options
)
INSTALL(FILES
- ${CMAKE_CURRENT_BINARY_DIR}/template-arguments
- ${CMAKE_CURRENT_BINARY_DIR}/Make.global_options
+ ${CMAKE_BINARY_DIR}/${DEAL_II_COMMON_RELDIR}/template-arguments
+ ${CMAKE_BINARY_DIR}/${DEAL_II_COMMON_RELDIR}/Make.global_options
DESTINATION ${DEAL_II_COMMON_RELDIR}
COMPONENT compat_files
)
ENDIF()
+#
+# Finally, add a target to create the "binary" file in
+# ${DEAL_II_PROJECT_CONFIG_RELDIR} and add it to the library target:
+#
+ADD_CUSTOM_TARGET(setup_build_dir ALL
+ COMMAND ${CMAKE_COMMAND} -E touch
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/binary
+ COMMAND ${CMAKE_COMMAND} -E touch
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake
+ )
+ADD_DEPENDENCIES(library setup_build_dir)
+
+#
+# And a script to remove it upon installation from the install prefix:
+#
+INSTALL(CODE
+ "
+ FILE(REMOVE ${CMAKE_INSTALL_PREFIX}/${DEAL_II_PROJECT_CONFIG_RELDIR}/binary)
+ EXECUTE_PROCESS(COMMAND ${CMAKE_COMMAND} -E touch
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}Config.cmake
+ )
+ "
+ COMPONENT library
+ )
GET_FILENAME_COMPONENT(_path "${_path}" PATH)
ENDWHILE()
+#
+# Is this project reside in a build directory or at an installed location?
+#
+
+IF(EXISTS ${DEAL_II_PATH}/${DEAL_II_PROJECT_CONFIG_RELDIR}/binary)
+ SET(DEAL_II_BUILD_DIR TRUE)
+ENDIF()
+
#
# Print a message after inclusion of this file:
#
SET(DEAL_II_PROJECT_CONFIG_INCLUDED TRUE)
IF(NOT ${DEAL_II_PACKAGE_NAME}_FIND_QUIETLY)
- MESSAGE(STATUS
- "Using the ${DEAL_II_PACKAGE_NAME}-${DEAL_II_PACKAGE_VERSION} installation found at ${DEAL_II_PATH}"
- )
+ IF(DEAL_II_BUILD_DIR)
+ MESSAGE(STATUS
+ "Using the ${DEAL_II_PACKAGE_NAME}-${DEAL_II_PACKAGE_VERSION} build directory found at ${DEAL_II_PATH}"
+ )
+ ELSE()
+ MESSAGE(STATUS
+ "Using the ${DEAL_II_PACKAGE_NAME}-${DEAL_II_PACKAGE_VERSION} installation found at ${DEAL_II_PATH}"
+ )
+ ENDIF()
ENDIF()
#
# Information about include directories and libraries
#
-
-SET(DEAL_II_INCLUDE_DIRS "@CONFIG_INCLUDE_DIRS@")
+IF(DEAL_II_BUILD_DIR)
+ SET(DEAL_II_INCLUDE_DIRS "@CONFIG_INCLUDE_DIRS_BINARY@")
+ELSE()
+ SET(DEAL_II_INCLUDE_DIRS "@CONFIG_INCLUDE_DIRS_INSTALLATION@")
+ENDIF()
# Full list of libraries for the debug target:
SET(DEAL_II_LIBRARIES_DEBUG "@CONFIG_LIBRARIES_DEBUG@")
# Information about library targets
#
-SET(DEAL_II_TARGET_CONFIG "${DEAL_II_PATH}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}Targets.cmake")
+IF(DEAL_II_BUILD_DIR)
+ SET(DEAL_II_TARGET_CONFIG "${DEAL_II_PATH}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}BuildTargets.cmake")
+ELSE()
+ SET(DEAL_II_TARGET_CONFIG "${DEAL_II_PATH}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}Targets.cmake")
+ENDIF()
# The Debug target:
SET(DEAL_II_TARGET_DEBUG "@CONFIG_TARGET_DEBUG@")
#
-D = @CMAKE_INSTALL_PREFIX@
+ifeq ($(wildcard $(D)/@DEAL_II_PROJECT_CONFIG_RELDIR@/binary),)
+ D = @MAKEFILE_D_INSTALLATION@
+else
+ D = @MAKEFILE_D_BINARY@
+endif
+
CXX = @CMAKE_CXX_COMPILER@
CXX-ID = @CMAKE_CXX_COMPILER_ID@
CXX-VERSION = @CMAKE_CXX_COMPILER_VERSION@
# Include paths as command line flags
-INCLUDE = @MAKEFILE_INCLUDE_DIRS@
+ifeq ($(wildcard $(D)/@DEAL_II_PROJECT_CONFIG_RELDIR@/binary),)
+ INCLUDE = @MAKEFILE_INCLUDE_DIRS_INSTALLATION@
+else
+ INCLUDE = @MAKEFILE_INCLUDE_DIRS_BINARY@
+endif
# compiler flags for debug and optimized mode
# Define custom targets to easily switch the build type:
ADD_CUSTOM_TARGET(debug
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR}
+ COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
COMMENT "Switch CMAKE_BUILD_TYPE to Debug"
)
ADD_CUSTOM_TARGET(release
COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
+ COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
COMMENT "Switch CMAKE_BUILD_TYPE to Release"
)
MACRO(EXPAND_INSTANTIATIONS _target _inst_in_files)
- IF(CMAKE_CROSSCOMPILING)
- SET(_ei_exec native-expand_instantiations)
- ELSE()
- SET(_ei_exec expand_instantiations)
- ENDIF()
-
FOREACH (_inst_in_file ${_inst_in_files})
STRING(REGEX REPLACE "\\.in$" "" _inst_file "${_inst_in_file}" )
ADD_CUSTOM_COMMAND(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${_inst_file}
- DEPENDS ${_ei_exec}
- ${CMAKE_BINARY_DIR}/cmake/config/template-arguments
+ DEPENDS expand_instantiations
+ ${CMAKE_BINARY_DIR}/${DEAL_II_COMMON_RELDIR}/template-arguments
${CMAKE_CURRENT_SOURCE_DIR}/${_inst_in_file}
- COMMAND ${_ei_exec}
- ARGS ${CMAKE_BINARY_DIR}/cmake/config/template-arguments
+ COMMAND expand_instantiations
+ ARGS ${CMAKE_BINARY_DIR}/${DEAL_II_COMMON_RELDIR}/template-arguments
< ${CMAKE_CURRENT_SOURCE_DIR}/${_inst_in_file}
> ${CMAKE_CURRENT_BINARY_DIR}/${_inst_file}
)
TRILINOS_CONFIG_FOUND
)
-MARK_AS_ADVANCED(TRILINOS_CONFIG_DIR)
+MARK_AS_ADVANCED(TRILINOS_CONFIG EPETRA_CONFIG_H)
IF(TRILINOS_FOUND)
##
## ---------------------------------------------------------------------
-ADD_EXECUTABLE(expand_instantiations expand_instantiations.cc)
IF(NOT CMAKE_CROSSCOMPILING)
- #
- # Export expand_instantiations for use in cross compilation:
- #
+ ADD_EXECUTABLE(expand_instantiations expand_instantiations.cc)
+
+ IF(DEAL_II_COMPONENT_COMPAT_FILES)
+ SET_TARGET_PROPERTIES(expand_instantiations
+ PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/common/scripts
+ )
+ INSTALL(TARGETS expand_instantiations
+ EXPORT ${DEAL_II_PROJECT_CONFIG_NAME}Targets
+ DESTINATION ${DEAL_II_COMMON_RELDIR}/scripts
+ COMPONENT compat_files
+ )
+ ENDIF()
+
EXPORT(TARGETS expand_instantiations
- FILE ${CMAKE_CURRENT_BINARY_DIR}/importExecutables.cmake
- NAMESPACE native- # Prefix a namespace
+ FILE ${CMAKE_BINARY_DIR}/importExecutables.cmake
)
ENDIF()
# Build and install the old dependency_resolution and report_features
# executables:
#
+
INCLUDE_DIRECTORIES(
BEFORE # Ensure deal.II include directories come first
- ${CMAKE_BINARY_DIR}/include/
+ #
+ # Reverse order due to BEFORE:
+ #
${CMAKE_SOURCE_DIR}/include/
+ ${CMAKE_BINARY_DIR}/include/
)
ADD_EXECUTABLE(make_dependencies make_dependencies.cc)
ADD_CUSTOM_TARGET(run_report_features COMMAND report_features)
ENDIF()
- INSTALL(TARGETS make_dependencies report_features expand_instantiations
+ SET_TARGET_PROPERTIES(make_dependencies report_features
+ PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/common/scripts
+ )
+
+ INSTALL(TARGETS make_dependencies report_features
EXPORT ${DEAL_II_PROJECT_CONFIG_NAME}Targets
DESTINATION ${DEAL_II_COMMON_RELDIR}/scripts
COMPONENT compat_files
${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/deal_ii_objects_${_build_lowercase}
)
ENDFOREACH()
+FILE(WRITE
+ ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/deal_ii_source_includes
+ ""
+ )
+
+#
+# Cleanup deal.IITargets.cmake in the build directory:
+#
+FILE(REMOVE
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}BuildTargets.cmake
+ )
########################################################################
SET(DEAL_II_NATIVE "DEAL_II_NATIVE-NOTFOUND" CACHE FILEPATH
"A pointer to a native deal.Ii build directory"
)
- INCLUDE(${DEAL_II_NATIVE}/cmake/scripts/importExecutables.cmake)
+ IF(DEAL_II_NATIVE MATCHES "-NOTFOUND")
+ MESSAGE(FATAL_ERROR
+ "Please set the CMake variable DEAL_II_NATIVE to a valid path that points to a native deal.II build directory"
+ )
+ ENDIF()
+ INCLUDE(${DEAL_II_NATIVE}/importExecutables.cmake)
ENDIF()
{
QString pattern_description = index.data(Qt::StatusTipRole).toString(); // load pattern description
// stored in the StatusLine
- QRegExp rx_string("\\b(Anything|MultipleSelection|List)\\b"),
+ QRegExp rx_string("\\b(Anything|MultipleSelection|List|Map)\\b"),
rx_filename("\\b(FileName)\\b"),
rx_dirname("\\b(DirectoryName)\\b"),
rx_integer("\\b(Integer)\\b"),
* The ParameterDelegate class implements special delegates for the QTreeWidget class used in the parameterGUI.
* The QTreeWidget class provides some different standard delegates for editing parameters shown in the
* tree structure. The ParameterDelegate class provides special editors for the different types of parameters defined in
- * the ParameterHandler class. For all parameter types based on strings as "Anything", "MultipleSelection" and
+ * the ParameterHandler class. For all parameter types based on strings as "Anything", "MultipleSelection" "Map" and
* "List" a simple line editor will be shown up. In the case of integer and double type parameters the editor is a spin box and for
* "Selection" type parameters a combo box will be shown up. For parameters of type "FileName" and "DirectoryName"
* the delegate shows a @ref BrowseLineEdit editor. The column of the tree structure with the parameter values has to be set
you may receive an error message of the form:</p>
<pre>
-CMake Warning at source/CMakeLists.txt:65 (ADD_LIBRARY):
- Cannot generate a safe runtime search path for target deal_II.g because
- files in some directories may conflict with libraries in implicit
- directories:
- runtime library [libtbb.so.2] in /usr/lib may be hidden by files in:
- /my/private/lib
+ CMake Warning at source/CMakeLists.txt:65 (ADD_LIBRARY):
+ Cannot generate a safe runtime search path for target deal_II.g because
+ files in some directories may conflict with libraries in implicit
+ directories:
+
+ runtime library [libtbb.so.2] in /usr/lib may be hidden by files in:
+ /my/private/lib
+
+ Some of these libraries may not be found correctly.
- Some of these libraries may not be found correctly.
</pre>
<p>This is not a problem of CMake or deal.II, but rather a general
<meta name="svn_id" content="$Id$">
<meta name="keywords" content="deal.II">
</head>
-
+
<body>
<h1>How to use CMake to configure your projects with <acronym>deal.II</acronym></h1>
<ol>
<li><a href="#cmakesimple.multiple">Adding multiple executable targets</a></li>
<li><a href="#cmakesimple.libs">Adding libraries and common source files</a></li>
+ <li><a href="#cmakesimple.build_type">Switching build types</a></li>
<li><a href="#cmakesimple.run">Adding a "run" target</a></li>
</ol>
</li>
<h2>Simple CMakeLists.txt</h2>
<p>
- In this section, we start out with a
- minimal <code>CMakeLists.txt</code> based on
- the <code>DEAL_II_SETUP_TARGET</code> macro. This method gives
- full control of what's happening and is easily extensible to
- more complex projects, as exemplified in the subsections here an
- later in the section on <a href="cmakeadvanced">advanced
- topics</a>. Here is a full example
+ In this section, we start out with a minimal <code>CMakeLists.txt</code>
+ based on the <code>DEAL_II_SETUP_TARGET</code> macro. This method gives
+ full control of what's happening and is easily extensible to more complex
+ projects, as exemplified in the subsections here and later in the section
+ on <a href="cmakeadvanced">advanced topics</a>. Here is a full example
(<a href="CMakeLists.txt.sample3" target="_top">plain text</a>
version):
-
+
<p class="todo"> Fix plain text after finalizing!</p>
<pre class="cmake">
FIND_PACKAGE(deal.II 8.0 REQUIRED
HINTS ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
-)
+ )
DEAL_II_INITIALIZE_CACHED_VARIABLES()
-PROJECT(myproject CXX)
+PROJECT(myproject)
ADD_EXECUTABLE(mycode mycode.cc)
DEAL_II_SETUP_TARGET(mycode)
</p>
<p>
- Next, we find our deal.II installation, in this case requiring at
- least version 8.0, which is the first version using
- CMake. The <code>HINTS</code> are a list of directories where the
- install directory of deal.II is likely to be found. Here, we check
- whether we are in a subdirectory (first and second level) of the
- deal.II installation and otherwise use the
- variable <code>DEAL_II_DIR</code>.
- <q class="todo">Matthias, what is the difference between the first
- and the last entry?</q> This list can be changed according to your
- preferences. After finding the deal.II project, we fetch its chached
- variables. You can inspect these for instance
- with <code>ccmake</code>.
+ Next, we find our deal.II installation with the help of the
+ <code>FIND_PACKAGE</code> command. In this case requiring at least
+ version 8.0, which is the first version using CMake. The
+ <code>HINTS</code> are a list of directories where the install directory
+ of deal.II is likely to be found. First, the location possibly defined in
+ the CMake variable <code>DEAL_II_DIR</code> is considered. After that, we
+ check whether we are in a subdirectory (first and second level) of the
+ deal.II installation and otherwise use the environment variable
+ <code>DEAL_II_DIR</code>. If all of these hints fail the default system
+ locations <code>/usr/</code> and <code>/usr/local/</code> are considered.
+ The list after <code>HINTS</code> can be changed according to your
+ preferences.
</p>
<p>
- Every <code>CMAkeLists.txt</code> must contain a project definition,
- which we do next. <q class="todo">Matthias, I have been happy
- without the CXX...?</q>
+ After finding the deal.II project, we fetch a set of cached variables
+ with the <code>DEAL_II_INITIALIZE_CACHED_VARIABLES()</code> macro. You
+ can inspect these for instance with <code>ccmake</code>.
+</p>
+
+<p class="todo"> link to detailed explanation of the macro below</p>
+
+<p>
+ Every <code>CMakeLists.txt</code> must contain a project definition,
+ which we do next.
</p>
<p>
Finally, the last two lines define the executable that is to be
- produced and its source code.
+ produced and its source code. The <code>DEAL_II_SETUP_TARGET</code> macro
+ will set up necessary incllude directories, compile defintions and the
+ link interface.
</p>
-<p class="todo"> Explain the macros or link to explanation below</p>
+<p class="todo"> link to detailed explanation of the macro below</p>
<a name="cmakesimple.multiple"></a>
<h3>Adding multiple executable targets</h3>
<pre class="cmake">
ADD_LIBRARY(mylib libsrc1.cc libsrc2.cc libsrc3.cc)
+DEAL_II_SETUP_TARGET(mylib)
ADD_EXECUTABLE(mycode mycode.cc)
DEAL_II_SETUP_TARGET(mycode)
DEAL_II_SETUP_TARGET(mycode2)
</pre>
-<p>You should be aware though that <code>common.cc</code> will be
-compiled for each target, not only once.</p>
+<p>You should be aware though that in this case <code>common.cc</code> will
+be compiled for each target, not only once. If you want to avoid this and
+still don't want to use a shared library or static archive, another option
+is to create an <code>OBJECT</code> "library":</p>
+
+<pre class="cmake">
+ADD_LIBRARY(common OBJECT common.cc)
+DEAL_II_SETUP_TARGET(common)
-<p class="todo"> Matthias, is this correct?</p>
+ADD_EXECUTABLE(mycode mycode.cc $<TARGET_OBJECTS:common>)
+DEAL_II_SETUP_TARGET(mycode)
+ADD_EXECUTABLE(mycode2 mycode2.cc $<TARGET_OBJECTS:common>)
+DEAL_II_SETUP_TARGET(mycode2)
+</pre>
+This will compile <code>common.cc</code> once for the object target
+<code>common</code> and link the resulting object file into the two
+executables.
+
+<a name="cmakesimple.build_type"></a>
+<h3>Switching build types</h3>
+
+<p> The build type is controlled via the variable
+<code>CMAKE_BUILD_TYPE</code>. If it is set to <code>Debug</code>
+executables and libraries specified in your <code>CMakeLists.txt</code>
+file will be compiled in debug mode and linked against the debug version of
+the deal.II library. Contrary <code>Release</code> will build in optimized
+mode and link against the optimized release version of deal.II. You can set
+<code>CMAKE_BUILD_TYPE</code> with the help of <code>ccmake</code> or via
+<code>cmake</code> on the command line: </p>
- <a name="cmakesimple.run"></a>
- <h3>Adding a "run" target</h3>
-
- <p>
- If you wish to have a "run" target for make, like in the deal.II
- tutorial, specify one this way (obviously, a single "run" target
- can only run a single executable):
- </p>
+<pre class="cmake">
+$ cmake -DCMAKE_BUILD_TYPE="Debug" .
+
+$ cmake -DCMAKE_BUILD_TYPE="Release" .
+</pre>
+
+Alternatively, you can specify custom targets to switch the build type and
+compile automatically:
+
+<pre class="cmake">
+ADD_CUSTOM_TARGET(debug
+ COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR}
+ COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
+ COMMENT "Switch CMAKE_BUILD_TYPE to Debug"
+ )
+
+ADD_CUSTOM_TARGET(release
+ COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
+ COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
+ COMMENT "Switch CMAKE_BUILD_TYPE to Release"
+ )
+</pre>
+With that, switching the build type and compiling the project can be done
+very conveniently via:
+<pre class="cmake">
+$ make debug
+
+$ make release
+</pre>
+
+<a name="cmakesimple.run"></a>
+<h3>Adding a "run" target</h3>
+
+<p> If you wish to have a "run" target for make, like in the deal.II
+tutorial, specify one this way (obviously, a single "run" target can only
+run a single executable): </p>
<pre class="cmake">
ADD_CUSTOM_TARGET(run COMMAND mycode
- COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration")
+ COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration"
+ )
</pre>
<a name="cmakeadvanced"></a>
<h2> Advanced <code>CMakeLists.txt</code></h2>
-<pre class="todo">
-
- A Subsection for each of the following:
-
- - control statements (if and foreach)
-
- - file glob to pick up sources
-
- - DEAL_II_SETUP_TARGET and DEAL_II_INITIALIZE_CACHED_VARIABLES
- revisited
-
- - add include dirs and compile definitions to a directory or target
-
- - provide a ./run folder and target
-
- - custom targets to switch between release and debug
-
- - installation
+<p class="todo">This section covers some advanced topics for a user
+<code>CMakeLIists.txt</code> file.</p>
+
+<a name="cmakeadvanced.control"></a>
+<h3>Control statements</h3>
+
+<p>Control statements in CMake take the following form:
+<pre class="cmake">
+IF(<expression>)
+ ...
+ENDIF()
+</pre>
+or in long form:
+<pre class="cmake">
+IF(<expression1>)
+ ...
+ELSEIF(<expression2>)
+ ...
+ELSE()
+ ...
+ENDIF()
+</pre>
+Please note the (somehow uncommon) empty, opening and closing brackets
+behind <code>ELSE()</code> and <code>ENDIF()</code>.
+<code><expression></code> can take a multitude of different forms,
+have a look at the
+<a href="http://cmake.org/cmake/help/v2.8.8/cmake.html">CMake
+documentation</a> for a complete list. Important examples are:
+<p class="todo">Link!</p>
+<pre class="cmake">
+IF(${variable})
+ - the body will be evaluated if the variable "variable" is defined and
+ synonymous to true, e.g. 1, TRUE, ON, YES (modulo case insensitivity)
+
+IF(variable MATCHES <regular expression>)
+ - the body will be evaluated if the variable "variable" is defined and
+ matches the specified regular expression
+
+IF("${variable}" STREQUAL "foobar")
+ - the body will be evaluated if both strings are equal. Note that
+ "${variable}" will be replaced by the content of the (string)
+ variable "variable"
+</pre>
+<p>An expression can be negated by prefixing <code>NOT</code>:</p>
+<pre class="cmake">
+IF(NOT <expression>)
+ ...
+ENDIF()
</pre>
-<a name="extfinding"></a>
-<h3> Finding the deal.II library </h3>
+<p>Loops are implemented with the help of <code>WHILE</code> and
+<code>FOR</code> statements. The former takes the same
+<code><expression></code> as the <code>IF</code> statement:</p>
+<pre class="cmake">
+WHILE(<expression>)
+ ...
+ENDWHILE()
+</pre>
+Given a variable <code>list</code> containing a list, the individual
+elements <code>element</code> can be accessed with a <code>FOREACH</code>
+statement:
+<pre class="cmake">
+FOREACH(element ${list})
+ ...
+ENDFOREACH
+</pre>
+Note: It is also possible to specify the list directly:
+<pre class="cmake">
+FOREACH(element foo bar baz)
+ # The variable element will iterate through foo, bar and baz.
+ENDFOREACH
+</pre>
-<p>
- Finding the <acronym>deal.II</acronym> library should be no more than
- <pre>
-
- FIND_PACKAGE(deal.II REQUIRED)
- </pre>
- in your CMakeLists.txt file. You may have to hint for the location
- of the <code>deal.IIConfig.cmake</code> file. Either by directly
- specifying <code>deal.II_DIR</code> to point to the path were the
- <code>deal.IIConfig.cmake</code> file is located:
-
-<p class="todo"> is deal.II_DIR still accurate? We use different
-above. BTW, that is an <b>UGLY</b> variable name.</p>
-
-<pre>
-cmake -Ddeal.II_DIR=/path/to/the/config/file <...>
+<a name="cmakeadvanced.globs"></a>
+<h3>File globs</h3>
+
+<p>A very common task is to pick up a list of source files from a
+directory. You can either manage a list of source files in
+<code>CMakeLists.txt</code> by hand, e.g. by manually updating all source
+files for a given target, or you can use a glob to automate this process.
+The following example will pick up every source file under
+<code>SOURCE_DIR/sources/</code> and add it to an executable:</p>
+
+<pre class="cmake">
+FILE(GLOB sources ${CMAKE_SOURCE_DIR}/source/*.cc)
+ADD_EXECUTABLE(mycode ${sources})
</pre>
- or by specifying a search path via <code>CMAKE_PREFIX_PATH</code>,
- e.g.
- <pre>
-
- cmake -DCMAKE_PREFIX_PATH=~/workspace/local
- </pre>
- In this case <code>deal.IIConfig.cmake</code> will be searched
- for in
- <pre>
-
- ~/workspace/local/
- ~/workspace/local/lib/cmake/deal.II/
- </pre>
-</p>
+
+<p>Please be aware of one caveat of this approach: Due to the fact that
+CMake is a <i>build system generator</i> the resulting build configuration
+(for make) has no way to detect whether a new source file was added (or
+removed) and that it has to call back to cmake. So, after adding a new
+source file you have to touch a <code>CMakeLists.txt</code> file or to run
+<code>cmake .</code> again by hand.
+
+
+<a name="cmakeadvanced.setup_target"></a>
+<h3><code>DEAL_II_SETUP_TARGET</code> revisited</h3>
+
+<p class="todo"> Well, write this! </p>
+
+
+<a name="cmakeadvanced.cached_variables"></a>
+<h3><code>DEAL_II_INITIALIZE_CACHED_VARIABLES</code> revisited</h3>
+
+<p class="todo"> Well, write this! </p>
+
+
+<a name="cmakeadvanced.properties"></a>
+<h3>Customizing include directories and compile definitions</h3>
+
+<p class="todo"> Well, write this! </p>
+
+
+<a name="cmakeadvanced.external_libraries"></a>
+<h3>External libraries</h3>
+
+<p class="todo"> Well, write this! </p>
+
+
+<a name="cmakeadvanced.run"></a>
+<h3>The "run" target revisited</h3>
+
+<p class="todo"> Well, write this! </p>
+
+
+<a name="cmakeadvanced.install"></a>
+<h3>Install a project</h3>
+
<a name="cmakesauto"></a>
<h2>Autopilot style CMakeLists.txt</h2>
# directories that contain example code fragments that are included (see
# the \include command).
-EXAMPLE_PATH = @CMAKE_BINARY_DIR@/doc/doxygen/tutorial
+EXAMPLE_PATH = @CMAKE_BINARY_DIR@/doc/doxygen/tutorial \
+ @CMAKE_SOURCE_DIR@/examples/doxygen
# If the value of the EXAMPLE_PATH tag contains directories, you can use the
# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
<ol>
+ <li>
+ New: Like the usual DoFHandler class, the hp::DoFHandler class now also
+ has a cache that makes operations such as <code>cell-@>get_dof_indices(...)</code>
+ faster. This should accelerate many parts of the library that deal with
+ hp finite elements.
+ <br>
+ (Wolfgang Bangerth, 2013/09/10)
+ </li>
+
+ <li>
+ New: parallel::distributed::Triangulation now supports periodic boundaries.
+ DoFTools::make_periodicity_constraints and similar functions are now working
+ on parallel::distributed::Triangulation objects.
+ <br>
+ (Tobin Isaac, Craig Michoski, Daniel Arndt, 2013/09/06)
+ </li>
+
<li>
New: It is now possible to compile and link deal.II against LLVM's libcxx. For
this, a few issues with C++ standard violations are resolved.
<h3>Specific improvements</h3>
<ol>
+ <li>
+ Enhancement: It is now possible to use the build directory directly without
+ the need to install first. For this, a second copy of all necessary project
+ configuration files (deal.IIConfig.cmake, ..., Make.global_options) are
+ generated and deployed in the build directory. (This is fully compatible with
+ the old possibility to install into the build dir.)
+ <br>
+ (Matthias Maier, 2013/09/15)
+ </li>
+
+ <li>
+ Fixed: DoFTools::extract_locally_*_dofs now instantiated for hp::DofHandler.
+ <br>
+ (Jean-Paul Pelteret, 2013/09/11)
+ </li>
+
+ <li>
+ Changed: distributed::parallel:BlockVector::operator= now allows importing
+ of ghost values like all other vector types. Also added some new constructors
+ for BlockVector and Vector using IndexSets to mirror the other linear algebra
+ classes.
+ <br>
+ (Timo Heister, 2013/09/04)
+ </li>
+
+ <li>
+ Fixed: VectorTools::compute_no_normal_flux_constraints had a bug that
+ only manifested on complex meshes. This is now fixed.
+ <br>
+ (Chih-Che Chueh, Wolfgang Bangerth, 2013/09/04)
+ </li>
+
<li>
New: All vector classes now have functions <code>extract_subvector_to()</code>
that allow extracting not just a single value but a whole set.
conductivity and a time scale of relevance to affect the evolution of the
earth's interior and surface structure.
+@note If you are interested in using the program as the basis for your own
+experiments, you will also want to take a look at its continuation in
+step-32. Furthermore, step-32 later was developed into the much larger open
+source code Aspect (see http://aspect.dealii.org/ ) that can solve realistic
+problems and that you may want to investigate before trying to morph step-31
+into something that can solve whatever you want to solve.
<h3>%Boundary and initial conditions</h3>
if we state the CFL condition as the requirement that the time step be
small enough so that the distance transport advects in each time step
is no longer than one <i>grid point</i> away (which for $Q_1$ elements
-is $h_K$, but for $Q_2$ elements is $h_K/2$). It turns out that $\beta$
+is $h_K$, but for $Q_2$ elements is $h_K/2$). It turns out that $\beta$
needs to be slightly larger for obtaining stable results also late in
the simulation at times larger than 60, so we actually choose it as
$\beta = 0.034$ in the code.
is, of course, to make it faster and/or increase the resolution of the
program, in particular in 3d. This is the topic of the step-32
tutorial program which will implement strategies to solve this problem in
-%parallel on a cluster.
+%parallel on a cluster. It is also the basis of the much larger open
+source code Aspect (see http://aspect.dealii.org/ ) that can solve realistic
+problems and that constitutes the further development of step-32.
Another direction would be to make the fluid flow more realistic. The program
was initially written to simulate various cases simulating the convection of
<b>
M. Kronbichler, T. Heister, W. Bangerth:
<i>High Accuracy Mantle Convection Simulation through Modern Numerical
- Methods</i>, Geophysical Journal International, 2012, 191, 12-29.
+ Methods</i>, Geophysical Journal International, 2012, 191, 12-29.
<a href="http://dx.doi.org/10.1111/j.1365-246X.2012.05609.x">[DOI]</a>
</b>
+
+The continuation of development of this program has led to the much larger
+open source code Aspect (see http://aspect.dealii.org/ ) which is much more
+flexible in solving many kinds of related problems.
</i>
here, namely the KellyErrorEstimator applied to the temperature, is
able to.
</ul>
+
+
+<h3> Possible extensions </h3>
+
+There are many ways to extend the current program. However, rather than
+discussing them here, let us point to the much larger open
+source code Aspect (see http://aspect.dealii.org/ ) that constitutes the
+further development of step-32 and that already includes many such possible
+extensions.
namespace EquationData
{
const double eta = 1e21; /* Pa s */
- const double kappa = 1e-6; /* m / s */
+ const double kappa = 1e-6; /* m^2 / s */
const double reference_density = 3300; /* kg / m^3 */
const double reference_temperature = 293; /* K */
const double expansion_coefficient = 2e-5; /* 1/K */
}
}
}
-
+
joint_solution.compress(VectorOperation::insert);
IndexSet locally_relevant_joint_dofs(joint_dof_handler.n_dofs());
+<br>
+
+<i>
+This program was contributed by Martin Kronbichler and Scott Miller.
+</i>
+
<a name="Intro"></a>
<h1>Introduction</h1>
One common argument against the use of discontinuous Galerkin elements
is the large number of globally coupled degrees of freedom that one
-must solve in an implicit system. For the FE_DGP_Monomial basis, each
-scalar solution component is represented by polynomials of degree p
-which yields $(1/dim!)*\prod_{i=1}^{dim}(k+i)$ degrees of freedom per
+must solve in an implicit system. This is because, unlike continuous finite
+elements, in typical discontinuous elements there is one degree of freedom at
+each vertex <i>for each of the adjacent elements</i>, rather than just one,
+and similarly for edges and faces. As another example,
+for the FE_DGP_Monomial basis, each
+scalar solution component is represented by polynomials of degree $p$
+which yields $(1/dim!)*\prod_{i=1}^{dim}(p+i)$ degrees of freedom per
element. Typically, all degrees of freedom in an element are coupled
to all of the degrees of freedom in the adjacent elements. The resulting
discrete equations yield very large linear systems very quickly, especially
<h4> Reducing the size of the linear system </h4>
To alleviate the computational cost of solving such large linear systems,
-the hybridizable discontinuous Galerkin (HDG) methodology has recently been
-developed by Cockburn and co-workers
-(<b>
- N.C. Nguyen and J. Peraire:
+the hybridizable discontinuous Galerkin (HDG) methodology was introduced
+by Cockburn and co-workers
+(see the references in the recent HDG overview article by
+ N.C. Nguyen and J. Peraire:
<i>Hybridizable discontinuous Galerkin methods for partial differential
equations in continuum mechanics</i>, Journal of Computational Physics,
2012, 231:18, 5955-5988.
- <a href="http://dx.doi.org/10.1016/j.jcp.2012.02.033">[DOI]</a>
-</b>).
+ <a href="http://dx.doi.org/10.1016/j.jcp.2012.02.033">[DOI]</a>).
The HDG method achieves
this goal by formulating the mathematical problem using Dirichlet-to-Neumann
mappings. The partial differential equations are first written as a first
order system, and each field is then discretized via a DG method. At this
-point the single-valued ``trace'' values on the skeleton of the
-mesh, i.e. element edges, are taken to be independent unknown quantities.
-The Dirichlet-to-Neumann map concept then permits the solution procedure:
+point the single-valued "trace" values on the skeleton of the
+mesh, i.e. element faces, are taken to be independent unknown quantities.
+The Dirichlet-to-Neumann map concept then permits the following solution procedure:
<ol>
- <li> Solve a Dirichlet problem on each element individually, where the local/interior
-element solutions no longer couple to neighboring elements. This is known as the
+ <li> Use local element interior data to enforce a Neumann condition on the
+skeleton of the triangulation. The global problem is then to solve for the
+trace values, which are the only globally coupled unknowns.
+ <li> Use the known skeleton values as Dirichlet data for solving local
+element-level solutions. This is known as the
'local solver', and is an <i>embarrassingly parallel</i> element-by-element
solution process.
- <li> Use the local element equations to construct the global problem for the
-trace values. These are the only globally coupled unknowns.
- <li> Solve for the element-local solutions using the computed trace values.
</ol>
The above procedure also has a linear algebra interpretation and is referred to
(D - C A^{-1} B) \Lambda &=& G - C A^{-1} F \\
A U &=& F - B \Lambda
@f}
-The steps in the Dirichlet-to-Neumann map concept hence correspond to
+The point is that the presence of $A^{-1}$ is not a problem because $A$ is a
+block diagonal matrix where each block corresponds to one cell and is
+therefore easy enough to invert.
+The coupling to other cells is introduced by the matrices
+<i>B</i> and <i>C</i> over the skeleton variable. The block-diagonality of
+<i>A</i> and the structure in <i>B</i> and <i>C</i> allow us to invert the
+matrix <i>A</i> element by element (the local solution of the Dirichlet
+problem) and subtract $CA^{-1}B$ from $D$. The steps in the Dirichlet-to-Neumann map concept hence correspond to
<ol>
<li> constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$,
<li> solving the Schur complement system for $\Lambda$, and
<li> solving the equation for <i>U</i> using the second equation which uses $\Lambda$.
</ol>
-The important ingredient from the linear algebra point of view is that the
-matrix <i>A</i> is block-diagonal with block size equal to the number of
-degrees of freedom of the interior DG variables which are always only related
-to a single cell. The coupling to other cells is introduced by the matrices
-<i>B</i> and <i>C</i> over the skeleton variable. The block-diagonality of
-<i>A</i> and the structure in <i>B</i> and <i>C</i> allow us to invert the
-matrix <i>A</i> element by element (the local solution of the Dirichlet
-problem) and subtract $CA^{-1}B$ from $D$.
<h4> Solution quality and rates of convergence</h4>
Another criticism of traditional DG methods is that the approximate fluxes
converge suboptimally. The local HDG solutions can be shown to converge
-as $\mathcal{O}(h^{p+1})$, i.e., at optimal order. Additionally, a super-convergence property can
-be used to post-process a new approximate solution that converges at the
-rate $\mathcal{O}(h^{p+2})$.
+as $\mathcal{O}(h^{p+1})$, i.e., at optimal order. Additionally, a
+super-convergence property can be used to post-process a new approximate
+solution that converges at the rate $\mathcal{O}(h^{p+2})$.
We multiply these equations by the weight functions $\mathbf{v}, w$
and integrate by parts over every element $K$ to obtain:
@f{eqnarray*}
- (\kappa^{-1} \mathbf{q}, \mathbf{v})_K - (u, \nabla\cdot\mathbf{v})_K
- + <\hat{u}, \mathbf{v}\cdot\mathbf{n}>_{\partial K} &=& 0, \\
- - (\mathbf{c} u + \mathbf{q}, \nabla w)_K
- + <(\hat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}, w>_{\partial K}
- &=& (f,w)_K.
+ (\mathbf{v}, \kappa^{-1} \mathbf{q})_K - (\nabla\cdot\mathbf{v}, u)_K
+ + \left<\mathbf{v}\cdot\mathbf{n}, \hat{u}\right>_{\partial K} &=& 0, \\
+ - (\nabla w, \mathbf{c} u + \mathbf{q})_K
+ + \left<(w, \hat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}\right>_{\partial K}
+ &=& (w,f)_K.
@f}
-The terms decorated with a $\hat{}$ denote the numerical traces (also commonly referred
+The terms decorated with a hat denote the numerical traces (also commonly referred
to as numerical fluxes). They are approximations
to the interior values on the boundary of the element. To ensure conservation,
-these terms must be single-valued on any given element edge $\partial K$.
+these terms must be single-valued on any given element edge $\partial K$ even
+though, with discontinuous shape functions, there may of course be multiple
+values coming from the cells adjacent to an interface.
We eliminate the numerical trace $\hat{\mathbf{q}}$ by using traces of the form:
@f{eqnarray*}
\hat{\mathbf{c} u}+\hat{\mathbf{q}} = \mathbf{c}\hat{u} + \mathbf{q}
+ \tau(u - \hat{u})\mathbf{n} \quad \text{ on } \partial K.
@f}
-The local stabilization parameter $\tau$ has effects on stability and accuracy of
-HDG solutions; see the literature for a further discussion.
+The local stabilization parameter $\tau$ has effects on stability and accuracy
+of HDG solutions; see the literature for a further discussion. A stabilization
+parameter of unity is reported to be the choice which gives best results. A
+stabilization parameter $\tau$ that tends to infinity prohibits jumps in the
+solution over the element boundaries, making the HDG solution approach the
+approximation of continuous finite elements. In the program below, we choose
+the stabilization parameter as
+@f{eqnarray*}
+ \tau = \frac{\kappa}{\ell} + |\mathbf{c} \cdot \mathbf{n}|
+@f}
+where we set the diffusion length scale to $\ell = \frac{1}{5}$.
The trace/skeleton variables in HDG methods are single-valued on element faces. As such,
they must strongly represent the Dirichlet data on $\partial\Omega_D$. We introduce
@f}
Eliminating $\hat{u}$ from the weak form in favor of $\lambda$, summing the elemental
-contributions across all elements in the triangulation, and enforcing the normal
-component of the numerical flux, we arrive at the final form of the problem:
-Find $(\mathbf{q}, u, \lambda)$ such that
+contributions across all elements in the triangulation, enforcing the normal
+component of the numerical flux, and integrating by parts
+on the equation weighted by $w$, we arrive at the final form of the problem:
+Find $(\mathbf{q}_h, u_h, \lambda_h) \in
+\mathcal{V}_h^p \times \mathcal{W}_h^p \times \mathcal{M}_h^p$ such that
@f{eqnarray*}
- (\kappa^{-1} \mathbf{q}, \mathbf{v})_{\mathcal{T}}
- - (u, \nabla\cdot\mathbf{v})_{\mathcal{T}}
- + <\lambda, \mathbf{v}\cdot\mathbf{n}>_{\partial\mathcal{T}}
+ (\mathbf{v}, \kappa^{-1} \mathbf{q}_h)_{\mathcal{T}}
+ - ( \nabla\cdot\mathbf{v}, u_h)_{\mathcal{T}}
+ + \left<\mathbf{v}\cdot\mathbf{n}, \lambda_h\right>_{\partial\mathcal{T}}
&=&
- - <g_D, \mathbf{v}\cdot\mathbf{n}>_{\partial\Omega_D},
+ - \left<\mathbf{v}\cdot\mathbf{n}, g_D\right>_{\partial\Omega_D},
+ \quad \forall \mathbf{v} \in \mathcal{V}_h^p,
\\
- - (\mathbf{c} u + \mathbf{q}, \nabla w)_{\mathcal{T}}
- + <(\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}, w>_{\partial \mathcal{T}}
+ (w, \mathbf{c}\nabla u_h + \nabla \cdot \mathbf{q}_h)_{\mathcal{T}}
+ + \left<w, \tau (u_h - \lambda_h)\right>_{\partial \mathcal{T}}
&=&
- (f, w)_{\mathcal{T}},
+ (w, f)_{\mathcal{T}},
+ \quad \forall w \in \mathcal{W}_h^p,
\\
- < \left[ \negthinspace \left[ (\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}
- \right] \negthinspace \right], \mu>_{\partial \mathcal{T}}
+ \left< \mu, \mathbf{c} \lambda_h\cdot \mathbf{n}
+ + \mathbf{q}_h\cdot \mathbf{n}
+ + \tau (u_h - \lambda_h)\right>_{\partial \mathcal{T}}
&=&
- <g_N, \mu>_{\partial\Omega_N}
+ \left<\mu, g_N\right>_{\partial\Omega_N},
+ \quad \forall \mu \in \mathcal{M}_h^p.
@f}
+The unknowns $(\mathbf{q}_h, u_h)$ are referred to as local variables; they are
+represented as standard DG variables. The unknown $\lambda_h$ is the skeleton
+variable which has support on the codimension-1 surfaces (faces) of the mesh.
+In the weak form given above, we can note the following coupling patterns:
+<ol>
+ <li> The matrix $A$ consists of local-local coupling terms. These arise when the
+ local weighting functions $(\mathbf{v}, w)$ multiply the local solution terms
+ $(\mathbf{q}_h, u_h)$.
+ <li> The matrix $B$ represents the local-face coupling. These are the terms
+ with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable
+ $\lambda_h$.
+ <li> The matrix $C$ represents the face-local coupling, which involves the
+ weighting function $\mu$ multiplying the local solutions $(\mathbf{q}_h, u_h)$.
+ <li> The matrix $D$ is the face-face coupling;
+ terms involve both $\mu$ and $\lambda_h$.
+</ol>
+
+<h4> Post-processing and super-convergence </h4>
<h3> Problem specific data </h3>
- WorkStream to parallelize local solvers. Workstream is already used in step-32, step-44.
- Reconstructing the trace
- Post-processing the solution for superconvergence
-- DataOutFaces: direct output of the global solution
\ No newline at end of file
+- DataOutFaces: direct output of the global solution
+
One can see the error reduction upon grid refinement, and for the cases where
-global refinement was performed, also the convergence rates can be seen. The
-quadratic convergence rates of Q1 elements in the $L_2$ norm for both the
-scalar variable and the gradient variable can clearly be seen, as is the cubic
-rate for the postprocessed scalar variable in the $L_2$ norm. Likewise, the
-scalar variable and gradient for Q3 elements converge at fourth order and the
-postprocessed scalar variable at fifth order.
+global refinement was performed, also the convergence rates. The quadratic
+convergence rates of Q1 elements in the $L_2$ norm for both the scalar
+variable and the gradient variable is apparent, as is the cubic rate for the
+postprocessed scalar variable in the $L_2$ norm. Note that is is a distinctive
+feature of an HDG solution. In typical continuous finite element, the gradient
+of the solution of order <i>p</i> converges at rate <i>p</i> only, as opposed
+to <i>p</i>+1 for the actual solution. Even though superconvergence results
+for finite elements are also available (e.g. superconvergent patch recovery
+first introduced by Zienkiewicz and Zhu), these are typically limited to
+structured meshes and other special cases. Likewise, the scalar variable and
+gradient for Q3 elements converge at fourth order and the postprocessed scalar
+variable at fifth order.
The same convergence rates are observed in 3d.
@code
improvement makes most sense.
<ol>
- <li> Better linear solvers: We use a GMRES iterative solver without
+ <li> Better linear solvers: We use a BiCGStab iterative solver without
preconditioner, where the number of iteration increases with increasing
problem size (the number of iterations for Q1 elements and global
- refinements start at 8 for the small sizes but increase up to 2592 for the
+ refinements start at 35 for the small sizes but increase up to 701 for the
largest size). To do better, one could for example use an algebraic
multigrid preconditioner from Trilinos. For diffusion-dominated problems as
the problem at hand with finer meshes, such a solver can be designed that
<li> Speed up assembly by pre-assembling parts that do not change from one
cell to another (those that do neither contain variable coefficients nor
- mapping-dependent terms.
+ mapping-dependent terms).
</ol>
// @sect3{Include files}
//
// Most of the deal.II include files have already been covered in previous
-// examples are are not commented on.
+// examples and are not commented on.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/tensor_function.h>
// the simulation.
#include <deal.II/numerics/data_out_faces.h>
+namespace Step51
+{
+
using namespace dealii;
// @sect3{Equation data}
* this->width);
}
-// @sect3{The Step51 HDG solver class}
+// @sect3{The HDG solver class}
// The HDG solution procedure follows closely that of step-7. The major
-// difference is the use of 3 different sets of <code>DoFHandler</code> and FE
+// difference is the use of three different sets of <code>DoFHandler</code> and FE
// objects, along with the <code>ChunkSparseMatrix</code> and the
// corresponding solutions vectors. We also use WorkStream to enable a
// multithreaded local solution process which exploits the embarrassingly
// parallel nature of the local solver. For WorkStream, we define the local
// operations on a cell and a copy function into the global matrix and
-// vector. We do this once for the assembly (which is run twice, once when we
+// vector. We do this both for the assembly (which is run twice, once when we
// generate the system matrix and once when we compute the element-interior
-// solutions from the skeleton values) and once for the postprocessing where
+// solutions from the skeleton values) and for the postprocessing where
// we extract a solution that converges at higher order.
template <int dim>
-class Step51
+class HDG
{
public:
enum RefinementMode
global_refinement, adaptive_refinement
};
- Step51 (const unsigned int degree,
+ HDG (const unsigned int degree,
const RefinementMode refinement_mode);
void run ();
private:
+// Data for the assembly and solution of the primal variables.
struct PerTaskData;
struct ScratchData;
+// Post-processing the solution to obtain $u^*$ is an element-by-element
+// procedure; as such, we do not need to assemble any global data and do
+// not declare any 'task data' for WorkStream to use.
struct PostProcessScratchData;
void setup_system ();
void assemble_system (const bool reconstruct_trace = false);
+ void solve ();
+ void postprocess ();
+
+ void refine_grid (const unsigned int cylce);
+ void output_results (const unsigned int cycle);
+
+// The following three functions are used by WorkStream to do the actual work of
+// the program.
void assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData &scratch,
PerTaskData &task_data);
+
void copy_local_to_global(const PerTaskData &data);
- void solve ();
- void postprocess ();
+
void postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
PostProcessScratchData &scratch,
unsigned int &empty_data);
- void refine_grid (const unsigned int cylce);
- void output_results (const unsigned int cycle);
+
Triangulation<dim> triangulation;
Vector<double> solution_local;
// The new finite element type and corresponding <code>DoFHandler</code> are
- // used for the global solution that couples the element-level local
+ // used for the global skeleton solution that couples the element-level local
// solutions.
FE_FaceQ<dim> fe;
DoFHandler<dim> dof_handler;
ConvergenceTable convergence_table;
};
+// @sect3{The HDG class implementation}
-
+// @sect4{Constructor} The constructor is similar to those in other examples,
+// with the exception of handling multiple <code>DoFHandler</code> and
+// <code>FiniteElement</code> objects. Note that we create a system of finite
+// elements for the local DG part, including the gradient/flux part and the
+// scalar part.
template <int dim>
-Step51<dim>::Step51 (const unsigned int degree,
+HDG<dim>::HDG (const unsigned int degree,
const RefinementMode refinement_mode) :
fe_local (FE_DGQ<dim>(degree), dim,
FE_DGQ<dim>(degree), 1),
-// First come the definition of the local data structures for the parallel
-// assembly. The first structure @p PerTaskData contains the local vector and
-// matrix that are written into the global matrix, whereas the ScratchData
-// contains all data that we need for the local assembly.
+// @sect4{HDG::setup_system}
+// The system for an HDG solution is setup in an analogous manner to most
+// of the other tutorial programs. We are careful to distribute dofs with
+// all of our <code>DoFHandler</code> objects. The @p solution and @p system_matrix
+// objects go with the global skeleton solution.
+template <int dim>
+void
+HDG<dim>::setup_system ()
+{
+ dof_handler_local.distribute_dofs(fe_local);
+ dof_handler.distribute_dofs(fe);
+ dof_handler_u_post.distribute_dofs(fe_u_post);
+
+ std::cout << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+
+ solution_local.reinit (dof_handler_local.n_dofs());
+ solution_u_post.reinit (dof_handler_u_post.n_dofs());
+
+ constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ typename FunctionMap<dim>::type boundary_functions;
+ Solution<dim> solution_function;
+ boundary_functions[0] = &solution_function;
+ VectorTools::project_boundary_values (dof_handler,
+ boundary_functions,
+ QGauss<dim-1>(fe.degree+1),
+ constraints);
+ constraints.close ();
+
+ {
+ CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, csp,
+ constraints, false);
+ sparsity_pattern.copy_from(csp, fe.dofs_per_face);
+ }
+ system_matrix.reinit (sparsity_pattern);
+}
+
+
+
+// @sect4{HDG::PerTaskData} Next come the definition of the local data
+// structures for the parallel assembly. The first structure @p PerTaskData
+// contains the local vector and matrix that are written into the global
+// matrix, whereas the ScratchData contains all data that we need for the
+// local assembly. There is one variable worth noting here, namely the boolean
+// variable @p trace_reconstruct. As mentioned introdution, we solve the HDG
+// system in two steps. First, we create a linear system for the skeleton
+// system where we condense the local part into it by $D-CA^{-1}B$. Then, we
+// solve for the local part using the skeleton solution. For these two steps,
+// we need the same matrices on the elements twice, which we want to compute
+// by two assembly steps. Since most of the code is similar, we do this with
+// the same function but only switch between the two based on a flag that we
+// set when starting the assembly. Since we need to pass this information on
+// to the local worker routines, we store it once in the task data.
template <int dim>
-struct Step51<dim>::PerTaskData
+struct HDG<dim>::PerTaskData
{
FullMatrix<double> cell_matrix;
Vector<double> cell_vector;
dof_indices(n_dofs),
trace_reconstruct(trace_reconstruct)
{}
-
- void reset()
- {
- cell_matrix = 0.0;
- cell_vector = 0.0;
- }
};
+// @sect4{HDG::ScratchData} @p ScratchData contains persistent data for each
+// thread within <code>WorkStream</code>. The <code>FEValues</code>, matrix,
+// and vector objects should be familiar by now. There are two objects that
+// need to be discussed: @p std::vector<std::vector<unsigned int> >
+// fe_local_support_on_face and @p std::vector<std::vector<unsigned int> >
+// fe_support_on_face. These are used to indicate whether or not the finite
+// elements chosen have support (non-zero values) on a given face of the
+// reference cell for the local part associated to @p fe_local and the
+// skeleton part @p f, which is why we can extract this information in the
+// constructor and store it once for all cells that we work on. Had we not
+// stored this information, we would be forced to assemble a large number of
+// zero terms on each cell, which would significantly slow the program.
template <int dim>
-struct Step51<dim>::ScratchData
+struct HDG<dim>::ScratchData
{
FEValues<dim> fe_values_local;
FEFaceValues<dim> fe_face_values_local;
RightHandSide<dim> right_hand_side;
const Solution<dim> exact_solution;
- // Full constructor
ScratchData(const FiniteElement<dim> &fe,
const FiniteElement<dim> &fe_local,
const QGauss<dim> &quadrature_formula,
}
}
- // Copy constructor
ScratchData(const ScratchData &sd)
:
fe_values_local (sd.fe_values_local.get_fe(),
fe_local_support_on_face(sd.fe_local_support_on_face),
fe_support_on_face(sd.fe_support_on_face)
{}
+};
- // We manually reset our matrices to zero in the assembly process,
- // since certain matrices are only used in the reconstruction process.
- // We therefore do not implement an methods in <code>reset()</code>, but
- // need to have it for the WorkStream interface.
- void reset() {}
-};
+// @sect4{HDG::PostProcessScratchData}
+// @p PostProcessScratchData contains the data used by <code>WorkStream</code>
+// when post-processing the local solution $u^*$. It is similar, but much
+// simpler, than @p ScratchData.
template <int dim>
-struct Step51<dim>::PostProcessScratchData
+struct HDG<dim>::PostProcessScratchData
{
FEValues<dim> fe_values_local;
FEValues<dim> fe_values;
Vector<double> cell_rhs;
Vector<double> cell_sol;
- // Full constructor
PostProcessScratchData(const FiniteElement<dim> &fe,
const FiniteElement<dim> &fe_local,
const QGauss<dim> &quadrature_formula,
cell_sol (fe.dofs_per_cell)
{}
- // Copy constructor
PostProcessScratchData(const PostProcessScratchData &sd)
:
fe_values_local (sd.fe_values_local.get_fe(),
cell_rhs (sd.cell_rhs),
cell_sol (sd.cell_sol)
{}
+};
- void reset()
- {
- cell_matrix = 0.;
- cell_rhs = 0.;
- cell_sol = 0.;
- }
-};
+// @sect4{HDG::copy_local_to_global}
+// If we are in the first step of the solution, i.e. @p trace_reconstruct=false,
+// then we assemble the global system.
template <int dim>
-void Step51<dim>::copy_local_to_global(const PerTaskData &data)
+void HDG<dim>::copy_local_to_global(const PerTaskData &data)
{
if (data.trace_reconstruct == false)
constraints.distribute_local_to_global (data.cell_matrix,
system_matrix, system_rhs);
}
-template <int dim>
-void
-Step51<dim>::setup_system ()
-{
- dof_handler_local.distribute_dofs(fe_local);
- dof_handler.distribute_dofs(fe);
- dof_handler_u_post.distribute_dofs(fe_u_post);
-
- std::cout << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
-
- solution.reinit (dof_handler.n_dofs());
- system_rhs.reinit (dof_handler.n_dofs());
-
- solution_local.reinit (dof_handler_local.n_dofs());
- solution_u_post.reinit (dof_handler_u_post.n_dofs());
-
- constraints.clear ();
- DoFTools::make_hanging_node_constraints (dof_handler, constraints);
- typename FunctionMap<dim>::type boundary_functions;
- Solution<dim> solution_function;
- boundary_functions[0] = &solution_function;
- VectorTools::project_boundary_values (dof_handler,
- boundary_functions,
- QGauss<dim-1>(fe.degree+1),
- constraints);
- constraints.close ();
-
- {
- CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dof_handler, csp,
- constraints, false);
- sparsity_pattern.copy_from(csp, fe.dofs_per_face);
- }
- system_matrix.reinit (sparsity_pattern);
-}
-
+// @sect4{HDG::assemble_system}
+// The @p assemble_system function is similar to <code>Step-32</code>, where
+// the quadrature formula and the update flags are set up, and then
+// <code>WorkStream</code> is used to do the work in a multi-threaded manner.
+// The @p trace_reconstruct input parameter is used to decide whether we are
+// solving for the local solution (true) or the global skeleton solution
+// (false).
template <int dim>
void
-Step51<dim>::assemble_system (const bool trace_reconstruct)
+HDG<dim>::assemble_system (const bool trace_reconstruct)
{
const QGauss<dim> quadrature_formula(fe.degree+1);
const QGauss<dim-1> face_quadrature_formula(fe.degree+1);
WorkStream::run(dof_handler.begin_active(),
dof_handler.end(),
*this,
- &Step51<dim>::assemble_system_one_cell,
- &Step51<dim>::copy_local_to_global,
+ &HDG<dim>::assemble_system_one_cell,
+ &HDG<dim>::copy_local_to_global,
scratch,
task_data);
}
+
+// @sect4{HDG::assemble_system_one_cell}
+// The real work of the HDG program is done by @p assemble_system_one_cell.
+// Assembling the local matrices $A, B, C$ is done here, along with the
+// local contributions of the global matrix $D$.
template <int dim>
void
-Step51<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData &scratch,
PerTaskData &task_data)
{
- // Construct iterator for dof_handler_local
+// Construct iterator for dof_handler_local for FEValues reinit function.
typename DoFHandler<dim>::active_cell_iterator
loc_cell (&triangulation,
cell->level(),
const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size();
const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size();
- // const unsigned int dofs_per_cell = scratch.fe_face_values.get_fe().dofs_per_cell;
const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell;
- // Choose stabilization parameter to be 5 * diffusion = 5
- const double tau_stab_diffusion = 5.;
-
const FEValuesExtractors::Vector fluxes (0);
const FEValuesExtractors::Scalar scalar (dim);
{
scratch.lf_matrix = 0;
scratch.fl_matrix = 0;
- task_data.reset();
+ task_data.cell_matrix = 0;
+ task_data.cell_vector = 0;
}
scratch.fe_values_local.reinit (loc_cell);
+ // We first compute the cell-interior contribution to @p ll_matrix matrix
+ // (referred to as matrix $A$ in the introduction) corresponding to
+ // local-local coupling, as well as the local right-hand-side vector. We
+ // store the values at each quadrature point for the basis functions, the
+ // right-hand-side value, and the convection velocity.
for (unsigned int q=0; q<n_q_points; ++q)
{
const double rhs_value
}
}
+ // Face terms are assembled on all faces of all elements. This is in
+ // contrast to more traditional DG methods, where each face is only visited
+ // once in the assembly procedure.
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
{
scratch.fe_face_values_local.reinit(loc_cell, face);
scratch.fe_face_values.reinit(cell, face);
+
+ // The already obtained $\hat{u}$ values are needed when solving for the
+ // local variables.
if (task_data.trace_reconstruct)
scratch.fe_face_values.get_function_values (solution, scratch.trace_values);
const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
const Tensor<1,dim> convection
= scratch.convection_velocity.value(quadrature_point);
- const double tau_stab = (tau_stab_diffusion +
+
+ // Here we compute the stabilization parameter discussed in the
+ // introduction: since the diffusion is one and the diffusion length
+ // scale is set to 1/5, it simply results in a contribution of 5 for
+ // the diffusion part and the magnitude of convection through the
+ // element boundary in a centered-like scheme for the convection
+ // part.
+ const double tau_stab = (5. +
std::abs(convection * normal));
+ // We store the non-zero flux and scalar values, making use of the
+ // support_on_face information we calculated in @p ScratchData.
for (unsigned int k=0; k<scratch.fe_local_support_on_face[face].size(); ++k)
{
const unsigned int kk=scratch.fe_local_support_on_face[face][k];
scratch.u_phi[k] = scratch.fe_face_values_local[scalar].value(kk,q);
}
+ // When @p trace_reconstruct=false, we are preparing assemble the
+ // system for the skeleton variable $\lambda$. If this is the case,
+ // we must assemble all local matrices associated with the problem:
+ // local-local, local-face, face-local, and face-face. The
+ // face-face matrix is stored as @p TaskData::cell_matrix, so that
+ // it can be assembled into the global system by @p
+ // copy_local_to_global.
if (!task_data.trace_reconstruct)
{
for (unsigned int k=0; k<scratch.fe_support_on_face[face].size(); ++k)
tau_stab) * scratch.u_phi[i])
* scratch.tr_phi[j]
) * JxW;
+
+ // Note the sign of the face-local matrix. We negate the
+ // sign during assembly here so that we can use the
+ // FullMatrix::mmult with addition when computing the
+ // Schur complement.
scratch.fl_matrix(jj,ii) -= (
(scratch.q_phi[i] * normal
+
}
}
+ // This last term adds the contribution of the term $\left<w,\tau
+ // u_h\right>_{\partial \mathcal T}$ to the local matrix. As opposed
+ // to the face matrices above, we need it in both assembly stages.
for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
for (unsigned int j=0; j<scratch.fe_local_support_on_face[face].size(); ++j)
{
scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
}
- // compute the local right hand side contributions from trace
+ // When @p trace_reconstruct=true, we are solving for the local
+ // solutions on an element by element basis. The local
+ // right-hand-side is calculated by replacing the basis functions @p
+ // tr_phi in the @p lf_matrix computation by the computed values @p
+ // trace_values. Of course, the sign of the matrix is now minus
+ // since we have moved everything to the other side of the equation.
if (task_data.trace_reconstruct)
for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
{
}
}
+ // Once assembly of all of the local contributions is complete, we must either:
+ // (1) assemble the global system, or (2) compute the local solution values and
+ // save them.
+ // In either case, the first step is to invert the local-local matrix.
scratch.ll_matrix.gauss_jordan();
+
+ // For (1), we compute the Schur complement and add it to the @p
+ // cell_matrix, matrix $D$ in the introduction.
if (task_data.trace_reconstruct == false)
{
scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix);
scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true);
cell->get_dof_indices(task_data.dof_indices);
}
+ // For (2), we are simply solving (ll_matrix).(solution_local) = (l_rhs).
+ // Hence, we multiply @p l_rhs by our already inverted local-local matrix
+ // and store the result using the <code>set_dof_values</code> function.
else
{
scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs);
+// @sect4{HDG::solve}
+// The skeleton solution is solved for by using a BiCGStab solver with
+// identity preconditioner.
template <int dim>
-void Step51<dim>::solve ()
+void HDG<dim>::solve ()
{
SolverControl solver_control (system_matrix.m()*10,
1e-11*system_rhs.l2_norm());
constraints.distribute(solution);
- // update local values
+ // Once we have solved for the skeleton solution,
+ // we can solve for the local solutions in an element-by-element
+ // fashion. We do this by re-using the same @p assemble_system function
+ // but switching @p trace_reconstruct to true.
assemble_system(true);
}
+ // @sect4{HDG::postprocess}
+
+ // The postprocess method serves two purposes. First, we want to construct a
+ // post-processed scalar variables in the element space of degree $p+1$ that
+ // we hope will converge at order $p+2$. This is again an element-by-element
+ // process and only involves the scalar solution as well as the gradient on
+ // the local cell. To do this, we introduce the already defined scratch data
+ // together with some update flags and run the work stream to do this in
+ // parallel.
+ //
+ // Secondly, we want to compute discretization errors just as we did in
+ // step-7. The overall procedure is similar with calls to
+ // VectorTools::integrate_difference. The difference is in how we compute
+ // the errors for the scalar variable and the gradient variable. In step-7,
+ // we did this by computing @p L2_norm or @p H1_seminorm
+ // contributions. Here, we have a DoFHandler with these two contributions
+ // computed and sorted by their vector component, <code>[0, dim)</code> for the
+ // gradient and @p dim for the scalar. To compute their value, we hence use
+ // a ComponentSelectFunction with either of them, together with the @p
+ // SolutionAndGradient class introduced above that contains the analytic
+ // parts of either of them. Eventually, we also compute the L2-error of the
+ // post-processed solution and add the results into the convergence table.
template <int dim>
void
-Step51<dim>::postprocess()
+HDG<dim>::postprocess()
{
- // construct post-processed solution with (hopefully) higher order of
- // accuracy
{
const QGauss<dim> quadrature_formula(fe_u_post.degree+1);
const UpdateFlags local_flags (update_values);
WorkStream::run(dof_handler_u_post.begin_active(),
dof_handler_u_post.end(),
- std_cxx1x::bind (&Step51<dim>::postprocess_one_cell,
+ std_cxx1x::bind (&HDG<dim>::postprocess_one_cell,
std_cxx1x::ref(*this),
std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
std_cxx1x::function<void(const unsigned int&)>(),
0U);
}
- // Compute some convergence rates, etc., and add to a table
Vector<float> difference_per_cell (triangulation.n_active_cells());
ComponentSelectFunction<dim> value_select (dim, dim+1);
convergence_table.add_value("val L2-post", post_error);
}
+
+
+ // @sect4{HDG::postprocess_one_cell}
+ //
+ // This is the actual work done for the postprocessing. According to the
+ // discussion in the introduction, we need to set up a system that projects
+ // the gradient part of the DG solution onto the gradient of the
+ // post-processed variable. Moreover, we need to set the average of the new
+ // post-processed variable to be equal the average of the scalar DG solution
+ // on the cell.
+ //
+ // More technically speaking, the projection of the gradient is a system
+ // that would potentially fills our @p dofs_per_cell times @p dofs_per_cell
+ // matrix but is singular (the sum of all rows would be zero because the
+ // constant function has zero gradient). Therefore, we take one row away and
+ // use it for imposing the average of the scalar value. We pick the first
+ // row for the scalar part, even though we could pick any row for $\mathcal
+ // Q_{-p}$ elements. However, had we used FE_DGP elements instead, the first
+ // row would correspond to the constant part already and deleting e.g. the
+ // last row would give us a singular system. This way, our program can also
+ // be used for those elements.
template <int dim>
void
-Step51<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
PostProcessScratchData &scratch,
unsigned int &)
{
- scratch.reset();
-
typename DoFHandler<dim>::active_cell_iterator
loc_cell (&triangulation,
cell->level(),
scratch.cell_rhs(0) = sum;
}
+ // Having assembled all terms, we can again go on and solve the linear
+ // system. We again invert the matrix and then multiply the inverse by the
+ // right hand side. An alternative (and more numerically stable) would have
+ // been to only factorize the matrix and apply the factorization.
scratch.cell_matrix.gauss_jordan();
scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs);
cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
}
+
+// @sect4{HDG::output_results}
+// We have 3 sets of results that we would like to output: the local solution,
+// the post-processed local solution, and the skeleton solution. The former 2
+// both `live' on element volumes, wheras the latter lives on codimention-1 surfaces
+// of the triangulation. Our @p output_results function writes all local solutions
+// to the same vtk file, even though they correspond to different <code>DoFHandler</code>
+// objects. The graphical output for the skeleton variable is done through
+// use of the <code>DataOutFaces</code> class.
template <int dim>
-void Step51<dim>::output_results (const unsigned int cycle)
+void HDG<dim>::output_results (const unsigned int cycle)
{
std::string filename;
switch (refinement_mode)
std::ofstream output (filename.c_str());
DataOut<dim> data_out;
+
+// We first define the names and types of the local solution,
+// and add the data to @p data_out.
std::vector<std::string> names (dim, "gradient");
names.push_back ("solution");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
component_interpretation
- (dim+1, DataComponentInterpretation::component_is_part_of_vector);
+ (dim+1, DataComponentInterpretation::component_is_part_of_vector);
component_interpretation[dim]
- = DataComponentInterpretation::component_is_scalar;
+ = DataComponentInterpretation::component_is_scalar;
data_out.add_data_vector (dof_handler_local, solution_local,
names, component_interpretation);
- // Post-processed solution: can now add more than 1 dof_handler to
- // the DataOut object!
+// The second data item we add is the post-processed solution.
+// In this case, it is a single scalar variable belonging to
+// a different DoFHandler.
std::vector<std::string> post_name(1,"u_post");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
post_comp_type(1, DataComponentInterpretation::component_is_scalar);
face_out += ".vtk";
std::ofstream face_output (face_out.c_str());
+// The <code>DataOutFaces</code> class works analagously to the <code>DataOut</code>
+// class when we have a <code>DoFHandler</code> that defines the solution on
+// the skeleton of the triangulation. We treat it as such here, and the code is
+// similar to that above.
DataOutFaces<dim> data_out_face(false);
std::vector<std::string> face_name(1,"lambda");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_out_face.write_vtk (face_output);
}
+// @sect4{HDG::refine_grid}
+// We implement two different refinement cases for HDG, just as in
+// <code>Step-7</code>: adaptive_refinement and global_refinement. The
+// global_refinement option recreates the entire triangulation every
+// time. This is because we want to use a finer sequence of meshes than what
+// we would get with one refinement step, namely 2, 3, 4, 6, 8, 12, 16, ...
+// elements per direction.
+// The adaptive_refinement mode uses the <code>KellyErrorEstimator</code> to
+// give a decent indication of the non-regular regions in the scalar local
+// solutions.
template <int dim>
-void Step51<dim>::refine_grid (const unsigned int cycle)
+void HDG<dim>::refine_grid (const unsigned int cycle)
{
if (cycle == 0)
{
endc = triangulation.end();
for (; cell!=endc; ++cell)
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12)
- ||
- (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12))
- cell->face(face)->set_boundary_indicator (1);
+ if (cell->face(face)->at_boundary())
+ if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12)
+ ||
+ (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12))
+ cell->face(face)->set_boundary_indicator (1);
}
-
-
-
-
+// @sect4{HDG::run}
+// The functionality here is basically the same as <code>Step-7</code>.
+// We loop over 10 cycles, refining the grid on each one. At the end,
+// convergence tables are created.
template <int dim>
-void Step51<dim>::run ()
+void HDG<dim>::run ()
{
for (unsigned int cycle=0; cycle<10; ++cycle)
{
convergence_table.set_precision("val L2-post", 3);
convergence_table.set_scientific("val L2-post", true);
+ // There is one minor change for the convergence table compared to step-7:
+ // Since we did not refine our mesh by a factor two in each cycle (but
+ // rather used the sequence 2, 3, 4, 6, 8, 12, ...), we need to tell the
+ // convergence rate evaluation about this. We do this by setting the number
+ // of cells as a reference column and additionally specifying the dimension
+ // of the problem, which gives the computation the necessary information for
+ // how much the mesh was refinement given a certain increase in the number
+ // of cells.
if (refinement_mode == global_refinement)
{
convergence_table
convergence_table.write_text(std::cout);
}
+} // end of namespace Step51
int main (int argc, char **argv)
{
<< "=============================================" << std::endl
<< std::endl;
- Step51<dim> hdg_problem (1, Step51<dim>::adaptive_refinement);
+ Step51::HDG<dim> hdg_problem (1, Step51::HDG<dim>::adaptive_refinement);
hdg_problem.run ();
std::cout << std::endl;
<< "===========================================" << std::endl
<< std::endl;
- Step51<dim> hdg_problem (1, Step51<dim>::global_refinement);
+ Step51::HDG<dim> hdg_problem (1, Step51::HDG<dim>::global_refinement);
hdg_problem.run ();
std::cout << std::endl;
<< "===========================================" << std::endl
<< std::endl;
- Step51<dim> hdg_problem (3, Step51<dim>::global_refinement);
+ Step51::HDG<dim> hdg_problem (3, Step51::HDG<dim>::global_refinement);
hdg_problem.run ();
std::cout << std::endl;
#include <deal.II/grid/tria.h>
#include <deal.II/base/std_cxx1x/function.h>
+#include <deal.II/base/std_cxx1x/tuple.h>
+#include <set>
#include <vector>
#include <list>
#include <utility>
const std::vector<types::global_dof_index> &
get_p4est_tree_to_coarse_cell_permutation() const;
+
+ /**
+ * Join faces in the p4est forest due to periodic boundary conditions.
+ *
+ * The entries in the std::vector should have the form
+ * std_cxx1x::tuple<cell1, face_no1, cell2, face_no2>.
+ *
+ * The vector can be filled by the function
+ * DoFTools::identify_periodic_face_pairs.
+ *
+ * @note Before this function can be used the triangulation has to be
+ * initialized and must not be refined.
+ * Calling this function more than once is possible, but not recommended:
+ * The function destroys and rebuilds the p4est forest each time it is called.
+ */
+ void
+ add_periodicity
+ (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >&);
+
+
private:
/**
* MPI communicator to be
* triangulation.
*/
typename dealii::internal::p4est::types<dim>::forest *parallel_forest;
+ /**
+ * A data structure that holds some
+ * information about the ghost cells of the triangulation.
+ */
+ typename dealii::internal::p4est::types<dim>::ghost *parallel_ghost;
/**
* A flag that indicates
*/
void attach_mesh_data();
+ /**
+ * fills a map that, for each vertex, lists all the processors whose
+ * subdomains are adjacent to that vertex. Used by
+ * DoFHandler::Policy::ParallelDistributed.
+ */
+ void
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors);
template <int, int> friend class dealii::internal::DoFHandler::Policy::ParallelDistributed;
};
* in hierarchical ordering is the ith deal cell starting
* from begin(0).
*/
- const std::vector<unsigned int> &
+ const std::vector<types::global_dof_index> &
get_p4est_tree_to_coarse_cell_permutation() const;
/**
* these variables at a
* couple places anyway.
*/
- std::vector<unsigned int> coarse_cell_to_p4est_tree_permutation;
- std::vector<unsigned int> p4est_tree_to_coarse_cell_permutation;
+ std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
+ std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
/**
* dummy settings
*/
Settings settings;
+ /**
+ * Like above, this method, which is only implemented for dim = 2 or 3,
+ * needs a stub because it is used in dof_handler_policy.cc
+ */
+ void
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors);
+
};
}
}
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_levels.h>
#include <deal.II/dofs/dof_faces.h>
-#include <deal.II/hp/dof_levels.h>
+#include <deal.II/hp/dof_level.h>
#include <deal.II/hp/dof_faces.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_iterator.templates.h>
template <int dim, int spacedim, bool lda>
static
void
- update_cell_dof_indices_cache (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim>, lda> &)
+ update_cell_dof_indices_cache (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim>, lda> &accessor)
{
-//TODO[WB]: should implement a dof indices cache for hp as well
+ // caches are only for cells with DoFs, i.e., for active ones
+ if (accessor.has_children())
+ return;
- // not implemented, but should also
- // not be called
- Assert (false, ExcNotImplemented());
+ const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
+
+ // make sure the cache is at least
+ // as big as we need it when
+ // writing to the last element of
+ // this cell
+ Assert (static_cast<unsigned int>(accessor.present_index)
+ <
+ accessor.dof_handler->levels[accessor.present_level]
+ ->cell_cache_offsets.size(),
+ ExcInternalError());
+ Assert (accessor.dof_handler->levels[accessor.present_level]
+ ->cell_cache_offsets[accessor.present_index]
+ <=
+ accessor.dof_handler->levels[accessor.present_level]
+ ->cell_dof_indices_cache.size(),
+ ExcInternalError());
+
+ std::vector<types::global_dof_index> dof_indices (dofs_per_cell);
+ static_cast<const dealii::DoFAccessor<dim,dealii::hp::DoFHandler<dim,spacedim>,lda> &>
+ (accessor).get_dof_indices (dof_indices, accessor.active_fe_index());
+
+ types::global_dof_index *next_dof_index
+ = &accessor.dof_handler->levels[accessor.present_level]
+ ->cell_dof_indices_cache[accessor.dof_handler->levels[accessor.present_level]
+ ->cell_cache_offsets[accessor.present_index]];
+ for (unsigned int i=0; i<dofs_per_cell; ++i, ++next_dof_index)
+ *next_dof_index = dof_indices[i];
}
+
+
/**
* Implement setting dof
* indices on a
}
- /**
- * A function that collects the
- * global indices of degrees of
- * freedom. This function works
- * for ::DoFHandler and all
- * template arguments and copies
- * the data out of the cache that
- * we hold for each cell.
- */
- template <int dim, int spacedim, bool lda>
- static
- void
- get_dof_indices (const DoFCellAccessor<DoFHandler<dim,spacedim>, lda> &accessor,
- std::vector<types::global_dof_index> &dof_indices)
- {
- AssertDimension (dof_indices.size(), accessor.get_fe().dofs_per_cell);
-
- Assert (!accessor.has_children(),
- ExcMessage ("Cell must be active."));
-
- types::global_dof_index *cache = &accessor.dof_handler->levels[accessor.level()]
- ->cell_dof_indices_cache[accessor.present_index * accessor.get_fe().dofs_per_cell];
- for (unsigned int i=0; i<accessor.get_fe().dofs_per_cell; ++i, ++cache)
- dof_indices[i] = *cache;
- }
-
- /**
- * Same function as above except
- * that it works for
- * hp::DoFHandler objects that do
- * not have a cache for the local
- * DoF indices.
- */
- template <int dim, int spacedim, bool lda>
- static
- void
- get_dof_indices (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim>, lda> &accessor,
- std::vector<types::global_dof_index> &dof_indices)
- {
- // no caching for
- // hp::DoFHandler
- // implemented
- typedef
- dealii::DoFAccessor<dim,dealii::hp::DoFHandler<dim,spacedim>, lda>
- DoFAccessor;
- accessor.DoFAccessor::get_dof_indices (dof_indices,
- accessor.active_fe_index());
- }
-
- /**
- * A function that collects the
- * values of degrees of freedom. This
- * function works for ::DoFHandler
- * and all template arguments and
- * uses the data from the cache of
- * indices that we hold for each
- * cell.
- */
- template <int dim, int spacedim, bool lda, class InputVector, typename ForwardIterator>
- static
- void
- get_dof_values (const DoFCellAccessor<DoFHandler<dim,spacedim>, lda> &accessor,
- const InputVector &values,
- ForwardIterator local_values_begin,
- ForwardIterator local_values_end)
- {
- typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim>, lda> BaseClass;
- Assert (static_cast<unsigned int>(local_values_end-local_values_begin)
- == accessor.get_fe().dofs_per_cell,
- typename BaseClass::ExcVectorDoesNotMatch());
- Assert (values.size() == accessor.get_dof_handler().n_dofs(),
- typename BaseClass::ExcVectorDoesNotMatch());
-
- Assert (!accessor.has_children(),
- ExcMessage ("Cell must be active."));
-
- types::global_dof_index *cache = &accessor.dof_handler->levels[accessor.level()]
- ->cell_dof_indices_cache[accessor.present_index *
- accessor.get_fe().dofs_per_cell];
-
- values.extract_subvector_to (cache,
- cache + accessor.get_fe().dofs_per_cell,
- local_values_begin);
- }
-
- /**
- * Same function as above except
- * that it works for
- * hp::DoFHandler objects that do
- * not have a cache for the local
- * DoF indices.
- */
- template <int dim, int spacedim, bool lda, class InputVector, typename ForwardIterator>
- static
- void
- get_dof_values (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim>, lda> &accessor,
- const InputVector &values,
- ForwardIterator local_values_begin,
- ForwardIterator local_values_end)
- {
- typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim>, lda> BaseClass;
-
- // no caching for hp::DoFHandler
- // implemented
- Assert (static_cast<unsigned int>(local_values_end-local_values_begin)
- == accessor.get_fe().dofs_per_cell,
- typename BaseClass::ExcVectorDoesNotMatch());
- const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
-
- std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
- get_dof_indices (accessor, local_dof_indices);
-
- values.extract_subvector_to (local_dof_indices.begin(),
- local_dof_indices.end(),
- local_values_begin);
- }
-
-
- /**
- * A function that collects the
- * values of degrees of freedom. This
- * function works for ::DoFHandler
- * and all template arguments and
- * uses the data from the cache of
- * indices that we hold for each
- * cell.
- */
- template <int dim, int spacedim, bool lda, class InputVector, typename ForwardIterator>
- static
- void
- get_dof_values (const DoFCellAccessor<DoFHandler<dim,spacedim>, lda> &accessor,
- const ConstraintMatrix &constraints,
- const InputVector &values,
- ForwardIterator local_values_begin,
- ForwardIterator local_values_end)
- {
- typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim>, lda> BaseClass;
- Assert (static_cast<unsigned int>(local_values_end-local_values_begin)
- == accessor.get_fe().dofs_per_cell,
- typename BaseClass::ExcVectorDoesNotMatch());
- Assert (values.size() == accessor.get_dof_handler().n_dofs(),
- typename BaseClass::ExcVectorDoesNotMatch());
-
- Assert (!accessor.has_children(),
- ExcMessage ("Cell must be active."));
-
- types::global_dof_index *cache = &accessor.dof_handler->levels[accessor.level()]
- ->cell_dof_indices_cache[accessor.present_index *
- accessor.get_fe().dofs_per_cell];
- constraints.get_dof_values(values, *cache, local_values_begin,
- local_values_end);
- }
-
- /**
- * Same function as above except
- * that it works for
- * hp::DoFHandler objects that do
- * not have a cache for the local
- * DoF indices.
- */
- template <int dim, int spacedim, bool lda, class InputVector, typename ForwardIterator>
- static
- void
- get_dof_values (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim>, lda> &accessor,
- const ConstraintMatrix &constraints,
- const InputVector &values,
- ForwardIterator local_values_begin,
- ForwardIterator local_values_end)
- {
- // no caching for hp::DoFHandler
- // implemented
- typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim>, lda> BaseClass;
- Assert (static_cast<unsigned int>(local_values_end-local_values_begin)
- == accessor.get_fe().dofs_per_cell,
- typename BaseClass::ExcVectorDoesNotMatch());
- const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
-
- std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
- get_dof_indices (accessor, local_dof_indices);
-
- constraints.get_dof_values (values, local_dof_indices.begin(),
- local_values_begin, local_values_end);
- }
-
-
- /**
- * Same set of functions as above
- * except that it sets rather than
- * gets values
- */
- template <int dim, int spacedim, bool lda, class OutputVector, typename number>
- static
- void
- set_dof_values (const DoFCellAccessor<DoFHandler<dim,spacedim>, lda> &accessor,
- const dealii::Vector<number> &local_values,
- OutputVector &values)
- {
- typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim>, lda> BaseClass;
- Assert (local_values.size() == accessor.get_fe().dofs_per_cell,
- typename BaseClass::ExcVectorDoesNotMatch());
- Assert (values.size() == accessor.get_dof_handler().n_dofs(),
- typename BaseClass::ExcVectorDoesNotMatch());
-
- Assert (!accessor.has_children(),
- ExcMessage ("Cell must be active."));
-
- types::global_dof_index *cache = &accessor.dof_handler->levels[accessor.level()]
- ->cell_dof_indices_cache[accessor.present_index * accessor.get_fe().dofs_per_cell];
- for (unsigned int i=0; i<accessor.get_fe().dofs_per_cell; ++i, ++cache)
- values(*cache) = local_values(i);
- }
-
-
-
- template <int dim, int spacedim, bool lda, class OutputVector, typename number>
- static
- void
- set_dof_values (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim>, lda> &accessor,
- const dealii::Vector<number> &local_values,
- OutputVector &values)
- {
- // no caching for hp::DoFHandler
- // implemented
- const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
-
- std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
- get_dof_indices (accessor, local_dof_indices);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- values(local_dof_indices[i]) = local_values(i);
- }
-
/**
* Do what the active_fe_index
ExcMessage ("Can't ask for DoF indices on artificial cells."));
AssertDimension (dof_indices.size(), this->get_fe().dofs_per_cell);
- dealii::internal::DoFCellAccessor::Implementation::get_dof_indices (*this, dof_indices);
+ const types::global_dof_index *cache
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+ for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
+ dof_indices[i] = *cache;
}
DoFCellAccessor<DH,lda>::get_dof_values (const InputVector &values,
Vector<number> &local_values) const
{
- Assert (this->is_artificial() == false,
- ExcMessage ("Can't ask for DoF indices on artificial cells."));
- AssertDimension (local_values.size(), this->get_fe().dofs_per_cell);
-
- dealii::internal::DoFCellAccessor::Implementation
- ::get_dof_values (*this, values, local_values.begin(), local_values.end());
+ get_dof_values (values, local_values.begin(), local_values.end());
}
{
Assert (this->is_artificial() == false,
ExcMessage ("Can't ask for DoF indices on artificial cells."));
- dealii::internal::DoFCellAccessor::Implementation
- ::get_dof_values (*this, values, local_values_begin, local_values_end);
+ Assert (!this->has_children(),
+ ExcMessage ("Cell must be active."));
+
+ Assert (static_cast<unsigned int>(local_values_end-local_values_begin)
+ == this->get_fe().dofs_per_cell,
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
+ Assert (values.size() == this->get_dof_handler().n_dofs(),
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
+
+ const types::global_dof_index *cache
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+
+ values.extract_subvector_to (cache,
+ cache + this->get_fe().dofs_per_cell,
+ local_values_begin);
}
{
Assert (this->is_artificial() == false,
ExcMessage ("Can't ask for DoF indices on artificial cells."));
- dealii::internal::DoFCellAccessor::Implementation
- ::get_dof_values (*this, constraints, values,
- local_values_begin, local_values_end);
+ Assert (!this->has_children(),
+ ExcMessage ("Cell must be active."));
+
+ Assert (static_cast<unsigned int>(local_values_end-local_values_begin)
+ == this->get_fe().dofs_per_cell,
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
+ Assert (values.size() == this->get_dof_handler().n_dofs(),
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
+
+
+ const types::global_dof_index *cache
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+
+ constraints.get_dof_values(values, *cache, local_values_begin,
+ local_values_end);
}
{
Assert (this->is_artificial() == false,
ExcMessage ("Can't ask for DoF indices on artificial cells."));
- dealii::internal::DoFCellAccessor::Implementation
- ::set_dof_values (*this, local_values, values);
+ Assert (!this->has_children(),
+ ExcMessage ("Cell must be active."));
+
+ Assert (static_cast<unsigned int>(local_values.size())
+ == this->get_fe().dofs_per_cell,
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
+ Assert (values.size() == this->get_dof_handler().n_dofs(),
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
+
+
+ const types::global_dof_index *cache
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+
+ for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
+ values(*cache) = local_values(i);
}
*/
DoFObjects<dim> dof_object;
+ /**
+ * Return a pointer to the beginning of the DoF indices cache
+ * for a given cell.
+ *
+ * @param obj_index The number of the cell we are looking at.
+ * @param dofs_per_cell The number of DoFs per cell for this cell.
+ * @return A pointer to the first DoF index for the current cell. The
+ * next dofs_per_cell indices are for the current cell.
+ */
+ const types::global_dof_index *
+ get_cell_cache_start (const unsigned int obj_index,
+ const unsigned int dofs_per_cell) const;
/**
* Determine an estimate for the
};
+
+ template <int dim>
+ inline
+ const types::global_dof_index *
+ DoFLevel<dim>::get_cell_cache_start (const unsigned int obj_index,
+ const unsigned int dofs_per_cell) const
+ {
+ Assert (obj_index*dofs_per_cell+dofs_per_cell
+ <=
+ cell_dof_indices_cache.size(),
+ ExcInternalError());
+
+ return &cell_dof_indices_cache[obj_index*dofs_per_cell];
+ }
+
+
+
template <int dim>
inline
std::size_t
* More information on the topic can be found in the
* @ref GlossFaceOrientation "glossary" article.
*
- * @note This function will not work for DoFHandler objects that are
- * built on a parallel::distributed::Triangulation object unless both
- * faces (or their children) are owned by the current processor.
+ * @note For DoFHandler objects that are built on a
+ * parallel::distributed::Triangulation object
+ * parallel::distributed::Triangulation::add_periodicity has to be called
+ * before.
*
* @author Matthias Maier, 2012
*/
* function will be used for which the respective flag was set in the
* component mask.
*
- * @note This function will not work for DoFHandler objects that are
- * built on a parallel::distributed::Triangulation object.
+ * @note For DoFHandler objects that are built on a
+ * parallel::distributed::Triangulation object
+ * parallel::distributed::Triangulation::add_periodicity has to be called
+ * before.
*
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*
* @p orthogonal_equality. This can be used to implement conditions such
* as $u(0,y)=u(1,y+1)$.
*
- * @note This function will not work for DoFHandler objects that are
- * built on a parallel::distributed::Triangulation object.
+ * @note For DoFHandler objects that are built on a
+ * parallel::distributed::Triangulation object
+ * parallel::distributed::Triangulation::add_periodicity has to be called
+ * before.
*
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*
* meshes with cells not in @ref GlossFaceOrientation
* "standard orientation".
*
- * @note This function will not work for DoFHandler objects that are
- * built on a parallel::distributed::Triangulation object.
+ * @note For DoFHandler objects that are built on a
+ * parallel::distributed::Triangulation object
+ * parallel::distributed::Triangulation::add_periodicity has to be called
+ * before.
*
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*/
* meshes with cells not in @ref GlossFaceOrientation
* "standard orientation".
*
- * @note This function will not work for DoFHandler objects that are
- * built on a parallel::distributed::Triangulation object.
+ * @note For DoFHandler objects that are built on a
+ * parallel::distributed::Triangulation object
+ * parallel::distributed::Triangulation::add_periodicity has to be called
+ * before.
*
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*/
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/std_cxx1x/tuple.h>
#include <bitset>
#include <list>
const int direction,
const dealii::Tensor<1,DH::space_dimension> &offset);
+ /**
+ * Add periodicity information to the @p periodicity_vector for the
+ * boundaries with boundary id @p b_id1 and @p b_id2 in cartesian
+ * direction @p direction.
+ *
+ * This function tries to match all faces belonging to the first
+ * boundary with faces belonging to the second boundary
+ * by comparing the center of the cell faces. To find the correct
+ * corresponding faces, the direction argument indicates in which
+ * cartesian direction periodicity should be set.
+ * ((0,1,2) -> (x,y,z)-direction)
+ *
+ * The output of this function can be used in
+ * parallel::distributed::Triangulation::add_periodicity
+ */
+ template<typename DH>
+ void
+ identify_periodic_face_pairs
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+ typename DH::cell_iterator, unsigned int> >
+ &periodicity_vector);
/**
#include <deal.II/base/config.h>
-#include <deal.II/dofs/dof_levels.h>
-#include <deal.II/hp/dof_objects.h>
#include <deal.II/hp/fe_collection.h>
#include <vector>
{
namespace hp
{
+ /**
+ * Store the indices of the degrees of freedom which are located on
+ * objects of dimension @p structdim < dim, i.e., for faces or edges
+ * of cells. This is opposed to the internal::hp::DoFLevels class
+ * that stores the DoF indices on cells.
+ *
+ * The things we store here is very similar to what is stored in the
+ * internal::DoFHandler::DoFObjects classes (see there for more
+ * information, in particular on the layout of the class hierarchy,
+ * and the use of file names).
+ *
+ * <h4>Offset computations</h4>
+ *
+ * For hp methods, not all cells may use the same finite element, and
+ * it is consequently more complicated to determine where the DoF
+ * indices for a given line, quad, or hex are stored. As described in
+ * the documentation of the internal::DoFHandler::DoFLevel class, we
+ * can compute the location of the first line DoF, for example, by
+ * calculating the offset as <code>line_index *
+ * dof_handler.get_fe().dofs_per_line</code>. This of course doesn't
+ * work any more if different lines may have different numbers of
+ * degrees of freedom associated with them. Consequently, rather than
+ * using this simple multiplication, the dofs array has an associated
+ * array dof_offsets. The data corresponding to a
+ * line then starts at index <code>line_dof_offsets[line_index]</code>
+ * within the <code>line_dofs</code> array.
+ *
+ *
+ * <h4>Multiple data sets per object</h4>
+ *
+ * If two adjacent cells use different finite elements, then
+ * the face that they share needs to store DoF indices for both
+ * involved finite elements. While faces therefore have to have at
+ * most two sets of DoF indices, it is easy to see that edges and
+ * vertices can have as many sets of DoF indices associated with them
+ * as there are adjacent cells.
+ *
+ * Consequently, for objects that have a lower dimensionality than
+ * cells, we have to store a map from the finite element index to the
+ * set of DoF indices associated. Since real sets are typically very
+ * inefficient to store, and since most of the time we expect the
+ * number of individual keys to be small (frequently, adjacent cells
+ * will have the same finite element, and only a single entry will
+ * exist in the map), what we do is instead to store a linked list. In
+ * this format, the first entry starting at position
+ * <code>lines.dofs[lines.dof_offsets[line_index]]</code> will denote
+ * the finite element index of the set of DoF indices following; after
+ * this set, we will store the finite element index of the second set
+ * followed by the corresponding DoF indices; and so on. Finally, when
+ * all finite element indices adjacent to this object have been
+ * covered, we write a -1 to indicate the end of the list.
+ *
+ * Access to this kind of data, as well as the distinction between
+ * cells and objects of lower dimensionality are encoded in the
+ * accessor functions, DoFObjects::set_dof_index() and
+ * DoFLevel::get_dof_index(). They are able to traverse this
+ * list and pick out or set a DoF index given the finite element index
+ * and its location within the set of DoFs corresponding to this
+ * finite element.
+ *
+ *
+ * @ingroup hp
+ * @author Tobias Leicht, 2006
+ */
+ template <int structdim>
+ class DoFIndicesOnFacesOrEdges
+ {
+ public:
+ /**
+ * Store the start index for
+ * the degrees of freedom of each
+ * object in the @p dofs array.
+ *
+ * The type we store is then obviously the type the @p dofs array
+ * uses for indexing.
+ */
+ std::vector<unsigned int> dof_offsets;
+
+ /**
+ * Store the global indices of
+ * the degrees of freedom. See
+ * DoFLevel() for detailed
+ * information.
+ */
+ std::vector<types::global_dof_index> dofs;
+
+ /**
+ * Set the global index of
+ * the @p local_index-th
+ * degree of freedom located
+ * on the object with number @p
+ * obj_index to the value
+ * given by @p global_index. The @p
+ * dof_handler argument is
+ * used to access the finite
+ * element that is to be used
+ * to compute the location
+ * where this data is stored.
+ *
+ * The third argument, @p
+ * fe_index, denotes which of
+ * the finite elements
+ * associated with this
+ * object we shall
+ * access. Refer to the
+ * general documentation of
+ * the internal::hp::DoFLevel
+ * class template for more
+ * information.
+ */
+ template <int dim, int spacedim>
+ void
+ set_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const types::global_dof_index global_index,
+ const unsigned int obj_level);
+
+ /**
+ * Return the global index of
+ * the @p local_index-th
+ * degree of freedom located
+ * on the object with number @p
+ * obj_index. The @p
+ * dof_handler argument is
+ * used to access the finite
+ * element that is to be used
+ * to compute the location
+ * where this data is stored.
+ *
+ * The third argument, @p
+ * fe_index, denotes which of
+ * the finite elements
+ * associated with this
+ * object we shall
+ * access. Refer to the
+ * general documentation of
+ * the internal::hp::DoFLevel
+ * class template for more
+ * information.
+ */
+ template <int dim, int spacedim>
+ types::global_dof_index
+ get_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const unsigned int obj_level) const;
+
+ /**
+ * Return the number of
+ * finite elements that are
+ * active on a given
+ * object. If this is a cell,
+ * the answer is of course
+ * one. If it is a face, the
+ * answer may be one or two,
+ * depending on whether the
+ * two adjacent cells use the
+ * same finite element or
+ * not. If it is an edge in
+ * 3d, the possible return
+ * value may be one or any
+ * other value larger than
+ * that.
+ *
+ * If the object is not part
+ * of an active cell, then no
+ * degrees of freedom have
+ * been distributed and zero
+ * is returned.
+ */
+ template <int dim, int spacedim>
+ unsigned int
+ n_active_fe_indices (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index) const;
+
+ /**
+ * Return the fe_index of the
+ * n-th active finite element
+ * on this object.
+ */
+ template <int dim, int spacedim>
+ types::global_dof_index
+ nth_active_fe_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n) const;
+
+ /**
+ * Check whether a given
+ * finite element index is
+ * used on the present
+ * object or not.
+ */
+ template <int dim, int spacedim>
+ bool
+ fe_index_is_active (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int obj_level) const;
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object.
+ */
+ std::size_t memory_consumption () const;
+ };
+
+
/**
* These classes are similar to the internal::hp::DoFLevel classes. We here store
* @author Tobias Leicht, 2006
*/
template<int dim>
- class DoFFaces;
+ class DoFIndicesOnFaces;
/**
* @author Tobias Leicht, 2006
*/
template<>
- class DoFFaces<1>
+ class DoFIndicesOnFaces<1>
{
public:
/**
* @author Tobias Leicht, 2006
*/
template<>
- class DoFFaces<2>
+ class DoFIndicesOnFaces<2>
{
public:
/**
* Indices of DoFs on the lines that bound cells.
*/
- internal::hp::DoFObjects<1> lines;
+ internal::hp::DoFIndicesOnFacesOrEdges<1> lines;
/**
* Determine an estimate for the
* @author Tobias Leicht, 2006
*/
template<>
- class DoFFaces<3>
+ class DoFIndicesOnFaces<3>
{
public:
/**
* Indices of DoFs on the lines that form the edges of cells.
*/
- internal::hp::DoFObjects<1> lines;
+ internal::hp::DoFIndicesOnFacesOrEdges<1> lines;
/**
* Indices of DoFs on the quads that bound cells.
*/
- internal::hp::DoFObjects<2> quads;
+ internal::hp::DoFIndicesOnFacesOrEdges<2> quads;
/**
* Determine an estimate for the
std::size_t memory_consumption () const;
};
+
+ // --------------------- inline and template functions ------------------
+
+ template <int structdim>
+ template <int dim, int spacedim>
+ inline
+ types::global_dof_index
+ DoFIndicesOnFacesOrEdges<structdim>::
+ get_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const unsigned int obj_level) const
+ {
+ Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
+ ExcMessage ("You need to specify a FE index when working "
+ "with hp DoFHandlers"));
+ Assert (&dof_handler != 0,
+ ExcMessage ("No DoFHandler is specified for this iterator"));
+ Assert (&dof_handler.get_fe() != 0,
+ ExcMessage ("No finite element collection is associated with "
+ "this DoFHandler"));
+ Assert (fe_index < dof_handler.get_fe().size(),
+ ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
+ Assert (local_index <
+ dof_handler.get_fe()[fe_index].template n_dofs_per_object<structdim>(),
+ ExcIndexRange(local_index, 0,
+ dof_handler.get_fe()[fe_index]
+ .template n_dofs_per_object<structdim>()));
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const types::global_dof_index starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ while (true)
+ {
+ Assert (*pointer != numbers::invalid_dof_index,
+ ExcInternalError());
+ if (*pointer == fe_index)
+ return *(pointer + 1 + local_index);
+ else
+ pointer += static_cast<types::global_dof_index>(
+ dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>() + 1);
+ }
+ }
+
+
+
+ template <int structdim>
+ template <int dim, int spacedim>
+ inline
+ void
+ DoFIndicesOnFacesOrEdges<structdim>::
+ set_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const types::global_dof_index global_index,
+ const unsigned int obj_level)
+ {
+ Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
+ ExcMessage ("You need to specify a FE index when working "
+ "with hp DoFHandlers"));
+ Assert (&dof_handler != 0,
+ ExcMessage ("No DoFHandler is specified for this iterator"));
+ Assert (&dof_handler.get_fe() != 0,
+ ExcMessage ("No finite element collection is associated with "
+ "this DoFHandler"));
+ Assert (fe_index < dof_handler.get_fe().size(),
+ ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
+ Assert (local_index <
+ dof_handler.get_fe()[fe_index].template n_dofs_per_object<structdim>(),
+ ExcIndexRange(local_index, 0,
+ dof_handler.get_fe()[fe_index]
+ .template n_dofs_per_object<structdim>()));
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const types::global_dof_index starting_offset = dof_offsets[obj_index];
+ types::global_dof_index *pointer = &dofs[starting_offset];
+ while (true)
+ {
+ Assert (*pointer != numbers::invalid_dof_index,
+ ExcInternalError());
+ if (*pointer == fe_index)
+ {
+ *(pointer + 1 + local_index) = global_index;
+ return;
+ }
+ else
+ pointer += dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>() + 1;
+ }
+ }
+
+
+
+ template <int structdim>
+ template <int dim, int spacedim>
+ inline
+ unsigned int
+ DoFIndicesOnFacesOrEdges<structdim>::
+ n_active_fe_indices (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index) const
+ {
+ Assert (&dof_handler != 0,
+ ExcMessage ("No DoFHandler is specified for this iterator"));
+ Assert (&dof_handler.get_fe() != 0,
+ ExcMessage ("No finite element collection is associated with "
+ "this DoFHandler"));
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ if (dof_offsets[obj_index] == numbers::invalid_dof_index)
+ return 0;
+
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with this
+ // object. hop along the list of index sets until we find the
+ // one with the correct fe_index, and then poke into that
+ // part. trigger an exception if we can't find a set for this
+ // particular fe_index
+ const unsigned int starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ unsigned int counter = 0;
+ while (true)
+ {
+ if (*pointer == numbers::invalid_dof_index)
+ // end of list reached
+ return counter;
+ else
+ {
+ ++counter;
+ pointer += dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>() + 1;
+ }
+ }
+ }
+
+
+
+ template <int structdim>
+ template <int dim, int spacedim>
+ inline
+ types::global_dof_index
+ DoFIndicesOnFacesOrEdges<structdim>::
+ nth_active_fe_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_level,
+ const unsigned int obj_index,
+ const unsigned int n) const
+ {
+ Assert (&dof_handler != 0,
+ ExcMessage ("No DoFHandler is specified for this iterator"));
+ Assert (&dof_handler.get_fe() != 0,
+ ExcMessage ("No finite element collection is associated with "
+ "this DoFHandler"));
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ Assert (n < n_active_fe_indices(dof_handler, obj_index),
+ ExcIndexRange (n, 0,
+ n_active_fe_indices(dof_handler, obj_index)));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const unsigned int starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ unsigned int counter = 0;
+ while (true)
+ {
+ Assert (*pointer != numbers::invalid_dof_index,
+ ExcInternalError());
+
+ const unsigned int fe_index = *pointer;
+
+ Assert (fe_index < dof_handler.get_fe().size(),
+ ExcInternalError());
+
+ if (counter == n)
+ return fe_index;
+
+ ++counter;
+ pointer += dof_handler.get_fe()[fe_index]
+ .template n_dofs_per_object<structdim>() + 1;
+ }
+ }
+
+
+
+ template <int structdim>
+ template <int dim, int spacedim>
+ inline
+ bool
+ DoFIndicesOnFacesOrEdges<structdim>::
+ fe_index_is_active (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
+ const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int obj_level) const
+ {
+ Assert (&dof_handler != 0,
+ ExcMessage ("No DoFHandler is specified for this iterator"));
+ Assert (&dof_handler.get_fe() != 0,
+ ExcMessage ("No finite element collection is associated with "
+ "this DoFHandler"));
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, static_cast<unsigned int>(dof_offsets.size())));
+ Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
+ ExcMessage ("You need to specify a FE index when working "
+ "with hp DoFHandlers"));
+ Assert (fe_index < dof_handler.get_fe().size(),
+ ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const types::global_dof_index starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ while (true)
+ {
+ if (*pointer == numbers::invalid_dof_index)
+ // end of list reached
+ return false;
+ else if (*pointer == fe_index)
+ return true;
+ else
+ pointer += static_cast<types::global_dof_index>(
+ dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>()+1);
+ }
+ }
+
+
} // namespace hp
} // namespace internal
{
namespace hp
{
- template <int> class DoFLevel;
- template <int> class DoFFaces;
- template <int> class DoFObjects;
+ class DoFLevel;
+ template <int> class DoFIndicesOnFaces;
+ template <int> class DoFIndicesOnFacesOrEdges;
namespace DoFHandler
{
* <tt>levels[]</tt> tree of
* the Triangulation objects.
*/
- std::vector<dealii::internal::hp::DoFLevel<dim>*> levels;
+ std::vector<dealii::internal::hp::DoFLevel*> levels;
/**
* Space to store the DoF
* <tt>faces</tt> pointer of
* the Triangulation objects.
*/
- dealii::internal::hp::DoFFaces<dim> *faces;
+ dealii::internal::hp::DoFIndicesOnFaces<dim> *faces;
/**
* A structure that contains all
* the functions that set and
* retrieve vertex dof indices.
*/
- template <int> friend class dealii::internal::hp::DoFLevel;
- template <int> friend class dealii::internal::hp::DoFObjects;
+ template <int> friend class dealii::internal::hp::DoFIndicesOnFacesOrEdges;
friend struct dealii::internal::hp::DoFHandler::Implementation;
};
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2003 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__hp_dof_level_h
+#define __deal2__hp_dof_level_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+
+#include <vector>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace hp
+{
+ template <int, int> class DoFHandler;
+ template <int, int> class FECollection;
+}
+
+
+namespace internal
+{
+ namespace hp
+ {
+ namespace DoFHandler
+ {
+ struct Implementation;
+ }
+ }
+ namespace DoFCellAccessor
+ {
+ struct Implementation;
+ }
+}
+
+
+namespace internal
+{
+ namespace hp
+ {
+ /**
+ * This is the class that stores the degrees of freedom on cells in a hp
+ * hierarchy. Compared to faces and edges, the task here is simple since
+ * each cell can only have a single active finite element index. Consequently,
+ * all we need is one long array with DoF indices and one array of offsets
+ * where each cell's indices start within the array of indices. This is in
+ * contrast to the DoFObjects class where each face or edge may have more than
+ * one associated finite element with corresponding degrees of freedom.
+ *
+ * The data stored here is represented by three arrays
+ * - The @p active_fe_indices array stores for each cell which
+ * finite element is used on this cell. Since some cells are not active
+ * on the current level, some entries in this array may represent an
+ * invalid value.
+ * - The @p dof_indices array stores for each active cell on the current
+ * level the dofs that associated with the <i>interior</i> of the cell,
+ * i.e., the @p dofs_per_line dofs associated with the line in 1d, and
+ * @p dofs_per_quad and @p dofs_per_hex in 2d and 3d. These numbers are
+ * in general smaller than @p dofs_per_cell.
+ * - The @p dof_offsets array stores, for each cell, the starting point
+ * of the dof indices corresponding to this cell in the @p dof_indices
+ * array. This is analogous to how we store data in compressed row storage
+ * for sparse matrices. For cells that are not active on the current level,
+ * we store an invalid value for the starting index.
+ *
+ * <h3>Compression</h3>
+ *
+ * It is common for the indices stored in @p dof_indices for one cell to be
+ * numbered consecutively. For example, using the standard numbering (without
+ * renumbering DoFs), the quad dofs on the first cell of a mesh when using a
+ * $Q_3$ element will be numbered <tt>12, 13, 14, 15</tt>. This allows for
+ * compression if we only store the first entry and have some way to mark
+ * the DoFs on this object as compressed. Here, compression means that we
+ * know that subsequent DoF indices can be obtained from the previous ones
+ * by just incrementing them by one -- in other words, we use a variant of doing
+ * run-length encoding. The way to do this is that we
+ * use positive FE indices for uncompressed sets of DoFs and if a set of indices
+ * is compressed, then we instead store the FE index in binary complement (which
+ * we can identify by looking at the sign bit when interpreting the number as a
+ * signed one). There are two functions, compress_data() and uncompress_data()
+ * that convert between the two possible representations.
+ *
+ * Note that compression is not always possible. For example, if one renumbered
+ * the example above using DoFRenumbering::downstream with $(1,0)^T$ as
+ * direction, then they would likely be numbered <tt>12, 14, 13, 15</tt>, which
+ * can not be compressed using run-length encoding.
+ */
+ class DoFLevel
+ {
+ private:
+ /**
+ * The type in which we store the offsets into the dof_indices array.
+ */
+ typedef unsigned int offset_type;
+
+ /**
+ * The type in which we store the active FE index.
+ */
+ typedef unsigned short int active_fe_index_type;
+
+ /**
+ * A signed type that matches the type in which we store the active FE
+ * index. We use this in computing binary complements.
+ */
+ typedef signed short int signed_active_fe_index_type;
+
+ /**
+ * Indices specifying the finite element of hp::FECollection to
+ * use for the different cells on the current level. The vector
+ * stores one element per cell since the active_fe_index is
+ * unique for cells.
+ *
+ * If a cell is not active on the level corresponding to the
+ * current object (i.e., it has children on higher levels) then
+ * it does not have an associated fe index and we store
+ * an invalid fe index marker instead.
+ */
+ std::vector<active_fe_index_type> active_fe_indices;
+
+ /**
+ * Store the start index for the degrees of freedom of each
+ * object in the @p dof_indices array. If the cell corresponding to
+ * a particular index in this array is not active on this level,
+ * then we do not store any DoFs for it. In that case, the offset
+ * we store here must be an invalid number and indeed we store
+ * <code>(std::vector<types::global_dof_index>::size_type)(-1)</code>
+ * for it.
+ *
+ * The type we store is then obviously the type the @p dof_indices array
+ * uses for indexing.
+ */
+ std::vector<offset_type> dof_offsets;
+
+ /**
+ * Store the global indices of the degrees of freedom.
+ * information. The dof_offsets field determines where each
+ * (active) cell's data is stored.
+ */
+ std::vector<types::global_dof_index> dof_indices;
+
+ /**
+ * The offsets for each cell of the cache that holds all DoF indices.
+ */
+ std::vector<offset_type> cell_cache_offsets;
+
+ /**
+ * Cache for the DoF indices
+ * on cells. The size of this
+ * array equals the sum over all cells of selected_fe[active_fe_index[cell]].dofs_per_cell.
+ */
+ std::vector<types::global_dof_index> cell_dof_indices_cache;
+
+ public:
+
+ /**
+ * Set the global index of
+ * the @p local_index-th
+ * degree of freedom located
+ * on the object with number @p
+ * obj_index to the value
+ * given by @p global_index. The @p
+ * dof_handler argument is
+ * used to access the finite
+ * element that is to be used
+ * to compute the location
+ * where this data is stored.
+ *
+ * The third argument, @p
+ * fe_index, denotes which of
+ * the finite elements
+ * associated with this
+ * object we shall
+ * access. Refer to the
+ * general documentation of
+ * the internal::hp::DoFLevel
+ * class template for more
+ * information.
+ */
+ void
+ set_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const types::global_dof_index global_index);
+
+ /**
+ * Return the global index of
+ * the @p local_index-th
+ * degree of freedom located
+ * on the object with number @p
+ * obj_index. The @p
+ * dof_handler argument is
+ * used to access the finite
+ * element that is to be used
+ * to compute the location
+ * where this data is stored.
+ *
+ * The third argument, @p
+ * fe_index, denotes which of
+ * the finite elements
+ * associated with this
+ * object we shall
+ * access. Refer to the
+ * general documentation of
+ * the internal::hp::DoFLevel
+ * class template for more
+ * information.
+ */
+ types::global_dof_index
+ get_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index) const;
+
+ /**
+ * Return the fe_index of the
+ * active finite element
+ * on this object.
+ */
+ unsigned int
+ active_fe_index (const unsigned int obj_index) const;
+
+ /**
+ * Check whether a given
+ * finite element index is
+ * used on the present
+ * object or not.
+ */
+ bool
+ fe_index_is_active (const unsigned int obj_index,
+ const unsigned int fe_index) const;
+
+ /**
+ * Set the fe_index of the
+ * active finite element
+ * on this object.
+ */
+ void
+ set_active_fe_index (const unsigned int obj_index,
+ const unsigned int fe_index);
+
+ /**
+ * Return a pointer to the beginning of the DoF indices cache
+ * for a given cell.
+ *
+ * @param obj_index The number of the cell we are looking at.
+ * @param dofs_per_cell The number of DoFs per cell for this cell. This
+ * is not used for the hp case but necessary to keep the interface
+ * the same as for the non-hp case.
+ * @return A pointer to the first DoF index for the current cell. The
+ * next dofs_per_cell indices are for the current cell.
+ */
+ const types::global_dof_index *
+ get_cell_cache_start (const unsigned int obj_index,
+ const unsigned int dofs_per_cell) const;
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object.
+ */
+ std::size_t memory_consumption () const;
+
+ private:
+ /**
+ * Compress the arrays that store dof indices by using a variant
+ * of run-length encoding. See the general documentation of this
+ * class for more information.
+ *
+ * @param fe_collection The object that can tell us how many
+ * degrees of freedom each of the finite elements has that we
+ * store in this object.
+ */
+ template <int dim, int spacedim>
+ void compress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection);
+
+ /**
+ * Uncompress the arrays that store dof indices by using a variant
+ * of run-length encoding. See the general documentation of this
+ * class for more information.
+ *
+ * @param fe_collection The object that can tell us how many
+ * degrees of freedom each of the finite elements has that we
+ * store in this object.
+ */
+ template <int dim, int spacedim>
+ void uncompress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection);
+
+ /**
+ * Make hp::DoFHandler and its auxiliary class a friend since it
+ * is the class that needs to create these data structures.
+ */
+ template <int, int> friend class dealii::hp::DoFHandler;
+ friend struct dealii::internal::hp::DoFHandler::Implementation;
+ friend struct dealii::internal::DoFCellAccessor::Implementation;
+ };
+
+
+ // -------------------- template functions --------------------------------
+
+ inline
+ types::global_dof_index
+ DoFLevel::
+ get_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index) const
+ {
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (fe_index == active_fe_indices[obj_index],
+ ExcMessage ("FE index does not match that of the present cell"));
+
+ // see if the dof_indices array has been compressed for this
+ // particular cell
+ if ((signed_active_fe_index_type)active_fe_indices[obj_index]>=0)
+ return dof_indices[dof_offsets[obj_index]+local_index];
+ else
+ return dof_indices[dof_offsets[obj_index]]+local_index;
+ }
+
+
+
+ inline
+ void
+ DoFLevel::
+ set_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const types::global_dof_index global_index)
+ {
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+ Assert ((signed_active_fe_index_type)active_fe_indices[obj_index]>=0,
+ ExcMessage ("This function can no longer be called after compressing the dof_indices array"));
+ Assert (fe_index == active_fe_indices[obj_index],
+ ExcMessage ("FE index does not match that of the present cell"));
+ dof_indices[dof_offsets[obj_index]+local_index] = global_index;
+ }
+
+
+
+ inline
+ unsigned int
+ DoFLevel::
+ active_fe_index (const unsigned int obj_index) const
+ {
+ Assert (obj_index < active_fe_indices.size(),
+ ExcIndexRange (obj_index, 0, active_fe_indices.size()));
+
+ if (((signed_active_fe_index_type)active_fe_indices[obj_index]) >= 0)
+ return active_fe_indices[obj_index];
+ else
+ return (active_fe_index_type)~(signed_active_fe_index_type)active_fe_indices[obj_index];
+ }
+
+
+
+ inline
+ bool
+ DoFLevel::
+ fe_index_is_active (const unsigned int obj_index,
+ const unsigned int fe_index) const
+ {
+ return (fe_index == active_fe_index(obj_index));
+ }
+
+
+ inline
+ void
+ DoFLevel::
+ set_active_fe_index (const unsigned int obj_index,
+ const unsigned int fe_index)
+ {
+ Assert (obj_index < active_fe_indices.size(),
+ ExcIndexRange (obj_index, 0, active_fe_indices.size()));
+
+ active_fe_indices[obj_index] = fe_index;
+ }
+
+
+
+ inline
+ const types::global_dof_index *
+ DoFLevel::get_cell_cache_start (const unsigned int obj_index,
+ const unsigned int dofs_per_cell) const
+ {
+ Assert (obj_index < cell_cache_offsets.size(),
+ ExcInternalError());
+ Assert (cell_cache_offsets[obj_index]+dofs_per_cell
+ <=
+ cell_dof_indices_cache.size(),
+ ExcInternalError());
+
+ return &cell_dof_indices_cache[cell_cache_offsets[obj_index]];
+ }
+ } // namespace hp
+
+} // namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-// $Id$
-//
-// Copyright (C) 2006 - 2013 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-#ifndef __deal2__hp_dof_objects_h
-#define __deal2__hp_dof_objects_h
-
-#include <deal.II/base/config.h>
-#include <deal.II/hp/fe_collection.h>
-#include <deal.II/hp/dof_handler.h>
-
-#include <vector>
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace internal
-{
- namespace hp
- {
-
- /**
- * Store the indices of the degrees of freedom which are located on
- * objects of dimension @p structdim < dim, i.e., for faces or edges
- * of cells. This is opposed to the internal::hp::DoFLevels class
- * that stores the DoF indices on cells.
- *
- * The things we store here is very similar to what is stored in the
- * internal::DoFHandler::DoFObjects classes (see there for more
- * information, in particular on the layout of the class hierarchy,
- * and the use of file names).
- *
- * <h4>Offset computations</h4>
- *
- * For hp methods, not all cells may use the same finite element, and
- * it is consequently more complicated to determine where the DoF
- * indices for a given line, quad, or hex are stored. As described in
- * the documentation of the internal::DoFHandler::DoFLevel class, we
- * can compute the location of the first line DoF, for example, by
- * calculating the offset as <code>line_index *
- * dof_handler.get_fe().dofs_per_line</code>. This of course doesn't
- * work any more if different lines may have different numbers of
- * degrees of freedom associated with them. Consequently, rather than
- * using this simple multiplication, the dofs array has an associated
- * array dof_offsets. The data corresponding to a
- * line then starts at index <code>line_dof_offsets[line_index]</code>
- * within the <code>line_dofs</code> array.
- *
- *
- * <h4>Multiple data sets per object</h4>
- *
- * If two adjacent cells use different finite elements, then
- * the face that they share needs to store DoF indices for both
- * involved finite elements. While faces therefore have to have at
- * most two sets of DoF indices, it is easy to see that edges and
- * vertices can have as many sets of DoF indices associated with them
- * as there are adjacent cells.
- *
- * Consequently, for objects that have a lower dimensionality than
- * cells, we have to store a map from the finite element index to the
- * set of DoF indices associated. Since real sets are typically very
- * inefficient to store, and since most of the time we expect the
- * number of individual keys to be small (frequently, adjacent cells
- * will have the same finite element, and only a single entry will
- * exist in the map), what we do is instead to store a linked list. In
- * this format, the first entry starting at position
- * <code>lines.dofs[lines.dof_offsets[line_index]]</code> will denote
- * the finite element index of the set of DoF indices following; after
- * this set, we will store the finite element index of the second set
- * followed by the corresponding DoF indices; and so on. Finally, when
- * all finite element indices adjacent to this object have been
- * covered, we write a -1 to indicate the end of the list.
- *
- * Access to this kind of data, as well as the distinction between
- * cells and objects of lower dimensionality are encoded in the
- * accessor functions, DoFObjects::set_dof_index() and
- * DoFLevel::get_dof_index(). They are able to traverse this
- * list and pick out or set a DoF index given the finite element index
- * and its location within the set of DoFs corresponding to this
- * finite element.
- *
- *
- * @ingroup hp
- * @author Tobias Leicht, 2006
- */
- template <int structdim>
- class DoFObjects
- {
- public:
- /**
- * Store the start index for
- * the degrees of freedom of each
- * object in the @p dofs array.
- *
- * The type we store is then obviously the type the @p dofs array
- * uses for indexing.
- */
- std::vector<std::vector<types::global_dof_index>::size_type> dof_offsets;
-
- /**
- * Store the global indices of
- * the degrees of freedom. See
- * DoFLevel() for detailed
- * information.
- */
- std::vector<types::global_dof_index> dofs;
-
- /**
- * Set the global index of
- * the @p local_index-th
- * degree of freedom located
- * on the object with number @p
- * obj_index to the value
- * given by @p global_index. The @p
- * dof_handler argument is
- * used to access the finite
- * element that is to be used
- * to compute the location
- * where this data is stored.
- *
- * The third argument, @p
- * fe_index, denotes which of
- * the finite elements
- * associated with this
- * object we shall
- * access. Refer to the
- * general documentation of
- * the internal::hp::DoFLevel
- * class template for more
- * information.
- */
- template <int dim, int spacedim>
- void
- set_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const unsigned int local_index,
- const types::global_dof_index global_index,
- const unsigned int obj_level);
-
- /**
- * Return the global index of
- * the @p local_index-th
- * degree of freedom located
- * on the object with number @p
- * obj_index. The @p
- * dof_handler argument is
- * used to access the finite
- * element that is to be used
- * to compute the location
- * where this data is stored.
- *
- * The third argument, @p
- * fe_index, denotes which of
- * the finite elements
- * associated with this
- * object we shall
- * access. Refer to the
- * general documentation of
- * the internal::hp::DoFLevel
- * class template for more
- * information.
- */
- template <int dim, int spacedim>
- types::global_dof_index
- get_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const unsigned int local_index,
- const unsigned int obj_level) const;
-
- /**
- * Return the number of
- * finite elements that are
- * active on a given
- * object. If this is a cell,
- * the answer is of course
- * one. If it is a face, the
- * answer may be one or two,
- * depending on whether the
- * two adjacent cells use the
- * same finite element or
- * not. If it is an edge in
- * 3d, the possible return
- * value may be one or any
- * other value larger than
- * that.
- *
- * If the object is not part
- * of an active cell, then no
- * degrees of freedom have
- * been distributed and zero
- * is returned.
- */
- template <int dim, int spacedim>
- unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index) const;
-
- /**
- * Return the fe_index of the
- * n-th active finite element
- * on this object.
- */
- template <int dim, int spacedim>
- types::global_dof_index
- nth_active_fe_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n) const;
-
- /**
- * Check whether a given
- * finite element index is
- * used on the present
- * object or not.
- */
- template <int dim, int spacedim>
- bool
- fe_index_is_active (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const unsigned int obj_level) const;
-
- /**
- * Determine an estimate for the
- * memory consumption (in bytes)
- * of this object.
- */
- std::size_t memory_consumption () const;
- };
-
-
- // --------------------- inline and template functions ------------------
-
- template <int structdim>
- template <int dim, int spacedim>
- inline
- types::global_dof_index
- DoFObjects<structdim>::
- get_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const unsigned int local_index,
- const unsigned int obj_level) const
- {
- Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
- ExcMessage ("You need to specify a FE index when working "
- "with hp DoFHandlers"));
- Assert (&dof_handler != 0,
- ExcMessage ("No DoFHandler is specified for this iterator"));
- Assert (&dof_handler.get_fe() != 0,
- ExcMessage ("No finite element collection is associated with "
- "this DoFHandler"));
- Assert (fe_index < dof_handler.get_fe().size(),
- ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
- Assert (local_index <
- dof_handler.get_fe()[fe_index].template n_dofs_per_object<structdim>(),
- ExcIndexRange(local_index, 0,
- dof_handler.get_fe()[fe_index]
- .template n_dofs_per_object<structdim>()));
- Assert (obj_index < dof_offsets.size(),
- ExcIndexRange (obj_index, 0, dof_offsets.size()));
-
- // make sure we are on an
- // object for which DoFs have
- // been allocated at all
- Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
- ExcMessage ("You are trying to access degree of freedom "
- "information for an object on which no such "
- "information is available"));
-
- Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
-
- // there may be multiple finite elements associated with
- // this object. hop along the list of index sets until we
- // find the one with the correct fe_index, and then poke
- // into that part. trigger an exception if we can't find a
- // set for this particular fe_index
- const types::global_dof_index starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- while (true)
- {
- Assert (*pointer != numbers::invalid_dof_index,
- ExcInternalError());
- if (*pointer == fe_index)
- return *(pointer + 1 + local_index);
- else
- pointer += static_cast<types::global_dof_index>(
- dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<structdim>() + 1);
- }
- }
-
-
-
- template <int structdim>
- template <int dim, int spacedim>
- inline
- void
- DoFObjects<structdim>::
- set_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const unsigned int local_index,
- const types::global_dof_index global_index,
- const unsigned int obj_level)
- {
- Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
- ExcMessage ("You need to specify a FE index when working "
- "with hp DoFHandlers"));
- Assert (&dof_handler != 0,
- ExcMessage ("No DoFHandler is specified for this iterator"));
- Assert (&dof_handler.get_fe() != 0,
- ExcMessage ("No finite element collection is associated with "
- "this DoFHandler"));
- Assert (fe_index < dof_handler.get_fe().size(),
- ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
- Assert (local_index <
- dof_handler.get_fe()[fe_index].template n_dofs_per_object<structdim>(),
- ExcIndexRange(local_index, 0,
- dof_handler.get_fe()[fe_index]
- .template n_dofs_per_object<structdim>()));
- Assert (obj_index < dof_offsets.size(),
- ExcIndexRange (obj_index, 0, dof_offsets.size()));
-
- // make sure we are on an
- // object for which DoFs have
- // been allocated at all
- Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
- ExcMessage ("You are trying to access degree of freedom "
- "information for an object on which no such "
- "information is available"));
-
- Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
-
- // there may be multiple finite elements associated with
- // this object. hop along the list of index sets until we
- // find the one with the correct fe_index, and then poke
- // into that part. trigger an exception if we can't find a
- // set for this particular fe_index
- const types::global_dof_index starting_offset = dof_offsets[obj_index];
- types::global_dof_index *pointer = &dofs[starting_offset];
- while (true)
- {
- Assert (*pointer != numbers::invalid_dof_index,
- ExcInternalError());
- if (*pointer == fe_index)
- {
- *(pointer + 1 + local_index) = global_index;
- return;
- }
- else
- pointer += dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<structdim>() + 1;
- }
- }
-
-
-
- template <int structdim>
- template <int dim, int spacedim>
- inline
- unsigned int
- DoFObjects<structdim>::
- n_active_fe_indices (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index) const
- {
- Assert (&dof_handler != 0,
- ExcMessage ("No DoFHandler is specified for this iterator"));
- Assert (&dof_handler.get_fe() != 0,
- ExcMessage ("No finite element collection is associated with "
- "this DoFHandler"));
- Assert (obj_index < dof_offsets.size(),
- ExcIndexRange (obj_index, 0, dof_offsets.size()));
-
- // make sure we are on an
- // object for which DoFs have
- // been allocated at all
- if (dof_offsets[obj_index] == numbers::invalid_dof_index)
- return 0;
-
- Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
-
- // there may be multiple finite elements associated with this
- // object. hop along the list of index sets until we find the
- // one with the correct fe_index, and then poke into that
- // part. trigger an exception if we can't find a set for this
- // particular fe_index
- const unsigned int starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- unsigned int counter = 0;
- while (true)
- {
- if (*pointer == numbers::invalid_dof_index)
- // end of list reached
- return counter;
- else
- {
- ++counter;
- pointer += dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<structdim>() + 1;
- }
- }
- }
-
-
-
- template <int structdim>
- template <int dim, int spacedim>
- inline
- types::global_dof_index
- DoFObjects<structdim>::
- nth_active_fe_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_level,
- const unsigned int obj_index,
- const unsigned int n) const
- {
- Assert (&dof_handler != 0,
- ExcMessage ("No DoFHandler is specified for this iterator"));
- Assert (&dof_handler.get_fe() != 0,
- ExcMessage ("No finite element collection is associated with "
- "this DoFHandler"));
- Assert (obj_index < dof_offsets.size(),
- ExcIndexRange (obj_index, 0, dof_offsets.size()));
-
- // make sure we are on an
- // object for which DoFs have
- // been allocated at all
- Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
- ExcMessage ("You are trying to access degree of freedom "
- "information for an object on which no such "
- "information is available"));
-
- Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
-
- Assert (n < n_active_fe_indices(dof_handler, obj_index),
- ExcIndexRange (n, 0,
- n_active_fe_indices(dof_handler, obj_index)));
-
- // there may be multiple finite elements associated with
- // this object. hop along the list of index sets until we
- // find the one with the correct fe_index, and then poke
- // into that part. trigger an exception if we can't find a
- // set for this particular fe_index
- const unsigned int starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- unsigned int counter = 0;
- while (true)
- {
- Assert (*pointer != numbers::invalid_dof_index,
- ExcInternalError());
-
- const unsigned int fe_index = *pointer;
-
- Assert (fe_index < dof_handler.get_fe().size(),
- ExcInternalError());
-
- if (counter == n)
- return fe_index;
-
- ++counter;
- pointer += dof_handler.get_fe()[fe_index]
- .template n_dofs_per_object<structdim>() + 1;
- }
- }
-
-
-
- template <int structdim>
- template <int dim, int spacedim>
- inline
- bool
- DoFObjects<structdim>::
- fe_index_is_active (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
- const unsigned int obj_index,
- const unsigned int fe_index,
- const unsigned int obj_level) const
- {
- Assert (&dof_handler != 0,
- ExcMessage ("No DoFHandler is specified for this iterator"));
- Assert (&dof_handler.get_fe() != 0,
- ExcMessage ("No finite element collection is associated with "
- "this DoFHandler"));
- Assert (obj_index < dof_offsets.size(),
- ExcIndexRange (obj_index, 0, static_cast<unsigned int>(dof_offsets.size())));
- Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
- ExcMessage ("You need to specify a FE index when working "
- "with hp DoFHandlers"));
- Assert (fe_index < dof_handler.get_fe().size(),
- ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
-
- // make sure we are on an
- // object for which DoFs have
- // been allocated at all
- Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
- ExcMessage ("You are trying to access degree of freedom "
- "information for an object on which no such "
- "information is available"));
-
- Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
-
- // there may be multiple finite elements associated with
- // this object. hop along the list of index sets until we
- // find the one with the correct fe_index, and then poke
- // into that part. trigger an exception if we can't find a
- // set for this particular fe_index
- const types::global_dof_index starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- while (true)
- {
- if (*pointer == numbers::invalid_dof_index)
- // end of list reached
- return false;
- else if (*pointer == fe_index)
- return true;
- else
- pointer += static_cast<types::global_dof_index>(
- dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<structdim>()+1);
- }
- }
-
- }
-}
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
*/
BlockVector (const std::vector<size_type> &block_sizes);
+ /**
+ * Construct a block vector with an IndexSet for the local range
+ * and ghost entries for each block.
+ */
+ BlockVector (const std::vector<IndexSet> &local_ranges,
+ const std::vector<IndexSet> &ghost_indices,
+ const MPI_Comm communicator);
+
+ /**
+ * Same as above but the ghost indicies are assumed to be empty.
+ */
+ BlockVector (const std::vector<IndexSet> &local_ranges,
+ const MPI_Comm communicator);
+
/**
* Destructor. Clears memory.
*/
}
+ template <typename Number>
+ inline
+ BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
+ const std::vector<IndexSet> &ghost_indices,
+ const MPI_Comm communicator)
+ {
+ std::vector<size_type> sizes(local_ranges.size());
+ for (unsigned int i=0; i<local_ranges.size(); ++i)
+ sizes[i] = local_ranges[i].size();
+
+ this->block_indices.reinit(sizes);
+ this->components.resize(this->n_blocks());
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ this->block(i).reinit(local_ranges[i], ghost_indices[i], communicator);
+ }
+
+
+ template <typename Number>
+ inline
+ BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
+ const MPI_Comm communicator)
+ {
+ std::vector<size_type> sizes(local_ranges.size());
+ for (unsigned int i=0; i<local_ranges.size(); ++i)
+ sizes[i] = local_ranges[i].size();
+
+ this->block_indices.reinit(sizes);
+ this->components.resize(this->n_blocks());
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ this->block(i).reinit(local_ranges[i], communicator);
+ }
+
+
template <typename Number>
inline
BlockVector<Number> &
BlockVector<Number>::operator = (const BlockVector &v)
{
- reinit (v, true);
- BaseClass::operator = (v);
+ // we only allow assignment to vectors with the same number of blocks
+ // or to an empty BlockVector
+ Assert (this->n_blocks() == 0 || this->n_blocks() == v.n_blocks(),
+ ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
+
+ if (this->n_blocks() != v.n_blocks())
+ reinit(v.n_blocks(), true);
+
+ for (size_type i=0; i<this->n_blocks(); ++i)
+ this->components[i] = v.block(i);
+
+ this->collect_sizes();
return *this;
}
const IndexSet &ghost_indices,
const MPI_Comm communicator);
+ /**
+ * Same constructor as above but without any ghost indices.
+ */
+ Vector (const IndexSet &local_range,
+ const MPI_Comm communicator);
+
/**
* Create the vector based on the parallel partitioning described in @p
* partitioner. The input argument is a shared pointer, which store the
const IndexSet &ghost_indices,
const MPI_Comm communicator);
+ /**
+ * Same as above, but without ghost entries.
+ */
+ void reinit (const IndexSet &local_range,
+ const MPI_Comm communicator);
+
/**
* Initialize the vector given to the parallel partitioning described in
* @p partitioner. The input argument is a shared pointer, which store
+ template <typename Number>
+ inline
+ Vector<Number>::Vector (const IndexSet &local_range,
+ const MPI_Comm communicator)
+ :
+ allocated_size (0),
+ val (0),
+ import_data (0),
+ vector_view (0, static_cast<Number *>(0))
+ {
+ IndexSet ghost_indices(local_range.size());
+ reinit (local_range, ghost_indices, communicator);
+ }
+
+
+
template <typename Number>
inline
Vector<Number>::Vector (const size_type size)
reinit (c, true);
}
vector_view = c.vector_view;
+ update_ghost_values();
return *this;
}
const IndexSet &ghost_indices,
const MPI_Comm communicator)
{
- // set up parallel partitioner with index sets
- // and communicator
+ // set up parallel partitioner with index sets and communicator
+ std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> new_partitioner
+ (new Utilities::MPI::Partitioner (locally_owned_indices,
+ ghost_indices, communicator));
+ reinit (new_partitioner);
+ }
+
+
+
+ template <typename Number>
+ void
+ Vector<Number>::reinit (const IndexSet &locally_owned_indices,
+ const MPI_Comm communicator)
+ {
+ // set up parallel partitioner with index sets and communicator
+ IndexSet ghost_indices(locally_owned_indices.size());
std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> new_partitioner
(new Utilities::MPI::Partitioner (locally_owned_indices,
ghost_indices, communicator));
{
namespace Vector
{
- typedef types::global_dof_index size_type;
-
- template<typename T>
+ template <typename T>
void set_subrange (const T s,
- const size_type begin,
- const size_type end,
+ const typename dealii::Vector<T>::size_type begin,
+ const typename dealii::Vector<T>::size_type end,
dealii::Vector<T> &dst)
{
if (s == T())
std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s);
}
- template<typename T>
- void copy_subrange (const size_type begin,
- const size_type end,
+
+ template <typename T>
+ void copy_subrange (const typename dealii::Vector<T>::size_type begin,
+ const typename dealii::Vector<T>::size_type end,
const dealii::Vector<T> &src,
dealii::Vector<T> &dst)
{
(end-begin)*sizeof(T));
}
- template<typename T, typename U>
- void copy_subrange (const size_type begin,
- const size_type end,
+
+ template <typename T, typename U>
+ void copy_subrange (const typename dealii::Vector<T>::size_type begin,
+ const typename dealii::Vector<T>::size_type end,
const dealii::Vector<T> &src,
dealii::Vector<U> &dst)
{
*p = *q;
}
+
template<typename T, typename U>
- void copy_subrange_wrap (const size_type begin,
- const size_type end,
+ void copy_subrange_wrap (const typename dealii::Vector<T>::size_type begin,
+ const typename dealii::Vector<T>::size_type end,
const dealii::Vector<T> &src,
dealii::Vector<U> &dst)
{
copy_subrange (begin, end, src, dst);
}
+
template <typename T, typename U>
void copy_vector (const dealii::Vector<T> &src,
dealii::Vector<U> &dst)
if (PointerComparison::equal(&src, &dst))
return;
- const size_type vec_size = src.size();
- const size_type dst_size = dst.size();
+ const typename dealii::Vector<T>::size_type vec_size = src.size();
+ const typename dealii::Vector<U>::size_type dst_size = dst.size();
if (dst_size != vec_size)
dst.reinit (vec_size, true);
if (vec_size>internal::Vector::minimum_parallel_grain_size)
// The code returns the result as the last argument in order to make
// spawning tasks simpler and use automatic template deduction.
template <typename Operation, typename Number, typename Number2,
- typename ResultType>
+ typename ResultType, typename size_type>
void accumulate (const Operation &op,
const Number *X,
const Number2 *Y,
ResultType r0, r1, r2, r3;
Threads::TaskGroup<> task_group;
task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
- ResultType>,
+ ResultType,size_type>,
op, X, Y, power, new_size, r0);
task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
- ResultType>,
+ ResultType,size_type>,
op, X+new_size, Y+new_size, power,
new_size, r1);
task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
- ResultType>,
+ ResultType,size_type>,
op, X+2*new_size, Y+2*new_size, power,
new_size, r2);
task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
- ResultType>,
+ ResultType,size_type>,
op, X+3*new_size, Y+3*new_size, power,
vec_size-3*new_size, r3);
task_group.join_all();
* using a Q1-mapping for cell
* boundaries to call the other
* point_difference() function.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
void point_difference (const DoFHandler<dim,spacedim> &dof,
* this function uses an
* arbitrary mapping to evaluate
* the difference.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
void point_difference (const Mapping<dim, spacedim> &mapping,
* using a Q1-mapping for cell
* boundaries to call the other
* point_difference() function.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
void
/**
* Same as above for hp.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
void
* "Possibilities for extensions" part of
* the results section of @ref step_3
* "step-3".
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
double
/**
* Same as above for hp.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
double
* function of the same name,
* this function uses an arbitrary
* mapping to evaluate the difference.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
void
/**
* Same as above for hp.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
void
* function of the same name,
* this function uses an arbitrary
* mapping to evaluate the difference.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
double
/**
* Same as above for hp.
+ *
+ * @note If the cell in which the point is found
+ * is not locally owned, an exception of type
+ * VectorTools<dim, InVector, spacedim>::ExcPointNotAvailableHere
+ * is thrown.
*/
template <int dim, class InVector, int spacedim>
double
* Exception
*/
DeclException0 (ExcNonInterpolatingFE);
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcPointNotAvailableHere);
}
std::vector<double> &dof_values,
std::vector<bool> &dofs_processed)
{
- const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->line (line)->diameter ();
+ const double tol = 0.5 * cell->face (face)->line (line)->diameter () / cell->get_fe ().degree;
const unsigned int dim = 3;
const unsigned int spacedim = 3;
{
case 2:
{
+ const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree;
std::vector<Point<dim> >
tangentials (fe_values.n_quadrature_points);
= reference_quadrature_points[q_point];
shifted_reference_point_1 (face_coordinate_direction[face])
- += 1e-13;
+ += tol;
shifted_reference_point_2 (face_coordinate_direction[face])
- -= 1e-13;
+ -= tol;
tangentials[q_point]
- = 2e13
- * (fe_values.get_mapping ()
+ = (fe_values.get_mapping ()
.transform_unit_to_real_cell (cell,
shifted_reference_point_1)
-
fe_values.get_mapping ()
.transform_unit_to_real_cell (cell,
- shifted_reference_point_2));
+ shifted_reference_point_2))
+ / tol;
tangentials[q_point]
/= std::sqrt (tangentials[q_point].square ());
// Compute the degrees
// if the degree of
// freedom is not
// already constrained.
- const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->diameter ();
-
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
&& !(constraints.is_constrained (face_dof_indices[dof])))
{
constraints.add_line (face_dof_indices[dof]);
- if (std::abs (dof_values[dof]) > tol)
+ if (std::abs (dof_values[dof]) > 1e-13)
constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
}
}
// Store the computed
// values in the global
// vector.
-
cell->face (face)->get_dof_indices (face_dof_indices,
cell->active_fe_index ());
- const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter ();
-
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
&& !(constraints.is_constrained (face_dof_indices[dof])))
{
constraints.add_line (face_dof_indices[dof]);
- if (std::abs (dof_values[dof]) > tol)
+ if (std::abs (dof_values[dof]) > 1e-13)
constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
}
}
cell->face (face)->get_dof_indices (face_dof_indices,
cell->active_fe_index ());
- const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->diameter ();
-
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
&& !(constraints.is_constrained (face_dof_indices[dof])))
{
constraints.add_line (face_dof_indices[dof]);
- if (std::abs (dof_values[dof]) > tol)
+ if (std::abs (dof_values[dof]) > 1e-13)
constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
}
}
cell->face (face)->get_dof_indices (face_dof_indices,
cell->active_fe_index ());
- const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter ();
-
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
&& !(constraints.is_constrained (face_dof_indices[dof])))
{
constraints.add_line (face_dof_indices[dof]);
- if (std::abs (dof_values[dof]) > tol)
+ if (std::abs (dof_values[dof]) > 1e-13)
constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
}
}
hp::FEFaceValues<dim,spacedim> x_fe_face_values (mapping_collection,
fe_collection,
face_quadrature_collection,
- update_q_points);
+ update_q_points |
+ update_normal_vectors);
// have a map that stores normal vectors for each vector-dof tuple we want
// to constrain. since we can get at the same vector dof tuple more than
// to check that they are small *relative to something
// else*). we do this and then normalize the length of the
// vector back to one, just to be on the safe side
+ //
+ // one more point: we would like to use the "real" normal
+ // vector here, as provided by the boundary description
+ // and as opposed to what we get from the FEValues object.
+ // we do this in the immediately next line, but as is
+ // obvious, the boundary only has a vague idea which side
+ // of a cell it is on -- indicated by the face number. in
+ // other words, it may provide the inner or outer normal.
+ // by and large, there is no harm from this, since the
+ // tangential vector we compute is still the same. however,
+ // we do average over normal vectors from adjacent cells
+ // and if they have recorded normal vectors from the inside
+ // once and from the outside the other time, then this
+ // averaging is going to run into trouble. as a consequence
+ // we ask the mapping after all for its normal vector,
+ // but we only ask it so that we can possibly correct the
+ // sign of the normal vector provided by the boundary
+ // if they should point in different directions. this is the
+ // case in tests/deal.II/no_flux_11.
Point<dim> normal_vector
= (cell->face(face_no)->get_boundary()
.normal_vector (cell->face(face_no),
fe_values.quadrature_point(i)));
+ if (normal_vector * fe_values.normal_vector(i) < 0)
+ normal_vector *= -1;
Assert (std::fabs(normal_vector.norm() - 1) < 1e-14,
ExcInternalError());
for (unsigned int d=0; d<dim; ++d)
if (std::fabs(normal_vector[d]) < 1e-13)
normal_vector[d] = 0;
normal_vector /= normal_vector.norm();
-
+
// now enter the (dofs,(normal_vector,cell)) entry into
// the map
dof_to_normals_map
double exponent = exponent_1;
const unsigned int n_components = dof.get_fe().n_components();
- const bool fe_is_system = (n_components != 1);
if (weight!=0)
{
const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+ AssertThrow(cell_point.first->is_locally_owned(),
+ ExcPointNotAvailableHere());
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
ExcInternalError());
cell_point
= GridTools::find_active_cell_around_point (mapping, dof, point);
+ AssertThrow(cell_point.first->is_locally_owned(),
+ ExcPointNotAvailableHere());
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
ExcInternalError());
const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+ AssertThrow(cell_point.first->is_locally_owned(),
+ ExcPointNotAvailableHere());
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
ExcInternalError());
const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+ AssertThrow(cell_point.first->is_locally_owned(),
+ ExcPointNotAvailableHere());
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
ExcInternalError());
const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+ AssertThrow(cell_point.first->is_locally_owned(),
+ ExcPointNotAvailableHere());
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
ExcInternalError());
deal_ii_objects_${build_lowercase}
)
ADD_LIBRARY(${DEAL_II_BASE_NAME}${DEAL_II_${build}_SUFFIX}
- base/dummy.cc # Workaround for a bug in the Xcode generator
+ dummy.cc # Workaround for a bug in the Xcode generator
${deal_ii_objects_${build_lowercase}}
)
COMPILE_DEFINITIONS "${DEAL_II_DEFINITIONS};${DEAL_II_DEFINITIONS_${build}}"
COMPILE_FLAGS "${DEAL_II_CXX_FLAGS_${build}}"
INSTALL_NAME_DIR "${CMAKE_INSTALL_PREFIX}/${DEAL_II_LIBRARY_RELDIR}"
+ ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+ LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+ RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
)
TARGET_LINK_LIBRARIES(${DEAL_II_BASE_NAME}${DEAL_II_${build}_SUFFIX}
${DEAL_II_EXTERNAL_LIBRARIES}
)
+ FILE(MAKE_DIRECTORY
+ ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}
+ )
+ EXPORT(TARGETS ${DEAL_II_BASE_NAME}${DEAL_II_${build}_SUFFIX}
+ FILE ${CMAKE_BINARY_DIR}/${DEAL_II_PROJECT_CONFIG_RELDIR}/${DEAL_II_PROJECT_CONFIG_NAME}BuildTargets.cmake
+ APPEND
+ )
+
INSTALL(TARGETS ${DEAL_II_BASE_NAME}${DEAL_II_${build}_SUFFIX}
EXPORT ${DEAL_II_PROJECT_CONFIG_NAME}Targets
DESTINATION ${DEAL_II_LIBRARY_RELDIR}
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/algorithms/*.h
- )
-
SET(_src
operator.cc
timestep_control.cc
operator.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/algorithms/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_algorithms OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_algorithms "${_inst}")
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/base/*.h
- )
-
SET(_src
auto_derivative_function.cc
conditional_ostream.cc
config.cc
convergence_table.cc
data_out_base.cc
- # do NOT include dummy.cc
event.cc
exceptions.cc
flow_function.cc
polynomials_adini.cc
polynomials_bdm.cc
polynomials_nedelec.cc
- polynomials_piecewise.cc
polynomial_space.cc
polynomials_p.cc
+ polynomials_piecewise.cc
polynomials_raviart_thomas.cc
quadrature.cc
quadrature_lib.cc
utilities.cc
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/base/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_base OBJECT ${_src} ${_header}
${CMAKE_BINARY_DIR}/include/deal.II/base/config.h
)
-
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/distributed/*.h
- )
-
SET(_src
grid_refinement.cc
solution_transfer.cc
tria.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/distributed/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_distributed OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_distributed "${_inst}")
-
# include <p4est_vtk.h>
# include <p4est_ghost.h>
# include <p4est_communication.h>
+# include <p4est_iterate.h>
# include <p8est_bits.h>
# include <p8est_extended.h>
# include <p8est_vtk.h>
# include <p8est_ghost.h>
# include <p8est_communication.h>
+# include <p8est_iterate.h>
#endif
#include <algorithm>
namespace p4est
{
/**
- * A structure whose explicit
- * specializations contain
- * pointers to the relevant
- * p4est_* and p8est_*
- * functions. Using this
- * structure, for example by
- * saying
- * functions<dim>::quadrant_compare,
- * we can write code in a
- * dimension independent way,
- * either calling
- * p4est_quadrant_compare or
- * p8est_quadrant_compare,
- * depending on template
- * argument.
+ * A structure whose explicit specializations contain pointers to the
+ * relevant p4est_* and p8est_* functions. Using this structure, for
+ * example by saying functions<dim>::quadrant_compare, we can write code
+ * in a dimension independent way, either calling p4est_quadrant_compare
+ * or p8est_quadrant_compare, depending on template argument.
*/
template <int dim> struct functions;
int (&quadrant_is_ancestor) (const types<2>::quadrant *q1,
const types<2>::quadrant *q2);
+ static
+ int (&quadrant_ancestor_id) (const types<2>::quadrant *q,
+ int level);
+
static
int (&comm_find_owner) (types<2>::forest *p4est,
const types<2>::locidx which_tree,
types<2>::topidx num_corners,
types<2>::topidx num_vtt);
+ static
+ void (&connectivity_join_faces) (types<2>::connectivity *conn,
+ types<2>::topidx tree_left,
+ types<2>::topidx tree_right,
+ int face_left,
+ int face_right,
+ int orientation);
+
+
+
static
void (&connectivity_destroy) (p4est_connectivity_t *connectivity);
void (&connectivity_save) (const char *filename,
types<2>::connectivity *connectivity);
+ static
+ int (&connectivity_is_valid) (types<2>::connectivity *connectivity);
+
static
types<2>::connectivity *(&connectivity_load) (const char *filename,
long *length);
const types<2>::quadrant *q2)
= p4est_quadrant_is_ancestor;
+ int (&functions<2>::quadrant_ancestor_id) (const types<2>::quadrant *q,
+ int level)
+ = p4est_quadrant_ancestor_id;
+
int (&functions<2>::comm_find_owner) (types<2>::forest *p4est,
const types<2>::locidx which_tree,
const types<2>::quadrant *q,
types<2>::topidx num_vtt)
= p4est_connectivity_new;
+#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+ void (&functions<2>::connectivity_join_faces) (types<2>::connectivity *conn,
+ types<2>::topidx tree_left,
+ types<2>::topidx tree_right,
+ int face_left,
+ int face_right,
+ int orientation)
+ = p4est_connectivity_join_faces;
+#endif
+
void (&functions<2>::connectivity_destroy) (p4est_connectivity_t *connectivity)
= p4est_connectivity_destroy;
types<2>::connectivity *connectivity)
= p4est_connectivity_save;
+ int (&functions<2>::connectivity_is_valid) (types<2>::connectivity
+ *connectivity)
+ = p4est_connectivity_is_valid;
+
types<2>::connectivity *
(&functions<2>::connectivity_load) (const char *filename,
long *length)
int (&quadrant_is_ancestor) (const types<3>::quadrant *q1,
const types<3>::quadrant *q2);
+ static
+ int (&quadrant_ancestor_id) (const types<3>::quadrant *q,
+ int level);
+
static
int (&comm_find_owner) (types<3>::forest *p4est,
const types<3>::locidx which_tree,
types<3>::topidx num_corners,
types<3>::topidx num_ctt);
+ static
+ void (&connectivity_join_faces) (types<3>::connectivity *conn,
+ types<3>::topidx tree_left,
+ types<3>::topidx tree_right,
+ int face_left,
+ int face_right,
+ int orientation);
+
static
void (&connectivity_destroy) (p8est_connectivity_t *connectivity);
void (&connectivity_save) (const char *filename,
types<3>::connectivity *connectivity);
+ static
+ int (&connectivity_is_valid) (types<3>::connectivity *connectivity);
+
static
types<3>::connectivity *(&connectivity_load) (const char *filename,
long *length);
const types<3>::quadrant *q2)
= p8est_quadrant_is_ancestor;
+ int (&functions<3>::quadrant_ancestor_id) (const types<3>::quadrant *q,
+ int level)
+ = p8est_quadrant_ancestor_id;
+
int (&functions<3>::comm_find_owner) (types<3>::forest *p4est,
const types<3>::locidx which_tree,
const types<3>::quadrant *q,
void (&functions<3>::connectivity_destroy) (p8est_connectivity_t *connectivity)
= p8est_connectivity_destroy;
+#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+ void (&functions<3>::connectivity_join_faces) (types<3>::connectivity *conn,
+ types<3>::topidx tree_left,
+ types<3>::topidx tree_right,
+ int face_left,
+ int face_right,
+ int orientation)
+ = p8est_connectivity_join_faces;
+#endif
+
types<3>::forest *(&functions<3>::new_forest) (MPI_Comm mpicomm,
types<3>::connectivity *connectivity,
types<3>::locidx min_quadrants,
types<3>::connectivity *connectivity)
= p8est_connectivity_save;
+ int (&functions<3>::connectivity_is_valid) (types<3>::connectivity
+ *connectivity)
+ = p8est_connectivity_is_valid;
+
types<3>::connectivity *
(&functions<3>::connectivity_load) (const char *filename,
long *length)
{
return functions<dim>::quadrant_is_ancestor(&q1, &q2);
}
+
+ /**
+ * This struct templatizes the p4est iterate structs and function
+ * prototypes, which are used to execute callback functions for faces,
+ * edges, and corners that require local neighborhood information, i.e.
+ * the neighboring cells */
+ template <int dim> struct iter;
+
+ template <> struct iter<2>
+ {
+ typedef p4est_iter_corner_info_t corner_info;
+ typedef p4est_iter_corner_side_t corner_side;
+ typedef p4est_iter_corner_t corner_iter;
+ typedef p4est_iter_face_info_t face_info;
+ typedef p4est_iter_face_side_t face_side;
+ typedef p4est_iter_face_t face_iter;
+ };
+
+ template <> struct iter<3>
+ {
+ typedef p8est_iter_corner_info_t corner_info;
+ typedef p8est_iter_corner_side_t corner_side;
+ typedef p8est_iter_corner_t corner_iter;
+ typedef p8est_iter_edge_info_t edge_info;
+ typedef p8est_iter_edge_side_t edge_side;
+ typedef p8est_iter_edge_t edge_iter;
+ typedef p8est_iter_face_info_t face_info;
+ typedef p8est_iter_face_side_t face_side;
+ typedef p8est_iter_face_t face_iter;
+ };
+
}
}
/**
- * Set all vertex and cell
- * related information in the
- * p4est connectivity structure.
+ * Set all vertex and cell related information in the p4est connectivity
+ * structure.
*/
template <int dim, int spacedim>
void
const bool set_vertex_info,
typename internal::p4est::types<dim>::connectivity *connectivity)
{
- // copy the vertices into the
- // connectivity structure. the
- // triangulation exports the array of
- // vertices, but some of the entries are
- // sometimes unused; this shouldn't be
- // the case for a newly created
- // triangulation, but make sure
+ // copy the vertices into the connectivity structure. the triangulation
+ // exports the array of vertices, but some of the entries are sometimes
+ // unused; this shouldn't be the case for a newly created triangulation,
+ // but make sure
//
- // note that p4est stores coordinates as
- // a triplet of values even in 2d
+ // note that p4est stores coordinates as a triplet of values even in 2d
Assert (triangulation.get_used_vertices().size() ==
triangulation.get_vertices().size(),
ExcInternalError());
triangulation.get_vertices()[v][2]);
}
- // next store the tree_to_vertex indices
- // (each tree is here only a single cell
- // in the coarse mesh). p4est requires
- // vertex numbering in clockwise
+ // next store the tree_to_vertex indices (each tree is here only a single
+ // cell in the coarse mesh). p4est requires vertex numbering in clockwise
// orientation
//
- // while we're at it, also copy the
- // neighborship information between cells
+ // while we're at it, also copy the neighborship information between cells
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
connectivity->tree_to_corner[index*GeometryInfo<dim>::vertices_per_cell+v] = cell->vertex_index(v);
}
- // neighborship information. if a
- // cell is at a boundary, then enter
+ // neighborship information. if a cell is at a boundary, then enter
// the index of the cell itself here
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary() == false)
connectivity->tree_to_tree[index*GeometryInfo<dim>::faces_per_cell + f]
= coarse_cell_to_p4est_tree_permutation[cell->index()];
- // fill tree_to_face, which is
- // essentially neighbor_to_neighbor;
- // however, we have to remap the
- // resulting face number as well
+ // fill tree_to_face, which is essentially neighbor_to_neighbor;
+ // however, we have to remap the resulting face number as well
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary() == false)
{
case 3:
{
/*
- * The values for
- * tree_to_face are in
- * 0..23 where ttf % 6
- * gives the face
- * number and ttf / 6
- * the face orientation
- * code. The
- * orientation is
- * determined as
- * follows. Let
- * my_face and
- * other_face be the
- * two face numbers of
- * the connecting trees
- * in 0..5. Then the
- * first face vertex of
- * the lower of my_face
- * and other_face
- * connects to a face
- * vertex numbered 0..3
- * in the higher of
- * my_face and
- * other_face. The
- * face orientation is
- * defined as this
- * number. If my_face
- * == other_face,
- * treating either of
- * both faces as the
- * lower one leads to
- * the same result.
+ * The values for tree_to_face are in 0..23 where ttf % 6
+ * gives the face number and ttf / 6 the face orientation
+ * code. The orientation is determined as follows. Let
+ * my_face and other_face be the two face numbers of the
+ * connecting trees in 0..5. Then the first face vertex of
+ * the lower of my_face and other_face connects to a face
+ * vertex numbered 0..3 in the higher of my_face and
+ * other_face. The face orientation is defined as this
+ * number. If my_face == other_face, treating either of
+ * both faces as the lower one leads to the same result.
*/
connectivity->tree_to_face[index*6 + f]
unsigned int v = 0;
- // global vertex index of
- // vertex 0 on face of
- // cell with smaller
- // local face index
+ // global vertex index of vertex 0 on face of cell with
+ // smaller local face index
unsigned int g_idx =
cell_list[smaller_idx]->vertex_index(
GeometryInfo<dim>::face_to_cell_vertices(
cell_list[smaller_idx]->face_rotation(face_idx_list[smaller_idx]))
);
- // loop over vertices on
- // face from other cell
- // and compare global
- // vertex numbers
+ // loop over vertices on face from other cell and compare
+ // global vertex numbers
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
{
unsigned int idx
const typename internal::p4est::types<dim>::forest &forest,
const types::subdomain_id my_subdomain)
{
- // check if this cell exists in
- // the local p4est cell
+ // check if this cell exists in the local p4est cell
if (sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
&p4est_cell,
internal::p4est::functions<dim>::quadrant_compare)
}
else
{
- // no, cell not found in
- // local part of
- // p4est. this means that
- // the local part is more
- // refined than the current
- // cell. if this cell has
- // no children of its own,
- // we need to refine it,
- // and if it does already
- // have children then loop
- // over all children and
- // see if they are locally
- // available as well
+ // no, cell not found in local part of p4est. this means that the
+ // local part is more refined than the current cell. if this cell has
+ // no children of its own, we need to refine it, and if it does
+ // already have children then loop over all children and see if they
+ // are locally available as well
if (dealii_cell->has_children () == false)
dealii_cell->set_refine_flag ();
else
&p4est_child[c])
== false)
{
- // no, this child
- // is locally not
- // available in
- // the p4est.
- // delete all
- // its children
- // but, because this
- // may not be successfull,
- // make sure to mark all
- // children recursively as
- // not local.
+ // no, this child is locally not available in the p4est.
+ // delete all its children but, because this may not be
+ // successfull, make sure to mark all children recursively
+ // as not local.
delete_all_children<dim,spacedim> (dealii_cell->child(c));
dealii_cell->child(c)
->recursively_set_subdomain_id(numbers::artificial_subdomain_id);
}
else
{
- // at least some
- // part of the
- // tree rooted in
- // this child is
- // locally
- // available
+ // at least some part of the tree rooted in this child is
+ // locally available
match_tree_recursively<dim,spacedim> (tree,
dealii_cell->child(c),
p4est_child[c],
template <int dim, int spacedim>
void
- match_quadrant_recursively (const typename internal::p4est::types<dim>::quadrant &this_quadrant,
- const typename Triangulation<dim,spacedim>::cell_iterator &dealii_cell,
- const typename internal::p4est::types<dim>::quadrant &ghost_quadrant,
- const typename internal::p4est::types<dim>::forest &forest,
- unsigned int ghost_owner)
+ match_quadrant (const dealii::Triangulation<dim,spacedim> *tria,
+ unsigned int dealii_index,
+ typename internal::p4est::types<dim>::quadrant &ghost_quadrant,
+ unsigned int ghost_owner)
{
- if (internal::p4est::functions<dim>::quadrant_is_equal(&this_quadrant, &ghost_quadrant))
- {
- // this is the ghostcell
+ int i, child_id;
+ int l = ghost_quadrant.level;
- if (dealii_cell->has_children())
- delete_all_children<dim,spacedim> (dealii_cell);
- else
+ for (i = 0; i < l; i++)
+ {
+ typename Triangulation<dim,spacedim>::cell_iterator cell (tria, i, dealii_index);
+ if (cell->has_children () == false)
{
- dealii_cell->clear_coarsen_flag();
- dealii_cell->set_subdomain_id(ghost_owner);
+ cell->clear_coarsen_flag();
+ cell->set_refine_flag ();
+ return;
}
- return;
- }
- if (! internal::p4est::functions<dim>::quadrant_is_ancestor ( &this_quadrant, &ghost_quadrant))
- {
- return;
+ child_id = internal::p4est::functions<dim>::quadrant_ancestor_id (&ghost_quadrant, i + 1);
+ dealii_index = cell->child_index(child_id);
}
- if (dealii_cell->has_children () == false)
- {
- dealii_cell->clear_coarsen_flag();
- dealii_cell->set_refine_flag ();
- return;
- }
-
- typename internal::p4est::types<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
- for (unsigned int c=0;
- c<GeometryInfo<dim>::max_children_per_cell; ++c)
- switch (dim)
- {
- case 2:
- P4EST_QUADRANT_INIT(&p4est_child[c]);
- break;
- case 3:
- P8EST_QUADRANT_INIT(&p4est_child[c]);
- break;
- default:
- Assert (false, ExcNotImplemented());
- }
-
-
- internal::p4est::functions<dim>::
- quadrant_childrenv (&this_quadrant, p4est_child);
-
- for (unsigned int c=0;
- c<GeometryInfo<dim>::max_children_per_cell; ++c)
-
+ typename Triangulation<dim,spacedim>::cell_iterator cell (tria, l, dealii_index);
+ if (cell->has_children())
+ delete_all_children<dim,spacedim> (cell);
+ else
{
- match_quadrant_recursively<dim,spacedim> (p4est_child[c], dealii_cell->child(c), ghost_quadrant, forest, ghost_owner);
+ cell->clear_coarsen_flag();
+ cell->set_subdomain_id(ghost_owner);
}
-
-
}
{
//this cell got refined
- //attach to the first child, because
- // we can only attach to active
+ //attach to the first child, because we can only attach to active
// quadrants
typename internal::p4est::types<dim>::quadrant
p4est_child[GeometryInfo<dim>::max_children_per_cell];
}
else
{
- //it's children got coarsened into
- //this cell
+ //it's children got coarsened into this cell
typename internal::p4est::types<dim>::quadrant *q;
q = static_cast<typename internal::p4est::types<dim>::quadrant *> (
sc_array_index (const_cast<sc_array_t *>(&tree.quadrants), idx)
quadrant_overlaps_tree (const_cast<typename internal::p4est::types<dim>::tree *>(&tree),
&p4est_cell)
== false))
- // this quadrant and none of
- // it's children belong to us.
+ // this quadrant and none of it's children belong to us.
return;
/**
- * A data structure that we use to store
- * which cells (indicated by
- * internal::p4est::types<dim>::quadrant objects) shall be
- * refined and which shall be coarsened.
+ * A data structure that we use to store which cells (indicated by
+ * internal::p4est::types<dim>::quadrant objects) shall be refined and which
+ * shall be coarsened.
*/
template <int dim, int spacedim>
class RefineAndCoarsenList
typename internal::p4est::types<dim>::forest &forest);
/**
- * A callback function that we
- * pass to the p4est data
- * structures when a forest is
- * to be refined. The p4est
- * functions call it back with
- * a tree (the index of the
- * tree that grows out of a
- * given coarse cell) and a
- * refinement path from that
- * coarse cell to a
- * terminal/leaf cell. The
- * function returns whether the
- * corresponding cell in the
- * deal.II triangulation has
- * the refined flag set.
+ * A callback function that we pass to the p4est data structures when a
+ * forest is to be refined. The p4est functions call it back with a tree
+ * (the index of the tree that grows out of a given coarse cell) and a
+ * refinement path from that coarse cell to a terminal/leaf cell. The
+ * function returns whether the corresponding cell in the deal.II
+ * triangulation has the refined flag set.
*/
static
int
typename internal::p4est::types<dim>::quadrant *quadrant);
/**
- * Same as the refine_callback
- * function, but return whether
- * all four of the given
- * children of a non-terminal
- * cell are to be coarsened
- * away.
+ * Same as the refine_callback function, but return whether all four of
+ * the given children of a non-terminal cell are to be coarsened away.
*/
static
int
template <int dim, int spacedim>
RefineAndCoarsenList<dim,spacedim>::
RefineAndCoarsenList (const Triangulation<dim,spacedim> &triangulation,
- const std::vector<types::global_dof_index> &p4est_tree_to_coarse_cell_permutation,
- const types::subdomain_id my_subdomain,
+ const std::vector<types::global_dof_index> &p4est_tree_to_coarse_cell_permutation,
+ const types::subdomain_id my_subdomain,
typename internal::p4est::types<dim>::forest &forest)
:
forest(forest)
{
- // count how many flags are set and
- // allocate that much memory
+ // count how many flags are set and allocate that much memory
unsigned int n_refine_flags = 0,
n_coarsen_flags = 0;
for (typename Triangulation<dim,spacedim>::active_cell_iterator
coarsen_list.reserve (n_coarsen_flags);
- // now build the lists of cells that
- // are flagged. note that p4est will
- // traverse its cells in the order in
- // which trees appear in the
- // forest. this order is not the same
- // as the order of coarse cells in the
- // deal.II Triangulation because we
- // have translated everything by the
- // coarse_cell_to_p4est_tree_permutation
- // permutation. in order to make sure
- // that the output array is already in
- // the correct order, traverse our
- // coarse cells in the same order in
- // which p4est will:
+ // now build the lists of cells that are flagged. note that p4est will
+ // traverse its cells in the order in which trees appear in the
+ // forest. this order is not the same as the order of coarse cells in the
+ // deal.II Triangulation because we have translated everything by the
+ // coarse_cell_to_p4est_tree_permutation permutation. in order to make
+ // sure that the output array is already in the correct order, traverse
+ // our coarse cells in the same order in which p4est will:
for (unsigned int c=0; c<triangulation.n_cells(0); ++c)
{
unsigned int coarse_cell_index =
Assert(coarsen_list.size() == n_coarsen_flags,
ExcInternalError());
- // make sure that our ordering in fact
- // worked
+ // make sure that our ordering in fact worked
for (unsigned int i=1; i<refine_list.size(); ++i)
Assert (refine_list[i].p.which_tree >=
refine_list[i-1].p.which_tree,
RefineAndCoarsenList<dim,spacedim> *this_object
= reinterpret_cast<RefineAndCoarsenList<dim,spacedim>*>(forest->user_pointer);
- // if there are no more cells in our
- // list the current cell can't be
+ // if there are no more cells in our list the current cell can't be
// flagged for refinement
if (this_object->current_refine_pointer == this_object->refine_list.end())
return false;
this_object->current_refine_pointer->p.which_tree,
ExcInternalError());
- // if p4est hasn't yet reached the tree
- // of the next flagged cell the current
- // cell can't be flagged for refinement
+ // if p4est hasn't yet reached the tree of the next flagged cell the
+ // current cell can't be flagged for refinement
if (coarse_cell_index <
this_object->current_refine_pointer->p.which_tree)
return false;
- // now we're in the right tree in the
- // forest
+ // now we're in the right tree in the forest
Assert (coarse_cell_index <=
this_object->current_refine_pointer->p.which_tree,
ExcInternalError());
- // make sure that the p4est loop over
- // cells hasn't gotten ahead of our own
+ // make sure that the p4est loop over cells hasn't gotten ahead of our own
// pointer
Assert (internal::p4est::functions<dim>::
quadrant_compare (quadrant,
&*this_object->current_refine_pointer) <= 0,
ExcInternalError());
- // now, if the p4est cell is one in the
- // list, it is supposed to be refined
+ // now, if the p4est cell is one in the list, it is supposed to be refined
if (internal::p4est::functions<dim>::
quadrant_is_equal (quadrant, &*this_object->current_refine_pointer))
{
RefineAndCoarsenList<dim,spacedim> *this_object
= reinterpret_cast<RefineAndCoarsenList<dim,spacedim>*>(forest->user_pointer);
- // if there are no more cells in our
- // list the current cell can't be
+ // if there are no more cells in our list the current cell can't be
// flagged for coarsening
if (this_object->current_coarsen_pointer ==
this_object->coarsen_list.end())
this_object->current_coarsen_pointer->p.which_tree,
ExcInternalError());
- // if p4est hasn't yet reached the tree
- // of the next flagged cell the current
- // cell can't be flagged for coarsening
+ // if p4est hasn't yet reached the tree of the next flagged cell the
+ // current cell can't be flagged for coarsening
if (coarse_cell_index <
this_object->current_coarsen_pointer->p.which_tree)
return false;
- // now we're in the right tree in the
- // forest
+ // now we're in the right tree in the forest
Assert (coarse_cell_index <=
this_object->current_coarsen_pointer->p.which_tree,
ExcInternalError());
- // make sure that the p4est loop over
- // cells hasn't gotten ahead of our own
+ // make sure that the p4est loop over cells hasn't gotten ahead of our own
// pointer
Assert (internal::p4est::functions<dim>::
quadrant_compare (children[0],
&*this_object->current_coarsen_pointer) <= 0,
ExcInternalError());
- // now, if the p4est cell is one in the
- // list, it is supposed to be coarsened
+ // now, if the p4est cell is one in the list, it is supposed to be
+ // coarsened
if (internal::p4est::functions<dim>::
quadrant_is_equal (children[0],
&*this_object->current_coarsen_pointer))
// move current pointer one up
++this_object->current_coarsen_pointer;
- // note that the next 3 cells in
- // our list need to correspond to
- // the other siblings of the cell
- // we have just found
+ // note that the next 3 cells in our list need to correspond to the
+ // other siblings of the cell we have just found
for (unsigned int c=1; c<GeometryInfo<dim>::max_children_per_cell; ++c)
{
Assert (internal::p4est::functions<dim>::
{
Singleton ()
{
- // ensure that the
- // initialization
- // code is run only
- // once, even if we
- // link with 1d, 2d,
- // and 3d libraries
+ // ensure that the initialization code is run only once, even if we
+ // link with 1d, 2d, and 3d libraries
static bool initialized = false;
if (initialized == false)
if (deinitialized == false)
{
- // p4est has no
- // p4est_finalize
- // function
+ // p4est has no p4est_finalize function
sc_finalize ();
deinitialized = true;
};
public:
- // do run the initialization code, at least the first time around we
- // get to this function
+ // do run the initialization code, at least the first time around we get
+ // to this function
static void do_initialize ()
{
static Singleton singleton;
const typename dealii::Triangulation<dim,spacedim>::MeshSmoothing smooth_grid,
const Settings settings_)
:
- // do not check
- // for distorted
- // cells
+ // do not check for distorted cells
dealii::Triangulation<dim,spacedim>
(smooth_grid,
false),
n_attached_datas(0),
n_attached_deserialize(0)
{
- // initialize p4est. do this in a separate function since it has
- // to happen only once, even if we have triangulation objects
- // for several different space dimensions
+ // initialize p4est. do this in a separate function since it has to
+ // happen only once, even if we have triangulation objects for several
+ // different space dimensions
dealii::internal::p4est::InitFinalize::do_initialize ();
+ parallel_ghost = 0;
+
number_cache.n_locally_owned_active_cells
.resize (Utilities::MPI::n_mpi_processes (mpi_communicator));
}
Assert (parallel_forest == 0, ExcInternalError());
Assert (refinement_in_progress == false, ExcInternalError());
- // get rid of the unique
- // communicator used here again
+ // get rid of the unique communicator used here again
MPI_Comm_free (&mpi_communicator);
}
}
catch (const typename dealii::Triangulation<dim,spacedim>::DistortedCellList &)
{
- // the underlying
- // triangulation should not
- // be checking for
- // distorted cells
+ // the underlying triangulation should not be checking for distorted
+ // cells
AssertThrow (false, ExcInternalError());
}
- // note that now we have some content in
- // the p4est objects and call the
- // functions that do the actual work
- // (which are dimension dependent, so
+ // note that now we have some content in the p4est objects and call the
+ // functions that do the actual work (which are dimension dependent, so
// separate)
triangulation_has_content = true;
}
catch (const typename Triangulation<dim>::DistortedCellList &)
{
- // the underlying
- // triangulation should not
- // be checking for
- // distorted cells
+ // the underlying triangulation should not be checking for distorted
+ // cells
AssertThrow (false, ExcInternalError());
}
{
triangulation_has_content = false;
+ if (parallel_ghost != 0)
+ {
+ dealii::internal::p4est::functions<dim>::ghost_destroy (parallel_ghost);
+ parallel_ghost = 0;
+ }
+
if (parallel_forest != 0)
{
dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
Triangulation<dim,spacedim>::
load(const char *filename)
{
+ if (parallel_ghost != 0) {
+ dealii::internal::p4est::functions<dim>::ghost_destroy (parallel_ghost);
+ parallel_ghost = 0;
+ }
dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
parallel_forest = 0;
dealii::internal::p4est::functions<dim>::connectivity_destroy (connectivity);
Assert (this->n_cells(0) > 0, ExcInternalError());
Assert (this->n_levels() == 1, ExcInternalError());
- // data structures that counts how many
- // cells touch each vertex
- // (vertex_touch_count), and which cells
- // touch a given vertex (together with
- // the local numbering of that vertex
- // within the cells that touch it)
+ // data structures that counts how many cells touch each vertex
+ // (vertex_touch_count), and which cells touch a given vertex (together
+ // with the local numbering of that vertex within the cells that touch
+ // it)
std::vector<unsigned int> vertex_touch_count;
std::vector<
std::list<
vertex_touch_count.end(),
0);
- // now create a connectivity
- // object with the right sizes
- // for all arrays. set vertex
- // information only in debug
- // mode (saves a few bytes in
- // optimized mode)
+ // now create a connectivity object with the right sizes for all
+ // arrays. set vertex information only in debug mode (saves a few bytes
+ // in optimized mode)
const bool set_vertex_info
#ifdef DEBUG
= true
Assert (p4est_connectivity_is_valid (connectivity) == 1,
ExcInternalError());
- // now create a forest out of the
- // connectivity data structure
+ // now create a forest out of the connectivity data structure
parallel_forest
= dealii::internal::p4est::functions<2>::
new_forest (mpi_communicator,
- // TODO: This is a verbatim copy of the 2,2
- // case. However, we can't just specialize the
- // dim template argument, but let spacedim open
+ // TODO: This is a verbatim copy of the 2,2 case. However, we can't just
+ // specialize the dim template argument, but let spacedim open
template <>
void
Triangulation<2,3>::copy_new_triangulation_to_p4est (dealii::internal::int2type<2>)
Assert (this->n_cells(0) > 0, ExcInternalError());
Assert (this->n_levels() == 1, ExcInternalError());
- // data structures that counts how many
- // cells touch each vertex
- // (vertex_touch_count), and which cells
- // touch a given vertex (together with
- // the local numbering of that vertex
- // within the cells that touch it)
+ // data structures that counts how many cells touch each vertex
+ // (vertex_touch_count), and which cells touch a given vertex (together
+ // with the local numbering of that vertex within the cells that touch
+ // it)
std::vector<unsigned int> vertex_touch_count;
std::vector<
std::list<
vertex_touch_count.end(),
0);
- // now create a connectivity
- // object with the right sizes
- // for all arrays. set vertex
- // information only in debug
- // mode (saves a few bytes in
- // optimized mode)
+ // now create a connectivity object with the right sizes for all
+ // arrays. set vertex information only in debug mode (saves a few bytes
+ // in optimized mode)
const bool set_vertex_info
#ifdef DEBUG
= true
Assert (p4est_connectivity_is_valid (connectivity) == 1,
ExcInternalError());
- // now create a forest out of the
- // connectivity data structure
+ // now create a forest out of the connectivity data structure
parallel_forest
= dealii::internal::p4est::functions<2>::
new_forest (mpi_communicator,
Assert (this->n_cells(0) > 0, ExcInternalError());
Assert (this->n_levels() == 1, ExcInternalError());
- // data structures that counts how many
- // cells touch each vertex
- // (vertex_touch_count), and which cells
- // touch a given vertex (together with
- // the local numbering of that vertex
- // within the cells that touch it)
+ // data structures that counts how many cells touch each vertex
+ // (vertex_touch_count), and which cells touch a given vertex (together
+ // with the local numbering of that vertex within the cells that touch
+ // it)
std::vector<unsigned int> vertex_touch_count;
std::vector<
std::list<
num_ett = std::accumulate (edge_touch_count.begin(),
edge_touch_count.end(),
0);
- // now create a connectivity object with
- // the right sizes for all arrays
-
+
+ // now create a connectivity object with the right sizes for all arrays
const bool set_vertex_info
#ifdef DEBUG
= true
Assert (p8est_connectivity_is_valid (connectivity) == 1,
ExcInternalError());
- // now create a forest out of the
- // connectivity data structure
+ // now create a forest out of the connectivity data structure
parallel_forest
= dealii::internal::p4est::functions<3>::
new_forest (mpi_communicator,
void
Triangulation<dim,spacedim>::copy_local_forest_to_triangulation ()
{
- // disable mesh smoothing for
- // recreating the deal.II
- // triangulation, otherwise we might
- // not be able to reproduce the p4est
- // mesh exactly. We restore the
- // original smoothing at the end of
- // this function. Note that the
- // smoothing flag is used in the normal
+ // disable mesh smoothing for recreating the deal.II triangulation,
+ // otherwise we might not be able to reproduce the p4est mesh
+ // exactly. We restore the original smoothing at the end of this
+ // function. Note that the smoothing flag is used in the normal
// refinement process.
typename Triangulation<dim,spacedim>::MeshSmoothing
- save_smooth = this->smooth_grid;
- this->smooth_grid = dealii::Triangulation<dim,spacedim>::none;
- bool mesh_changed = false;
+ save_smooth = this->smooth_grid;
+
+ // We will refine manually to match the p4est further down, which
+ // obeys a level difference of 2 at each vertex (see the balance call
+ // to p4est). We can disable this here so we store fewer artificial
+ // cells (in some cases). For geometric multigrid it turns out that
+ // we will miss level cells at shared vertices if we ignore this.
+ // See tests/mpi/mg_06.
+ if (settings & construct_multigrid_hierarchy)
+ this->smooth_grid = dealii::Triangulation<dim,spacedim>::limit_level_difference_at_vertices;
+ else
+ this->smooth_grid = dealii::Triangulation<dim,spacedim>::none;
+ bool mesh_changed = false;
// remove all deal.II refinements. Note that we could skip this and
// start from our current state, because the algorithm later coarsens as
}
catch (const typename Triangulation<dim, spacedim>::DistortedCellList &)
{
- // the underlying
- // triangulation should not
- // be checking for
+ // the underlying triangulation should not be checking for
// distorted cells
AssertThrow (false, ExcInternalError());
}
// query p4est for the ghost cells
- typename dealii::internal::p4est::types<dim>::ghost *ghostlayer;
- ghostlayer = dealii::internal::p4est::functions<dim>::ghost_new (parallel_forest,
+ if (parallel_ghost != 0)
+ {
+ dealii::internal::p4est::functions<dim>::ghost_destroy (parallel_ghost);
+ parallel_ghost = 0;
+ }
+ parallel_ghost = dealii::internal::p4est::functions<dim>::ghost_new (parallel_forest,
(dim == 2
?
typename dealii::internal::p4est::types<dim>::
- balance_type(P4EST_BALANCE_CORNER)
+ balance_type(P4EST_CONNECT_CORNER)
:
typename dealii::internal::p4est::types<dim>::
- balance_type(P8EST_BALANCE_CORNER)));
+ balance_type(P8EST_CONNECT_CORNER)));
- Assert (ghostlayer, ExcInternalError());
+ Assert (parallel_ghost, ExcInternalError());
// set all cells to artificial. we will later set it to the correct
cell != this->end(0);
++cell)
{
- // if this processor
- // stores no part of the
- // forest that comes out
- // of this coarse grid
- // cell, then we need to
- // delete all children of
- // this cell (the coarse
- // grid cell remains)
+ // if this processor stores no part of the forest that comes out
+ // of this coarse grid cell, then we need to delete all children
+ // of this cell (the coarse grid cell remains)
if (tree_exists_locally<dim,spacedim>(parallel_forest,
coarse_cell_to_p4est_tree_permutation[cell->index()])
== false)
else
{
-
- // this processor
- // stores at least a
- // part of the tree
- // that comes out of
- // this cell.
-
+ // this processor stores at least a part of the tree that
+ // comes out of this cell.
typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
typename dealii::internal::p4est::types<dim>::tree *tree =
}
}
- // check mesh for ghostcells,
- // refine as necessary.
- // iterate over every ghostquadrant,
- // find corresponding deal coarsecell
- // and recurse.
+ // check mesh for ghostcells, refine as necessary. iterate over
+ // every ghostquadrant, find corresponding deal coarsecell and
+ // recurse.
typename dealii::internal::p4est::types<dim>::quadrant *quadr;
unsigned int ghost_owner=0;
+ typename dealii::internal::p4est::types<dim>::topidx ghost_tree=0;
- for (unsigned int g_idx=0; g_idx<ghostlayer->ghosts.elem_count; ++g_idx)
+ for (unsigned int g_idx=0; g_idx<parallel_ghost->ghosts.elem_count; ++g_idx)
{
- while (g_idx >= (unsigned int)ghostlayer->proc_offsets[ghost_owner+1])
+ while (g_idx >= (unsigned int)parallel_ghost->proc_offsets[ghost_owner+1])
++ghost_owner;
+ while (g_idx >= (unsigned int)parallel_ghost->tree_offsets[ghost_tree+1])
+ ++ghost_tree;
quadr = static_cast<typename dealii::internal::p4est::types<dim>::quadrant *>
- ( sc_array_index(&ghostlayer->ghosts, g_idx) );
+ ( sc_array_index(¶llel_ghost->ghosts, g_idx) );
unsigned int coarse_cell_index =
- p4est_tree_to_coarse_cell_permutation[quadr->p.piggy3.which_tree];
-
- const typename Triangulation<dim,spacedim>::cell_iterator
- cell (this, 0U, coarse_cell_index);
+ p4est_tree_to_coarse_cell_permutation[ghost_tree];
- typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
- dealii::internal::p4est::init_coarse_quadrant<dim> (p4est_coarse_cell);
-
- match_quadrant_recursively<dim,spacedim> (p4est_coarse_cell, cell, *quadr,
- *parallel_forest, ghost_owner);
+ match_quadrant<dim,spacedim> (this, coarse_cell_index, *quadr, ghost_owner);
}
- // fix all the flags to make
- // sure we have a consistent
- // mesh
+ // fix all the flags to make sure we have a consistent mesh
this->prepare_coarsening_and_refinement ();
// see if any flags are still set
break;
}
- // actually do the refinement
- // but prevent the refinement
- // hook below from taking
- // over
+ // actually do the refinement but prevent the refinement hook below
+ // from taking over
const bool saved_refinement_in_progress = refinement_in_progress;
refinement_in_progress = true;
}
catch (const typename Triangulation<dim,spacedim>::DistortedCellList &)
{
- // the underlying
- // triangulation should not
- // be checking for
+ // the underlying triangulation should not be checking for
// distorted cells
AssertThrow (false, ExcInternalError());
}
++num_ghosts;
}
- Assert( num_ghosts == ghostlayer->ghosts.elem_count, ExcInternalError());
+ Assert( num_ghosts == parallel_ghost->ghosts.elem_count, ExcInternalError());
#endif
// we need this information for all our ancestors and the same-level neighbors of our own cells (=level ghosts)
if (settings & construct_multigrid_hierarchy)
{
- // step 1: We set our own ids all the way down and all the others to -1. Note that we do not fill
- // other cells we could figure out the same way, because we might accidentally set an id for a
- // cell that is not a ghost cell.
+ // step 1: We set our own ids all the way down and all the others to
+ // -1. Note that we do not fill other cells we could figure out the
+ // same way, because we might accidentally set an id for a cell that
+ // is not a ghost cell.
for (unsigned int lvl=this->n_levels(); lvl>0;)
{
--lvl;
}
}
- //step 2: make sure all the neighbors to our level_cells exist. Need to look up in p4est...
+ //step 2: make sure all the neighbors to our level_cells exist. Need
+ //to look up in p4est...
std::vector<std::vector<bool> > marked_vertices(this->n_levels());
for (unsigned int lvl=0; lvl<this->n_levels(); ++lvl)
{
if (cell->child(c)->level_subdomain_id()==this->locally_owned_subdomain())
{
- //at least one of the children belongs to us, so make sure we set the level subdomain id
+ //at least one of the children belongs to us, so
+ //make sure we set the level subdomain id
types::subdomain_id mark = numbers::artificial_subdomain_id;
mark = cell->child(0)->level_subdomain_id();
Assert(mark != numbers::artificial_subdomain_id, ExcInternalError()); //we should know the child(0)
}
}
+ //step 4: Special case: on each level we need all the face neighbors
+ // of our own level cells these are normally on the same level,
+ // unless the neighbor is active and coarser. It can end up on a
+ // different processor. Luckily, the level_subdomain_id can be
+ // figured out without communication, because the cell is active
+ // (and so level_subdomain_id=subdomain_id):
+ for (typename Triangulation<dim,spacedim>::cell_iterator cell = this->begin(); cell!=this->end(); ++cell)
+ {
+ if (cell->level_subdomain_id()!=this->locally_owned_subdomain())
+ continue;
+
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ {
+ if (cell->face(f)->at_boundary())
+ continue;
+ if (cell->neighbor(f)->level() < cell->level()
+ &&
+ cell->neighbor(f)->level_subdomain_id()!=this->locally_owned_subdomain())
+ {
+ Assert(cell->neighbor(f)->active(), ExcInternalError());
+ Assert(cell->neighbor(f)->subdomain_id() != numbers::artificial_subdomain_id, ExcInternalError());
+ Assert(cell->neighbor(f)->level_subdomain_id() == numbers::artificial_subdomain_id
+ || cell->neighbor(f)->level_subdomain_id() == cell->neighbor(f)->subdomain_id(), ExcInternalError());
+ cell->neighbor(f)->set_level_subdomain_id(cell->neighbor(f)->subdomain_id());
+ }
+ }
+ }
}
- // check that our local copy has
- // exactly as many cells as the
- // p4est original (at least if we
- // are on only one processor);
- // for parallel computations, we
- // want to check that we have at
- // least as many as p4est stores
- // locally (in the future we
- // should check that we have
- // exactly as many non-artificial
- // cells as
- // parallel_forest->local_num_quadrants)
+ // check that our local copy has exactly as many cells as the p4est
+ // original (at least if we are on only one processor); for parallel
+ // computations, we want to check that we have at least as many as p4est
+ // stores locally (in the future we should check that we have exactly as
+ // many non-artificial cells as parallel_forest->local_num_quadrants)
{
const unsigned int total_local_cells = this->n_active_cells();
total_local_cells,
ExcInternalError());
- // count the number of owned, active
- // cells and compare with p4est.
+ // count the number of owned, active cells and compare with p4est.
unsigned int n_owned = 0;
for (typename Triangulation<dim,spacedim>::active_cell_iterator
cell = this->begin_active();
}
- dealii::internal::p4est::functions<dim>::ghost_destroy (ghostlayer);
-
this->smooth_grid = save_smooth;
}
void
Triangulation<dim, spacedim>::execute_coarsening_and_refinement ()
{
- // first make sure that recursive
- // calls are handled correctly
+ // first make sure that recursive calls are handled correctly
if (refinement_in_progress == true)
{
dealii::Triangulation<dim,spacedim>::execute_coarsening_and_refinement ();
}
}
- // now do the work we're
- // supposed to do when we are
- // in charge
+ // now do the work we're supposed to do when we are in charge
refinement_in_progress = true;
this->prepare_coarsening_and_refinement ();
- // make sure all flags are
- // cleared on cells we don't
- // own, since nothing good can
- // come of that if they are
- // still around
+ // make sure all flags are cleared on cells we don't own, since nothing
+ // good can come of that if they are still around
for (typename Triangulation<dim,spacedim>::active_cell_iterator
cell = this->begin_active();
cell != this->end(); ++cell)
}
- // count how many cells will be refined
- // and coarsened, and allocate that much
- // memory
+ // count how many cells will be refined and coarsened, and allocate that
+ // much memory
RefineAndCoarsenList<dim,spacedim>
refine_and_coarsen_list (*this,
p4est_tree_to_coarse_cell_permutation,
my_subdomain,
*parallel_forest);
- // copy refine and coarsen flags into
- // p4est and execute the refinement and
- // coarsening. this uses the
- // refine_and_coarsen_list just built,
- // which is communicated to the callback
- // functions through the
+ // copy refine and coarsen flags into p4est and execute the refinement
+ // and coarsening. this uses the refine_and_coarsen_list just built,
+ // which is communicated to the callback functions through the
// user_pointer
Assert (parallel_forest->user_pointer == this,
ExcInternalError());
parallel_forest->user_pointer = &refine_and_coarsen_list;
+ if (parallel_ghost != 0)
+ {
+ dealii::internal::p4est::functions<dim>::ghost_destroy (parallel_ghost);
+ parallel_ghost = 0;
+ }
dealii::internal::p4est::functions<dim>::
refine (parallel_forest, /* refine_recursive */ false,
&RefineAndCoarsenList<dim,spacedim>::refine_callback,
&RefineAndCoarsenList<dim,spacedim>::coarsen_callback,
/*init_callback=*/NULL);
- // make sure all cells in the lists have
- // been consumed
+ // make sure all cells in the lists have been consumed
Assert (refine_and_coarsen_list.pointers_are_at_end(),
ExcInternalError());
(dim == 2
?
typename dealii::internal::p4est::types<dim>::
- balance_type(P4EST_BALANCE_FULL)
+ balance_type(P4EST_CONNECT_FULL)
:
typename dealii::internal::p4est::types<dim>::
- balance_type(P8EST_BALANCE_FULL)),
+ balance_type(P8EST_CONNECT_FULL)),
/*init_callback=*/NULL);
- // before repartitioning the mesh let
- // others attach mesh related info
- // (such as SolutionTransfer data) to
- // the p4est
+ // before repartitioning the mesh let others attach mesh related info
+ // (such as SolutionTransfer data) to the p4est
attach_mesh_data();
- // partition the new mesh between
- // all processors
+ // partition the new mesh between all processors
dealii::internal::p4est::functions<dim>::
partition (parallel_forest,
/* prepare coarsening */ 1,
/* weight_callback */ NULL);
- // finally copy back from local
- // part of tree to deal.II
- // triangulation. before doing
- // so, make sure there are no
- // refine or coarsen flags
- // pending
+ // finally copy back from local part of tree to deal.II
+ // triangulation. before doing so, make sure there are no refine or
+ // coarsen flags pending
for (typename Triangulation<dim,spacedim>::active_cell_iterator
cell = this->begin_active();
cell != this->end(); ++cell)
}
catch (const typename Triangulation<dim>::DistortedCellList &)
{
- // the underlying
- // triangulation should not
- // be checking for
- // distorted cells
+ // the underlying triangulation should not be checking for distorted
+ // cells
AssertThrow (false, ExcInternalError());
}
number_cache.n_global_active_cells
= std::accumulate (number_cache.n_locally_owned_active_cells.begin(),
number_cache.n_locally_owned_active_cells.end(),
- /* ensure sum is
- computed with
- correct data
- type:*/
+ /* ensure sum is computed with correct data type:*/
static_cast<types::global_dof_index>(0));
number_cache.n_global_levels = Utilities::MPI::max(this->n_levels(), mpi_communicator);
}
dealii::internal::p4est::init_coarse_quadrant<dim> (p4est_coarse_cell);
- // parent_cell is not correct here,
- // but is only used in a refined
+ // parent_cell is not correct here, but is only used in a refined
// cell
post_mesh_data_recursively<dim, spacedim> (*tree,
cell,
}
// important: only remove data if we are not in the deserialization
- // process. There, each SolutionTransfer registers and unpacks
- // before the next one does this, so n_attached_datas is only 1 here.
- // This would destroy the saved data before the second SolutionTransfer
- // can get it. This created a bug that is documented in
+ // process. There, each SolutionTransfer registers and unpacks before
+ // the next one does this, so n_attached_datas is only 1 here. This
+ // would destroy the saved data before the second SolutionTransfer can
+ // get it. This created a bug that is documented in
// tests/mpi/p4est_save_03 with more than one SolutionTransfer.
if (!n_attached_datas && n_attached_deserialize == 0)
{
template <int dim, int spacedim>
- const std::vector<unsigned int> &
+ const std::vector<types::global_dof_index> &
Triangulation<dim, spacedim>::get_p4est_tree_to_coarse_cell_permutation() const
{
return p4est_tree_to_coarse_cell_permutation;
}
+ namespace
+ {
+ /**
+ * This is the callback data structure used to fill
+ * vertices_with_ghost_neighbors via the p4est_iterate tool
+ */
+ template <int dim, int spacedim>
+ struct find_ghosts
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim> *triangulation;
+ sc_array_t * subids;
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ *vertices_with_ghost_neighbors;
+ };
+ /** At a corner (vertex), determine if any of the neighboring cells are
+ * ghosts. If there are, find out their subdomain ids, and if this is a
+ * local vertex, then add these subdomain ids to the map
+ * vertices_with_ghost_neighbors of that index
+ */
+ template <int dim, int spacedim>
+ void
+ find_ghosts_corner
+ (typename dealii::internal::p4est::iter<dim>::corner_info * info,
+ void *user_data)
+ {
+ int i, j;
+ int nsides = info->sides.elem_count;
+ typename dealii::internal::p4est::iter<dim>::corner_side * sides =
+ (typename dealii::internal::p4est::iter<dim>::corner_side *)
+ (info->sides.array);
+ struct find_ghosts<dim,spacedim> *fg = static_cast<struct find_ghosts<dim,spacedim> *>(user_data);
+ sc_array_t *subids = fg->subids;
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim> *triangulation = fg->triangulation;
+ int nsubs;
+ dealii::types::subdomain_id *subdomain_ids;
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors;
+
+ subids->elem_count = 0;
+ for (i = 0; i < nsides; i++) {
+ if (sides[i].is_ghost) {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad));
+ Assert (cell->is_ghost(), ExcMessage ("ghost quad did not find ghost cell"));
+ dealii::types::subdomain_id *subid =
+ static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
+ *subid = cell->subdomain_id();
+ }
+ }
+
+ if (!subids->elem_count) {
+ return;
+ }
+
+ nsubs = (int) subids->elem_count;
+ subdomain_ids = (dealii::types::subdomain_id *) (subids->array);
+
+ for (i = 0; i < nsides; i++) {
+ if (!sides[i].is_ghost) {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad));
+
+ Assert (!cell->is_ghost(), ExcMessage ("local quad found ghost cell"));
+
+ for (j = 0; j < nsubs; j++) {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(sides[i].corner)]
+ .insert (subdomain_ids[j]);
+ }
+ }
+ }
+
+ subids->elem_count = 0;
+ }
+
+ /** Similar to find_ghosts_corner, but for the hanging vertex in the
+ * middle of an edge
+ */
+ template <int dim, int spacedim>
+ void
+ find_ghosts_edge
+ (typename dealii::internal::p4est::iter<dim>::edge_info * info,
+ void *user_data)
+ {
+ int i, j, k;
+ int nsides = info->sides.elem_count;
+ typename dealii::internal::p4est::iter<dim>::edge_side * sides =
+ (typename dealii::internal::p4est::iter<dim>::edge_side *)
+ (info->sides.array);
+ struct find_ghosts<dim,spacedim> *fg = static_cast<struct find_ghosts<dim,spacedim> *>(user_data);
+ sc_array_t *subids = fg->subids;
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim> *triangulation = fg->triangulation;
+ int nsubs;
+ dealii::types::subdomain_id *subdomain_ids;
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors;
+
+ subids->elem_count = 0;
+ for (i = 0; i < nsides; i++) {
+ if (sides[i].is_hanging) {
+ for (j = 0; j < 2; j++) {
+ if (sides[i].is.hanging.is_ghost[j]) {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+ dealii::types::subdomain_id *subid =
+ static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
+ *subid = cell->subdomain_id();
+ }
+ }
+ }
+ }
+
+ if (!subids->elem_count) {
+ return;
+ }
+
+ nsubs = (int) subids->elem_count;
+ subdomain_ids = (dealii::types::subdomain_id *) (subids->array);
+
+ for (i = 0; i < nsides; i++) {
+ if (sides[i].is_hanging) {
+ for (j = 0; j < 2; j++) {
+ if (!sides[i].is.hanging.is_ghost[j]) {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+
+ for (k = 0; k < nsubs; k++) {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_edge_corners[sides[i].edge][1^j])]
+ .insert (subdomain_ids[k]);
+ }
+ }
+ }
+ }
+ }
+
+ subids->elem_count = 0;
+ }
+
+ /** Similar to find_ghosts_corner, but for the hanging vertex in the
+ * middle of a face
+ */
+ template <int dim, int spacedim>
+ void
+ find_ghosts_face
+ (typename dealii::internal::p4est::iter<dim>::face_info * info,
+ void *user_data)
+ {
+ int i, j, k;
+ int nsides = info->sides.elem_count;
+ typename dealii::internal::p4est::iter<dim>::face_side * sides =
+ (typename dealii::internal::p4est::iter<dim>::face_side *)
+ (info->sides.array);
+ struct find_ghosts<dim,spacedim> *fg = static_cast<struct find_ghosts<dim,spacedim> *>(user_data);
+ sc_array_t *subids = fg->subids;
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim> *triangulation = fg->triangulation;
+ int nsubs;
+ dealii::types::subdomain_id *subdomain_ids;
+ std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors;
+ int limit = (dim == 2) ? 2 : 4;
+
+ subids->elem_count = 0;
+ for (i = 0; i < nsides; i++) {
+ if (sides[i].is_hanging) {
+ for (j = 0; j < limit; j++) {
+ if (sides[i].is.hanging.is_ghost[j]) {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+ dealii::types::subdomain_id *subid =
+ static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
+ *subid = cell->subdomain_id();
+ }
+ }
+ }
+ }
+
+ if (!subids->elem_count) {
+ return;
+ }
+
+ nsubs = (int) subids->elem_count;
+ subdomain_ids = (dealii::types::subdomain_id *) (subids->array);
+
+ for (i = 0; i < nsides; i++) {
+ if (sides[i].is_hanging) {
+ for (j = 0; j < limit; j++) {
+ if (!sides[i].is.hanging.is_ghost[j]) {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+
+ for (k = 0; k < nsubs; k++) {
+ if (dim == 2) {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(p4est_face_corners[sides[i].face][(limit - 1)^j])]
+ .insert (subdomain_ids[k]);
+ }
+ else {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_face_corners[sides[i].face][(limit - 1)^j])]
+ .insert (subdomain_ids[k]);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ subids->elem_count = 0;
+ }
+ }
+
+ template <>
+ void
+ Triangulation<1,1>::
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+ template <>
+ void
+ Triangulation<1,2>::
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+ template <>
+ void
+ Triangulation<1,3>::
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+ /**
+ * Determine the neighboring subdomains that are adjacent to each vertex.
+ * This is achieved via the p4est_iterate tool
+ */
+ template <>
+ void
+ Triangulation<2,2>::
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
+ {
+ struct find_ghosts<2,2> fg;
+
+ fg.subids = sc_array_new (sizeof (dealii::types::subdomain_id));
+ fg.triangulation = this;
+ fg.vertices_with_ghost_neighbors = &(vertices_with_ghost_neighbors);
+
+ p4est_iterate (this->parallel_forest, this->parallel_ghost, static_cast<void *>(&fg),
+ NULL, find_ghosts_face<2,2>, find_ghosts_corner<2,2>);
+
+ sc_array_destroy (fg.subids);
+ }
+
+ /**
+ * Determine the neighboring subdomains that are adjacent to each vertex.
+ * This is achieved via the p4est_iterate tool
+ */
+ template <>
+ void
+ Triangulation<2,3>::
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
+ {
+ struct find_ghosts<2,3> fg;
+
+ fg.subids = sc_array_new (sizeof (dealii::types::subdomain_id));
+ fg.triangulation = this;
+ fg.vertices_with_ghost_neighbors = &(vertices_with_ghost_neighbors);
+
+ p4est_iterate (this->parallel_forest, this->parallel_ghost, static_cast<void *>(&fg),
+ NULL, find_ghosts_face<2,3>, find_ghosts_corner<2,3>);
+
+ sc_array_destroy (fg.subids);
+ }
+
+ /**
+ * Determine the neighboring subdomains that are adjacent to each vertex.
+ * This is achieved via the p8est_iterate tool
+ */
+ template <>
+ void
+ Triangulation<3,3>::
+ fill_vertices_with_ghost_neighbors
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
+ {
+ struct find_ghosts<3,3> fg;
+
+ fg.subids = sc_array_new (sizeof (dealii::types::subdomain_id));
+ fg.triangulation = this;
+ fg.vertices_with_ghost_neighbors = &(vertices_with_ghost_neighbors);
+
+ p8est_iterate (this->parallel_forest, this->parallel_ghost, static_cast<void *>(&fg),
+ NULL, find_ghosts_face<3,3>, find_ghosts_edge<3,3>, find_ghosts_corner<3,3>);
+
+ sc_array_destroy (fg.subids);
+ }
template <int dim, int spacedim>
MPI_Comm
}
+ template <int dim, int spacedim>
+ void
+ Triangulation<dim,spacedim>::add_periodicity
+ (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >&
+ periodicity_vector)
+ {
+#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+ Assert (triangulation_has_content == true,
+ ExcMessage ("The triangulation is empty!"));
+ Assert (this->n_levels() == 1,
+ ExcMessage ("The triangulation is refined!"));
+
+ typename std::vector
+ <typename std_cxx1x::tuple
+ <cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >::const_iterator periodic_tuple;
+ periodic_tuple = periodicity_vector.begin();
+
+ typename std::vector
+ <typename std_cxx1x::tuple
+ <cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >::const_iterator periodic_end;
+ periodic_end = periodicity_vector.end();
+
+ for (; periodic_tuple<periodic_end; ++periodic_tuple)
+ {
+ const cell_iterator first_cell=std_cxx1x::get<0>(*periodic_tuple);
+ const cell_iterator second_cell=std_cxx1x::get<2>(*periodic_tuple);
+ const unsigned int face_right=std_cxx1x::get<3>(*periodic_tuple);
+ const unsigned int face_left=std_cxx1x::get<1>(*periodic_tuple);
+
+ //respective cells of the matching faces in p4est
+ const unsigned int tree_left
+ = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
+ first_cell)];
+ const unsigned int tree_right
+ = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
+ second_cell)];
+
+ dealii::internal::p4est::functions<dim>::
+ connectivity_join_faces (connectivity,
+ tree_left,
+ tree_right,
+ face_left,
+ face_right,
+ /* orientation */ 0);
+ /* The orientation parameter above describes the difference between
+ * the cell on the left and the cell on the right would number of the
+ * corners of the face. In the periodic domains most users will want,
+ * this orientation will be 0, i.e. the two cells would number the
+ * corners the same way. More exotic periodicity, like moebius strips
+ * or converting an unstructured quad/hex mesh into a periodic domain,
+ * are not supported right now, and undefined behavior will occur if
+ * users try to make them periodic. This may be addressed at a later
+ * date.
+ */
+ }
+
+
+ Assert(dealii::internal::p4est::functions<dim>::connectivity_is_valid
+ (connectivity) == 1,
+ ExcInternalError());
+
+ // now create a forest out of the connectivity data structure
+ dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
+ parallel_forest
+ = dealii::internal::p4est::functions<dim>::
+ new_forest (mpi_communicator,
+ connectivity,
+ /* minimum initial number of quadrants per tree */ 0,
+ /* minimum level of upfront refinement */ 0,
+ /* use uniform upfront refinement */ 1,
+ /* user_data_size = */ 0,
+ /* user_data_constructor = */ NULL,
+ /* user_pointer */ this);
+
+#else
+ Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
+#endif
+ }
+
template <int dim, int spacedim>
std::size_t
}
catch (const typename dealii::Triangulation<dim,spacedim>::DistortedCellList &)
{
- // the underlying
- // triangulation should not
- // be checking for
- // distorted cells
+ // the underlying triangulation should not be checking for distorted
+ // cells
AssertThrow (false, ExcInternalError());
}
- // note that now we have some content in
- // the p4est objects and call the
- // functions that do the actual work
- // (which are dimension dependent, so
+ // note that now we have some content in the p4est objects and call the
+ // functions that do the actual work (which are dimension dependent, so
// separate)
triangulation_has_content = true;
}
catch (const typename Triangulation<dim>::DistortedCellList &)
{
- // the underlying
- // triangulation should not
- // be checking for
- // distorted cells
+ // the underlying triangulation should not be checking for distorted
+ // cells
AssertThrow (false, ExcInternalError());
}
Triangulation<dim,spacedim>::
attach_mesh_data()
{
- // determine size of memory in bytes to
- // attach to each cell. This needs to
- // be constant because of p4est.
+ // determine size of memory in bytes to attach to each cell. This needs
+ // to be constant because of p4est.
if (attached_data_size==0)
{
Assert(n_attached_datas==0, ExcInternalError());
}
}
+ template <int dim, int spacedim>
+ typename dealii::Triangulation<dim,spacedim>::cell_iterator
+ cell_from_quad
+ (typename dealii::parallel::distributed::Triangulation<dim,spacedim> *triangulation,
+ typename dealii::internal::p4est::types<dim>::topidx treeidx,
+ typename dealii::internal::p4est::types<dim>::quadrant &quad)
+ {
+ int i, l = quad.level;
+ int child_id;
+ const std::vector<types::global_dof_index> perm = triangulation->get_p4est_tree_to_coarse_cell_permutation ();
+ types::global_dof_index dealii_index = perm[treeidx];
+
+ for (i = 0; i < l; i++) {
+ typename dealii::Triangulation<dim,spacedim>::cell_iterator cell (triangulation, i, dealii_index);
+ child_id = internal::p4est::functions<dim>::quadrant_ancestor_id (&quad, i + 1);
+ Assert (cell->has_children (), ExcMessage ("p4est quadrant does not correspond to a cell!"));
+ dealii_index = cell->child_index(child_id);
+ }
+
+ typename dealii::Triangulation<dim,spacedim>::cell_iterator out_cell (triangulation, l, dealii_index);
+
+ return out_cell;
+ }
+
- // TODO: again problems with specialization in
- // only one template argument
+ // TODO: again problems with specialization in only one template argument
template <>
Triangulation<1,1>::Triangulation (MPI_Comm)
{
}
template <>
- const std::vector<unsigned int> &
+ const std::vector<types::global_dof_index> &
Triangulation<1,1>::get_p4est_tree_to_coarse_cell_permutation() const
{
- static std::vector<unsigned int> a;
+ static std::vector<types::global_dof_index> a;
return a;
}
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/dofs/*.h
- )
-
SET(_src
block_info.cc
dof_accessor.cc
SET(_inst
block_info.inst.in
- dof_accessor.inst.in
dof_accessor_get.inst.in
+ dof_accessor.inst.in
dof_accessor_set.inst.in
dof_handler.inst.in
dof_handler_policy.inst.in
dof_objects.inst.in
dof_renumbering.inst.in
- dof_tools.inst.in
dof_tools_constraints.inst.in
+ dof_tools.inst.in
dof_tools_sparsity.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/dofs/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_dofs OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_dofs "${_inst}")
-
}
}
+ //additionally (multigrid only), we can have the case that children of our neighbor
+ //have us as a neighbor. In this case we and the children are active.
+ if (dealii_cell->active())
+ {
+ for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+ {
+ if (dealii_cell->at_boundary(f))
+ continue;
+ typename DoFHandler<dim,spacedim>::level_cell_iterator neighbor = dealii_cell->neighbor(f);
+ if (!neighbor->has_children())
+ continue;
+
+ for (unsigned int subface=0; subface<GeometryInfo<dim>::max_children_per_face; ++subface)
+ {
+ typename DoFHandler<dim,spacedim>::level_cell_iterator child = dealii_cell->neighbor_child_on_subface(f,subface);
+ dealii::types::subdomain_id dest = child->subdomain_id();
+ Assert(dest != dealii::numbers::artificial_subdomain_id, ExcInternalError());
+ if (dest != tria.locally_owned_subdomain())
+ send_to.insert(dest);
+ }
+
+ }
+
+ }
+
+
+ // send if we have something to send
if (send_to.size() > 0)
{
// this cell's dof_indices
MPI_Get_count(&status, MPI_BYTE, &len);
receive.resize(len);
+#ifdef DEBUG
+ Assert(senders.find(status.MPI_SOURCE)!=senders.end(), ExcInternalError());
+#endif
+
char *ptr = &receive[0];
MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
tr->get_communicator(), &status);
unsigned int sent=needs_to_get_cells.size();
unsigned int recv=senders.size();
- MPI_Reduce(&sent, &sum_send, 1, MPI_UNSIGNED, MPI_SUM, 0, tr->get_communicator());
- MPI_Reduce(&recv, &sum_recv, 1, MPI_UNSIGNED, MPI_SUM, 0, tr->get_communicator());
+ MPI_Reduce(&sent, &sum_send, 1, MPI_UNSIGNED, MPI_SUM, 56345, tr->get_communicator());
+ MPI_Reduce(&recv, &sum_recv, 1, MPI_UNSIGNED, MPI_SUM, 56346, tr->get_communicator());
Assert(sum_send==sum_recv, ExcInternalError());
}
#endif
if (!cell->is_artificial())
cell->set_user_flag();
- //mark the vertices we are interested
- //in, i.e. belonging to own and marked cells
- const std::vector<bool> locally_active_vertices
- = mark_locally_active_vertices (*tr);
-
// add each ghostcells'
// subdomain to the vertex and
// keep track of interesting
// neighbors
std::map<unsigned int, std::set<dealii::types::subdomain_id> >
vertices_with_ghost_neighbors;
- for (typename DoFHandler<dim,spacedim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell)
- if (cell->is_ghost ())
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- if (locally_active_vertices[cell->vertex_index(v)])
- vertices_with_ghost_neighbors[cell->vertex_index(v)]
- .insert (cell->subdomain_id());
+
+ tr->fill_vertices_with_ghost_neighbors (vertices_with_ghost_neighbors);
/* Send and receive cells. After this,
for (typename DoFHandler<dim,spacedim>::level_cell_iterator
cell = dof_handler.begin(level);
cell != dof_handler.end(level); ++cell)
- if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id && cell->level_subdomain_id() != tr->locally_owned_subdomain())
+ if (cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id
+ && cell->level_subdomain_id() != tr->locally_owned_subdomain())
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
if (locally_active_vertices[cell->vertex_index(v)])
vertices_with_ghost_neighbors[cell->vertex_index(v)]
local_dof_indices.end(),
DoFHandler<dim,spacedim>::invalid_dof_index))
{
- Assert(false, ExcMessage ("not all dofs got distributed!"));
+ Assert(false, ExcMessage ("not all DoFs got distributed!"));
}
}
}
for (cell = dof_handler.begin_active(); cell != endc; ++cell)
if (!cell->is_artificial())
cell->set_user_flag();
- //mark the vertices we are interested
- //in, i.e. belonging to own and marked cells
- const std::vector<bool> locally_active_vertices
- = mark_locally_active_vertices (*tr);
// add each ghostcells'
// subdomain to the vertex and
// neighbors
std::map<unsigned int, std::set<dealii::types::subdomain_id> >
vertices_with_ghost_neighbors;
- for (typename DoFHandler<dim,spacedim>::active_cell_iterator
- cell = dof_handler.begin_active();
- cell != dof_handler.end(); ++cell)
- if (cell->is_ghost ())
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- if (locally_active_vertices[cell->vertex_index(v)])
- vertices_with_ghost_neighbors[cell->vertex_index(v)]
- .insert (cell->subdomain_id());
+
+ tr->fill_vertices_with_ghost_neighbors (vertices_with_ghost_neighbors);
// Send and receive cells. After this, only
// the local cells are marked, that received
(const DoFHandler<deal_II_dimension> & dof_handler,
IndexSet & dof_set);
+template
+void
+DoFTools::extract_locally_owned_dofs<hp::DoFHandler<deal_II_dimension> >
+(const hp::DoFHandler<deal_II_dimension> & dof_handler,
+ IndexSet & dof_set);
+
template
void
DoFTools::extract_locally_active_dofs<DoFHandler<deal_II_dimension> >
(const DoFHandler<deal_II_dimension> & dof_handler,
IndexSet & dof_set);
+template
+void
+DoFTools::extract_locally_active_dofs<hp::DoFHandler<deal_II_dimension> >
+(const hp::DoFHandler<deal_II_dimension> & dof_handler,
+ IndexSet & dof_set);
template
void
DoFTools::extract_locally_relevant_dofs<DoFHandler<deal_II_dimension> >
(const DoFHandler<deal_II_dimension> & dof_handler,
IndexSet & dof_set);
+
+template
+void
+DoFTools::extract_locally_relevant_dofs<hp::DoFHandler<deal_II_dimension> >
+(const hp::DoFHandler<deal_II_dimension> & dof_handler,
+ IndexSet & dof_set);
template
void
face_1->get_dof_indices(dofs_1, face_1_index);
face_2->get_dof_indices(dofs_2, face_2_index);
+ for (unsigned int i=0; i < dofs_per_face; i++) {
+ if (dofs_1[i] == numbers::invalid_dof_index ||
+ dofs_2[i] == numbers::invalid_dof_index) {
+ /* If either of these faces have no indices, stop. This is so
+ * that there is no attempt to match artificial cells of
+ * parallel distributed triangulations.
+ *
+ * While it seems like we ought to be able to avoid even calling
+ * set_periodicity_constraints for artificial faces, this
+ * situation can arise when a face that is being made periodic
+ * is only partially touched by the local subdomain.
+ * make_periodicity_constraints will be called recursively even
+ * for the section of the face that is not touched by the local
+ * subdomain.
+ *
+ * Until there is a better way to determine if the cells that
+ * neighbor a face are artificial, we simply test to see if the
+ * face does not have a valid dof initialization.
+ */
+ return;
+ }
+ }
+
// Well, this is a hack:
//
// There is no
Assert (0<=direction && direction<space_dim,
ExcIndexRange (direction, 0, space_dim));
-#if defined(DEBUG) && defined(DEAL_II_WITH_P4EST)
- // Check whether we run on a non parallel mesh or on a
- // parallel::distributed::Triangulation in serial
- {
- typedef parallel::distributed::Triangulation<DH::dimension,DH::space_dimension> PTRIA;
- const PTRIA *ptria_p = dynamic_cast<const PTRIA *> (&dof_handler.get_tria());
- Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1),
- ExcMessage ("This function can not be used with distributed triangulations."
- "See the documentation for more information."));
- }
-#endif
-
Assert (b_id1 != b_id2,
ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
"different to denote different boundaries."));
Assert(dim == space_dim,
ExcNotImplemented());
-#if defined(DEBUG) && defined(DEAL_II_WITH_P4EST)
- // Check whether we run on a non parallel mesh or on a
- // parallel::distributed::Triangulation in serial
- {
- typedef typename parallel::distributed::Triangulation<DH::dimension,DH::space_dimension> PTRIA;
- const PTRIA *ptria_p = dynamic_cast<const PTRIA *> (&dof_handler.get_tria());
- Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1),
- ExcMessage ("This function can not be used with distributed triangulations."
- "See the documentation for more information."));
- }
-#endif
-
typedef typename DH::face_iterator FaceIterator;
typedef std::map<FaceIterator, FaceIterator> FaceMap;
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/fe/*.h
- )
-
SET(_src
block_mask.cc
component_mask.cc
fe_nothing.cc
fe_poly.cc
fe_poly_tensor.cc
- fe_q.cc
fe_q_base.cc
+ fe_q.cc
fe_q_dg0.cc
fe_q_hierarchical.cc
fe_q_iso_q1.cc
fe_nothing.inst.in
fe_poly.inst.in
fe_poly_tensor.inst.in
- fe_q_base.inst.in
+ fe_q_base.inst.in
fe_q_dg0.inst.in
fe_q_hierarchical.inst.in
fe_q.inst.in
mapping_q.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/fe/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_fe OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_fe "${_inst}")
-
// this
for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
- lohvs[unit_point].end(),0.) - 1)<1e-12*this->degree*this->degree,
+ lohvs[unit_point].end(),0.) - 1)<1e-13*this->degree,
ExcInternalError());
}
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/grid/*.h
- )
-
SET(_src
grid_generator.cc
grid_in.cc
grid_tools.inst.in
intergrid_map.inst.in
manifold.inst.in
- manifold_lib.inst.in
tria_accessor.inst.in
tria_boundary.inst.in
tria_boundary_lib.inst.in
tria_objects.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/grid/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_grid OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_grid "${_inst}")
-
default:
Assert (false, ExcNotImplemented());
return "";
- };
+ }
}
for (unsigned int i=0; i<vertices.size(); ++i)
if (vertex_used[i])
- {
out << '\t' << vertices[i] << '\n';
- };
// write cells or faces
const bool write_cells = dx_flags.write_cells;
for (unsigned int d=spacedim+1; d<=3; ++d)
out << " 0"; // fill with zeroes
out << std::endl;
- };
+ }
// Write cells preamble
out << "$ENDNOD" << std::endl
++vertex)
out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
out << std::endl;
- };
+ }
// write faces and lines with
// non-zero boundary indicator
<< "# For a description of the UCD format see the AVS Developer's guide."
<< '\n'
<< "#" << '\n';
- };
+ }
// start with ucd data
out << n_vertices << ' '
for (unsigned int d=spacedim+1; d<=3; ++d)
out << " 0"; // fill with zeroes
out << '\n';
- };
+ }
// write cells. Enumerate cells
// consecutively, starting with 1
break;
default:
Assert (false, ExcNotImplemented());
- };
+ }
// it follows a list of the
// vertices of each cell. in 1d
++vertex)
out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex])+1 << ' ';
out << '\n';
- };
+ }
// write faces and lines with
// non-zero boundary indicator
break;
default:
Assert (false, ExcNotImplemented ());
- };
+ }
out << "\n";
// (iii) write vertex ordering
break;
default:
Assert (false, ExcNotImplemented ());
- };
+ }
// (iv) write a list of vertices of cells
out << "\n#"
break;
default:
Assert (false, ExcNotImplemented());
- };
+ }
out << static_cast<unsigned int>(face->boundary_indicator())
<< ' '
<< static_cast<unsigned int>(face->boundary_indicator())
out << '\n';
++index;
- };
+ }
}
break;
default:
Assert (false, ExcNotImplemented());
- };
+ }
// note: vertex numbers are 1-base
for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
out << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1 << ' ';
out << '\n';
++index;
- };
+ }
}
{
Assert(false, ExcInternalError());
break;
- };
+ }
case 2:
{
line_list.push_back (LineEntry(p0, p1,
face->user_flag_set(),
cell->level()));
- };
- };
- };
+ }
+ }
+ }
break;
- };
+ }
case 3:
{
y_max = std::max (y_max, line->second(1));
max_level = std::max (max_level, line->level);
- };
+ }
// scale in x-direction such that
// in the output 0 <= x <= 300.
"/MCshow { currentpoint stroke m\n"
"exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
<< '\n';
- };
+ }
out << "%%EndProlog" << '\n'
<< '\n';
// set fine lines
out << eps_flags_base.line_width << " setlinewidth" << '\n';
- };
+ }
// now write the lines
const Point<2> offset(x_min, y_min);
out << ")] "
<< "] -6 MCshow"
<< '\n';
- };
- };
+ }
+ }
// and the vertex numbers
if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
<< ")] "
<< "] -6 MCshow"
<< '\n';
- };
- };
+ }
+ }
out << "showpage" << '\n';
#include <deal.II/base/std_cxx1x/array.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary.h>
}
+ template<typename DH>
+ void
+ identify_periodic_face_pairs
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+ typename DH::cell_iterator, unsigned int> >
+ &periodicity_vector)
+ {
+ static const unsigned int dim = DH::dimension;
+ static const unsigned int spacedim = DH::space_dimension;
+
+ Assert (0<=direction && direction<spacedim,
+ ExcIndexRange(direction, 0, spacedim));
+
+ Assert (b_id1 != b_id2,
+ ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
+ "different to denote different boundaries."));
+
+ std::map<std::pair<typename DH::cell_iterator,
+ unsigned int>,
+ Point<spacedim> > face_locations;
+
+ // Collect faces with boundary_indicator b_id1
+ typename DH::cell_iterator cell = dof_handler.begin();
+ typename DH::cell_iterator endc = dof_handler.end();
+ for (; cell != endc; ++cell)
+ for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+ if(cell->face(f)->boundary_indicator() == b_id1)
+ {
+ Point<spacedim> face_center (cell->face(f)->center());
+ face_center(direction)=0.;
+ const std::pair<typename DH::cell_iterator, unsigned int>
+ cell_face_pair = std::make_pair(cell, f);
+ face_locations[cell_face_pair]=face_center;
+ }
+
+ // Match faces with boundary_indicator b_id2 to the ones in
+ //face_locations
+ cell = dof_handler.begin();
+ for (; cell != endc; ++cell)
+ for(unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+ if(cell->face(f)->boundary_indicator() == b_id2)
+ {
+ typename std::map<std::pair<typename DH::cell_iterator,
+ unsigned int>,
+ Point<spacedim> >::iterator p
+ = face_locations.begin();
+
+ Point<spacedim> center2 (cell->face(f)->center());
+ center2(direction)=0.;
+
+ for (; p != face_locations.end(); ++p)
+ {
+ if (center2.distance(p->second) < 1e-4*cell->face(f)->diameter())
+ {
+ const std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+ typename DH::cell_iterator, unsigned int>
+ periodic_tuple (p->first.first,
+ p->first.second,
+ cell,
+ f);
+
+ periodicity_vector.push_back(periodic_tuple);
+
+ face_locations.erase(p);
+ break;
+ }
+ Assert (p != face_locations.end(),
+ ExcMessage ("No corresponding face was found!"));
+ }
+ }
+
+ Assert (face_locations.size() == 0,
+ ExcMessage ("There are unmatched faces!"));
+ }
+
} /* namespace GridTools */
int,
const Tensor<1,deal_II_space_dimension> &);
+ template
+ void
+ identify_periodic_face_pairs
+ (const X &,
+ const types::boundary_id,
+ const types::boundary_id,
+ const unsigned int,
+ std::vector<std_cxx1x::tuple<X::cell_iterator, unsigned int,
+ X::cell_iterator, unsigned int> > &);
+
#endif
\}
#endif
}
+for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+ #if deal_II_dimension >= 2
+ namespace GridTools \{
+ template
+ void
+ identify_periodic_face_pairs
+ (const parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+ const types::boundary_id,
+ const types::boundary_id,
+ const unsigned int,
+ std::vector<std_cxx1x::tuple<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator,
+ unsigned int,
+ parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator,
+ unsigned int> > &);
+ \}
+ #endif
+#endif
+}
+
+
FlatManifold<spacedim>::FlatManifold ()
{}
-
template <int spacedim>
Point<spacedim>
FlatManifold<spacedim>::
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/hp/*.h
- )
-
SET(_src
dof_faces.cc
dof_handler.cc
- dof_levels.cc
- dof_objects.cc
+ dof_level.cc
fe_collection.cc
fe_values.cc
mapping_collection.cc
mapping_collection.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/hp/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_hp OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_hp "${_inst}")
-
{
namespace hp
{
+// ---------------------- DoFObjects ----------------------------
+ template <int structdim>
std::size_t
- DoFFaces<1>::memory_consumption () const
+ DoFIndicesOnFacesOrEdges<structdim>::memory_consumption () const
+ {
+ return (MemoryConsumption::memory_consumption (dofs) +
+ MemoryConsumption::memory_consumption (dof_offsets));
+ }
+
+
+ // explicit instantiations
+ template
+ std::size_t
+ DoFIndicesOnFacesOrEdges<1>::memory_consumption () const;
+
+ template
+ std::size_t
+ DoFIndicesOnFacesOrEdges<2>::memory_consumption () const;
+
+ template
+ std::size_t
+ DoFIndicesOnFacesOrEdges<3>::memory_consumption () const;
+
+
+// ---------------------- DoFFaces ----------------------------
+
+ std::size_t
+ DoFIndicesOnFaces<1>::memory_consumption () const
{
return 0;
}
std::size_t
- DoFFaces<2>::memory_consumption () const
+ DoFIndicesOnFaces<2>::memory_consumption () const
{
return MemoryConsumption::memory_consumption (lines);
}
std::size_t
- DoFFaces<3>::memory_consumption () const
+ DoFIndicesOnFaces<3>::memory_consumption () const
{
- return (MemoryConsumption::memory_consumption (quads) +
+ return (MemoryConsumption::memory_consumption (lines) +
MemoryConsumption::memory_consumption (quads) );
}
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/std_cxx1x/bind.h>
#include <deal.II/hp/dof_handler.h>
-#include <deal.II/hp/dof_objects.h>
-#include <deal.II/hp/dof_levels.h>
+#include <deal.II/hp/dof_level.h>
#include <deal.II/hp/dof_faces.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/grid/tria_accessor.h>
// active_fe_indices field which
// we have to backup before
{
- std::vector<std::vector<unsigned int> >
+ std::vector<std::vector<DoFLevel::active_fe_index_type> >
active_fe_backup(dof_handler.levels.size ());
for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
std::swap (dof_handler.levels[level]->active_fe_indices,
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
- dof_handler.levels.push_back (new internal::hp::DoFLevel<dim>);
+ dof_handler.levels.push_back (new internal::hp::DoFLevel);
std::swap (active_fe_backup[level],
dof_handler.levels[level]->active_fe_indices);
}
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
dof_handler.levels[level]->dof_offsets
- = std::vector<std::vector<types::global_dof_index>::size_type> (
+ = std::vector<DoFLevel::offset_type> (
dof_handler.tria->n_raw_lines(level),
- (std::vector<types::global_dof_index>::size_type)(-1));
+ (DoFLevel::offset_type)(-1));
+ dof_handler.levels[level]->cell_cache_offsets
+ = std::vector<DoFLevel::offset_type> (
+ dof_handler.tria->n_raw_lines(level),
+ (DoFLevel::offset_type)(-1));
types::global_dof_index next_free_dof = 0;
+ types::global_dof_index cache_size = 0;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell=dof_handler.begin_active(level);
cell!=dof_handler.end_active(level); ++cell)
{
dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
next_free_dof += cell->get_fe().dofs_per_line;
+
+ dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
+ cache_size += cell->get_fe().dofs_per_cell;
}
- dof_handler.levels[level]->dofs
+ dof_handler.levels[level]->dof_indices
= std::vector<types::global_dof_index> (next_free_dof,
DoFHandler<dim,spacedim>::invalid_dof_index);
+ dof_handler.levels[level]->cell_dof_indices_cache
+ = std::vector<types::global_dof_index> (cache_size,
+ DoFHandler<dim,spacedim>::invalid_dof_index);
}
// safety check: make sure that
if (!cell->has_children())
counter += cell->get_fe().dofs_per_line;
- Assert (dof_handler.levels[level]->dofs.size() == counter,
+ Assert (dof_handler.levels[level]->dof_indices.size() == counter,
ExcInternalError());
Assert (static_cast<unsigned int>
(std::count (dof_handler.levels[level]->dof_offsets.begin(),
dof_handler.levels[level]->dof_offsets.end(),
- (std::vector<types::global_dof_index>::size_type)(-1)))
+ (DoFLevel::offset_type)(-1)))
==
dof_handler.tria->n_raw_lines(level) - dof_handler.tria->n_active_lines(level),
ExcInternalError());
// active_fe_indices field which
// we have to backup before
{
- std::vector<std::vector<unsigned int> >
+ std::vector<std::vector<DoFLevel::active_fe_index_type> >
active_fe_backup(dof_handler.levels.size ());
for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
std::swap (dof_handler.levels[level]->active_fe_indices,
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
- dof_handler.levels.push_back (new internal::hp::DoFLevel<dim>);
+ dof_handler.levels.push_back (new internal::hp::DoFLevel);
std::swap (active_fe_backup[level],
dof_handler.levels[level]->active_fe_indices);
}
- dof_handler.faces = new internal::hp::DoFFaces<2>;
+ dof_handler.faces = new internal::hp::DoFIndicesOnFaces<2>;
}
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
dof_handler.levels[level]->dof_offsets
- = std::vector<std::vector<types::global_dof_index>::size_type> (
+ = std::vector<DoFLevel::offset_type> (
dof_handler.tria->n_raw_quads(level),
- (std::vector<types::global_dof_index>::size_type)(-1));
+ (DoFLevel::offset_type)(-1));
+ dof_handler.levels[level]->cell_cache_offsets
+ = std::vector<DoFLevel::offset_type> (
+ dof_handler.tria->n_raw_quads(level),
+ (DoFLevel::offset_type)(-1));
types::global_dof_index next_free_dof = 0;
+ types::global_dof_index cache_size = 0;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell=dof_handler.begin_active(level);
cell!=dof_handler.end_active(level); ++cell)
{
dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
next_free_dof += cell->get_fe().dofs_per_quad;
+
+ dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
+ cache_size += cell->get_fe().dofs_per_cell;
}
- dof_handler.levels[level]->dofs
+ dof_handler.levels[level]->dof_indices
= std::vector<types::global_dof_index> (next_free_dof,
DoFHandler<dim,spacedim>::invalid_dof_index);
+ dof_handler.levels[level]->cell_dof_indices_cache
+ = std::vector<types::global_dof_index> (cache_size,
+ DoFHandler<dim,spacedim>::invalid_dof_index);
}
// safety check: make sure that
if (!cell->has_children())
counter += cell->get_fe().dofs_per_quad;
- Assert (dof_handler.levels[level]->dofs.size() == counter,
+ Assert (dof_handler.levels[level]->dof_indices.size() == counter,
ExcInternalError());
Assert (static_cast<unsigned int>
(std::count (dof_handler.levels[level]->dof_offsets.begin(),
dof_handler.levels[level]->dof_offsets.end(),
- (std::vector<types::global_dof_index>::size_type)(-1)))
+ (DoFLevel::offset_type)(-1)))
==
dof_handler.tria->n_raw_quads(level) - dof_handler.tria->n_active_quads(level),
ExcInternalError());
// then allocate as much space as
// we need and prime the linked
// list for lines (see the
- // description in hp::DoFLevels)
+ // description in hp::DoFLevel)
// with the indices we will
// need. note that our task is
// more complicated since two
// active ones will have a
// non-invalid value later on
dof_handler.faces->lines.dof_offsets
- = std::vector<std::vector<types::global_dof_index>::size_type>
- (dof_handler.tria->n_raw_lines(),
- (std::vector<types::global_dof_index>::size_type)(-1));
+ = std::vector<unsigned int> (dof_handler.tria->n_raw_lines(),
+ (unsigned int)(-1));
dof_handler.faces->lines.dofs
= std::vector<types::global_dof_index> (n_line_slots,
DoFHandler<dim,spacedim>::invalid_dof_index);
// active_fe_indices field which
// we have to backup before
{
- std::vector<std::vector<unsigned int> >
+ std::vector<std::vector<DoFLevel::active_fe_index_type> >
active_fe_backup(dof_handler.levels.size ());
for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
std::swap (dof_handler.levels[level]->active_fe_indices,
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
- dof_handler.levels.push_back (new internal::hp::DoFLevel<dim>);
+ dof_handler.levels.push_back (new internal::hp::DoFLevel);
std::swap (active_fe_backup[level],
dof_handler.levels[level]->active_fe_indices);
}
- dof_handler.faces = new internal::hp::DoFFaces<3>;
+ dof_handler.faces = new internal::hp::DoFIndicesOnFaces<3>;
}
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
dof_handler.levels[level]->dof_offsets
- = std::vector<std::vector<types::global_dof_index>::size_type>
- (dof_handler.tria->n_raw_hexs(level),
- (std::vector<types::global_dof_index>::size_type)(-1));
+ = std::vector<DoFLevel::offset_type> (
+ dof_handler.tria->n_raw_hexs(level),
+ (DoFLevel::offset_type)(-1));
+ dof_handler.levels[level]->cell_cache_offsets
+ = std::vector<DoFLevel::offset_type> (
+ dof_handler.tria->n_raw_hexs(level),
+ (DoFLevel::offset_type)(-1));
types::global_dof_index next_free_dof = 0;
+ types::global_dof_index cache_size = 0;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell=dof_handler.begin_active(level);
cell!=dof_handler.end_active(level); ++cell)
{
dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
next_free_dof += cell->get_fe().dofs_per_hex;
+
+ dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
+ cache_size += cell->get_fe().dofs_per_cell;
}
- dof_handler.levels[level]->dofs
+ dof_handler.levels[level]->dof_indices
= std::vector<types::global_dof_index> (next_free_dof,
DoFHandler<dim,spacedim>::invalid_dof_index);
+ dof_handler.levels[level]->cell_dof_indices_cache
+ = std::vector<types::global_dof_index> (cache_size,
+ DoFHandler<dim,spacedim>::invalid_dof_index);
}
// safety check: make sure that
if (!cell->has_children())
counter += cell->get_fe().dofs_per_hex;
- Assert (dof_handler.levels[level]->dofs.size() == counter,
+ Assert (dof_handler.levels[level]->dof_indices.size() == counter,
ExcInternalError());
Assert (static_cast<unsigned int>
(std::count (dof_handler.levels[level]->dof_offsets.begin(),
dof_handler.levels[level]->dof_offsets.end(),
- (std::vector<types::global_dof_index>::size_type)(-1)))
+ (DoFLevel::offset_type)(-1)))
==
dof_handler.tria->n_raw_hexs(level) - dof_handler.tria->n_active_hexs(level),
ExcInternalError());
// then allocate as much space as
// we need and prime the linked
// list for quad (see the
- // description in hp::DoFLevels)
+ // description in hp::DoFLevel)
// with the indices we will
// need. note that our task is
// more complicated since two
if (true)
{
dof_handler.faces->quads.dof_offsets
- = std::vector<std::vector<types::global_dof_index>::size_type>
+ = std::vector<unsigned int>
(dof_handler.tria->n_raw_quads(),
- (std::vector<types::global_dof_index>::size_type)(-1));
+ (unsigned int)(-1));
dof_handler.faces->quads.dofs
= std::vector<types::global_dof_index> (n_quad_slots,
DoFHandler<dim,spacedim>::invalid_dof_index);
= std::vector<IndexSet> (1,
number_cache.locally_owned_dofs);
+ // update the cache used for cell dof indices and compress the data on the levels
+ for (active_cell_iterator cell = begin_active();
+ cell != end(); ++cell)
+ cell->update_cell_dof_indices_cache ();
+
+ for (unsigned int level=0; level<levels.size(); ++level)
+ levels[level]->compress_data (*finite_elements);
+
// finally restore the user flags
const_cast<Triangulation<dim,spacedim> &>(*tria).load_user_flags(user_flags);
}
}
#endif
+ // uncompress the internal storage scheme of dofs on cells
+ // so that we can access dofs in turns
+ for (unsigned int level=0; level<levels.size(); ++level)
+ levels[level]->uncompress_data (*finite_elements);
+
+ // do the renumbering
renumber_dofs_internal (new_numbers, dealii::internal::int2type<dim>());
+
+ // update the cache used for cell dof indices
+ for (active_cell_iterator cell = begin_active();
+ cell != end(); ++cell)
+ cell->update_cell_dof_indices_cache ();
+
+ // now re-compress the dof indices
+ for (unsigned int level=0; level<levels.size(); ++level)
+ levels[level]->compress_data (*finite_elements);
}
// Create sufficiently many
// hp::DoFLevels.
while (levels.size () < tria->n_levels ())
- levels.push_back (new dealii::internal::hp::DoFLevel<dim>);
+ levels.push_back (new dealii::internal::hp::DoFLevel);
// then make sure that on each
// level we have the appropriate
{
Assert (has_children.size () == levels.size (), ExcInternalError ());
- // In each refinement at most one new
- // new level may appear. If that happened
- // it is appended to the DoFHandler
- // levels.
- if (levels.size () < tria->n_levels ())
- levels.push_back (new dealii::internal::hp::DoFLevel<dim>);
+ // Normally only one level is added, but if this Triangulation
+ // is created by copy_triangulation, it can be more than one level.
+ while (levels.size () < tria->n_levels ())
+ levels.push_back (new dealii::internal::hp::DoFLevel);
// Coarsening can lead to the loss
// of levels. Hence remove them.
levels.pop_back ();
}
+ Assert(levels.size () == tria->n_levels (), ExcInternalError());
+
// Resize active_fe_indices
// vectors. use zero indicator to
// extend
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2003 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/hp/dof_level.h>
+#include <deal.II/hp/fe_collection.h>
+#include <iostream>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+ namespace hp
+ {
+ template <int dim, int spacedim>
+ void
+ DoFLevel::compress_data (const dealii::hp::FECollection<dim,spacedim> &fe_collection)
+ {
+ return;
+
+ if (dof_offsets.size() == 0 || dof_indices.size()==0)
+ return;
+
+ // in a first run through, count how many new slots we need in the
+ // dof_indices array after compression. note that the 'cell'
+ // counter is incremented inside the loop
+ unsigned int new_size = 0;
+ for (unsigned int cell=0; cell<dof_offsets.size(); )
+ // see if this cell is active on the current level
+ if (dof_offsets[cell] != (offset_type)(-1))
+ {
+ // find the next cell active on this level
+ unsigned int next_cell = cell+1;
+ while ((next_cell<dof_offsets.size()) &&
+ (dof_offsets[next_cell] == (offset_type)(-1)))
+ ++next_cell;
+
+ const unsigned int next_offset = (next_cell < dof_offsets.size() ?
+ dof_offsets[next_cell] :
+ dof_indices.size());
+
+ Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
+ ExcInternalError());
+
+ // see if the range of dofs for this cell can be compressed and if so
+ // how many slots we have to store for them
+ if (next_offset > dof_offsets[cell])
+ {
+ bool compressible = true;
+ for (unsigned int j=dof_offsets[cell]+1; j<next_offset; ++j)
+ if (dof_indices[j] != dof_indices[j-1]+1)
+ {
+ compressible = false;
+ break;
+ }
+ if (compressible == true)
+ new_size += 1;
+ else
+ new_size += (next_offset-dof_offsets[cell]);
+ }
+
+ // then move on to the next cell
+ cell = next_cell;
+ }
+ else
+ ++cell;
+
+ // now allocate the new array and copy into it whatever we need
+ std::vector<types::global_dof_index> new_dof_indices;
+ new_dof_indices.reserve(new_size);
+ for (unsigned int cell=0; cell<dof_offsets.size(); )
+ // see if this cell is active on the current level
+ if (dof_offsets[cell] != (offset_type)(-1))
+ {
+ // find the next cell active on this level
+ unsigned int next_cell = cell+1;
+ while ((next_cell<dof_offsets.size()) &&
+ (dof_offsets[next_cell] == (offset_type)(-1)))
+ ++next_cell;
+
+ const unsigned int next_offset = (next_cell < dof_offsets.size() ?
+ dof_offsets[next_cell] :
+ dof_indices.size());
+
+ Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
+ ExcInternalError());
+
+ // see if the range of dofs for this cell can be compressed and if so
+ // how many slots we have to store for them
+ if (next_offset > dof_offsets[cell])
+ {
+ bool compressible = true;
+ for (unsigned int j=dof_offsets[cell]+1; j<next_offset; ++j)
+ if (dof_indices[j] != dof_indices[j-1]+1)
+ {
+ compressible = false;
+ break;
+ }
+
+ // if this cell is compressible, then copy the first index and mark this
+ // in the dof_offsets array
+ if (compressible == true)
+ {
+ new_dof_indices.push_back (dof_indices[dof_offsets[cell]]);
+
+ // make sure that the current active_fe_index indicates
+ // that this entry hasn't been compressed yet
+ Assert ((signed_active_fe_index_type)active_fe_indices[cell] >= 0, ExcInternalError());
+
+ // then mark the compression
+ active_fe_indices[cell] = (active_fe_index_type)~(signed_active_fe_index_type)active_fe_indices[cell];
+ }
+ else
+ for (unsigned int i=dof_offsets[cell]; i<next_offset; ++i)
+ new_dof_indices.push_back (dof_indices[i]);
+ }
+
+ // then move on to the next cell
+ cell = next_cell;
+ }
+ else
+ ++cell;
+
+ // finally swap old and new content
+ Assert (new_dof_indices.size() == new_size, ExcInternalError());
+ dof_indices.swap (new_dof_indices);
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ DoFLevel::uncompress_data(const dealii::hp::FECollection<dim,spacedim> &fe_collection)
+ {
+ return;
+
+ if (dof_offsets.size() == 0 || dof_indices.size()==0)
+ return;
+
+ // in a first run through, count how many new slots we need in the
+ // dof_indices array after uncompression.
+ unsigned int new_size = 0;
+ for (unsigned int cell=0; cell<dof_offsets.size(); ++cell)
+ if (dof_offsets[cell] != (offset_type)(-1))
+ {
+ // we know now that the slot for this cell is used. extract the
+ // active_fe_index for it and see how many entries we need
+ new_size += fe_collection[active_fe_index(cell)].template n_dofs_per_object<dim>();
+ }
+
+ // now allocate the new array and copy into it whatever we need
+ std::vector<types::global_dof_index> new_dof_indices;
+ new_dof_indices.reserve(new_size);
+ std::vector<offset_type> new_dof_offsets (dof_offsets.size(), (offset_type)(-1));
+ for (unsigned int cell=0; cell<dof_offsets.size(); )
+ // see if this cell is active on the current level
+ if (dof_offsets[cell] != (offset_type)(-1))
+ {
+ // find the next cell active on this level
+ unsigned int next_cell = cell+1;
+ while ((next_cell<dof_offsets.size()) &&
+ (dof_offsets[next_cell] == (offset_type)(-1)))
+ ++next_cell;
+
+ const unsigned int next_offset = (next_cell < dof_offsets.size() ?
+ dof_offsets[next_cell] :
+ dof_indices.size());
+
+ // set offset for this cell
+ new_dof_offsets[cell] = new_dof_indices.size();
+
+ // see if we need to uncompress this set of dofs
+ if ((signed_active_fe_index_type)active_fe_indices[cell]>=0)
+ {
+ // apparently not. simply copy them
+ Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
+ ExcInternalError());
+ for (unsigned int i=dof_offsets[cell]; i<next_offset; ++i)
+ new_dof_indices.push_back (dof_indices[i]);
+ }
+ else
+ {
+ // apparently so. uncompress
+ Assert (next_offset-dof_offsets[cell] == 1,
+ ExcInternalError());
+ for (unsigned int i=0; i<fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(); ++i)
+ new_dof_indices.push_back (dof_indices[dof_offsets[cell]]+i);
+ }
+
+ // then move on to the next cell
+ cell = next_cell;
+ }
+ else
+ ++cell;
+
+ // verify correct size, then swap arrays
+ Assert (new_dof_indices.size() == new_size, ExcInternalError());
+ dof_indices.swap (new_dof_indices);
+ dof_offsets.swap (new_dof_offsets);
+ }
+
+
+ std::size_t
+ DoFLevel::memory_consumption () const
+ {
+ return (MemoryConsumption::memory_consumption (active_fe_indices) +
+ MemoryConsumption::memory_consumption (dof_indices) +
+ MemoryConsumption::memory_consumption (dof_offsets) +
+ MemoryConsumption::memory_consumption (cell_cache_offsets) +
+ MemoryConsumption::memory_consumption(cell_dof_indices_cache));
+ }
+
+
+ // explicit instantiations
+ template void DoFLevel::compress_data(const dealii::hp::FECollection<1,1> &);
+ template void DoFLevel::compress_data(const dealii::hp::FECollection<1,2> &);
+ template void DoFLevel::compress_data(const dealii::hp::FECollection<1,3> &);
+ template void DoFLevel::compress_data(const dealii::hp::FECollection<2,2> &);
+ template void DoFLevel::compress_data(const dealii::hp::FECollection<2,3> &);
+ template void DoFLevel::compress_data(const dealii::hp::FECollection<3,3> &);
+
+ template void DoFLevel::uncompress_data(const dealii::hp::FECollection<1,1> &);
+ template void DoFLevel::uncompress_data(const dealii::hp::FECollection<1,2> &);
+ template void DoFLevel::uncompress_data(const dealii::hp::FECollection<1,3> &);
+ template void DoFLevel::uncompress_data(const dealii::hp::FECollection<2,2> &);
+ template void DoFLevel::uncompress_data(const dealii::hp::FECollection<2,3> &);
+ template void DoFLevel::uncompress_data(const dealii::hp::FECollection<3,3> &);
+ }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+++ /dev/null
-// ---------------------------------------------------------------------
-// $Id$
-//
-// Copyright (C) 2003 - 2013 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/hp/dof_levels.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace internal
-{
- namespace hp
- {
- template <int dim>
- std::size_t
- DoFLevel<dim>::memory_consumption () const
- {
- return (MemoryConsumption::memory_consumption (active_fe_indices) +
- MemoryConsumption::memory_consumption (dofs) +
- MemoryConsumption::memory_consumption (dof_offsets));
- }
-
-
- // explicit instantiations
- template std::size_t DoFLevel<1>::memory_consumption () const;
- template std::size_t DoFLevel<2>::memory_consumption () const;
- template std::size_t DoFLevel<3>::memory_consumption () const;
- }
-}
-
-DEAL_II_NAMESPACE_CLOSE
+++ /dev/null
-// ---------------------------------------------------------------------
-// $Id$
-//
-// Copyright (C) 2006 - 2013 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/hp/dof_objects.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace internal
-{
- namespace hp
- {
- template <int structdim>
- std::size_t
- DoFObjects<structdim>::memory_consumption () const
- {
- return (MemoryConsumption::memory_consumption (dofs) +
- MemoryConsumption::memory_consumption (dof_offsets));
- }
-
-
- // explicit instantiations
- template
- std::size_t
- DoFObjects<1>::memory_consumption () const;
-
- template
- std::size_t
- DoFObjects<2>::memory_consumption () const;
-
- template
- std::size_t
- DoFObjects<3>::memory_consumption () const;
- }
-}
-
-DEAL_II_NAMESPACE_CLOSE
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/integrators/*.h
+SET(_src
)
-SET(_src
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/integrators/*.h
)
DEAL_II_ADD_LIBRARY(obj_integrators OBJECT ${_src} ${_header})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/lac/*.h
- )
-
SET(_src
block_matrix_array.cc
block_sparse_matrix.cc
vector.cc
vector_memory.cc
vector_view.cc
-)
+ )
SET(_inst
block_sparse_matrix.inst.in
vector.inst.in
vector_memory.inst.in
vector_view.inst.in
-)
+ )
+
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/lac/*.h
+ )
DEAL_II_ADD_LIBRARY(obj_lac OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_lac "${_inst}")
-
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/matrix_free/*.h
- )
-
SET(_src
matrix_free.cc
)
matrix_free.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/matrix_free/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_matrix_free OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_matrix_free "${_inst}")
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/meshworker/*.h
- )
-
SET(_src
mesh_worker.cc
mesh_worker_info.cc
mesh_worker_vector_selector.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/meshworker/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_meshworker OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_meshworker "${_inst}")
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/multigrid/*.h
- )
-
SET(_src
mg_base.cc
mg_dof_handler.cc
multigrid.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/multigrid/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_multigrid OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_multigrid "${_inst}")
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/numerics/*.h
- )
-
SET(_src
data_out.cc
data_out_dof_data.cc
data_postprocessor.cc
derivative_approximation.cc
dof_output_operator.cc
- error_estimator.cc
error_estimator_1d.cc
+ error_estimator.cc
fe_field_function.cc
- matrix_tools.cc
histogram.cc
+ matrix_tools.cc
point_value_history.cc
solution_transfer.cc
solution_transfer_inst2.cc
)
SET(_inst
- data_out.inst.in
data_out_dof_data.inst.in
data_out_faces.inst.in
+ data_out.inst.in
data_out_rotation.inst.in
data_out_stack.inst.in
data_postprocessor.inst.in
derivative_approximation.inst.in
dof_output_operator.inst.in
- error_estimator.inst.in
error_estimator_1d.inst.in
+ error_estimator.inst.in
fe_field_function.inst.in
matrix_tools.inst.in
point_value_history.inst.in
vector_tools_rhs.inst.in
)
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/numerics/*.h
+ )
+
DEAL_II_ADD_LIBRARY(obj_numerics OBJECT ${_src} ${_header} ${_inst})
EXPAND_INSTANTIATIONS(obj_numerics "${_inst}")
-
template <typename number>
-inline
bool Histogram::logarithmic_less (const number n1,
const number n2)
{
values[i].end(),
logarithmic_less_function),
logarithmic_less_function);
- };
+ }
break;
- };
+ }
default:
Assert (false, ExcInternalError());
- };
+ }
// move right bound arbitrarily if
// necessary. sometimes in logarithmic
projects / dealii-svn.git / commitdiff