Invalid = static_cast<std::uint8_t>(-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.
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 <int dim>
+ 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<dim>::faces_per_cell);
+ return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
+ }
+ else if (reference_cell == Type::Tri)
+ {
+ AssertIndexRange(face_no, 3);
+ static const std::array<Tensor<1, dim>, 3> table = {
+ {Point<dim>(1, 0),
+ Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
+ Point<dim>(0, -1)}};
+
+ return table[face_no];
+ }
+ else if (reference_cell == Type::Tet)
+ {
+ AssertIndexRange(face_no, 4);
+ static const std::array<std::array<Tensor<1, dim>, 2>, 4> table = {
+ {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
+ {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
+ {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
+ {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
+ +std::pow(1.0 / 3.0, 1.0 / 4.0),
+ 0),
+ Point<dim>(-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<std::array<Tensor<1, dim>, 2>, 5> table = {
+ {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
+ {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
+ {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
+ Point<dim>(0, 0, 1)}},
+ {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
+ {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}}}};
+
+ return table[face_no][i];
+ }
+ else if (reference_cell == Type::Pyramid)
+ {
+ AssertIndexRange(face_no, 5);
+ static const std::array<std::array<Tensor<1, dim>, 2>, 5> table = {
+ {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
+ {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
+ Point<dim>(0, 1, 0)}},
+ {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
+ Point<dim>(0, 1, 0)}},
+ {{Point<dim>(1, 0, 0),
+ Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
+ {{Point<dim>(1, 0, 0),
+ Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
+
+ return table[face_no][i];
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
namespace internal
{
std::vector<Tensor<1, spacedim>>(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<dim>::face_indices())
+ unit_tangentials[i].resize(n_original_q_points);
+ std::fill(unit_tangentials[i].begin(),
+ unit_tangentials[i].end(),
+ ReferenceCell::unit_tangential_vectors<dim>(
+ 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<dim>::unit_tangential_vectors[i][0]);
- if (dim > 2)
- {
- unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
- .resize(n_original_q_points);
- std::fill(
- unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
- .begin(),
- unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
- .end(),
- GeometryInfo<dim>::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<dim>(
+ 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());
- }
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/grid/reference_cell.h>
+
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+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<dim>(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);
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff