Reference documentation for deal.II version Git 418c2197d0 2021-10-27 08:19:48 +0200
\(\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\}}\)
function_spherical.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2018 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 
24 namespace 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  {
322  AssertThrow(false, ExcNotImplemented());
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  {
333  AssertThrow(false, ExcNotImplemented());
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  {
344  AssertThrow(false, ExcNotImplemented());
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 double value(const Point< dim > &point, const unsigned int component=0) const override
virtual std::array< double, 6 > shessian(const std::array< double, dim > &sp, const unsigned int component) const
Spherical(const Point< dim > &center=Point< dim >(), const unsigned int n_components=1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
static ::ExceptionBase & ExcDivideByZero()
virtual std::array< double, dim > sgradient(const std::array< double, dim > &sp, const unsigned int component) const
virtual SymmetricTensor< 2, dim > hessian(const Point< dim > &p, const unsigned int component=0) const override
virtual double svalue(const std::array< double, dim > &sp, const unsigned int component) const
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
std::array< double, dim > to_spherical(const Point< dim > &point)
#define Assert(cond, exc)
Definition: exceptions.h:1461
const Tensor< 1, dim > coordinate_system_offset
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
virtual std::size_t memory_consumption() const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
static ::ExceptionBase & ExcNotImplemented()