From ab83714d1218572fe8606da057991945b8025765 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Fri, 7 Jan 2022 22:31:25 +0100 Subject: [PATCH] Add some more tests --- tests/meshworker/scratch_data_02_hp.cc | 180 ++++++++++++ ...hp.with_muparser=on.with_umfpack=on.output | 2 + tests/meshworker/scratch_data_03_hp.cc | 268 ++++++++++++++++++ ...hp.with_muparser=on.with_umfpack=on.output | 2 + 4 files changed, 452 insertions(+) create mode 100644 tests/meshworker/scratch_data_02_hp.cc create mode 100644 tests/meshworker/scratch_data_02_hp.with_muparser=on.with_umfpack=on.output create mode 100644 tests/meshworker/scratch_data_03_hp.cc create mode 100644 tests/meshworker/scratch_data_03_hp.with_muparser=on.with_umfpack=on.output diff --git a/tests/meshworker/scratch_data_02_hp.cc b/tests/meshworker/scratch_data_02_hp.cc new file mode 100644 index 0000000000..4a33aee1c1 --- /dev/null +++ b/tests/meshworker/scratch_data_02_hp.cc @@ -0,0 +1,180 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2018 - 2021 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. + * + * --------------------------------------------------------------------- + */ + +// Check mesh_loop + ScratchData for Laplacian (conforming) +// - hp variant + +#include +#include +#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" + +using namespace MeshWorker; + +template +void +test() +{ + Triangulation tria; + hp::FECollection fe(FE_Q(1), + FE_Q(2)); + DoFHandler dh(tria); + + FunctionParser rhs_function("1"); + FunctionParser boundary_function("0"); + + AffineConstraints constraints; + + GridGenerator::hyper_cube(tria); + tria.refine_global(5); + tria.execute_coarsening_and_refinement(); + + const auto bounding_box = GridTools::compute_bounding_box(tria); + for (auto &cell : dh.active_cell_iterators()) + { + bool use_alt_index = true; + for (unsigned int d = 0; d < spacedim; ++d) + if (cell->center()[d] < bounding_box.center()[d]) + { + use_alt_index = false; + break; + } + + if (use_alt_index) + cell->set_active_fe_index(1); + else + cell->set_active_fe_index(0); + } + + dh.distribute_dofs(fe); + + DoFTools::make_hanging_node_constraints(dh, constraints); + VectorTools::interpolate_boundary_values(dh, + 0, + boundary_function, + constraints); + + constraints.close(); + + SparsityPattern sparsity; + + { + DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dh, dsp); + sparsity.copy_from(dsp); + } + + SparseMatrix matrix; + matrix.reinit(sparsity); + + Vector solution(dh.n_dofs()); + Vector rhs(dh.n_dofs()); + + hp::QCollection cell_quad(QGauss(3), QGauss(4)); + hp::QCollection face_quad(QGauss(3), QGauss(4)); + + UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + UpdateFlags face_flags = update_JxW_values; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<1, 1, 1>; + + ScratchData scratch(fe, cell_quad, cell_flags, face_quad, face_flags); + CopyData copy(0); + + auto cell = dh.begin_active(); + auto endc = dh.end(); + + using Iterator = decltype(cell); + + auto cell_worker = + [&rhs_function](const Iterator &cell, ScratchData &s, CopyData &c) { + const auto &fev = s.reinit(cell); + const auto &JxW = s.get_JxW_values(); + const auto &p = s.get_quadrature_points(); + + c.reinit(s.n_dofs_per_cell()); + c.local_dof_indices[0] = s.get_local_dof_indices(); + + for (unsigned int q = 0; q < p.size(); ++q) + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + { + c.matrices[0](i, j) += + fev.shape_grad(i, q) * fev.shape_grad(j, q) * JxW[q]; + } + c.vectors[0](i) += + fev.shape_value(i, q) * rhs_function.value(p[q]) * JxW[q]; + } + }; + + auto copier = [&constraints, &matrix, &rhs](const CopyData &c) { + constraints.distribute_local_to_global( + c.matrices[0], c.vectors[0], c.local_dof_indices[0], matrix, rhs); + }; + + mesh_loop(cell, endc, cell_worker, copier, scratch, copy, assemble_own_cells); + + SparseDirectUMFPACK inv; + inv.initialize(matrix); + + inv.vmult(solution, rhs); + constraints.distribute(solution); + + deallog << "Linfty norm of solution: " << solution.linfty_norm() << std::endl; +} + + +int +main() +{ + initlog(); + + test<2, 2>(); +} diff --git a/tests/meshworker/scratch_data_02_hp.with_muparser=on.with_umfpack=on.output b/tests/meshworker/scratch_data_02_hp.with_muparser=on.with_umfpack=on.output new file mode 100644 index 0000000000..b5ad838ee8 --- /dev/null +++ b/tests/meshworker/scratch_data_02_hp.with_muparser=on.with_umfpack=on.output @@ -0,0 +1,2 @@ + +DEAL::Linfty norm of solution: 0.0737281 diff --git a/tests/meshworker/scratch_data_03_hp.cc b/tests/meshworker/scratch_data_03_hp.cc new file mode 100644 index 0000000000..8f9807fc74 --- /dev/null +++ b/tests/meshworker/scratch_data_03_hp.cc @@ -0,0 +1,268 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2018 - 2021 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. + * + * --------------------------------------------------------------------- + */ + +// Solve Laplacian using SIPG + mesh_loop + ScratchData +// - hp variant + +#include +#include +#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" + +using namespace MeshWorker; + +template +void +test() +{ + Triangulation tria; + hp::FECollection fe(FE_DGQ(1), + FE_DGQ(2)); + DoFHandler dh(tria); + + FunctionParser rhs_function("1"); + FunctionParser boundary_function("0"); + + AffineConstraints constraints; + constraints.close(); + + GridGenerator::hyper_cube(tria); + tria.refine_global(5); + tria.execute_coarsening_and_refinement(); + + const auto bounding_box = GridTools::compute_bounding_box(tria); + for (auto &cell : dh.active_cell_iterators()) + { + bool use_alt_index = true; + for (unsigned int d = 0; d < spacedim; ++d) + if (cell->center()[d] < bounding_box.center()[d]) + { + use_alt_index = false; + break; + } + + if (use_alt_index) + cell->set_active_fe_index(1); + else + cell->set_active_fe_index(0); + } + + dh.distribute_dofs(fe); + + SparsityPattern sparsity; + + { + DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dh, dsp); + sparsity.copy_from(dsp); + } + + SparseMatrix matrix; + matrix.reinit(sparsity); + + Vector solution(dh.n_dofs()); + Vector rhs(dh.n_dofs()); + + hp::QCollection cell_quad(QGauss(3), QGauss(4)); + hp::QCollection face_quad(QGauss(3), QGauss(4)); + + UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + UpdateFlags face_flags = update_values | update_gradients | + update_quadrature_points | update_normal_vectors | + update_JxW_values; + + // Stabilization for SIPG + double gamma = 100; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<1 + GeometryInfo::faces_per_cell, + 1, + 1 + GeometryInfo::faces_per_cell>; + + ScratchData scratch(fe, cell_quad, cell_flags, face_quad, face_flags); + CopyData copy(0); + + auto cell = dh.begin_active(); + auto endc = dh.end(); + + using Iterator = decltype(cell); + + auto cell_worker = + [&rhs_function](const Iterator &cell, ScratchData &s, CopyData &c) { + const auto &fev = s.reinit(cell); + const auto &JxW = s.get_JxW_values(); + const auto &p = s.get_quadrature_points(); + + c.reinit(0, s.n_dofs_per_cell()); + c.local_dof_indices[0] = s.get_local_dof_indices(); + + for (unsigned int q = 0; q < p.size(); ++q) + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + { + c.matrices[0](i, j) += + fev.shape_grad(i, q) * fev.shape_grad(j, q) * JxW[q]; + } + c.vectors[0](i) += + fev.shape_value(i, q) * rhs_function.value(p[q]) * JxW[q]; + } + }; + + auto boundary_worker = [gamma, &boundary_function](const Iterator & cell, + const unsigned int &f, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f); + const auto &JxW = s.get_JxW_values(); + const auto &p = s.get_quadrature_points(); + const auto &n = s.get_normal_vectors(); + + for (unsigned int q = 0; q < p.size(); ++q) + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + { + c.matrices[0](i, j) += + (-fev.shape_grad(i, q) * n[q] * fev.shape_value(j, q) + + -fev.shape_grad(j, q) * n[q] * fev.shape_value(i, q) + + gamma / cell->face(f)->diameter() * fev.shape_value(i, q) * + fev.shape_value(j, q)) * + JxW[q]; + } + c.vectors[0](i) += + ((gamma / cell->face(f)->diameter() * fev.shape_value(i, q) - + fev.shape_grad(i, q) * n[q]) * + boundary_function.value(p[q])) * + JxW[q]; + } + }; + + auto face_worker = [gamma](const Iterator & cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator & ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, ncell, f, sf); + const auto &JxW = s.get_JxW_values(); + const auto &nfev = s.reinit_neighbor(cell, ncell, nf, nsf); + + c.reinit(f + 1, s.n_dofs_per_cell(), s.n_neighbor_dofs_per_cell()); + c.local_dof_indices[f + 1] = s.get_neighbor_dof_indices(); + + const auto &p = s.get_quadrature_points(); + const auto &n = s.get_normal_vectors(); + const auto &nn = s.get_neighbor_normal_vectors(); + + const double gh = gamma / cell->face(f)->diameter(); + + for (unsigned int q = 0; q < p.size(); ++q) + { + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + c.matrices[0](i, j) += + (-.5 * fev.shape_grad(i, q) * n[q] * fev.shape_value(j, q) + + -.5 * fev.shape_value(i, q) * n[q] * fev.shape_grad(j, q) + + gh * fev.shape_value(i, q) * fev.shape_value(j, q)) * + JxW[q]; + + for (unsigned int j = 0; j < nfev.dofs_per_cell; ++j) + c.matrices[f + 1](i, j) += + (-.5 * fev.shape_grad(i, q) * nn[q] * nfev.shape_value(j, q) + + -.5 * fev.shape_value(i, q) * n[q] * nfev.shape_grad(j, q) - + gh * fev.shape_value(i, q) * nfev.shape_value(j, q)) * + JxW[q]; + } + } + }; + + auto copier = [&constraints, &matrix, &rhs](const CopyData &c) { + constraints.distribute_local_to_global( + c.matrices[0], c.vectors[0], c.local_dof_indices[0], matrix, rhs); + + // We don't necessarily visit all faces, so we have to filter + // some of them out here. + for (const unsigned int f : GeometryInfo::face_indices()) + if (c.local_dof_indices[f + 1].size() > 0) + constraints.distribute_local_to_global(c.matrices[f + 1], + c.local_dof_indices[0], + c.local_dof_indices[f + 1], + matrix); + }; + + mesh_loop(cell, + endc, + cell_worker, + copier, + scratch, + copy, + assemble_own_cells | assemble_boundary_faces | + assemble_own_interior_faces_both, + boundary_worker, + face_worker); + + SparseDirectUMFPACK inv; + inv.initialize(matrix); + + inv.vmult(solution, rhs); + constraints.distribute(solution); + + deallog << "Linfty norm of solution " << solution.linfty_norm() << std::endl; +} + + +int +main() +{ + initlog(); + test<2, 2>(); +} diff --git a/tests/meshworker/scratch_data_03_hp.with_muparser=on.with_umfpack=on.output b/tests/meshworker/scratch_data_03_hp.with_muparser=on.with_umfpack=on.output new file mode 100644 index 0000000000..1b2777af16 --- /dev/null +++ b/tests/meshworker/scratch_data_03_hp.with_muparser=on.with_umfpack=on.output @@ -0,0 +1,2 @@ + +DEAL::Linfty norm of solution 0.0737281 -- 2.39.5