From 3d39add1d98dda2218e907120535b31c973d81a4 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 18 Sep 2021 21:58:34 +0200 Subject: [PATCH] New test case --- tests/matrix_free/matrix_vector_faces_33.cc | 126 ++++++++++++++++++ .../matrix_free/matrix_vector_faces_33.output | 13 ++ 2 files changed, 139 insertions(+) create mode 100644 tests/matrix_free/matrix_vector_faces_33.cc create mode 100644 tests/matrix_free/matrix_vector_faces_33.output diff --git a/tests/matrix_free/matrix_vector_faces_33.cc b/tests/matrix_free/matrix_vector_faces_33.cc new file mode 100644 index 0000000000..c792e5ff7d --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_33.cc @@ -0,0 +1,126 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2020 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. +// +// --------------------------------------------------------------------- + + + +// tests matrix-free face evaluation when only boundary faces but not inner +// faces are enabled. Otherwise the same test as matrix_vector_faces_05 (FE_Q, +// Laplacian, weak imposition of Dirichlet boundary condition) + +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_faces_common.h" + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(5 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + // std::cout << "Number of cells: " << + // dof.get_triangulation().n_active_cells() << std::endl; std::cout << "Number + // of degrees of freedom: " << dof.n_dofs() << std::endl; std::cout << "Number + // of constraints: " << constraints.n_constraints() << std::endl; + + MappingQ mapping(dof.get_fe().degree + 1); + + Vector in(dof.n_dofs()), out(dof.n_dofs()); + Vector out_dist(out); + + // Set random seed for reproducibility + Testing::srand(42); + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + if (constraints.is_constrained(i)) + continue; + const double entry = Testing::rand() / (double)RAND_MAX; + in(i) = entry; + } + + constexpr unsigned int n_q_points_1d = fe_degree + 1; + + // assemble sparse matrix with MeshWorker + SparsityPattern sparsity; + SparseMatrix matrix; + { + DynamicSparsityPattern d_sparsity(dof.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dof, d_sparsity); + sparsity.copy_from(d_sparsity); + } + matrix.reinit(sparsity); + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = + update_values | update_gradients | update_jacobians; + info_box.add_update_flags_all(update_flags); + info_box.initialize_gauss_quadrature(n_q_points_1d, + n_q_points_1d, + n_q_points_1d); + info_box.initialize(dof.get_fe(), mapping); + + MeshWorker::DoFInfo dof_info(dof); + + MeshWorker::Assembler::MatrixSimple> assembler; + assembler.initialize(matrix); + + MatrixIntegrator integrator; + MeshWorker::integration_loop( + dof.begin_active(), dof.end(), dof_info, info_box, integrator, assembler); + + matrix.vmult(out, in); + + // zero constrained dofs + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + if (constraints.is_constrained(i)) + out(i) = 0; + + MatrixFree mf_data; + const QGauss<1> quad(n_q_points_1d > 0 ? n_q_points_1d : + dof.get_fe().degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 3; + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + mf_data.reinit(mapping, dof, constraints, quad, data); + + // Check that only boundary faces are set up as requested + Assert(mf_data.n_inner_face_batches() == 0, ExcInternalError()); + Assert(mf_data.n_boundary_face_batches() > 0, ExcInternalError()); + + MatrixFreeTest, 1> mf( + mf_data); + mf.vmult(out_dist, in); + + out_dist -= out; + const double diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << std::endl; + + deallog << std::endl; +} diff --git a/tests/matrix_free/matrix_vector_faces_33.output b/tests/matrix_free/matrix_vector_faces_33.output new file mode 100644 index 0000000000..52f68efa16 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_33.output @@ -0,0 +1,13 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: -- 2.39.5