From c2553dd9c2d62e5fc5130ba8d57ab300e1240842 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 22 Jul 2019 18:31:02 +0200 Subject: [PATCH] Add tests for the new level feature --- tests/mappings/mapping_fe_field_04.cc | 112 ++++++++++++ tests/mappings/mapping_fe_field_04.output | 7 + tests/mappings/mapping_fe_field_05.cc | 114 ++++++++++++ tests/mappings/mapping_fe_field_05.output | 7 + tests/mappings/mapping_fe_field_06.cc | 172 ++++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=3.output | 23 +++ 6 files changed, 435 insertions(+) create mode 100644 tests/mappings/mapping_fe_field_04.cc create mode 100644 tests/mappings/mapping_fe_field_04.output create mode 100644 tests/mappings/mapping_fe_field_05.cc create mode 100644 tests/mappings/mapping_fe_field_05.output create mode 100644 tests/mappings/mapping_fe_field_06.cc create mode 100644 tests/mappings/mapping_fe_field_06.with_mpi=true.with_p4est=true.mpirun=3.output diff --git a/tests/mappings/mapping_fe_field_04.cc b/tests/mappings/mapping_fe_field_04.cc new file mode 100644 index 0000000000..af2f1ec978 --- /dev/null +++ b/tests/mappings/mapping_fe_field_04.cc @@ -0,0 +1,112 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 MappingFEField when initialized on all multigrid levels + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + Triangulation tria( + Triangulation::limit_level_difference_at_vertices); + GridGenerator::hyper_cube(tria); + + tria.refine_global(2); + + FESystem fe(FE_Q(2), spacedim); + DoFHandler dh(tria); + + dh.distribute_dofs(fe); + dh.distribute_mg_dofs(); + + deallog << "dim, spacedim: " << dim << ", " << spacedim << std::endl + << "cells: " << tria.n_active_cells() << ", dofs: " << dh.n_dofs() + << std::endl; + + // Create a Mapping + LinearAlgebra::distributed::Vector map_vector(dh.n_dofs()); + VectorTools::get_position_vector(dh, map_vector); + MGLevelObject> level_vectors( + 0, tria.n_global_levels() - 1); + for (unsigned int level = 0; level < tria.n_global_levels(); ++level) + level_vectors[level].reinit(dh.n_dofs(level)); + + MGTransferMatrixFree transfer; + transfer.build(dh); + transfer.interpolate_to_mg(dh, level_vectors, map_vector); + MappingFEField, + DoFHandler> + mapping(dh, level_vectors); + + QGauss quad(1); + FEValues fe_values_ref(fe, + quad, + update_jacobians | update_quadrature_points); + FEValues fe_values(mapping, + fe, + quad, + update_jacobians | update_quadrature_points); + + for (const auto &cell : tria.cell_iterators()) + { + fe_values_ref.reinit(cell); + fe_values.reinit(cell); + + if (fe_values_ref.quadrature_point(0).distance( + fe_values.quadrature_point(0)) > 1e-12) + deallog << "Mapped point should be " + << fe_values_ref.quadrature_point(0) << " and is " + << fe_values.quadrature_point(0) << std::endl; + Tensor<2, dim> jac_ref = fe_values_ref.jacobian(0), + jac = fe_values.jacobian(0); + if ((jac_ref - jac).norm() > 1e-12) + deallog << "Jacobian should be " << jac_ref << " and is " << jac + << std::endl; + } + deallog << "OK" << std::endl; +} + +using namespace dealii; +int +main() +{ + initlog(); + test<2, 2>(); + test<3, 3>(); +} diff --git a/tests/mappings/mapping_fe_field_04.output b/tests/mappings/mapping_fe_field_04.output new file mode 100644 index 0000000000..a9e8f81573 --- /dev/null +++ b/tests/mappings/mapping_fe_field_04.output @@ -0,0 +1,7 @@ + +DEAL::dim, spacedim: 2, 2 +DEAL::cells: 16, dofs: 162 +DEAL::OK +DEAL::dim, spacedim: 3, 3 +DEAL::cells: 64, dofs: 2187 +DEAL::OK diff --git a/tests/mappings/mapping_fe_field_05.cc b/tests/mappings/mapping_fe_field_05.cc new file mode 100644 index 0000000000..0a9cf0026b --- /dev/null +++ b/tests/mappings/mapping_fe_field_05.cc @@ -0,0 +1,114 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 MappingFEField when initialized on all multigrid levels + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + Triangulation tria( + Triangulation::limit_level_difference_at_vertices); + GridGenerator::hyper_ball(tria); + + tria.refine_global(1); + + FESystem fe(FE_Q(2), spacedim); + DoFHandler dh(tria); + + dh.distribute_dofs(fe); + dh.distribute_mg_dofs(); + + deallog << "dim, spacedim: " << dim << ", " << spacedim << std::endl + << "cells: " << tria.n_active_cells() << ", dofs: " << dh.n_dofs() + << std::endl; + + // Create a Mapping + LinearAlgebra::distributed::Vector map_vector(dh.n_dofs()); + VectorTools::get_position_vector(dh, map_vector); + MGLevelObject> level_vectors( + 0, tria.n_global_levels() - 1); + for (unsigned int level = 0; level < tria.n_global_levels(); ++level) + level_vectors[level].reinit(dh.n_dofs(level)); + + MGTransferMatrixFree transfer; + transfer.build(dh); + transfer.interpolate_to_mg(dh, level_vectors, map_vector); + MappingFEField, + DoFHandler> + mapping(dh, level_vectors); + MappingQGeneric mapping_ref(fe.degree); + + QGauss quad(1); + FEValues fe_values_ref(mapping_ref, + fe, + quad, + update_jacobians | update_quadrature_points); + FEValues fe_values(mapping, + fe, + quad, + update_jacobians | update_quadrature_points); + + for (const auto &cell : tria.cell_iterators()) + { + fe_values_ref.reinit(cell); + fe_values.reinit(cell); + + if (fe_values_ref.quadrature_point(0).distance( + fe_values.quadrature_point(0)) > 1e-12) + deallog << "Mapped point should be " + << fe_values_ref.quadrature_point(0) << " and is " + << fe_values.quadrature_point(0) << std::endl; + Tensor<2, dim> jac_ref = fe_values_ref.jacobian(0), + jac = fe_values.jacobian(0); + if ((jac_ref - jac).norm() > 1e-12) + deallog << "Jacobian should be " << jac_ref << " and is " << jac + << std::endl; + } + deallog << "OK" << std::endl; +} + +using namespace dealii; +int +main() +{ + initlog(); + test<2, 2>(); + test<3, 3>(); +} diff --git a/tests/mappings/mapping_fe_field_05.output b/tests/mappings/mapping_fe_field_05.output new file mode 100644 index 0000000000..eca0ce18d5 --- /dev/null +++ b/tests/mappings/mapping_fe_field_05.output @@ -0,0 +1,7 @@ + +DEAL::dim, spacedim: 2, 2 +DEAL::cells: 20, dofs: 178 +DEAL::OK +DEAL::dim, spacedim: 3, 3 +DEAL::cells: 56, dofs: 1551 +DEAL::OK diff --git a/tests/mappings/mapping_fe_field_06.cc b/tests/mappings/mapping_fe_field_06.cc new file mode 100644 index 0000000000..4a7ec1698e --- /dev/null +++ b/tests/mappings/mapping_fe_field_06.cc @@ -0,0 +1,172 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 MappingFEField when initialized manually on levels on a distributed +// triangulation + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + parallel::distributed::Triangulation tria( + MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation:: + construct_multigrid_hierarchy); + GridGenerator::hyper_ball(tria); + + tria.refine_global(2); + if (tria.begin_active()->is_locally_owned()) + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + FESystem fe(FE_Q(2), spacedim); + DoFHandler dh(tria); + + dh.distribute_dofs(fe); + dh.distribute_mg_dofs(); + + deallog << "dim, spacedim: " << dim << ", " << spacedim << std::endl + << "cells: " << tria.n_active_cells() << ", dofs: " << dh.n_dofs() + << std::endl; + + // Create a Mapping + std::vector> level_vectors( + tria.n_global_levels()); + MappingQGeneric mapping_ref(fe.degree); + FEValues fe_values_setup(mapping_ref, + dh.get_fe(), + Quadrature( + dh.get_fe().get_unit_support_points()), + update_quadrature_points); + for (unsigned int level = 0; level < tria.n_global_levels(); ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dh, level, relevant_dofs); + level_vectors[level].reinit(dh.locally_owned_mg_dofs(level), + relevant_dofs, + tria.get_communicator()); + std::vector dof_indices(fe.dofs_per_cell); + for (const auto &cell : dh.mg_cell_iterators_on_level(level)) + if (cell->level_subdomain_id() != numbers::artificial_subdomain_id) + { + fe_values_setup.reinit(cell); + cell->get_active_or_mg_dof_indices(dof_indices); + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + { + const unsigned int coordinate_direction = + fe.system_to_component_index(i).first; + const Point point = fe_values_setup.quadrature_point(i); + level_vectors[level](dof_indices[i]) = + point[coordinate_direction]; + } + } + level_vectors[level].update_ghost_values(); + } + + MappingFEField, + DoFHandler> + mapping(dh, level_vectors); + + QGauss quad(1); + FEValues fe_values_ref(mapping_ref, + fe, + quad, + update_jacobians | update_quadrature_points); + FEValues fe_values(mapping, + fe, + quad, + update_jacobians | update_quadrature_points); + + for (const auto &cell : tria.cell_iterators()) + if (cell->level_subdomain_id() != numbers::artificial_subdomain_id) + { + fe_values_ref.reinit(cell); + fe_values.reinit(cell); + + if (fe_values_ref.quadrature_point(0).distance( + fe_values.quadrature_point(0)) > 1e-12) + deallog << "Mapped point should be " + << fe_values_ref.quadrature_point(0) << " and is " + << fe_values.quadrature_point(0) << std::endl; + Tensor<2, dim> jac_ref = fe_values_ref.jacobian(0), + jac = fe_values.jacobian(0); + if ((jac_ref - jac).norm() > 1e-12) + deallog << "Jacobian should be " << jac_ref << " and is " << jac + << std::endl; + } + + // shift the Euler vectors and check whether the result is still correct + for (unsigned int level = 0; level < tria.n_global_levels(); ++level) + level_vectors[level].add(1.1); + + for (const auto &cell : tria.cell_iterators()) + if (cell->level_subdomain_id() != numbers::artificial_subdomain_id) + { + fe_values_ref.reinit(cell); + fe_values.reinit(cell); + + Point shift; + for (unsigned int d = 0; d < dim; ++d) + shift[d] = 1.1; + if ((fe_values_ref.quadrature_point(0) + shift) + .distance(fe_values.quadrature_point(0)) > 1e-12) + deallog << "Mapped point should be " + << fe_values_ref.quadrature_point(0) + shift << " and is " + << fe_values.quadrature_point(0) << std::endl; + Tensor<2, dim> jac_ref = fe_values_ref.jacobian(0), + jac = fe_values.jacobian(0); + if ((jac_ref - jac).norm() > 1e-12) + deallog << "Jacobian should be " << jac_ref << " and is " << jac + << std::endl; + } + deallog << "OK" << std::endl; +} + +using namespace dealii; +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, + argv, + testing_max_num_threads()); + MPILogInitAll log; + test<2, 2>(); + test<3, 3>(); +} diff --git a/tests/mappings/mapping_fe_field_06.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/mappings/mapping_fe_field_06.with_mpi=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..6048744c3e --- /dev/null +++ b/tests/mappings/mapping_fe_field_06.with_mpi=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,23 @@ + +DEAL:0::dim, spacedim: 2, 2 +DEAL:0::cells: 65, dofs: 674 +DEAL:0::OK +DEAL:0::dim, spacedim: 3, 3 +DEAL:0::cells: 392, dofs: 11451 +DEAL:0::OK + +DEAL:1::dim, spacedim: 2, 2 +DEAL:1::cells: 74, dofs: 674 +DEAL:1::OK +DEAL:1::dim, spacedim: 3, 3 +DEAL:1::cells: 448, dofs: 11451 +DEAL:1::OK + + +DEAL:2::dim, spacedim: 2, 2 +DEAL:2::cells: 65, dofs: 674 +DEAL:2::OK +DEAL:2::dim, spacedim: 3, 3 +DEAL:2::cells: 420, dofs: 11451 +DEAL:2::OK + -- 2.39.5