-
Notifications
You must be signed in to change notification settings - Fork 707
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Functions for calculating angles between vectors.
- Loading branch information
1 parent
9a23ab9
commit b802714
Showing
2 changed files
with
137 additions
and
3 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,133 @@ | ||
// --------------------------------------------------------------------- | ||
// | ||
// 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_relations_h | ||
#define dealii_relations_h | ||
|
||
#include <deal.II/base/config.h> | ||
|
||
#include <deal.II/base/tensor.h> | ||
|
||
#include <cmath> | ||
|
||
DEAL_II_NAMESPACE_OPEN | ||
|
||
|
||
namespace Physics | ||
{ | ||
/** | ||
* Functions to determine relations between vectors. | ||
*/ | ||
namespace Relations | ||
{ | ||
/** | ||
* 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 <int spacedim, typename Number> | ||
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 describes the relation with respect to the vector @p axis and will | ||
* be in the range $[-\pi, \pi]$. | ||
* | ||
* The vector @p axis needs to be a unit vector and 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 <int spacedim, typename Number> | ||
Number | ||
angle(const Tensor<1, spacedim, Number> &a, | ||
const Tensor<1, spacedim, Number> &b, | ||
const Tensor<1, spacedim, Number> &axis); | ||
} // namespace Relations | ||
} // namespace Physics | ||
|
||
|
||
|
||
#ifndef DOXYGEN | ||
|
||
|
||
|
||
template <int spacedim, typename Number> | ||
inline Number | ||
Physics::Relations::angle(const Tensor<1, spacedim, Number> &a, | ||
const Tensor<1, spacedim, Number> &b) | ||
{ | ||
Assert(a.norm() > 1.e-12 && b.norm() > 1.e-12, | ||
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<Number>::epsilon()) | ||
argument = std::copysign(1., argument); | ||
|
||
return std::acos(argument); | ||
} | ||
|
||
|
||
|
||
template <int spacedim, typename Number> | ||
inline Number | ||
Physics::Relations::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 && std::abs(axis * b) < 1.e-12, | ||
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters