Reference documentation for deal.II version 9.6.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
function_spherical.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 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
17#include <deal.II/base/point.h>
18
19#include <algorithm>
20#include <cmath>
21
23namespace Functions
24{
25 // other implementations/notes:
26 // https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/geometry/euclidean/threed/SphericalCoordinates.java
27 // http://mathworld.wolfram.com/SphericalCoordinates.html
28
29 /*derivation of Hessian in Maxima as function of tensor products of unit
30 vectors:
31
32 depends(ur,[theta,phi]);
33 depends(utheta,theta);
34 depends(uphi,[theta,phi]);
35 depends(f,[r,theta,phi]);
36 declare([f,r,theta,phi], scalar)@f$
37 dotscrules: true;
38 grads(a):=ur.diff(a,r)+(1/r)*uphi.diff(a,phi)+(1/(r*sin(phi)))*utheta.diff(a,theta);
39
40
41 H : factor(grads(grads(f)));
42 H2: subst([diff(ur,theta)=sin(phi)*utheta,
43 diff(utheta,theta)=-cos(phi)*uphi-sin(phi)*ur,
44 diff(uphi,theta)=cos(phi)*utheta,
45 diff(ur,phi)=uphi,
46 diff(uphi,phi)=-ur],
47 H);
48 H3: trigsimp(fullratsimp(H2));
49
50
51 srules : [diff(f,r)=sg0,
52 diff(f,theta)=sg1,
53 diff(f,phi)=sg2,
54 diff(f,r,2)=sh0,
55 diff(f,theta,2)=sh1,
56 diff(f,phi,2)=sh2,
57 diff(f,r,1,theta,1)=sh3,
58 diff(f,r,1,phi,1)=sh4,
59 diff(f,theta,1,phi,1)=sh5,
60 cos(phi)=cos_phi,
61 cos(theta)=cos_theta,
62 sin(phi)=sin_phi,
63 sin(theta)=sin_theta
64 ]@f$
65
66 c_utheta2 : distrib(subst(srules, ratcoeff(expand(H3), utheta.utheta)));
67 c_utheta_ur : (subst(srules, ratcoeff(expand(H3), utheta.ur)));
68 (subst(srules, ratcoeff(expand(H3), ur.utheta))) - c_utheta_ur;
69 c_utheta_uphi : (subst(srules, ratcoeff(expand(H3), utheta.uphi)));
70 (subst(srules, ratcoeff(expand(H3), uphi.utheta))) - c_utheta_uphi;
71 c_ur2 : (subst(srules, ratcoeff(expand(H3), ur.ur)));
72 c_ur_uphi : (subst(srules, ratcoeff(expand(H3), ur.uphi)));
73 (subst(srules, ratcoeff(expand(H3), uphi.ur))) - c_ur_uphi;
74 c_uphi2 : (subst(srules, ratcoeff(expand(H3), uphi.uphi)));
75
76
77 where (used later to do tensor products):
78
79 ur : [cos(theta)*sin(phi), sin(theta)*sin(phi), cos(phi)];
80 utheta : [-sin(theta), cos(theta), 0];
81 uphi : [cos(theta)*cos(phi), sin(theta)*cos(phi), -sin(phi)];
82
83 with the following proof of substitution rules above:
84
85 -diff(ur,theta)+sin(phi)*utheta;
86 trigsimp(-diff(utheta,theta)-cos(phi)*uphi-sin(phi)*ur);
87 -diff(uphi,theta)+cos(phi)*utheta;
88 -diff(ur,phi)+uphi;
89 -diff(uphi,phi)-ur;
90
91 */
92
93 namespace
94 {
98 template <int dim>
99 void
100 set_unit_vectors(const double cos_theta,
101 const double sin_theta,
102 const double cos_phi,
103 const double sin_phi,
104 Tensor<1, dim> &unit_r,
105 Tensor<1, dim> &unit_theta,
106 Tensor<1, dim> &unit_phi)
107 {
108 unit_r[0] = cos_theta * sin_phi;
109 unit_r[1] = sin_theta * sin_phi;
110 unit_r[2] = cos_phi;
111
112 unit_theta[0] = -sin_theta;
113 unit_theta[1] = cos_theta;
114 unit_theta[2] = 0.;
115
116 unit_phi[0] = cos_theta * cos_phi;
117 unit_phi[1] = sin_theta * cos_phi;
118 unit_phi[2] = -sin_phi;
119 }
120
121
125 template <int dim>
126 void
127 add_outer_product(SymmetricTensor<2, dim> &out,
128 const double val,
129 const Tensor<1, dim> &in1,
130 const Tensor<1, dim> &in2)
131 {
132 if (val != 0.)
133 for (unsigned int i = 0; i < dim; ++i)
134 for (unsigned int j = i; j < dim; ++j)
135 out[i][j] += (in1[i] * in2[j] + in1[j] * in2[i]) * val;
136 }
137
141 template <int dim>
142 void
143 add_outer_product(SymmetricTensor<2, dim> &out,
144 const double val,
145 const Tensor<1, dim> &in)
146 {
147 if (val != 0.)
148 for (unsigned int i = 0; i < dim; ++i)
149 for (unsigned int j = i; j < dim; ++j)
150 out[i][j] += val * in[i] * in[j];
151 }
152 } // namespace
153
154
155
156 template <int dim>
158 const unsigned int n_components)
159 : Function<dim>(n_components)
160 , coordinate_system_offset(p)
161 {
162 AssertThrow(dim == 3, ExcNotImplemented());
163 }
164
165
166
167 template <int dim>
168 double
170 const unsigned int component) const
171 {
172 const Point<dim> p = p_ - coordinate_system_offset;
173 const std::array<double, dim> sp =
175 return svalue(sp, component);
176 }
177
178
179
180 template <int dim>
183 const unsigned int /*component*/) const
184
185 {
187 return {};
188 }
189
190
191
192 template <>
194 Spherical<3>::gradient(const Point<3> &p_, const unsigned int component) const
195 {
196 constexpr int dim = 3;
197 const Point<dim> p = p_ - coordinate_system_offset;
198 const std::array<double, dim> sp =
200 const std::array<double, dim> sg = sgradient(sp, component);
201
202 // somewhat backwards, but we need cos/sin's for unit vectors
203 const double cos_theta = std::cos(sp[1]);
204 const double sin_theta = std::sin(sp[1]);
205 const double cos_phi = std::cos(sp[2]);
206 const double sin_phi = std::sin(sp[2]);
207
208 Tensor<1, dim> unit_r, unit_theta, unit_phi;
209 set_unit_vectors(
210 cos_theta, sin_theta, cos_phi, sin_phi, unit_r, unit_theta, unit_phi);
211
212 Tensor<1, dim> res;
213
214 if (sg[0] != 0.)
215 {
216 res += unit_r * sg[0];
217 }
218
219 if (sg[1] * sin_phi != 0.)
220 {
221 Assert(sp[0] != 0., ExcDivideByZero());
222 res += unit_theta * sg[1] / (sp[0] * sin_phi);
223 }
224
225 if (sg[2] != 0.)
226 {
227 Assert(sp[0] != 0., ExcDivideByZero());
228 res += unit_phi * sg[2] / sp[0];
229 }
230
231 return res;
232 }
233
234
235
236 template <int dim>
239 const unsigned int /*component*/) const
240 {
242 return {};
243 }
244
245
246
247 template <>
249 Spherical<3>::hessian(const Point<3> &p_, const unsigned int component) const
250
251 {
252 constexpr int dim = 3;
253 const Point<dim> p = p_ - coordinate_system_offset;
254 const std::array<double, dim> sp =
256 const std::array<double, dim> sg = sgradient(sp, component);
257 const std::array<double, 6> sh = shessian(sp, component);
258
259 // somewhat backwards, but we need cos/sin's for unit vectors
260 const double cos_theta = std::cos(sp[1]);
261 const double sin_theta = std::sin(sp[1]);
262 const double cos_phi = std::cos(sp[2]);
263 const double sin_phi = std::sin(sp[2]);
264 const double r = sp[0];
265
266 Tensor<1, dim> unit_r, unit_theta, unit_phi;
267 set_unit_vectors(
268 cos_theta, sin_theta, cos_phi, sin_phi, unit_r, unit_theta, unit_phi);
269
270 const double sin_phi2 = sin_phi * sin_phi;
271 const double r2 = r * r;
272 Assert(r != 0., ExcDivideByZero());
273
274 const double c_utheta2 =
275 sg[0] / r + ((sin_phi != 0.) ? (cos_phi * sg[2]) / (r2 * sin_phi) +
276 sh[1] / (r2 * sin_phi2) :
277 0.);
278 const double c_utheta_ur =
279 ((sin_phi != 0.) ? (r * sh[3] - sg[1]) / (r2 * sin_phi) : 0.);
280 const double c_utheta_uphi =
281 ((sin_phi != 0.) ? (sh[5] * sin_phi - cos_phi * sg[1]) / (r2 * sin_phi2) :
282 0.);
283 const double c_ur2 = sh[0];
284 const double c_ur_uphi = (r * sh[4] - sg[2]) / r2;
285 const double c_uphi2 = (sh[2] + r * sg[0]) / r2;
286
287 // go through each tensor product
289
290 add_outer_product(res, c_utheta2, unit_theta);
291
292 add_outer_product(res, c_utheta_ur, unit_theta, unit_r);
293
294 add_outer_product(res, c_utheta_uphi, unit_theta, unit_phi);
295
296 add_outer_product(res, c_ur2, unit_r);
297
298 add_outer_product(res, c_ur_uphi, unit_r, unit_phi);
299
300 add_outer_product(res, c_uphi2, unit_phi);
301
302 return res;
303 }
304
305
306
307 template <int dim>
308 std::size_t
310 {
311 return sizeof(Spherical<dim>);
312 }
313
314
315
316 template <int dim>
317 double
318 Spherical<dim>::svalue(const std::array<double, dim> & /* sp */,
319 const unsigned int /*component*/) const
320 {
322 return 0.;
323 }
324
325
326
327 template <int dim>
328 std::array<double, dim>
329 Spherical<dim>::sgradient(const std::array<double, dim> & /* sp */,
330 const unsigned int /*component*/) const
331 {
333 return std::array<double, dim>();
334 }
335
336
337
338 template <int dim>
339 std::array<double, 6>
340 Spherical<dim>::shessian(const std::array<double, dim> & /* sp */,
341 const unsigned int /*component*/) const
342 {
344 return std::array<double, 6>();
345 }
346
347
348
349 // explicit instantiations
350 template class Spherical<1>;
351 template class Spherical<2>;
352 template class Spherical<3>;
353
354} // namespace Functions
355
virtual SymmetricTensor< 2, dim > hessian(const Point< dim > &p, const unsigned int component=0) const override
Spherical(const Point< dim > &center=Point< dim >(), const unsigned int n_components=1)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
virtual std::array< double, dim > sgradient(const std::array< double, dim > &sp, const unsigned int component) const
virtual double value(const Point< dim > &point, const unsigned int component=0) const override
virtual std::size_t memory_consumption() const override
virtual std::array< double, 6 > shessian(const std::array< double, dim > &sp, const unsigned int component) const
virtual double svalue(const std::array< double, dim > &sp, const unsigned int component) const
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDivideByZero()
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
std::array< double, dim > to_spherical(const Point< dim > &point)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)