From: Jean-Paul Pelteret Date: Thu, 30 Apr 2020 21:58:42 +0000 (+0200) Subject: Add another test for mesh_loop X-Git-Tag: v9.2.0-rc1~123^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d734aa0dc43f0a889a25e52ddfeea1304f39b1c4;p=dealii.git Add another test for mesh_loop --- diff --git a/tests/meshworker/mesh_loop_07.cc b/tests/meshworker/mesh_loop_07.cc new file mode 100644 index 0000000000..0b3b71cf85 --- /dev/null +++ b/tests/meshworker/mesh_loop_07.cc @@ -0,0 +1,169 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Check that mesh loop returns the correct values for integral over +// volumes and boundaries + +#include + +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + + +template +void +run() +{ + deallog << "Dim: " << dim << std::endl; + + const FE_Q fe(1); + const QGauss cell_quadrature(fe.degree + 1); + const QGauss face_quadrature(fe.degree + 1); + + Triangulation triangulation; + GridGenerator::subdivided_hyper_cube(triangulation, 4, 0.0, 1.0); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<1, 1, 1>; + using CellIteratorType = decltype(dof_handler.begin_active()); + + // Volume integral + { + ScratchData scratch(dof_handler.get_fe(), + cell_quadrature, + update_JxW_values); + CopyData copy(1); + + auto cell_worker = [](const CellIteratorType &cell, + ScratchData & scratch_data, + CopyData & copy_data) { + const auto &fe_values = scratch_data.reinit(cell); + copy_data = 0; + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + { + copy_data.vectors[0][0] += 1.0 * fe_values.JxW(q_point); + } + }; + + double volume = 0.0; + auto copier = [&volume](const CopyData ©_data) { + volume += copy_data.vectors[0][0]; + }; + + MeshWorker::mesh_loop(dof_handler.active_cell_iterators(), + cell_worker, + copier, + scratch, + copy, + MeshWorker::assemble_own_cells); + deallog << "Volume: " << volume << std::endl; + + double reference_volume = 0.0; + for (auto &cell : dof_handler.active_cell_iterators()) + reference_volume += cell->measure(); + + Assert(std::abs(volume - reference_volume) < 1e-6, + ExcMessage("Volumes do not match. Reference value: " + + Utilities::to_string(reference_volume) + + "; Calculated value: " + Utilities::to_string(volume))); + } + + // Boundary integral + { + ScratchData scratch(dof_handler.get_fe(), + cell_quadrature, + update_default, + face_quadrature, + update_JxW_values); + CopyData copy(1); + + std::function + empty_cell_worker; + + auto boundary_worker = [](const CellIteratorType &cell, + const unsigned int face, + ScratchData & scratch_data, + CopyData & copy_data) { + const auto &fe_face_values = scratch_data.reinit(cell, face); + + for (unsigned int q_point = 0; + q_point < fe_face_values.n_quadrature_points; + ++q_point) + { + copy_data.vectors[0][0] += 1.0 * fe_face_values.JxW(q_point); + } + }; + + double area = 0.0; + auto copier = [&area](const CopyData ©_data) { + // *** + // The problem is observed here. For some faces, the copy_data that + // is passed in is reset, even though there is always work being + // done in the boundary_worker. + // *** + area += copy_data.vectors[0][0]; + }; + + MeshWorker::mesh_loop(dof_handler.active_cell_iterators(), + empty_cell_worker, + copier, + scratch, + copy, + MeshWorker::assemble_boundary_faces, + boundary_worker); + deallog << "Area: " << area << std::endl; + + double reference_area = 0.0; + for (auto &cell : dof_handler.active_cell_iterators()) + for (const unsigned int face : GeometryInfo::face_indices()) + if (cell->face(face)->at_boundary()) + { + reference_area += cell->face(face)->measure(); + } + + Assert(std::abs(area - reference_area) < 1e-6, + ExcMessage("Areas do not match. Reference value: " + + Utilities::to_string(reference_area) + + "; Calculated value: " + Utilities::to_string(area))); + } + + deallog << "OK" << std::endl; +} + + +int +main() +{ + initlog(); + + run<2>(); + run<3>(); + + deallog << "OK" << std::endl; +} diff --git a/tests/meshworker/mesh_loop_07.output b/tests/meshworker/mesh_loop_07.output new file mode 100644 index 0000000000..fccc9bdb0c --- /dev/null +++ b/tests/meshworker/mesh_loop_07.output @@ -0,0 +1,10 @@ + +DEAL::Dim: 2 +DEAL::Volume: 1.00000 +DEAL::Area: 4.00000 +DEAL::OK +DEAL::Dim: 3 +DEAL::Volume: 1.00000 +DEAL::Area: 6.00000 +DEAL::OK +DEAL::OK