--- /dev/null
+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.
+<br>
+(Simon Sticko, 2020/01/28)
template <typename CellIterator>
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 <typename DoFHandlerType>
+ inline bool
+ always_couple_on_faces(
+ const typename DoFHandlerType::active_cell_iterator &,
+ const unsigned int)
+ {
+ return true;
+ }
+ } // namespace internal
+} // namespace DoFTools
+
#endif
/**
* 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<dim>::active_cell_iterator &cell,
+ * const unsigned int face_index) {
+ * const Point<dim> &face_center = cell->face(face_index)->center();
+ * return 0 < face_center[0];
+ * };
*/
template <typename DoFHandlerType,
typename SparsityPatternType,
typename number>
void
- make_flux_sparsity_pattern(const DoFHandlerType & dof,
- SparsityPatternType & sparsity,
- const AffineConstraints<number> &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<number> &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<DoFHandlerType>);
/**
* Create the sparsity pattern for boundary matrices. See the general
{
// implementation of the same function in namespace DoFTools for
// non-hp DoFHandlers
- template <typename DoFHandlerType,
+ template <int dim,
+ int spacedim,
typename SparsityPatternType,
typename number>
void
- make_flux_sparsity_pattern(const DoFHandlerType & dof,
- SparsityPatternType & sparsity,
- const AffineConstraints<number> &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<dim, spacedim> &dof,
+ SparsityPatternType & sparsity,
+ const AffineConstraints<number> &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<dim, spacedim>::active_cell_iterator &,
+ const unsigned int)> &face_has_flux_coupling)
{
- const FiniteElement<DoFHandlerType::dimension,
- DoFHandlerType::space_dimension> &fe = dof.get_fe();
+ const FiniteElement<dim, spacedim> &fe = dof.get_fe();
std::vector<types::global_dof_index> dofs_on_this_cell(
fe.dofs_per_cell);
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<DoFHandlerType::dimension>::faces_per_cell);
+ Table<2, bool> support_on_face(fe.dofs_per_cell,
+ GeometryInfo<dim>::faces_per_cell);
for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
- for (unsigned int f = 0;
- f < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
- ++f)
+ for (unsigned int f = 0; f < GeometryInfo<dim>::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
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<dim, spacedim>::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())) &&
bool_int_dof_mask);
// Loop over all interior neighbors
for (unsigned int face_n = 0;
- face_n <
- GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
+ face_n < GeometryInfo<dim>::faces_per_cell;
++face_n)
{
- const typename DoFHandlerType::face_iterator cell_face =
- cell->face(face_n);
+ const typename DoFHandler<dim, spacedim>::face_iterator
+ cell_face = cell->face(face_n);
const bool periodic_neighbor =
cell->has_periodic_neighbor(face_n);
}
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<dim, spacedim>::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
// 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);
sub_nr != cell_face->n_children();
++sub_nr)
{
- const typename DoFHandlerType::level_cell_iterator
- sub_neighbor =
+ const typename DoFHandler<dim, spacedim>::
+ level_cell_iterator sub_neighbor =
periodic_neighbor ?
cell->periodic_neighbor_child_on_subface(
face_n, sub_nr) :
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<dim, spacedim>::
+ 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
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
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<DoFHandlerType>);
}
+
+
template <typename DoFHandlerType,
typename SparsityPatternType,
typename number>
void
- make_flux_sparsity_pattern(const DoFHandlerType & dof,
- SparsityPatternType & sparsity,
- const AffineConstraints<number> &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<number> &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
"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
const bool,
const Table<2, Coupling> &,
const Table<2, Coupling> &,
- const types::subdomain_id);
+ const types::subdomain_id,
+ const std::function<bool(
+ const typename DoFHandler<deal_II_dimension>::active_cell_iterator &,
+ const unsigned int)> &);
template void
DoFTools::make_flux_sparsity_pattern<DoFHandler<deal_II_dimension>, SP>(
const bool,
const Table<2, Coupling> &,
const Table<2, Coupling> &,
- const types::subdomain_id);
+ const types::subdomain_id,
+ const std::function<bool(const typename hp::DoFHandler<
+ deal_II_dimension>::active_cell_iterator &,
+ const unsigned int)> &);
#if deal_II_dimension < 3
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/dof_handler.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#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 <int dim, class DoFHandlerType>
+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<dim> ¢er = cell->face(face_index)->center();
+ return std::abs(center[0]) < 1e-3;
+ return true;
+ };
+
+ AffineConstraints<double> 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 <int dim>
+void
+test_with_both_dof_handlers(const Triangulation<dim> &triangulation)
+{
+ const FE_Q<dim> element(1);
+
+ deallog << "dealii::DoFHandler" << std::endl;
+ {
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(element);
+ create_and_print_pattern<dim>(dof_handler);
+ }
+
+ deallog << std::endl;
+
+ deallog << "hp::DoFHandler" << std::endl;
+ {
+ hp::FECollection<dim> fe_collection(element);
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe_collection);
+ create_and_print_pattern<dim>(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<dim> triangulation;
+ create_3_cell_triangulation(triangulation);
+
+ test_with_both_dof_handlers(triangulation);
+ return 0;
+}
--- /dev/null
+
+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]
projects / dealii.git / commitdiff