From dfd5f531b6f1e2b1376fa215c53a5dcf39a8f6e2 Mon Sep 17 00:00:00 2001 From: Sebastian Kinnewig Date: Fri, 8 Dec 2023 12:28:47 +0100 Subject: [PATCH] Add a convergence test that verifies that hanging nodes on FE_NedelecSZ elements are handled correctly. --- source/dofs/dof_tools_constraints.cc | 29 +- ...fe_nedelec_sz_hanging_nodes_convergence.cc | 447 ++++++++++++++++++ ...edelec_sz_hanging_nodes_convergence.output | 5 + 3 files changed, 467 insertions(+), 14 deletions(-) create mode 100644 tests/fe/fe_nedelec_sz_hanging_nodes_convergence.cc create mode 100644 tests/fe/fe_nedelec_sz_hanging_nodes_convergence.output diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index bb9907c83a..92f101f19d 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -1594,19 +1594,19 @@ namespace DoFTools // internal lines for (unsigned int line = 0; line < n_internal_lines_on_children; ++line) - { - unsigned int row_start = line * fe.n_dofs_per_line(); - unsigned int line_mother = line / 2; - unsigned int row_mother = - (line_mother * 2) * fe.n_dofs_per_line(); - for (unsigned int row = 0; row < fe.n_dofs_per_line(); - ++row) - for (unsigned int i = 0; - i < n_lines_on_mother * fe.n_dofs_per_line(); - ++i) - constraints_matrix[row + row_start][i] = - fe.constraints()(row + row_mother, i); - } + { + unsigned int row_start = line * fe.n_dofs_per_line(); + unsigned int line_mother = line / 2; + unsigned int row_mother = + (line_mother * 2) * fe.n_dofs_per_line(); + for (unsigned int row = 0; row < fe.n_dofs_per_line(); + ++row) + for (unsigned int i = 0; + i < n_lines_on_mother * fe.n_dofs_per_line(); + ++i) + constraints_matrix[row + row_start][i] = + fe.constraints()(row + row_mother, i); + } for (unsigned int line = 0; line < n_internal_lines_on_children; ++line) @@ -1617,7 +1617,8 @@ namespace DoFTools (line_mother * 2) * fe.n_dofs_per_line(); for (unsigned int row = 0; row < fe.n_dofs_per_line(); ++row) - for (unsigned int i = n_lines_on_mother * fe.n_dofs_per_line(); + for (unsigned int i = + n_lines_on_mother * fe.n_dofs_per_line(); i < dofs_on_mother.size(); ++i) constraints_matrix[row + row_start][i] = diff --git a/tests/fe/fe_nedelec_sz_hanging_nodes_convergence.cc b/tests/fe/fe_nedelec_sz_hanging_nodes_convergence.cc new file mode 100644 index 0000000000..43615783f0 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_hanging_nodes_convergence.cc @@ -0,0 +1,447 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +namespace ConvergenceTest +{ + + template + class SolutionValues : public Function + { + public: + SolutionValues() + : Function(dim) + {} + + void + vector_value(const Point &p, Vector &values) const override; + + void + vector_value_list(const std::vector> &points, + std::vector> &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::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 &values) const + { + values(0) = sin(p[1]); + values(1) = sin(p[0]); + } + + template <> + void + SolutionValues<3>::vector_value(const Point<3> &p, + Vector &values) const + { + values(0) = sin(p[1]); + values(1) = sin(p[2]); + values(2) = 0.0; + } + + + + template + class MaxwellProblem + { + public: + MaxwellProblem(unsigned int refinements, + unsigned int poly_degree, + unsigned int n_iterations); + + void + run(); + + private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(double radius); + double + output_error(); + + Triangulation triangulation; + DoFHandler dof_handler; + FE_NedelecSZ fe; + + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution, system_rhs; + + const unsigned int refinements; + const unsigned int poly_degree; + const unsigned int n_iterations; + }; + + template + MaxwellProblem::MaxwellProblem(const unsigned int refinements, + const unsigned int poly_degree, + const unsigned int n_iterations) + : dof_handler(triangulation) + , fe(poly_degree) + , refinements(refinements) + , poly_degree(poly_degree) + , n_iterations(n_iterations) + {} + + template <> + void + MaxwellProblem<2>::make_grid() + { + const unsigned int dim = 2; + + const Point center(0.0, 0.0); + const double radius = 1.00; + + GridGenerator::hyper_ball_balanced(triangulation, center, radius); + + triangulation.refine_global(refinements); + } + + template <> + void + MaxwellProblem<3>::make_grid() + { + const double radius = 1.00; + const double x_thickness = 0.75; + + GridGenerator::cylinder(triangulation, radius, x_thickness / 2.0); + + for (auto &cell : triangulation.active_cell_iterators()) + for (auto &face : cell->face_iterators()) + if (face->at_boundary()) + face->set_boundary_id(0); + + triangulation.refine_global(refinements); + } + + template + void + MaxwellProblem::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(), + 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 + void + MaxwellProblem::assemble_system() + { + system_matrix = 0; + system_rhs = 0; + + const unsigned int curl_dim = (dim == 2) ? 1 : 3; + + // choose the quadrature formulas + QGauss 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 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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + std::vector 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> phi_i(n_q_points); + std::vector> 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 + void + MaxwellProblem::solve() + { + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); + constraints.distribute(solution); + } + + template + void + MaxwellProblem::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 + double + MaxwellProblem::output_error() + { + SolutionValues exact_solution; + Vector diff_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(dof_handler, + solution, + exact_solution, + diff_per_cell, + QGauss(poly_degree + 2), + VectorTools::L2_norm); + return diff_per_cell.l2_norm(); + } + + template + void + MaxwellProblem::run() + { + deallog << "Testing for dim = " << dim + << ", polynomial_degree p = " << poly_degree << std::endl; + std::vector L2_error(n_iterations); + bool passed = false; + + for (unsigned int cycle = 0; cycle < n_iterations; ++cycle) + { + if (cycle == 0) + make_grid(); + else + refine_grid(0.25); + + setup_system(); + assemble_system(); + solve(); + L2_error[cycle] = output_error(); + + // check that we are close to the analytical solution + passed = (std::fabs(L2_error[cycle]) < 1e-3); + + if (!passed) + { + deallog << "FAILED" << std::endl; + deallog << "L2 Error = " << L2_error << "(threshold = 1e-3)" + << std::endl; + } + + + // check for convergence + if (cycle > 0) + { + passed = (L2_error[cycle - 1] >= L2_error[cycle]); + + if (!passed) + { + deallog << "FAILED" << std::endl; + deallog << "Convergence error in step " << cycle << std::endl; + } + } + } + + if (passed) + deallog << "OK" << std::endl; + } +} // namespace ConvergenceTest + +int +main() +{ + initlog(); + + using namespace ConvergenceTest; + + const unsigned int refinements = 1; + const unsigned int poly_degree = 2; +#ifdef DEBUG + const unsigned int n_iterations = 2; +#else + const unsigned int n_iterations = 3; +#endif + + MaxwellProblem<2> maxwell_2d(refinements, poly_degree, n_iterations); + maxwell_2d.run(); + + MaxwellProblem<3> maxwell_3d(refinements, poly_degree, n_iterations); + maxwell_3d.run(); + + return 0; +} diff --git a/tests/fe/fe_nedelec_sz_hanging_nodes_convergence.output b/tests/fe/fe_nedelec_sz_hanging_nodes_convergence.output new file mode 100644 index 0000000000..8a79e42d27 --- /dev/null +++ b/tests/fe/fe_nedelec_sz_hanging_nodes_convergence.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 -- 2.39.5