From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Mon, 14 Mar 2022 19:10:46 +0000 (+0100)
Subject: Add ScalarType as template parameter to MeshWorker::CopyData
X-Git-Tag: v9.4.0-rc1~360^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d4fc1d44af1397d08f7c4564d2b6d8ab32e48a3e;p=dealii.git

Add ScalarType as template parameter to MeshWorker::CopyData

The `matrices` and `vectors` can be made to store entries of type double, std::complex<double>, float, etc.
---

diff --git a/include/deal.II/meshworker/copy_data.h b/include/deal.II/meshworker/copy_data.h
index 1945437804..69f857b0bf 100644
--- a/include/deal.II/meshworker/copy_data.h
+++ b/include/deal.II/meshworker/copy_data.h
@@ -44,10 +44,13 @@ namespace MeshWorker
    * - @tparam n_matrices: Size of the array of matrices
    * - @tparam n_vectors: size of the array of vectors
    * - @tparam n_dof_indices: size of the array of local dof indices
+   * - @tparam ScalarType: The data type stored by the vectors and matrices.
+   *           This must be a real or complex floating point number.
    */
-  template <int n_matrices    = 1,
-            int n_vectors     = n_matrices,
-            int n_dof_indices = n_matrices>
+  template <int n_matrices      = 1,
+            int n_vectors       = n_matrices,
+            int n_dof_indices   = n_matrices,
+            typename ScalarType = double>
   struct CopyData
   {
     /**
@@ -74,8 +77,8 @@ namespace MeshWorker
     /**
      * Deep copy constructor.
      */
-    CopyData(const CopyData<n_matrices, n_vectors, n_dof_indices> &other) =
-      default;
+    CopyData(const CopyData<n_matrices, n_vectors, n_dof_indices, ScalarType>
+               &other) = default;
 
     /**
      * Reinitialize everything the same @p size. This is usually the number of
@@ -119,12 +122,12 @@ namespace MeshWorker
     /**
      * An array of local matrices.
      */
-    std::array<FullMatrix<double>, n_matrices> matrices;
+    std::array<FullMatrix<ScalarType>, n_matrices> matrices;
 
     /**
      * An array of local vectors.
      */
-    std::array<Vector<double>, n_vectors> vectors;
+    std::array<Vector<ScalarType>, n_vectors> vectors;
 
     /**
      * An array of local degrees of freedom indices.
@@ -138,8 +141,11 @@ namespace MeshWorker
   //
   // Template definitions
   //
-  template <int n_matrices, int n_vectors, int n_dof_indices>
-  CopyData<n_matrices, n_vectors, n_dof_indices>::CopyData(
+  template <int n_matrices,
+            int n_vectors,
+            int n_dof_indices,
+            typename ScalarType>
+  CopyData<n_matrices, n_vectors, n_dof_indices, ScalarType>::CopyData(
     const unsigned int size)
   {
     reinit(size);
@@ -147,8 +153,11 @@ namespace MeshWorker
 
 
 
-  template <int n_matrices, int n_vectors, int n_dof_indices>
-  CopyData<n_matrices, n_vectors, n_dof_indices>::CopyData(
+  template <int n_matrices,
+            int n_vectors,
+            int n_dof_indices,
+            typename ScalarType>
+  CopyData<n_matrices, n_vectors, n_dof_indices, ScalarType>::CopyData(
     const ndarray<unsigned int, n_matrices, 2> &   matrix_sizes,
     const std::array<unsigned int, n_vectors> &    vector_sizes,
     const std::array<unsigned int, n_dof_indices> &dof_indices_sizes)
@@ -165,9 +174,12 @@ namespace MeshWorker
 
 
 
-  template <int n_matrices, int n_vectors, int n_dof_indices>
+  template <int n_matrices,
+            int n_vectors,
+            int n_dof_indices,
+            typename ScalarType>
   void
-  CopyData<n_matrices, n_vectors, n_dof_indices>::reinit(
+  CopyData<n_matrices, n_vectors, n_dof_indices, ScalarType>::reinit(
     const unsigned int size)
   {
     for (auto &m : matrices)
@@ -180,9 +192,12 @@ namespace MeshWorker
 
 
 
-  template <int n_matrices, int n_vectors, int n_dof_indices>
+  template <int n_matrices,
+            int n_vectors,
+            int n_dof_indices,
+            typename ScalarType>
   void
-  CopyData<n_matrices, n_vectors, n_dof_indices>::reinit(
+  CopyData<n_matrices, n_vectors, n_dof_indices, ScalarType>::reinit(
     const unsigned int index,
     const unsigned int size)
   {
@@ -191,9 +206,12 @@ namespace MeshWorker
 
 
 
-  template <int n_matrices, int n_vectors, int n_dof_indices>
+  template <int n_matrices,
+            int n_vectors,
+            int n_dof_indices,
+            typename ScalarType>
   void
-  CopyData<n_matrices, n_vectors, n_dof_indices>::reinit(
+  CopyData<n_matrices, n_vectors, n_dof_indices, ScalarType>::reinit(
     const unsigned int index,
     const unsigned int size_rows,
     const unsigned int size_columns)
diff --git a/tests/meshworker/mesh_loop_09a.cc b/tests/meshworker/mesh_loop_09a.cc
new file mode 100644
index 0000000000..950b788a57
--- /dev/null
+++ b/tests/meshworker/mesh_loop_09a.cc
@@ -0,0 +1,143 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 - 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// test CopyData for complex values: vector
+// This test is based off meshworker/scratch_data_04.cc.
+
+#include <deal.II/base/iterator_range.h>
+
+#include <deal.II/fe/fe_q.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 <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <cmath>
+#include <complex>
+#include <fstream>
+#include <iostream>
+#include <unordered_map>
+
+#include "../tests.h"
+
+using namespace MeshWorker;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  Triangulation<dim, spacedim> tria;
+  FE_Q<dim, spacedim>          fe(1);
+  DoFHandler<dim, spacedim>    dh(tria);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(1);
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  dh.distribute_dofs(fe);
+
+  QGauss<dim>     quad(3);
+  QGauss<dim - 1> face_quad(3);
+
+  UpdateFlags cell_flags = update_JxW_values;
+  UpdateFlags face_flags = update_JxW_values;
+
+  using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
+  using CopyData    = MeshWorker::CopyData<0, 1, 0, std::complex<double>>;
+
+  ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags);
+  CopyData    copy(3); // store measures
+
+  using Iterator = decltype(dh.begin_active());
+
+  auto                       measures = std::make_tuple(0.0, 0.0, 0.0);
+  const std::complex<double> i{0, 1};
+
+
+  auto cell_worker = [&i](const Iterator &cell, ScratchData &s, CopyData &c) {
+    const auto &fev = s.reinit(cell);
+    const auto &JxW = s.get_JxW_values();
+    for (auto w : JxW)
+      c.vectors[0][0] += w * i;
+  };
+
+  auto boundary_worker = [&i](const Iterator &    cell,
+                              const unsigned int &f,
+                              ScratchData &       s,
+                              CopyData &          c) {
+    const auto &fev = s.reinit(cell, f);
+    const auto &JxW = s.get_JxW_values();
+    for (auto w : JxW)
+      c.vectors[0][1] += w * i;
+  };
+
+  auto face_worker = [&i](const Iterator &    cell,
+                          const unsigned int &f,
+                          const unsigned int &sf,
+                          const Iterator &    ncell,
+                          const unsigned int &nf,
+                          const unsigned int &nsf,
+                          ScratchData &       s,
+                          CopyData &          c) {
+    const auto &fev  = s.reinit(cell, f, sf);
+    const auto &nfev = s.reinit_neighbor(ncell, nf, nsf);
+
+    const auto &JxW  = s.get_JxW_values();
+    const auto &nJxW = s.get_neighbor_JxW_values();
+    for (auto w : JxW)
+      c.vectors[0][2] += w * i;
+  };
+
+  auto copier = [&measures](const CopyData &c) {
+    std::get<0>(measures) += std::imag(c.vectors[0][0]);
+    std::get<1>(measures) += std::imag(c.vectors[0][1]);
+    std::get<2>(measures) += std::imag(c.vectors[0][2]);
+  };
+
+  mesh_loop(dh.active_cell_iterators(),
+            cell_worker,
+            copier,
+            scratch,
+            copy,
+            assemble_own_cells | assemble_boundary_faces |
+              assemble_own_interior_faces_once,
+            boundary_worker,
+            face_worker);
+
+  deallog << "Testing <" << dim << ',' << spacedim << '>' << std::endl;
+
+  deallog << "Volume: " << std::get<0>(measures) << std::endl;
+
+  deallog << "Boundary surface: " << std::get<1>(measures) << std::endl;
+
+  deallog << "Interior surface: " << std::get<2>(measures) << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test<2, 2>();
+  test<2, 3>();
+  test<3, 3>();
+}
diff --git a/tests/meshworker/mesh_loop_09a.output b/tests/meshworker/mesh_loop_09a.output
new file mode 100644
index 0000000000..b9c1ecd3d3
--- /dev/null
+++ b/tests/meshworker/mesh_loop_09a.output
@@ -0,0 +1,13 @@
+
+DEAL::Testing <2,2>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 3.00000
+DEAL::Testing <2,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 3.00000
+DEAL::Testing <3,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 6.00000
+DEAL::Interior surface: 3.75000
diff --git a/tests/meshworker/mesh_loop_09b.cc b/tests/meshworker/mesh_loop_09b.cc
new file mode 100644
index 0000000000..8ab2d471a1
--- /dev/null
+++ b/tests/meshworker/mesh_loop_09b.cc
@@ -0,0 +1,143 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 - 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// test CopyData for complex values: matrix
+// This test is based off meshworker/scratch_data_04.cc.
+
+#include <deal.II/base/iterator_range.h>
+
+#include <deal.II/fe/fe_q.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 <deal.II/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <cmath>
+#include <complex>
+#include <fstream>
+#include <iostream>
+#include <unordered_map>
+
+#include "../tests.h"
+
+using namespace MeshWorker;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  Triangulation<dim, spacedim> tria;
+  FE_Q<dim, spacedim>          fe(1);
+  DoFHandler<dim, spacedim>    dh(tria);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(1);
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  dh.distribute_dofs(fe);
+
+  QGauss<dim>     quad(3);
+  QGauss<dim - 1> face_quad(3);
+
+  UpdateFlags cell_flags = update_JxW_values;
+  UpdateFlags face_flags = update_JxW_values;
+
+  using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
+  using CopyData    = MeshWorker::CopyData<1, 0, 0, std::complex<double>>;
+
+  ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags);
+  CopyData    copy(3); // store measures
+
+  using Iterator = decltype(dh.begin_active());
+
+  auto                       measures = std::make_tuple(0.0, 0.0, 0.0);
+  const std::complex<double> i{0, 1};
+
+
+  auto cell_worker = [&i](const Iterator &cell, ScratchData &s, CopyData &c) {
+    const auto &fev = s.reinit(cell);
+    const auto &JxW = s.get_JxW_values();
+    for (auto w : JxW)
+      c.matrices[0][0][0] += w * i;
+  };
+
+  auto boundary_worker = [&i](const Iterator &    cell,
+                              const unsigned int &f,
+                              ScratchData &       s,
+                              CopyData &          c) {
+    const auto &fev = s.reinit(cell, f);
+    const auto &JxW = s.get_JxW_values();
+    for (auto w : JxW)
+      c.matrices[0][0][1] += w * i;
+  };
+
+  auto face_worker = [&i](const Iterator &    cell,
+                          const unsigned int &f,
+                          const unsigned int &sf,
+                          const Iterator &    ncell,
+                          const unsigned int &nf,
+                          const unsigned int &nsf,
+                          ScratchData &       s,
+                          CopyData &          c) {
+    const auto &fev  = s.reinit(cell, f, sf);
+    const auto &nfev = s.reinit_neighbor(ncell, nf, nsf);
+
+    const auto &JxW  = s.get_JxW_values();
+    const auto &nJxW = s.get_neighbor_JxW_values();
+    for (auto w : JxW)
+      c.matrices[0][0][2] += w * i;
+  };
+
+  auto copier = [&measures](const CopyData &c) {
+    std::get<0>(measures) += std::imag(c.matrices[0][0][0]);
+    std::get<1>(measures) += std::imag(c.matrices[0][0][1]);
+    std::get<2>(measures) += std::imag(c.matrices[0][0][2]);
+  };
+
+  mesh_loop(dh.active_cell_iterators(),
+            cell_worker,
+            copier,
+            scratch,
+            copy,
+            assemble_own_cells | assemble_boundary_faces |
+              assemble_own_interior_faces_once,
+            boundary_worker,
+            face_worker);
+
+  deallog << "Testing <" << dim << ',' << spacedim << '>' << std::endl;
+
+  deallog << "Volume: " << std::get<0>(measures) << std::endl;
+
+  deallog << "Boundary surface: " << std::get<1>(measures) << std::endl;
+
+  deallog << "Interior surface: " << std::get<2>(measures) << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test<2, 2>();
+  test<2, 3>();
+  test<3, 3>();
+}
diff --git a/tests/meshworker/mesh_loop_09b.output b/tests/meshworker/mesh_loop_09b.output
new file mode 100644
index 0000000000..b9c1ecd3d3
--- /dev/null
+++ b/tests/meshworker/mesh_loop_09b.output
@@ -0,0 +1,13 @@
+
+DEAL::Testing <2,2>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 3.00000
+DEAL::Testing <2,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 4.00000
+DEAL::Interior surface: 3.00000
+DEAL::Testing <3,3>
+DEAL::Volume: 1.00000
+DEAL::Boundary surface: 6.00000
+DEAL::Interior surface: 3.75000