]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add NM::MappingInfo::AD to be able to reinit with global weights 15292/head
authorJohannes Heinz <johannes.heinz@tuwien.ac.at>
Thu, 1 Jun 2023 15:44:45 +0000 (17:44 +0200)
committerJohannes Heinz <johannes.heinz@tuwien.ac.at>
Wed, 7 Jun 2023 10:58:23 +0000 (12:58 +0200)
doc/news/changes/minor/20230601Heinz [new file with mode: 0644]
include/deal.II/non_matching/mapping_info.h
tests/non_matching/mapping_info_03.cc [new file with mode: 0644]
tests/non_matching/mapping_info_03.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20230601Heinz b/doc/news/changes/minor/20230601Heinz
new file mode 100644 (file)
index 0000000..af0099c
--- /dev/null
@@ -0,0 +1,3 @@
+Improved: Make it possible to fill JxW in NonMatching::MappingInfo from quadratures which contain JxW as weights.
+<br>
+(Johannes Heinz, 2023/06/01)
index 0f745240f39cb85082bba3f04e15ae5dee081454..85c7c512bdcc4322cd013367cec2b8ec86b0a420 100644 (file)
@@ -207,6 +207,27 @@ namespace NonMatching
     using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
       Number>::vectorized_value_type;
 
+    /**
+     * Collects the options which can be used to specify the behavior during
+     * reinitialization.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor which sets the default arguments.
+       */
+      AdditionalData(const bool use_global_weights = false)
+        : use_global_weights(use_global_weights)
+      {}
+
+      /**
+       * During initialization, assume that the Quadrature object contains
+       * global weights as, e.g., obtained by
+       * QSimplex::compute_affine_transformation().
+       */
+      bool use_global_weights;
+    };
+
     /**
      * Constructor.
      *
@@ -215,8 +236,13 @@ namespace NonMatching
      * @param update_flags Specify the quantities to be computed by the mapping
      * during the call of reinit(). These update flags are also handed to a
      * FEEvaluation object if you construct it with this MappingInfo object.
+     *
+     * @param additional_data Additional data for the class to specify the
+     * behavior during reinitialization.
      */
-    MappingInfo(const Mapping<dim> &mapping, const UpdateFlags update_flags);
+    MappingInfo(const Mapping<dim> & mapping,
+                const UpdateFlags    update_flags,
+                const AdditionalData additional_data = AdditionalData());
 
     /**
      * Compute the mapping information for the incoming cell and unit
@@ -448,10 +474,11 @@ namespace NonMatching
      * Store the requested mapping data.
      */
     void
-    store_mapping_data(const unsigned int unit_points_index_offset,
-                       const unsigned int n_q_points,
-                       const unsigned int n_q_points_unvectorized,
-                       const MappingData &mapping_data);
+    store_mapping_data(const unsigned int         unit_points_index_offset,
+                       const unsigned int         n_q_points,
+                       const unsigned int         n_q_points_unvectorized,
+                       const MappingData &        mapping_data,
+                       const std::vector<double> &weights);
 
     /**
      * Compute the compressed cell index.
@@ -523,6 +550,11 @@ namespace NonMatching
      */
     UpdateFlags update_flags_mapping;
 
+    /**
+     * AdditionalData for this object.
+     */
+    const AdditionalData additional_data;
+
     /**
      * Stores the index offset into the arrays @p JxW_values, @p jacobians,
      * @p inverse_jacobians and @p normal_vectors.
@@ -602,11 +634,13 @@ namespace NonMatching
 
   template <int dim, int spacedim, typename Number>
   MappingInfo<dim, spacedim, Number>::MappingInfo(
-    const Mapping<dim> &mapping,
-    const UpdateFlags   update_flags)
+    const Mapping<dim> & mapping,
+    const UpdateFlags    update_flags,
+    const AdditionalData additional_data)
     : mapping(&mapping)
     , update_flags(update_flags)
     , update_flags_mapping(update_default)
+    , additional_data(additional_data)
   {
     // translate update flags
     if (update_flags & update_jacobians)
@@ -704,7 +738,8 @@ namespace NonMatching
     store_mapping_data(0,
                        n_q_points_data,
                        n_q_points_unvectorized[0],
-                       mapping_data);
+                       mapping_data,
+                       quadrature.get_weights());
 
     state = State::single_cell;
     is_reinitialized();
@@ -814,7 +849,8 @@ namespace NonMatching
         store_mapping_data(data_index_offsets[cell_index],
                            n_q_points_data,
                            n_q_points_unvectorized[cell_index],
-                           mapping_data);
+                           mapping_data,
+                           quadrature_vector[cell_index].get_weights());
 
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -837,6 +873,11 @@ namespace NonMatching
     const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
     const unsigned int                                 n_unfiltered_cells)
   {
+    Assert(
+      additional_data.use_global_weights == false,
+      ExcMessage(
+        "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
+
     do_cell_index_compression =
       n_unfiltered_cells != numbers::invalid_unsigned_int;
 
@@ -911,7 +952,8 @@ namespace NonMatching
         store_mapping_data(data_index_offsets[cell_index],
                            n_q_points_data,
                            n_q_points_unvectorized[cell_index],
-                           mapping_data);
+                           mapping_data,
+                           quadrature_vector[cell_index].get_weights());
 
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -1051,7 +1093,8 @@ namespace NonMatching
             store_mapping_data(data_index_offsets[current_face_index],
                                n_q_points_data,
                                n_q_points_unvectorized[current_face_index],
-                               mapping_data);
+                               mapping_data,
+                               quadrature_on_face.get_weights());
           }
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -1219,7 +1262,8 @@ namespace NonMatching
     const unsigned int              unit_points_index_offset,
     const unsigned int              n_q_points,
     const unsigned int              n_q_points_unvectorized,
-    const MappingInfo::MappingData &mapping_data)
+    const MappingInfo::MappingData &mapping_data,
+    const std::vector<double> &     weights)
   {
     const unsigned int n_lanes =
       dealii::internal::VectorizedArrayTrait<Number>::width();
@@ -1244,9 +1288,19 @@ namespace NonMatching
                     inverse_jacobians[offset][s][d], v) =
                     mapping_data.inverse_jacobians[q * n_lanes + v][s][d];
             if (update_flags_mapping & UpdateFlags::update_JxW_values)
-              dealii::internal::VectorizedArrayTrait<Number>::get(
-                JxW_values[offset], v) =
-                mapping_data.JxW_values[q * n_lanes + v];
+              {
+                if (additional_data.use_global_weights)
+                  {
+                    dealii::internal::VectorizedArrayTrait<Number>::get(
+                      JxW_values[offset], v) = weights[q * n_lanes + v];
+                  }
+                else
+                  {
+                    dealii::internal::VectorizedArrayTrait<Number>::get(
+                      JxW_values[offset], v) =
+                      mapping_data.JxW_values[q * n_lanes + v];
+                  }
+              }
             if (update_flags_mapping & UpdateFlags::update_normal_vectors)
               for (unsigned int s = 0; s < spacedim; ++s)
                 dealii::internal::VectorizedArrayTrait<Number>::get(
diff --git a/tests/non_matching/mapping_info_03.cc b/tests/non_matching/mapping_info_03.cc
new file mode 100644 (file)
index 0000000..cda1208
--- /dev/null
@@ -0,0 +1,150 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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 NonMatching::MappingInfo if reinit with quadratures which contain JxW as
+ * weights.
+ */
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/non_matching/mapping_info.h>
+
+#include "../tests.h"
+
+void
+test()
+{
+  constexpr unsigned int dim          = 2;
+  constexpr unsigned int spacedim     = 2;
+  constexpr unsigned int degree       = 1;
+  constexpr unsigned int n_components = 1;
+  using Number                        = double;
+
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::subdivided_hyper_cube(tria, 4);
+
+  FE_Q<dim>     fe(degree);
+  MappingQ<dim> mapping(degree);
+
+
+  deallog << "Check JxW faces..." << std::endl;
+  {
+    NonMatching::MappingInfo<dim>::AdditionalData additional_data;
+    additional_data.use_global_weights = true;
+    NonMatching::MappingInfo<dim> mapping_info(mapping,
+                                               update_JxW_values,
+                                               additional_data);
+
+    // 1) build vector of quadratures
+    std::vector<std::vector<Quadrature<dim - 1>>> quad_vec;
+    // prescribe JxW (there is no meaning in the actual values, they just have
+    // to stay the same when fetched with FEPointEvaluation)
+    double JxW = 1.0;
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        std::vector<Quadrature<dim - 1>> quad;
+        for (auto f : cell->face_indices())
+          {
+            dealii::QGauss<dim - 1> face_quadrature(degree + 1);
+            std::vector<double> weights(face_quadrature.get_weights().size());
+            for (auto &w : weights)
+              {
+                w = JxW;
+                JxW += 1.0;
+              }
+
+            quad.emplace_back(
+              Quadrature<dim - 1>(face_quadrature.get_points(), weights));
+          }
+        quad_vec.push_back(quad);
+      }
+
+    // 2) reinit mapping info
+    mapping_info.reinit_faces(tria.active_cell_iterators(), quad_vec);
+
+    FEPointEvaluation<n_components, dim, spacedim, Number> fe_point_eval(
+      mapping_info, fe);
+
+    // 3) print JxW
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        for (auto f : cell->face_indices())
+          {
+            fe_point_eval.reinit(cell->active_cell_index(), f);
+            for (unsigned int q : fe_point_eval.quadrature_point_indices())
+              deallog << fe_point_eval.JxW(q) << std::endl;
+          }
+      }
+  }
+
+  deallog << "\n\nCheck JxW cells..." << std::endl;
+  {
+    NonMatching::MappingInfo<dim>::AdditionalData additional_data;
+    additional_data.use_global_weights = true;
+    NonMatching::MappingInfo<dim> mapping_info(mapping,
+                                               update_JxW_values,
+                                               additional_data);
+
+
+    // 1) build vector of quadratures
+    std::vector<Quadrature<dim>> quad_vec;
+    // prescribe JxW (there is no meaning in the actual values, they just have
+    // to stay the same when fetched with FEPointEvaluation)
+    double JxW = 1.0;
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        dealii::QGauss<dim> cell_quadrature(degree + 1);
+        std::vector<double> weights(cell_quadrature.get_weights().size());
+        for (auto &w : weights)
+          {
+            w = JxW;
+            JxW += 1.0;
+          }
+
+        quad_vec.emplace_back(
+          Quadrature<dim>(cell_quadrature.get_points(), weights));
+      }
+
+    // 2) reinit mapping info
+    mapping_info.reinit_cells(tria.active_cell_iterators(), quad_vec);
+    FEPointEvaluation<n_components, dim, spacedim, Number> fe_point_eval(
+      mapping_info, fe);
+
+    // 3) print JxW
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        fe_point_eval.reinit(cell->active_cell_index());
+        for (unsigned int q : fe_point_eval.quadrature_point_indices())
+          deallog << fe_point_eval.JxW(q) << std::endl;
+      }
+  }
+
+
+  deallog << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  initlog();
+  test();
+}
diff --git a/tests/non_matching/mapping_info_03.output b/tests/non_matching/mapping_info_03.output
new file mode 100644 (file)
index 0000000..7c1ee20
--- /dev/null
@@ -0,0 +1,198 @@
+
+DEAL::Check JxW faces...
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL::6.00000
+DEAL::7.00000
+DEAL::8.00000
+DEAL::9.00000
+DEAL::10.0000
+DEAL::11.0000
+DEAL::12.0000
+DEAL::13.0000
+DEAL::14.0000
+DEAL::15.0000
+DEAL::16.0000
+DEAL::17.0000
+DEAL::18.0000
+DEAL::19.0000
+DEAL::20.0000
+DEAL::21.0000
+DEAL::22.0000
+DEAL::23.0000
+DEAL::24.0000
+DEAL::25.0000
+DEAL::26.0000
+DEAL::27.0000
+DEAL::28.0000
+DEAL::29.0000
+DEAL::30.0000
+DEAL::31.0000
+DEAL::32.0000
+DEAL::33.0000
+DEAL::34.0000
+DEAL::35.0000
+DEAL::36.0000
+DEAL::37.0000
+DEAL::38.0000
+DEAL::39.0000
+DEAL::40.0000
+DEAL::41.0000
+DEAL::42.0000
+DEAL::43.0000
+DEAL::44.0000
+DEAL::45.0000
+DEAL::46.0000
+DEAL::47.0000
+DEAL::48.0000
+DEAL::49.0000
+DEAL::50.0000
+DEAL::51.0000
+DEAL::52.0000
+DEAL::53.0000
+DEAL::54.0000
+DEAL::55.0000
+DEAL::56.0000
+DEAL::57.0000
+DEAL::58.0000
+DEAL::59.0000
+DEAL::60.0000
+DEAL::61.0000
+DEAL::62.0000
+DEAL::63.0000
+DEAL::64.0000
+DEAL::65.0000
+DEAL::66.0000
+DEAL::67.0000
+DEAL::68.0000
+DEAL::69.0000
+DEAL::70.0000
+DEAL::71.0000
+DEAL::72.0000
+DEAL::73.0000
+DEAL::74.0000
+DEAL::75.0000
+DEAL::76.0000
+DEAL::77.0000
+DEAL::78.0000
+DEAL::79.0000
+DEAL::80.0000
+DEAL::81.0000
+DEAL::82.0000
+DEAL::83.0000
+DEAL::84.0000
+DEAL::85.0000
+DEAL::86.0000
+DEAL::87.0000
+DEAL::88.0000
+DEAL::89.0000
+DEAL::90.0000
+DEAL::91.0000
+DEAL::92.0000
+DEAL::93.0000
+DEAL::94.0000
+DEAL::95.0000
+DEAL::96.0000
+DEAL::97.0000
+DEAL::98.0000
+DEAL::99.0000
+DEAL::100.000
+DEAL::101.000
+DEAL::102.000
+DEAL::103.000
+DEAL::104.000
+DEAL::105.000
+DEAL::106.000
+DEAL::107.000
+DEAL::108.000
+DEAL::109.000
+DEAL::110.000
+DEAL::111.000
+DEAL::112.000
+DEAL::113.000
+DEAL::114.000
+DEAL::115.000
+DEAL::116.000
+DEAL::117.000
+DEAL::118.000
+DEAL::119.000
+DEAL::120.000
+DEAL::121.000
+DEAL::122.000
+DEAL::123.000
+DEAL::124.000
+DEAL::125.000
+DEAL::126.000
+DEAL::127.000
+DEAL::128.000
+DEAL::
+
+Check JxW cells...
+DEAL::1.00000
+DEAL::2.00000
+DEAL::3.00000
+DEAL::4.00000
+DEAL::5.00000
+DEAL::6.00000
+DEAL::7.00000
+DEAL::8.00000
+DEAL::9.00000
+DEAL::10.0000
+DEAL::11.0000
+DEAL::12.0000
+DEAL::13.0000
+DEAL::14.0000
+DEAL::15.0000
+DEAL::16.0000
+DEAL::17.0000
+DEAL::18.0000
+DEAL::19.0000
+DEAL::20.0000
+DEAL::21.0000
+DEAL::22.0000
+DEAL::23.0000
+DEAL::24.0000
+DEAL::25.0000
+DEAL::26.0000
+DEAL::27.0000
+DEAL::28.0000
+DEAL::29.0000
+DEAL::30.0000
+DEAL::31.0000
+DEAL::32.0000
+DEAL::33.0000
+DEAL::34.0000
+DEAL::35.0000
+DEAL::36.0000
+DEAL::37.0000
+DEAL::38.0000
+DEAL::39.0000
+DEAL::40.0000
+DEAL::41.0000
+DEAL::42.0000
+DEAL::43.0000
+DEAL::44.0000
+DEAL::45.0000
+DEAL::46.0000
+DEAL::47.0000
+DEAL::48.0000
+DEAL::49.0000
+DEAL::50.0000
+DEAL::51.0000
+DEAL::52.0000
+DEAL::53.0000
+DEAL::54.0000
+DEAL::55.0000
+DEAL::56.0000
+DEAL::57.0000
+DEAL::58.0000
+DEAL::59.0000
+DEAL::60.0000
+DEAL::61.0000
+DEAL::62.0000
+DEAL::63.0000
+DEAL::64.0000
+DEAL::

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.