From: Peter Munch Date: Tue, 10 Oct 2023 06:00:15 +0000 (+0200) Subject: Fix ReferenceCell::standard_vs_true_line_orientation() for simplices X-Git-Tag: relicensing~415^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16119%2Fhead;p=dealii.git Fix ReferenceCell::standard_vs_true_line_orientation() for simplices --- diff --git a/doc/news/changes/minor/20231010Munch b/doc/news/changes/minor/20231010Munch new file mode 100644 index 0000000000..01a1ed7593 --- /dev/null +++ b/doc/news/changes/minor/20231010Munch @@ -0,0 +1,4 @@ +Fixed: The querying of line orientations in the case +of tetrahedra has been fixed. +
+(Peter Munch, David Wells, 2023/10/10) diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 11bcb633b6..8167ab01b1 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -522,6 +522,7 @@ public: */ bool standard_vs_true_line_orientation(const unsigned int line, + const unsigned int face, const unsigned char face_orientation, const bool line_orientation) const; @@ -2598,6 +2599,7 @@ ReferenceCell::n_face_orientations(const unsigned int face_no) const inline bool ReferenceCell::standard_vs_true_line_orientation( const unsigned int line, + const unsigned int face, const unsigned char combined_face_orientation, const bool line_orientation) const { @@ -2610,6 +2612,26 @@ ReferenceCell::standard_vs_true_line_orientation( return (line_orientation == bool_table[line / 2][combined_face_orientation]); } + else if (*this == ReferenceCells::Tetrahedron) + { + static constexpr unsigned int X = numbers::invalid_unsigned_int; + static constexpr dealii::ndarray combined_lines{ + {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}}; + + const auto combined_line = combined_lines[face][line]; + + Assert(combined_line != X, + ExcMessage( + "This function can only be called for following face-line " + "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),")); + + static constexpr dealii::ndarray bool_table{ + {{{false, true, false, true, false, true}}, + {{true, false, true, false, true, false}}}}; + + return (line_orientation == + bool_table[combined_line][combined_face_orientation]); + } else // TODO: This might actually be wrong for some of the other // kinds of objects. We should check this diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index 846fe46c5b..e3a5a0061d 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -854,6 +854,7 @@ namespace internal // Then query how that line is oriented within that face: return accessor.reference_cell().standard_vs_true_line_orientation( line_index, + face_index, combined_face_orientation(accessor, face_index), accessor.quad(face_index)->line_orientation(line_within_face_index)); } @@ -1170,13 +1171,16 @@ namespace internal const auto quad = cell.quad(f); const std::array my_orientations{ {ref_cell.standard_vs_true_line_orientation( - 0, orientation, quad->line_orientation(my_indices[0])), + 0, f, orientation, quad->line_orientation(my_indices[0])), ref_cell.standard_vs_true_line_orientation( - 1, orientation, quad->line_orientation(my_indices[1])), + 1, f, orientation, quad->line_orientation(my_indices[1])), ref_cell.standard_vs_true_line_orientation( - 2, orientation, quad->line_orientation(my_indices[2])), + 2, f, orientation, quad->line_orientation(my_indices[2])), ref_cell.standard_vs_true_line_orientation( - 3, orientation, quad->line_orientation(my_indices[3]))}}; + 3, + f, + orientation, + quad->line_orientation(my_indices[3]))}}; for (unsigned int l = 0; l < 4; ++l) line_orientations[4 * (f - 4) + l] = my_orientations[l]; } @@ -1193,9 +1197,12 @@ namespace internal const auto quad = cell.quad(f); const std::array my_orientations{ {ref_cell.standard_vs_true_line_orientation( - 0, orientation, quad->line_orientation(my_indices[0])), + 0, f, orientation, quad->line_orientation(my_indices[0])), ref_cell.standard_vs_true_line_orientation( - 1, orientation, quad->line_orientation(my_indices[1]))}}; + 1, + f, + orientation, + quad->line_orientation(my_indices[1]))}}; line_orientations[8 + f] = my_orientations[0]; line_orientations[10 + f] = my_orientations[1]; } @@ -1214,27 +1221,33 @@ namespace internal ref_cell.standard_to_real_face_line(2, 1, orientations[1]), ref_cell.standard_to_real_face_line(1, 2, orientations[2])}}; line_orientations[0] = ref_cell.standard_vs_true_line_orientation( + 0, 0, orientations[0], cell.quad(0)->line_orientation(my_indices[0])); line_orientations[1] = ref_cell.standard_vs_true_line_orientation( 1, + 0, orientations[0], cell.quad(0)->line_orientation(my_indices[1])); line_orientations[2] = ref_cell.standard_vs_true_line_orientation( 2, + 0, orientations[0], cell.quad(0)->line_orientation(my_indices[2])); line_orientations[3] = ref_cell.standard_vs_true_line_orientation( + 1, 1, orientations[1], cell.quad(1)->line_orientation(my_indices[3])); line_orientations[4] = ref_cell.standard_vs_true_line_orientation( 2, + 1, orientations[1], cell.quad(1)->line_orientation(my_indices[4])); line_orientations[5] = ref_cell.standard_vs_true_line_orientation( 1, + 2, orientations[2], cell.quad(2)->line_orientation(my_indices[5])); } diff --git a/tests/simplex/orientation_04.cc b/tests/simplex/orientation_04.cc new file mode 100644 index 0000000000..fac1116f44 --- /dev/null +++ b/tests/simplex/orientation_04.cc @@ -0,0 +1,187 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2022 - 2022 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. + * + * --------------------------------------------------------------------- + */ + +// Test a mesh with two tetrahedra for all possible orientations. Similar to +// orientation_02 but also checks that line orientations are correct. + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +std::tuple +create_triangulation(const std::vector> &vertices_, + const std::vector> &cell_data_, + const unsigned int face_n, + const unsigned int n_permuations, + Triangulation<3> &tria) +{ + const ReferenceCell ref_cell = ReferenceCells::Tetrahedron; + auto vertices = vertices_; + auto cell_data = cell_data_; + + Point<3> extra_vertex; + for (unsigned int i = 0; i < 3; ++i) + extra_vertex += ref_cell.template vertex<3>(ref_cell.face_to_cell_vertices( + face_n, i, ReferenceCell::default_combined_face_orientation())); + + extra_vertex /= 3.0; + extra_vertex += ref_cell.template unit_normal_vectors<3>(face_n); + + vertices.push_back(extra_vertex); + + cell_data.emplace_back(); + cell_data.back().vertices.resize(0); + for (unsigned int i = 0; i < 3; ++i) + cell_data.back().vertices.push_back(ref_cell.face_to_cell_vertices( + face_n, i, ref_cell.default_combined_face_orientation())); + cell_data.back().vertices.push_back(ref_cell.n_vertices()); + std::sort(cell_data.back().vertices.begin(), cell_data.back().vertices.end()); + + unsigned int permutation_n = 0; + do + { + tria.clear(); + tria.create_triangulation(vertices, cell_data, SubCellData()); + ++permutation_n; + } + while ((permutation_n < n_permuations) && + std::next_permutation(cell_data.back().vertices.begin(), + cell_data.back().vertices.end())); + + const auto cell = tria.begin(); + + const auto face = cell->face(face_n); + + auto ncell = tria.begin(); + ncell++; + ncell->face(face_n); + + unsigned int nf = 0; + for (; nf < ref_cell.n_faces(); ++nf) + if (ncell->face(nf) == face) + break; + + return {nf, ncell->combined_face_orientation(nf)}; +} + + + +void +test() +{ + const unsigned int dim = 3; + + double previous_error = 1.0; + + for (unsigned int f = 0; f < 4; ++f) + { + for (unsigned int r = 0; r < 24; ++r) + { + unsigned int orientation = r; + unsigned int face_no = f; + + Triangulation tria; + + Triangulation<3> dummy; + GridGenerator::reference_cell(dummy, ReferenceCells::Tetrahedron); + + auto vertices = dummy.get_vertices(); + + std::vector> cells; + + { + CellData<3> cell; + cell.vertices = {0, 1, 2, 3}; + cell.material_id = 0; + cells.push_back(cell); + } + + std::tie(face_no, orientation) = + create_triangulation(vertices, cells, face_no, r, tria); + + bool success = true; + + auto cell = tria.begin(); + cell++; + + std::vector verticess; + + for (const auto v : cell->vertex_indices()) + verticess.emplace_back(cell->vertex_index(v)); + + for (unsigned int ll = 0; ll < 3; ++ll) + { + const unsigned int l = + cell->reference_cell().face_to_cell_lines(face_no, ll, 1); + + const auto orientation_exp = cell->line_orientation(l); + + std::pair p0; + p0.first = + verticess[cell->reference_cell().line_to_cell_vertices(l, 0)]; + p0.second = + verticess[cell->reference_cell().line_to_cell_vertices(l, 1)]; + + std::pair p1; + p1.first = cell->line(l)->vertex_index(0); + p1.second = cell->line(l)->vertex_index(1); + + if (orientation_exp == false) + std::swap(p1.first, p1.second); + + success &= (p0 == p1); + } + + if (success) + deallog << "x "; + else + deallog << "o "; + } + } + + deallog << std::endl; +} + +int +main() +{ + initlog(); + + test(); +} diff --git a/tests/simplex/orientation_04.output b/tests/simplex/orientation_04.output new file mode 100644 index 0000000000..133d3507c9 --- /dev/null +++ b/tests/simplex/orientation_04.output @@ -0,0 +1,2 @@ + +DEAL::x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x