From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 25 Feb 2020 18:51:46 +0000 (+0100)
Subject: Add test case
X-Git-Tag: v9.2.0-rc1~432^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=20c95de33474c3d56eb31462ad98a0ae92d12c67;p=dealii.git

Add test case
---

diff --git a/tests/matrix_free/matrix_vector_faces_28.cc b/tests/matrix_free/matrix_vector_faces_28.cc
new file mode 100644
index 0000000000..8cfe4fb36d
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_faces_28.cc
@@ -0,0 +1,174 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that MatrixFree::update_mapping() works as expected. The test is
+// otherwise similar to matrix_vector_faces_03.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe_field.h>
+
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "create_mesh.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_faces_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  create_mesh(tria);
+
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  typename Triangulation<dim>::active_cell_iterator cell, endc;
+  cell = tria.begin_active();
+  endc = tria.end();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 0.5)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  // if (fe_degree == 1)
+  //  tria.refine_global(1);
+  cell = tria.begin_active();
+  for (unsigned int i = 0; i < 9 - 3 * dim; ++i)
+    {
+      cell                 = tria.begin_active();
+      endc                 = tria.end();
+      unsigned int counter = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (counter % (7 - i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  deallog << "Degree of element: " << fe_degree << std::endl;
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  FESystem<dim>   fe_grid(FE_Q<dim>(dof.get_fe().degree), dim);
+  DoFHandler<dim> dof_grid(tria);
+  dof_grid.distribute_dofs(fe_grid);
+  Vector<double> euler;
+  euler.reinit(dof_grid.n_dofs());
+  const ComponentMask mask(dim, true);
+  VectorTools::get_position_vector(dof_grid, euler, mask);
+
+  MappingFEField<dim> mapping(dof_grid, euler, mask);
+
+  Vector<double> in(dof.n_dofs()), out(dof.n_dofs());
+  Vector<double> out_dist(out);
+
+  // Set random seed for reproducibility
+  Testing::srand(42);
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      const double entry = Testing::rand() / (double)RAND_MAX;
+      in(i)              = entry;
+    }
+
+  auto create_reference_result = [&]() {
+    // assemble sparse matrix with MeshWorker
+    SparsityPattern      sparsity;
+    SparseMatrix<double> matrix;
+    {
+      DynamicSparsityPattern d_sparsity(dof.n_dofs());
+      DoFTools::make_flux_sparsity_pattern(dof, d_sparsity);
+      sparsity.copy_from(d_sparsity);
+    }
+    const unsigned int n_q_points_1d = dof.get_fe().degree + 1;
+    matrix.reinit(sparsity);
+    MeshWorker::IntegrationInfoBox<dim> info_box;
+    UpdateFlags                         update_flags =
+      update_values | update_gradients | update_jacobians;
+    info_box.add_update_flags_all(update_flags);
+    info_box.initialize_gauss_quadrature(n_q_points_1d,
+                                         n_q_points_1d,
+                                         n_q_points_1d);
+    info_box.initialize(dof.get_fe(), mapping);
+
+    MeshWorker::DoFInfo<dim> dof_info(dof);
+
+    MeshWorker::Assembler::MatrixSimple<SparseMatrix<double>> assembler;
+    assembler.initialize(matrix);
+
+    MatrixIntegrator<dim> integrator;
+    MeshWorker::integration_loop<dim, dim>(
+      dof.begin_active(), dof.end(), dof_info, info_box, integrator, assembler);
+
+    matrix.vmult(out, in);
+  };
+
+  MatrixFree<dim>                          mf_data;
+  const QGauss<1>                          quad(dof.get_fe().degree + 1);
+  typename MatrixFree<dim>::AdditionalData data;
+  data.tasks_parallel_scheme = MatrixFree<dim>::AdditionalData::none;
+  data.tasks_block_size      = 3;
+  data.mapping_update_flags_inner_faces =
+    (update_gradients | update_JxW_values);
+  data.mapping_update_flags_boundary_faces =
+    (update_gradients | update_JxW_values);
+
+  mf_data.reinit(mapping, dof, constraints, quad, data);
+
+  MatrixFreeTest<dim, fe_degree, fe_degree + 1, double, Vector<double>, 1> mf(
+    mf_data);
+  mf.vmult(out_dist, in);
+  create_reference_result();
+
+  out_dist -= out;
+  double diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+  deallog << "Norm of difference:           " << diff_norm << std::endl;
+
+  // deform the geometry slightly
+  euler(0) += 0.001;
+  euler(euler.size() - 1) += 0.001;
+
+  mf.vmult(out_dist, in);
+  create_reference_result();
+
+  out_dist -= out;
+  diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+  deallog << "Norm of difference no update: " << diff_norm << std::endl;
+
+  mf_data.update_mapping(mapping);
+  mf.vmult(out_dist, in);
+  out_dist -= out;
+  diff_norm = out_dist.linfty_norm() / out.linfty_norm();
+  deallog << "Norm of difference w update:  " << diff_norm << std::endl;
+}
diff --git a/tests/matrix_free/matrix_vector_faces_28.output b/tests/matrix_free/matrix_vector_faces_28.output
new file mode 100644
index 0000000000..11ef5b8048
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_faces_28.output
@@ -0,0 +1,17 @@
+
+DEAL:2d::Degree of element: 1
+DEAL:2d::Norm of difference:           4.87e-16
+DEAL:2d::Norm of difference no update: 0.0179
+DEAL:2d::Norm of difference w update:  4.87e-16
+DEAL:2d::Degree of element: 2
+DEAL:2d::Norm of difference:           1.92e-15
+DEAL:2d::Norm of difference no update: 0.00672
+DEAL:2d::Norm of difference w update:  2.50e-15
+DEAL:3d::Degree of element: 1
+DEAL:3d::Norm of difference:           1.88e-16
+DEAL:3d::Norm of difference no update: 5.11e-05
+DEAL:3d::Norm of difference w update:  1.88e-16
+DEAL:3d::Degree of element: 2
+DEAL:3d::Norm of difference:           1.06e-16
+DEAL:3d::Norm of difference no update: 5.82e-05
+DEAL:3d::Norm of difference w update:  1.06e-16