From 1a935f405dae11aa293bc0aa32f44ea6d442ce5f Mon Sep 17 00:00:00 2001 From: Niklas Wik Date: Thu, 16 Jun 2022 15:15:18 +0200 Subject: [PATCH] Tests for hanging nodes --- tests/matrix_free/matrix_vector_rt_03.cc | 74 +++++++++++++++++++ tests/matrix_free/matrix_vector_rt_03.output | 29 ++++++++ tests/matrix_free/matrix_vector_rt_common.h | 24 +++--- tests/matrix_free/matrix_vector_rt_face_03.cc | 70 ++++++++++++++++++ .../matrix_vector_rt_face_03.output | 22 ++++++ .../matrix_vector_rt_face_common.h | 37 ++++++---- 6 files changed, 230 insertions(+), 26 deletions(-) create mode 100644 tests/matrix_free/matrix_vector_rt_03.cc create mode 100644 tests/matrix_free/matrix_vector_rt_03.output create mode 100644 tests/matrix_free/matrix_vector_rt_face_03.cc create mode 100644 tests/matrix_free/matrix_vector_rt_face_03.output diff --git a/tests/matrix_free/matrix_vector_rt_03.cc b/tests/matrix_free/matrix_vector_rt_03.cc new file mode 100644 index 0000000000..d817e027a7 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_03.cc @@ -0,0 +1,74 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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 test is the same as matrix_vector_rt_01.cc but with hanging nodes. + +#include + +#include "../tests.h" + +#include "matrix_vector_rt_common.h" + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + if (dim < 3 || fe_degree < 2) + tria.refine_global(2); + else + tria.refine_global(1); + typename Triangulation::active_cell_iterator cell, endc; + cell = tria.begin_active(), endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + cell = tria.begin_active(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 0.2) + 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(); + cell = tria.begin_active(); + for (unsigned int i = 0; i < 10 - 3 * dim; ++i) + { + cell = tria.begin_active(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_RaviartThomasNodal fe(fe_degree - 1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + constraints.close(); + + deallog << "Using " << dof.get_fe().get_name() << std::endl; + deallog << "Number of cells: " << dof.get_triangulation().n_active_cells() + << std::endl; + deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl + << std::endl; + do_test(dof, constraints, TestType::values_gradients); +} diff --git a/tests/matrix_free/matrix_vector_rt_03.output b/tests/matrix_free/matrix_vector_rt_03.output new file mode 100644 index 0000000000..f327e829b4 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_03.output @@ -0,0 +1,29 @@ + +DEAL:2d::Using FE_RaviartThomasNodal<2>(1) +DEAL:2d::Number of cells: 481 +DEAL:2d::Number of degrees of freedom: 4548 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 1.13503e-13 +DEAL:2d:: +DEAL:2d::Using FE_RaviartThomasNodal<2>(2) +DEAL:2d::Number of cells: 481 +DEAL:2d::Number of degrees of freedom: 9708 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 2.05175e-13 +DEAL:2d:: +DEAL:3d::Using FE_RaviartThomasNodal<3>(1) +DEAL:3d::Number of cells: 85 +DEAL:3d::Number of degrees of freedom: 2392 +DEAL:3d:: +DEAL:3d::Testing Values and Gradients +DEAL:3d::Norm of difference: 1.98105e-14 +DEAL:3d:: +DEAL:3d::Using FE_RaviartThomasNodal<3>(2) +DEAL:3d::Number of cells: 85 +DEAL:3d::Number of degrees of freedom: 7677 +DEAL:3d:: +DEAL:3d::Testing Values and Gradients +DEAL:3d::Norm of difference: 3.46605e-14 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_rt_common.h b/tests/matrix_free/matrix_vector_rt_common.h index 90ea0b41c8..40d7022da1 100644 --- a/tests/matrix_free/matrix_vector_rt_common.h +++ b/tests/matrix_free/matrix_vector_rt_common.h @@ -173,10 +173,11 @@ do_test(const DoFHandler & dof, MatrixFree mf_data; { const QGaussLobatto<1> quad(fe_degree + 2); + const MappingQ mapping(fe_degree + 2); typename MatrixFree::AdditionalData data; data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; data.mapping_update_flags = update_gradients | update_hessians; - mf_data.reinit(dof, constraints, quad, data); + mf_data.reinit(mapping, dof, constraints, quad, data); } // create vector with random entries @@ -191,23 +192,27 @@ do_test(const DoFHandler & dof, continue; initial_condition[i] = random_value(); } + constraints.distribute(initial_condition); // MatrixFree solution MatrixFreeTest mf(mf_data, test_type); mf.test_functions(solution, initial_condition); + // Evaluation with FEValues SparsityPattern sp; SparseMatrix system_matrix; DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs()); - DoFTools::make_sparsity_pattern(dof, dsp); + DoFTools::make_sparsity_pattern(dof, dsp, constraints, false); sp.copy_from(dsp); system_matrix.reinit(sp); - FEValues fe_val(dof.get_fe(), + const MappingQ mapping(fe_degree + 2); + FEValues fe_val(mapping, + dof.get_fe(), Quadrature(mf_data.get_quadrature(0)), - update_values | update_gradients | update_JxW_values); - + update_values | update_gradients | update_piola | + update_JxW_values); const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell; std::vector local_dof_indices(dofs_per_cell); @@ -248,17 +253,16 @@ do_test(const DoFHandler & dof, } } cell->get_dof_indices(local_dof_indices); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - system_matrix.add(local_dof_indices[i], - local_dof_indices[j], - local_matrix(i, j)); + constraints.distribute_local_to_global(local_matrix, + local_dof_indices, + system_matrix); } Vector ref(solution.size()); // Compute reference system_matrix.vmult(ref, initial_condition); + constraints.set_zero(ref); ref -= solution; diff --git a/tests/matrix_free/matrix_vector_rt_face_03.cc b/tests/matrix_free/matrix_vector_rt_face_03.cc new file mode 100644 index 0000000000..66736efe8b --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_03.cc @@ -0,0 +1,70 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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 test does the same as matrix_vector_rt_face_01.cc but also uses hanging +// nodes. + +#include + +#include "../tests.h" + +#include "matrix_vector_rt_face_common.h" + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + typename Triangulation::active_cell_iterator cell, endc; + cell = tria.begin_active(), endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + cell = tria.begin_active(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 0.2) + 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(); + cell = tria.begin_active(); + for (unsigned int i = 0; i < 10 - 3 * dim; ++i) + { + cell = tria.begin_active(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_RaviartThomasNodal fe(fe_degree - 1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + constraints.close(); + + deallog << "Using " << dof.get_fe().get_name() << std::endl; + deallog << "Number of cells: " << tria.n_active_cells() << std::endl; + deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl + << std::endl; + do_test(dof, constraints, TestType::values_gradients); + // do_test(dof, constraints, TestType::divergence); +} diff --git a/tests/matrix_free/matrix_vector_rt_face_03.output b/tests/matrix_free/matrix_vector_rt_face_03.output new file mode 100644 index 0000000000..31d0141aee --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_03.output @@ -0,0 +1,22 @@ + +DEAL:2d::Using FE_RaviartThomasNodal<2>(1) +DEAL:2d::Number of cells: 85 +DEAL:2d::Number of degrees of freedom: 798 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 1.46726e-13 +DEAL:2d:: +DEAL:2d::Using FE_RaviartThomasNodal<2>(2) +DEAL:2d::Number of cells: 85 +DEAL:2d::Number of degrees of freedom: 1707 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 1.47883e-13 +DEAL:2d:: +DEAL:3d::Using FE_RaviartThomasNodal<3>(1) +DEAL:3d::Number of cells: 22 +DEAL:3d::Number of degrees of freedom: 672 +DEAL:3d:: +DEAL:3d::Testing Values and Gradients +DEAL:3d::Norm of difference: 3.31436e-14 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_rt_face_common.h b/tests/matrix_free/matrix_vector_rt_face_common.h index dd65b7f391..39a3ffa46d 100644 --- a/tests/matrix_free/matrix_vector_rt_face_common.h +++ b/tests/matrix_free/matrix_vector_rt_face_common.h @@ -225,17 +225,21 @@ do_test(const DoFHandler & dof, { deallog << "Testing " << enum_to_string(test_type) << std::endl; + constexpr unsigned int n_q_points = fe_degree + 2; + constexpr unsigned int m_degree = fe_degree + 2; + + const MappingQ mapping(m_degree); MatrixFree mf_data; { - const QGaussLobatto<1> quad(fe_degree + 2); - const MappingQ mapping(fe_degree); + const QGaussLobatto<1> quad(n_q_points); typename MatrixFree::AdditionalData data; data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; - data.mapping_update_flags = update_gradients | update_JxW_values; + data.mapping_update_flags = + (update_gradients | update_JxW_values | update_hessians); data.mapping_update_flags_inner_faces = - (update_gradients | update_JxW_values); + (update_gradients | update_JxW_values | update_hessians); data.mapping_update_flags_boundary_faces = - (update_gradients | update_JxW_values); + (update_gradients | update_JxW_values | update_hessians); mf_data.reinit(mapping, dof, constraints, quad, data); } @@ -251,23 +255,25 @@ do_test(const DoFHandler & dof, continue; initial_condition[i] = random_value(); } + constraints.distribute(initial_condition); MatrixFreeTest mf(mf_data, test_type); mf.test_functions(solution, initial_condition); + // Evaluation with FEFaceValues SparsityPattern sp; SparseMatrix system_matrix; DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs()); - DoFTools::make_sparsity_pattern(dof, dsp); + DoFTools::make_sparsity_pattern(dof, dsp, constraints, false); sp.copy_from(dsp); system_matrix.reinit(sp); - FEFaceValues fe_val( - dof.get_fe(), - QGaussLobatto(mf_data.get_quadrature(0).size()), - update_values | update_gradients | update_JxW_values); - + FEFaceValues fe_val(mapping, + dof.get_fe(), + QGaussLobatto(n_q_points), + update_values | update_gradients | + update_JxW_values | update_piola); const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell; std::vector local_dof_indices(dofs_per_cell); @@ -308,17 +314,16 @@ do_test(const DoFHandler & dof, } } cell->get_dof_indices(local_dof_indices); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - system_matrix.add(local_dof_indices[i], - local_dof_indices[j], - local_matrix(i, j)); + constraints.distribute_local_to_global(local_matrix, + local_dof_indices, + system_matrix); } Vector ref(solution.size()); // Compute reference system_matrix.vmult(ref, initial_condition); + constraints.set_zero(ref); ref -= solution; -- 2.39.5