From d47c7c477427d8745d8f4dab3669efdb86a289df Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Mon, 29 Nov 2021 18:49:53 -0700 Subject: [PATCH] Functions for calculating angles between vectors. --- doc/news/changes/minor/20211129Fehling | 4 + include/deal.II/physics/vector_relations.h | 140 +++++++++++++++++++++ source/grid/manifold_lib.cc | 8 +- 3 files changed, 149 insertions(+), 3 deletions(-) create mode 100644 doc/news/changes/minor/20211129Fehling create mode 100644 include/deal.II/physics/vector_relations.h diff --git a/doc/news/changes/minor/20211129Fehling b/doc/news/changes/minor/20211129Fehling new file mode 100644 index 0000000000..fa639e56b0 --- /dev/null +++ b/doc/news/changes/minor/20211129Fehling @@ -0,0 +1,4 @@ +New: Namespace Physics::VectorRelations features functions to compute +angles between (spatial) vectors. +
+(Marc Fehling, 2021/11/29) diff --git a/include/deal.II/physics/vector_relations.h b/include/deal.II/physics/vector_relations.h new file mode 100644 index 0000000000..efc263447c --- /dev/null +++ b/include/deal.II/physics/vector_relations.h @@ -0,0 +1,140 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_vector_relations_h +#define dealii_vector_relations_h + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace Physics +{ + /** + * Functions to compute relations between spatial vectors. + */ + namespace VectorRelations + { + /** + * Calculate the angle $\theta$ between two vectors @p a and @p b. The returned + * angle will be in the range $[0, \pi]$. + * + * This function uses the geometric definition of the scalar product. + * @f[ + * \vec{a} \cdot \vec{b} = \|\vec{a}\| \|\vec{b}\| \cos(\theta) + * @f] + */ + template + Number + angle(const Tensor<1, spacedim, Number> &a, + const Tensor<1, spacedim, Number> &b); + + /** + * Calculate the angle $\theta$ between two vectors @p a and @p b, where both + * vectors are located in a plane described by a normal vector @p axis. + * + * The angle computed by this function corresponds to the rotation angle + * that would transform the vector @p a into the vector @p b around the vector + * @p axis. Thus, contrary to the function above, we get a @em signed angle + * which will be in the range $[-\pi, \pi]$. + * + * The vector @p axis needs to be a unit vector and be perpendicular to both + * vectors @p a and @p b. + * + * This function uses the geometric definitions of both the scalar and cross + * product. + * @f{align*}{ + * \vec{a} \cdot \vec{b} &= \|\vec{a}\| \|\vec{b}\| \cos(\theta) \\ + * \vec{a} \times \vec{b} &= \|\vec{a}\| \|\vec{b}\| \sin(\theta) \vec{n} + * @f} + * We can create the tangent of the angle using both products. + * @f[ + * \tan{\theta} + * = \frac{\sin(\theta)}{\cos(theta)} + * = \frac{(\vec{a} \times \vec{b}) \cdot \vec{n}}{\vec{a} \cdot \vec{b}} + * @f] + * + * @note Only applicable for three-dimensional vectors `spacedim == 3`. + */ + template + Number + signed_angle(const Tensor<1, spacedim, Number> &a, + const Tensor<1, spacedim, Number> &b, + const Tensor<1, spacedim, Number> &axis); + } // namespace VectorRelations +} // namespace Physics + + + +#ifndef DOXYGEN + + + +template +inline Number +Physics::VectorRelations::angle(const Tensor<1, spacedim, Number> &a, + const Tensor<1, spacedim, Number> &b) +{ + const Number a_norm = a.norm(); + const Number b_norm = b.norm(); + Assert(a_norm > 1.e-12 * b_norm && a_norm > 1.e-12 * b_norm, + ExcMessage("Both vectors need to be non-zero!")); + + Number argument = (a * b) / a_norm / b_norm; + + // std::acos returns nan if argument is out of domain [-1,+1]. + // if argument slightly overshoots these bounds, set it to the bound. + // allow for 8*eps as a tolerance. + if ((1. - std::abs(argument)) < 8. * std::numeric_limits::epsilon()) + argument = std::copysign(1., argument); + + return std::acos(argument); +} + + + +template +inline Number +Physics::VectorRelations::signed_angle(const Tensor<1, spacedim, Number> &a, + const Tensor<1, spacedim, Number> &b, + const Tensor<1, spacedim, Number> &axis) +{ + Assert(spacedim == 3, + ExcMessage("This function can only be used with spacedim==3!")); + + Assert(std::abs(axis.norm() - 1.) < 1.e-12, + ExcMessage("The axial vector is not a unit vector.")); + Assert(std::abs(axis * a) < 1.e-12 * b.norm() && + std::abs(axis * b) < 1.e-12 * a.norm(), + ExcMessage("The vectors are not perpendicular to the axial vector.")); + + const Number dot = a * b; + const Number det = axis * cross_product_3d(a, b); + + return std::atan2(det, dot); +} + + + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 052718befa..bc78bc58b5 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -26,6 +26,8 @@ #include +#include + DEAL_II_DISABLE_EXTRA_DIAGNOSTICS #include DEAL_II_ENABLE_EXTRA_DIAGNOSTICS @@ -1129,9 +1131,9 @@ CylindricalManifold::pull_back( // Then compute the angle between the projection direction and // another vector orthogonal to the direction vector. - const double dot = normal_direction * p_diff; - const double det = direction * cross_product_3d(normal_direction, p_diff); - const double phi = std::atan2(det, dot); + const double phi = Physics::VectorRelations::signed_angle(normal_direction, + p_diff, + /*axis=*/direction); // Return distance from the axis, angle and signed distance on the axis. return Point<3>(p_diff.norm(), phi, lambda); -- 2.39.5