Reference documentation for deal.II version GIT relicensing-756-gc37f912904 2024-05-25 08:20:01+00:00
Searching...
No Matches
vector_relations.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
4// Copyright (C) 2021 - 2022 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_vector_relations_h
16#define dealii_vector_relations_h
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/base/tensor.h>
21
22#include <cmath>
23#include <limits>
24
26
27
28namespace Physics
29{
33 namespace VectorRelations
34 {
44 template <int spacedim, typename Number>
45 Number
48
76 template <int spacedim, typename Number>
77 Number
80 const Tensor<1, spacedim, Number> &axis);
81 } // namespace VectorRelations
82} // namespace Physics
83
84
85
86#ifndef DOXYGEN
87
88
89
90template <int spacedim, typename Number>
91inline Number
94{
95 const Number a_norm = a.norm();
96 const Number b_norm = b.norm();
97 Assert(a_norm > 1.e-12 * b_norm && a_norm > 1.e-12 * b_norm,
98 ExcMessage("Both vectors need to be non-zero!"));
99
100 Number argument = (a * b) / a_norm / b_norm;
101
102 // std::acos returns nan if argument is out of domain [-1,+1].
103 // if argument slightly overshoots these bounds, set it to the bound.
104 // allow for 8*eps as a tolerance.
105 if ((1. - std::abs(argument)) < 8. * std::numeric_limits<Number>::epsilon())
106 argument = std::copysign(1., argument);
107
108 return std::acos(argument);
109}
110
111
112
113template <int spacedim, typename Number>
114inline Number
117 const Tensor<1, spacedim, Number> &axis)
118{
119 Assert(spacedim == 3,
120 ExcMessage("This function can only be used with spacedim==3!"));
121
122 Assert(std::abs(axis.norm() - 1.) < 1.e-12,
123 ExcMessage("The axial vector is not a unit vector."));
124 Assert(std::abs(axis * a) < 1.e-12 * a.norm() &&
125 std::abs(axis * b) < 1.e-12 * b.norm(),
126 ExcMessage("The vectors are not perpendicular to the axial vector."));
127
128 const Number dot = a * b;
129 const Number det = axis * cross_product_3d(a, b);
130
131 return std::atan2(det, dot);
132}
133
134
135
136#endif // DOXYGEN
137
139
140#endif
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number signed_angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b, const Tensor< 1, spacedim, Number > &axis)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > acos(const ::VectorizedArray< Number, width > &x)