From f16f94d2ff712d15f9675240ca75c37c291617e3 Mon Sep 17 00:00:00 2001 From: Sebastian Kinnewig Date: Mon, 15 Jan 2024 17:18:27 +0100 Subject: [PATCH] Add another test to verify that inhomogeneous boundary conditions with hanging nodes work as expected. --- ...hanging_nodes_on_inhomogeneous_boundary.cc | 413 ++++++++++++++++++ ...ing_nodes_on_inhomogeneous_boundary.output | 5 + 2 files changed, 418 insertions(+) create mode 100644 tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.cc create mode 100644 tests/fe/fe_nedelec_sz_hanging_nodes_on_inhomogeneous_boundary.output 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 +#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); + + 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; + }; + + template + MaxwellProblem::MaxwellProblem(const unsigned int refinements, + const unsigned int poly_degree) + : dof_handler(triangulation) + , fe(poly_degree) + , refinements(refinements) + , poly_degree(poly_degree) + {} + + template + void + MaxwellProblem::make_grid() + { + GridGenerator::hyper_ball_balanced(triangulation); + + static const SphericalManifold 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 + 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; + 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 -- 2.39.5