From 6f55414f972df267f4791cfda2c198b22ecf25f7 Mon Sep 17 00:00:00 2001 From: Johannes Heinz Date: Wed, 25 Jan 2023 14:38:02 +0100 Subject: [PATCH] Get mapped vertices of face --- doc/news/changes/minor/20230125Heinz | 3 + include/deal.II/fe/mapping.h | 13 +++ source/fe/mapping.cc | 28 +++++++ .../mappings/mapping_get_vertices_on_face.cc | 84 +++++++++++++++++++ .../mapping_get_vertices_on_face.output | 7 ++ 5 files changed, 135 insertions(+) create mode 100644 doc/news/changes/minor/20230125Heinz create mode 100644 tests/mappings/mapping_get_vertices_on_face.cc create mode 100644 tests/mappings/mapping_get_vertices_on_face.output diff --git a/doc/news/changes/minor/20230125Heinz b/doc/news/changes/minor/20230125Heinz new file mode 100644 index 0000000000..c4bed5b02a --- /dev/null +++ b/doc/news/changes/minor/20230125Heinz @@ -0,0 +1,3 @@ +Improved: Mapping can now return vertices on a face of a cell. +
+(Johannes Heinz, 2023/01/25) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 9c2568abc6..99c3fc2109 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -347,6 +347,19 @@ public: get_vertices( const typename Triangulation::cell_iterator &cell) const; + /** + * Return the mapped vertices of a face. + * + * Same as above but working on a given face of a cell. + * + * @param[in] cell The cell containing the face. + * @param[in] face_no The number of the face within the cell. + */ + boost::container::small_vector, + GeometryInfo::vertices_per_face> + get_vertices(const typename Triangulation::cell_iterator &cell, + const unsigned int face_no) const; + /** * Return the mapped center of a cell. * diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index dc87105a03..f4f7fdfec4 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -52,6 +52,34 @@ Mapping::get_vertices( +template +boost::container::small_vector, + GeometryInfo::vertices_per_face> +Mapping::get_vertices( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no) const +{ + boost::container::small_vector, + GeometryInfo::vertices_per_face> + face_vertices; + + const auto &cell_vertices = get_vertices(cell); + const auto &reference_cell = cell->reference_cell(); + const auto face_orientation = cell->combined_face_orientation(face_no); + + for (const unsigned int v : + reference_cell.face_reference_cell(face_no).vertex_indices()) + { + face_vertices.push_back( + cell_vertices[reference_cell.face_to_cell_vertices( + face_no, v, face_orientation)]); + } + + return face_vertices; +} + + + template Point Mapping::get_center( diff --git a/tests/mappings/mapping_get_vertices_on_face.cc b/tests/mappings/mapping_get_vertices_on_face.cc new file mode 100644 index 0000000000..53ef12bbfa --- /dev/null +++ b/tests/mappings/mapping_get_vertices_on_face.cc @@ -0,0 +1,84 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 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. +// +// --------------------------------------------------------------------- + + +// Checks if vertices on face are given in correct oder + +#include + +#include +#include + +#include "../tests.h" + + +template +std::vector> +get_vertices(const typename Triangulation::cell_iterator &cell, + const typename Triangulation::face_iterator &face, + const Mapping & mapping) +{ + std::vector> vertices(face->n_vertices()); + + for (unsigned int i = 0; i < face->n_vertices(); ++i) + { + vertices[i] = mapping.transform_unit_to_real_cell( + cell, mapping.transform_real_to_unit_cell(cell, face->vertex(i))); + } + return vertices; +} + +template +void +test_vertex_order() +{ + deallog << "dim=" << dim << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube(tria, -1.0, 1.0); + + MappingQ mapping(1); + + for (const auto &cell : tria.active_cell_iterators()) + for (unsigned int f = 0; f < cell->n_faces(); ++f) + { + auto const &vertices_detour = + get_vertices(cell, cell->face(f), mapping); + auto const &vertices = mapping.get_vertices(cell, f); + + AssertDimension(vertices_detour.size(), vertices.size()); + + for (unsigned int i = 0; i < vertices.size(); ++i) + { + auto const diff = vertices[i] - vertices_detour[i]; + AssertThrow(diff.norm() < 1e-12, ExcMessage("Vertices differ!")); + } + } + + deallog << "OK" << std::endl; +} + +int +main() +{ + initlog(); + deallog << std::setprecision(12); + + test_vertex_order<1>(); + test_vertex_order<2>(); + test_vertex_order<3>(); + + return 0; +} diff --git a/tests/mappings/mapping_get_vertices_on_face.output b/tests/mappings/mapping_get_vertices_on_face.output new file mode 100644 index 0000000000..4afcfebbaf --- /dev/null +++ b/tests/mappings/mapping_get_vertices_on_face.output @@ -0,0 +1,7 @@ + +DEAL::dim=1 +DEAL::OK +DEAL::dim=2 +DEAL::OK +DEAL::dim=3 +DEAL::OK -- 2.39.5