]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix HangingNodes::setup_constraints() for multiple components 12702/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 24 Aug 2021 08:34:09 +0000 (10:34 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 25 Aug 2021 17:11:18 +0000 (19:11 +0200)
include/deal.II/matrix_free/hanging_nodes_internal.h
tests/matrix_free/matrix_vector_stokes_noflux2.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_stokes_noflux2.output [new file with mode: 0644]

index 88a1c7c1e92efe006fa04c6fb5cf551cdf22ffcb..225b654d047403eeb3124f2ccdd24f987412bf5b 100644 (file)
@@ -346,18 +346,22 @@ namespace internal
             auto &mask = masks[comp];
             mask       = ConstraintKinds::unconstrained;
 
-            const auto &fe = cell->get_fe().base_element(base_element_index);
+            const auto &fe_base =
+              cell->get_fe().base_element(base_element_index);
 
-            if (dim == 1 || dynamic_cast<const FE_Q<dim> *>(&fe) == nullptr)
+            if (dim == 1 ||
+                dynamic_cast<const FE_Q<dim> *>(&fe_base) == nullptr)
               continue;
 
-            const unsigned int fe_degree = fe.tensor_degree();
+            const unsigned int fe_degree = fe_base.tensor_degree();
             const unsigned int n_dofs_1d = fe_degree + 1;
             const unsigned int dofs_per_face =
               Utilities::fixed_power<dim - 1>(n_dofs_1d);
 
             std::vector<types::global_dof_index> neighbor_dofs_all(
               idx_offset.back());
+            std::vector<types::global_dof_index> neighbor_dofs_all_temp(
+              idx_offset.back());
 
             std::vector<types::global_dof_index> neighbor_dofs(dofs_per_face);
 
@@ -659,20 +663,24 @@ namespace internal
                                   ExcInternalError();
 
                                 // Copy the unconstrained values
-                                neighbor_dofs.resize(n_dofs_1d * n_dofs_1d *
-                                                     n_dofs_1d);
                                 DoFCellAccessor<dim, dim, false>(
                                   &neighbor_cell->get_triangulation(),
                                   neighbor_cell->level(),
                                   neighbor_cell->index(),
                                   &cell->get_dof_handler())
-                                  .get_dof_indices(neighbor_dofs);
+                                  .get_dof_indices(neighbor_dofs_all);
                                 // If the vector is distributed, we need to
                                 // transform the global indices to local ones.
                                 if (partitioner)
-                                  for (auto &index : neighbor_dofs)
+                                  for (auto &index : neighbor_dofs_all)
                                     index = partitioner->global_to_local(index);
 
+                                for (unsigned int i = 0;
+                                     i < neighbor_dofs_all_temp.size();
+                                     ++i)
+                                  neighbor_dofs_all_temp[i] =
+                                    neighbor_dofs_all[lexicographic_mapping[i]];
+
                                 for (unsigned int i = 0; i < n_dofs_1d; ++i)
                                   {
                                     // Get local dof index along line
@@ -680,14 +688,12 @@ namespace internal
                                       line_dof_idx(local_line, i, n_dofs_1d);
 
                                     dof_indices[idx + idx_offset[comp]] =
-                                      neighbor_dofs
-                                        [lexicographic_mapping
-                                           [fe.component_to_system_index(
-                                             comp,
-                                             line_dof_idx(
-                                               local_line_neighbor,
-                                               flipped ? fe_degree - i : i,
-                                               n_dofs_1d))]];
+                                      neighbor_dofs_all_temp
+                                        [line_dof_idx(local_line_neighbor,
+                                                      flipped ? fe_degree - i :
+                                                                i,
+                                                      n_dofs_1d) +
+                                         idx_offset[comp]];
                                   }
 
                                 // Stop looping over edge neighbors
diff --git a/tests/matrix_free/matrix_vector_stokes_noflux2.cc b/tests/matrix_free/matrix_vector_stokes_noflux2.cc
new file mode 100644 (file)
index 0000000..c3d2ea5
--- /dev/null
@@ -0,0 +1,340 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// matrix-vector products by comparing with the result of deal.II sparse
+// matrix. The implementation is similar to matrix_vector_stokes_noflux.cc
+// but tests edge hanging-node constraints.
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/block_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 <complex>
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, int degree_p, typename VectorType>
+class MatrixFreeTest
+{
+public:
+  using CellIterator = typename DoFHandler<dim>::active_cell_iterator;
+  using Number       = double;
+
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  void
+  local_apply(const MatrixFree<dim, Number> &              data,
+              VectorType &                                 dst,
+              const VectorType &                           src,
+              const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    using vector_t = VectorizedArray<Number>;
+    FEEvaluation<dim, degree_p + 1, degree_p + 2, dim, Number> velocity(data,
+                                                                        0);
+    FEEvaluation<dim, degree_p, degree_p + 2, 1, Number> pressure(data, 1);
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        velocity.reinit(cell);
+        velocity.read_dof_values(src.block(0));
+        velocity.evaluate(EvaluationFlags::gradients);
+        pressure.reinit(cell);
+        pressure.read_dof_values(src.block(1));
+        pressure.evaluate(EvaluationFlags::values);
+
+        for (unsigned int q = 0; q < velocity.n_q_points; ++q)
+          {
+            SymmetricTensor<2, dim, vector_t> sym_grad_u =
+              velocity.get_symmetric_gradient(q);
+            vector_t pres = pressure.get_value(q);
+            vector_t div  = -trace(sym_grad_u);
+            pressure.submit_value(div, q);
+
+            // subtract p * I
+            for (unsigned int d = 0; d < dim; ++d)
+              sym_grad_u[d][d] -= pres;
+
+            velocity.submit_symmetric_gradient(sym_grad_u, q);
+          }
+
+        velocity.integrate(EvaluationFlags::gradients);
+        velocity.distribute_local_to_global(dst.block(0));
+        pressure.integrate(EvaluationFlags::values);
+        pressure.distribute_local_to_global(dst.block(1));
+      }
+  }
+
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    dst = 0;
+    data.cell_loop(&MatrixFreeTest<dim, degree_p, VectorType>::local_apply,
+                   this,
+                   dst,
+                   src);
+  };
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation);
+  triangulation.refine_global();
+  for (const auto &cell : triangulation.active_cell_iterators())
+    if (cell->index() > 0)
+      cell->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement();
+
+  MappingQ<dim>   mapping(3);
+  FE_Q<dim>       fe_u_scal(fe_degree + 1);
+  FESystem<dim>   fe_u(fe_u_scal, dim);
+  FE_Q<dim>       fe_p(fe_degree);
+  FESystem<dim>   fe(fe_u_scal, dim, fe_p, 1);
+  DoFHandler<dim> dof_handler_u(triangulation);
+  DoFHandler<dim> dof_handler_p(triangulation);
+  DoFHandler<dim> dof_handler(triangulation);
+
+  MatrixFree<dim, double> mf_data;
+
+  AffineConstraints<double> constraints, constraints_u, constraints_p;
+
+  BlockSparsityPattern      sparsity_pattern;
+  BlockSparseMatrix<double> system_matrix;
+
+  BlockVector<double> solution;
+  BlockVector<double> system_rhs;
+  BlockVector<double> mf_solution;
+
+  dof_handler.distribute_dofs(fe);
+  dof_handler_u.distribute_dofs(fe_u);
+  dof_handler_p.distribute_dofs(fe_p);
+  std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
+  stokes_sub_blocks[dim] = 1;
+  DoFRenumbering::component_wise(dof_handler, stokes_sub_blocks);
+
+  std::set<types::boundary_id> no_normal_flux_boundaries;
+  no_normal_flux_boundaries.insert(0);
+  no_normal_flux_boundaries.insert(1);
+  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+  VectorTools::compute_no_normal_flux_constraints(
+    dof_handler, 0, no_normal_flux_boundaries, constraints, mapping);
+  constraints.close();
+  DoFTools::make_hanging_node_constraints(dof_handler_u, constraints_u);
+  VectorTools::compute_no_normal_flux_constraints(
+    dof_handler_u, 0, no_normal_flux_boundaries, constraints_u, mapping);
+  constraints_u.close();
+  DoFTools::make_hanging_node_constraints(dof_handler_p, constraints_p);
+  constraints_p.close();
+
+  const std::vector<types::global_dof_index> dofs_per_block =
+    DoFTools::count_dofs_per_fe_block(dof_handler, stokes_sub_blocks);
+
+  // std::cout << "Number of active cells: "
+  //          << triangulation.n_active_cells()
+  //          << std::endl
+  //          << "Number of degrees of freedom: "
+  //          << dof_handler.n_dofs()
+  //          << " (" << n_u << '+' << n_p << ')'
+  //          << std::endl;
+
+  {
+    BlockDynamicSparsityPattern csp(2, 2);
+
+    for (unsigned int d = 0; d < 2; ++d)
+      for (unsigned int e = 0; e < 2; ++e)
+        csp.block(d, e).reinit(dofs_per_block[d], dofs_per_block[e]);
+
+    csp.collect_sizes();
+
+    DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false);
+    sparsity_pattern.copy_from(csp);
+  }
+
+  system_matrix.reinit(sparsity_pattern);
+
+  // this is from step-22
+  {
+    QGauss<dim> quadrature_formula(fe_degree + 2);
+
+    FEValues<dim> fe_values(mapping,
+                            fe,
+                            quadrature_formula,
+                            update_values | update_JxW_values |
+                              update_gradients);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    const FEValuesExtractors::Vector velocities(0);
+    const FEValuesExtractors::Scalar pressure(dim);
+
+    std::vector<SymmetricTensor<2, dim>> phi_grads_u(dofs_per_cell);
+    std::vector<double>                  div_phi_u(dofs_per_cell);
+    std::vector<double>                  phi_p(dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator cell =
+                                                     dof_handler.begin_active(),
+                                                   endc = dof_handler.end();
+    for (; cell != endc; ++cell)
+      {
+        fe_values.reinit(cell);
+        local_matrix = 0;
+
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            for (unsigned int k = 0; k < dofs_per_cell; ++k)
+              {
+                phi_grads_u[k] = fe_values[velocities].symmetric_gradient(k, q);
+                div_phi_u[k]   = fe_values[velocities].divergence(k, q);
+                phi_p[k]       = fe_values[pressure].value(k, q);
+              }
+
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              {
+                for (unsigned int j = 0; j <= i; ++j)
+                  {
+                    local_matrix(i, j) +=
+                      (phi_grads_u[i] * phi_grads_u[j] -
+                       div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
+                      fe_values.JxW(q);
+                  }
+              }
+          }
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
+            local_matrix(i, j) = local_matrix(j, i);
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global(local_matrix,
+                                               local_dof_indices,
+                                               system_matrix);
+      }
+  }
+
+
+  solution.reinit(2);
+  for (unsigned int d = 0; d < 2; ++d)
+    solution.block(d).reinit(dofs_per_block[d]);
+  solution.collect_sizes();
+
+  system_rhs.reinit(solution);
+  mf_solution.reinit(solution);
+
+  // fill system_rhs with random numbers
+  for (unsigned int j = 0; j < system_rhs.block(0).size(); ++j)
+    if (constraints_u.is_constrained(j) == false)
+      {
+        const double val       = -1 + 2. * random_value<double>();
+        system_rhs.block(0)(j) = val;
+      }
+  for (unsigned int j = 0; j < system_rhs.block(1).size(); ++j)
+    if (constraints_p.is_constrained(j) == false)
+      {
+        const double val       = -1 + 2. * random_value<double>();
+        system_rhs.block(1)(j) = val;
+      }
+
+  // setup matrix-free structure
+  {
+    std::vector<const DoFHandler<dim> *> dofs;
+    dofs.push_back(&dof_handler_u);
+    dofs.push_back(&dof_handler_p);
+    std::vector<const AffineConstraints<double> *> constraints;
+    constraints.push_back(&constraints_u);
+    constraints.push_back(&constraints_p);
+    QGauss<1> quad(fe_degree + 2);
+    // no parallelism
+    mf_data.reinit(mapping,
+                   dofs,
+                   constraints,
+                   quad,
+                   typename MatrixFree<dim>::AdditionalData(
+                     MatrixFree<dim>::AdditionalData::none));
+  }
+
+  system_matrix.vmult(solution, system_rhs);
+
+  MatrixFreeTest<dim, fe_degree, BlockVector<double>> mf(mf_data);
+  mf.vmult(mf_solution, system_rhs);
+
+  // Verification
+  mf_solution -= solution;
+  const double error    = mf_solution.linfty_norm();
+  const double relative = solution.linfty_norm();
+  deallog << "Verification fe degree " << fe_degree << ": " << error / relative
+          << std::endl
+          << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(3);
+
+  {
+    deallog << std::endl << "Test with doubles" << std::endl << std::endl;
+    deallog.push("2d");
+    test<2, 1>();
+    test<2, 2>();
+    test<2, 3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_stokes_noflux2.output b/tests/matrix_free/matrix_vector_stokes_noflux2.output
new file mode 100644 (file)
index 0000000..0434059
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::
+DEAL::Test with doubles
+DEAL::
+DEAL:2d::Verification fe degree 1: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 2: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 3: 0
+DEAL:2d::
+DEAL:3d::Verification fe degree 1: 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.