]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add MatrixFreeTools::categorize_accoring_boundary_ids_for_ecl 10871/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 1 Sep 2020 10:50:06 +0000 (12:50 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 4 Sep 2020 07:49:49 +0000 (09:49 +0200)
include/deal.II/matrix_free/tools.h
tests/matrix_free/categorize_by_boundary_ids_01.cc [new file with mode: 0644]
tests/matrix_free/categorize_by_boundary_ids_01.output [new file with mode: 0644]

index cfffa41214ede84890f75a564a1a9424677dccf5..ffcdd5b790d131edb20e861bd33521ea80e26c51 100644 (file)
@@ -18,6 +18,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/grid/tria.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -27,7 +28,83 @@ DEAL_II_NAMESPACE_OPEN
  * evaluation.
  */
 namespace MatrixFreeTools
-{} // namespace MatrixFreeTools
+{
+  /**
+   * Modify @p additional_data so that cells are categorized
+   * according to their boundary IDs, making face integrals in the case of
+   * cell-centric loop simpler.
+   */
+  template <int dim, typename AdditionalData>
+  void
+  categorize_by_boundary_ids(
+    const Triangulation<dim> &tria,
+    AdditionalData &          additional_data,
+    const unsigned int        level = numbers::invalid_unsigned_int)
+  {
+    bool is_mg = (level != numbers::invalid_unsigned_int);
+
+    // ... create list for the category of each cell
+    if (is_mg)
+      additional_data.cell_vectorization_category.resize(
+        std::distance(tria.begin(level), tria.end(level)));
+    else
+      additional_data.cell_vectorization_category.resize(tria.n_active_cells());
+
+    // ... set up scaling factor
+    std::vector<unsigned int> factors(dim * 2);
+
+    const auto bids = tria.get_boundary_ids();
+
+    std::map<unsigned int, unsigned int> bid_map;
+    for (const auto id : tria.get_boundary_ids())
+      bid_map[id] = bid_map.size() + 1;
+
+    {
+      unsigned int n_bids = bids.size() + 1;
+      int          offset = 1;
+      for (unsigned int i = 0; i < dim * 2; i++, offset = offset * n_bids)
+        factors[i] = offset;
+    }
+
+    auto to_category = [&](auto &cell) {
+      unsigned int c_num = 0;
+      for (unsigned int i = 0; i < dim * 2; i++)
+        {
+          auto &face = *cell->face(i);
+          if (face.at_boundary())
+            c_num += factors[i] * bid_map[face.boundary_id()];
+        }
+      return c_num;
+    };
+
+    if (!is_mg)
+      {
+        for (auto cell = tria.begin_active(); cell != tria.end(); ++cell)
+          {
+            if (cell->is_locally_owned())
+              additional_data
+                .cell_vectorization_category[cell->active_cell_index()] =
+                to_category(cell);
+          }
+      }
+    else
+      {
+        for (auto cell = tria.begin(level); cell != tria.end(level); ++cell)
+          {
+            if (cell->is_locally_owned_on_level())
+              additional_data.cell_vectorization_category[cell->index()] =
+                to_category(cell);
+          }
+      }
+
+    // ... finalize set up of matrix_free
+    additional_data.hold_all_faces_to_owned_cells        = true;
+    additional_data.cell_vectorization_categories_strict = true;
+    additional_data.mapping_update_flags_faces_by_cells =
+      additional_data.mapping_update_flags_inner_faces |
+      additional_data.mapping_update_flags_boundary_faces;
+  }
+} // namespace MatrixFreeTools
 
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/matrix_free/categorize_by_boundary_ids_01.cc b/tests/matrix_free/categorize_by_boundary_ids_01.cc
new file mode 100644 (file)
index 0000000..e0e582a
--- /dev/null
@@ -0,0 +1,128 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests MatrixFreeTools::categorize_by_boundary_ids()
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(const unsigned int n_refinements)
+{
+  const unsigned int fe_degree = 1;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+
+  tria.refine_global(n_refinements);
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  QGauss<1>                                        quadrature(fe_degree + 1);
+  typename MatrixFree<dim, double>::AdditionalData additional_data;
+  additional_data.mapping_update_flags                = update_values;
+  additional_data.mapping_update_flags_inner_faces    = update_values;
+  additional_data.mapping_update_flags_boundary_faces = update_values;
+  additional_data.tasks_parallel_scheme =
+    MatrixFree<dim, double>::AdditionalData::none;
+
+  AffineConstraints<double> dummy;
+  dummy.close();
+
+  const auto process = [&]() {
+    bool result = true;
+
+    MatrixFree<dim, double> data;
+    data.reinit(dof_handler, dummy, quadrature, additional_data);
+
+    using VectorType = Vector<double>;
+
+    VectorType vec;
+    data.initialize_dof_vector(vec);
+
+    data.template cell_loop<VectorType, VectorType>(
+      [&](const auto &, auto &, const auto &, const auto cell_range) {
+        for (unsigned int cell = cell_range.first; cell < cell_range.second;
+             ++cell)
+          {
+            for (unsigned int face = 0;
+                 face < GeometryInfo<dim>::faces_per_cell;
+                 ++face)
+              {
+                bool temp = true;
+
+                for (unsigned int v = 1; v < data.n_components_filled(cell);
+                     ++v)
+                  temp &= (data.get_faces_by_cells_boundary_id(cell, face)[0] ==
+                           data.get_faces_by_cells_boundary_id(cell, face)[v]);
+
+                result &= temp;
+              }
+          }
+      },
+      vec,
+      vec);
+
+    return result;
+  };
+
+  AssertDimension(process(), false); // not categorize
+
+  MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data);
+
+  AssertDimension(process(), true); // categorized
+
+  deallog << "OK!" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  {
+    deallog.push("2D");
+    test<2>(3);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("3D");
+    test<3>(3);
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/categorize_by_boundary_ids_01.output b/tests/matrix_free/categorize_by_boundary_ids_01.output
new file mode 100644 (file)
index 0000000..417aae8
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:2D::OK!
+DEAL:3D::OK!

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.