Merge pull request #12952 from bangerth/filters

Provide C++20-style operator| to filter iterator ranges.

doc/news/changes/minor/20211115Bangerth-c

@@ -0,0 +1,42 @@
+New: In the same spirit as provided by the C++20 [range
+adaptors](https://en.cppreference.com/w/cpp/ranges) feature, it is now
+possible to "filter" ranges over the active cells of Triangulation
+and DoFHandler objects.
+<br>
+This allows to replace:
+@code
+  DoFHandler<dim> dof_handler;
+  ...
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      if (cell->is_locally_owned())
+        {
+          fe_values.reinit (cell);
+          ...do the local integration on 'cell'...;
+        }
+    }
+@endcode
+by:
+@code
+  DoFHandler<dim> dof_handler;
+  ...
+  const auto filtered_iterators_range =
+    filter_iterators();
+  for (const auto &cell :
+          dof_handler.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+    {
+      fe_values.reinit (cell);
+      ...do the local integration on 'cell'...;
+    }
+@endcode
+Here, the `operator|` is to be interpreted in the same way as is done in
+the [range adaptors](https://en.cppreference.com/w/cpp/ranges) feature
+that is part of [C++20](https://en.wikipedia.org/wiki/C%2B%2B20). It has
+the same meaning as the `|` symbol on the command line: It takes what is
+on its left as its inputs, and filters and transforms to produce some
+output. In the example above, it "filters" all of the active cell iterators
+and removes those that do not satisfy the predicate -- that is, it produces
+a range of iterators that only contains those cells that are both active
+and locally owned.
+<br>
+(Wolfgang Bangerth, 2021/11/03)

include/deal.II/grid/filtered_iterator.h

@@ -856,7 +856,7 @@ namespace internal
  * @ingroup CPP11
  */
 template <typename BaseIterator, typename Predicate>
-IteratorRange<FilteredIterator<BaseIterator>>
+inline IteratorRange<FilteredIterator<BaseIterator>>
 filter_iterators(IteratorRange<BaseIterator> i, const Predicate &p)
 {
   FilteredIterator<BaseIterator> fi(p, *(i.begin()));
@@ -917,6 +917,57 @@ filter_iterators(IteratorRange<BaseIterator> i,
 }
 
 
+
+/**
+ * Filter the  given range of iterators using a predicate. This allows to
+ * replace:
+ * @code
+ *   DoFHandler<dim> dof_handler;
+ *   ...
+ *   for (const auto &cell : dof_handler.active_cell_iterators())
+ *     {
+ *       if (cell->is_locally_owned())
+ *         {
+ *           fe_values.reinit (cell);
+ *           ...do the local integration on 'cell'...;
+ *         }
+ *     }
+ * @endcode
+ * by:
+ * @code
+ *   DoFHandler<dim> dof_handler;
+ *   ...
+ *   const auto filtered_iterators_range =
+ *     filter_iterators();
+ *   for (const auto &cell :
+ *           dof_handler.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+ *     {
+ *       fe_values.reinit (cell);
+ *       ...do the local integration on 'cell'...;
+ *     }
+ * @endcode
+ * Here, the `operator|` is to be interpreted in the same way as is done in
+ * the [range adaptors](https://en.cppreference.com/w/cpp/ranges) feature
+ * that is part of [C++20](https://en.wikipedia.org/wiki/C%2B%2B20). It has
+ * the same meaning as the `|` symbol on the command line: It takes what is
+ * on its left as its inputs, and filters and transforms to produce some
+ * output. In the example above, it "filters" all of the active cell iterators
+ * and removes those that do not satisfy the predicate -- that is, it produces
+ * a range of iterators that only contains those cells that are both active
+ * and locally owned.
+ *
+ * @relatesalso FilteredIterator
+ * @ingroup CPP11
+ */
+template <typename BaseIterator, typename Predicate>
+inline IteratorRange<FilteredIterator<BaseIterator>>
+operator|(IteratorRange<BaseIterator> i, const Predicate &p)
+{
+  return filter_iterators(i, p);
+}
+
+
+
 /* ------------------ Inline functions and templates ------------ */
 
 

source/particles/generators.cc

@@ -205,23 +205,19 @@ namespace Particles
       particle_handler.reserve(particle_handler.n_locally_owned_particles() +
                                n_particles_to_generate);
 
-      for (const auto &cell : triangulation.active_cell_iterators())
+      for (const auto &cell : triangulation.active_cell_iterators() |
+                                IteratorFilters::LocallyOwnedCell())
         {
-          if (cell->is_locally_owned())
+          for (const auto &reference_location : particle_reference_locations)
             {
-              for (const auto &reference_location :
-                   particle_reference_locations)
-                {
-                  const Point<spacedim> position_real =
-                    mapping.transform_unit_to_real_cell(cell,
-                                                        reference_location);
-
-                  particle_handler.insert_particle(position_real,
-                                                   reference_location,
-                                                   particle_index,
-                                                   cell);
-                  ++particle_index;
-                }
+              const Point<spacedim> position_real =
+                mapping.transform_unit_to_real_cell(cell, reference_location);
+
+              particle_handler.insert_particle(position_real,
+                                               reference_location,
+                                               particle_index,
+                                               cell);
+              ++particle_index;
             }
         }
 
@@ -398,23 +394,23 @@ namespace Particles
             // between their weight and the local weight integral
             types::particle_index particles_created = 0;
 
-            for (const auto &cell : triangulation.active_cell_iterators())
-              if (cell->is_locally_owned())
-                {
-                  const types::particle_index cumulative_particles_to_create =
-                    std::llround(
-                      static_cast<double>(n_local_particles) *
-                      cumulative_cell_weights[cell->active_cell_index()] /
-                      local_weight_integral);
-
-                  // Compute particles for this cell as difference between
-                  // number of particles that should be created including this
-                  // cell minus the number of particles already created.
-                  particles_per_cell[cell->active_cell_index()] =
-                    cumulative_particles_to_create - particles_created;
-                  particles_created +=
-                    particles_per_cell[cell->active_cell_index()];
-                }
+            for (const auto &cell : triangulation.active_cell_iterators() |
+                                      IteratorFilters::LocallyOwnedCell())
+              {
+                const types::particle_index cumulative_particles_to_create =
+                  std::llround(
+                    static_cast<double>(n_local_particles) *
+                    cumulative_cell_weights[cell->active_cell_index()] /
+                    local_weight_integral);
+
+                // Compute particles for this cell as difference between
+                // number of particles that should be created including this
+                // cell minus the number of particles already created.
+                particles_per_cell[cell->active_cell_index()] =
+                  cumulative_particles_to_create - particles_created;
+                particles_created +=
+                  particles_per_cell[cell->active_cell_index()];
+              }
           }
       }
 
@@ -424,22 +420,22 @@ namespace Particles
                                  n_local_particles);
         unsigned int current_particle_index = start_particle_id;
 
-        for (const auto &cell : triangulation.active_cell_iterators())
-          if (cell->is_locally_owned())
-            {
-              for (unsigned int i = 0;
-                   i < particles_per_cell[cell->active_cell_index()];
-                   ++i)
-                {
-                  random_particle_in_cell_insert(cell,
-                                                 current_particle_index,
-                                                 random_number_generator,
-                                                 particle_handler,
-                                                 mapping);
-
-                  ++current_particle_index;
-                }
-            }
+        for (const auto &cell : triangulation.active_cell_iterators() |
+                                  IteratorFilters::LocallyOwnedCell())
+          {
+            for (unsigned int i = 0;
+                 i < particles_per_cell[cell->active_cell_index()];
+                 ++i)
+              {
+                random_particle_in_cell_insert(cell,
+                                               current_particle_index,
+                                               random_number_generator,
+                                               particle_handler,
+                                               mapping);
+
+                ++current_particle_index;
+              }
+          }
 
         particle_handler.update_cached_numbers();
       }
@@ -503,18 +499,14 @@ namespace Particles
                                  triangulation.n_active_cells());
 
       //       Loop through cells and gather gauss points
-      for (const auto &cell : triangulation.active_cell_iterators())
+      for (const auto &cell : triangulation.active_cell_iterators() |
+                                IteratorFilters::LocallyOwnedCell())
         {
-          if (cell->is_locally_owned())
+          for (const auto &reference_location : particle_reference_locations)
             {
-              for (const auto &reference_location :
-                   particle_reference_locations)
-                {
-                  const Point<spacedim> position_real =
-                    mapping.transform_unit_to_real_cell(cell,
-                                                        reference_location);
-                  points_to_generate.push_back(position_real);
-                }
+              const Point<spacedim> position_real =
+                mapping.transform_unit_to_real_cell(cell, reference_location);
+              points_to_generate.push_back(position_real);
             }
         }
       particle_handler.insert_global_particles(points_to_generate,

tests/grid/filtered_iterator_04.cc

@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 
-// check filtered iterators using multiple predicate
+// check filtered iterators using multiple predicates
 
 #include <deal.II/grid/filtered_iterator.h>
 #include <deal.II/grid/grid_generator.h>

tests/grid/filtered_iterator_04_operator.cc

@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check filtered iterators using multiple predicates
+//
+// Compared to filtered_iterator_04, this test checks the availability
+// of operator| to create the filter.
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <set>
+
+#include "../tests.h"
+
+using active_cell_iterator = Triangulation<2>::active_cell_iterator;
+
+void
+test()
+{
+  Triangulation<2> tria;
+  GridGenerator::hyper_cube(tria, -1, 1);
+  tria.refine_global(1);
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global(2);
+
+  // we now have a number of cells,
+  // flag them with some subdomain
+  // ids based on their position, in
+  // particular we take the quadrant
+  // (octant)
+  active_cell_iterator cell = tria.begin_active(), endc = tria.end();
+  for (unsigned int i = 0; cell != endc; ++cell)
+    {
+      unsigned int subdomain = i % 3;
+
+      cell->set_subdomain_id(subdomain);
+      ++i;
+    };
+
+  // Count the cells that are on the boundary and have a subdomain_id of 0
+  std::set<active_cell_iterator> cell_set;
+  for (cell = tria.begin_active(); cell != endc; ++cell)
+    if ((cell->subdomain_id() == 0) && (cell->at_boundary()))
+      cell_set.insert(cell);
+
+
+  unsigned int n_filtered_cells = 0;
+  for (auto filtered_cell : tria.active_cell_iterators() |
+                              IteratorFilters::AtBoundary() |
+                              IteratorFilters::SubdomainEqualTo(0))
+    {
+      AssertThrow(cell_set.count(filtered_cell) == 1,
+                  ExcMessage("Wrong cell filtered."));
+      ++n_filtered_cells;
+    }
+  AssertThrow(n_filtered_cells == cell_set.size(),
+              ExcMessage("Filtered cells missing."));
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(4);
+
+  test();
+
+  deallog << "OK" << std::endl;
+  ;
+
+  return 0;
+}

tests/grid/filtered_iterator_04_operator.output

@@ -0,0 +1,2 @@
+
+DEAL::OK