Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
vector_relations.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2021 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_vector_relations_h
17#define dealii_vector_relations_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/tensor.h>
22
23#include <cmath>
24#include <limits>
25
27
28
29namespace Physics
30{
34 namespace VectorRelations
35 {
45 template <int spacedim, typename Number>
46 Number
49
77 template <int spacedim, typename Number>
78 Number
81 const Tensor<1, spacedim, Number> &axis);
82 } // namespace VectorRelations
83} // namespace Physics
84
85
86
87#ifndef DOXYGEN
88
89
90
91template <int spacedim, typename Number>
92inline Number
95{
96 const Number a_norm = a.norm();
97 const Number b_norm = b.norm();
98 Assert(a_norm > 1.e-12 * b_norm && a_norm > 1.e-12 * b_norm,
99 ExcMessage("Both vectors need to be non-zero!"));
100
101 Number argument = (a * b) / a_norm / b_norm;
102
103 // std::acos returns nan if argument is out of domain [-1,+1].
104 // if argument slightly overshoots these bounds, set it to the bound.
105 // allow for 8*eps as a tolerance.
106 if ((1. - std::abs(argument)) < 8. * std::numeric_limits<Number>::epsilon())
107 argument = std::copysign(1., argument);
108
109 return std::acos(argument);
110}
111
112
113
114template <int spacedim, typename Number>
115inline Number
118 const Tensor<1, spacedim, Number> &axis)
119{
120 Assert(spacedim == 3,
121 ExcMessage("This function can only be used with spacedim==3!"));
122
123 Assert(std::abs(axis.norm() - 1.) < 1.e-12,
124 ExcMessage("The axial vector is not a unit vector."));
125 Assert(std::abs(axis * a) < 1.e-12 * a.norm() &&
126 std::abs(axis * b) < 1.e-12 * b.norm(),
127 ExcMessage("The vectors are not perpendicular to the axial vector."));
128
129 const Number dot = a * b;
130 const Number det = axis * cross_product_3d(a, b);
131
132 return std::atan2(det, dot);
133}
134
135
136
137#endif // DOXYGEN
138
140
141#endif
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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 > &)