From: Marco Feder <marco.feder@sissa.it>
Date: Wed, 3 May 2023 08:52:55 +0000 (+0200)
Subject: fix empty quad in MappingInfo::reinit_cells()
X-Git-Tag: v9.5.0-rc1~269^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15167%2Fhead;p=dealii.git

fix empty quad in MappingInfo::reinit_cells()
---

diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h
index 7a106c18ac..6312d22c24 100644
--- a/include/deal.II/non_matching/mapping_info.h
+++ b/include/deal.II/non_matching/mapping_info.h
@@ -402,7 +402,7 @@ namespace NonMatching
     std::vector<Quadrature<dim>> quadrature_vector(n_cells);
     for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
       quadrature_vector[cell_index] =
-        Quadrature<dim>(quadrature_vector[cell_index].get_points());
+        Quadrature<dim>(unit_points_vector[cell_index]);
 
     reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
   }
diff --git a/tests/non_matching/mapping_info_02.cc b/tests/non_matching/mapping_info_02.cc
new file mode 100644
index 0000000000..fb0ca11143
--- /dev/null
+++ b/tests/non_matching/mapping_info_02.cc
@@ -0,0 +1,120 @@
+// ---------------------------------------------------------------------
+//
+// 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 the NonMatching::MappingInfo class together with FEPointEvaluation when
+ * `reinit_cells` is called with a vector of vector of qpoints, rather than a
+ * vector of quadratures.
+ */
+
+#include <deal.II/base/function_signed_distance.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/non_matching/fe_values.h>
+#include <deal.II/non_matching/mapping_info.h>
+#include <deal.II/non_matching/mesh_classifier.h>
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+void
+test_vec_of_qpoints()
+{
+  constexpr unsigned int dim    = 2;
+  constexpr unsigned int degree = 1;
+
+  FE_Q<dim> fe_q(degree);
+
+  Triangulation<dim> tria;
+
+  MappingQ<dim> mapping(degree);
+
+  GridGenerator::subdivided_hyper_cube(tria, 4);
+
+  DoFHandler<dim> dof_handler(tria);
+
+  dof_handler.distribute_dofs(fe_q);
+
+  Functions::SignedDistance::Sphere<dim> level_set;
+
+  Vector<double> level_set_vec(dof_handler.n_dofs());
+
+  VectorTools::interpolate(dof_handler, level_set, level_set_vec);
+
+  NonMatching::MeshClassifier<dim> mesh_classifier(dof_handler, level_set_vec);
+  mesh_classifier.reclassify();
+
+  hp::QCollection<1> q_collection((QGauss<1>(degree)));
+
+  NonMatching::DiscreteQuadratureGenerator<dim> quadrature_generator(
+    q_collection, dof_handler, level_set_vec);
+
+  // FEPointEvaluation
+  NonMatching::MappingInfo<dim> mapping_info_cell(
+    mapping, update_values | update_gradients | update_JxW_values);
+
+  std::vector<std::vector<Point<dim>>> unit_points_vector;
+  for (const auto &cell : tria.active_cell_iterators())
+    {
+      quadrature_generator.generate(cell);
+      unit_points_vector.emplace_back(
+        quadrature_generator.get_inside_quadrature().get_points());
+    }
+
+  mapping_info_cell.reinit_cells(tria.active_cell_iterators(),
+                                 unit_points_vector);
+
+  Vector<double>      src(dof_handler.n_dofs());
+  std::vector<double> solution_values_in(fe_q.dofs_per_cell);
+
+  for (auto &v : src)
+    v = random_value<double>();
+  FEPointEvaluation<1, dim, dim, double> fe_point_cell(mapping_info_cell, fe_q);
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      cell->get_dof_values(src,
+                           solution_values_in.begin(),
+                           solution_values_in.end());
+
+      fe_point_cell.reinit(cell->active_cell_index());
+      fe_point_cell.evaluate(solution_values_in,
+                             EvaluationFlags::values |
+                               EvaluationFlags::gradients);
+      for (const auto q : fe_point_cell.quadrature_point_indices())
+        {
+          deallog << fe_point_cell.get_value(q) << std::endl;
+          deallog << fe_point_cell.get_gradient(q) << std::endl;
+        }
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  initlog();
+  test_vec_of_qpoints();
+}
diff --git a/tests/non_matching/mapping_info_02.output b/tests/non_matching/mapping_info_02.output
new file mode 100644
index 0000000000..28c8091576
--- /dev/null
+++ b/tests/non_matching/mapping_info_02.output
@@ -0,0 +1,33 @@
+
+DEAL::0.704027
+DEAL::-0.860928 0.693937
+DEAL::0.575505
+DEAL::-0.167248 -0.620078
+DEAL::0.553163
+DEAL::-0.0114928 -0.562178
+DEAL::0.488064
+DEAL::-0.543415 1.43809
+DEAL::0.671952
+DEAL::0.333629 -0.950543
+DEAL::0.497412
+DEAL::-1.72995 -0.00467201
+DEAL::0.460992
+DEAL::1.43859 -0.175191
+DEAL::0.678670
+DEAL::0.449138 -0.138598
+DEAL::0.664544
+DEAL::-0.258019 0.891278
+DEAL::0.586666
+DEAL::-0.365003 0.718707
+DEAL::0.436108
+DEAL::-0.302437 -0.0385718
+DEAL::0.471293
+DEAL::-1.75106 0.226205
+DEAL::0.451239
+DEAL::1.78818 -1.47462
+DEAL::0.473064
+DEAL::-0.171464 -2.58544
+DEAL::0.512506
+DEAL::0.0738393 -1.94578
+DEAL::0.515943
+DEAL::-1.71419 -0.781521