From c337e9dc44f8026b3b54a6882e66d783cfd5cfb3 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Tue, 14 Jan 2020 11:29:15 +0100 Subject: [PATCH] make_flux_sparsity_pattern with flux couplings on a subset of the faces Add an optional std::function input parameter to make_flux_sparsity_pattern which describe over which faces there should be flux couplings in the created sparsity pattern. internal::make_flux_sparsity_pattern has one implementation for a generic DoFHandlerType and a specialization for hp::DoFHandler. To circument a bug where MSVC can't deduce which of these internal functions should be called, change DoFHandlerType to dealii::DoFHandler in the generic function. To avoid another bug where MSVC can't construct a templated class, call one of the overloaded versions of make_flux_sparsity_pattern with always_couple_on_faces (even if this is already the default value). --- doc/news/changes/minor/20200128SimonSticko | 5 + include/deal.II/dofs/dof_tools.h | 58 ++++++++- source/dofs/dof_tools_sparsity.cc | 113 +++++++++------- source/dofs/dof_tools_sparsity.inst.in | 10 +- tests/dofs/dof_tools_25.cc | 143 +++++++++++++++++++++ tests/dofs/dof_tools_25.output | 12 ++ 6 files changed, 288 insertions(+), 53 deletions(-) create mode 100644 doc/news/changes/minor/20200128SimonSticko create mode 100644 tests/dofs/dof_tools_25.cc create mode 100644 tests/dofs/dof_tools_25.output diff --git a/doc/news/changes/minor/20200128SimonSticko b/doc/news/changes/minor/20200128SimonSticko new file mode 100644 index 0000000000..9081573487 --- /dev/null +++ b/doc/news/changes/minor/20200128SimonSticko @@ -0,0 +1,5 @@ +Add an optional std::function parameter to make_flux_sparsity_pattern which can +be used to specify over which faces there should be a flux coupling in the +created sparsity pattern. +
+(Simon Sticko, 2020/01/28) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index b5a51ae3cf..f83c001f6a 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -68,6 +68,28 @@ namespace GridTools template struct PeriodicFacePair; } + +namespace DoFTools +{ + namespace internal + { + /* + * Default value of the face_has_flux_coupling parameter of + * make_flux_sparsity_pattern. Defined here (instead of using a default + * lambda in the parameter list) to avoid a bug in gcc where the same lambda + * gets defined multiple times. + */ + template + inline bool + always_couple_on_faces( + const typename DoFHandlerType::active_cell_iterator &, + const unsigned int) + { + return true; + } + } // namespace internal +} // namespace DoFTools + #endif /** @@ -775,18 +797,40 @@ namespace DoFTools * components of a finite element are continuous and some discontinuous, * allowing constraints to be imposed on the continuous part while also * building the flux terms needed for the discontinuous part. + * + * The optional @param face_has_flux_coupling can be used to specify on which + * faces flux couplings occur. This allows for creating a sparser pattern when + * using a bilinear form where flux terms only appear on a subset of the faces + * in the triangulation. By default flux couplings are added over all internal + * faces. @param face_has_flux_coupling should be a function that takes an + * active_cell_iterator and a face index and should return true if there is a + * flux coupling over the face. When using the ::dealii::DoFHandler we could, + * for example, use + * + * @code + * auto face_has_flux_coupling = + * [](const typename DoFHandler::active_cell_iterator &cell, + * const unsigned int face_index) { + * const Point &face_center = cell->face(face_index)->center(); + * return 0 < face_center[0]; + * }; */ template void - make_flux_sparsity_pattern(const DoFHandlerType & dof, - SparsityPatternType & sparsity, - const AffineConstraints &constraints, - const bool keep_constrained_dofs, - const Table<2, Coupling> &couplings, - const Table<2, Coupling> &face_couplings, - const types::subdomain_id subdomain_id); + make_flux_sparsity_pattern( + const DoFHandlerType & dof, + SparsityPatternType & sparsity, + const AffineConstraints &constraints, + const bool keep_constrained_dofs, + const Table<2, Coupling> & couplings, + const Table<2, Coupling> & face_couplings, + const types::subdomain_id subdomain_id, + const std::function< + bool(const typename DoFHandlerType::active_cell_iterator &, + const unsigned int)> &face_has_flux_coupling = + &internal::always_couple_on_faces); /** * Create the sparsity pattern for boundary matrices. See the general diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc index b3c05e9e93..a4c3ddc967 100644 --- a/source/dofs/dof_tools_sparsity.cc +++ b/source/dofs/dof_tools_sparsity.cc @@ -759,20 +759,24 @@ namespace DoFTools { // implementation of the same function in namespace DoFTools for // non-hp DoFHandlers - template void - make_flux_sparsity_pattern(const DoFHandlerType & dof, - SparsityPatternType & sparsity, - const AffineConstraints &constraints, - const bool keep_constrained_dofs, - const Table<2, Coupling> &int_mask, - const Table<2, Coupling> &flux_mask, - const types::subdomain_id subdomain_id) + make_flux_sparsity_pattern( + const DoFHandler &dof, + SparsityPatternType & sparsity, + const AffineConstraints &constraints, + const bool keep_constrained_dofs, + const Table<2, Coupling> & int_mask, + const Table<2, Coupling> & flux_mask, + const types::subdomain_id subdomain_id, + const std::function< + bool(const typename DoFHandler::active_cell_iterator &, + const unsigned int)> &face_has_flux_coupling) { - const FiniteElement &fe = dof.get_fe(); + const FiniteElement &fe = dof.get_fe(); std::vector dofs_on_this_cell( fe.dofs_per_cell); @@ -783,13 +787,10 @@ namespace DoFTools int_dof_mask = dof_couplings_from_component_couplings(fe, int_mask), flux_dof_mask = dof_couplings_from_component_couplings(fe, flux_mask); - Table<2, bool> support_on_face( - fe.dofs_per_cell, - GeometryInfo::faces_per_cell); + Table<2, bool> support_on_face(fe.dofs_per_cell, + GeometryInfo::faces_per_cell); for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) - for (unsigned int f = 0; - f < GeometryInfo::faces_per_cell; - ++f) + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) support_on_face(i, f) = fe.has_support_on_face(i, f); // Convert the int_dof_mask to bool_int_dof_mask so we can pass it @@ -801,8 +802,9 @@ namespace DoFTools if (int_dof_mask(i, j) != none) bool_int_dof_mask(i, j) = true; - typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); for (; cell != endc; ++cell) if (((subdomain_id == numbers::invalid_subdomain_id) || (subdomain_id == cell->subdomain_id())) && @@ -816,12 +818,11 @@ namespace DoFTools bool_int_dof_mask); // Loop over all interior neighbors for (unsigned int face_n = 0; - face_n < - GeometryInfo::faces_per_cell; + face_n < GeometryInfo::faces_per_cell; ++face_n) { - const typename DoFHandlerType::face_iterator cell_face = - cell->face(face_n); + const typename DoFHandler::face_iterator + cell_face = cell->face(face_n); const bool periodic_neighbor = cell->has_periodic_neighbor(face_n); @@ -846,8 +847,11 @@ namespace DoFTools } else { - typename DoFHandlerType::level_cell_iterator neighbor = - cell->neighbor_or_periodic_neighbor(face_n); + if (!face_has_flux_coupling(cell, face_n)) + continue; + + typename DoFHandler::level_cell_iterator + neighbor = cell->neighbor_or_periodic_neighbor(face_n); // If the cells are on the same level (and both are // active, locally-owned cells) then only add to the // sparsity pattern if the current cell is 'greater' in @@ -886,7 +890,7 @@ namespace DoFTools // if (neighbor->has_children()) section below. We need to // do this since we otherwise iterate over the children of // the face, which are always 0 in 1D. - if (DoFHandlerType::dimension == 1) + if (dim == 1) while (neighbor->has_children()) neighbor = neighbor->child(face_n == 0 ? 1 : 0); @@ -896,8 +900,8 @@ namespace DoFTools sub_nr != cell_face->n_children(); ++sub_nr) { - const typename DoFHandlerType::level_cell_iterator - sub_neighbor = + const typename DoFHandler:: + level_cell_iterator sub_neighbor = periodic_neighbor ? cell->periodic_neighbor_child_on_subface( face_n, sub_nr) : @@ -1069,7 +1073,11 @@ namespace DoFTools const bool keep_constrained_dofs, const Table<2, Coupling> & int_mask, const Table<2, Coupling> & flux_mask, - const types::subdomain_id subdomain_id) + const types::subdomain_id subdomain_id, + const std::function< + bool(const typename dealii::hp::DoFHandler:: + active_cell_iterator &, + const unsigned int)> &face_has_flux_coupling) { // while the implementation above is quite optimized and caches a // lot of data (see e.g. the int/flux_dof_mask tables), this is no @@ -1156,6 +1164,9 @@ namespace DoFTools level_cell_iterator neighbor = cell->neighbor_or_periodic_neighbor(face); + if (!face_has_flux_coupling(cell, face)) + continue; + // Like the non-hp case: If the cells are on the same // level (and both are active, locally-owned cells) then // only add to the sparsity pattern if the current cell @@ -1336,26 +1347,34 @@ namespace DoFTools const bool keep_constrained_dofs = true; - make_flux_sparsity_pattern(dof, - sparsity, - dummy, - keep_constrained_dofs, - int_mask, - flux_mask, - subdomain_id); + make_flux_sparsity_pattern( + dof, + sparsity, + dummy, + keep_constrained_dofs, + int_mask, + flux_mask, + subdomain_id, + internal::always_couple_on_faces); } + + template void - make_flux_sparsity_pattern(const DoFHandlerType & dof, - SparsityPatternType & sparsity, - const AffineConstraints &constraints, - const bool keep_constrained_dofs, - const Table<2, Coupling> &int_mask, - const Table<2, Coupling> &flux_mask, - const types::subdomain_id subdomain_id) + make_flux_sparsity_pattern( + const DoFHandlerType & dof, + SparsityPatternType & sparsity, + const AffineConstraints &constraints, + const bool keep_constrained_dofs, + const Table<2, Coupling> & int_mask, + const Table<2, Coupling> & flux_mask, + const types::subdomain_id subdomain_id, + const std::function< + bool(const typename DoFHandlerType::active_cell_iterator &, + const unsigned int)> &face_has_flux_coupling) { // do the error checking and frame code here, and then pass on to more // specialized functions in the internal namespace @@ -1390,16 +1409,22 @@ namespace DoFTools "associated DoF handler objects, asking for any subdomain other " "than the locally owned one does not make sense.")); + Assert( + face_has_flux_coupling, + ExcMessage( + "The function which specifies if a flux coupling occurs over a given " + "face is empty.")); + internal::make_flux_sparsity_pattern(dof, sparsity, constraints, keep_constrained_dofs, int_mask, flux_mask, - subdomain_id); + subdomain_id, + face_has_flux_coupling); } - } // end of namespace DoFTools diff --git a/source/dofs/dof_tools_sparsity.inst.in b/source/dofs/dof_tools_sparsity.inst.in index f73035578f..cce112d35e 100644 --- a/source/dofs/dof_tools_sparsity.inst.in +++ b/source/dofs/dof_tools_sparsity.inst.in @@ -194,7 +194,10 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS; const bool, const Table<2, Coupling> &, const Table<2, Coupling> &, - const types::subdomain_id); + const types::subdomain_id, + const std::function::active_cell_iterator &, + const unsigned int)> &); template void DoFTools::make_flux_sparsity_pattern, SP>( @@ -220,7 +223,10 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS; const bool, const Table<2, Coupling> &, const Table<2, Coupling> &, - const types::subdomain_id); + const types::subdomain_id, + const std::function::active_cell_iterator &, + const unsigned int)> &); #if deal_II_dimension < 3 diff --git a/tests/dofs/dof_tools_25.cc b/tests/dofs/dof_tools_25.cc new file mode 100644 index 0000000000..7f67b758c2 --- /dev/null +++ b/tests/dofs/dof_tools_25.cc @@ -0,0 +1,143 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test the version of make_flux_sparsity_pattern that takes an additional +// std::function describing over which faces flux couplings occur. +// +// Set up a 1D-triangulation with 3 cells (left, middle and right) where the +// face between left and middle is located at x=0: +// +// |-L-|-M-|-R-| +// x=0 +// +// Distribute FE_Q<1>(1) elements over this grid. Then call +// make_flux_sparsity_pattern with an std::function specifying that there should +// only be a flux coupling between the left and middle cell. Verify that this +// is the case for the created sparsity pattern. + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + +// Call the make_flux_sparsity_pattern with a lambda specifying that only the +// face at x = 0 should have a flux coupling. Print the constructed sparsity +// pattern to deallog. +template +void +create_and_print_pattern(const DoFHandlerType &dof_handler) +{ + DynamicSparsityPattern dynamic_pattern(dof_handler.n_dofs()); + + auto face_has_flux_coupling = + [](const typename DoFHandlerType::active_cell_iterator &cell, + const unsigned int face_index) { + // Only add a flux coupling if the face is at x = 0. + const Point ¢er = cell->face(face_index)->center(); + return std::abs(center[0]) < 1e-3; + return true; + }; + + AffineConstraints constraints; + constraints.close(); + + const bool keep_constrained_dofs = true; + + // Set all couplings to be connected. + const unsigned int n_components = 1; + Table<2, DoFTools::Coupling> couplings(n_components, n_components); + couplings[0][0] = DoFTools::Coupling::always; + Table<2, DoFTools::Coupling> face_couplings(n_components, n_components); + face_couplings[0][0] = DoFTools::Coupling::always; + + const types::subdomain_id subdomain_id = 0; + + DoFTools::make_flux_sparsity_pattern(dof_handler, + dynamic_pattern, + constraints, + keep_constrained_dofs, + couplings, + face_couplings, + subdomain_id, + face_has_flux_coupling); + + + dynamic_pattern.print(deallog.get_file_stream()); +} + + + +// The implementation is different depending on which DoFHandler is used. Test +// both of them. +template +void +test_with_both_dof_handlers(const Triangulation &triangulation) +{ + const FE_Q element(1); + + deallog << "dealii::DoFHandler" << std::endl; + { + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(element); + create_and_print_pattern(dof_handler); + } + + deallog << std::endl; + + deallog << "hp::DoFHandler" << std::endl; + { + hp::FECollection fe_collection(element); + hp::DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe_collection); + create_and_print_pattern(dof_handler); + } +} + + + +// Create the 3-cell-triangulation described at the top. +void create_3_cell_triangulation(Triangulation<1> &triangulation) +{ + const unsigned int n_elements = 3; + const double left = -1; + const double right = 2; + + GridGenerator::subdivided_hyper_cube(triangulation, n_elements, left, right); +} + + + +int +main() +{ + initlog(); + + const int dim = 1; + Triangulation triangulation; + create_3_cell_triangulation(triangulation); + + test_with_both_dof_handlers(triangulation); + return 0; +} diff --git a/tests/dofs/dof_tools_25.output b/tests/dofs/dof_tools_25.output new file mode 100644 index 0000000000..c2353b3664 --- /dev/null +++ b/tests/dofs/dof_tools_25.output @@ -0,0 +1,12 @@ + +DEAL::dealii::DoFHandler +[0,0,1,2] +[1,0,1,2] +[2,0,1,2,3] +[3,2,3] +DEAL:: +DEAL::hp::DoFHandler +[0,0,1,2] +[1,0,1,2] +[2,0,1,2,3] +[3,2,3] -- 2.39.5