From 3c50b4b63efd99ae0c029165f19fd4ffe90043d3 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 14 Mar 2020 22:50:07 +0100 Subject: [PATCH] One new test, one included in existing one --- tests/matrix_free/compress_mapping.cc | 105 +++++++++++++++++ tests/matrix_free/compress_mapping.output | 8 ++ tests/matrix_free/mapping_info_01.cc | 134 ++++++++++++++++++++++ tests/matrix_free/mapping_info_01.output | 42 +++++++ 4 files changed, 289 insertions(+) create mode 100644 tests/matrix_free/mapping_info_01.cc create mode 100644 tests/matrix_free/mapping_info_01.output diff --git a/tests/matrix_free/compress_mapping.cc b/tests/matrix_free/compress_mapping.cc index cc9cf159cf..9fb09bf2b2 100644 --- a/tests/matrix_free/compress_mapping.cc +++ b/tests/matrix_free/compress_mapping.cc @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -176,6 +177,108 @@ test_parallelogram() +template +class DeformedManifold : public ChartManifold +{ +public: + DeformedManifold() = default; + + virtual std::unique_ptr> + clone() const + { + return std_cxx14::make_unique>(); + } + + virtual Point + push_forward(const Point &chart_point) const + { + Point result = chart_point; + result[0] = std::tanh(2.0 * result[0]) / std::tanh(2.0); + return result; + } + + virtual Point + pull_back(const Point &space_point) const + { + Point result = space_point; + result[0] = 0.5 * std::atanh(result[0] * std::tanh(2.0)); + return result; + } +}; + + + +template +void +test_deformed_cube() +{ + deallog << "Deformed hyper cube" << std::endl; + Triangulation tria; + GridGenerator::hyper_cube(tria, -1., 1.); + tria.set_all_manifold_ids(0); + tria.set_manifold(0, DeformedManifold()); + tria.refine_global(6 - dim); + + FE_Q fe(2); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + constraints.close(); + + const QGauss<1> quad(3); + MatrixFree mf; + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = + update_gradients | update_normal_vectors; + data.mapping_update_flags_boundary_faces = + update_gradients | update_normal_vectors; + + mf.reinit(dof, constraints, quad, data); + const unsigned int n_macro_cells = mf.n_macro_cells(); + Assert(n_macro_cells > 1, ExcInternalError()); + + { + std::vector n_cell_types(4, 0); + for (unsigned int i = 0; i < n_macro_cells; ++i) + n_cell_types[mf.get_mapping_info().get_cell_type(i)]++; + + // should have all Cartesian type and no other cell type + AssertDimension(n_cell_types[0], n_macro_cells); + + // should have as many different Jacobians as we have cell batches in x + // direction; as the mesh is a cube, we can easily calculate it + deallog << "Number of different Jacobians: " + << mf.get_mapping_info().cell_data[0].jacobians[0].size() / 2 + << std::endl; + } + + // check again, now using a mapping that displaces points + { + MappingQGeneric mapping(3); + mf.reinit(mapping, dof, constraints, quad, data); + + std::vector n_cell_types(4, 0); + for (unsigned int i = 0; i < n_macro_cells; ++i) + n_cell_types[mf.get_mapping_info().get_cell_type(i)]++; + + // should have all general type and no other cell type + AssertDimension(n_cell_types[3], n_macro_cells); + + // should have as many different Jacobians as we have cell batches in x + // direction times the number of quadrature points; as the mesh is a cube, + // we can easily calculate it + deallog << "Number of different Jacobians: " + << mf.get_mapping_info().cell_data[0].jacobians[0].size() + << std::endl; + } + + deallog << "OK" << std::endl; +} + + + int main() { @@ -188,11 +291,13 @@ main() test<2>(); test_cube<2>(); test_parallelogram<2>(); + test_deformed_cube<2>(); deallog.pop(); deallog.push("3d"); test<3>(); test_cube<3>(); test_parallelogram<3>(); + test_deformed_cube<3>(); deallog.pop(); } } diff --git a/tests/matrix_free/compress_mapping.output b/tests/matrix_free/compress_mapping.output index 9ced1ac068..ec4b7b3083 100644 --- a/tests/matrix_free/compress_mapping.output +++ b/tests/matrix_free/compress_mapping.output @@ -5,9 +5,17 @@ DEAL:2d::Hyper cube DEAL:2d::OK DEAL:2d::Parallelogram DEAL:2d::OK +DEAL:2d::Deformed hyper cube +DEAL:2d::Number of different Jacobians: 8 +DEAL:2d::Number of different Jacobians: 72 +DEAL:2d::OK DEAL:3d::General mesh DEAL:3d::OK DEAL:3d::Hyper cube DEAL:3d::OK DEAL:3d::Parallelogram DEAL:3d::OK +DEAL:3d::Deformed hyper cube +DEAL:3d::Number of different Jacobians: 4 +DEAL:3d::Number of different Jacobians: 108 +DEAL:3d::OK diff --git a/tests/matrix_free/mapping_info_01.cc b/tests/matrix_free/mapping_info_01.cc new file mode 100644 index 0000000000..1afa066ca7 --- /dev/null +++ b/tests/matrix_free/mapping_info_01.cc @@ -0,0 +1,134 @@ +// --------------------------------------------------------------------- +// +// 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Let MatrixFree initialize the geometry for ridiculously high numbers of +// quadrature points and compute the volume and surface of a hyper ball. In an +// initial implementation we used to consume a lot of memory and take a long +// time, whereas the new state is relatively cheap. + +#include + +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +template +void +test(const unsigned int degree) +{ + Triangulation tria; + GridGenerator::hyper_ball(tria, Point(), 1.); + + FE_DGQ fe(1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + MappingQGeneric mapping(degree); + + MatrixFree mf_data; + const QGauss<1> quad(degree); + typename MatrixFree::AdditionalData data; + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + mf_data.reinit( + mapping, + std::vector *>({{&dof}}), + std::vector *>({{&constraints}}), + std::vector>({{QGauss<1>(std::max(degree / 2, 1U)), + QGauss<1>(degree + 1), + QGauss<1>(3 * degree / 2)}}), + data); + + deallog << std::setw(5) << degree; + for (unsigned int index = 0; index < 3; ++index) + { + double volume = 0; + FEEvaluation fe_eval(mf_data, 0, index); + for (unsigned int cell = 0; cell < mf_data.n_cell_batches(); ++cell) + { + fe_eval.reinit(cell); + VectorizedArray local = 0; + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + local += fe_eval.JxW(q); + for (unsigned int v = 0; + v < mf_data.n_active_entries_per_cell_batch(cell); + ++v) + volume += local[v]; + } + deallog << std::setw(11) + << std::abs((dim == 2 ? numbers::PI : 4. / 3. * numbers::PI) - + volume); + } + deallog << " "; + for (unsigned int index = 0; index < 3; ++index) + { + double area = 0; + FEFaceEvaluation fe_eval(mf_data, true, 0, index); + for (unsigned int face = mf_data.n_inner_face_batches(); + face < + mf_data.n_boundary_face_batches() + mf_data.n_inner_face_batches(); + ++face) + { + fe_eval.reinit(face); + VectorizedArray local = 0; + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + local += fe_eval.JxW(q); + for (unsigned int v = 0; + v < mf_data.n_active_entries_per_face_batch(face); + ++v) + area += local[v]; + } + deallog << std::setw(11) << std::abs(2. * (dim - 1) * numbers::PI - area); + } + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + + deallog << std::setprecision(4); + deallog << "2D volume area" << std::endl; + deallog << "deg err deg/2 err deg+1 err 3deg/2 " + << "err deg/2 err deg+1 err 3deg/2" << std::endl; + for (unsigned int q = 1; q < 20; ++q) + test<2>(q); + deallog << std::endl; + + deallog << "3D volume area" << std::endl; + deallog << "deg err deg/2 err deg+1 err 3deg/2 " + << "err deg/2 err deg+1 err 3deg/2" << std::endl; + for (unsigned int q = 1; q < 17; ++q) + test<3>(q); + deallog << std::endl; + + return 0; +} diff --git a/tests/matrix_free/mapping_info_01.output b/tests/matrix_free/mapping_info_01.output new file mode 100644 index 0000000000..885e6fcddc --- /dev/null +++ b/tests/matrix_free/mapping_info_01.output @@ -0,0 +1,42 @@ + +DEAL::2D volume area +DEAL::deg err deg/2 err deg+1 err 3deg/2 err deg/2 err deg+1 err 3deg/2 +DEAL::1 1.142 1.142 1.142 0.6263 0.6263 0.6263 +DEAL::2 0.02998 0.03702 0.03702 0.6263 0.03493 0.03493 +DEAL::3 0.2367 0.0005893 0.0005893 0.003915 0.0001972 0.0001972 +DEAL::4 0.004112 5.461e-06 5.461e-06 0.008015 5.006e-06 4.602e-06 +DEAL::5 0.006630 3.288e-08 3.288e-08 8.248e-05 8.603e-09 1.122e-08 +DEAL::6 3.905e-05 1.386e-10 1.386e-10 4.224e-05 1.254e-10 1.136e-10 +DEAL::7 5.266e-05 4.334e-13 4.370e-13 3.522e-07 9.948e-14 1.368e-13 +DEAL::8 1.580e-07 8.882e-16 3.553e-15 1.179e-07 1.776e-15 2.665e-15 +DEAL::9 1.964e-07 6.661e-15 1.776e-15 8.182e-10 3.553e-15 2.665e-15 +DEAL::10 3.578e-10 3.553e-15 2.665e-15 2.036e-10 0.000 8.882e-16 +DEAL::11 4.246e-10 2.665e-15 3.553e-15 1.210e-12 1.776e-15 1.776e-15 +DEAL::12 5.302e-13 8.438e-15 1.332e-15 2.327e-13 8.882e-16 1.776e-15 +DEAL::13 5.844e-13 4.441e-15 1.776e-15 7.105e-15 1.776e-15 8.882e-16 +DEAL::14 8.882e-16 1.776e-15 1.155e-14 8.882e-16 1.776e-15 8.882e-16 +DEAL::15 7.105e-15 5.329e-15 7.994e-15 6.217e-15 8.882e-15 4.441e-15 +DEAL::16 1.865e-14 1.243e-14 4.885e-15 8.882e-16 0.000 1.776e-15 +DEAL::17 1.243e-14 1.332e-15 1.332e-14 8.882e-16 8.882e-16 1.776e-15 +DEAL::18 4.441e-16 4.885e-15 9.326e-15 2.665e-15 1.776e-15 8.882e-16 +DEAL::19 3.997e-15 5.773e-15 7.994e-15 4.441e-15 0.000 2.665e-15 +DEAL:: +DEAL::3D volume area +DEAL::deg err deg/2 err deg+1 err 3deg/2 err deg/2 err deg+1 err 3deg/2 +DEAL::1 2.747 2.649 2.747 4.566 4.566 4.566 +DEAL::2 0.1217 0.1120 0.1120 0.5664 0.2217 0.2217 +DEAL::3 0.8702 0.01703 0.01703 3.949 0.02024 0.02024 +DEAL::4 0.03972 0.0005453 0.0005447 0.03028 0.001079 0.001038 +DEAL::5 0.1274 0.0001358 0.0001358 0.3096 5.472e-05 6.758e-05 +DEAL::6 0.002403 5.127e-06 5.121e-06 0.006959 1.007e-05 9.541e-06 +DEAL::7 0.01271 1.081e-06 1.080e-06 0.03144 1.443e-07 2.279e-08 +DEAL::8 0.0001101 4.117e-08 4.111e-08 0.001177 8.077e-08 7.626e-08 +DEAL::9 0.001234 9.165e-09 9.161e-09 0.002962 7.696e-09 5.659e-09 +DEAL::10 7.936e-07 3.567e-10 3.561e-10 0.0001491 6.945e-10 6.492e-10 +DEAL::11 0.0001194 7.771e-11 7.768e-11 0.0002851 1.123e-10 9.018e-11 +DEAL::12 7.006e-07 3.033e-12 3.001e-12 1.740e-05 5.846e-12 5.437e-12 +DEAL::13 1.151e-05 6.777e-13 6.688e-13 2.715e-05 1.450e-12 1.204e-12 +DEAL::14 1.345e-07 3.553e-14 9.770e-15 1.934e-06 4.619e-14 3.908e-14 +DEAL::15 1.112e-06 1.688e-14 4.885e-14 2.602e-06 2.842e-14 3.197e-14 +DEAL::16 1.865e-08 3.464e-14 1.332e-14 2.088e-07 1.421e-14 1.066e-14 +DEAL:: -- 2.39.5