From 18ddf29df715d9aae7780ad6cdcb4c4d35616ee0 Mon Sep 17 00:00:00 2001 From: heltai Date: Mon, 16 Sep 2013 21:53:47 +0000 Subject: [PATCH] Merged from trunk git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30750 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/bundled/CMakeLists.txt | 6 + deal.II/cmake/config/CMakeLists.txt | 122 +- deal.II/cmake/config/Config.cmake.in | 33 +- deal.II/cmake/config/Make.global_options.in | 13 +- .../macro_deal_ii_invoke_autopilot.cmake | 2 + .../macros/macro_expand_instantiations.cmake | 14 +- deal.II/cmake/modules/FindTRILINOS.cmake | 2 +- deal.II/cmake/scripts/CMakeLists.txt | 32 +- deal.II/cmake/setup_finalize.cmake | 11 + deal.II/cmake/setup_post_project_call.cmake | 7 +- .../parameter_gui/parameter_delegate.cpp | 2 +- .../parameter_gui/parameter_delegate.h | 2 +- deal.II/doc/development/cmake.html | 16 +- deal.II/doc/development/cmakelists.html | 318 +++-- deal.II/doc/doxygen/options.dox.in | 3 +- deal.II/doc/news/changes.h | 49 + deal.II/examples/step-31/doc/intro.dox | 6 + deal.II/examples/step-31/doc/results.dox | 6 +- deal.II/examples/step-32/doc/intro.dox | 6 +- deal.II/examples/step-32/doc/results.dox | 9 + deal.II/examples/step-32/step-32.cc | 4 +- deal.II/examples/step-51/doc/intro.dox | 148 +- deal.II/examples/step-51/doc/results.dox | 24 +- deal.II/examples/step-51/step-51.cc | 424 ++++-- deal.II/include/deal.II/distributed/tria.h | 52 +- .../deal.II/dofs/dof_accessor.templates.h | 341 ++--- deal.II/include/deal.II/dofs/dof_levels.h | 29 + deal.II/include/deal.II/dofs/dof_tools.h | 31 +- deal.II/include/deal.II/grid/grid_tools.h | 26 + deal.II/include/deal.II/hp/dof_faces.h | 520 ++++++- deal.II/include/deal.II/hp/dof_handler.h | 13 +- deal.II/include/deal.II/hp/dof_level.h | 432 ++++++ deal.II/include/deal.II/hp/dof_objects.h | 540 ------- .../deal.II/lac/parallel_block_vector.h | 63 +- deal.II/include/deal.II/lac/parallel_vector.h | 29 + .../deal.II/lac/parallel_vector.templates.h | 18 +- .../include/deal.II/lac/vector.templates.h | 42 +- .../include/deal.II/numerics/vector_tools.h | 55 + .../deal.II/numerics/vector_tools.templates.h | 67 +- deal.II/source/CMakeLists.txt | 13 +- deal.II/source/algorithms/CMakeLists.txt | 8 +- deal.II/source/base/CMakeLists.txt | 12 +- deal.II/source/distributed/CMakeLists.txt | 9 +- deal.II/source/distributed/tria.cc | 1260 +++++++++++------ deal.II/source/dofs/CMakeLists.txt | 13 +- deal.II/source/dofs/dof_handler_policy.cc | 69 +- deal.II/source/dofs/dof_tools.inst.in | 17 + deal.II/source/dofs/dof_tools_constraints.cc | 47 +- deal.II/source/{base => }/dummy.cc | 0 deal.II/source/fe/CMakeLists.txt | 13 +- deal.II/source/fe/mapping_q.cc | 2 +- deal.II/source/grid/CMakeLists.txt | 10 +- deal.II/source/grid/grid_out.cc | 52 +- deal.II/source/grid/grid_tools.cc | 80 ++ deal.II/source/grid/grid_tools.inst.in | 32 + deal.II/source/grid/manifold.cc | 1 - deal.II/source/hp/CMakeLists.txt | 12 +- deal.II/source/hp/dof_faces.cc | 34 +- deal.II/source/hp/dof_handler.cc | 134 +- deal.II/source/hp/dof_level.cc | 244 ++++ deal.II/source/hp/dof_levels.cc | 43 - deal.II/source/hp/dof_objects.cc | 50 - deal.II/source/integrators/CMakeLists.txt | 6 +- deal.II/source/lac/CMakeLists.txt | 13 +- deal.II/source/matrix_free/CMakeLists.txt | 8 +- deal.II/source/meshworker/CMakeLists.txt | 8 +- deal.II/source/multigrid/CMakeLists.txt | 8 +- deal.II/source/numerics/CMakeLists.txt | 17 +- deal.II/source/numerics/histogram.cc | 7 +- 69 files changed, 3770 insertions(+), 1969 deletions(-) create mode 100644 deal.II/include/deal.II/hp/dof_level.h delete mode 100644 deal.II/include/deal.II/hp/dof_objects.h rename deal.II/source/{base => }/dummy.cc (100%) create mode 100644 deal.II/source/hp/dof_level.cc delete mode 100644 deal.II/source/hp/dof_levels.cc delete mode 100644 deal.II/source/hp/dof_objects.cc diff --git a/deal.II/bundled/CMakeLists.txt b/deal.II/bundled/CMakeLists.txt index 846b614f2a..897b4632e1 100644 --- a/deal.II/bundled/CMakeLists.txt +++ b/deal.II/bundled/CMakeLists.txt @@ -35,6 +35,9 @@ IF(FEATURE_BOOST_BUNDLED_CONFIGURED) 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) @@ -88,6 +91,9 @@ IF(FEATURE_THREADS_BUNDLED_CONFIGURED) 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() diff --git a/deal.II/cmake/config/CMakeLists.txt b/deal.II/cmake/config/CMakeLists.txt index 9bcda3e751..d640217ef4 100644 --- a/deal.II/cmake/config/CMakeLists.txt +++ b/deal.II/cmake/config/CMakeLists.txt @@ -14,13 +14,32 @@ ## ## --------------------------------------------------------------------- +# +# 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 ) @@ -30,7 +49,13 @@ CONFIGURE_FILE( # # ######################################################################## -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 @@ -93,24 +118,37 @@ FOREACH(_build ${DEAL_II_BUILD_TYPES}) ) 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" ) @@ -123,16 +161,15 @@ FOREACH(var ${res}) 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 ) @@ -150,13 +187,6 @@ IF(DEAL_II_COMPONENT_COMPAT_FILES) # 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} @@ -268,16 +298,62 @@ IF(DEAL_II_COMPONENT_COMPAT_FILES) 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 + ) diff --git a/deal.II/cmake/config/Config.cmake.in b/deal.II/cmake/config/Config.cmake.in index 150793f5cc..53417f2bdf 100644 --- a/deal.II/cmake/config/Config.cmake.in +++ b/deal.II/cmake/config/Config.cmake.in @@ -66,6 +66,14 @@ WHILE(NOT "${_path}" STREQUAL "") 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: # @@ -73,9 +81,15 @@ ENDWHILE() 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() @@ -138,8 +152,11 @@ SET(DEAL_II_STATIC_EXECUTABLE "@DEAL_II_STATIC_EXECUTABLE@") # # 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@") @@ -155,7 +172,11 @@ SET(DEAL_II_LIBRARIES "@CONFIG_LIBRARIES@") # 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@") diff --git a/deal.II/cmake/config/Make.global_options.in b/deal.II/cmake/config/Make.global_options.in index cc00835e74..3ebd2f232a 100644 --- a/deal.II/cmake/config/Make.global_options.in +++ b/deal.II/cmake/config/Make.global_options.in @@ -23,7 +23,12 @@ # -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@ @@ -82,7 +87,11 @@ lib-deal2.g = @MAKEFILE_TARGETS_DEBUG@ # 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 diff --git a/deal.II/cmake/macros/macro_deal_ii_invoke_autopilot.cmake b/deal.II/cmake/macros/macro_deal_ii_invoke_autopilot.cmake index a895702cd6..4391d4bcf0 100644 --- a/deal.II/cmake/macros/macro_deal_ii_invoke_autopilot.cmake +++ b/deal.II/cmake/macros/macro_deal_ii_invoke_autopilot.cmake @@ -85,11 +85,13 @@ MACRO(DEAL_II_INVOKE_AUTOPILOT) # 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" ) diff --git a/deal.II/cmake/macros/macro_expand_instantiations.cmake b/deal.II/cmake/macros/macro_expand_instantiations.cmake index 5516eb01de..684860d475 100644 --- a/deal.II/cmake/macros/macro_expand_instantiations.cmake +++ b/deal.II/cmake/macros/macro_expand_instantiations.cmake @@ -35,22 +35,16 @@ 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} ) diff --git a/deal.II/cmake/modules/FindTRILINOS.cmake b/deal.II/cmake/modules/FindTRILINOS.cmake index 2ff349ebac..f36c514ed0 100644 --- a/deal.II/cmake/modules/FindTRILINOS.cmake +++ b/deal.II/cmake/modules/FindTRILINOS.cmake @@ -163,7 +163,7 @@ FIND_PACKAGE_HANDLE_STANDARD_ARGS(TRILINOS DEFAULT_MSG TRILINOS_CONFIG_FOUND ) -MARK_AS_ADVANCED(TRILINOS_CONFIG_DIR) +MARK_AS_ADVANCED(TRILINOS_CONFIG EPETRA_CONFIG_H) IF(TRILINOS_FOUND) diff --git a/deal.II/cmake/scripts/CMakeLists.txt b/deal.II/cmake/scripts/CMakeLists.txt index 5ac93c43cc..96539b7835 100644 --- a/deal.II/cmake/scripts/CMakeLists.txt +++ b/deal.II/cmake/scripts/CMakeLists.txt @@ -14,14 +14,22 @@ ## ## --------------------------------------------------------------------- -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() @@ -30,10 +38,14 @@ IF(DEAL_II_COMPONENT_COMPAT_FILES) # 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) @@ -81,7 +93,11 @@ IF(DEAL_II_COMPONENT_COMPAT_FILES) 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 diff --git a/deal.II/cmake/setup_finalize.cmake b/deal.II/cmake/setup_finalize.cmake index a651931f34..0584cc7696 100644 --- a/deal.II/cmake/setup_finalize.cmake +++ b/deal.II/cmake/setup_finalize.cmake @@ -72,6 +72,17 @@ FOREACH(_build ${DEAL_II_BUILD_TYPES}) ${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 + ) ######################################################################## diff --git a/deal.II/cmake/setup_post_project_call.cmake b/deal.II/cmake/setup_post_project_call.cmake index 690acf5d5d..60aad894df 100644 --- a/deal.II/cmake/setup_post_project_call.cmake +++ b/deal.II/cmake/setup_post_project_call.cmake @@ -49,5 +49,10 @@ IF(CMAKE_CROSSCOMPILING) 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() diff --git a/deal.II/contrib/parameter_gui/parameter_delegate.cpp b/deal.II/contrib/parameter_gui/parameter_delegate.cpp index 1add7f46e0..1dd33b59cc 100644 --- a/deal.II/contrib/parameter_gui/parameter_delegate.cpp +++ b/deal.II/contrib/parameter_gui/parameter_delegate.cpp @@ -102,7 +102,7 @@ namespace dealii { 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"), diff --git a/deal.II/contrib/parameter_gui/parameter_delegate.h b/deal.II/contrib/parameter_gui/parameter_delegate.h index 155fa8edaa..3c41359266 100644 --- a/deal.II/contrib/parameter_gui/parameter_delegate.h +++ b/deal.II/contrib/parameter_gui/parameter_delegate.h @@ -39,7 +39,7 @@ namespace dealii * 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 diff --git a/deal.II/doc/development/cmake.html b/deal.II/doc/development/cmake.html index 86c663bfb9..208e2617b7 100644 --- a/deal.II/doc/development/cmake.html +++ b/deal.II/doc/development/cmake.html @@ -419,15 +419,17 @@ you may receive an error message of the form:

-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.
 

This is not a problem of CMake or deal.II, but rather a general diff --git a/deal.II/doc/development/cmakelists.html b/deal.II/doc/development/cmakelists.html index 44fb38ae3f..fe34b34891 100644 --- a/deal.II/doc/development/cmakelists.html +++ b/deal.II/doc/development/cmakelists.html @@ -10,7 +10,7 @@ - +

How to use CMake to configure your projects with deal.II

@@ -20,6 +20,7 @@
  1. Adding multiple executable targets
  2. Adding libraries and common source files
  3. +
  4. Switching build types
  5. Adding a "run" target
@@ -58,16 +59,14 @@

Simple CMakeLists.txt

- In this section, we start out with a - minimal CMakeLists.txt based on - the DEAL_II_SETUP_TARGET 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 advanced - topics. Here is a full example + In this section, we start out with a minimal CMakeLists.txt + based on the DEAL_II_SETUP_TARGET 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 advanced topics. Here is a full example (plain text version): - +

Fix plain text after finalizing!

@@ -75,10 +74,10 @@ CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
 
 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)
@@ -92,33 +91,42 @@ DEAL_II_SETUP_TARGET(mycode)
 

- Next, we find our deal.II installation, in this case requiring at - least version 8.0, which is the first version using - CMake. The HINTS 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 DEAL_II_DIR. - Matthias, what is the difference between the first - and the last entry? 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 ccmake. + Next, we find our deal.II installation with the help of the + FIND_PACKAGE command. In this case requiring at least + version 8.0, which is the first version using CMake. The + HINTS 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 DEAL_II_DIR 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 + DEAL_II_DIR. If all of these hints fail the default system + locations /usr/ and /usr/local/ are considered. + The list after HINTS can be changed according to your + preferences.

- Every CMAkeLists.txt must contain a project definition, - which we do next. Matthias, I have been happy - without the CXX...? + After finding the deal.II project, we fetch a set of cached variables + with the DEAL_II_INITIALIZE_CACHED_VARIABLES() macro. You + can inspect these for instance with ccmake. +

+ +

link to detailed explanation of the macro below

+ +

+ Every CMakeLists.txt must contain a project definition, + which we do next.

Finally, the last two lines define the executable that is to be - produced and its source code. + produced and its source code. The DEAL_II_SETUP_TARGET macro + will set up necessary incllude directories, compile defintions and the + link interface.

-

Explain the macros or link to explanation below

+

link to detailed explanation of the macro below

Adding multiple executable targets

@@ -152,6 +160,7 @@ with GLOB.
 ADD_LIBRARY(mylib libsrc1.cc libsrc2.cc libsrc3.cc)
+DEAL_II_SETUP_TARGET(mylib)
 
 ADD_EXECUTABLE(mycode mycode.cc)
 DEAL_II_SETUP_TARGET(mycode)
@@ -174,85 +183,214 @@ ADD_EXECUTABLE(mycode2 mycode2.cc common.cc)
 DEAL_II_SETUP_TARGET(mycode2)
 
-

You should be aware though that common.cc will be -compiled for each target, not only once.

+

You should be aware though that in this case common.cc 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 OBJECT "library":

+ +
+ADD_LIBRARY(common OBJECT common.cc)
+DEAL_II_SETUP_TARGET(common)
 
-

Matthias, is this correct?

+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) +
+This will compile common.cc once for the object target +common and link the resulting object file into the two +executables. + + +

Switching build types

+ +

The build type is controlled via the variable +CMAKE_BUILD_TYPE. If it is set to Debug +executables and libraries specified in your CMakeLists.txt +file will be compiled in debug mode and linked against the debug version of +the deal.II library. Contrary Release will build in optimized +mode and link against the optimized release version of deal.II. You can set +CMAKE_BUILD_TYPE with the help of ccmake or via +cmake on the command line:

- -

Adding a "run" target

- -

- 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): -

+
+$ cmake -DCMAKE_BUILD_TYPE="Debug" .
+
+$ cmake -DCMAKE_BUILD_TYPE="Release" .
+
+ +Alternatively, you can specify custom targets to switch the build type and +compile automatically: + +
+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"
+  )
+
+With that, switching the build type and compiling the project can be done +very conveniently via: +
+$ make debug
+
+$ make release
+
+ + +

Adding a "run" target

+ +

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):

 ADD_CUSTOM_TARGET(run COMMAND mycode
-  COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration")
+  COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration"
+  )
 

Advanced CMakeLists.txt

-
-  
-  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
+

This section covers some advanced topics for a user +CMakeLIists.txt file.

+ + +

Control statements

+ +

Control statements in CMake take the following form: +

+IF(<expression>)
+  ...
+ENDIF()
+
+or in long form: +
+IF(<expression1>)
+  ...
+ELSEIF(<expression2>)
+  ...
+ELSE()
+  ...
+ENDIF()
+
+Please note the (somehow uncommon) empty, opening and closing brackets +behind ELSE() and ENDIF(). +<expression> can take a multitude of different forms, +have a look at the +CMake +documentation for a complete list. Important examples are: +

Link!

+
+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"
+
+

An expression can be negated by prefixing NOT:

+
+IF(NOT <expression>)
+  ...
+ENDIF()
 
- -

Finding the deal.II library

+

Loops are implemented with the help of WHILE and +FOR statements. The former takes the same +<expression> as the IF statement:

+
+WHILE(<expression>)
+  ...
+ENDWHILE()
+
+Given a variable list containing a list, the individual +elements element can be accessed with a FOREACH +statement: +
+FOREACH(element ${list})
+  ...
+ENDFOREACH
+
+Note: It is also possible to specify the list directly: +
+FOREACH(element foo bar baz)
+  # The variable element will iterate through foo, bar and baz.
+ENDFOREACH
+
-

- Finding the deal.II library should be no more than -

-    
-    FIND_PACKAGE(deal.II REQUIRED)
-  
- in your CMakeLists.txt file. You may have to hint for the location - of the deal.IIConfig.cmake file. Either by directly - specifying deal.II_DIR to point to the path were the - deal.IIConfig.cmake file is located: - -

is deal.II_DIR still accurate? We use different -above. BTW, that is an UGLY variable name.

- -
    
-cmake -Ddeal.II_DIR=/path/to/the/config/file <...>
 
+
+

File globs

+ +

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 +CMakeLists.txt 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 +SOURCE_DIR/sources/ and add it to an executable:

+ +
+FILE(GLOB sources ${CMAKE_SOURCE_DIR}/source/*.cc)
+ADD_EXECUTABLE(mycode ${sources})
 
- or by specifying a search path via CMAKE_PREFIX_PATH, - e.g. -
-    
-    cmake -DCMAKE_PREFIX_PATH=~/workspace/local
-  
- In this case deal.IIConfig.cmake will be searched - for in -
-    
-    ~/workspace/local/
-    ~/workspace/local/lib/cmake/deal.II/
-  
-

+ +

Please be aware of one caveat of this approach: Due to the fact that +CMake is a build system generator 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 CMakeLists.txt file or to run +cmake . again by hand. + + + +

DEAL_II_SETUP_TARGET revisited

+ +

Well, write this!

+ + + +

DEAL_II_INITIALIZE_CACHED_VARIABLES revisited

+ +

Well, write this!

+ + + +

Customizing include directories and compile definitions

+ +

Well, write this!

+ + + +

External libraries

+ +

Well, write this!

+ + + +

The "run" target revisited

+ +

Well, write this!

+ + + +

Install a project

+

Autopilot style CMakeLists.txt

diff --git a/deal.II/doc/doxygen/options.dox.in b/deal.II/doc/doxygen/options.dox.in index b57c0c6291..e30740fbdd 100644 --- a/deal.II/doc/doxygen/options.dox.in +++ b/deal.II/doc/doxygen/options.dox.in @@ -477,7 +477,8 @@ EXCLUDE_PATTERNS = *.templates.h # 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 diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c918de3ea1..89ca970443 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -53,6 +53,23 @@ inconvenience this causes.
    +
  1. + New: Like the usual DoFHandler class, the hp::DoFHandler class now also + has a cache that makes operations such as cell-@>get_dof_indices(...) + faster. This should accelerate many parts of the library that deal with + hp finite elements. +
    + (Wolfgang Bangerth, 2013/09/10) +
  2. + +
  3. + New: parallel::distributed::Triangulation now supports periodic boundaries. + DoFTools::make_periodicity_constraints and similar functions are now working + on parallel::distributed::Triangulation objects. +
    + (Tobin Isaac, Craig Michoski, Daniel Arndt, 2013/09/06) +
  4. +
  5. 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. @@ -68,6 +85,38 @@ inconvenience this causes.

    Specific improvements

      +
    1. + 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.) +
      + (Matthias Maier, 2013/09/15) +
    2. + +
    3. + Fixed: DoFTools::extract_locally_*_dofs now instantiated for hp::DofHandler. +
      + (Jean-Paul Pelteret, 2013/09/11) +
    4. + +
    5. + 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. +
      + (Timo Heister, 2013/09/04) +
    6. + +
    7. + Fixed: VectorTools::compute_no_normal_flux_constraints had a bug that + only manifested on complex meshes. This is now fixed. +
      + (Chih-Che Chueh, Wolfgang Bangerth, 2013/09/04) +
    8. +
    9. New: All vector classes now have functions extract_subvector_to() that allow extracting not just a single value but a whole set. diff --git a/deal.II/examples/step-31/doc/intro.dox b/deal.II/examples/step-31/doc/intro.dox index 9f11f04b81..7676de331c 100644 --- a/deal.II/examples/step-31/doc/intro.dox +++ b/deal.II/examples/step-31/doc/intro.dox @@ -117,6 +117,12 @@ shorter than for the same amount of heat to be distributed by thermal 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.

      %Boundary and initial conditions

      diff --git a/deal.II/examples/step-31/doc/results.dox b/deal.II/examples/step-31/doc/results.dox index 36fe1bc8c2..acccb42087 100644 --- a/deal.II/examples/step-31/doc/results.dox +++ b/deal.II/examples/step-31/doc/results.dox @@ -512,7 +512,7 @@ half the size for the $Q_1$ element, a fact that may not be surprising 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 grid point 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. @@ -593,7 +593,9 @@ There are various ways to extend the current program. Of particular interest 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 diff --git a/deal.II/examples/step-32/doc/intro.dox b/deal.II/examples/step-32/doc/intro.dox index da3dbb3875..25317c67cb 100644 --- a/deal.II/examples/step-32/doc/intro.dox +++ b/deal.II/examples/step-32/doc/intro.dox @@ -17,9 +17,13 @@ The work discussed here is also presented in the following publication: M. Kronbichler, T. Heister, W. Bangerth: High Accuracy Mantle Convection Simulation through Modern Numerical - Methods, Geophysical Journal International, 2012, 191, 12-29. + Methods, Geophysical Journal International, 2012, 191, 12-29. [DOI] + +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. diff --git a/deal.II/examples/step-32/doc/results.dox b/deal.II/examples/step-32/doc/results.dox index 1bc42080b4..ca0d584e14 100644 --- a/deal.II/examples/step-32/doc/results.dox +++ b/deal.II/examples/step-32/doc/results.dox @@ -471,3 +471,12 @@ solution we are really interested in better than the criterion used here, namely the KellyErrorEstimator applied to the temperature, is able to. + + +

      Possible extensions

      + +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. diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 2e9dda22f7..6819a0b1a6 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -113,7 +113,7 @@ namespace Step32 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 */ @@ -3328,7 +3328,7 @@ namespace Step32 } } } - + joint_solution.compress(VectorOperation::insert); IndexSet locally_relevant_joint_dofs(joint_dof_handler.n_dofs()); diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox index 1c14fa37c9..4aedf0fe90 100644 --- a/deal.II/examples/step-51/doc/intro.dox +++ b/deal.II/examples/step-51/doc/intro.dox @@ -1,3 +1,9 @@ +
      + + +This program was contributed by Martin Kronbichler and Scott Miller. + +

      Introduction

      @@ -8,9 +14,13 @@ discontinuous Galkerin method for the convection-diffusion equation. 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 for each of the adjacent elements, 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 @@ -18,30 +28,29 @@ for systems of equations in dim=2 or dim=3.

      Reducing the size of the linear system

      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 -( - 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: Hybridizable discontinuous Galerkin methods for partial differential equations in continuum mechanics, Journal of Computational Physics, 2012, 231:18, 5955-5988. - [DOI] -). + [DOI]). 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:
        -
      1. 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 +
      2. 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. +
      3. Use the known skeleton values as Dirichlet data for solving local +element-level solutions. This is known as the 'local solver', and is an embarrassingly parallel element-by-element solution process. -
      4. Use the local element equations to construct the global problem for the -trace values. These are the only globally coupled unknowns. -
      5. Solve for the element-local solutions using the computed trace values.
      The above procedure also has a linear algebra interpretation and is referred to @@ -60,28 +69,27 @@ approach similar to step-20, which results in the following two steps: (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 +B and C over the skeleton variable. The block-diagonality of +A and the structure in B and C allow us to invert the +matrix A 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
      1. constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$,
      2. solving the Schur complement system for $\Lambda$, and
      3. solving the equation for U using the second equation which uses $\Lambda$.
      -The important ingredient from the linear algebra point of view is that the -matrix A 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 -B and C over the skeleton variable. The block-diagonality of -A and the structure in B and C allow us to invert the -matrix A element by element (the local solution of the Dirichlet -problem) and subtract $CA^{-1}B$ from $D$.

      Solution quality and rates of convergence

      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})$. @@ -119,24 +127,35 @@ the above equation as the first order system: 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 @@ -150,27 +169,51 @@ a new variable $\lambda$ such that @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}} &=& - - _{\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_{\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}} &=& - _{\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: +
        +
      1. 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)$. +
      2. The matrix $B$ represents the local-face coupling. These are the terms + with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable + $\lambda_h$. +
      3. The matrix $C$ represents the face-local coupling, which involves the + weighting function $\mu$ multiplying the local solutions $(\mathbf{q}_h, u_h)$. +
      4. The matrix $D$ is the face-face coupling; + terms involve both $\mu$ and $\lambda_h$. +
      + +

      Post-processing and super-convergence

      Problem specific data

      @@ -204,4 +247,5 @@ Implementation notes: - 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 + diff --git a/deal.II/examples/step-51/doc/results.dox b/deal.II/examples/step-51/doc/results.dox index 96eabb1724..5cd183d641 100644 --- a/deal.II/examples/step-51/doc/results.dox +++ b/deal.II/examples/step-51/doc/results.dox @@ -48,12 +48,18 @@ cells dofs val L2 grad L2 val L2-post 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 p converges at rate p only, as opposed +to p+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 @@ -168,10 +174,10 @@ runtime of the program. This also gives a clear indication of where an improvement makes most sense.
        -
      1. Better linear solvers: We use a GMRES iterative solver without +
      2. 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 @@ -180,5 +186,5 @@ improvement makes most sense.
      3. 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).
      diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index 93e756143e..d385f69061 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -22,7 +22,7 @@ // @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 #include #include @@ -80,6 +80,9 @@ // the simulation. #include +namespace Step51 +{ + using namespace dealii; // @sect3{Equation data} @@ -295,21 +298,21 @@ double RightHandSide::value (const Point &p, * 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 DoFHandler and FE +// difference is the use of three different sets of DoFHandler and FE // objects, along with the ChunkSparseMatrix 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 -class Step51 +class HDG { public: enum RefinementMode @@ -317,30 +320,41 @@ public: 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::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::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 triangulation; @@ -352,7 +366,7 @@ private: Vector solution_local; // The new finite element type and corresponding DoFHandler 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 fe; DoFHandler dof_handler; @@ -391,10 +405,15 @@ private: 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 DoFHandler and +// FiniteElement 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 -Step51::Step51 (const unsigned int degree, +HDG::HDG (const unsigned int degree, const RefinementMode refinement_mode) : fe_local (FE_DGQ(degree), dim, FE_DGQ(degree), 1), @@ -408,12 +427,67 @@ Step51::Step51 (const unsigned int degree, -// 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 DoFHandler objects. The @p solution and @p system_matrix +// objects go with the global skeleton solution. +template +void +HDG::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::type boundary_functions; + Solution solution_function; + boundary_functions[0] = &solution_function; + VectorTools::project_boundary_values (dof_handler, + boundary_functions, + QGauss(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 -struct Step51::PerTaskData +struct HDG::PerTaskData { FullMatrix cell_matrix; Vector cell_vector; @@ -427,18 +501,24 @@ struct Step51::PerTaskData 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 WorkStream. The FEValues, matrix, +// and vector objects should be familiar by now. There are two objects that +// need to be discussed: @p std::vector > +// fe_local_support_on_face and @p std::vector > +// 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 -struct Step51::ScratchData +struct HDG::ScratchData { FEValues fe_values_local; FEFaceValues fe_face_values_local; @@ -465,7 +545,6 @@ struct Step51::ScratchData RightHandSide right_hand_side; const Solution exact_solution; - // Full constructor ScratchData(const FiniteElement &fe, const FiniteElement &fe_local, const QGauss &quadrature_formula, @@ -507,7 +586,6 @@ struct Step51::ScratchData } } - // Copy constructor ScratchData(const ScratchData &sd) : fe_values_local (sd.fe_values_local.get_fe(), @@ -534,17 +612,16 @@ struct Step51::ScratchData 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 reset(), but - // need to have it for the WorkStream interface. - void reset() {} -}; +// @sect4{HDG::PostProcessScratchData} +// @p PostProcessScratchData contains the data used by WorkStream +// when post-processing the local solution $u^*$. It is similar, but much +// simpler, than @p ScratchData. template -struct Step51::PostProcessScratchData +struct HDG::PostProcessScratchData { FEValues fe_values_local; FEValues fe_values; @@ -556,7 +633,6 @@ struct Step51::PostProcessScratchData Vector cell_rhs; Vector cell_sol; - // Full constructor PostProcessScratchData(const FiniteElement &fe, const FiniteElement &fe_local, const QGauss &quadrature_formula, @@ -572,7 +648,6 @@ struct Step51::PostProcessScratchData cell_sol (fe.dofs_per_cell) {} - // Copy constructor PostProcessScratchData(const PostProcessScratchData &sd) : fe_values_local (sd.fe_values_local.get_fe(), @@ -587,18 +662,15 @@ struct Step51::PostProcessScratchData 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 -void Step51::copy_local_to_global(const PerTaskData &data) +void HDG::copy_local_to_global(const PerTaskData &data) { if (data.trace_reconstruct == false) constraints.distribute_local_to_global (data.cell_matrix, @@ -607,49 +679,18 @@ void Step51::copy_local_to_global(const PerTaskData &data) system_matrix, system_rhs); } -template -void -Step51::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::type boundary_functions; - Solution solution_function; - boundary_functions[0] = &solution_function; - VectorTools::project_boundary_values (dof_handler, - boundary_functions, - QGauss(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 Step-32, where +// the quadrature formula and the update flags are set up, and then +// WorkStream 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 void -Step51::assemble_system (const bool trace_reconstruct) +HDG::assemble_system (const bool trace_reconstruct) { const QGauss quadrature_formula(fe.degree+1); const QGauss face_quadrature_formula(fe.degree+1); @@ -675,20 +716,25 @@ Step51::assemble_system (const bool trace_reconstruct) WorkStream::run(dof_handler.begin_active(), dof_handler.end(), *this, - &Step51::assemble_system_one_cell, - &Step51::copy_local_to_global, + &HDG::assemble_system_one_cell, + &HDG::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 void -Step51::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, +HDG::assemble_system_one_cell (const typename DoFHandler::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::active_cell_iterator loc_cell (&triangulation, cell->level(), @@ -698,12 +744,8 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce 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); @@ -713,10 +755,16 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce { 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::assemble_system_one_cell (const typename DoFHandler::active_ce } } + // 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::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); @@ -762,9 +816,18 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce const Point 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::assemble_system_one_cell (const typename DoFHandler::active_ce 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::assemble_system_one_cell (const typename DoFHandler::active_ce 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 + @@ -823,6 +898,9 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce } } + // This last term adds the contribution of the term $\left_{\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::assemble_system_one_cell (const typename DoFHandler::active_ce 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::assemble_system_one_cell (const typename DoFHandler::active_ce } } + // 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); @@ -852,6 +942,9 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce 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 set_dof_values function. else { scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); @@ -861,8 +954,11 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce +// @sect4{HDG::solve} +// The skeleton solution is solved for by using a BiCGStab solver with +// identity preconditioner. template -void Step51::solve () +void HDG::solve () { SolverControl solver_control (system_matrix.m()*10, 1e-11*system_rhs.l2_norm()); @@ -878,18 +974,41 @@ void Step51::solve () 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, [0, dim) 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 void -Step51::postprocess() +HDG::postprocess() { - // construct post-processed solution with (hopefully) higher order of - // accuracy { const QGauss quadrature_formula(fe_u_post.degree+1); const UpdateFlags local_flags (update_values); @@ -903,7 +1022,7 @@ Step51::postprocess() WorkStream::run(dof_handler_u_post.begin_active(), dof_handler_u_post.end(), - std_cxx1x::bind (&Step51::postprocess_one_cell, + std_cxx1x::bind (&HDG::postprocess_one_cell, std_cxx1x::ref(*this), std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), std_cxx1x::function(), @@ -911,7 +1030,6 @@ Step51::postprocess() 0U); } - // Compute some convergence rates, etc., and add to a table Vector difference_per_cell (triangulation.n_active_cells()); ComponentSelectFunction value_select (dim, dim+1); @@ -950,14 +1068,33 @@ Step51::postprocess() 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 void -Step51::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, +HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, PostProcessScratchData &scratch, unsigned int &) { - scratch.reset(); - typename DoFHandler::active_cell_iterator loc_cell (&triangulation, cell->level(), @@ -1009,14 +1146,27 @@ Step51::postprocess_one_cell (const typename DoFHandler::active_cell_i 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 DoFHandler +// objects. The graphical output for the skeleton variable is done through +// use of the DataOutFaces class. template -void Step51::output_results (const unsigned int cycle) +void HDG::output_results (const unsigned int cycle) { std::string filename; switch (refinement_mode) @@ -1040,18 +1190,22 @@ void Step51::output_results (const unsigned int cycle) std::ofstream output (filename.c_str()); DataOut data_out; + +// We first define the names and types of the local solution, +// and add the data to @p data_out. std::vector names (dim, "gradient"); names.push_back ("solution"); std::vector 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 post_name(1,"u_post"); std::vector post_comp_type(1, DataComponentInterpretation::component_is_scalar); @@ -1066,6 +1220,10 @@ void Step51::output_results (const unsigned int cycle) face_out += ".vtk"; std::ofstream face_output (face_out.c_str()); +// The DataOutFaces class works analagously to the DataOut +// class when we have a DoFHandler 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 data_out_face(false); std::vector face_name(1,"lambda"); std::vector @@ -1080,10 +1238,20 @@ void Step51::output_results (const unsigned int cycle) data_out_face.write_vtk (face_output); } +// @sect4{HDG::refine_grid} +// We implement two different refinement cases for HDG, just as in +// Step-7: 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 KellyErrorEstimator to +// give a decent indication of the non-regular regions in the scalar local +// solutions. template -void Step51::refine_grid (const unsigned int cycle) +void HDG::refine_grid (const unsigned int cycle) { if (cycle == 0) { @@ -1139,18 +1307,19 @@ void Step51::refine_grid (const unsigned int cycle) endc = triangulation.end(); for (; cell!=endc; ++cell) for (unsigned int face=0; face::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 Step-7. +// We loop over 10 cycles, refining the grid on each one. At the end, +// convergence tables are created. template -void Step51::run () +void HDG::run () { for (unsigned int cycle=0; cycle<10; ++cycle) { @@ -1173,6 +1342,14 @@ void Step51::run () 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 @@ -1185,6 +1362,7 @@ void Step51::run () convergence_table.write_text(std::cout); } +} // end of namespace Step51 int main (int argc, char **argv) { @@ -1203,7 +1381,7 @@ int main (int argc, char **argv) << "=============================================" << std::endl << std::endl; - Step51 hdg_problem (1, Step51::adaptive_refinement); + Step51::HDG hdg_problem (1, Step51::HDG::adaptive_refinement); hdg_problem.run (); std::cout << std::endl; @@ -1214,7 +1392,7 @@ int main (int argc, char **argv) << "===========================================" << std::endl << std::endl; - Step51 hdg_problem (1, Step51::global_refinement); + Step51::HDG hdg_problem (1, Step51::HDG::global_refinement); hdg_problem.run (); std::cout << std::endl; @@ -1225,7 +1403,7 @@ int main (int argc, char **argv) << "===========================================" << std::endl << std::endl; - Step51 hdg_problem (3, Step51::global_refinement); + Step51::HDG hdg_problem (3, Step51::HDG::global_refinement); hdg_problem.run (); std::cout << std::endl; diff --git a/deal.II/include/deal.II/distributed/tria.h b/deal.II/include/deal.II/distributed/tria.h index 6275739a61..95c34afc1d 100644 --- a/deal.II/include/deal.II/distributed/tria.h +++ b/deal.II/include/deal.II/distributed/tria.h @@ -25,7 +25,9 @@ #include #include +#include +#include #include #include #include @@ -697,6 +699,27 @@ namespace parallel const std::vector & 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. + * + * 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 >&); + + private: /** * MPI communicator to be @@ -759,6 +782,11 @@ namespace parallel * triangulation. */ typename dealii::internal::p4est::types::forest *parallel_forest; + /** + * A data structure that holds some + * information about the ghost cells of the triangulation. + */ + typename dealii::internal::p4est::types::ghost *parallel_ghost; /** * A flag that indicates @@ -892,6 +920,15 @@ namespace parallel */ 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 > + &vertices_with_ghost_neighbors); template friend class dealii::internal::DoFHandler::Policy::ParallelDistributed; }; @@ -946,7 +983,7 @@ namespace parallel * in hierarchical ordering is the ith deal cell starting * from begin(0). */ - const std::vector & + const std::vector & get_p4est_tree_to_coarse_cell_permutation() const; /** @@ -970,8 +1007,8 @@ namespace parallel * these variables at a * couple places anyway. */ - std::vector coarse_cell_to_p4est_tree_permutation; - std::vector p4est_tree_to_coarse_cell_permutation; + std::vector coarse_cell_to_p4est_tree_permutation; + std::vector p4est_tree_to_coarse_cell_permutation; /** * dummy settings @@ -990,6 +1027,15 @@ namespace parallel */ 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 > + &vertices_with_ghost_neighbors); + }; } } diff --git a/deal.II/include/deal.II/dofs/dof_accessor.templates.h b/deal.II/include/deal.II/dofs/dof_accessor.templates.h index 0a6ce9538c..6ccc484d85 100644 --- a/deal.II/include/deal.II/dofs/dof_accessor.templates.h +++ b/deal.II/include/deal.II/dofs/dof_accessor.templates.h @@ -24,7 +24,7 @@ #include #include #include -#include +#include #include #include #include @@ -2556,15 +2556,44 @@ namespace internal template static void - update_cell_dof_indices_cache (const DoFCellAccessor, lda> &) + update_cell_dof_indices_cache (const DoFCellAccessor, 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(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 dof_indices (dofs_per_cell); + static_cast,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 - static - void - get_dof_indices (const DoFCellAccessor, lda> &accessor, - std::vector &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 - static - void - get_dof_indices (const DoFCellAccessor, lda> &accessor, - std::vector &dof_indices) - { - // no caching for - // hp::DoFHandler - // implemented - typedef - dealii::DoFAccessor, 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 - static - void - get_dof_values (const DoFCellAccessor, lda> &accessor, - const InputVector &values, - ForwardIterator local_values_begin, - ForwardIterator local_values_end) - { - typedef dealii::DoFAccessor, lda> BaseClass; - Assert (static_cast(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 - static - void - get_dof_values (const DoFCellAccessor, lda> &accessor, - const InputVector &values, - ForwardIterator local_values_begin, - ForwardIterator local_values_end) - { - typedef dealii::DoFAccessor, lda> BaseClass; - - // no caching for hp::DoFHandler - // implemented - Assert (static_cast(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 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 - static - void - get_dof_values (const DoFCellAccessor, lda> &accessor, - const ConstraintMatrix &constraints, - const InputVector &values, - ForwardIterator local_values_begin, - ForwardIterator local_values_end) - { - typedef dealii::DoFAccessor, lda> BaseClass; - Assert (static_cast(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 - static - void - get_dof_values (const DoFCellAccessor, 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, lda> BaseClass; - Assert (static_cast(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 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 - static - void - set_dof_values (const DoFCellAccessor, lda> &accessor, - const dealii::Vector &local_values, - OutputVector &values) - { - typedef dealii::DoFAccessor, 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 - static - void - set_dof_values (const DoFCellAccessor, lda> &accessor, - const dealii::Vector &local_values, - OutputVector &values) - { - // no caching for hp::DoFHandler - // implemented - const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell; - - std::vector local_dof_indices (dofs_per_cell); - get_dof_indices (accessor, local_dof_indices); - - for (unsigned int i=0; i::get_dof_indices (std::vector & 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; iget_fe().dofs_per_cell; ++i, ++cache) + dof_indices[i] = *cache; } @@ -3531,12 +3332,7 @@ void DoFCellAccessor::get_dof_values (const InputVector &values, Vector &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()); } @@ -3551,8 +3347,22 @@ DoFCellAccessor::get_dof_values (const InputVector &values, { 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(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); } @@ -3568,9 +3378,22 @@ DoFCellAccessor::get_dof_values (const ConstraintMatrix &constraints, { 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(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); } @@ -3584,8 +3407,22 @@ DoFCellAccessor::set_dof_values (const Vector &local_values, { 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(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; iget_fe().dofs_per_cell; ++i, ++cache) + values(*cache) = local_values(i); } diff --git a/deal.II/include/deal.II/dofs/dof_levels.h b/deal.II/include/deal.II/dofs/dof_levels.h index a153fa1c7b..c657388fe9 100644 --- a/deal.II/include/deal.II/dofs/dof_levels.h +++ b/deal.II/include/deal.II/dofs/dof_levels.h @@ -90,6 +90,18 @@ namespace internal */ DoFObjects 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 @@ -108,6 +120,23 @@ namespace internal }; + + template + 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*dofs_per_cell+dofs_per_cell + <= + cell_dof_indices_cache.size(), + ExcInternalError()); + + return &cell_dof_indices_cache[obj_index*dofs_per_cell]; + } + + + template inline std::size_t diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index b7e1d21a7c..e4d237faf3 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -1078,9 +1078,10 @@ namespace DoFTools * 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 */ @@ -1131,8 +1132,10 @@ namespace DoFTools * 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" * @@ -1157,8 +1160,10 @@ namespace DoFTools * @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" * @@ -1190,8 +1195,10 @@ namespace DoFTools * 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" */ @@ -1217,8 +1224,10 @@ namespace DoFTools * 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" */ diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index eabd288760..83e40cf654 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -1153,6 +1154,31 @@ namespace GridTools 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 + 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 > + &periodicity_vector); /** diff --git a/deal.II/include/deal.II/hp/dof_faces.h b/deal.II/include/deal.II/hp/dof_faces.h index 2937bf902d..0f3f2cbb0e 100644 --- a/deal.II/include/deal.II/hp/dof_faces.h +++ b/deal.II/include/deal.II/hp/dof_faces.h @@ -19,8 +19,6 @@ #include -#include -#include #include #include @@ -38,6 +36,218 @@ 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). + * + *

      Offset computations

      + * + * 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 line_index * + * dof_handler.get_fe().dofs_per_line. 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 line_dof_offsets[line_index] + * within the line_dofs array. + * + * + *

      Multiple data sets per object

      + * + * 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 + * lines.dofs[lines.dof_offsets[line_index]] 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 + 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 dof_offsets; + + /** + * Store the global indices of + * the degrees of freedom. See + * DoFLevel() for detailed + * information. + */ + std::vector 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 + void + set_dof_index (const dealii::hp::DoFHandler &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 + types::global_dof_index + get_dof_index (const dealii::hp::DoFHandler &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 + unsigned int + n_active_fe_indices (const dealii::hp::DoFHandler &dof_handler, + const unsigned int obj_index) const; + + /** + * Return the fe_index of the + * n-th active finite element + * on this object. + */ + template + types::global_dof_index + nth_active_fe_index (const dealii::hp::DoFHandler &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 + bool + fe_index_is_active (const dealii::hp::DoFHandler &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 @@ -64,7 +274,7 @@ namespace internal * @author Tobias Leicht, 2006 */ template - class DoFFaces; + class DoFIndicesOnFaces; /** @@ -75,7 +285,7 @@ namespace internal * @author Tobias Leicht, 2006 */ template<> - class DoFFaces<1> + class DoFIndicesOnFaces<1> { public: /** @@ -93,13 +303,13 @@ namespace internal * @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 @@ -117,18 +327,18 @@ namespace internal * @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 @@ -138,6 +348,298 @@ namespace internal std::size_t memory_consumption () const; }; + + // --------------------- inline and template functions ------------------ + + template + template + inline + types::global_dof_index + DoFIndicesOnFacesOrEdges:: + get_dof_index (const dealii::hp::DoFHandler &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::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(), + ExcIndexRange(local_index, 0, + dof_handler.get_fe()[fe_index] + .template n_dofs_per_object())); + 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( + dof_handler.get_fe()[*pointer] + .template n_dofs_per_object() + 1); + } + } + + + + template + template + inline + void + DoFIndicesOnFacesOrEdges:: + set_dof_index (const dealii::hp::DoFHandler &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::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(), + ExcIndexRange(local_index, 0, + dof_handler.get_fe()[fe_index] + .template n_dofs_per_object())); + 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() + 1; + } + } + + + + template + template + inline + unsigned int + DoFIndicesOnFacesOrEdges:: + n_active_fe_indices (const dealii::hp::DoFHandler &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() + 1; + } + } + } + + + + template + template + inline + types::global_dof_index + DoFIndicesOnFacesOrEdges:: + nth_active_fe_index (const dealii::hp::DoFHandler &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() + 1; + } + } + + + + template + template + inline + bool + DoFIndicesOnFacesOrEdges:: + fe_index_is_active (const dealii::hp::DoFHandler &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(dof_offsets.size()))); + Assert ((fe_index != dealii::hp::DoFHandler::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( + dof_handler.get_fe()[*pointer] + .template n_dofs_per_object()+1); + } + } + + } // namespace hp } // namespace internal diff --git a/deal.II/include/deal.II/hp/dof_handler.h b/deal.II/include/deal.II/hp/dof_handler.h index d4ff6b16f7..4ee8bab130 100644 --- a/deal.II/include/deal.II/hp/dof_handler.h +++ b/deal.II/include/deal.II/hp/dof_handler.h @@ -38,9 +38,9 @@ namespace internal { namespace hp { - template class DoFLevel; - template class DoFFaces; - template class DoFObjects; + class DoFLevel; + template class DoFIndicesOnFaces; + template class DoFIndicesOnFacesOrEdges; namespace DoFHandler { @@ -790,7 +790,7 @@ namespace hp * levels[] tree of * the Triangulation objects. */ - std::vector*> levels; + std::vector levels; /** * Space to store the DoF @@ -799,7 +799,7 @@ namespace hp * faces pointer of * the Triangulation objects. */ - dealii::internal::hp::DoFFaces *faces; + dealii::internal::hp::DoFIndicesOnFaces *faces; /** * A structure that contains all @@ -900,8 +900,7 @@ namespace hp * the functions that set and * retrieve vertex dof indices. */ - template friend class dealii::internal::hp::DoFLevel; - template friend class dealii::internal::hp::DoFObjects; + template friend class dealii::internal::hp::DoFIndicesOnFacesOrEdges; friend struct dealii::internal::hp::DoFHandler::Implementation; }; diff --git a/deal.II/include/deal.II/hp/dof_level.h b/deal.II/include/deal.II/hp/dof_level.h new file mode 100644 index 0000000000..80c4e24b3f --- /dev/null +++ b/deal.II/include/deal.II/hp/dof_level.h @@ -0,0 +1,432 @@ +// --------------------------------------------------------------------- +// $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 +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace hp +{ + template class DoFHandler; + template 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 interior 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. + * + *

      Compression

      + * + * 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 12, 13, 14, 15. 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 12, 14, 13, 15, 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_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 + * (std::vector::size_type)(-1) + * for it. + * + * The type we store is then obviously the type the @p dof_indices array + * uses for indexing. + */ + std::vector 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 dof_indices; + + /** + * The offsets for each cell of the cache that holds all DoF indices. + */ + std::vector 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 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 + void compress_data (const dealii::hp::FECollection &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 + void uncompress_data (const dealii::hp::FECollection &fe_collection); + + /** + * Make hp::DoFHandler and its auxiliary class a friend since it + * is the class that needs to create these data structures. + */ + template 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 diff --git a/deal.II/include/deal.II/hp/dof_objects.h b/deal.II/include/deal.II/hp/dof_objects.h deleted file mode 100644 index 59ec451d53..0000000000 --- a/deal.II/include/deal.II/hp/dof_objects.h +++ /dev/null @@ -1,540 +0,0 @@ -// --------------------------------------------------------------------- -// $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 -#include -#include - -#include - -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). - * - *

      Offset computations

      - * - * 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 line_index * - * dof_handler.get_fe().dofs_per_line. 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 line_dof_offsets[line_index] - * within the line_dofs array. - * - * - *

      Multiple data sets per object

      - * - * 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 - * lines.dofs[lines.dof_offsets[line_index]] 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 - 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::size_type> dof_offsets; - - /** - * Store the global indices of - * the degrees of freedom. See - * DoFLevel() for detailed - * information. - */ - std::vector 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 - void - set_dof_index (const dealii::hp::DoFHandler &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 - types::global_dof_index - get_dof_index (const dealii::hp::DoFHandler &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 - unsigned int - n_active_fe_indices (const dealii::hp::DoFHandler &dof_handler, - const unsigned int obj_index) const; - - /** - * Return the fe_index of the - * n-th active finite element - * on this object. - */ - template - types::global_dof_index - nth_active_fe_index (const dealii::hp::DoFHandler &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 - bool - fe_index_is_active (const dealii::hp::DoFHandler &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 - template - inline - types::global_dof_index - DoFObjects:: - get_dof_index (const dealii::hp::DoFHandler &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::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(), - ExcIndexRange(local_index, 0, - dof_handler.get_fe()[fe_index] - .template n_dofs_per_object())); - 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( - dof_handler.get_fe()[*pointer] - .template n_dofs_per_object() + 1); - } - } - - - - template - template - inline - void - DoFObjects:: - set_dof_index (const dealii::hp::DoFHandler &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::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(), - ExcIndexRange(local_index, 0, - dof_handler.get_fe()[fe_index] - .template n_dofs_per_object())); - 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() + 1; - } - } - - - - template - template - inline - unsigned int - DoFObjects:: - n_active_fe_indices (const dealii::hp::DoFHandler &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() + 1; - } - } - } - - - - template - template - inline - types::global_dof_index - DoFObjects:: - nth_active_fe_index (const dealii::hp::DoFHandler &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() + 1; - } - } - - - - template - template - inline - bool - DoFObjects:: - fe_index_is_active (const dealii::hp::DoFHandler &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(dof_offsets.size()))); - Assert ((fe_index != dealii::hp::DoFHandler::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( - dof_handler.get_fe()[*pointer] - .template n_dofs_per_object()+1); - } - } - - } -} - -DEAL_II_NAMESPACE_CLOSE - -#endif diff --git a/deal.II/include/deal.II/lac/parallel_block_vector.h b/deal.II/include/deal.II/lac/parallel_block_vector.h index d0f5387851..747d12a33e 100644 --- a/deal.II/include/deal.II/lac/parallel_block_vector.h +++ b/deal.II/include/deal.II/lac/parallel_block_vector.h @@ -125,6 +125,20 @@ namespace parallel */ BlockVector (const std::vector &block_sizes); + /** + * Construct a block vector with an IndexSet for the local range + * and ghost entries for each block. + */ + BlockVector (const std::vector &local_ranges, + const std::vector &ghost_indices, + const MPI_Comm communicator); + + /** + * Same as above but the ghost indicies are assumed to be empty. + */ + BlockVector (const std::vector &local_ranges, + const MPI_Comm communicator); + /** * Destructor. Clears memory. */ @@ -356,6 +370,41 @@ namespace parallel } + template + inline + BlockVector::BlockVector (const std::vector &local_ranges, + const std::vector &ghost_indices, + const MPI_Comm communicator) + { + std::vector sizes(local_ranges.size()); + for (unsigned int i=0; iblock_indices.reinit(sizes); + this->components.resize(this->n_blocks()); + + for (unsigned int i=0; in_blocks(); ++i) + this->block(i).reinit(local_ranges[i], ghost_indices[i], communicator); + } + + + template + inline + BlockVector::BlockVector (const std::vector &local_ranges, + const MPI_Comm communicator) + { + std::vector sizes(local_ranges.size()); + for (unsigned int i=0; iblock_indices.reinit(sizes); + this->components.resize(this->n_blocks()); + + for (unsigned int i=0; in_blocks(); ++i) + this->block(i).reinit(local_ranges[i], communicator); + } + + template inline @@ -454,8 +503,18 @@ namespace parallel BlockVector & BlockVector::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; in_blocks(); ++i) + this->components[i] = v.block(i); + + this->collect_sizes(); return *this; } diff --git a/deal.II/include/deal.II/lac/parallel_vector.h b/deal.II/include/deal.II/lac/parallel_vector.h index 6d966154a0..7e1d853e52 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.h +++ b/deal.II/include/deal.II/lac/parallel_vector.h @@ -154,6 +154,12 @@ namespace parallel 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 @@ -204,6 +210,12 @@ namespace parallel 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 @@ -1002,6 +1014,22 @@ namespace parallel + template + inline + Vector::Vector (const IndexSet &local_range, + const MPI_Comm communicator) + : + allocated_size (0), + val (0), + import_data (0), + vector_view (0, static_cast(0)) + { + IndexSet ghost_indices(local_range.size()); + reinit (local_range, ghost_indices, communicator); + } + + + template inline Vector::Vector (const size_type size) @@ -1074,6 +1102,7 @@ namespace parallel reinit (c, true); } vector_view = c.vector_view; + update_ghost_values(); return *this; } diff --git a/deal.II/include/deal.II/lac/parallel_vector.templates.h b/deal.II/include/deal.II/lac/parallel_vector.templates.h index 23f0734d03..334966c5f0 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.templates.h +++ b/deal.II/include/deal.II/lac/parallel_vector.templates.h @@ -144,8 +144,22 @@ namespace parallel 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 new_partitioner + (new Utilities::MPI::Partitioner (locally_owned_indices, + ghost_indices, communicator)); + reinit (new_partitioner); + } + + + + template + void + Vector::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 new_partitioner (new Utilities::MPI::Partitioner (locally_owned_indices, ghost_indices, communicator)); diff --git a/deal.II/include/deal.II/lac/vector.templates.h b/deal.II/include/deal.II/lac/vector.templates.h index 8aa37a6396..b6181ffdf3 100644 --- a/deal.II/include/deal.II/lac/vector.templates.h +++ b/deal.II/include/deal.II/lac/vector.templates.h @@ -327,12 +327,10 @@ namespace internal { namespace Vector { - typedef types::global_dof_index size_type; - - template + template void set_subrange (const T s, - const size_type begin, - const size_type end, + const typename dealii::Vector::size_type begin, + const typename dealii::Vector::size_type end, dealii::Vector &dst) { if (s == T()) @@ -341,9 +339,10 @@ namespace internal std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s); } - template - void copy_subrange (const size_type begin, - const size_type end, + + template + void copy_subrange (const typename dealii::Vector::size_type begin, + const typename dealii::Vector::size_type end, const dealii::Vector &src, dealii::Vector &dst) { @@ -351,9 +350,10 @@ namespace internal (end-begin)*sizeof(T)); } - template - void copy_subrange (const size_type begin, - const size_type end, + + template + void copy_subrange (const typename dealii::Vector::size_type begin, + const typename dealii::Vector::size_type end, const dealii::Vector &src, dealii::Vector &dst) { @@ -364,15 +364,17 @@ namespace internal *p = *q; } + template - void copy_subrange_wrap (const size_type begin, - const size_type end, + void copy_subrange_wrap (const typename dealii::Vector::size_type begin, + const typename dealii::Vector::size_type end, const dealii::Vector &src, dealii::Vector &dst) { copy_subrange (begin, end, src, dst); } + template void copy_vector (const dealii::Vector &src, dealii::Vector &dst) @@ -380,8 +382,8 @@ namespace internal if (PointerComparison::equal(&src, &dst)) return; - const size_type vec_size = src.size(); - const size_type dst_size = dst.size(); + const typename dealii::Vector::size_type vec_size = src.size(); + const typename dealii::Vector::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) @@ -578,7 +580,7 @@ namespace internal // The code returns the result as the last argument in order to make // spawning tasks simpler and use automatic template deduction. template + typename ResultType, typename size_type> void accumulate (const Operation &op, const Number *X, const Number2 *Y, @@ -679,18 +681,18 @@ namespace internal ResultType r0, r1, r2, r3; Threads::TaskGroup<> task_group; task_group += Threads::new_task(&accumulate, + ResultType,size_type>, op, X, Y, power, new_size, r0); task_group += Threads::new_task(&accumulate, + ResultType,size_type>, op, X+new_size, Y+new_size, power, new_size, r1); task_group += Threads::new_task(&accumulate, + ResultType,size_type>, op, X+2*new_size, Y+2*new_size, power, new_size, r2); task_group += Threads::new_task(&accumulate, + ResultType,size_type>, op, X+3*new_size, Y+3*new_size, power, vec_size-3*new_size, r3); task_group.join_all(); diff --git a/deal.II/include/deal.II/numerics/vector_tools.h b/deal.II/include/deal.II/numerics/vector_tools.h index 5cb317c497..801e080db9 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.h +++ b/deal.II/include/deal.II/numerics/vector_tools.h @@ -2030,6 +2030,11 @@ namespace VectorTools * 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::ExcPointNotAvailableHere + * is thrown. */ template void point_difference (const DoFHandler &dof, @@ -2054,6 +2059,11 @@ namespace VectorTools * 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::ExcPointNotAvailableHere + * is thrown. */ template void point_difference (const Mapping &mapping, @@ -2077,6 +2087,11 @@ namespace VectorTools * 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::ExcPointNotAvailableHere + * is thrown. */ template void @@ -2087,6 +2102,11 @@ namespace VectorTools /** * Same as above for hp. + * + * @note If the cell in which the point is found + * is not locally owned, an exception of type + * VectorTools::ExcPointNotAvailableHere + * is thrown. */ template void @@ -2112,6 +2132,11 @@ namespace VectorTools * "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::ExcPointNotAvailableHere + * is thrown. */ template double @@ -2121,6 +2146,11 @@ namespace VectorTools /** * Same as above for hp. + * + * @note If the cell in which the point is found + * is not locally owned, an exception of type + * VectorTools::ExcPointNotAvailableHere + * is thrown. */ template double @@ -2142,6 +2172,11 @@ namespace VectorTools * 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::ExcPointNotAvailableHere + * is thrown. */ template void @@ -2153,6 +2188,11 @@ namespace VectorTools /** * Same as above for hp. + * + * @note If the cell in which the point is found + * is not locally owned, an exception of type + * VectorTools::ExcPointNotAvailableHere + * is thrown. */ template void @@ -2174,6 +2214,11 @@ namespace VectorTools * 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::ExcPointNotAvailableHere + * is thrown. */ template double @@ -2184,6 +2229,11 @@ namespace VectorTools /** * Same as above for hp. + * + * @note If the cell in which the point is found + * is not locally owned, an exception of type + * VectorTools::ExcPointNotAvailableHere + * is thrown. */ template double @@ -2301,6 +2351,11 @@ namespace VectorTools * Exception */ DeclException0 (ExcNonInterpolatingFE); + + /** + * Exception + */ + DeclException0 (ExcPointNotAvailableHere); } diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 0a906013a5..92670d6ac5 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -2787,7 +2787,7 @@ namespace VectorTools std::vector &dof_values, std::vector &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; @@ -2991,6 +2991,7 @@ namespace VectorTools { case 2: { + const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree; std::vector > tangentials (fe_values.n_quadrature_points); @@ -3024,18 +3025,18 @@ namespace VectorTools = 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 @@ -3412,15 +3413,13 @@ namespace VectorTools // 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]); } } @@ -3507,19 +3506,16 @@ namespace VectorTools // 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]); } } @@ -3611,15 +3607,13 @@ namespace VectorTools 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]); } } @@ -3705,15 +3699,13 @@ namespace VectorTools 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]); } } @@ -4250,7 +4242,8 @@ namespace VectorTools hp::FEFaceValues 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 @@ -4347,17 +4340,38 @@ namespace VectorTools // 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 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::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); @@ -5412,6 +5427,8 @@ namespace VectorTools cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); @@ -5449,6 +5466,8 @@ namespace VectorTools const std::pair::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); @@ -5485,6 +5504,8 @@ namespace VectorTools const std::pair::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); @@ -5520,6 +5541,8 @@ namespace VectorTools const std::pair::active_cell_iterator, Point > cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); diff --git a/deal.II/source/CMakeLists.txt b/deal.II/source/CMakeLists.txt index bab0c9796e..5b21083aac 100644 --- a/deal.II/source/CMakeLists.txt +++ b/deal.II/source/CMakeLists.txt @@ -63,7 +63,7 @@ FOREACH(build ${DEAL_II_BUILD_TYPES}) 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}} ) @@ -82,6 +82,9 @@ FOREACH(build ${DEAL_II_BUILD_TYPES}) 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} @@ -89,6 +92,14 @@ FOREACH(build ${DEAL_II_BUILD_TYPES}) ${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} diff --git a/deal.II/source/algorithms/CMakeLists.txt b/deal.II/source/algorithms/CMakeLists.txt index c77a4e71cf..b9a7e70f1b 100644 --- a/deal.II/source/algorithms/CMakeLists.txt +++ b/deal.II/source/algorithms/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -29,5 +25,9 @@ SET(_inst 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}") diff --git a/deal.II/source/base/CMakeLists.txt b/deal.II/source/base/CMakeLists.txt index d60e54902c..892456ad89 100644 --- a/deal.II/source/base/CMakeLists.txt +++ b/deal.II/source/base/CMakeLists.txt @@ -16,17 +16,12 @@ 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 @@ -53,9 +48,9 @@ SET(_src 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 @@ -71,7 +66,10 @@ SET(_src 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 ) - diff --git a/deal.II/source/distributed/CMakeLists.txt b/deal.II/source/distributed/CMakeLists.txt index 9648b6ae59..4bb6e252cd 100644 --- a/deal.II/source/distributed/CMakeLists.txt +++ b/deal.II/source/distributed/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -32,6 +28,9 @@ SET(_inst 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}") - diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc index d4b8767c96..ce4cbc6e97 100644 --- a/deal.II/source/distributed/tria.cc +++ b/deal.II/source/distributed/tria.cc @@ -32,12 +32,14 @@ # include # include # include +# include # include # include # include # include # include +# include #endif #include @@ -56,21 +58,11 @@ namespace internal namespace p4est { /** - * A structure whose explicit - * specializations contain - * pointers to the relevant - * p4est_* and p8est_* - * functions. Using this - * structure, for example by - * saying - * functions::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::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 struct functions; @@ -104,6 +96,10 @@ namespace internal 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, @@ -116,6 +112,16 @@ namespace internal 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); @@ -171,6 +177,9 @@ namespace internal 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); @@ -233,6 +242,10 @@ namespace internal 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, @@ -245,6 +258,16 @@ namespace internal 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; @@ -301,6 +324,10 @@ namespace internal 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) @@ -365,6 +392,10 @@ namespace internal 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, @@ -379,6 +410,14 @@ namespace internal 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); @@ -434,6 +473,9 @@ namespace internal 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); @@ -497,6 +539,10 @@ namespace internal 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, @@ -514,6 +560,16 @@ namespace internal 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, @@ -567,6 +623,10 @@ namespace internal 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) @@ -668,6 +728,37 @@ namespace internal { return functions::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 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; + }; + } } @@ -725,9 +816,8 @@ namespace /** - * 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 void @@ -740,16 +830,12 @@ namespace const bool set_vertex_info, typename internal::p4est::types::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()); @@ -769,14 +855,11 @@ namespace 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::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); @@ -792,8 +875,7 @@ namespace connectivity->tree_to_corner[index*GeometryInfo::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::faces_per_cell; ++f) if (cell->face(f)->at_boundary() == false) @@ -803,10 +885,8 @@ namespace connectivity->tree_to_tree[index*GeometryInfo::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::faces_per_cell; ++f) if (cell->face(f)->at_boundary() == false) { @@ -822,37 +902,16 @@ namespace 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] @@ -873,10 +932,8 @@ namespace 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::face_to_cell_vertices( @@ -887,10 +944,8 @@ namespace 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::vertices_per_face; ++i) { unsigned int idx @@ -1072,8 +1127,7 @@ namespace const typename internal::p4est::types::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(&tree.quadrants), &p4est_cell, internal::p4est::functions::quadrant_compare) @@ -1086,19 +1140,11 @@ namespace } 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 @@ -1131,29 +1177,18 @@ namespace &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 (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 (tree, dealii_cell->child(c), p4est_child[c], @@ -1167,66 +1202,36 @@ namespace template void - match_quadrant_recursively (const typename internal::p4est::types::quadrant &this_quadrant, - const typename Triangulation::cell_iterator &dealii_cell, - const typename internal::p4est::types::quadrant &ghost_quadrant, - const typename internal::p4est::types::forest &forest, - unsigned int ghost_owner) + match_quadrant (const dealii::Triangulation *tria, + unsigned int dealii_index, + typename internal::p4est::types::quadrant &ghost_quadrant, + unsigned int ghost_owner) { - if (internal::p4est::functions::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 (dealii_cell); - else + for (i = 0; i < l; i++) + { + typename Triangulation::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::quadrant_is_ancestor ( &this_quadrant, &ghost_quadrant)) - { - return; + child_id = internal::p4est::functions::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::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - for (unsigned int c=0; - c::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:: - quadrant_childrenv (&this_quadrant, p4est_child); - - for (unsigned int c=0; - c::max_children_per_cell; ++c) - + typename Triangulation::cell_iterator cell (tria, l, dealii_index); + if (cell->has_children()) + delete_all_children (cell); + else { - match_quadrant_recursively (p4est_child[c], dealii_cell->child(c), ghost_quadrant, forest, ghost_owner); + cell->clear_coarsen_flag(); + cell->set_subdomain_id(ghost_owner); } - - } @@ -1313,8 +1318,7 @@ namespace { //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::quadrant p4est_child[GeometryInfo::max_children_per_cell]; @@ -1371,8 +1375,7 @@ namespace } else { - //it's children got coarsened into - //this cell + //it's children got coarsened into this cell typename internal::p4est::types::quadrant *q; q = static_cast::quadrant *> ( sc_array_index (const_cast(&tree.quadrants), idx) @@ -1410,8 +1413,7 @@ namespace quadrant_overlaps_tree (const_cast::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; @@ -1494,10 +1496,9 @@ namespace /** - * A data structure that we use to store - * which cells (indicated by - * internal::p4est::types::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::quadrant objects) shall be refined and which + * shall be coarsened. */ template class RefineAndCoarsenList @@ -1509,21 +1510,12 @@ namespace typename internal::p4est::types::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 @@ -1532,12 +1524,8 @@ namespace typename internal::p4est::types::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 @@ -1578,14 +1566,13 @@ namespace template RefineAndCoarsenList:: RefineAndCoarsenList (const Triangulation &triangulation, - const std::vector &p4est_tree_to_coarse_cell_permutation, - const types::subdomain_id my_subdomain, + const std::vector &p4est_tree_to_coarse_cell_permutation, + const types::subdomain_id my_subdomain, typename internal::p4est::types::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::active_cell_iterator @@ -1606,20 +1593,13 @@ namespace 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= refine_list[i-1].p.which_tree, @@ -1718,8 +1697,7 @@ namespace RefineAndCoarsenList *this_object = reinterpret_cast*>(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; @@ -1728,29 +1706,25 @@ namespace 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:: 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:: quadrant_is_equal (quadrant, &*this_object->current_refine_pointer)) { @@ -1774,8 +1748,7 @@ namespace RefineAndCoarsenList *this_object = reinterpret_cast*>(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()) @@ -1785,29 +1758,26 @@ namespace 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:: 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:: quadrant_is_equal (children[0], &*this_object->current_coarsen_pointer)) @@ -1815,10 +1785,8 @@ namespace // 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::max_children_per_cell; ++c) { Assert (internal::p4est::functions:: @@ -1849,12 +1817,8 @@ namespace internal { 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) @@ -1874,9 +1838,7 @@ namespace internal if (deinitialized == false) { - // p4est has no - // p4est_finalize - // function + // p4est has no p4est_finalize function sc_finalize (); deinitialized = true; @@ -1885,8 +1847,8 @@ namespace internal }; 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; @@ -1920,9 +1882,7 @@ namespace parallel const typename dealii::Triangulation::MeshSmoothing smooth_grid, const Settings settings_) : - // do not check - // for distorted - // cells + // do not check for distorted cells dealii::Triangulation (smooth_grid, false), @@ -1938,11 +1898,13 @@ namespace parallel 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)); } @@ -1960,8 +1922,7 @@ namespace parallel 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); } @@ -1982,17 +1943,13 @@ namespace parallel } catch (const typename dealii::Triangulation::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; @@ -2006,10 +1963,8 @@ namespace parallel } catch (const typename Triangulation::DistortedCellList &) { - // the underlying - // triangulation should not - // be checking for - // distorted cells + // the underlying triangulation should not be checking for distorted + // cells AssertThrow (false, ExcInternalError()); } @@ -2024,6 +1979,12 @@ namespace parallel { triangulation_has_content = false; + if (parallel_ghost != 0) + { + dealii::internal::p4est::functions::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } + if (parallel_forest != 0) { dealii::internal::p4est::functions::destroy (parallel_forest); @@ -2121,6 +2082,10 @@ namespace parallel Triangulation:: load(const char *filename) { + if (parallel_ghost != 0) { + dealii::internal::p4est::functions::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } dealii::internal::p4est::functions::destroy (parallel_forest); parallel_forest = 0; dealii::internal::p4est::functions::connectivity_destroy (connectivity); @@ -2198,12 +2163,10 @@ namespace parallel 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 vertex_touch_count; std::vector< std::list< @@ -2218,12 +2181,9 @@ namespace parallel 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 @@ -2249,8 +2209,7 @@ namespace parallel 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, @@ -2265,9 +2224,8 @@ namespace parallel - // 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>) @@ -2276,12 +2234,10 @@ namespace parallel 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 vertex_touch_count; std::vector< std::list< @@ -2296,12 +2252,9 @@ namespace parallel 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 @@ -2327,8 +2280,7 @@ namespace parallel 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, @@ -2351,12 +2303,10 @@ namespace parallel 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 vertex_touch_count; std::vector< std::list< @@ -2384,9 +2334,8 @@ namespace parallel 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 @@ -2483,8 +2432,7 @@ namespace parallel 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, @@ -2503,20 +2451,26 @@ namespace parallel void Triangulation::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::MeshSmoothing - save_smooth = this->smooth_grid; - this->smooth_grid = dealii::Triangulation::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::limit_level_difference_at_vertices; + else + this->smooth_grid = dealii::Triangulation::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 @@ -2545,9 +2499,7 @@ namespace parallel } catch (const typename Triangulation::DistortedCellList &) { - // the underlying - // triangulation should not - // be checking for + // the underlying triangulation should not be checking for // distorted cells AssertThrow (false, ExcInternalError()); } @@ -2557,17 +2509,21 @@ namespace parallel // query p4est for the ghost cells - typename dealii::internal::p4est::types::ghost *ghostlayer; - ghostlayer = dealii::internal::p4est::functions::ghost_new (parallel_forest, + if (parallel_ghost != 0) + { + dealii::internal::p4est::functions::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } + parallel_ghost = dealii::internal::p4est::functions::ghost_new (parallel_forest, (dim == 2 ? typename dealii::internal::p4est::types:: - balance_type(P4EST_BALANCE_CORNER) + balance_type(P4EST_CONNECT_CORNER) : typename dealii::internal::p4est::types:: - 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 @@ -2585,14 +2541,9 @@ namespace parallel 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(parallel_forest, coarse_cell_to_p4est_tree_permutation[cell->index()]) == false) @@ -2604,13 +2555,8 @@ namespace parallel 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::quadrant p4est_coarse_cell; typename dealii::internal::p4est::types::tree *tree = @@ -2625,38 +2571,30 @@ namespace parallel } } - // 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::quadrant *quadr; unsigned int ghost_owner=0; + typename dealii::internal::p4est::types::topidx ghost_tree=0; - for (unsigned int g_idx=0; g_idxghosts.elem_count; ++g_idx) + for (unsigned int g_idx=0; g_idxghosts.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::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::cell_iterator - cell (this, 0U, coarse_cell_index); + p4est_tree_to_coarse_cell_permutation[ghost_tree]; - typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; - dealii::internal::p4est::init_coarse_quadrant (p4est_coarse_cell); - - match_quadrant_recursively (p4est_coarse_cell, cell, *quadr, - *parallel_forest, ghost_owner); + match_quadrant (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 @@ -2671,10 +2609,8 @@ namespace parallel 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; @@ -2684,9 +2620,7 @@ namespace parallel } catch (const typename Triangulation::DistortedCellList &) { - // the underlying - // triangulation should not - // be checking for + // the underlying triangulation should not be checking for // distorted cells AssertThrow (false, ExcInternalError()); } @@ -2710,7 +2644,7 @@ namespace parallel ++num_ghosts; } - Assert( num_ghosts == ghostlayer->ghosts.elem_count, ExcInternalError()); + Assert( num_ghosts == parallel_ghost->ghosts.elem_count, ExcInternalError()); #endif @@ -2720,9 +2654,10 @@ namespace parallel // 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; @@ -2739,7 +2674,8 @@ namespace parallel } } - //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 > marked_vertices(this->n_levels()); for (unsigned int lvl=0; lvln_levels(); ++lvl) @@ -2784,7 +2720,8 @@ namespace parallel { 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) @@ -2795,24 +2732,44 @@ namespace parallel } } + //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::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::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(); @@ -2825,8 +2782,7 @@ namespace parallel 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::active_cell_iterator cell = this->begin_active(); @@ -2841,8 +2797,6 @@ namespace parallel } - dealii::internal::p4est::functions::ghost_destroy (ghostlayer); - this->smooth_grid = save_smooth; } @@ -2852,8 +2806,7 @@ namespace parallel void Triangulation::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::execute_coarsening_and_refinement (); @@ -2884,17 +2837,12 @@ namespace parallel } } - // 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::active_cell_iterator cell = this->begin_active(); cell != this->end(); ++cell) @@ -2905,26 +2853,27 @@ namespace parallel } - // 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 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::ghost_destroy (parallel_ghost); + parallel_ghost = 0; + } dealii::internal::p4est::functions:: refine (parallel_forest, /* refine_recursive */ false, &RefineAndCoarsenList::refine_callback, @@ -2934,8 +2883,7 @@ namespace parallel &RefineAndCoarsenList::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()); @@ -2949,32 +2897,26 @@ namespace parallel (dim == 2 ? typename dealii::internal::p4est::types:: - balance_type(P4EST_BALANCE_FULL) + balance_type(P4EST_CONNECT_FULL) : typename dealii::internal::p4est::types:: - 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:: 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::active_cell_iterator cell = this->begin_active(); cell != this->end(); ++cell) @@ -2989,10 +2931,8 @@ namespace parallel } catch (const typename Triangulation::DistortedCellList &) { - // the underlying - // triangulation should not - // be checking for - // distorted cells + // the underlying triangulation should not be checking for distorted + // cells AssertThrow (false, ExcInternalError()); } @@ -3037,10 +2977,7 @@ namespace parallel 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(0)); number_cache.n_global_levels = Utilities::MPI::max(this->n_levels(), mpi_communicator); } @@ -3145,8 +3082,7 @@ namespace parallel dealii::internal::p4est::init_coarse_quadrant (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 (*tree, cell, @@ -3164,10 +3100,10 @@ namespace parallel } // 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) { @@ -3185,13 +3121,313 @@ namespace parallel template - const std::vector & + const std::vector & Triangulation::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 + struct find_ghosts + { + typename dealii::parallel::distributed::Triangulation *triangulation; + sc_array_t * subids; + std::map > + *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 + void + find_ghosts_corner + (typename dealii::internal::p4est::iter::corner_info * info, + void *user_data) + { + int i, j; + int nsides = info->sides.elem_count; + typename dealii::internal::p4est::iter::corner_side * sides = + (typename dealii::internal::p4est::iter::corner_side *) + (info->sides.array); + struct find_ghosts *fg = static_cast *>(user_data); + sc_array_t *subids = fg->subids; + typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; + int nsubs; + dealii::types::subdomain_id *subdomain_ids; + std::map > + *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::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(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::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 + void + find_ghosts_edge + (typename dealii::internal::p4est::iter::edge_info * info, + void *user_data) + { + int i, j, k; + int nsides = info->sides.elem_count; + typename dealii::internal::p4est::iter::edge_side * sides = + (typename dealii::internal::p4est::iter::edge_side *) + (info->sides.array); + struct find_ghosts *fg = static_cast *>(user_data); + sc_array_t *subids = fg->subids; + typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; + int nsubs; + dealii::types::subdomain_id *subdomain_ids; + std::map > + *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::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); + dealii::types::subdomain_id *subid = + static_cast(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::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 + void + find_ghosts_face + (typename dealii::internal::p4est::iter::face_info * info, + void *user_data) + { + int i, j, k; + int nsides = info->sides.elem_count; + typename dealii::internal::p4est::iter::face_side * sides = + (typename dealii::internal::p4est::iter::face_side *) + (info->sides.array); + struct find_ghosts *fg = static_cast *>(user_data); + sc_array_t *subids = fg->subids; + typename dealii::parallel::distributed::Triangulation *triangulation = fg->triangulation; + int nsubs; + dealii::types::subdomain_id *subdomain_ids; + std::map > + *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::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j])); + dealii::types::subdomain_id *subid = + static_cast(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::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 > + &vertices_with_ghost_neighbors) + { + Assert (false, ExcNotImplemented()); + } + + template <> + void + Triangulation<1,2>:: + fill_vertices_with_ghost_neighbors + (std::map > + &vertices_with_ghost_neighbors) + { + Assert (false, ExcNotImplemented()); + } + + template <> + void + Triangulation<1,3>:: + fill_vertices_with_ghost_neighbors + (std::map > + &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 > + &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(&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 > + &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(&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 > + &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(&fg), + NULL, find_ghosts_face<3,3>, find_ghosts_edge<3,3>, find_ghosts_corner<3,3>); + + sc_array_destroy (fg.subids); + } template MPI_Comm @@ -3201,6 +3437,88 @@ namespace parallel } + template + void + Triangulation::add_periodicity + (const std::vector >& + 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 + >::const_iterator periodic_tuple; + periodic_tuple = periodicity_vector.begin(); + + typename std::vector + >::const_iterator periodic_end; + periodic_end = periodicity_vector.end(); + + for (; periodic_tuple(*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:: + 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::connectivity_is_valid + (connectivity) == 1, + ExcInternalError()); + + // now create a forest out of the connectivity data structure + dealii::internal::p4est::functions::destroy (parallel_forest); + parallel_forest + = dealii::internal::p4est::functions:: + 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 std::size_t @@ -3251,17 +3569,13 @@ namespace parallel } catch (const typename dealii::Triangulation::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; @@ -3298,10 +3612,8 @@ namespace parallel } catch (const typename Triangulation::DistortedCellList &) { - // the underlying - // triangulation should not - // be checking for - // distorted cells + // the underlying triangulation should not be checking for distorted + // cells AssertThrow (false, ExcInternalError()); } @@ -3314,9 +3626,8 @@ namespace parallel Triangulation:: 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()); @@ -3360,10 +3671,33 @@ namespace parallel } } + template + typename dealii::Triangulation::cell_iterator + cell_from_quad + (typename dealii::parallel::distributed::Triangulation *triangulation, + typename dealii::internal::p4est::types::topidx treeidx, + typename dealii::internal::p4est::types::quadrant &quad) + { + int i, l = quad.level; + int child_id; + const std::vector 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::cell_iterator cell (triangulation, i, dealii_index); + child_id = internal::p4est::functions::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::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) { @@ -3414,10 +3748,10 @@ namespace parallel } template <> - const std::vector & + const std::vector & Triangulation<1,1>::get_p4est_tree_to_coarse_cell_permutation() const { - static std::vector a; + static std::vector a; return a; } diff --git a/deal.II/source/dofs/CMakeLists.txt b/deal.II/source/dofs/CMakeLists.txt index e3bf9d5631..5f2c031c9d 100644 --- a/deal.II/source/dofs/CMakeLists.txt +++ b/deal.II/source/dofs/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -38,18 +34,21 @@ SET(_src 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}") - diff --git a/deal.II/source/dofs/dof_handler_policy.cc b/deal.II/source/dofs/dof_handler_policy.cc index 30a9c96ecb..3db4b2104a 100644 --- a/deal.II/source/dofs/dof_handler_policy.cc +++ b/deal.II/source/dofs/dof_handler_policy.cc @@ -1212,6 +1212,33 @@ namespace internal } } + //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::faces_per_cell;++f) + { + if (dealii_cell->at_boundary(f)) + continue; + typename DoFHandler::level_cell_iterator neighbor = dealii_cell->neighbor(f); + if (!neighbor->has_children()) + continue; + + for (unsigned int subface=0; subface::max_children_per_face; ++subface) + { + typename DoFHandler::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 @@ -1806,6 +1833,10 @@ namespace internal 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); @@ -1861,8 +1892,8 @@ namespace internal 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 @@ -2058,25 +2089,14 @@ namespace internal 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 locally_active_vertices - = mark_locally_active_vertices (*tr); - // add each ghostcells' // subdomain to the vertex and // keep track of interesting // neighbors std::map > vertices_with_ghost_neighbors; - for (typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (cell->is_ghost ()) - for (unsigned int v=0; v::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, @@ -2320,7 +2340,8 @@ namespace internal for (typename DoFHandler::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::vertices_per_cell; ++v) if (locally_active_vertices[cell->vertex_index(v)]) vertices_with_ghost_neighbors[cell->vertex_index(v)] @@ -2364,7 +2385,7 @@ namespace internal local_dof_indices.end(), DoFHandler::invalid_dof_index)) { - Assert(false, ExcMessage ("not all dofs got distributed!")); + Assert(false, ExcMessage ("not all DoFs got distributed!")); } } } @@ -2546,10 +2567,6 @@ namespace internal 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 locally_active_vertices - = mark_locally_active_vertices (*tr); // add each ghostcells' // subdomain to the vertex and @@ -2557,14 +2574,8 @@ namespace internal // neighbors std::map > vertices_with_ghost_neighbors; - for (typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (cell->is_ghost ()) - for (unsigned int v=0; v::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 diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in index f4121a9d63..9309ae6bc2 100644 --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@ -203,18 +203,35 @@ DoFTools::extract_locally_owned_dofs > (const DoFHandler & dof_handler, IndexSet & dof_set); +template +void +DoFTools::extract_locally_owned_dofs > +(const hp::DoFHandler & dof_handler, + IndexSet & dof_set); + template void DoFTools::extract_locally_active_dofs > (const DoFHandler & dof_handler, IndexSet & dof_set); +template +void +DoFTools::extract_locally_active_dofs > +(const hp::DoFHandler & dof_handler, + IndexSet & dof_set); template void DoFTools::extract_locally_relevant_dofs > (const DoFHandler & dof_handler, IndexSet & dof_set); + +template +void +DoFTools::extract_locally_relevant_dofs > +(const hp::DoFHandler & dof_handler, + IndexSet & dof_set); template void diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc index 6e21770c30..402557bee3 100644 --- a/deal.II/source/dofs/dof_tools_constraints.cc +++ b/deal.II/source/dofs/dof_tools_constraints.cc @@ -1710,6 +1710,29 @@ namespace DoFTools 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 @@ -1983,18 +2006,6 @@ namespace DoFTools Assert (0<=direction && direction PTRIA; - const PTRIA *ptria_p = dynamic_cast (&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.")); @@ -2077,18 +2088,6 @@ namespace DoFTools 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 PTRIA; - const PTRIA *ptria_p = dynamic_cast (&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 FaceMap; diff --git a/deal.II/source/base/dummy.cc b/deal.II/source/dummy.cc similarity index 100% rename from deal.II/source/base/dummy.cc rename to deal.II/source/dummy.cc diff --git a/deal.II/source/fe/CMakeLists.txt b/deal.II/source/fe/CMakeLists.txt index 5fa933789b..57b9e686f8 100644 --- a/deal.II/source/fe/CMakeLists.txt +++ b/deal.II/source/fe/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -37,8 +33,8 @@ SET(_src 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 @@ -72,7 +68,7 @@ SET(_inst 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 @@ -96,6 +92,9 @@ SET(_inst 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}") - diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc index 26a26c7b25..c653133dfd 100644 --- a/deal.II/source/fe/mapping_q.cc +++ b/deal.II/source/fe/mapping_q.cc @@ -558,7 +558,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const // this for (unsigned int unit_point=0; unit_pointdegree*this->degree, + lohvs[unit_point].end(),0.) - 1)<1e-13*this->degree, ExcInternalError()); } diff --git a/deal.II/source/grid/CMakeLists.txt b/deal.II/source/grid/CMakeLists.txt index 7c20b596ce..18e30fe8c7 100644 --- a/deal.II/source/grid/CMakeLists.txt +++ b/deal.II/source/grid/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -48,7 +44,6 @@ SET(_inst 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 @@ -56,6 +51,9 @@ SET(_inst 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}") - diff --git a/deal.II/source/grid/grid_out.cc b/deal.II/source/grid/grid_out.cc index 1a48929745..5a204ef522 100644 --- a/deal.II/source/grid/grid_out.cc +++ b/deal.II/source/grid/grid_out.cc @@ -488,7 +488,7 @@ GridOut::default_suffix (const OutputFormat output_format) default: Assert (false, ExcNotImplemented()); return ""; - }; + } } @@ -686,9 +686,7 @@ void GridOut::write_dx (const Triangulation &tria, for (unsigned int i=0; i &tria, for (unsigned int d=spacedim+1; d<=3; ++d) out << " 0"; // fill with zeroes out << std::endl; - }; + } // Write cells preamble out << "$ENDNOD" << std::endl @@ -1008,7 +1006,7 @@ void GridOut::write_msh (const Triangulation &tria, ++vertex) out << cell->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; out << std::endl; - }; + } // write faces and lines with // non-zero boundary indicator @@ -1065,7 +1063,7 @@ void GridOut::write_ucd (const Triangulation &tria, << "# For a description of the UCD format see the AVS Developer's guide." << '\n' << "#" << '\n'; - }; + } // start with ucd data out << n_vertices << ' ' @@ -1088,7 +1086,7 @@ void GridOut::write_ucd (const Triangulation &tria, for (unsigned int d=spacedim+1; d<=3; ++d) out << " 0"; // fill with zeroes out << '\n'; - }; + } // write cells. Enumerate cells // consecutively, starting with 1 @@ -1112,7 +1110,7 @@ void GridOut::write_ucd (const Triangulation &tria, break; default: Assert (false, ExcNotImplemented()); - }; + } // it follows a list of the // vertices of each cell. in 1d @@ -1133,7 +1131,7 @@ void GridOut::write_ucd (const Triangulation &tria, ++vertex) out << cell->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; out << '\n'; - }; + } // write faces and lines with // non-zero boundary indicator @@ -2284,7 +2282,7 @@ void GridOut::write_mathgl (const Triangulation &tria, break; default: Assert (false, ExcNotImplemented ()); - }; + } out << "\n"; // (iii) write vertex ordering @@ -2308,7 +2306,7 @@ void GridOut::write_mathgl (const Triangulation &tria, break; default: Assert (false, ExcNotImplemented ()); - }; + } // (iv) write a list of vertices of cells out << "\n#" @@ -2553,7 +2551,7 @@ void GridOut::write_msh_faces (const Triangulation &tria, break; default: Assert (false, ExcNotImplemented()); - }; + } out << static_cast(face->boundary_indicator()) << ' ' << static_cast(face->boundary_indicator()) @@ -2565,7 +2563,7 @@ void GridOut::write_msh_faces (const Triangulation &tria, out << '\n'; ++index; - }; + } } @@ -2709,14 +2707,14 @@ void GridOut::write_ucd_faces (const Triangulation &tria, break; default: Assert (false, ExcNotImplemented()); - }; + } // note: vertex numbers are 1-base for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) out << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1 << ' '; out << '\n'; ++index; - }; + } } @@ -3272,7 +3270,7 @@ namespace internal { Assert(false, ExcInternalError()); break; - }; + } case 2: { @@ -3390,12 +3388,12 @@ namespace internal line_list.push_back (LineEntry(p0, p1, face->user_flag_set(), cell->level())); - }; - }; - }; + } + } + } break; - }; + } case 3: { @@ -3498,7 +3496,7 @@ namespace internal 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. @@ -3585,14 +3583,14 @@ namespace internal "/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); @@ -3637,8 +3635,8 @@ namespace internal out << ")] " << "] -6 MCshow" << '\n'; - }; - }; + } + } // and the vertex numbers if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true)) @@ -3672,8 +3670,8 @@ namespace internal << ")] " << "] -6 MCshow" << '\n'; - }; - }; + } + } out << "showpage" << '\n'; diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index a49ede55e2..eeeb0dc6ad 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -2540,6 +2541,85 @@ next_cell: } + template + 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 > + &periodicity_vector) + { + static const unsigned int dim = DH::dimension; + static const unsigned int spacedim = DH::space_dimension; + + Assert (0<=direction && direction, + Point > 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::faces_per_cell; ++f) + if(cell->face(f)->boundary_indicator() == b_id1) + { + Point face_center (cell->face(f)->center()); + face_center(direction)=0.; + const std::pair + 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::faces_per_cell;++f) + if(cell->face(f)->boundary_indicator() == b_id2) + { + typename std::map, + Point >::iterator p + = face_locations.begin(); + + Point 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 + 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 */ diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index 4411897224..fdaecb93c5 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -271,9 +271,41 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II 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 > &); + #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 &, + const types::boundary_id, + const types::boundary_id, + const unsigned int, + std::vector::cell_iterator, + unsigned int, + parallel::distributed::Triangulation::cell_iterator, + unsigned int> > &); + \} + #endif +#endif +} + + diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc index 77e93e2d5d..374c4a3a4d 100644 --- a/deal.II/source/grid/manifold.cc +++ b/deal.II/source/grid/manifold.cc @@ -102,7 +102,6 @@ template FlatManifold::FlatManifold () {} - template Point FlatManifold:: diff --git a/deal.II/source/hp/CMakeLists.txt b/deal.II/source/hp/CMakeLists.txt index 140accecb5..d360552acc 100644 --- a/deal.II/source/hp/CMakeLists.txt +++ b/deal.II/source/hp/CMakeLists.txt @@ -16,15 +16,10 @@ 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 @@ -37,6 +32,9 @@ SET(_inst 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}") - diff --git a/deal.II/source/hp/dof_faces.cc b/deal.II/source/hp/dof_faces.cc index 280ab65a9a..f7a531311f 100644 --- a/deal.II/source/hp/dof_faces.cc +++ b/deal.II/source/hp/dof_faces.cc @@ -23,9 +23,35 @@ namespace internal { namespace hp { +// ---------------------- DoFObjects ---------------------------- + template std::size_t - DoFFaces<1>::memory_consumption () const + DoFIndicesOnFacesOrEdges::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; } @@ -33,7 +59,7 @@ namespace internal std::size_t - DoFFaces<2>::memory_consumption () const + DoFIndicesOnFaces<2>::memory_consumption () const { return MemoryConsumption::memory_consumption (lines); } @@ -41,9 +67,9 @@ namespace internal 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) ); } diff --git a/deal.II/source/hp/dof_handler.cc b/deal.II/source/hp/dof_handler.cc index 8a32c2022a..d6afadf43c 100644 --- a/deal.II/source/hp/dof_handler.cc +++ b/deal.II/source/hp/dof_handler.cc @@ -18,8 +18,7 @@ #include #include #include -#include -#include +#include #include #include #include @@ -539,7 +538,7 @@ namespace internal // active_fe_indices field which // we have to backup before { - std::vector > + std::vector > active_fe_backup(dof_handler.levels.size ()); for (unsigned int level = 0; levelactive_fe_indices, @@ -553,7 +552,7 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { - dof_handler.levels.push_back (new internal::hp::DoFLevel); + dof_handler.levels.push_back (new internal::hp::DoFLevel); std::swap (active_fe_backup[level], dof_handler.levels[level]->active_fe_indices); } @@ -578,11 +577,16 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { dof_handler.levels[level]->dof_offsets - = std::vector::size_type> ( + = std::vector ( dof_handler.tria->n_raw_lines(level), - (std::vector::size_type)(-1)); + (DoFLevel::offset_type)(-1)); + dof_handler.levels[level]->cell_cache_offsets + = std::vector ( + 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::active_cell_iterator cell=dof_handler.begin_active(level); cell!=dof_handler.end_active(level); ++cell) @@ -590,11 +594,17 @@ namespace internal { 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 (next_free_dof, DoFHandler::invalid_dof_index); + dof_handler.levels[level]->cell_dof_indices_cache + = std::vector (cache_size, + DoFHandler::invalid_dof_index); } // safety check: make sure that @@ -614,12 +624,12 @@ namespace internal 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 (std::count (dof_handler.levels[level]->dof_offsets.begin(), dof_handler.levels[level]->dof_offsets.end(), - (std::vector::size_type)(-1))) + (DoFLevel::offset_type)(-1))) == dof_handler.tria->n_raw_lines(level) - dof_handler.tria->n_active_lines(level), ExcInternalError()); @@ -654,7 +664,7 @@ namespace internal // active_fe_indices field which // we have to backup before { - std::vector > + std::vector > active_fe_backup(dof_handler.levels.size ()); for (unsigned int level = 0; levelactive_fe_indices, @@ -668,11 +678,11 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { - dof_handler.levels.push_back (new internal::hp::DoFLevel); + 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>; } @@ -695,11 +705,16 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { dof_handler.levels[level]->dof_offsets - = std::vector::size_type> ( + = std::vector ( dof_handler.tria->n_raw_quads(level), - (std::vector::size_type)(-1)); + (DoFLevel::offset_type)(-1)); + dof_handler.levels[level]->cell_cache_offsets + = std::vector ( + 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::active_cell_iterator cell=dof_handler.begin_active(level); cell!=dof_handler.end_active(level); ++cell) @@ -707,11 +722,17 @@ namespace internal { 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 (next_free_dof, DoFHandler::invalid_dof_index); + dof_handler.levels[level]->cell_dof_indices_cache + = std::vector (cache_size, + DoFHandler::invalid_dof_index); } // safety check: make sure that @@ -731,12 +752,12 @@ namespace internal 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 (std::count (dof_handler.levels[level]->dof_offsets.begin(), dof_handler.levels[level]->dof_offsets.end(), - (std::vector::size_type)(-1))) + (DoFLevel::offset_type)(-1))) == dof_handler.tria->n_raw_quads(level) - dof_handler.tria->n_active_quads(level), ExcInternalError()); @@ -750,7 +771,7 @@ namespace internal // 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 @@ -847,9 +868,8 @@ namespace internal // active ones will have a // non-invalid value later on dof_handler.faces->lines.dof_offsets - = std::vector::size_type> - (dof_handler.tria->n_raw_lines(), - (std::vector::size_type)(-1)); + = std::vector (dof_handler.tria->n_raw_lines(), + (unsigned int)(-1)); dof_handler.faces->lines.dofs = std::vector (n_line_slots, DoFHandler::invalid_dof_index); @@ -1020,7 +1040,7 @@ namespace internal // active_fe_indices field which // we have to backup before { - std::vector > + std::vector > active_fe_backup(dof_handler.levels.size ()); for (unsigned int level = 0; levelactive_fe_indices, @@ -1034,11 +1054,11 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { - dof_handler.levels.push_back (new internal::hp::DoFLevel); + 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>; } @@ -1061,11 +1081,16 @@ namespace internal for (unsigned int level=0; leveln_levels(); ++level) { dof_handler.levels[level]->dof_offsets - = std::vector::size_type> - (dof_handler.tria->n_raw_hexs(level), - (std::vector::size_type)(-1)); + = std::vector ( + dof_handler.tria->n_raw_hexs(level), + (DoFLevel::offset_type)(-1)); + dof_handler.levels[level]->cell_cache_offsets + = std::vector ( + 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::active_cell_iterator cell=dof_handler.begin_active(level); cell!=dof_handler.end_active(level); ++cell) @@ -1073,11 +1098,17 @@ namespace internal { 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 (next_free_dof, DoFHandler::invalid_dof_index); + dof_handler.levels[level]->cell_dof_indices_cache + = std::vector (cache_size, + DoFHandler::invalid_dof_index); } // safety check: make sure that @@ -1097,12 +1128,12 @@ namespace internal 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 (std::count (dof_handler.levels[level]->dof_offsets.begin(), dof_handler.levels[level]->dof_offsets.end(), - (std::vector::size_type)(-1))) + (DoFLevel::offset_type)(-1))) == dof_handler.tria->n_raw_hexs(level) - dof_handler.tria->n_active_hexs(level), ExcInternalError()); @@ -1116,7 +1147,7 @@ namespace internal // 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 @@ -1217,9 +1248,9 @@ namespace internal if (true) { dof_handler.faces->quads.dof_offsets - = std::vector::size_type> + = std::vector (dof_handler.tria->n_raw_quads(), - (std::vector::size_type)(-1)); + (unsigned int)(-1)); dof_handler.faces->quads.dofs = std::vector (n_quad_slots, DoFHandler::invalid_dof_index); @@ -2691,6 +2722,14 @@ namespace hp = std::vector (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; levelcompress_data (*finite_elements); + // finally restore the user flags const_cast &>(*tria).load_user_flags(user_flags); } @@ -2727,7 +2766,22 @@ namespace hp } #endif + // uncompress the internal storage scheme of dofs on cells + // so that we can access dofs in turns + for (unsigned int level=0; leveluncompress_data (*finite_elements); + + // do the renumbering renumber_dofs_internal (new_numbers, dealii::internal::int2type()); + + // 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; levelcompress_data (*finite_elements); } @@ -3088,7 +3142,7 @@ namespace hp // Create sufficiently many // hp::DoFLevels. while (levels.size () < tria->n_levels ()) - levels.push_back (new dealii::internal::hp::DoFLevel); + levels.push_back (new dealii::internal::hp::DoFLevel); // then make sure that on each // level we have the appropriate @@ -3161,12 +3215,10 @@ namespace hp { 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); + // 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. @@ -3176,6 +3228,8 @@ namespace hp levels.pop_back (); } + Assert(levels.size () == tria->n_levels (), ExcInternalError()); + // Resize active_fe_indices // vectors. use zero indicator to // extend diff --git a/deal.II/source/hp/dof_level.cc b/deal.II/source/hp/dof_level.cc new file mode 100644 index 0000000000..fc11565a33 --- /dev/null +++ b/deal.II/source/hp/dof_level.cc @@ -0,0 +1,244 @@ +// --------------------------------------------------------------------- +// $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 +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + namespace hp + { + template + void + DoFLevel::compress_data (const dealii::hp::FECollection &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(), + 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 new_dof_indices; + new_dof_indices.reserve(new_size); + for (unsigned int cell=0; cell(), + 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= 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 + void + DoFLevel::uncompress_data(const dealii::hp::FECollection &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(); + } + + // now allocate the new array and copy into it whatever we need + std::vector new_dof_indices; + new_dof_indices.reserve(new_size); + std::vector new_dof_offsets (dof_offsets.size(), (offset_type)(-1)); + for (unsigned int cell=0; cell=0) + { + // apparently not. simply copy them + Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object(), + ExcInternalError()); + for (unsigned int i=dof_offsets[cell]; i(); ++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 diff --git a/deal.II/source/hp/dof_levels.cc b/deal.II/source/hp/dof_levels.cc deleted file mode 100644 index 6e28511809..0000000000 --- a/deal.II/source/hp/dof_levels.cc +++ /dev/null @@ -1,43 +0,0 @@ -// --------------------------------------------------------------------- -// $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 -#include - -DEAL_II_NAMESPACE_OPEN - -namespace internal -{ - namespace hp - { - template - std::size_t - DoFLevel::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 diff --git a/deal.II/source/hp/dof_objects.cc b/deal.II/source/hp/dof_objects.cc deleted file mode 100644 index 361e15da2a..0000000000 --- a/deal.II/source/hp/dof_objects.cc +++ /dev/null @@ -1,50 +0,0 @@ -// --------------------------------------------------------------------- -// $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 -#include - -DEAL_II_NAMESPACE_OPEN - -namespace internal -{ - namespace hp - { - template - std::size_t - DoFObjects::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 diff --git a/deal.II/source/integrators/CMakeLists.txt b/deal.II/source/integrators/CMakeLists.txt index 45ba6b7b61..ba7be39101 100644 --- a/deal.II/source/integrators/CMakeLists.txt +++ b/deal.II/source/integrators/CMakeLists.txt @@ -16,11 +16,11 @@ 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}) diff --git a/deal.II/source/lac/CMakeLists.txt b/deal.II/source/lac/CMakeLists.txt index bc2e40555c..ebdcde19de 100644 --- a/deal.II/source/lac/CMakeLists.txt +++ b/deal.II/source/lac/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -79,7 +75,7 @@ SET(_src vector.cc vector_memory.cc vector_view.cc -) + ) SET(_inst block_sparse_matrix.inst.in @@ -99,8 +95,11 @@ SET(_inst 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}") - diff --git a/deal.II/source/matrix_free/CMakeLists.txt b/deal.II/source/matrix_free/CMakeLists.txt index ee5c8bc58a..b3fe95a237 100644 --- a/deal.II/source/matrix_free/CMakeLists.txt +++ b/deal.II/source/matrix_free/CMakeLists.txt @@ -16,10 +16,6 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) -FILE(GLOB _header - ${CMAKE_SOURCE_DIR}/include/deal.II/matrix_free/*.h - ) - SET(_src matrix_free.cc ) @@ -28,5 +24,9 @@ SET(_inst 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}") diff --git a/deal.II/source/meshworker/CMakeLists.txt b/deal.II/source/meshworker/CMakeLists.txt index f58b0fcf63..7819483220 100644 --- a/deal.II/source/meshworker/CMakeLists.txt +++ b/deal.II/source/meshworker/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -31,5 +27,9 @@ SET(_inst 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}") diff --git a/deal.II/source/multigrid/CMakeLists.txt b/deal.II/source/multigrid/CMakeLists.txt index 40c223f30d..36738e2b50 100644 --- a/deal.II/source/multigrid/CMakeLists.txt +++ b/deal.II/source/multigrid/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -40,5 +36,9 @@ SET(_inst 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}") diff --git a/deal.II/source/numerics/CMakeLists.txt b/deal.II/source/numerics/CMakeLists.txt index df05c2eb6c..3f8b5fa7af 100644 --- a/deal.II/source/numerics/CMakeLists.txt +++ b/deal.II/source/numerics/CMakeLists.txt @@ -16,10 +16,6 @@ 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 @@ -29,11 +25,11 @@ SET(_src 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 @@ -51,16 +47,16 @@ SET(_src ) 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 @@ -76,6 +72,9 @@ SET(_inst 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}") - diff --git a/deal.II/source/numerics/histogram.cc b/deal.II/source/numerics/histogram.cc index c1c89e903a..ac21b650a4 100644 --- a/deal.II/source/numerics/histogram.cc +++ b/deal.II/source/numerics/histogram.cc @@ -24,7 +24,6 @@ DEAL_II_NAMESPACE_OPEN template -inline bool Histogram::logarithmic_less (const number n1, const number n2) { @@ -120,14 +119,14 @@ void Histogram::evaluate (const std::vector > &values, values[i].end(), logarithmic_less_function), logarithmic_less_function); - }; + } break; - }; + } default: Assert (false, ExcInternalError()); - }; + } // move right bound arbitrarily if // necessary. sometimes in logarithmic -- 2.39.5