Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50:00+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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
geometric_utilities.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2024 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
19#include <deal.II/base/point.h>
20
21#include <array>
22#include <cmath>
23#include <limits>
24
25
27
28
29namespace GeometricUtilities
30{
31 namespace Coordinates
32 {
34 double,
35 << "The radius <" << arg1 << "> can not be negative.");
36
38 double,
39 << "The azimuth angle <" << arg1 << "> is not in [0,2Pi).");
40
42 double,
43 << "The polar angle <" << arg1 << "> is not in [0,Pi].");
44
45
46 template <int dim>
47 std::array<double, dim>
48 to_spherical(const Point<dim> &position)
49 {
50 std::array<double, dim> scoord{};
51 Assert(dim > 1, ExcNotImplemented());
52
53 // radius
54 if constexpr (dim > 1)
55 {
56 scoord[0] = position.norm();
57 // azimuth angle \theta:
58 scoord[1] = std::atan2(position[1], position[0]);
59 // correct to [0,2*pi)
60 if (scoord[1] < 0.0)
61 scoord[1] += 2.0 * numbers::PI;
62 }
63
64 // polar angle \phi:
65 if constexpr (dim == 3)
66 {
67 // acos returns the angle in the range [0,\pi]
68 if (scoord[0] > std::numeric_limits<double>::min())
69 scoord[2] = std::acos(position[2] / scoord[0]);
70 else
71 scoord[2] = 0.0;
72 }
73 return scoord;
74 }
75
76 template <std::size_t dim>
78 from_spherical(const std::array<double, dim> &scoord)
79 {
80 Point<dim> ccoord;
81
82 Assert(scoord[0] >= 0., NegativeRadius(scoord[0]));
83
84 Assert(scoord[1] >= 0. && scoord[1] < 2. * numbers::PI,
85 SphericalAzimuth(scoord[1]));
86
87 switch (dim)
88 {
89 case 2:
90 {
91 ccoord[0] = scoord[0] * std::cos(scoord[1]);
92 ccoord[1] = scoord[0] * std::sin(scoord[1]);
93 break;
94 }
95 case 3:
96 {
97 Assert(scoord[2] >= 0. && scoord[2] <= numbers::PI,
98 SphericalPolar(scoord[2]));
99
100 ccoord[0] = scoord[0] * std::sin(scoord[2]) * std::cos(scoord[1]);
101 ccoord[1] = scoord[0] * std::sin(scoord[2]) * std::sin(scoord[1]);
102 ccoord[2] = scoord[0] * std::cos(scoord[2]);
103 break;
104 }
105 default:
107 break;
108 }
109
110 return ccoord;
111 }
112
113 // explicit instantiations
114#include "base/geometric_utilities.inst"
115
116 } // namespace Coordinates
117} // namespace GeometricUtilities
118
Definition point.h:113
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & SphericalAzimuth(double arg1)
static ::ExceptionBase & NegativeRadius(double arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & SphericalPolar(double arg1)
#define DeclException1(Exception1, type1, outsequence)
Point< dim > from_spherical(const std::array< double, dim > &scoord)
std::array< double, dim > to_spherical(const Point< dim > &point)
constexpr double PI
Definition numbers.h:240
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > acos(const ::VectorizedArray< Number, width > &x)