From 4286f162cf62fbe814b9577eee3652214ca86cf5 Mon Sep 17 00:00:00 2001 From: peterrum Date: Sat, 27 Jul 2019 10:45:56 +0200 Subject: [PATCH] Enable periodicity for 1D --- doc/news/changes/minor/20190727PeterMunch | 3 + include/deal.II/dofs/dof_accessor.templates.h | 15 +- source/dofs/dof_tools_constraints.inst.in | 2 - source/grid/tria.cc | 256 ++++++++++-------- tests/bits/periodicity_01.cc | 15 +- tests/bits/periodicity_01.output | 4 + tests/grid/periodicity_01.cc | 3 +- tests/grid/periodicity_01.output | 1 + 8 files changed, 175 insertions(+), 124 deletions(-) create mode 100644 doc/news/changes/minor/20190727PeterMunch diff --git a/doc/news/changes/minor/20190727PeterMunch b/doc/news/changes/minor/20190727PeterMunch new file mode 100644 index 0000000000..7aa7d1c1f6 --- /dev/null +++ b/doc/news/changes/minor/20190727PeterMunch @@ -0,0 +1,3 @@ +New: Enable periodicity for 1D triangulations (only for dealii::DoFHandler). +
+(Peter Munch, 2019/07/27) diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index f8a4b7a690..be8285f8ba 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -35,6 +35,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -2565,8 +2566,11 @@ inline unsigned int DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access>:: n_active_fe_indices() const { - Assert(false, ExcNotImplemented()); - return numbers::invalid_unsigned_int; + Assert((std::is_same, + dealii::DoFHandler<1, spacedim>>::value == true), + ExcNotImplemented()); + + return 1; } @@ -2578,8 +2582,11 @@ inline unsigned int DoFAccessor<0, DoFHandlerType<1, spacedim>, level_dof_access>:: nth_active_fe_index(const unsigned int /*n*/) const { - Assert(false, ExcNotImplemented()); - return numbers::invalid_unsigned_int; + Assert((std::is_same, + dealii::DoFHandler<1, spacedim>>::value == true), + ExcNotImplemented()); + + return 0; } diff --git a/source/dofs/dof_tools_constraints.inst.in b/source/dofs/dof_tools_constraints.inst.in index 685a98d022..931ef21f26 100644 --- a/source/dofs/dof_tools_constraints.inst.in +++ b/source/dofs/dof_tools_constraints.inst.in @@ -49,7 +49,6 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) AffineConstraints> &); #endif -#if deal_II_dimension != 1 template void DoFTools::make_periodicity_constraints( const DH::face_iterator &, const DH::face_iterator &, @@ -83,7 +82,6 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) const int, AffineConstraints &, const ComponentMask &); -#endif } for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 31615cb9a6..17efe915c0 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -950,7 +950,12 @@ namespace Assert(face_1->at_boundary() && face_2->at_boundary(), ExcMessage("Periodic faces must be on the boundary")); - Assert(std::abs(cell_1->level() - cell_2->level()) < 2, ExcInternalError()); + // Check if the requirement that each edge can only have at most one hanging + // node, and as a consequence neighboring cells can differ by at most + // one refinement level is enforced. In 1d, there are no hanging nodes and + // so neighboring cells can differ by more than one refinement level. + Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2, + ExcInternalError()); // insert periodic face pair for both cells using CellFace = @@ -968,127 +973,24 @@ namespace Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError()); periodic_face_map.insert(periodic_faces); - // A lookup table on how to go through the child cells depending on the - // orientation: - // see Documentation of GeometryInfo for details - - static const int lookup_table_2d[2][2] = - // flip: + if (dim == 1) { - {0, 1}, // false - {1, 0} // true - }; - - static const int lookup_table_3d[2][2][2][4] = - // orientation flip rotation - {{{ - {0, 2, 1, 3}, // false false false - {2, 3, 0, 1} // false false true - }, - { - {3, 1, 2, 0}, // false true false - {1, 0, 3, 2} // false true true - }}, - {{ - {0, 1, 2, 3}, // true false false - {1, 3, 0, 2} // true false true - }, - { - {3, 2, 1, 0}, // true true false - {2, 0, 3, 1} // true true true - }}}; - - if (cell_1->has_children()) - { - if (cell_2->has_children()) + if (cell_1->has_children()) { - // In the case that both faces have children, we loop over all - // children and apply update_periodic_face_map_recursively - // recursively: - - Assert(face_1->n_children() == - GeometryInfo::max_children_per_face && - face_2->n_children() == - GeometryInfo::max_children_per_face, - ExcNotImplemented()); - - for (unsigned int i = 0; - i < GeometryInfo::max_children_per_face; - ++i) + if (cell_2->has_children()) { - // Lookup the index for the second face - unsigned int j = 0; - switch (dim) - { - case 2: - j = lookup_table_2d[face_flip][i]; - break; - case 3: - j = lookup_table_3d[face_orientation][face_flip] - [face_rotation][i]; - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - - // find subcell ids that belong to the subface indices - unsigned int child_cell_1 = - GeometryInfo::child_cell_on_face( - cell_1->refinement_case(), - n_face_1, - i, - cell_1->face_orientation(n_face_1), - cell_1->face_flip(n_face_1), - cell_1->face_rotation(n_face_1), - face_1->refinement_case()); - unsigned int child_cell_2 = - GeometryInfo::child_cell_on_face( - cell_2->refinement_case(), - n_face_2, - j, - cell_2->face_orientation(n_face_2), - cell_2->face_flip(n_face_2), - cell_2->face_rotation(n_face_2), - face_2->refinement_case()); - - Assert(cell_1->child(child_cell_1)->face(n_face_1) == - face_1->child(i), - ExcInternalError()); - Assert(cell_2->child(child_cell_2)->face(n_face_2) == - face_2->child(j), - ExcInternalError()); - - // precondition: subcell has the same orientation as cell (so - // that the face numbers coincide) recursive call update_periodic_face_map_recursively( - cell_1->child(child_cell_1), - cell_2->child(child_cell_2), + cell_1->child(n_face_1), + cell_2->child(n_face_2), n_face_1, n_face_2, orientation, periodic_face_map); } - } - else // only face_1 has children - { - for (unsigned int i = 0; - i < GeometryInfo::max_children_per_face; - ++i) + else // only face_1 has children { - // find subcell ids that belong to the subface indices - unsigned int child_cell_1 = - GeometryInfo::child_cell_on_face( - cell_1->refinement_case(), - n_face_1, - i, - cell_1->face_orientation(n_face_1), - cell_1->face_flip(n_face_1), - cell_1->face_rotation(n_face_1), - face_1->refinement_case()); - - // recursive call update_periodic_face_map_recursively( - cell_1->child(child_cell_1), + cell_1->child(n_face_1), cell_2, n_face_1, n_face_2, @@ -1097,6 +999,138 @@ namespace } } } + else // dim == 2 || dim == 3 + { + // A lookup table on how to go through the child cells depending on the + // orientation: + // see Documentation of GeometryInfo for details + + static const int lookup_table_2d[2][2] = + // flip: + { + {0, 1}, // false + {1, 0} // true + }; + + static const int lookup_table_3d[2][2][2][4] = + // orientation flip rotation + {{{ + {0, 2, 1, 3}, // false false false + {2, 3, 0, 1} // false false true + }, + { + {3, 1, 2, 0}, // false true false + {1, 0, 3, 2} // false true true + }}, + {{ + {0, 1, 2, 3}, // true false false + {1, 3, 0, 2} // true false true + }, + { + {3, 2, 1, 0}, // true true false + {2, 0, 3, 1} // true true true + }}}; + + if (cell_1->has_children()) + { + if (cell_2->has_children()) + { + // In the case that both faces have children, we loop over all + // children and apply update_periodic_face_map_recursively + // recursively: + + Assert(face_1->n_children() == + GeometryInfo::max_children_per_face && + face_2->n_children() == + GeometryInfo::max_children_per_face, + ExcNotImplemented()); + + for (unsigned int i = 0; + i < GeometryInfo::max_children_per_face; + ++i) + { + // Lookup the index for the second face + unsigned int j = 0; + switch (dim) + { + case 2: + j = lookup_table_2d[face_flip][i]; + break; + case 3: + j = lookup_table_3d[face_orientation][face_flip] + [face_rotation][i]; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + + // find subcell ids that belong to the subface indices + unsigned int child_cell_1 = + GeometryInfo::child_cell_on_face( + cell_1->refinement_case(), + n_face_1, + i, + cell_1->face_orientation(n_face_1), + cell_1->face_flip(n_face_1), + cell_1->face_rotation(n_face_1), + face_1->refinement_case()); + unsigned int child_cell_2 = + GeometryInfo::child_cell_on_face( + cell_2->refinement_case(), + n_face_2, + j, + cell_2->face_orientation(n_face_2), + cell_2->face_flip(n_face_2), + cell_2->face_rotation(n_face_2), + face_2->refinement_case()); + + Assert(cell_1->child(child_cell_1)->face(n_face_1) == + face_1->child(i), + ExcInternalError()); + Assert(cell_2->child(child_cell_2)->face(n_face_2) == + face_2->child(j), + ExcInternalError()); + + // precondition: subcell has the same orientation as cell + // (so that the face numbers coincide) recursive call + update_periodic_face_map_recursively( + cell_1->child(child_cell_1), + cell_2->child(child_cell_2), + n_face_1, + n_face_2, + orientation, + periodic_face_map); + } + } + else // only face_1 has children + { + for (unsigned int i = 0; + i < GeometryInfo::max_children_per_face; + ++i) + { + // find subcell ids that belong to the subface indices + unsigned int child_cell_1 = + GeometryInfo::child_cell_on_face( + cell_1->refinement_case(), + n_face_1, + i, + cell_1->face_orientation(n_face_1), + cell_1->face_flip(n_face_1), + cell_1->face_rotation(n_face_1), + face_1->refinement_case()); + + // recursive call + update_periodic_face_map_recursively( + cell_1->child(child_cell_1), + cell_2, + n_face_1, + n_face_2, + orientation, + periodic_face_map); + } + } + } + } } diff --git a/tests/bits/periodicity_01.cc b/tests/bits/periodicity_01.cc index 5af7e592b7..e381160fc0 100644 --- a/tests/bits/periodicity_01.cc +++ b/tests/bits/periodicity_01.cc @@ -43,11 +43,12 @@ test() Triangulation triangulation; std::vector repetitions(dim, 1); repetitions[0] = 2; - GridGenerator::subdivided_hyper_rectangle(triangulation, - repetitions, - Point(), - (dim == 2 ? Point(2, 1) : - Point(2, 1, 1))); + GridGenerator::subdivided_hyper_rectangle( + triangulation, + repetitions, + Point(), + (dim == 1 ? Point(2) : + (dim == 2 ? Point(2, 1) : Point(2, 1, 1)))); FE_Q fe(1); DoFHandler dof_handler(triangulation); @@ -57,6 +58,8 @@ test() DoFTools::make_periodicity_constraints(dof_handler.begin(0)->face(0), (++dof_handler.begin(0))->face(1), cm); + + deallog << "dim " << std::to_string(dim) << ":" << std::endl; cm.print(deallog.get_file_stream()); } @@ -66,7 +69,7 @@ int main() { initlog(); - + test<1>(); test<2>(); test<3>(); return 0; diff --git a/tests/bits/periodicity_01.output b/tests/bits/periodicity_01.output index eeb0da966c..97d6535c96 100644 --- a/tests/bits/periodicity_01.output +++ b/tests/bits/periodicity_01.output @@ -1,6 +1,10 @@ +DEAL::dim 1: + 2 0: 1.00000 +DEAL::dim 2: 4 0: 1.00000 5 2: 1.00000 +DEAL::dim 3: 8 0: 1.00000 9 2: 1.00000 10 4: 1.00000 diff --git a/tests/grid/periodicity_01.cc b/tests/grid/periodicity_01.cc index 588e79f458..cc5edd7229 100644 --- a/tests/grid/periodicity_01.cc +++ b/tests/grid/periodicity_01.cc @@ -67,7 +67,7 @@ test() unsigned int neighbor_count = 0; for (const auto &cell : tria.active_cell_iterators()) - for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) if (cell->has_periodic_neighbor(f)) ++neighbor_count; @@ -85,6 +85,7 @@ int main() { initlog(); + test<1>(); test<2>(); test<3>(); return 0; diff --git a/tests/grid/periodicity_01.output b/tests/grid/periodicity_01.output index cf04d7bf87..aee6245e1e 100644 --- a/tests/grid/periodicity_01.output +++ b/tests/grid/periodicity_01.output @@ -1,3 +1,4 @@ +DEAL::Found 2 faces with periodic neighbor DEAL::Found 48 faces with periodic neighbor DEAL::Found 1280 faces with periodic neighbor -- 2.39.5