for (const auto &cell : tria.cell_iterators_on_level(0))
{
// Save cell data
- CellData<dim> cell_data;
+ CellData<dim> cell_data(cell->n_vertices());
for (const unsigned int cell_vertex_n : cell->vertex_indices())
{
Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <vector>
+
+#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 <int dim>
+void
+refine(const std::vector<unsigned int> &n_refinements, Triangulation<dim> &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<unsigned int>(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 <int dim>
+void
+set_active_fe_indices(const std::vector<unsigned int> &fe_indices,
+ DoFHandler<dim> & 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 <int dim>
+void
+print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
+{
+ std::vector<types::global_dof_index> 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 <int dim>
+void
+print_dof_points(const DoFHandler<dim> &dofh)
+{
+ hp::MappingCollection<dim> mapping;
+ for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i)
+ mapping.push_back(MappingFE<dim>(dofh.get_fe(i)));
+
+ std::vector<Point<dim>> 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 <int dim>
+void
+test(const std::vector<unsigned int> & n_refinements,
+ const std::vector<unsigned int> & fe_indices,
+ const hp::FECollection<dim> & fe_collection,
+ const std::function<void(Triangulation<dim> &)> &grid_generator)
+{
+ // setup grid
+ Triangulation<dim> 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<dim> 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<double> constraints;
+ DoFTools::make_hanging_node_constraints(dofh, constraints);
+ constraints.close();
+ constraints.print(deallog.get_file_stream());
+
+ deallog << "OK" << std::endl;
+}
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/mapping_fe.h>
-
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include "../tests.h"
-
-// ----- diagnostics -----
-
-template <int dim>
-void
-print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
-{
- std::vector<types::global_dof_index> 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 <int dim>
-void
-print_dof_points(const DoFHandler<dim> &dofh)
-{
- std::vector<Point<dim>> points(dofh.n_dofs());
- DoFTools::map_dofs_to_support_points(MappingFE<dim>(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 <int dim>
-void
-test()
-{
- // setup grid
- Triangulation<dim> 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<dim> dofh(tria);
- dofh.distribute_dofs(Simplex::FE_P<dim>(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<double> constraints;
- DoFTools::make_hanging_node_constraints(dofh, constraints);
- constraints.print(deallog.get_file_stream());
-
- deallog << "OK" << std::endl;
-}
+#include "hanging_nodes.h"
int
initlog();
deallog.push("2d");
- test<2>();
+ {
+ const unsigned int dim = 2;
+
+ const auto subdivided_hyper_cube_with_simplices =
+ [](Triangulation<dim> &tria) {
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+ };
+
+ test<dim>({1, 0},
+ {0, 0},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(1)),
+ subdivided_hyper_cube_with_simplices);
+ test<dim>({1, 0},
+ {0, 0},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(2)),
+ subdivided_hyper_cube_with_simplices);
+ }
deallog.pop();
}
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
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/mapping_fe.h>
-
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/hp/fe_collection.h>
-#include <deal.II/hp/mapping_collection.h>
#include <deal.II/lac/affine_constraints.h>
#include "../tests.h"
-
-// ----- diagnostics -----
-
-template <int dim>
-void
-print_dof_indices_on_faces(const DoFHandler<dim> &dofh)
-{
- std::vector<types::global_dof_index> 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 <int dim>
-void
-print_dof_points(const DoFHandler<dim> &dofh)
-{
- hp::MappingCollection<dim> mapping;
- for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i)
- mapping.push_back(MappingFE<dim>(dofh.get_fe(i)));
-
- std::vector<Point<dim>> 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 <int dim>
-void
-test(const hp::FECollection<dim> &fes)
-{
- // setup grid
- Triangulation<dim> 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<dim> 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<double> constraints;
- DoFTools::make_hanging_node_constraints(dofh, constraints);
- constraints.print(deallog.get_file_stream());
-
- deallog << "OK" << std::endl;
-}
+#include "hanging_nodes.h"
int
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<dim> &tria) {
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+ };
+
+ test<dim>({0, 0},
+ {0, 1},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(2),
+ Simplex::FE_P<dim>(1)),
+ subdivided_hyper_cube_with_simplices);
+ test<dim>({0, 0},
+ {0, 1},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(1),
+ Simplex::FE_P<dim>(2)),
+ subdivided_hyper_cube_with_simplices);
+ }
deallog.pop();
}
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+
+#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<dim> &tria) {
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+ };
+
+ test<dim>({1, 0},
+ {0, 1},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(1),
+ Simplex::FE_P<dim>(2)),
+ subdivided_hyper_cube_with_simplices);
+ test<dim>({1, 0},
+ {0, 1},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(2),
+ Simplex::FE_P<dim>(1)),
+ subdivided_hyper_cube_with_simplices);
+ }
+ deallog.pop();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/fe/fe_q.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include <vector>
+
+#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<dim> &tria) {
+ GridGenerator::cube_and_pyramid(tria, 1);
+ };
+
+ for (unsigned int d = 1; d <= 2; ++d)
+ {
+ deallog << "degree: " << d << std::endl;
+ test<dim>({0, 1},
+ {0, 1},
+ hp::FECollection<dim>(FE_Q<dim>(d), Simplex::FE_P<dim>(d)),
+ cube_and_pyramid);
+ test<dim>({1, 0},
+ {0, 1},
+ hp::FECollection<dim>(FE_Q<dim>(d), Simplex::FE_P<dim>(d)),
+ cube_and_pyramid);
+ test<dim>({0, 1},
+ {1, 0},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(d), FE_Q<dim>(d)),
+ cube_and_pyramid);
+ test<dim>({1, 0},
+ {1, 0},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(d), FE_Q<dim>(d)),
+ cube_and_pyramid);
+ }
+ }
+ deallog.pop();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/fe/fe_q.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include <vector>
+
+#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<dim> &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<dim>({0, 0},
+ {0, 1},
+ hp::FECollection<dim>(FE_Q<dim>(q), Simplex::FE_P<dim>(p)),
+ cube_and_pyramid);
+ test<dim>({0, 0},
+ {1, 0},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(p), FE_Q<dim>(q)),
+ cube_and_pyramid);
+ }
+ }
+ deallog.pop();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/fe/fe_q.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/simplex/fe_lib.h>
+
+#include <vector>
+
+#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<dim> &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<dim>({0, 1},
+ {0, 1},
+ hp::FECollection<dim>(FE_Q<dim>(q), Simplex::FE_P<dim>(p)),
+ cube_and_pyramid);
+ test<dim>({1, 0},
+ {0, 1},
+ hp::FECollection<dim>(FE_Q<dim>(q), Simplex::FE_P<dim>(p)),
+ cube_and_pyramid);
+ test<dim>({0, 1},
+ {1, 0},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(p), FE_Q<dim>(q)),
+ cube_and_pyramid);
+ test<dim>({1, 0},
+ {1, 0},
+ hp::FECollection<dim>(Simplex::FE_P<dim>(p), FE_Q<dim>(q)),
+ cube_and_pyramid);
+ }
+ }
+ deallog.pop();
+}
--- /dev/null
+
+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
AssertThrow(false, ExcNotImplemented())
}
}
+
+
+
+ /**
+ * Adds a simplex cell to face @p face_no of a hyper_cube cell.
+ */
+ template <int dim, int spacedim>
+ void
+ cube_and_pyramid(Triangulation<dim, spacedim> &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<dim, spacedim> 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<Point<spacedim>> vertices;
+ Point<spacedim> 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<spacedim> tip = midpoint;
+ tip(coordinate) += (face_no % 2 == 1) ? 1. : -1.;
+ vertices.push_back(tip);
+
+ CellData<dim> simplex(vertices.size());
+ std::iota(simplex.vertices.begin(), simplex.vertices.end(), 0);
+
+ Triangulation<dim, spacedim> tria_simplex;
+ tria_simplex.create_triangulation(vertices, {simplex}, SubCellData());
+
+ GridGenerator::merge_triangulations(tria_cube, tria_simplex, tria);
+ }
} // namespace GridGenerator
} // namespace dealii