Reference documentation for deal.II version GIT 00e6fe71c9 2023-06-04 19:35:01+00:00
\(\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\}}\)
vector_relations.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 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 
29 namespace 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 
91 template <int spacedim, typename Number>
92 inline 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 
114 template <int spacedim, typename Number>
115 inline 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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression atan2(const Expression &y, const Expression &x)
Expression acos(const Expression &x)
Expression copysign(const Expression &value, const Expression &sign)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
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)