From: Sebastian Kinnewig <kinnewig@ifam.uni-hannover.de>
Date: Mon, 15 Jan 2024 16:18:27 +0000 (+0100)
Subject: Add another test to verify that inhomogeneous boundary conditions with hanging nodes... 
X-Git-Tag: relicensing~125^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f16f94d2ff712d15f9675240ca75c37c291617e3;p=dealii.git

Add another test to verify that inhomogeneous boundary conditions with hanging nodes work as expected.
---

diff --git a/tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.cc b/tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.cc
new file mode 100644
index 0000000000..21e640703c
--- /dev/null
+++ b/tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.cc
@@ -0,0 +1,413 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 - 2024 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 convergence test verifies that hanging nodes on FE_NedelecSZ
+// elements are handled correctly. The orientation of edges and faces is
+// automatically adapted in the presence of hanging edges and hanging
+// faces. Furthermore, the function make_hanging_node_constraints()
+// considers the orientation of the underlying edges and faces to prevent
+// sign conflicts.
+//
+// This test solves the real-valued Maxwell equation in 2D and 3D:
+//
+// curl(curl(E)) - E = 0,
+//
+// where we consider Dirichlet boundary data
+// such that we know the solution E.
+// In 2D:
+// E = (sin(y), sin(x))
+//
+// In 3D:
+// E = (sin(y), sin(z), 0)
+//
+// In the first step, we solve the Maxwell equation on a coarse grid.
+// After that, we refine the center of the domain so the grid contains
+// hanging edges (and faces in 3D). We expect that the L2 difference
+// between the numerical solution and the exact solution stays the same
+// or gets smaller in each step, i.e., with a finer and finer grid, we
+// converge more and more to the exact solution.
+// The geometry is chosen in such a way that it covers the most important
+// cases and is as small as possible.
+// To cover all aspects of make_hanging_node_constraints(), we have to
+// consider at least cubic base functions (that corresponds to
+// polynomial_degree = 2).
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nedelec_sz.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/vector_tools_boundary.h>
+
+#include <cmath>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace ConvergenceTest
+{
+
+  template <int dim>
+  class SolutionValues : public Function<dim>
+  {
+  public:
+    SolutionValues()
+      : Function<dim>(dim)
+    {}
+
+    void
+    vector_value(const Point<dim> &p, Vector<double> &values) const override;
+
+    void
+    vector_value_list(const std::vector<Point<dim>> &points,
+                      std::vector<Vector<double>>   &value_list) const override
+    {
+      Assert(value_list.size() == points.size(),
+             ExcDimensionMismatch(value_list.size(), points.size()));
+
+      for (unsigned int p = 0; p < points.size(); p++)
+        {
+          SolutionValues<dim>::vector_value(points[p], value_list[p]);
+        }
+    }
+  };
+
+  // here we are using artificial boundary values,
+  // the benefit is, that we know the exact solution for the
+  // electric field E, i.e. this is also the exact solution
+  // to the partial differential equation we aim to solve
+  template <>
+  void
+  SolutionValues<2>::vector_value(const Point<2> &p,
+                                  Vector<double> &values) const
+  {
+    values(0) = sin(p[1]);
+    values(1) = sin(p[0]);
+  }
+
+  template <>
+  void
+  SolutionValues<3>::vector_value(const Point<3> &p,
+                                  Vector<double> &values) const
+  {
+    values(0) = sin(p[1]);
+    values(1) = sin(p[2]);
+    values(2) = 0.0;
+  }
+
+
+
+  template <int dim>
+  class MaxwellProblem
+  {
+  public:
+    MaxwellProblem(unsigned int refinements, unsigned int poly_degree);
+
+    void
+    run();
+
+  private:
+    void
+    make_grid();
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid(double radius);
+    double
+    output_error();
+
+    Triangulation<dim> triangulation;
+    DoFHandler<dim>    dof_handler;
+    FE_NedelecSZ<dim>  fe;
+
+    AffineConstraints<double> constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double>       solution, system_rhs;
+
+    const unsigned int refinements;
+    const unsigned int poly_degree;
+  };
+
+  template <int dim>
+  MaxwellProblem<dim>::MaxwellProblem(const unsigned int refinements,
+                                      const unsigned int poly_degree)
+    : dof_handler(triangulation)
+    , fe(poly_degree)
+    , refinements(refinements)
+    , poly_degree(poly_degree)
+  {}
+
+  template <int dim>
+  void
+  MaxwellProblem<dim>::make_grid()
+  {
+    GridGenerator::hyper_ball_balanced(triangulation);
+
+    static const SphericalManifold<dim> spherical_manifold;
+    triangulation.set_manifold(0, spherical_manifold);
+
+    triangulation.refine_global(refinements);
+
+    // Refine one cell at the boundary,
+    // such that the resulting grid contains
+    // hanging nodes.
+    for (auto &cell : triangulation.cell_iterators())
+      {
+        if (!cell->is_active())
+          continue;
+
+        if (cell->at_boundary())
+          {
+            cell->set_refine_flag();
+            break;
+          }
+      }
+    triangulation.execute_coarsening_and_refinement();
+  }
+
+  template <int dim>
+  void
+  MaxwellProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp);
+
+    constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+    VectorTools::project_boundary_values_curl_conforming_l2(
+      dof_handler,
+      0 /* vector component*/,
+      SolutionValues<dim>(),
+      0 /* boundary id*/,
+      constraints);
+    constraints.close();
+
+    constraints.condense(dsp);
+
+    sparsity_pattern.copy_from(dsp);
+    system_matrix.reinit(sparsity_pattern);
+    system_rhs.reinit(dof_handler.n_dofs());
+    solution.reinit(dof_handler.n_dofs());
+  }
+
+  template <int dim>
+  void
+  MaxwellProblem<dim>::assemble_system()
+  {
+    system_matrix = 0;
+    system_rhs    = 0;
+
+    const unsigned int curl_dim = (dim == 2) ? 1 : 3;
+
+    // choose the quadrature formulas
+    QGauss<dim> quadrature_formula(fe.degree + 2);
+
+    // get the number of quadrature points and dofs
+    const unsigned int n_q_points    = quadrature_formula.size(),
+                       dofs_per_cell = fe.dofs_per_cell;
+
+    // set update flags
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    // Extractors for the real part
+    const FEValuesExtractors::Vector E_re(0);
+
+    // create the local left hand side and right hand side
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    // loop over all cells
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        if (cell->is_locally_owned() == false)
+          continue;
+
+        // initialize values:
+        cell_matrix = 0;
+        cell_rhs    = 0;
+        fe_values.reinit(cell);
+
+        for (const unsigned int i : fe_values.dof_indices())
+          {
+            // only compute this once
+            std::vector<Tensor<1, dim>>      phi_i(n_q_points);
+            std::vector<Tensor<1, curl_dim>> curl_phi_i(n_q_points);
+            for (unsigned int q_point = 0; q_point < n_q_points; q_point++)
+              {
+                phi_i[q_point]      = fe_values[E_re].value(i, q_point);
+                curl_phi_i[q_point] = fe_values[E_re].curl(i, q_point);
+              }
+
+            // we use here, that the problem is symmetrical
+            for (unsigned int j = i; j < dofs_per_cell; j++)
+              {
+                double mass_part = 0;
+                double curl_part = 0;
+
+                for (unsigned int q_point = 0; q_point < n_q_points; q_point++)
+                  {
+                    Tensor<1, dim> phi_j = fe_values[E_re].value(j, q_point);
+                    Tensor<1, curl_dim> curl_phi_j =
+                      fe_values[E_re].curl(j, q_point);
+
+                    curl_part +=
+                      curl_phi_i[q_point] * curl_phi_j * fe_values.JxW(q_point);
+
+                    mass_part +=
+                      phi_i[q_point] * phi_j * fe_values.JxW(q_point);
+                  }
+
+                double mass_term  = curl_part - mass_part;
+                cell_matrix(i, j) = mass_term;
+                cell_matrix(j, i) = mass_term;
+              }
+          }
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global(
+          cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+      }
+  }
+
+
+  template <int dim>
+  void
+  MaxwellProblem<dim>::solve()
+  {
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(system_matrix);
+    A_direct.vmult(solution, system_rhs);
+    constraints.distribute(solution);
+  }
+
+  template <int dim>
+  void
+  MaxwellProblem<dim>::refine_grid(const double radius)
+  {
+    for (auto &cell : triangulation.cell_iterators())
+      {
+        if (!cell->is_active())
+          continue;
+
+        double distance;
+        if (dim == 3)
+          distance = std::sqrt(std::pow(cell->center()[1], 2) +
+                               std::pow(cell->center()[2], 2));
+        else if (dim == 2)
+          distance = std::sqrt(std::pow(cell->center()[0], 2) +
+                               std::pow(cell->center()[1], 2));
+
+        if (distance < radius)
+          cell->set_refine_flag();
+      }
+
+    triangulation.execute_coarsening_and_refinement();
+  }
+
+  template <int dim>
+  double
+  MaxwellProblem<dim>::output_error()
+  {
+    SolutionValues<dim> exact_solution;
+    Vector<double>      diff_per_cell(triangulation.n_active_cells());
+    VectorTools::integrate_difference(dof_handler,
+                                      solution,
+                                      exact_solution,
+                                      diff_per_cell,
+                                      QGauss<dim>(poly_degree + 2),
+                                      VectorTools::L2_norm);
+    return diff_per_cell.l2_norm();
+  }
+
+  template <int dim>
+  void
+  MaxwellProblem<dim>::run()
+  {
+    deallog << "Testing for dim = " << dim
+            << ", polynomial_degree p = " << poly_degree << std::endl;
+    make_grid();
+    setup_system();
+    assemble_system();
+    solve();
+
+    // compute the error
+    double L2_error = output_error();
+
+    // check that we are close to the analytical solution
+    bool passed = (std::fabs(L2_error) < 1e-2);
+
+    if (passed)
+      deallog << "OK" << std::endl;
+    else
+      {
+        deallog << "FAILED" << std::endl;
+        deallog << "L2 Error = " << L2_error << "(threshold = 1e-2)"
+                << std::endl;
+      }
+  }
+} // namespace ConvergenceTest
+
+int
+main()
+{
+  initlog();
+
+  using namespace ConvergenceTest;
+
+  const unsigned int refinements = 1;
+  const unsigned int poly_degree = 2;
+
+  MaxwellProblem<2> maxwell_2d(refinements, poly_degree);
+  maxwell_2d.run();
+
+  MaxwellProblem<3> maxwell_3d(refinements, poly_degree);
+  maxwell_3d.run();
+
+  return 0;
+}
diff --git a/tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.output b/tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.output
new file mode 100644
index 0000000000..8a79e42d27
--- /dev/null
+++ b/tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.output
@@ -0,0 +1,5 @@
+
+DEAL::Testing for dim = 2, polynomial_degree p = 2
+DEAL::OK
+DEAL::Testing for dim = 3, polynomial_degree p = 2
+DEAL::OK