From 16c179da49eda913ac7ffe788653cb92bc8db1c4 Mon Sep 17 00:00:00 2001
From: heltai
Date: Tue, 19 Nov 2013 10:03:16 +0000
Subject: [PATCH] Merged from trunk
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@31704 0785d39b-7218-0410-832d-ea1e28bc413d
---
deal.II/bundled/functionparser/CMakeLists.txt | 3 +
deal.II/cmake/macros/macro_add_test.cmake | 23 +++-
deal.II/cmake/macros/macro_pickup_tests.cmake | 1 -
.../cmake/macros/macro_register_feature.cmake | 9 +-
deal.II/cmake/scripts/run_testsuite.cmake | 1 -
deal.II/cmake/setup_testsuite.cmake | 31 +++--
deal.II/doc/developers/testsuite.html | 18 ++-
deal.II/doc/news/changes.h | 43 ++++++-
deal.II/doc/users/Config.sample | 1 +
deal.II/include/deal.II/base/graph_coloring.h | 59 ++++++---
deal.II/include/deal.II/dofs/dof_handler.h | 12 +-
.../include/deal.II/dofs/dof_handler_policy.h | 22 ++++
deal.II/include/deal.II/grid/grid_generator.h | 2 +-
.../deal.II/lac/trilinos_sparse_matrix.h | 2 +-
.../deal.II/lac/trilinos_vector_base.h | 31 +++++
.../include/deal.II/numerics/vector_tools.h | 41 +++---
.../deal.II/numerics/vector_tools.templates.h | 86 ++++++++++++-
deal.II/source/dofs/dof_handler_policy.cc | 17 ++-
deal.II/source/dofs/dof_renumbering.inst.in | 23 ++--
deal.II/source/fe/fe_tools_interpolate.cc | 118 +++++++++++++----
deal.II/source/fe/mapping_q1.cc | 119 ++++++++++++------
deal.II/source/grid/grid_generator.cc | 40 +++++-
deal.II/source/grid/grid_generator.inst.in | 11 ++
deal.II/source/grid/grid_out.cc | 1 +
deal.II/source/lac/chunk_sparsity_pattern.cc | 9 ++
.../numerics/vector_tools_interpolate.inst.in | 4 +-
deal.II/tests/CMakeLists.txt | 17 ++-
deal.II/tests/quick_tests/CMakeLists.txt | 9 +-
deal.II/tests/quick_tests/step-petsc.cc | 2 +-
deal.II/tests/quick_tests/tbb.cc | 2 +-
30 files changed, 598 insertions(+), 159 deletions(-)
diff --git a/deal.II/bundled/functionparser/CMakeLists.txt b/deal.II/bundled/functionparser/CMakeLists.txt
index 9fb68b95dc..eaf39a348c 100644
--- a/deal.II/bundled/functionparser/CMakeLists.txt
+++ b/deal.II/bundled/functionparser/CMakeLists.txt
@@ -26,3 +26,6 @@ DEAL_II_ADD_LIBRARY(obj_functionparser OBJECT
fpoptimizer.cc
)
+IF(DEAL_II_WITH_THREADS)
+ DEAL_II_ADD_DEFINITIONS(obj_functionparser "FP_USE_THREAD_SAFE_EVAL")
+ENDIF()
diff --git a/deal.II/cmake/macros/macro_add_test.cmake b/deal.II/cmake/macros/macro_add_test.cmake
index fc95f00172..8c3843a2fb 100644
--- a/deal.II/cmake/macros/macro_add_test.cmake
+++ b/deal.II/cmake/macros/macro_add_test.cmake
@@ -62,12 +62,29 @@ MACRO(DEAL_II_ADD_TEST _category _test_name _comparison_file)
#
SET(_configuration)
- IF(_file MATCHES "debug")
+ IF(_file MATCHES "\\.debug\\.")
SET(_configuration DEBUG)
- ELSEIF(_file MATCHES "release")
+ ELSEIF(_file MATCHES "\\.release\\.")
SET(_configuration RELEASE)
ENDIF()
+ #
+ # A "binary" in the output file indicates binary output. In this case we
+ # have to switch to plain diff instead of (possibly) numdiff, which can
+ # only work on plain text files.
+ #
+ # TODO: The indirection with ${${_test_diff_variable}} is necessary to
+ # avoid quoting issues with the command line :-/ - come up with a fix for
+ # that
+ #
+
+ SET(_test_diff_variable TEST_DIFF)
+ IF(_file MATCHES "\\.binary\\.")
+ IF(NOT DIFF_EXECUTABLE MATCHES "-NOTFOUND")
+ SET(_test_diff_variable DIFF_EXECUTABLE)
+ ENDIF()
+ ENDIF()
+
#
# Determine whether the test should be run with mpirun:
#
@@ -191,7 +208,7 @@ MACRO(DEAL_II_ADD_TEST _category _test_name _comparison_file)
COMMAND rm -f ${_test_directory}/failing_diff
COMMAND touch ${_test_directory}/diff
COMMAND
- ${TEST_DIFF}
+ ${${_test_diff_variable}} # see comment in line 72
${_test_directory}/output
${_comparison_file}
> ${_test_directory}/diff
diff --git a/deal.II/cmake/macros/macro_pickup_tests.cmake b/deal.II/cmake/macros/macro_pickup_tests.cmake
index 1e840dee4a..45b69eddce 100644
--- a/deal.II/cmake/macros/macro_pickup_tests.cmake
+++ b/deal.II/cmake/macros/macro_pickup_tests.cmake
@@ -79,7 +79,6 @@ Invalid feature constraint \"${_match}\" in file
\"${_comparison}\":
The feature \"DEAL_II_${_feature}\" does not exist.\n"
)
- #SET(_define_test FALSE)
ENDIF()
ENDFOREACH()
diff --git a/deal.II/cmake/macros/macro_register_feature.cmake b/deal.II/cmake/macros/macro_register_feature.cmake
index a4408fe8d9..4184a22c27 100644
--- a/deal.II/cmake/macros/macro_register_feature.cmake
+++ b/deal.II/cmake/macros/macro_register_feature.cmake
@@ -23,8 +23,7 @@
# This macro will
#
# - add ${feature}_INCLUDE_DIRS and ${feature}_INCLUDE_PATH
-# to the list of (internal) include
-# directories
+# to the list of (internal) include directories
# - and if ${feature}_ADD_TO_USER_INCLUDE_DIRS is defined also to
# DEAL_II_USER_INCLUDE_DIRS
#
@@ -40,12 +39,9 @@
MACRO(REGISTER_FEATURE _feature)
- MESSAGE(STATUS "Register enabled external feature ${_feature}:")
-
# variables for include directories:
FOREACH(_var ${_feature}_INCLUDE_DIRS ${_feature}_INCLUDE_PATH)
IF(DEFINED ${_var})
- MESSAGE(STATUS " ${_var}")
INCLUDE_DIRECTORIES(${${_var}})
IF(${_feature}_ADD_TO_USER_INCLUDE_DIRS)
LIST(APPEND DEAL_II_USER_INCLUDE_DIRS ${${_var}})
@@ -56,7 +52,6 @@ MACRO(REGISTER_FEATURE _feature)
# variables for linker flags:
FOREACH(_var ${_feature}_LINKER_FLAGS ${_feature}_LINK_FLAGS)
IF(DEFINED ${_var})
- MESSAGE(STATUS " ${_var}")
ADD_FLAGS(DEAL_II_LINKER_FLAGS "${${_var}}")
ENDIF()
ENDFOREACH()
@@ -64,13 +59,11 @@ MACRO(REGISTER_FEATURE _feature)
# variables for compiler flags:
FOREACH(_var ${_feature}_CXX_FLAGS ${_feature}_COMPILE_FLAGS)
IF(DEFINED ${_var})
- MESSAGE(STATUS " ${_var}")
ADD_FLAGS(CMAKE_CXX_FLAGS "${${_var}}")
ENDIF()
ENDFOREACH()
IF(DEFINED ${_feature}_LIBRARIES)
- MESSAGE(STATUS " ${_feature}_LIBRARIES")
#
# Add ${_feature}_LIBRARIES to
# DEAL_II_EXTERNAL_LIBRARIES
diff --git a/deal.II/cmake/scripts/run_testsuite.cmake b/deal.II/cmake/scripts/run_testsuite.cmake
index f4b1444450..e952fb81cb 100644
--- a/deal.II/cmake/scripts/run_testsuite.cmake
+++ b/deal.II/cmake/scripts/run_testsuite.cmake
@@ -93,7 +93,6 @@
# TEST_TIME_LIMIT
# TEST_PICKUP_REGEX
# TEST_OVERRIDE_LOCATION
-# NUMDIFF_DIR
#
# For details, consult the ./README file.
#
diff --git a/deal.II/cmake/setup_testsuite.cmake b/deal.II/cmake/setup_testsuite.cmake
index 5aa5eb2718..e40a167517 100644
--- a/deal.II/cmake/setup_testsuite.cmake
+++ b/deal.II/cmake/setup_testsuite.cmake
@@ -65,6 +65,7 @@ ENDFOREACH()
#
# We need deal.II and Perl as external packages:
#
+
FIND_PACKAGE(deal.II 8.0 REQUIRED
HINTS ${DEAL_II_BINARY_DIR} ${DEAL_II_DIR}
)
@@ -75,8 +76,11 @@ FIND_PACKAGE(Perl REQUIRED)
#
# We need a diff tool, preferably numdiff:
#
+
FIND_PROGRAM(DIFF_EXECUTABLE
NAMES diff
+ HINTS ${DIFF_DIR}
+ PATH_SUFFIXES bin
)
FIND_PROGRAM(NUMDIFF_EXECUTABLE
@@ -87,37 +91,38 @@ FIND_PROGRAM(NUMDIFF_EXECUTABLE
MARK_AS_ADVANCED(DIFF_EXECUTABLE NUMDIFF_EXECUTABLE)
-IF( NUMDIFF_EXECUTABLE MATCHES "-NOTFOUND"
- AND DIFF_EXECUTABLE MATCHES "-NOTFOUND" )
- MESSAGE(FATAL_ERROR
- "Could not find diff or numdiff. One of those are required for running the testsuite."
- )
-ENDIF()
-
IF("${TEST_DIFF}" STREQUAL "")
+ #
+ # No TEST_DIFF is set, specify one:
+ #
+
IF(NOT NUMDIFF_EXECUTABLE MATCHES "-NOTFOUND")
- SET(TEST_DIFF ${NUMDIFF_EXECUTABLE} -a 1e-6 -s ' \\t\\n:')
+ SET(TEST_DIFF ${NUMDIFF_EXECUTABLE} -a 1e-6 -r 1e-8 -s ' \\t\\n:')
+ ELSEIF(NOT DIFF_EXECUTABLE MATCHES "-NOTFOUND")
+ SET(TEST_DIFF ${DIFF_EXECUTABLE})
ELSE()
- SET(TEST_DIFF ${DIFF_EXECUTABLE})
+ MESSAGE(FATAL_ERROR
+ "Could not find diff or numdiff. One of those are required for running the testsuite.\n"
+ "Please specify TEST_DIFF by hand."
+ )
ENDIF()
-ELSE()
- # TODO: I have no idea how to prepare a custom string comming possibly
- # through two layers of command line into a list...
- SEPARATE_ARGUMENTS(TEST_DIFF ${TEST_DIFF})
ENDIF()
#
# Set a default time limit of 600 seconds:
#
+
SET_IF_EMPTY(TEST_TIME_LIMIT 600)
#
# And finally, enable testing:
#
+
ENABLE_TESTING()
#
# A custom target that does absolutely nothing. It is used in the main
# project to trigger a "make rebuild_cache" if necessary.
#
+
ADD_CUSTOM_TARGET(regenerate)
diff --git a/deal.II/doc/developers/testsuite.html b/deal.II/doc/developers/testsuite.html
index 85f72404a6..c2691aa449 100644
--- a/deal.II/doc/developers/testsuite.html
+++ b/deal.II/doc/developers/testsuite.html
@@ -62,6 +62,7 @@
Restricting tests for build configurations
Restricting tests for feature configurations
Running tests with MPI
+ Tests with binary output
Changing condition for success
Adding new tests
@@ -425,7 +426,7 @@
just category/test.output
:
- category/test.[with_<feature>=<on|off>.]*[mpirun=<x>.][expect=<y>.][<debug|release>.]output
+ category/test.[with_<feature>=<on|off>.]*[mpirun=<x>.][expect=<y>.][binary.][<debug|release>.]output
Normally, a test will be set up so that it runs twice, once in debug and
once in release configuration.
@@ -501,6 +502,21 @@
files for different numbers of MPI processes.
+
+ Tests with binary output
+
+ If a test produces binary output add binary
to the
+ output file to indicate this:
+
+
+ category/test.binary.output
+
+ The testsuite ensures that a diff tool suitable for comparing binary
+ output files is used instead of the default diff tool, which (as in
+ the case of numdiff
) might be unable to compare binary
+ files.
+
+
Changing condition for success
diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index b5997b9a33..80b7a85a8f 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -211,6 +211,45 @@ inconvenience this causes.
Specific improvements
+ - Fixed: GridGenerator::parallelogram was not instantiated properly
+ on intel compilers.
+
+ (Timo Heister, 2013/11/17)
+
+
+ -
+ Fixed: MappingQ1::transform_real_to_unit_cell() could fail in some
+ cases with very elongated and twisted cells. This should now be fixed
+ with an algorithm that uses a better method of computing the Newton
+ convergence.
+
+ (Wolfgang Bangerth, 2013/11/17)
+
+
+ -
+ Fixed: VectorTools::compute_no_normal_flux_constraints had a bug that
+ only appeared in rare cases at vertices of the domain if one adjacent
+ cell had two boundary indicators selected for no normal flux and another
+ had only one. This is now fixed.
+
+ (Wolfgang Bangerth, 2013/11/17)
+
+
+ - Fixed: dealii::FETools::interpolation_difference was
+ not working for TrilinosWrappers::MPI::Vectors with ghost
+ entries. The TrilinosWrappers::VectorBase class has now a
+ get_mpi_communicator method similar to the PETSc vector
+ classes.
+
+ (Martin Steigemann, Martin Kronbichler, 2013/11/17)
+
+
+ - Fixed: Bundled fparser is now compiled with FP_USE_THREAD_SAFE_EVAL in
+ case of enabled threading support so that it is thread safe.
+
+ (Matthias Maier, reported by Francesco Cattoglio 2013/11/16)
+
+
- Fixed: The CellData class now has a default constructor that
sets the material and boundary indicators to zero. This fixes certain
internal errors with the Intel ICC compiler.
@@ -248,7 +287,7 @@ inconvenience this causes.
-
New: parallel::distributed::BlockVector has now methods update_ghost_values,
- compress, set_out_ghosts, and has_ghost_elements that do the respective
+ compress, zero_out_ghosts, and has_ghost_elements that do the respective
operation on each block of parallel::distributed::Vector.
(Martin Kronbichler, 2013/10/18)
@@ -344,7 +383,7 @@ inconvenience this causes.
New: All vector classes now have functions extract_subvector_to()
that allow extracting not just a single value but a whole set.
- (Fahad Alrasched, 2013/09/02)
+ (Fahad Alrashed, 2013/09/02)
-
diff --git a/deal.II/doc/users/Config.sample b/deal.II/doc/users/Config.sample
index 06b915700b..26faa6c15b 100644
--- a/deal.II/doc/users/Config.sample
+++ b/deal.II/doc/users/Config.sample
@@ -447,6 +447,7 @@
# SET(PETSC_VERSION_SUBMINOR "z" CACHE STRING "")
# SET(PETSC_VERSION_PATCH "p" CACHE STRING "")
# SET(PETSC_WITH_MPI_UNI FALSE CACHE BOOL "")
+# SET(PETSC_WITH_COMPLEX FALSE CACHE BOOL "")
#
diff --git a/deal.II/include/deal.II/base/graph_coloring.h b/deal.II/include/deal.II/base/graph_coloring.h
index 5de886009b..86601c87b5 100644
--- a/deal.II/include/deal.II/base/graph_coloring.h
+++ b/deal.II/include/deal.II/base/graph_coloring.h
@@ -39,6 +39,42 @@ namespace GraphColoring
{
namespace internal
{
+ /**
+ * Given two sets of indices that are assumed to be sorted, determine
+ * whether they will have a nonempty intersection. The actual intersection is
+ * not computed.
+ * @param indices1 A set of indices, assumed sorted.
+ * @param indices2 A set of indices, assumed sorted.
+ * @return Whether the two sets of indices do have a nonempty intersection.
+ */
+ inline
+ bool
+ have_nonempty_intersection (const std::vector &indices1,
+ const std::vector &indices2)
+ {
+ // we assume that both arrays are sorted, so we can walk
+ // them in lockstep and see if we encounter an index that's
+ // in both arrays. once we reach the end of either array,
+ // we know that there is no intersection
+ std::vector::const_iterator
+ p = indices1.begin(),
+ q = indices2.begin();
+ while ((p != indices1.end()) && (q != indices2.end()))
+ {
+ if (*p < *q)
+ ++p;
+ else if (*p > *q)
+ ++q;
+ else
+ // conflict found!
+ return true;
+ }
+
+ // no conflict found!
+ return false;
+ }
+
+
/**
* Create a partitioning of the given range of iterators using a simplified
* version of the Cuthill-McKee algorithm (Breadth First Search algorithm).
@@ -105,7 +141,9 @@ namespace GraphColoring
std::vector new_zone;
for (; previous_zone_it!=previous_zone_end; ++previous_zone_it)
{
- std::vector conflict_indices = get_conflict_indices(*previous_zone_it);
+ const std::vector
+ conflict_indices = get_conflict_indices(*previous_zone_it);
+
const unsigned int n_conflict_indices(conflict_indices.size());
for (unsigned int i=0; i > conflict_indices(partition_size);
std::vector > graph(partition_size);
- // Get the conflict indices associated to each iterator. The conflict_indices have to be sorted so
- // set_intersection can be used later.
+ // Get the conflict indices associated to each iterator. The conflict_indices have to
+ // be sorted so we can more easily find conflicts later on
for (unsigned int i=0; i conflict_indices_intersection;
- std::vector::iterator intersection_it;
for (unsigned int i=0; i::iterator degrees_it;
diff --git a/deal.II/include/deal.II/dofs/dof_handler.h b/deal.II/include/deal.II/dofs/dof_handler.h
index 93d765a670..a6a1d88142 100644
--- a/deal.II/include/deal.II/dofs/dof_handler.h
+++ b/deal.II/include/deal.II/dofs/dof_handler.h
@@ -1264,7 +1264,7 @@ void DoFHandler::save (Archive &ar,
// identifies the FE and the policy
unsigned int n_cells = tria->n_cells();
std::string fe_name = selected_fe->get_name();
- std::string policy_name = typeid(*policy).name();
+ std::string policy_name = policy->to_string();
ar &n_cells &fe_name &policy_name;
}
@@ -1305,9 +1305,13 @@ void DoFHandler::load (Archive &ar,
AssertThrow (fe_name == selected_fe->get_name(),
ExcMessage ("The finite element associated with this DoFHandler does not match "
"the one that was associated with the DoFHandler previously stored."));
- AssertThrow (policy_name == typeid(*policy).name(),
- ExcMessage ("The policy associated with this DoFHandler does not match "
- "the one that was associated with the DoFHandler previously stored."));
+ AssertThrow (policy_name == policy->to_string(),
+ ExcMessage (std::string ("The policy currently associated with this DoFHandler (")
+ + policy->to_string()
+ +std::string(") does not match the one that was associated with the "
+ "DoFHandler previously stored (")
+ + policy_name
+ + ")."));
}
diff --git a/deal.II/include/deal.II/dofs/dof_handler_policy.h b/deal.II/include/deal.II/dofs/dof_handler_policy.h
index 42eecf8afa..39a07a544f 100644
--- a/deal.II/include/deal.II/dofs/dof_handler_policy.h
+++ b/deal.II/include/deal.II/dofs/dof_handler_policy.h
@@ -64,6 +64,13 @@ namespace internal
*/
virtual ~PolicyBase ();
+ /**
+ * Convert type of policy to a string. This is used as a
+ * workaround in serialization because the intel compiler can
+ * not handle typeid correctly.
+ */
+ virtual std::string to_string () const = 0;
+
/**
* Distribute degrees of freedom on
* the object given as last argument.
@@ -108,6 +115,14 @@ namespace internal
virtual
NumberCache
distribute_dofs (dealii::DoFHandler &dof_handler) const;
+
+ /**
+ * Convert type of policy to a string. This is used as a
+ * workaround in serialization because the intel compiler can
+ * not handle typeid correctly.
+ */
+ virtual std::string to_string () const;
+
/**
* Distribute multigrid DoFs.
@@ -147,6 +162,13 @@ namespace internal
NumberCache
distribute_dofs (dealii::DoFHandler &dof_handler) const;
+ /**
+ * Convert type of policy to a string. This is used as a
+ * workaround in serialization because the intel compiler can
+ * not handle typeid correctly.
+ */
+ virtual std::string to_string () const;
+
/**
* Distribute multigrid DoFs.
*/
diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h
index a3f4021278..7307374abe 100644
--- a/deal.II/include/deal.II/grid/grid_generator.h
+++ b/deal.II/include/deal.II/grid/grid_generator.h
@@ -205,7 +205,7 @@ namespace GridGenerator
template
void
parallelogram (Triangulation &tria,
- const Point (&corners) [dim],
+ const Point corners[dim],
const bool colorize=false);
/**
diff --git a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
index 6f3327f85b..15c8731d6c 100644
--- a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
@@ -3404,7 +3404,7 @@ namespace TrilinosWrappers
}
- // If the calling matrix owns the row to which we want to add values, we
+ // If the calling processor owns the row to which we want to add values, we
// can directly call the Epetra_CrsMatrix input function, which is much
// faster than the Epetra_FECrsMatrix function.
if (row_partitioner().MyGID(static_cast(row)) == true)
diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h
index a721cdce8c..d5e2217ad9 100644
--- a/deal.II/include/deal.II/lac/trilinos_vector_base.h
+++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h
@@ -1051,6 +1051,13 @@ namespace TrilinosWrappers
* consumption in bytes.
*/
std::size_t memory_consumption () const;
+
+ /**
+ * Return a reference to the MPI
+ * communicator object in use with this
+ * object.
+ */
+ const MPI_Comm &get_mpi_communicator () const;
//@}
/**
@@ -2185,6 +2192,30 @@ namespace TrilinosWrappers
}
+
+ inline
+ const MPI_Comm &
+ VectorBase::get_mpi_communicator () const
+ {
+ static MPI_Comm comm;
+
+#ifdef DEAL_II_WITH_MPI
+
+ const Epetra_MpiComm *mpi_comm
+ = dynamic_cast(&vector->Map().Comm());
+ comm = mpi_comm->Comm();
+
+#else
+
+ comm = MPI_COMM_SELF;
+
+#endif
+
+ return comm;
+ }
+
+
+
#endif // DOXYGEN
}
diff --git a/deal.II/include/deal.II/numerics/vector_tools.h b/deal.II/include/deal.II/numerics/vector_tools.h
index d8c6be07ff..daf77a045d 100644
--- a/deal.II/include/deal.II/numerics/vector_tools.h
+++ b/deal.II/include/deal.II/numerics/vector_tools.h
@@ -1404,7 +1404,7 @@ namespace VectorTools
/**
- * Compute the constraints that correspond to boundary conditions of the
+ * This function computes the constraints that correspond to boundary conditions of the
* form $\vec n \cdot \vec u=0$, i.e. no normal flux if $\vec u$ is a
* vector-valued quantity. These conditions have exactly the form handled by
* the ConstraintMatrix class, so instead of creating a map between boundary
@@ -1563,8 +1563,19 @@ namespace VectorTools
* There are cases where one cell contributes two tangential directions and
* another one only one; for example, this would happen if both top and
* front faces of the left cell belong to the boundary selected whereas only
- * the top face of the right cell belongs to it. This case is not currently
- * implemented.
+ * the top face of the right cell belongs to it, maybe indicating the the entire
+ * front part of the domain is a smooth manifold whereas the top really forms
+ * two separate manifolds that meet in a ridge, and that no-flux boundary
+ * conditions are only desired on the front manifold and the right one on top.
+ * In cases like these, it's difficult to define what should happen. The
+ * current implementation simply ignores the one contribution from the
+ * cell that only contributes one normal vector. In the example shown, this
+ * is acceptable because the normal vector for the front face of the left
+ * cell is the same as the normal vector provided by the front face of
+ * the right cell (the surface is planar) but it would be a problem if the
+ * front manifold would be curved. Regardless, it is unclear how one would
+ * proceed in this case and ignoring the single cell is likely the best
+ * one can do.
*
*
*
Results
@@ -2032,7 +2043,7 @@ namespace VectorTools
* point_difference() function.
*
* @note If the cell in which the point is found
- * is not locally owned, an exception of type
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2059,9 +2070,9 @@ 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
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2089,7 +2100,7 @@ namespace VectorTools
* point_difference() function.
*
* @note If the cell in which the point is found
- * is not locally owned, an exception of type
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2104,7 +2115,7 @@ 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
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2134,7 +2145,7 @@ namespace VectorTools
* "step-3".
*
* @note If the cell in which the point is found
- * is not locally owned, an exception of type
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2146,9 +2157,9 @@ 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
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2174,7 +2185,7 @@ namespace VectorTools
* mapping to evaluate the difference.
*
* @note If the cell in which the point is found
- * is not locally owned, an exception of type
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2190,7 +2201,7 @@ 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
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2216,7 +2227,7 @@ namespace VectorTools
* mapping to evaluate the difference.
*
* @note If the cell in which the point is found
- * is not locally owned, an exception of type
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
@@ -2231,7 +2242,7 @@ 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
+ * is not locally owned, an exception of type
* VectorTools::ExcPointNotAvailableHere
* is thrown.
*/
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 92670d6ac5..93e6847b74 100644
--- a/deal.II/include/deal.II/numerics/vector_tools.templates.h
+++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h
@@ -2404,6 +2404,16 @@ namespace VectorTools
};
+ template
+ std::ostream & operator << (std::ostream &out,
+ const VectorDoFTuple &vdt)
+ {
+ for (unsigned int d=0; d vector_dofs;
vector_dofs.dof_indices[0] = face_dofs[i];
+ Assert(first_vector_component+dim<=fe.n_components(),
+ ExcMessage("Error: the finite element does not have enough components "
+ "to define a normal direction."));
+
for (unsigned int k=0; kcenter() << std::endl
+ << " normal=" << normal_vector << std::endl;
+#endif
}
}
@@ -4405,6 +4425,18 @@ namespace VectorTools
if (p == dof_to_normals_map.end())
same_dof_range[1] = dof_to_normals_map.end();
+#ifdef DEBUG_NO_NORMAL_FLUX
+ std::cout << "For dof indices <" << p->first << ">, found the following normals"
+ << std::endl;
+ for (typename DoFToNormalsMap::const_iterator
+ q = same_dof_range[0];
+ q != same_dof_range[1]; ++q)
+ std::cout << " " << q->second.first
+ << " from cell " << q->second.second
+ << std::endl;
+#endif
+
+
// now compute the reverse mapping: for each of the cells that
// contributed to the current set of vector dofs, add up the normal
// vectors. the values of the map are pairs of normal vectors and
@@ -4438,9 +4470,18 @@ namespace VectorTools
= std::make_pair ((old_normal * old_count + q->second.first) / (old_count + 1),
old_count + 1);
}
-
Assert (cell_to_normals_map.size() >= 1, ExcInternalError());
+#ifdef DEBUG_NO_NORMAL_FLUX
+ std::cout << " cell_to_normals_map:" << std::endl;
+ for (typename CellToNormalsMap::const_iterator
+ x = cell_to_normals_map.begin();
+ x != cell_to_normals_map.end(); ++x)
+ std::cout << " " << x->first << " -> ("
+ << x->second.first << ',' << x->second.second << ')'
+ << std::endl;
+#endif
+
// count the maximum number of contributions from each cell
unsigned int max_n_contributions_per_cell = 1;
for (typename CellToNormalsMap::const_iterator
@@ -4463,6 +4504,11 @@ namespace VectorTools
// encountered this vector dof exactly once while looping over all
// their faces. as stated in the documentation, this is the case
// where we want to simply average over all normal vectors
+ //
+ // the typical case is in 2d where multiple cells meet at one
+ // vertex sitting on the boundary. same in 3d for a vertex that
+ // is associated with only one of the boundary indicators passed
+ // to this function
case 1:
{
@@ -4572,7 +4618,22 @@ namespace VectorTools
// normal vectors it has contributed. we currently only implement
// if this is dim-1 for all cells (if a single cell has
// contributed dim, or if all adjacent cells have contributed 1
- // normal vector, this is already handled above)
+ // normal vector, this is already handled above).
+ //
+ // we only implement the case that all cells contribute
+ // dim-1 because we assume that we are following an edge
+ // of the domain (think: we are looking at a vertex
+ // located on one of the edges of a refined cube where the
+ // boundary indicators of the two adjacent faces of the
+ // cube are both listed in the set of boundary indicators
+ // passed to this function). in that case, all cells along
+ // that edge of the domain are assumed to have contributed
+ // dim-1 normal vectors. however, there are cases where
+ // this assumption is not justified (see the lengthy
+ // explanation in test no_flux_12.cc) and in those cases
+ // we simply ignore the cell that contributes only
+ // once. this is also discussed at length in the
+ // documentation of this function.
//
// for each contributing cell compute the tangential vector that
// remains unconstrained
@@ -4582,7 +4643,22 @@ namespace VectorTools
contribution != cell_contributions.end();
++contribution)
{
- Assert (contribution->second.size() == dim-1, ExcNotImplemented());
+#ifdef DEBUG_NO_NORMAL_FLUX
+ std::cout << " Treating edge case with dim-1 contributions." << std::endl
+ << " Looking at cell " << contribution->first
+ << " which has contributed these normal vectors:"
+ << std::endl;
+ for (typename std::list >::const_iterator
+ t = contribution->second.begin();
+ t != contribution->second.end();
+ ++t)
+ std::cout << " " << *t << std::endl;
+#endif
+
+ // as mentioned above, simply ignore cells that only
+ // contribute once
+ if (contribution->second.size() < dim-1)
+ continue;
Tensor<1,dim> normals[dim-1];
{
@@ -4744,7 +4820,7 @@ namespace VectorTools
// can find constrained ones
unsigned int n_total_constraints_found = 0;
for (typename std::set,
- PointComparator >::const_iterator it=vector_dofs.begin();
+ PointComparator >::const_iterator it=vector_dofs.begin();
it!=vector_dofs.end(); ++it)
{
unsigned int n_constraints = 0;
diff --git a/deal.II/source/dofs/dof_handler_policy.cc b/deal.II/source/dofs/dof_handler_policy.cc
index 3db4b2104a..4e8a60fc77 100644
--- a/deal.II/source/dofs/dof_handler_policy.cc
+++ b/deal.II/source/dofs/dof_handler_policy.cc
@@ -866,8 +866,16 @@ namespace internal
PolicyBase::~PolicyBase ()
{}
-
/* --------------------- class Sequential ---------------- */
+
+ template
+ std::string
+ Sequential::to_string () const
+ {
+ return "Policy::Sequential<"+Utilities::int_to_string(dim)+
+ ","+Utilities::int_to_string(spacedim)+">";
+ }
+
template
@@ -1926,6 +1934,13 @@ namespace internal
#endif // DEAL_II_WITH_P4EST
+ template
+ std::string
+ ParallelDistributed::to_string () const
+ {
+ return "Policy::ParallelDistributed<"+Utilities::int_to_string(dim)+
+ ","+Utilities::int_to_string(spacedim)+">";
+ }
template
diff --git a/deal.II/source/dofs/dof_renumbering.inst.in b/deal.II/source/dofs/dof_renumbering.inst.in
index d35eeb5dcd..2f977adbee 100644
--- a/deal.II/source/dofs/dof_renumbering.inst.in
+++ b/deal.II/source/dofs/dof_renumbering.inst.in
@@ -17,7 +17,7 @@
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
{
-#if deal_II_dimension <= deal_II_space_dimension
+#if deal_II_dimension <= deal_II_space_dimension
namespace DoFRenumbering
\{
namespace boost
@@ -56,6 +56,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS
void subdomain_wise >
(DoFHandler &);
+ template
+ void
+ compute_subdomain_wise (std::vector &new_dof_indices,
+ const DoFHandler &dof_handler);
+
\} // namespace DoFRenumbering
#endif
}
@@ -64,7 +69,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS
// where applicable and move to codimension cases above also when applicable
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
{
-#if deal_II_dimension == deal_II_space_dimension
+#if deal_II_dimension == deal_II_space_dimension
namespace DoFRenumbering
\{
namespace boost
@@ -125,10 +130,10 @@ namespace DoFRenumbering
compute_minimum_degree (std::vector &,
const hp::DoFHandler &, bool, bool);
\}
-
-
+
+
// non-boost functions:
-
+
template
void Cuthill_McKee >
(hp::DoFHandler&,
@@ -281,7 +286,7 @@ namespace DoFRenumbering
template
void downstream
- (DoFHandler&,
+ (DoFHandler&,
const unsigned int,
const Point&,
const bool);
@@ -352,7 +357,7 @@ namespace DoFRenumbering
const hp::DoFHandler &,
const std::vector &);
-
+
template
void subdomain_wise >
(hp::DoFHandler &);
@@ -374,7 +379,3 @@ namespace DoFRenumbering
\} // namespace DoFRenumbering
#endif
}
-
-
-
-
diff --git a/deal.II/source/fe/fe_tools_interpolate.cc b/deal.II/source/fe/fe_tools_interpolate.cc
index 8eef0c9e12..a1983255a1 100644
--- a/deal.II/source/fe/fe_tools_interpolate.cc
+++ b/deal.II/source/fe/fe_tools_interpolate.cc
@@ -93,15 +93,15 @@ namespace FETools
// if u1 is a parallel distributed
// PETSc vector, we check the local
// size of u1 for safety
- Assert(dynamic_cast(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.n_locally_owned_dofs()));
};
if (dynamic_cast(&u2) != 0)
if (dynamic_cast*>(&dof2) != 0)
{
- Assert(dynamic_cast(&u2)->local_size() == dof2.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u2)->local_size(), dof2.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u2)->local_size() == dof2.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u2)->local_size(), dof2.n_locally_owned_dofs()));
};
#endif
// allocate vectors at maximal
@@ -254,15 +254,15 @@ namespace FETools
// if u1 is a parallel distributed
// PETSc vector, we check the local
// size of u1 for safety
- Assert(dynamic_cast(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.n_locally_owned_dofs()));
};
if (dynamic_cast(&u1_interpolated) != 0)
if (dynamic_cast*>(&dof1) != 0)
{
- Assert(dynamic_cast(&u1_interpolated)->local_size() == dof1.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u1_interpolated)->local_size(), dof1.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u1_interpolated)->local_size() == dof1.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u1_interpolated)->local_size(), dof1.n_locally_owned_dofs()));
};
#endif
@@ -334,15 +334,15 @@ namespace FETools
// if u1 is a parallel distributed
// PETSc vector, we check the local
// size of u1 for safety
- Assert(dynamic_cast(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.n_locally_owned_dofs()));
};
if (dynamic_cast(&u1_interpolated) != 0)
if (dynamic_cast*>(&dof1) != 0)
{
- Assert(dynamic_cast(&u1_interpolated)->local_size() == dof1.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u1_interpolated)->local_size(), dof1.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u1_interpolated)->local_size() == dof1.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u1_interpolated)->local_size(), dof1.n_locally_owned_dofs()));
};
#endif
@@ -449,12 +449,41 @@ namespace FETools
DoFTools::extract_locally_relevant_dofs (dof2,
dof2_locally_relevant_dofs);
- PETScWrappers::MPI::Vector u2_out (u1.get_mpi_communicator(),
- dof2_locally_owned_dofs);
+ PETScWrappers::MPI::Vector u2_out (dof2_locally_owned_dofs,
+ u1.get_mpi_communicator());
interpolate(dof1, u1, dof2, constraints2, u2_out);
- PETScWrappers::MPI::Vector u2 (u1.get_mpi_communicator(),
- dof2_locally_owned_dofs,
- dof2_locally_relevant_dofs);
+ PETScWrappers::MPI::Vector u2 (dof2_locally_owned_dofs,
+ dof2_locally_relevant_dofs,
+ u1.get_mpi_communicator());
+ u2 = u2_out;
+ interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
+ }
+#endif
+
+ // special version for Trilinos
+#ifdef DEAL_II_USE_TRILINOS
+ template
+ void back_interpolate (const DoFHandler &dof1,
+ const ConstraintMatrix &constraints1,
+ const TrilinosWrappers::MPI::Vector &u1,
+ const DoFHandler &dof2,
+ const ConstraintMatrix &constraints2,
+ TrilinosWrappers::MPI::Vector &u1_interpolated)
+ {
+ // if u1 is a parallel distributed Trilinos vector, we create a
+ // vector u2 with based on the sets of locally owned and relevant
+ // dofs of dof2
+ IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs();
+ IndexSet dof2_locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs (dof2,
+ dof2_locally_relevant_dofs);
+
+ TrilinosWrappers::MPI::Vector u2_out (dof2_locally_owned_dofs,
+ u1.get_mpi_communicator());
+ interpolate(dof1, u1, dof2, constraints2, u2_out);
+ TrilinosWrappers::MPI::Vector u2 (dof2_locally_owned_dofs,
+ dof2_locally_relevant_dofs,
+ u1.get_mpi_communicator());
u2 = u2_out;
interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
}
@@ -537,15 +566,15 @@ namespace FETools
// if u1 is a parallel distributed
// PETSc vector, we check the local
// size of u1 for safety
- Assert(dynamic_cast(&u1)->local_size() == dof1.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u1)->local_size() == dof1.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u1)->local_size(), dof1.n_locally_owned_dofs()));
};
if (dynamic_cast(&u1_difference) != 0)
if (dynamic_cast*>(&dof1) != 0)
{
- Assert(dynamic_cast(&u1_difference)->local_size() == dof1.locally_owned_dofs().n_elements(),
- ExcDimensionMismatch(dynamic_cast(&u1_difference)->local_size(), dof1.locally_owned_dofs().n_elements()));
+ Assert(dynamic_cast(&u1_difference)->local_size() == dof1.n_locally_owned_dofs(),
+ ExcDimensionMismatch(dynamic_cast(&u1_difference)->local_size(), dof1.n_locally_owned_dofs()));
};
#endif
@@ -598,6 +627,50 @@ namespace FETools
+ namespace internal
+ {
+ namespace
+ {
+ template
+ void interpolation_difference (const DoFHandler &dof1,
+ const ConstraintMatrix &constraints1,
+ const InVector &u1,
+ const DoFHandler &dof2,
+ const ConstraintMatrix &constraints2,
+ OutVector &u1_difference)
+ {
+ back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
+ u1_difference.sadd(-1, u1);
+ }
+
+ // special version for Trilinos
+#ifdef DEAL_II_USE_TRILINOS
+ template
+ void interpolation_difference (const DoFHandler &dof1,
+ const ConstraintMatrix &constraints1,
+ const TrilinosWrappers::MPI::Vector &u1,
+ const DoFHandler &dof2,
+ const ConstraintMatrix &constraints2,
+ TrilinosWrappers::MPI::Vector &u1_difference)
+ {
+ back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
+
+ // Trilinos vectors with and without ghost entries are very different
+ // and we cannot use the sadd function directly, so we have to create
+ // a completely distributed vector first and copy the local entries
+ // from the vector with ghost entries
+ TrilinosWrappers::MPI::Vector u1_completely_distributed (u1_difference.vector_partitioner ());
+
+ u1_completely_distributed = u1;
+
+ u1_difference.sadd(-1, u1_completely_distributed);
+ }
+#endif
+ }
+ }
+
+
+
template
void interpolation_difference(const DoFHandler &dof1,
const ConstraintMatrix &constraints1,
@@ -615,8 +688,7 @@ namespace FETools
interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference);
else
{
- back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
- u1_difference.sadd(-1, u1);
+ internal::interpolation_difference(dof1, constraints1, u1, dof2, constraints2, u1_difference);
}
}
diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc
index 530a199602..7c708daea1 100644
--- a/deal.II/source/fe/mapping_q1.cc
+++ b/deal.II/source/fe/mapping_q1.cc
@@ -1684,33 +1684,69 @@ transform_real_to_unit_cell_internal
// Newton iteration to solve
- // f(x)=p(x)-p=0
- // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+ // f(x)=p(x)-p=0
+ // where we are looking for 'x' and p(x) is the forward transformation
+ // from unit to real cell. We solve this using a Newton iteration
+ // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+ // The start value is set to be the linear approximation to the cell
- // The start value was set to be the
- // linear approximation to the cell
-
- // The shape values and derivatives
- // of the mapping at this point are
+ // The shape values and derivatives of the mapping at this point are
// previously computed.
- // For the <2,3> case there is a
- // template specialization.
-
- // f(x)
-
Point p_unit = initial_p_unit;
compute_shapes(std::vector > (1, p_unit), mdata);
- Point p_real(transform_unit_to_real_cell_internal(mdata));
+ Point p_real = transform_unit_to_real_cell_internal(mdata);
Point f = p_real-p;
- const double eps = 1.e-12*cell->diameter();
+ // early out if we already have our point
+ if (f.square() < 1e-24 * cell->diameter() * cell->diameter())
+ return p_unit;
+
+ // we need to compare the position of the computed p(x) against the given
+ // point 'p'. We will terminate the iteration and return 'x' if they are
+ // less than eps apart. The question is how to choose eps -- or, put maybe
+ // more generally: in which norm we want these 'p' and 'p(x)' to be eps
+ // apart.
+ //
+ // the question is difficult since we may have to deal with very elongated
+ // cells where we may achieve 1e-12*h for the distance of these two points
+ // in the 'long' direction, but achieving this tolerance in the 'short'
+ // direction of the cell may not be possible
+ //
+ // what we do instead is then to terminate iterations if
+ // \| p(x) - p \|_A < eps
+ // where the A-norm is somehow induced by the transformation of the cell.
+ // in particular, we want to measure distances relative to the sizes of
+ // the cell in its principal directions.
+ //
+ // to define what exactly A should be, note that to first order we have
+ // the following (assuming that x* is the solution of the problem, i.e.,
+ // p(x*)=p):
+ // p(x) - p = p(x) - p(x*)
+ // = -grad p(x) * (x*-x) + higher order terms
+ // This suggest to measure with a norm that corresponds to
+ // A = {[grad p(x]^T [grad p(x)]}^{-1}
+ // because then
+ // \| p(x) - p \|_A \approx \| x - x* \|
+ // Consequently, we will try to enforce that
+ // \| p(x) - p \|_A = \| f \| <= eps
+ //
+ // Note that using this norm is a bit dangerous since the norm changes
+ // in every iteration (A isn't fixed by depends on xk). However, if the
+ // cell is not too deformed (it may be stretched, but not twisted) then
+ // the mapping is almost linear and A is indeed constant or nearly so.
+ const double eps = 1.e-11;
const unsigned int newton_iteration_limit = 20;
unsigned int newton_iteration=0;
- while (f.square()>eps*eps)
+ double last_f_weighted_norm = -1./0.;
+ do
{
+#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
+ std::cout << "Newton iteration " << newton_iteration << std::endl;
+#endif
+
// f'(x)
Tensor<2,spacedim> df;
for (unsigned int k=0; k delta;
- contract (delta, invert(df), static_cast&>(f));
+ Tensor<2,spacedim> df_inverse = invert(df);
+ contract (delta, df_inverse, static_cast&>(f));
+#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
+ std::cout << " delta=" << delta << std::endl;
+#endif
+
// do a line search
double step_length = 1;
do
{
- // update of p_unit. The
- // spacedimth component of
- // transformed point is simply
- // ignored in codimension one
- // case. When this component is
- // not zero, then we are
- // projecting the point to the
- // surface or curve identified
- // by the cell.
+ // update of p_unit. The spacedim-th component of transformed point
+ // is simply ignored in codimension one case. When this component is
+ // not zero, then we are projecting the point to the surface or
+ // curve identified by the cell.
Point p_unit_trial = p_unit;
for (unsigned int i=0; i p_real_trial = transform_unit_to_real_cell_internal(mdata);
const Point f_trial = p_real_trial-p;
+#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
+ std::cout << " step_length=" << step_length << std::endl
+ << " ||f || =" << f.norm() << std::endl
+ << " ||f*|| =" << f_trial.norm() << std::endl
+ << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl;
+#endif
+
// see if we are making progress with the current step length
// and if not, reduce it by a factor of two and try again
+ //
+ // strictly speaking, we should probably use the same norm as we use
+ // for the outer algorithm. in practice, line search is just a
+ // crutch to find a "reasonable" step length, and so using the l2
+ // norm is probably just fine
if (f_trial.norm() < f.norm())
{
p_real = p_real_trial;
@@ -1764,27 +1812,20 @@ transform_real_to_unit_cell_internal
else if (step_length > 0.05)
step_length /= 2;
else
- goto failure;
+ AssertThrow (false,
+ (typename Mapping::ExcTransformationFailed()));
}
while (true);
++newton_iteration;
if (newton_iteration > newton_iteration_limit)
- goto failure;
+ AssertThrow (false,
+ (typename Mapping::ExcTransformationFailed()));
+ last_f_weighted_norm = (df_inverse * f).norm();
}
-
+ while (last_f_weighted_norm > eps);
+
return p_unit;
-
- // if we get to the following label, then we have either run out
- // of Newton iterations, or the line search has not converged.
- // in either case, we need to give up, so throw an exception that
- // can then be caught
-failure:
- AssertThrow (false, (typename Mapping::ExcTransformationFailed()));
-
- // ...the compiler wants us to return something, though we can
- // of course never get here...
- return Point();
}
diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc
index 298de0aac3..b129a9d711 100644
--- a/deal.II/source/grid/grid_generator.cc
+++ b/deal.II/source/grid/grid_generator.cc
@@ -629,11 +629,30 @@ namespace GridGenerator
}
+
+ template<>
+ void
+ parallelogram (Triangulation<3> &tria,
+ const Point<3> corners[3],
+ const bool colorize)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+ template<>
+ void
+ parallelogram (Triangulation<1> &tria,
+ const Point<1> corners[1],
+ const bool colorize)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
// Implementation for 2D only
template<>
void
parallelogram (Triangulation<2> &tria,
- const Point<2> (&corners)[2],
+ const Point<2> corners[2],
const bool colorize)
{
std::vector > vertices (GeometryInfo<2>::vertices_per_cell);
@@ -655,6 +674,25 @@ namespace GridGenerator
}
+ template<>
+ void
+ parallelogram (Triangulation<1> &tria,
+ const Tensor<2,1> &corners,
+ const bool colorize)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+ template<>
+ void
+ parallelogram (Triangulation<3> &tria,
+ const Tensor<2,3> &corners,
+ const bool colorize)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+
template<>
void
parallelogram (Triangulation<2> &tria,
diff --git a/deal.II/source/grid/grid_generator.inst.in b/deal.II/source/grid/grid_generator.inst.in
index 5ce5fcc0f4..ff179bd53e 100644
--- a/deal.II/source/grid/grid_generator.inst.in
+++ b/deal.II/source/grid/grid_generator.inst.in
@@ -50,6 +50,17 @@ namespace GridGenerator \{
Triangulation &,
const unsigned int, const double, const double);
+ template
+ void
+ parallelogram (Triangulation &,
+ const Tensor<2,deal_II_dimension> &,
+ const bool );
+ template
+ void
+ parallelogram (Triangulation &,
+ const Point [deal_II_dimension],
+ const bool);
+
template void
subdivided_hyper_rectangle
(Triangulation &,
diff --git a/deal.II/source/grid/grid_out.cc b/deal.II/source/grid/grid_out.cc
index 5a204ef522..c110dcfdc8 100644
--- a/deal.II/source/grid/grid_out.cc
+++ b/deal.II/source/grid/grid_out.cc
@@ -368,6 +368,7 @@ namespace GridOutFlags
label_cell_index(label_cell_index),
label_material_id(label_material_id),
label_subdomain_id(label_subdomain_id),
+ label_level_subdomain_id(false),
draw_colorbar(draw_colorbar),
draw_legend(draw_legend)
{}
diff --git a/deal.II/source/lac/chunk_sparsity_pattern.cc b/deal.II/source/lac/chunk_sparsity_pattern.cc
index 57aaa9cf6c..4a2d34d402 100644
--- a/deal.II/source/lac/chunk_sparsity_pattern.cc
+++ b/deal.II/source/lac/chunk_sparsity_pattern.cc
@@ -755,6 +755,15 @@ void ChunkSparsityPattern::copy_from (const Com
const size_type ,
const bool);
template
+void ChunkSparsityPattern::copy_from (const CompressedSparsityPattern &,
+ const size_type);
+template
+void ChunkSparsityPattern::copy_from (const CompressedSetSparsityPattern &,
+ const size_type);
+template
+void ChunkSparsityPattern::copy_from (const CompressedSimpleSparsityPattern &,
+ const size_type);
+template
void ChunkSparsityPattern::create_from
(const unsigned int,
const unsigned int,
diff --git a/deal.II/source/numerics/vector_tools_interpolate.inst.in b/deal.II/source/numerics/vector_tools_interpolate.inst.in
index f877d0f4f0..e44dd85c4e 100644
--- a/deal.II/source/numerics/vector_tools_interpolate.inst.in
+++ b/deal.II/source/numerics/vector_tools_interpolate.inst.in
@@ -22,7 +22,7 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
template
void interpolate
- (const Mapping&,
+ (const Mapping::dimension,DoFHandler::space_dimension>&,
const DoFHandler&,
const Function&,
VEC&);
@@ -120,5 +120,3 @@ for (deal_II_dimension : DIMENSIONS)
\}
}
-
-
diff --git a/deal.II/tests/CMakeLists.txt b/deal.II/tests/CMakeLists.txt
index d594c3c3d4..74a5cfa2f5 100644
--- a/deal.II/tests/CMakeLists.txt
+++ b/deal.II/tests/CMakeLists.txt
@@ -23,6 +23,8 @@
# of complexity :-(
#
+MESSAGE(STATUS "")
+MESSAGE(STATUS "Setup testsuite with TEST_DIR ${TEST_DIR}")
ADD_SUBDIRECTORY(quick_tests)
#
@@ -57,6 +59,7 @@ SET(_options)
LIST(APPEND _options -DDEAL_II_SOURCE_DIR=${CMAKE_SOURCE_DIR})
LIST(APPEND _options -DDEAL_II_BINARY_DIR=${CMAKE_BINARY_DIR})
FOREACH(_var
+ DIFF_DIR
NUMDIFF_DIR
TEST_DIFF
TEST_OVERRIDE_LOCATION
@@ -70,6 +73,13 @@ FOREACH(_var
ENDIF()
ENDFOREACH()
+# Be as quiet as possible:
+IF(CMAKE_SYSTEM_NAME MATCHES "Windows")
+ SET(_shoo)
+ELSE()
+ SET(_shoo >/dev/null)
+ENDIF()
+
#
# Glob together a list of all subfolders to set up:
#
@@ -99,7 +109,7 @@ FOREACH(_category ${_categories})
${CMAKE_CURRENT_BINARY_DIR}/${_category}
COMMAND cd ${CMAKE_CURRENT_BINARY_DIR}/${_category} &&
${CMAKE_COMMAND} -G${CMAKE_GENERATOR} ${_options} ${_category_dir}
- > /dev/null # Shoo!
+ ${_shoo}
DEPENDS ${_category_dir}
COMMENT "Processing tests/${_category}"
)
@@ -111,14 +121,14 @@ FOREACH(_category ${_categories})
ADD_CUSTOM_TARGET(regen_tests_${_category}
COMMAND [ ! -d ${_category} ] || ${CMAKE_COMMAND}
--build ${CMAKE_CURRENT_BINARY_DIR}/${_category} --target regenerate
- > /dev/null # Shoo!
+ ${_shoo}
)
ADD_DEPENDENCIES(regen_tests regen_tests_${_category})
ADD_CUSTOM_TARGET(clean_tests_${_category}
COMMAND [ ! -d ${_category} ] || ${CMAKE_COMMAND}
--build ${CMAKE_CURRENT_BINARY_DIR}/${_category} --target clean
- > /dev/null # Shoo!
+ ${_shoo}
)
ADD_DEPENDENCIES(clean_tests clean_tests_${_category})
@@ -152,3 +162,4 @@ FOREACH(_category ${_categories})
ENDFOREACH()
MESSAGE(STATUS "Regenerating testsuite subprojects - Done")
+MESSAGE(STATUS "")
diff --git a/deal.II/tests/quick_tests/CMakeLists.txt b/deal.II/tests/quick_tests/CMakeLists.txt
index 727911cb15..a01e681928 100644
--- a/deal.II/tests/quick_tests/CMakeLists.txt
+++ b/deal.II/tests/quick_tests/CMakeLists.txt
@@ -21,7 +21,6 @@ ENABLE_TESTING()
# Use the first available build type (this prefers debug mode if available):
LIST(GET DEAL_II_BUILD_TYPES 0 _mybuild)
-MESSAGE(STATUS "")
MESSAGE(STATUS "Setup quick_tests in ${_mybuild} mode")
# define a macro to set up a quick test:
@@ -45,10 +44,10 @@ MACRO(make_quicktest test_basename build_name mpi_run)
)
ENDIF()
- ADD_CUSTOM_TARGET(${_target}.run
- DEPENDS ${_target}
- DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/${_target}-OK
- COMMAND echo "${_target}: PASSED.")
+ ADD_CUSTOM_TARGET(${_target}.run
+ DEPENDS ${_target}
+ DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/${_target}-OK
+ COMMAND echo "${_target}: PASSED.")
ADD_TEST(NAME ${_target}
COMMAND ${CMAKE_COMMAND} -DTRGT=${_target}.run -DTEST=${_target}
diff --git a/deal.II/tests/quick_tests/step-petsc.cc b/deal.II/tests/quick_tests/step-petsc.cc
index 9290f55a11..21586bf98a 100644
--- a/deal.II/tests/quick_tests/step-petsc.cc
+++ b/deal.II/tests/quick_tests/step-petsc.cc
@@ -1,5 +1,5 @@
/* ---------------------------------------------------------------------
- * $Id: step-petsc.cc $
+ * $Id$
*
* Copyright (C) 2013 by the deal.II authors
*
diff --git a/deal.II/tests/quick_tests/tbb.cc b/deal.II/tests/quick_tests/tbb.cc
index a30467bec6..654da479bf 100644
--- a/deal.II/tests/quick_tests/tbb.cc
+++ b/deal.II/tests/quick_tests/tbb.cc
@@ -1,5 +1,5 @@
// ---------------------------------------------------------------------
-// $Id: affinity.cc 31527 2013-11-03 09:58:45Z maier $
+// $Id$
//
// Copyright (C) 2013 by the deal.II authors
//
--
2.39.5