From: Peter Munch Date: Sun, 27 Dec 2020 06:34:18 +0000 (+0100) Subject: Introduce ReferenceCell::unit_tangential_vectors() X-Git-Tag: v9.3.0-rc1~708^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6ae859fd2ca38cf7b34aff481ce25f9f4ed1f821;p=dealii.git Introduce ReferenceCell::unit_tangential_vectors() --- diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 2d3d4580b1..2bfa0943b1 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -45,6 +45,31 @@ namespace ReferenceCell Invalid = static_cast(-1) }; + /** + * Return the dimension of the given reference-cell type @p type. + */ + inline unsigned int + get_dimension(const Type &type) + { + switch (type) + { + case Type::Vertex: + return 0; + case Type::Line: + return 1; + case Type::Tri: + case Type::Quad: + return 2; + case Type::Tet: + case Type::Pyramid: + case Type::Wedge: + case Type::Hex: + return 3; + default: + return numbers::invalid_unsigned_int; + } + } + /** * Return the correct simplex reference cell type for the given dimension * @p dim. @@ -374,7 +399,85 @@ namespace ReferenceCell return temp_; } + /* + * Return i-th unit tangential vector of a face of the reference cell. + * The vectors are arranged such that the + * cross product between the two vectors returns the unit normal vector. + */ + template + inline Tensor<1, dim> + unit_tangential_vectors(const Type & reference_cell, + const unsigned int face_no, + const unsigned int i) + { + AssertDimension(dim, get_dimension(reference_cell)); + AssertIndexRange(i, dim - 1); + if (reference_cell == get_hypercube(dim)) + { + AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + return GeometryInfo::unit_tangential_vectors[face_no][i]; + } + else if (reference_cell == Type::Tri) + { + AssertIndexRange(face_no, 3); + static const std::array, 3> table = { + {Point(1, 0), + Point(-std::sqrt(0.5), +std::sqrt(0.5)), + Point(0, -1)}}; + + return table[face_no]; + } + else if (reference_cell == Type::Tet) + { + AssertIndexRange(face_no, 4); + static const std::array, 2>, 4> table = { + {{{Point(0, 1, 0), Point(1, 0, 0)}}, + {{Point(1, 0, 0), Point(0, 0, 1)}}, + {{Point(0, 0, 1), Point(0, 1, 0)}}, + {{Point(-std::pow(1.0 / 3.0, 1.0 / 4.0), + +std::pow(1.0 / 3.0, 1.0 / 4.0), + 0), + Point(-std::pow(1.0 / 3.0, 1.0 / 4.0), + 0, + +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}}; + + return table[face_no][i]; + } + else if (reference_cell == Type::Wedge) + { + AssertIndexRange(face_no, 5); + static const std::array, 2>, 5> table = { + {{{Point(0, 1, 0), Point(1, 0, 0)}}, + {{Point(1, 0, 0), Point(0, 0, 1)}}, + {{Point(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0), + Point(0, 0, 1)}}, + {{Point(0, 0, 1), Point(0, 1, 0)}}, + {{Point(1, 0, 0), Point(0, 0, 1)}}}}; + + return table[face_no][i]; + } + else if (reference_cell == Type::Pyramid) + { + AssertIndexRange(face_no, 5); + static const std::array, 2>, 5> table = { + {{{Point(0, 1, 0), Point(1, 0, 0)}}, + {{Point(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)), + Point(0, 1, 0)}}, + {{Point(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)), + Point(0, 1, 0)}}, + {{Point(1, 0, 0), + Point(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}}, + {{Point(1, 0, 0), + Point(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}}; + + return table[face_no][i]; + } + + Assert(false, ExcNotImplemented()); + + return {}; + } namespace internal { diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 7921ebfd07..c004084ca2 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -151,259 +151,26 @@ MappingFE::InternalData::initialize_face( std::vector>(n_original_q_points)); // Compute tangentials to the unit cell. - if (this->fe.reference_cell_type() == ReferenceCell::get_hypercube(dim)) + const auto reference_cell_type = this->fe.reference_cell_type(); + const auto n_faces = + ReferenceCell::internal::Info::get_cell(reference_cell_type).n_faces(); + + for (unsigned int i = 0; i < n_faces; ++i) { - for (const unsigned int i : GeometryInfo::face_indices()) + unit_tangentials[i].resize(n_original_q_points); + std::fill(unit_tangentials[i].begin(), + unit_tangentials[i].end(), + ReferenceCell::unit_tangential_vectors( + reference_cell_type, i, 0)); + if (dim > 2) { - unit_tangentials[i].resize(n_original_q_points); - std::fill(unit_tangentials[i].begin(), - unit_tangentials[i].end(), - GeometryInfo::unit_tangential_vectors[i][0]); - if (dim > 2) - { - unit_tangentials[GeometryInfo::faces_per_cell + i] - .resize(n_original_q_points); - std::fill( - unit_tangentials[GeometryInfo::faces_per_cell + i] - .begin(), - unit_tangentials[GeometryInfo::faces_per_cell + i] - .end(), - GeometryInfo::unit_tangential_vectors[i][1]); - } + unit_tangentials[n_faces + i].resize(n_original_q_points); + std::fill(unit_tangentials[n_faces + i].begin(), + unit_tangentials[n_faces + i].end(), + ReferenceCell::unit_tangential_vectors( + reference_cell_type, i, 1)); } } - else if (this->fe.reference_cell_type() == ReferenceCell::Type::Tri) - { - Tensor<1, dim> t1; - constexpr int d0 = 0; - constexpr int d1 = 1 % dim; - - t1[d0] = 1; - t1[d1] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[0].emplace_back(t1); - t1[d0] = -std::sqrt(0.5); - t1[d1] = +std::sqrt(0.5); - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[1].emplace_back(t1); - t1[d0] = 0; - t1[d1] = -1; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[2].emplace_back(t1); - } - else if (this->fe.reference_cell_type() == ReferenceCell::Type::Tet) - { - Tensor<1, dim> t1; - constexpr int d0 = 0; - constexpr int d1 = 1 % dim; - constexpr int d2 = 2 % dim; - - t1[d0] = 0; - t1[d1] = 1; - t1[d2] = 0; // face 0 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[0].emplace_back(t1); - - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; // face 0 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[4].emplace_back(t1); - - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; // face 1 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[1].emplace_back(t1); - - t1[d0] = 0; - t1[d1] = 0; - t1[d2] = 1; // face 1 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[5].emplace_back(t1); - - t1[d0] = 0; - t1[d1] = 0; - t1[d2] = 1; // face 2 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[2].emplace_back(t1); - - t1[d0] = 0; - t1[d1] = 1; - t1[d2] = 0; // face 2 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[6].emplace_back(t1); - - t1[d0] = -std::pow(1.0 / 3.0, 1.0 / 4.0); - t1[d1] = +std::pow(1.0 / 3.0, 1.0 / 4.0); - t1[d2] = +0; // face 3 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[3].emplace_back(t1); - - t1[d0] = -std::pow(1.0 / 3.0, 1.0 / 4.0); - t1[d1] = +0; - t1[d2] = +std::pow(1.0 / 3.0, 1.0 / 4.0); // face 3 - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[7].emplace_back(t1); - } - else if (this->fe.reference_cell_type() == ReferenceCell::Type::Wedge) - { - Tensor<1, dim> t1; - constexpr int d0 = 0; - constexpr int d1 = 1 % dim; - constexpr int d2 = 2 % dim; - - // face 0 - t1[d0] = 0; - t1[d1] = 1; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[0].emplace_back(t1); - - // face 0 - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[5].emplace_back(t1); - - // face 1 - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[1].emplace_back(t1); - - // face 1 - t1[d0] = 0; - t1[d1] = 1; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[6].emplace_back(t1); - - // face 2 - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[2].emplace_back(t1); - - // face 2 - t1[d0] = 0; - t1[d1] = 0; - t1[d2] = 1; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[7].emplace_back(t1); - - // face 3 - t1[d0] = -1 / std::sqrt(2.0); - t1[d1] = +1 / std::sqrt(2.0); - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[3].emplace_back(t1); - - // face 3 - t1[d0] = 0; - t1[d1] = 0; - t1[d2] = 1; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[8].emplace_back(t1); - - // face 4 - t1[d0] = +0; - t1[d1] = +0; - t1[d2] = +1; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[4].emplace_back(t1); - - // face 4 - t1[d0] = +0; - t1[d1] = +1; - t1[d2] = +0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[9].emplace_back(t1); - } - else if (this->fe.reference_cell_type() == ReferenceCell::Type::Pyramid) - { - Tensor<1, dim> t1; - constexpr int d0 = 0; - constexpr int d1 = 1 % dim; - constexpr int d2 = 2 % dim; - - // face 0 - t1[d0] = 0; - t1[d1] = 1; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[0].emplace_back(t1); - - // face 0 - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[5].emplace_back(t1); - - // face 1 - t1[d0] = 1.0 / sqrt(2.0); - t1[d1] = 0; - t1[d2] = 1.0 / sqrt(2.0); - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[1].emplace_back(t1); - - // face 1 - t1[d0] = 0; - t1[d1] = 1; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[6].emplace_back(t1); - - // face 2 - t1[d0] = 1.0 / sqrt(2.0); - t1[d1] = 0; - t1[d2] = -1.0 / sqrt(2.0); - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[2].emplace_back(t1); - - // face 2 - t1[d0] = 0; - t1[d1] = 1; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[7].emplace_back(t1); - - // face 3 - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[3].emplace_back(t1); - - // face 3 - t1[d0] = 0; - t1[d1] = 1.0 / sqrt(2.0); - t1[d2] = 1.0 / sqrt(2.0); - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[8].emplace_back(t1); - - // face 4 - t1[d0] = 1; - t1[d1] = 0; - t1[d2] = 0; - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[4].emplace_back(t1); - - // face 4 - t1[d0] = 0; - t1[d1] = +1.0 / sqrt(2.0); - t1[d2] = -1.0 / sqrt(2.0); - for (unsigned int i = 0; i < n_original_q_points; i++) - unit_tangentials[9].emplace_back(t1); - } - else - { - Assert(false, ExcNotImplemented()); - } } } diff --git a/tests/simplex/unit_tangential_vectors_01.cc b/tests/simplex/unit_tangential_vectors_01.cc new file mode 100644 index 0000000000..660cb585f9 --- /dev/null +++ b/tests/simplex/unit_tangential_vectors_01.cc @@ -0,0 +1,55 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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 ReferenceCell::unit_tangential_vectors() for all reference-cell types. + + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test(const ReferenceCell::Type &reference_cell) +{ + for (const auto face_no : + ReferenceCell::internal::Info::get_cell(reference_cell).face_indices()) + for (unsigned int i = 0; i < dim - 1; ++i) + deallog << ReferenceCell::unit_tangential_vectors(reference_cell, + face_no, + i) + << std::endl; + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + + test<2>(ReferenceCell::Type::Tri); + test<2>(ReferenceCell::Type::Quad); + test<3>(ReferenceCell::Type::Tet); + test<3>(ReferenceCell::Type::Pyramid); + test<3>(ReferenceCell::Type::Wedge); + test<3>(ReferenceCell::Type::Hex); +} diff --git a/tests/simplex/unit_tangential_vectors_01.output b/tests/simplex/unit_tangential_vectors_01.output new file mode 100644 index 0000000000..46a477ea04 --- /dev/null +++ b/tests/simplex/unit_tangential_vectors_01.output @@ -0,0 +1,54 @@ + +DEAL::1.00000 0.00000 +DEAL::-0.707107 0.707107 +DEAL::0.00000 -1.00000 +DEAL:: +DEAL::0.00000 -1.00000 +DEAL::0.00000 1.00000 +DEAL::1.00000 0.00000 +DEAL::-1.00000 0.00000 +DEAL:: +DEAL::0.00000 1.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::0.00000 1.00000 0.00000 +DEAL::-0.759836 0.759836 0.00000 +DEAL::-0.759836 0.00000 0.759836 +DEAL:: +DEAL::0.00000 1.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.707107 0.00000 0.707107 +DEAL::0.00000 1.00000 0.00000 +DEAL::0.707107 0.00000 -0.707107 +DEAL::0.00000 1.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.00000 0.707107 0.707107 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.00000 0.707107 -0.707107 +DEAL:: +DEAL::0.00000 1.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::-0.707107 0.707107 0.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::0.00000 1.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL:: +DEAL::0.00000 -1.00000 0.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::0.00000 1.00000 0.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::0.00000 0.00000 -1.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.00000 0.00000 1.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::-1.00000 0.00000 0.00000 +DEAL::0.00000 1.00000 0.00000 +DEAL::1.00000 0.00000 0.00000 +DEAL::0.00000 1.00000 0.00000 +DEAL::