From: Marc Fehling Date: Mon, 1 Feb 2021 01:15:18 +0000 (-0700) Subject: simplex: verify hanging nodes on triangular and hybrid meshes. X-Git-Tag: v9.3.0-rc1~466^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=74b8171dd030a4b2de41955b8d913cde06974b91;p=dealii.git simplex: verify hanging nodes on triangular and hybrid meshes. --- diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 31bf8691d2..843c490b52 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -552,7 +552,7 @@ namespace GridTools for (const auto &cell : tria.cell_iterators_on_level(0)) { // Save cell data - CellData cell_data; + CellData cell_data(cell->n_vertices()); for (const unsigned int cell_vertex_n : cell->vertex_indices()) { Assert(cell->vertex_index(cell_vertex_n) < vertices.size(), diff --git a/tests/simplex/hanging_nodes.h b/tests/simplex/hanging_nodes.h new file mode 100644 index 0000000000..309ac2eb90 --- /dev/null +++ b/tests/simplex/hanging_nodes.h @@ -0,0 +1,238 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + + +// miscellaneous helper funcions for hanging_nodes_* tests + + +#include +#include + +#include + +#include +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + + +// ----- setup ----- + +/** + * On a Triangulation @p tria, refine those cells which are descendants of the + * same coarse cell by a specified amount of times. + * + * The index of an entry in @p n_refinements corresponds to the globally unique + * coarse cell id, while the entry itself describes the number of refinements. + */ +template +void +refine(const std::vector &n_refinements, Triangulation &tria) +{ + AssertDimension(n_refinements.size(), tria.n_cells(0)); + + const unsigned int max_level = + *std::max_element(n_refinements.begin(), n_refinements.end()); + + while (tria.n_levels() <= max_level) + { + for (const auto &cell : tria.active_cell_iterators()) + { + const unsigned int coarse_cell_id = cell->id().get_coarse_cell_id(); + + if (static_cast(cell->level()) < + n_refinements[coarse_cell_id]) + cell->set_refine_flag(); + } + + tria.execute_coarsening_and_refinement(); + } + +#ifdef DEBUG + for (const auto &cell : tria.active_cell_iterators()) + { + const unsigned int coarse_cell_id = cell->id().get_coarse_cell_id(); + + AssertDimension(cell->level(), n_refinements[coarse_cell_id]); + } +#endif +} + + +/** + * On a DoFHandler @p dofh, assign the same active FE index to all cells which + * are descendants of the same coarse cell. + * + * The index of an entry in @p fe_indices corresponds to the globally unique + * coarse cell id, while the entry itself is the active FE index. + */ +template +void +set_active_fe_indices(const std::vector &fe_indices, + DoFHandler & dofh) +{ + AssertDimension(fe_indices.size(), dofh.get_triangulation().n_cells(0)); + + for (const auto &cell : dofh.active_cell_iterators()) + { + const unsigned int coarse_cell_id = cell->id().get_coarse_cell_id(); + + cell->set_active_fe_index(fe_indices[coarse_cell_id]); + } +} + + + +// ----- diagnostics ----- + +/** + * Print the indices of degrees of freedom for each face on each cell in + * deallog. + * + * If a face is divided into subfaces on a neighboring cell which is refined, + * then also print the dof indices on the subfaces from the perspective of the + * neighboring cell. + */ +template +void +print_dof_indices_on_faces(const DoFHandler &dofh) +{ + std::vector dof_indices; + + for (const auto &cell : dofh.active_cell_iterators()) + { + const auto & fe = cell->get_fe(); + const unsigned int fe_index = cell->active_fe_index(); + + for (unsigned int f = 0; f < cell->n_faces(); ++f) + { + const auto &face = cell->face(f); + Assert(face->fe_index_is_active(fe_index), ExcInternalError()); + + const unsigned int n_dofs = + internal::DoFAccessorImplementation::Implementation::n_dof_indices( + *face, fe_index); + dof_indices.resize(n_dofs); + face->get_dof_indices(dof_indices, fe_index); + + deallog << "cell:" << cell->active_cell_index() << " face:" << f + << " dofs:"; + for (const auto &i : dof_indices) + deallog << i << " "; + deallog << std::endl; + + // also print dofs on subfaces if neighboring cell is refined. + // in this case, only one fe should be active on the subface. + for (unsigned int sf = 0; sf < face->n_children(); ++sf) + { + const auto & subface = face->child(sf); + const unsigned int subface_fe_index = + subface->nth_active_fe_index(0); + Assert(subface->n_active_fe_indices() == 1, ExcInternalError()); + + const unsigned int n_dofs = internal::DoFAccessorImplementation:: + Implementation::n_dof_indices(*subface, subface_fe_index); + dof_indices.resize(n_dofs); + subface->get_dof_indices(dof_indices, subface_fe_index); + + deallog << " subface:" << sf << " dofs:"; + for (const auto &i : dof_indices) + deallog << i << " "; + deallog << std::endl; + } + } + } +} + + +/** + * Print the index and physical coordinate of each degree of freedom in deallog. + */ +template +void +print_dof_points(const DoFHandler &dofh) +{ + hp::MappingCollection mapping; + for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i) + mapping.push_back(MappingFE(dofh.get_fe(i))); + + std::vector> points(dofh.n_dofs()); + DoFTools::map_dofs_to_support_points(mapping, dofh, points); + + for (unsigned int i = 0; i < dofh.n_dofs(); ++i) + deallog << "dof:" << i << " point:" << points[i] << std::endl; +} + + + +// ----- tests ----- + +/** + * Verify hanging node constraints on locally refined meshes. + * + * Creates a Triangulation based on @p grid_generator, refines the coarse cells + * @p n_refinements times, and assigns @p fe_indices to all descendants of the + * coarse cells based on the provided @p fe_collection. + * + * Makes hanging node constraints and prints them in deallog. + */ +template +void +test(const std::vector & n_refinements, + const std::vector & fe_indices, + const hp::FECollection & fe_collection, + const std::function &)> &grid_generator) +{ + // setup grid + Triangulation tria; + grid_generator(tria); + + AssertDimension(n_refinements.size(), tria.n_cells()); + AssertDimension(fe_indices.size(), tria.n_cells()); + + refine(n_refinements, tria); + +#if false + GridOut grid_out; + grid_out.write_vtk(tria, deallog.get_file_stream()); +#endif + + DoFHandler dofh(tria); + set_active_fe_indices(fe_indices, dofh); + dofh.distribute_dofs(fe_collection); + deallog << "ndofs: " << dofh.n_dofs() << std::endl; + +#if false + print_dof_points(dofh); + print_dof_indices_on_faces(dofh); +#endif + + // hanging node constraints + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dofh, constraints); + constraints.close(); + constraints.print(deallog.get_file_stream()); + + deallog << "OK" << std::endl; +} diff --git a/tests/simplex/hanging_nodes_01.cc b/tests/simplex/hanging_nodes_01.cc index 20762147c6..ac96c5eb37 100644 --- a/tests/simplex/hanging_nodes_01.cc +++ b/tests/simplex/hanging_nodes_01.cc @@ -31,8 +31,7 @@ #include #include -#include - +#include #include #include @@ -43,108 +42,7 @@ #include "../tests.h" - -// ----- diagnostics ----- - -template -void -print_dof_indices_on_faces(const DoFHandler &dofh) -{ - std::vector dof_indices; - - for (const auto &cell : dofh.active_cell_iterators()) - for (unsigned int f = 0; f < cell->n_faces(); ++f) - { - const auto &face = cell->face(f); - - if (face->has_children()) - { - for (unsigned int sf = 0; sf < face->n_children(); ++sf) - { - const auto &subface = face->child(sf); - Assert(subface->n_active_fe_indices() == 1, ExcInternalError()); - const unsigned int subface_fe_index = - subface->nth_active_fe_index(0); - const auto &subface_fe = dofh.get_fe(subface_fe_index); - - dof_indices.resize(subface_fe.n_dofs_per_face(f)); - subface->get_dof_indices(dof_indices, subface_fe_index); - - deallog << "cell:" << cell->active_cell_index() << " face:" << f - << " subface:" << sf << " dofs:"; - for (const auto &i : dof_indices) - deallog << i << " "; - deallog << std::endl; - } - } - else - { - Assert(face->n_active_fe_indices() == 1, ExcInternalError()); - const unsigned int face_fe_index = face->nth_active_fe_index(0); - const auto & face_fe = dofh.get_fe(face_fe_index); - - dof_indices.resize(face_fe.n_dofs_per_face(f)); - face->get_dof_indices(dof_indices, face_fe_index); - - deallog << "cell:" << cell->active_cell_index() << " face:" << f - << " dofs:"; - for (const auto &i : dof_indices) - deallog << i << " "; - deallog << std::endl; - } - } -} - - -template -void -print_dof_points(const DoFHandler &dofh) -{ - std::vector> points(dofh.n_dofs()); - DoFTools::map_dofs_to_support_points(MappingFE(dofh.get_fe()), - dofh, - points); - - for (unsigned int i = 0; i < dofh.n_dofs(); ++i) - deallog << "dof:" << i << " point:" << points[i] << std::endl; -} - - - -// ----- test ----- - -template -void -test() -{ - // setup grid - Triangulation tria; - GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); - - tria.begin_active()->set_refine_flag(); - tria.execute_coarsening_and_refinement(); - -#if false - GridOut grid_out; - grid_out.write_vtk(tria, deallog.get_file_stream()); -#endif - - DoFHandler dofh(tria); - dofh.distribute_dofs(Simplex::FE_P(1)); - deallog << "ndofs: " << dofh.n_dofs() << std::endl; - -#if false - print_dof_points(dofh); - print_dof_indices_on_faces(dofh); -#endif - - // hanging node constraints - AffineConstraints constraints; - DoFTools::make_hanging_node_constraints(dofh, constraints); - constraints.print(deallog.get_file_stream()); - - deallog << "OK" << std::endl; -} +#include "hanging_nodes.h" int @@ -153,6 +51,22 @@ main() initlog(); deallog.push("2d"); - test<2>(); + { + const unsigned int dim = 2; + + const auto subdivided_hyper_cube_with_simplices = + [](Triangulation &tria) { + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + }; + + test({1, 0}, + {0, 0}, + hp::FECollection(Simplex::FE_P(1)), + subdivided_hyper_cube_with_simplices); + test({1, 0}, + {0, 0}, + hp::FECollection(Simplex::FE_P(2)), + subdivided_hyper_cube_with_simplices); + } deallog.pop(); } diff --git a/tests/simplex/hanging_nodes_01.with_simplex_support=on.output b/tests/simplex/hanging_nodes_01.with_simplex_support=on.output index 8034a9e256..eb6b0d13ac 100644 --- a/tests/simplex/hanging_nodes_01.with_simplex_support=on.output +++ b/tests/simplex/hanging_nodes_01.with_simplex_support=on.output @@ -1,5 +1,14 @@ DEAL:2d::ndofs: 7 - 6 2: 0.500000 6 1: 0.500000 + 6 2: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 19 + 12 4: 1.00000 + 14 1: -0.125000 + 14 2: 0.375000 + 14 4: 0.750000 + 17 1: 0.375000 + 17 2: -0.125000 + 17 4: 0.750000 DEAL:2d::OK diff --git a/tests/simplex/hanging_nodes_02.cc b/tests/simplex/hanging_nodes_02.cc index 024848424e..838092d0f8 100644 --- a/tests/simplex/hanging_nodes_02.cc +++ b/tests/simplex/hanging_nodes_02.cc @@ -31,12 +31,11 @@ #include #include -#include - +#include +#include #include #include -#include #include @@ -45,86 +44,7 @@ #include "../tests.h" - -// ----- diagnostics ----- - -template -void -print_dof_indices_on_faces(const DoFHandler &dofh) -{ - std::vector dof_indices; - - for (const auto &cell : dofh.active_cell_iterators()) - for (unsigned int f = 0; f < cell->n_faces(); ++f) - { - const auto &face = cell->face(f); - - Assert(!face->has_children(), ExcInternalError()); - - const unsigned int fe_index = cell->active_fe_index(); - const auto & fe = cell->get_fe(); - - dof_indices.resize(fe.n_dofs_per_face(f)); - face->get_dof_indices(dof_indices, fe_index); - - deallog << "cell:" << cell->active_cell_index() << " face:" << f - << " dofs:"; - for (const auto &i : dof_indices) - deallog << i << " "; - deallog << std::endl; - } -} - - -template -void -print_dof_points(const DoFHandler &dofh) -{ - hp::MappingCollection mapping; - for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i) - mapping.push_back(MappingFE(dofh.get_fe(i))); - - std::vector> points(dofh.n_dofs()); - DoFTools::map_dofs_to_support_points(mapping, dofh, points); - - for (unsigned int i = 0; i < dofh.n_dofs(); ++i) - deallog << "dof:" << i << " point:" << points[i] << std::endl; -} - - -// ----- test ----- - -template -void -test(const hp::FECollection &fes) -{ - // setup grid - Triangulation tria; - GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); - -#if false - GridOut grid_out; - grid_out.write_vtk(tria, deallog.get_file_stream()); -#endif - - DoFHandler dofh(tria); - dofh.begin_active()->set_active_fe_index(1); - - dofh.distribute_dofs(fes); - deallog << "ndofs: " << dofh.n_dofs() << std::endl; - -#if false - print_dof_points(dofh); - print_dof_indices_on_faces(dofh); -#endif - - // hanging node constraints - AffineConstraints constraints; - DoFTools::make_hanging_node_constraints(dofh, constraints); - constraints.print(deallog.get_file_stream()); - - deallog << "OK" << std::endl; -} +#include "hanging_nodes.h" int @@ -133,7 +53,24 @@ main() initlog(); deallog.push("2d"); - test<2>(hp::FECollection<2>(Simplex::FE_P<2>(1), Simplex::FE_P<2>(2))); - test<2>(hp::FECollection<2>(Simplex::FE_P<2>(2), Simplex::FE_P<2>(1))); + { + const unsigned int dim = 2; + + const auto subdivided_hyper_cube_with_simplices = + [](Triangulation &tria) { + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + }; + + test({0, 0}, + {0, 1}, + hp::FECollection(Simplex::FE_P(2), + Simplex::FE_P(1)), + subdivided_hyper_cube_with_simplices); + test({0, 0}, + {0, 1}, + hp::FECollection(Simplex::FE_P(1), + Simplex::FE_P(2)), + subdivided_hyper_cube_with_simplices); + } deallog.pop(); } diff --git a/tests/simplex/hanging_nodes_02.with_simplex_support=on.output b/tests/simplex/hanging_nodes_02.with_simplex_support=on.output index 0bc8c9bb95..bcddbbf901 100644 --- a/tests/simplex/hanging_nodes_02.with_simplex_support=on.output +++ b/tests/simplex/hanging_nodes_02.with_simplex_support=on.output @@ -1,7 +1,7 @@ DEAL:2d::ndofs: 7 - 2 6: 0.500000 2 5: 0.500000 + 2 6: 0.500000 DEAL:2d::OK DEAL:2d::ndofs: 7 5 1: 0.500000 diff --git a/tests/simplex/hanging_nodes_03.cc b/tests/simplex/hanging_nodes_03.cc new file mode 100644 index 0000000000..c0e802b554 --- /dev/null +++ b/tests/simplex/hanging_nodes_03.cc @@ -0,0 +1,76 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + + +// verify hanging node constraints on locally hp-refined simplex mesh +// +// dofs will be enumerated as follows +// scenario 1: +// 9---1---0 +// |\ | +// | \ | +// 6--2,8 3 +// |\ |\ | +// | \| \| +// 4---5---7 + + +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include + +#include "../tests.h" + +#include "hanging_nodes.h" + + +int +main() +{ + initlog(); + + deallog.push("2d"); + { + const unsigned int dim = 2; + + const auto subdivided_hyper_cube_with_simplices = + [](Triangulation &tria) { + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + }; + + test({1, 0}, + {0, 1}, + hp::FECollection(Simplex::FE_P(1), + Simplex::FE_P(2)), + subdivided_hyper_cube_with_simplices); + test({1, 0}, + {0, 1}, + hp::FECollection(Simplex::FE_P(2), + Simplex::FE_P(1)), + subdivided_hyper_cube_with_simplices); + } + deallog.pop(); +} diff --git a/tests/simplex/hanging_nodes_03.with_simplex_support=on.output b/tests/simplex/hanging_nodes_03.with_simplex_support=on.output new file mode 100644 index 0000000000..589460aef2 --- /dev/null +++ b/tests/simplex/hanging_nodes_03.with_simplex_support=on.output @@ -0,0 +1,15 @@ + +DEAL:2d::ndofs: 10 + 2 7: 0.500000 + 2 9: 0.500000 + 8 7: 0.500000 + 8 9: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 16 + 9 1: 0.500000 + 9 2: 0.500000 + 11 1: 0.250000 + 11 2: 0.750000 + 14 1: 0.750000 + 14 2: 0.250000 +DEAL:2d::OK diff --git a/tests/simplex/hanging_nodes_hybrid_01.cc b/tests/simplex/hanging_nodes_hybrid_01.cc new file mode 100644 index 0000000000..6a98030110 --- /dev/null +++ b/tests/simplex/hanging_nodes_hybrid_01.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + + +// verify hanging node constraints on locally h-refined hybrid meshes +/* + * dofs will be enumerated as follows with d = 1 + * scenario 1: scenario 2: scenario 3: + * 2-------3 2-------3 7---8---9 + * | |\ | |\ | | |\ + * | | \ | | 6 | | | \ + * | | \ | |/|\ | | | \ + * | | 4 | 4 | 7 3---4---6 0 + * | | / | |\|/ | | | / + * | | / | | 5 | | | / + * | |/ | |/ | | |/ + * 0-------1 0-------1 1---2---5 + */ + + +#include + +#include + +#include + +#include + +#include "../tests.h" + +#include "hanging_nodes.h" +#include "simplex_grids.h" + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("2d"); + { + const unsigned int dim = 2; + + const auto cube_and_pyramid = [](Triangulation &tria) { + GridGenerator::cube_and_pyramid(tria, 1); + }; + + for (unsigned int d = 1; d <= 2; ++d) + { + deallog << "degree: " << d << std::endl; + test({0, 1}, + {0, 1}, + hp::FECollection(FE_Q(d), Simplex::FE_P(d)), + cube_and_pyramid); + test({1, 0}, + {0, 1}, + hp::FECollection(FE_Q(d), Simplex::FE_P(d)), + cube_and_pyramid); + test({0, 1}, + {1, 0}, + hp::FECollection(Simplex::FE_P(d), FE_Q(d)), + cube_and_pyramid); + test({1, 0}, + {1, 0}, + hp::FECollection(Simplex::FE_P(d), FE_Q(d)), + cube_and_pyramid); + } + } + deallog.pop(); +} diff --git a/tests/simplex/hanging_nodes_hybrid_01.with_simplex_support=on.output b/tests/simplex/hanging_nodes_hybrid_01.with_simplex_support=on.output new file mode 100644 index 0000000000..3ff5c89543 --- /dev/null +++ b/tests/simplex/hanging_nodes_hybrid_01.with_simplex_support=on.output @@ -0,0 +1,55 @@ + +DEAL:2d::degree: 1 +DEAL:2d::ndofs: 8 + 4 1: 0.500000 + 4 3: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 10 + 6 5: 0.500000 + 6 9: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 8 + 3 2: 0.500000 + 3 5: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 10 + 7 0: 0.500000 + 7 1: 0.500000 +DEAL:2d::OK +DEAL:2d::degree: 2 +DEAL:2d::ndofs: 22 + 9 5: 1.00000 + 11 1: 0.375000 + 11 3: -0.125000 + 11 5: 0.750000 + 15 1: -0.125000 + 15 3: 0.375000 + 15 5: 0.750000 +DEAL:2d::OK +DEAL:2d::ndofs: 29 + 14 1: 1.00000 + 15 1: 0.750000 + 15 13: 0.375000 + 15 25: -0.125000 + 26 1: 0.750000 + 26 13: -0.125000 + 26 25: 0.375000 +DEAL:2d::OK +DEAL:2d::ndofs: 22 + 8 3: 1.00000 + 10 3: 0.750000 + 10 7: 0.375000 + 10 13: -0.125000 + 15 3: 0.750000 + 15 7: -0.125000 + 15 13: 0.375000 +DEAL:2d::OK +DEAL:2d::ndofs: 29 + 15 3: 1.00000 + 16 0: 0.375000 + 16 1: -0.125000 + 16 3: 0.750000 + 26 0: -0.125000 + 26 1: 0.375000 + 26 3: 0.750000 +DEAL:2d::OK diff --git a/tests/simplex/hanging_nodes_hybrid_02.cc b/tests/simplex/hanging_nodes_hybrid_02.cc new file mode 100644 index 0000000000..cf7dc743d0 --- /dev/null +++ b/tests/simplex/hanging_nodes_hybrid_02.cc @@ -0,0 +1,79 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + + +// verify hanging node constraints on locally p-refined hybrid meshes +/* + * dofs will be enumerated as follows for degrees = {1, 2} + * scenario 1: scenario 2: + * 2-------3 1---5---8 + * | |\ | |\ + * | | 6 | | \ + * | | \ | | \ + * | 5 4 2 6 3 9 + * | | / | | / + * | | 7 | | / + * | |/ | |/ + * 0-------1 0---4---7 + */ + + +#include + +#include + +#include + +#include + +#include "../tests.h" + +#include "hanging_nodes.h" +#include "simplex_grids.h" + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("2d"); + { + const unsigned int dim = 2; + + const auto cube_and_pyramid = [](Triangulation &tria) { + GridGenerator::cube_and_pyramid(tria, 1); + }; + + for (unsigned int q = 1; q <= 4; ++q) + for (unsigned int p = 1; p <= 2; ++p) + { + if (q == p) + continue; + + deallog << "q_degree: " << q << ", p_degree: " << p << std::endl; + test({0, 0}, + {0, 1}, + hp::FECollection(FE_Q(q), Simplex::FE_P(p)), + cube_and_pyramid); + test({0, 0}, + {1, 0}, + hp::FECollection(Simplex::FE_P(p), FE_Q(q)), + cube_and_pyramid); + } + } + deallog.pop(); +} diff --git a/tests/simplex/hanging_nodes_hybrid_02.with_simplex_support=on.output b/tests/simplex/hanging_nodes_hybrid_02.with_simplex_support=on.output new file mode 100644 index 0000000000..1d51692dfd --- /dev/null +++ b/tests/simplex/hanging_nodes_hybrid_02.with_simplex_support=on.output @@ -0,0 +1,83 @@ + +DEAL:2d::q_degree: 1, p_degree: 2 +DEAL:2d::ndofs: 8 + 5 1: 0.500000 + 5 3: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 8 + 5 1: 0.500000 + 5 3: 0.500000 +DEAL:2d::OK +DEAL:2d::q_degree: 2, p_degree: 1 +DEAL:2d::ndofs: 10 + 3 7: 0.500000 + 3 8: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 10 + 3 7: 0.500000 + 3 8: 0.500000 +DEAL:2d::OK +DEAL:2d::q_degree: 3, p_degree: 1 +DEAL:2d::ndofs: 17 + 4 14: 0.723607 + 4 15: 0.276393 + 5 14: 0.276393 + 5 15: 0.723607 +DEAL:2d::OK +DEAL:2d::ndofs: 17 + 4 14: 0.723607 + 4 15: 0.276393 + 5 14: 0.276393 + 5 15: 0.723607 +DEAL:2d::OK +DEAL:2d::q_degree: 3, p_degree: 2 +DEAL:2d::ndofs: 20 + 4 14: 0.323607 + 4 15: -0.123607 + 4 17: 0.800000 + 5 14: -0.123607 + 5 15: 0.323607 + 5 17: 0.800000 +DEAL:2d::OK +DEAL:2d::ndofs: 20 + 4 14: 0.323607 + 4 15: -0.123607 + 4 17: 0.800000 + 5 14: -0.123607 + 5 15: 0.323607 + 5 17: 0.800000 +DEAL:2d::OK +DEAL:2d::q_degree: 4, p_degree: 1 +DEAL:2d::ndofs: 26 + 5 23: 0.827327 + 5 24: 0.172673 + 6 23: 0.500000 + 6 24: 0.500000 + 7 23: 0.172673 + 7 24: 0.827327 +DEAL:2d::OK +DEAL:2d::ndofs: 26 + 5 23: 0.827327 + 5 24: 0.172673 + 6 23: 0.500000 + 6 24: 0.500000 + 7 23: 0.172673 + 7 24: 0.827327 +DEAL:2d::OK +DEAL:2d::q_degree: 4, p_degree: 2 +DEAL:2d::ndofs: 28 + 5 22: 0.541613 + 5 23: -0.113041 + 5 25: 0.571429 + 6 22: -0.113041 + 6 23: 0.541613 + 6 25: 0.571429 +DEAL:2d::OK +DEAL:2d::ndofs: 28 + 5 22: 0.541613 + 5 23: -0.113041 + 5 25: 0.571429 + 6 22: -0.113041 + 6 23: 0.541613 + 6 25: 0.571429 +DEAL:2d::OK diff --git a/tests/simplex/hanging_nodes_hybrid_03.cc b/tests/simplex/hanging_nodes_hybrid_03.cc new file mode 100644 index 0000000000..63440a7801 --- /dev/null +++ b/tests/simplex/hanging_nodes_hybrid_03.cc @@ -0,0 +1,74 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + + +// verify hanging node constraints on locally hp-refined hybrid meshes + + +#include + +#include + +#include + +#include + +#include "../tests.h" + +#include "hanging_nodes.h" +#include "simplex_grids.h" + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("2d"); + { + const unsigned int dim = 2; + + const auto cube_and_pyramid = [](Triangulation &tria) { + GridGenerator::cube_and_pyramid(tria, 1); + }; + + for (unsigned int q = 1; q <= 4; ++q) + for (unsigned int p = 1; p <= 2; ++p) + { + if (q == p) + continue; + + deallog << "q_degree: " << q << ", p_degree: " << p << std::endl; + test({0, 1}, + {0, 1}, + hp::FECollection(FE_Q(q), Simplex::FE_P(p)), + cube_and_pyramid); + test({1, 0}, + {0, 1}, + hp::FECollection(FE_Q(q), Simplex::FE_P(p)), + cube_and_pyramid); + test({0, 1}, + {1, 0}, + hp::FECollection(Simplex::FE_P(p), FE_Q(q)), + cube_and_pyramid); + test({1, 0}, + {1, 0}, + hp::FECollection(Simplex::FE_P(p), FE_Q(q)), + cube_and_pyramid); + } + } + deallog.pop(); +} diff --git a/tests/simplex/hanging_nodes_hybrid_03.with_simplex_support=on.output b/tests/simplex/hanging_nodes_hybrid_03.with_simplex_support=on.output new file mode 100644 index 0000000000..cfc2558967 --- /dev/null +++ b/tests/simplex/hanging_nodes_hybrid_03.with_simplex_support=on.output @@ -0,0 +1,289 @@ + +DEAL:2d::q_degree: 1, p_degree: 2 +DEAL:2d::ndofs: 17 + 4 1: 0.500000 + 4 3: 0.500000 + 6 1: 0.750000 + 6 3: 0.250000 + 10 1: 0.250000 + 10 3: 0.750000 +DEAL:2d::OK +DEAL:2d::ndofs: 13 + 1 8: 0.500000 + 1 12: 0.500000 + 9 8: 0.500000 + 9 12: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 17 + 4 1: 0.500000 + 4 3: 0.500000 + 6 1: 0.750000 + 6 3: 0.250000 + 10 1: 0.250000 + 10 3: 0.750000 +DEAL:2d::OK +DEAL:2d::ndofs: 13 + 1 8: 0.500000 + 1 12: 0.500000 + 9 8: 0.500000 + 9 12: 0.500000 +DEAL:2d::OK +DEAL:2d::q_degree: 2, p_degree: 1 +DEAL:2d::ndofs: 13 + 3 7: 0.500000 + 3 10: 0.500000 + 8 7: 0.500000 + 8 10: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 26 + 12 0: 0.500000 + 12 1: 0.500000 + 13 0: 0.750000 + 13 1: 0.250000 + 23 0: 0.250000 + 23 1: 0.750000 +DEAL:2d::OK +DEAL:2d::ndofs: 13 + 3 7: 0.500000 + 3 10: 0.500000 + 8 7: 0.500000 + 8 10: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 26 + 12 0: 0.500000 + 12 1: 0.500000 + 13 0: 0.750000 + 13 1: 0.250000 + 23 0: 0.250000 + 23 1: 0.750000 +DEAL:2d::OK +DEAL:2d::q_degree: 3, p_degree: 1 +DEAL:2d::ndofs: 20 + 4 14: 0.723607 + 4 17: 0.276393 + 5 14: 0.276393 + 5 17: 0.723607 + 15 14: 0.500000 + 15 17: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 50 + 19 0: 0.500000 + 19 1: 0.500000 + 20 0: 0.861803 + 20 1: 0.138197 + 21 0: 0.638197 + 21 1: 0.361803 + 42 0: 0.361803 + 42 1: 0.638197 + 43 0: 0.138197 + 43 1: 0.861803 +DEAL:2d::OK +DEAL:2d::ndofs: 20 + 4 14: 0.723607 + 4 17: 0.276393 + 5 14: 0.276393 + 5 17: 0.723607 + 15 14: 0.500000 + 15 17: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 50 + 19 0: 0.500000 + 19 1: 0.500000 + 20 0: 0.861803 + 20 1: 0.138197 + 21 0: 0.638197 + 21 1: 0.361803 + 42 0: 0.361803 + 42 1: 0.638197 + 43 0: 0.138197 + 43 1: 0.861803 +DEAL:2d::OK +DEAL:2d::q_degree: 3, p_degree: 2 +DEAL:2d::ndofs: 29 + 5 4: 1.00000 + 5 14: -0.447214 + 5 20: 0.447214 + 15 4: 1.25000 + 15 14: -0.404508 + 15 20: 0.154508 + 17 4: 0.937500 + 17 14: 0.0716186 + 17 20: -0.00911863 + 22 4: 0.937500 + 22 14: -0.428381 + 22 20: 0.490881 +DEAL:2d::OK +DEAL:2d::ndofs: 53 + 22 3: 1.00000 + 23 0: 0.623607 + 23 1: -0.100000 + 23 3: 0.476393 + 24 0: 0.176393 + 24 1: -0.100000 + 24 3: 0.923607 + 45 0: -0.100000 + 45 1: 0.176393 + 45 3: 0.923607 + 46 0: -0.100000 + 46 1: 0.623607 + 46 3: 0.476393 +DEAL:2d::OK +DEAL:2d::ndofs: 29 + 5 4: 1.00000 + 5 14: -0.447214 + 5 20: 0.447214 + 15 4: 1.25000 + 15 14: -0.404508 + 15 20: 0.154508 + 17 4: 0.937500 + 17 14: 0.0716186 + 17 20: -0.00911863 + 22 4: 0.937500 + 22 14: -0.428381 + 22 20: 0.490881 +DEAL:2d::OK +DEAL:2d::ndofs: 53 + 22 3: 1.00000 + 23 0: 0.623607 + 23 1: -0.100000 + 23 3: 0.476393 + 24 0: 0.176393 + 24 1: -0.100000 + 24 3: 0.923607 + 45 0: -0.100000 + 45 1: 0.176393 + 45 3: 0.923607 + 46 0: -0.100000 + 46 1: 0.623607 + 46 3: 0.476393 +DEAL:2d::OK +DEAL:2d::q_degree: 4, p_degree: 1 +DEAL:2d::ndofs: 29 + 5 23: 0.827327 + 5 26: 0.172673 + 6 23: 0.500000 + 6 26: 0.500000 + 7 23: 0.172673 + 7 26: 0.827327 + 24 23: 0.500000 + 24 26: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 82 + 28 0: 0.500000 + 28 1: 0.500000 + 29 0: 0.913663 + 29 1: 0.0863366 + 30 0: 0.750000 + 30 1: 0.250000 + 31 0: 0.586337 + 31 1: 0.413663 + 67 0: 0.413663 + 67 1: 0.586337 + 68 0: 0.250000 + 68 1: 0.750000 + 69 0: 0.0863366 + 69 1: 0.913663 +DEAL:2d::OK +DEAL:2d::ndofs: 29 + 5 23: 0.827327 + 5 26: 0.172673 + 6 23: 0.500000 + 6 26: 0.500000 + 7 23: 0.172673 + 7 26: 0.827327 + 24 23: 0.500000 + 24 26: 0.500000 +DEAL:2d::OK +DEAL:2d::ndofs: 82 + 28 0: 0.500000 + 28 1: 0.500000 + 29 0: 0.913663 + 29 1: 0.0863366 + 30 0: 0.750000 + 30 1: 0.250000 + 31 0: 0.586337 + 31 1: 0.413663 + 67 0: 0.413663 + 67 1: 0.586337 + 68 0: 0.250000 + 68 1: 0.750000 + 69 0: 0.0863366 + 69 1: 0.913663 +DEAL:2d::OK +DEAL:2d::q_degree: 4, p_degree: 2 +DEAL:2d::ndofs: 38 + 6 5: 1.75000 + 6 23: -0.947822 + 6 29: 0.197822 + 7 5: 1.00000 + 7 23: -0.654654 + 7 29: 0.654654 + 24 5: 1.75000 + 24 23: -0.947822 + 24 29: 0.197822 + 26 5: 1.31250 + 26 23: -0.335866 + 26 29: 0.0233665 + 31 5: 1.31250 + 31 23: -0.835866 + 31 29: 0.523366 +DEAL:2d::OK +DEAL:2d::ndofs: 85 + 31 3: 1.00000 + 32 0: 0.755898 + 32 1: -0.0714286 + 32 3: 0.315530 + 33 0: 0.375000 + 33 1: -0.125000 + 33 3: 0.750000 + 34 0: 0.101245 + 34 1: -0.0714286 + 34 3: 0.970184 + 70 0: -0.0714286 + 70 1: 0.101245 + 70 3: 0.970184 + 71 0: -0.125000 + 71 1: 0.375000 + 71 3: 0.750000 + 72 0: -0.0714286 + 72 1: 0.755898 + 72 3: 0.315530 +DEAL:2d::OK +DEAL:2d::ndofs: 38 + 6 5: 1.75000 + 6 23: -0.947822 + 6 29: 0.197822 + 7 5: 1.00000 + 7 23: -0.654654 + 7 29: 0.654654 + 24 5: 1.75000 + 24 23: -0.947822 + 24 29: 0.197822 + 26 5: 1.31250 + 26 23: -0.335866 + 26 29: 0.0233665 + 31 5: 1.31250 + 31 23: -0.835866 + 31 29: 0.523366 +DEAL:2d::OK +DEAL:2d::ndofs: 85 + 31 3: 1.00000 + 32 0: 0.755898 + 32 1: -0.0714286 + 32 3: 0.315530 + 33 0: 0.375000 + 33 1: -0.125000 + 33 3: 0.750000 + 34 0: 0.101245 + 34 1: -0.0714286 + 34 3: 0.970184 + 70 0: -0.0714286 + 70 1: 0.101245 + 70 3: 0.970184 + 71 0: -0.125000 + 71 1: 0.375000 + 71 3: 0.750000 + 72 0: -0.0714286 + 72 1: 0.755898 + 72 3: 0.315530 +DEAL:2d::OK diff --git a/tests/simplex/simplex_grids.h b/tests/simplex/simplex_grids.h index 8d0e42dc1a..48e08c6697 100644 --- a/tests/simplex/simplex_grids.h +++ b/tests/simplex/simplex_grids.h @@ -376,5 +376,60 @@ namespace dealii AssertThrow(false, ExcNotImplemented()) } } + + + + /** + * Adds a simplex cell to face @p face_no of a hyper_cube cell. + */ + template + void + cube_and_pyramid(Triangulation &tria, + const unsigned int face_no = 1) + { + Assert(face_no % 2 == 1, + ExcMessage("Only works for odd face numbers. " + "GridReordering::reorder_cells() is not prepared for " + "simplices yet (uses GeometryInfo).")); + + // create cube + Triangulation tria_cube; + GridGenerator::hyper_cube(tria_cube); + const auto cube = tria_cube.begin_active(); + AssertIndexRange(face_no, cube->n_faces()); + + // extract vertices from specified face, store their midpoint + std::vector> vertices; + Point midpoint; + const auto shared_face = cube->face(face_no); + for (unsigned int v = 0; v < shared_face->n_vertices(); ++v) + { + const auto &vertex = shared_face->vertex(v); + vertices.push_back(vertex); + midpoint += vertex; + } + midpoint /= vertices.size(); + + // add simplex cell in coordinate direction + const unsigned int coordinate = face_no / 2; +#ifdef DEBUG + // all vertices should be in a plane + for (const auto &vertex : vertices) + Assert(midpoint(coordinate) == vertex(coordinate), ExcInternalError()); +#endif + + // add another vertex as tip of triangle/pyramid + Point tip = midpoint; + tip(coordinate) += (face_no % 2 == 1) ? 1. : -1.; + vertices.push_back(tip); + + CellData simplex(vertices.size()); + std::iota(simplex.vertices.begin(), simplex.vertices.end(), 0); + + Triangulation tria_simplex; + tria_simplex.create_triangulation(vertices, {simplex}, SubCellData()); + + GridGenerator::merge_triangulations(tria_cube, tria_simplex, tria); + } } // namespace GridGenerator } // namespace dealii