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