]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add new matrix-free test 18148/head
authorMartin Kronbichler <martin.kronbichler@rub.de>
Wed, 19 Feb 2025 07:31:30 +0000 (08:31 +0100)
committerMartin Kronbichler <martin.kronbichler@rub.de>
Tue, 25 Feb 2025 07:42:14 +0000 (08:42 +0100)
tests/matrix_free/integrate_functions_multife3.cc [new file with mode: 0644]
tests/matrix_free/integrate_functions_multife3.output [new file with mode: 0644]

diff --git a/tests/matrix_free/integrate_functions_multife3.cc b/tests/matrix_free/integrate_functions_multife3.cc
new file mode 100644 (file)
index 0000000..0034da2
--- /dev/null
@@ -0,0 +1,370 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2012 - 2025 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// similar to integrate_functions_multife, but using different quadrature
+// formula types QGauss and QGaussLobatto
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, int fe_degree, typename Number>
+class MatrixFreeTest
+{
+public:
+  using VectorType = std::vector<Vector<Number>>;
+
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in)
+    , fe_val0(data.get_dof_handler(0).get_fe(),
+              Quadrature<dim>(data.get_quadrature(0)),
+              update_values | update_gradients | update_JxW_values)
+    , fe_val01(data.get_dof_handler(0).get_fe(),
+               Quadrature<dim>(data.get_quadrature(1)),
+               update_values | update_gradients | update_JxW_values)
+    , fe_val1(data.get_dof_handler(1).get_fe(),
+              Quadrature<dim>(data.get_quadrature(1)),
+              update_values | update_gradients | update_JxW_values){};
+
+  void
+  operator()(const MatrixFree<dim, Number>               &data,
+             VectorType                                  &dst,
+             const VectorType                            &src,
+             const std::pair<unsigned int, unsigned int> &cell_range) const;
+
+  void
+  test_functions(VectorType &dst) const
+  {
+    for (unsigned int comp = 0; comp < dst.size(); ++comp)
+      dst[comp] = 0;
+    VectorType src_dummy;
+    data.cell_loop(&MatrixFreeTest<dim, fe_degree, Number>::operator(),
+                   this,
+                   dst,
+                   src_dummy);
+  };
+
+private:
+  const MatrixFree<dim, Number> &data;
+  mutable FEValues<dim>          fe_val0;
+  mutable FEValues<dim>          fe_val01;
+  mutable FEValues<dim>          fe_val1;
+};
+
+
+
+template <int dim, int fe_degree, typename Number>
+void
+MatrixFreeTest<dim, fe_degree, Number>::operator()(
+  const MatrixFree<dim, Number> &data,
+  std::vector<Vector<Number>>   &dst,
+  const std::vector<Vector<Number>> &,
+  const std::pair<unsigned int, unsigned int> &cell_range) const
+{
+  FEEvaluation<dim, fe_degree, fe_degree + 1, 1, Number> fe_eval0(data, 0, 0);
+  FEEvaluation<dim, fe_degree + 1, fe_degree + 2, 1, Number> fe_eval1(data,
+                                                                      1,
+                                                                      1);
+  FEEvaluation<dim, fe_degree, fe_degree + 2, 1, Number> fe_eval01(data, 0, 1);
+  const unsigned int n_q_points0    = fe_eval0.n_q_points;
+  const unsigned int n_q_points1    = fe_eval1.n_q_points;
+  const unsigned int dofs_per_cell0 = fe_eval0.dofs_per_cell;
+  const unsigned int dofs_per_cell1 = fe_eval1.dofs_per_cell;
+  AlignedVector<VectorizedArray<Number>> values0(n_q_points0);
+  AlignedVector<VectorizedArray<Number>> gradients0(dim * n_q_points0);
+  AlignedVector<VectorizedArray<Number>> values1(n_q_points1);
+  AlignedVector<VectorizedArray<Number>> gradients1(dim * n_q_points1);
+  std::vector<types::global_dof_index>   dof_indices0(dofs_per_cell0);
+  std::vector<types::global_dof_index>   dof_indices1(dofs_per_cell1);
+  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+    {
+      fe_eval0.reinit(cell);
+      fe_eval1.reinit(cell);
+      fe_eval01.reinit(cell);
+
+      // compare values with the ones the FEValues
+      // gives us. Those are seen as reference
+      for (unsigned int j = 0; j < data.n_active_entries_per_cell_batch(cell);
+           ++j)
+        {
+          // FE 0, Quad 0
+          // generate random numbers at quadrature
+          // points and test them with basis functions
+          // and their gradients
+          for (unsigned int q = 0; q < n_q_points0; ++q)
+            {
+              values0[q][j] = random_value<double>();
+              for (unsigned int d = 0; d < dim; ++d)
+                gradients0[q * dim + d][j] =
+                  -1. + 2. * (random_value<double>());
+            }
+          fe_val0.reinit(data.get_cell_iterator(cell, j, 0));
+          data.get_cell_iterator(cell, j, 0)->get_dof_indices(dof_indices0);
+
+          for (unsigned int i = 0; i < dofs_per_cell0; ++i)
+            {
+              double sum = 0.;
+              for (unsigned int q = 0; q < n_q_points0; ++q)
+                {
+                  sum +=
+                    values0[q][j] * fe_val0.shape_value(i, q) * fe_val0.JxW(q);
+                  for (unsigned int d = 0; d < dim; ++d)
+                    sum += (gradients0[q * dim + d][j] *
+                            fe_val0.shape_grad(i, q)[d] * fe_val0.JxW(q));
+                }
+              dst[0 + 1](dof_indices0[i]) += sum;
+            }
+
+          // FE 1, Quad 1
+          fe_val1.reinit(data.get_cell_iterator(cell, j, 1));
+          data.get_cell_iterator(cell, j, 1)->get_dof_indices(dof_indices1);
+
+          for (unsigned int q = 0; q < n_q_points1; ++q)
+            {
+              values1[q][j] = random_value<double>();
+              for (unsigned int d = 0; d < dim; ++d)
+                gradients1[q * dim + d][j] =
+                  -1. + 2. * (random_value<double>());
+            }
+          for (unsigned int i = 0; i < dofs_per_cell1; ++i)
+            {
+              double sum = 0.;
+              for (unsigned int q = 0; q < n_q_points1; ++q)
+                {
+                  sum +=
+                    values1[q][j] * fe_val1.shape_value(i, q) * fe_val1.JxW(q);
+                  for (unsigned int d = 0; d < dim; ++d)
+                    sum += (gradients1[q * dim + d][j] *
+                            fe_val1.shape_grad(i, q)[d] * fe_val1.JxW(q));
+                }
+              dst[2 + 1](dof_indices1[i]) += sum;
+            }
+
+          // FE 0, Quad 1
+          fe_val01.reinit(data.get_cell_iterator(cell, j, 0));
+          for (unsigned int i = 0; i < dofs_per_cell0; ++i)
+            {
+              double sum = 0.;
+              for (unsigned int q = 0; q < n_q_points1; ++q)
+                {
+                  sum += values1[q][j] * fe_val01.shape_value(i, q) *
+                         fe_val01.JxW(q);
+                  for (unsigned int d = 0; d < dim; ++d)
+                    sum += (gradients1[q * dim + d][j] *
+                            fe_val01.shape_grad(i, q)[d] * fe_val01.JxW(q));
+                }
+              dst[4 + 1](dof_indices0[i]) += sum;
+            }
+        }
+
+      // FE 0, Quad 0
+      for (unsigned int q = 0; q < n_q_points0; ++q)
+        {
+          fe_eval0.submit_value(values0[q], q);
+          Tensor<1, dim, VectorizedArray<Number>> submit;
+          for (unsigned int d = 0; d < dim; ++d)
+            submit[d] = gradients0[q * dim + d];
+          fe_eval0.submit_gradient(submit, q);
+        }
+      fe_eval0.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+      fe_eval0.distribute_local_to_global(dst[0]);
+
+      // FE 1, Quad 1
+      for (unsigned int q = 0; q < n_q_points1; ++q)
+        {
+          fe_eval1.submit_value(values1[q], q);
+          Tensor<1, dim, VectorizedArray<Number>> submit;
+          for (unsigned int d = 0; d < dim; ++d)
+            submit[d] = gradients1[q * dim + d];
+          fe_eval1.submit_gradient(submit, q);
+        }
+      fe_eval1.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+      fe_eval1.distribute_local_to_global(dst[2]);
+
+      // FE 0, Quad 1
+      for (unsigned int q = 0; q < n_q_points1; ++q)
+        {
+          fe_eval01.submit_value(values1[q], q);
+          Tensor<1, dim, VectorizedArray<Number>> submit;
+          for (unsigned int d = 0; d < dim; ++d)
+            submit[d] = gradients1[q * dim + d];
+          fe_eval01.submit_gradient(submit, q);
+        }
+      fe_eval01.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+      fe_eval01.distribute_local_to_global(dst[4]);
+    }
+}
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+  // create hyper ball geometry and refine some
+  // cells
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball(tria);
+
+  for (const auto &cell : tria.active_cell_iterators())
+    if (cell->center().norm() < 1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  for (const auto &cell : tria.active_cell_iterators())
+    if (cell->center().norm() < 0.2)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  if (dim < 3 || fe_degree < 2)
+    tria.refine_global(1);
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  for (unsigned int i = 0; i < 7 - 2 * dim; ++i)
+    {
+      unsigned int counter = 0;
+      for (const auto &cell : tria.active_cell_iterators())
+        {
+          if (counter % (7 - i) == 0)
+            cell->set_refine_flag();
+          ++counter;
+        }
+      tria.execute_coarsening_and_refinement();
+    }
+
+  FE_Q<dim>       fe0(fe_degree);
+  FE_Q<dim>       fe1(fe_degree + 1);
+  DoFHandler<dim> dof0(tria);
+  dof0.distribute_dofs(fe0);
+  DoFHandler<dim> dof1(tria);
+  dof1.distribute_dofs(fe1);
+
+  std::vector<const DoFHandler<dim> *> dof(2);
+  dof[0] = &dof0;
+  dof[1] = &dof1;
+
+  deallog << "Testing " << fe0.get_name() << " and " << fe1.get_name()
+          << std::endl;
+  // std::cout << "Number of cells: " << tria.n_active_cells() << std::endl;
+
+  std::vector<Vector<number>> dst(6);
+  dst[0].reinit(dof[0]->n_dofs());
+  dst[1].reinit(dst[0]);
+  dst[2].reinit(dof[1]->n_dofs());
+  dst[3].reinit(dst[2]);
+  dst[4].reinit(dst[0]);
+  dst[5].reinit(dst[0]);
+
+  std::vector<const AffineConstraints<double> *> constraints(2);
+  AffineConstraints<double>                      constraint0;
+  DoFTools::make_hanging_node_constraints(*dof[0], constraint0);
+  constraint0.close();
+  constraints[0] = &constraint0;
+  AffineConstraints<double> constraint1;
+  DoFTools::make_hanging_node_constraints(*dof[1], constraint1);
+  constraint1.close();
+  constraints[1] = &constraint1;
+
+  // std::cout << "Number of degrees of freedom FE 0: " << dof[0]->n_dofs() <<
+  // std::endl; std::cout << "Number of constraints FE 0: " <<
+  // constraints[0]->n_constraints() << std::endl; std::cout << "Number of
+  // degrees of freedom FE 1: " << dof[1]->n_dofs() << std::endl; std::cout <<
+  // "Number of constraints FE 1: " << constraints[1]->n_constraints() <<
+  // std::endl;
+
+  MatrixFree<dim, number> mf_data;
+  {
+    std::vector<Quadrature<1>> quad;
+    quad.push_back(QGauss<1>(fe_degree + 1));
+    quad.push_back(QGaussLobatto<1>(fe_degree + 2));
+    mf_data.reinit(MappingQ1<dim>{},
+                   dof,
+                   constraints,
+                   quad,
+                   typename MatrixFree<dim, number>::AdditionalData(
+                     MatrixFree<dim, number>::AdditionalData::none));
+  }
+
+  MatrixFreeTest<dim, fe_degree, number> mf(mf_data);
+  mf.test_functions(dst);
+
+  constraints[0]->condense(dst[1]);
+  constraints[1]->condense(dst[3]);
+  constraints[0]->condense(dst[5]);
+
+  dst[1] -= dst[0];
+  double diff_norm = dst[1].linfty_norm();
+  deallog << "FE 0, Quad 0; integration difference: " << diff_norm << std::endl;
+
+  dst[3] -= dst[2];
+  diff_norm = dst[3].linfty_norm();
+  deallog << "FE 1, Quad 1; integration difference: " << diff_norm << std::endl;
+
+  dst[5] -= dst[4];
+  diff_norm = dst[5].linfty_norm();
+  deallog << "FE 0, Quad 1; integration difference: " << diff_norm << std::endl
+          << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(3);
+
+  {
+    deallog.push("2d");
+    test<2, 1, double>();
+    test<2, 2, double>();
+    test<2, 3, double>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1, double>();
+    test<3, 2, double>();
+    deallog.pop();
+  }
+
+  {
+    deallog << std::endl << "Test with floats" << std::endl << std::endl;
+    deallog.push("2d");
+    test<2, 1, float>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1, float>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/integrate_functions_multife3.output b/tests/matrix_free/integrate_functions_multife3.output
new file mode 100644 (file)
index 0000000..3de4bbc
--- /dev/null
@@ -0,0 +1,39 @@
+
+DEAL:2d::Testing FE_Q<2>(1) and FE_Q<2>(2)
+DEAL:2d::FE 0, Quad 0; integration difference: 0
+DEAL:2d::FE 1, Quad 1; integration difference: 0
+DEAL:2d::FE 0, Quad 1; integration difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2) and FE_Q<2>(3)
+DEAL:2d::FE 0, Quad 0; integration difference: 0
+DEAL:2d::FE 1, Quad 1; integration difference: 0
+DEAL:2d::FE 0, Quad 1; integration difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3) and FE_Q<2>(4)
+DEAL:2d::FE 0, Quad 0; integration difference: 0
+DEAL:2d::FE 1, Quad 1; integration difference: 0
+DEAL:2d::FE 0, Quad 1; integration difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1) and FE_Q<3>(2)
+DEAL:3d::FE 0, Quad 0; integration difference: 0
+DEAL:3d::FE 1, Quad 1; integration difference: 0
+DEAL:3d::FE 0, Quad 1; integration difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2) and FE_Q<3>(3)
+DEAL:3d::FE 0, Quad 0; integration difference: 0
+DEAL:3d::FE 1, Quad 1; integration difference: 0
+DEAL:3d::FE 0, Quad 1; integration difference: 0
+DEAL:3d::
+DEAL::
+DEAL::Test with floats
+DEAL::
+DEAL:2d::Testing FE_Q<2>(1) and FE_Q<2>(2)
+DEAL:2d::FE 0, Quad 0; integration difference: 0
+DEAL:2d::FE 1, Quad 1; integration difference: 0
+DEAL:2d::FE 0, Quad 1; integration difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1) and FE_Q<3>(2)
+DEAL:3d::FE 0, Quad 0; integration difference: 0
+DEAL:3d::FE 1, Quad 1; integration difference: 0
+DEAL:3d::FE 0, Quad 1; integration difference: 0
+DEAL:3d::

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.