]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add spherical function 3856/head
authorDenis Davydov <davydden@gmail.com>
Mon, 30 Jan 2017 17:22:55 +0000 (18:22 +0100)
committerDenis Davydov <davydden@gmail.com>
Mon, 30 Jan 2017 22:04:00 +0000 (23:04 +0100)
doc/news/changes/minor/20170128DenisDavydov [new file with mode: 0644]
include/deal.II/base/function_spherical.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/function_spherical.cc [new file with mode: 0644]
tests/base/functions_spherical_01.cc [new file with mode: 0644]
tests/base/functions_spherical_01.output [new file with mode: 0644]
tests/base/functions_spherical_02.cc [new file with mode: 0644]
tests/base/functions_spherical_02.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170128DenisDavydov b/doc/news/changes/minor/20170128DenisDavydov
new file mode 100644 (file)
index 0000000..129c8d1
--- /dev/null
@@ -0,0 +1,4 @@
+New: Implement Functions::Spherical to represent a function given in spherical
+coordinates.
+<br>
+(Denis Davydov, 2017/01/28)
diff --git a/include/deal.II/base/function_spherical.h b/include/deal.II/base/function_spherical.h
new file mode 100644 (file)
index 0000000..1ed7726
--- /dev/null
@@ -0,0 +1,118 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__function_spherical_h
+#define dealii__function_spherical_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/std_cxx11/array.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions
+{
+
+  /**
+   * An abstract base class for a scalar-valued function $f=f(r,\theta,\phi)$
+   * defined in spherical coordinates. This class wraps transformation of values,
+   * gradients and hessians from spherical coordinates to the Cartesian coordinate
+   * system used by the Function base class.
+   * Therefore derived classes only need to implement those functions in
+   * spherical coordinates (specifically svalue(), sgradient() and
+   * shessian() ). The convention for angles is the same as in
+   * GeometricUtilities::Coordinates.
+   *
+   * @note This function is currently only implemented for dim==3 .
+   *
+   * @author Denis Davydov, 2017
+   */
+  template <int dim>
+  class Spherical : public Function<dim>
+  {
+  public:
+    /**
+     * Constructor which should be provided with @p center defining the origin
+     * of the coordinate system.
+     */
+    Spherical(const Point<dim> &center = Point<dim>());
+
+    /**
+     * Return the value of the function at the given point.
+     *
+     * This function converts the given point to spherical coordinates,
+     * calls svalue() with it, and returns the result.
+     */
+    virtual double value (const Point<dim> &point,
+                          const unsigned int component = 0) const;
+
+    /**
+     * Return the gradient with respect to the Cartesian coordinates at point @p p.
+     *
+     * This function converts the given point to spherical coordinates,
+     * calls sgradient() with it, and converts the result into Cartesian
+     * coordinates.
+     */
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                    const unsigned int  component = 0) const;
+
+    /**
+     * Return the Hessian with respect to the Cartesian coordinates at point @p p.
+     *
+     * This function converts the given point to spherical coordinates,
+     * calls sgradient and shessian() with it, and converts the result into
+     * Cartesian coordinates.
+     */
+    virtual SymmetricTensor<2,dim> hessian (const Point<dim> &p,
+                                            const unsigned int component=0) const;
+
+    std::size_t memory_consumption () const;
+
+  private:
+    /**
+     * Return the value at point @p sp. Here, @p sp is provided in spherical
+     * coordinates.
+     */
+    virtual double svalue(const std_cxx11::array<double, dim> &sp) const;
+
+    /**
+     * Return the gradient in spherical coordinates.
+     *
+     * The returned object should contain derivatives in the following order:
+     * $\{ f_{,r},\, f_{,\theta},\, f_{,\phi}\}$.
+     */
+    virtual std_cxx11::array<double, dim> sgradient(const std_cxx11::array<double, dim> &sp) const;
+
+    /**
+     * Return the Hessian in spherical coordinates.
+     *
+     * The returned object should contain derivatives in the following order:
+     * $\{ f_{,rr},\, f_{,\theta\theta},\, f_{,\phi\phi},\, f_{,r\theta},\, f_{,r\phi},\, f_{,\theta\phi}\}$.
+     */
+    virtual std_cxx11::array<double, 6> shessian (const std_cxx11::array<double, dim> &sp) const;
+
+    /**
+     * A vector from the origin to the center of spherical coordinate system.
+     */
+    const Tensor<1,dim> coordinate_system_offset;
+  };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
index 66ece18044b9a11e15be407da0551f7bb34ff522..0bff1fa5833d2860dca8619ba13c6fd6143bc0c8 100644 (file)
@@ -30,6 +30,7 @@ SET(_src
   function_lib.cc
   function_lib_cutoff.cc
   function_parser.cc
+  function_spherical.cc
   function_time.cc
   geometry_info.cc
   geometric_utilities.cc
diff --git a/source/base/function_spherical.cc b/source/base/function_spherical.cc
new file mode 100644 (file)
index 0000000..56b8117
--- /dev/null
@@ -0,0 +1,348 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/function_spherical.h>
+#include <deal.II/base/geometric_utilities.h>
+
+#include <cmath>
+#include <algorithm>
+
+DEAL_II_NAMESPACE_OPEN
+namespace Functions
+{
+
+  // other implementations/notes:
+  // https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/geometry/euclidean/threed/SphericalCoordinates.java
+  // http://mathworld.wolfram.com/SphericalCoordinates.html
+
+  /*derivation of Hessian in Maxima as function of tensor products of unit vectors:
+
+  depends(ur,[theta,phi]);
+  depends(utheta,theta);
+  depends(uphi,[theta,phi]);
+  depends(f,[r,theta,phi]);
+  declare([f,r,theta,phi], scalar)$
+  dotscrules: true;
+  grads(a):=ur.diff(a,r)+(1/r)*uphi.diff(a,phi)+(1/(r*sin(phi)))*utheta.diff(a,theta);
+
+
+  H : factor(grads(grads(f)));
+  H2: subst([diff(ur,theta)=sin(phi)*utheta,
+       diff(utheta,theta)=-cos(phi)*uphi-sin(phi)*ur,
+       diff(uphi,theta)=cos(phi)*utheta,
+       diff(ur,phi)=uphi,
+       diff(uphi,phi)=-ur],
+      H);
+  H3: trigsimp(fullratsimp(H2));
+
+
+  srules : [diff(f,r)=sg0,
+          diff(f,theta)=sg1,
+          diff(f,phi)=sg2,
+          diff(f,r,2)=sh0,
+          diff(f,theta,2)=sh1,
+          diff(f,phi,2)=sh2,
+          diff(f,r,1,theta,1)=sh3,
+          diff(f,r,1,phi,1)=sh4,
+          diff(f,theta,1,phi,1)=sh5,
+          cos(phi)=cos_phi,
+          cos(theta)=cos_theta,
+          sin(phi)=sin_phi,
+          sin(theta)=sin_theta
+        ]$
+
+  c_utheta2 : distrib(subst(srules, ratcoeff(expand(H3), utheta.utheta)));
+  c_utheta_ur : (subst(srules, ratcoeff(expand(H3), utheta.ur)));
+  (subst(srules, ratcoeff(expand(H3), ur.utheta))) - c_utheta_ur;
+  c_utheta_uphi : (subst(srules, ratcoeff(expand(H3), utheta.uphi)));
+  (subst(srules, ratcoeff(expand(H3), uphi.utheta))) - c_utheta_uphi;
+  c_ur2 : (subst(srules, ratcoeff(expand(H3), ur.ur)));
+  c_ur_uphi : (subst(srules, ratcoeff(expand(H3), ur.uphi)));
+  (subst(srules, ratcoeff(expand(H3), uphi.ur))) - c_ur_uphi;
+  c_uphi2 : (subst(srules, ratcoeff(expand(H3), uphi.uphi)));
+
+
+  where (used later to do tensor products):
+
+  ur     : [cos(theta)*sin(phi), sin(theta)*sin(phi), cos(phi)];
+  utheta : [-sin(theta), cos(theta), 0];
+  uphi   : [cos(theta)*cos(phi), sin(theta)*cos(phi), -sin(phi)];
+
+  with the following proof of substitution rules above:
+
+  -diff(ur,theta)+sin(phi)*utheta;
+  trigsimp(-diff(utheta,theta)-cos(phi)*uphi-sin(phi)*ur);
+  -diff(uphi,theta)+cos(phi)*utheta;
+  -diff(ur,phi)+uphi;
+  -diff(uphi,phi)-ur;
+
+   */
+
+  namespace
+  {
+
+    /**
+     * Evaluate unit vectors in spherical coordinates
+     */
+    template <int dim>
+    void set_unit_vectors(const double &cos_theta,
+                          const double &sin_theta,
+                          const double &cos_phi,
+                          const double &sin_phi,
+                          Tensor<1,dim> &unit_r,
+                          Tensor<1,dim> &unit_theta,
+                          Tensor<1,dim> &unit_phi)
+    {
+      unit_r[0]     = cos_theta * sin_phi;
+      unit_r[1]     = sin_theta * sin_phi;
+      unit_r[2]     = cos_phi;
+
+      unit_theta[0] = -sin_theta;
+      unit_theta[1] =  cos_theta;
+      unit_theta[2] = 0.;
+
+      unit_phi[0]   = cos_theta * cos_phi;
+      unit_phi[1]   = sin_theta * cos_phi;
+      unit_phi[2]   = -sin_phi;
+    }
+
+
+    /**
+     * calculates out[i][j] += v*(in1[i]*in2[j]+in1[j]*in2[i])
+     */
+    template <int dim>
+    void add_outer_product(SymmetricTensor<2,dim> &out,
+                           const double &val,
+                           const Tensor<1,dim> &in1,
+                           const Tensor<1,dim> &in2)
+    {
+      if (val != 0.)
+        for (unsigned int i = 0; i < dim; i++)
+          for (unsigned int j = i; j < dim; j++)
+            out[i][j] += (in1[i]*in2[j]+in1[j]*in2[i])*val;
+    }
+
+    /**
+     * calculates out[i][j] += v*in[i]in[j]
+     */
+    template <int dim>
+    void add_outer_product(SymmetricTensor<2,dim> &out,
+                           const double &val,
+                           const Tensor<1,dim> &in)
+    {
+      if (val != 0.)
+        for (unsigned int i = 0; i < dim; i++)
+          for (unsigned int j = i; j < dim; j++)
+            out[i][j] += val*in[i]*in[j];
+    }
+  }
+
+
+
+  template <int dim>
+  Spherical<dim>::Spherical(const Point<dim> &p)
+    :
+    Function<dim>(1),
+    coordinate_system_offset(p)
+  {
+    AssertThrow(dim==3,
+                ExcNotImplemented());
+  }
+
+
+
+  template <int dim>
+  double
+  Spherical<dim>::value (const Point<dim>   &p_,
+                         const unsigned int) const
+  {
+    const Point<dim> p = p_ - coordinate_system_offset;
+    const std_cxx11::array<double, dim> sp = GeometricUtilities::Coordinates::to_spherical(p);
+    return svalue(sp);
+  }
+
+
+
+  template <int dim>
+  Tensor<1,dim>
+  Spherical<dim>::gradient (const Point<dim>   &p_,
+                            const unsigned int) const
+  {
+    const Point<dim> p = p_ - coordinate_system_offset;
+    const std_cxx11::array<double, dim> sp = GeometricUtilities::Coordinates::to_spherical(p);
+    const std_cxx11::array<double, dim> sg = sgradient(sp);
+
+    // somewhat backwards, but we need cos/sin's for unit vectors
+    const double cos_theta = std::cos(sp[1]);
+    const double sin_theta = std::sin(sp[1]);
+    const double cos_phi   = std::cos(sp[2]);
+    const double sin_phi   = std::sin(sp[2]);
+
+    Tensor<1,dim> unit_r, unit_theta, unit_phi;
+    set_unit_vectors(cos_theta, sin_theta,
+                     cos_phi, sin_phi,
+                     unit_r, unit_theta, unit_phi);
+
+    Tensor<1,dim> res;
+
+    if (sg[0] != 0.)
+      {
+        res += unit_r * sg[0];
+      }
+
+    if (sg[1] * sin_phi != 0.)
+      {
+        Assert (sp[0] != 0.,
+                ExcDivideByZero());
+        res += unit_theta * sg[1] / (sp[0] * sin_phi);
+      }
+
+    if (sg[2] != 0.)
+      {
+        Assert (sp[0] != 0.,
+                ExcDivideByZero());
+        res += unit_phi * sg[2] / sp[0];
+      }
+
+    return res;
+  }
+
+
+
+  template <int dim>
+  SymmetricTensor<2,dim>
+  Spherical<dim>::hessian (const Point<dim> &p_,
+                           const unsigned int component) const
+  {
+    const Point<dim> p = p_ - coordinate_system_offset;
+    const std_cxx11::array<double, dim> sp = GeometricUtilities::Coordinates::to_spherical(p);
+    const std_cxx11::array<double, dim> sg = sgradient(sp);
+    const std_cxx11::array<double, 6>   sh = shessian(sp);
+
+    // somewhat backwards, but we need cos/sin's for unit vectors
+    const double cos_theta = std::cos(sp[1]);
+    const double sin_theta = std::sin(sp[1]);
+    const double cos_phi   = std::cos(sp[2]);
+    const double sin_phi   = std::sin(sp[2]);
+    const double r         = sp[0];
+
+    Tensor<1,dim> unit_r, unit_theta, unit_phi;
+    set_unit_vectors(cos_theta, sin_theta,
+                     cos_phi, sin_phi,
+                     unit_r, unit_theta, unit_phi);
+
+    const double sin_phi2 = sin_phi*sin_phi;
+    const double r2 = r*r;
+    Assert (r != 0.,
+            ExcDivideByZero());
+
+    const double c_utheta2     = sg[0]/r +
+                                 ((sin_phi!= 0.) ?
+                                  (cos_phi*sg[2])/(r2*sin_phi)+sh[1]/(r2*sin_phi2) :
+                                  0.);
+    const double c_utheta_ur   = ((sin_phi != 0.) ?
+                                  (r*sh[3]-sg[1])/(r2*sin_phi) :
+                                  0.);
+    const double c_utheta_uphi = ((sin_phi != 0.) ?
+                                  (sh[5]*sin_phi-cos_phi*sg[1])/(r2*sin_phi2) :
+                                  0.);
+    const double c_ur2         = sh[0];
+    const double c_ur_uphi     = (r*sh[4]-sg[2])/r2;
+    const double c_uphi2       = (sh[2]+r*sg[0])/r2;
+
+    // go through each tensor product
+    SymmetricTensor<2,dim> res;
+
+    add_outer_product(res,
+                      c_utheta2,
+                      unit_theta);
+
+    add_outer_product(res,
+                      c_utheta_ur,
+                      unit_theta,
+                      unit_r);
+
+    add_outer_product(res,
+                      c_utheta_uphi,
+                      unit_theta,
+                      unit_phi);
+
+    add_outer_product(res,
+                      c_ur2,
+                      unit_r);
+
+    add_outer_product(res,
+                      c_ur_uphi,
+                      unit_r,
+                      unit_phi);
+
+    add_outer_product(res,
+                      c_uphi2,
+                      unit_phi);
+
+    return res;
+  }
+
+
+
+  template <int dim>
+  std::size_t
+  Spherical<dim>::memory_consumption () const
+  {
+    return sizeof(Spherical<dim>);
+  }
+
+
+
+  template <int dim>
+  double
+  Spherical<dim>::svalue(const std_cxx11::array<double, dim> &sp) const
+  {
+    AssertThrow(false,
+                ExcNotImplemented());
+    return 0.;
+  }
+
+
+
+  template <int dim>
+  std_cxx11::array<double, dim>
+  Spherical<dim>::sgradient(const std_cxx11::array<double, dim> &sp) const
+  {
+    AssertThrow(false,
+                ExcNotImplemented());
+    return std_cxx11::array<double,dim>();
+  }
+
+
+
+  template <int dim>
+  std_cxx11::array<double, 6>
+  Spherical<dim>::shessian (const std_cxx11::array<double, dim> &sp) const
+  {
+    AssertThrow(false,
+                ExcNotImplemented());
+    return std_cxx11::array<double, 6>();
+  }
+
+
+
+// explicit instantiations
+  template class Spherical<3>;
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/base/functions_spherical_01.cc b/tests/base/functions_spherical_01.cc
new file mode 100644 (file)
index 0000000..8ce3bcf
--- /dev/null
@@ -0,0 +1,171 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test value, gradient and hessian of a Functions::Spherical as
+// compared to radially symmetric fucntion.
+
+#include "../tests.h"
+#include <cmath>
+#include <deal.II/base/function_spherical.h>
+#include <deal.II/base/geometric_utilities.h>
+
+
+// reference function f(x) = exp (-Z*r)
+template<int dim>
+class ExpFunc : public Function<dim>
+{
+public:
+  ExpFunc(const Point<dim> &origin,
+          const double     &Z)
+    : Function<dim>(1),
+      origin(origin),
+      Z(Z)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0) const
+  {
+    Tensor<1,dim> dist = point-origin;
+    const double r = dist.norm();
+    return std::exp(-Z*r);
+  }
+
+  virtual Tensor< 1, dim> gradient (const Point<dim > &p,
+                                    const unsigned int component=0) const
+  {
+    Tensor<1,dim> dist = p-origin;
+    const double r = dist.norm();
+    Assert (r>0.0, ExcMessage("r is not positive"));
+    dist/=r;
+    return -Z*std::exp(-Z*r)*dist;
+  }
+
+  virtual SymmetricTensor<2,dim> hessian (const Point<dim> &p,
+                                          const unsigned int component=0) const
+  {
+    Tensor<1,dim> dir = p-origin;
+    const double r = dir.norm();
+    Assert (r>0.0, ExcMessage("r is not positive"));
+    dir/=r;
+    SymmetricTensor<2,dim> dir_x_dir;
+    for (unsigned int i=0; i<dim; i++)
+      for (unsigned int j=i; j<dim; j++)
+        dir_x_dir[i][j] = dir[i] * dir[j];
+
+    return Z*std::exp(-Z*r)*( (Z+1.0/r)*dir_x_dir - unit_symmetric_tensor<dim>()/r );
+  }
+
+private:
+  const Point<dim> origin;
+  const double Z;
+};
+
+
+// same as above but using Functions::Spherical
+template<int dim>
+class ExpFunc2 : public Functions::Spherical<dim>
+{
+public:
+  ExpFunc2(const Point<dim> &origin,
+           const double     &Z)
+    : Functions::Spherical<dim>(origin),
+      Z(Z)
+  {}
+
+private:
+  virtual double svalue(const std_cxx11::array<double, dim> &sp) const
+  {
+    return std::exp(-Z*sp[0]);
+  }
+
+  virtual std_cxx11::array<double, dim> sgradient(const std_cxx11::array<double, dim> &sp) const
+  {
+    std_cxx11::array<double, dim> res;
+    res[0] = -Z*std::exp(-Z*sp[0]);
+    for (unsigned int i=1; i < dim; i++)
+      res[i] = 0.;
+    return res;
+  }
+
+  virtual std_cxx11::array<double, 6> shessian (const std_cxx11::array<double, dim> &sp) const
+  {
+    std_cxx11::array<double, 6> res;
+    res[0] = Z*Z*std::exp(-Z*sp[0]);
+    for (unsigned int i=1; i < 6; i++)
+      res[i] = 0.;
+    return res;
+  }
+
+  const double Z;
+};
+
+template <int dim>
+void check()
+{
+  Point<dim> center;
+  const double Z = 2.5;
+  center[1] = 2.0;
+  if (dim>2)
+    center[2] = -1.5;
+
+  ExpFunc<dim>  func(center, Z);
+  ExpFunc2<dim> func2(center, Z);
+
+  for (double r = 0.1; r < 10; r+= 0.35)
+    for (double theta = 0; theta < 2*numbers::PI; theta+= numbers::PI/3.)
+      for (double phi = 0.01; phi <= numbers::PI; phi+= numbers::PI/4.)
+        {
+          std_cxx11::array<double, dim> sp;
+          sp[0] = r;
+          sp[1] = theta;
+          sp[2] = phi;
+          Point<dim> p = GeometricUtilities::Coordinates::from_spherical(sp);
+          for (unsigned int i = 0; i < dim; i++)
+            p[i]+= center[i];
+
+          // check values:
+          const double v1 = func.value(p);
+          const double v2 = func2.value(p);
+          AssertThrow( std::fabs(v1 - v2) <= std::abs(v1)*1e-10,
+                       ExcInternalError());
+
+          // check gradients:
+          const Tensor<1,dim> g1 = func.gradient(p);
+          const Tensor<1,dim> g2 = func2.gradient(p);
+          const Tensor<1,dim> gd = g1-g2;
+          AssertThrow ( gd.norm() <= g1.norm() * 1e-10,
+                        ExcInternalError());
+
+          // check hessian:
+          const SymmetricTensor<2,dim> h1 = func.hessian(p);
+          const SymmetricTensor<2,dim> h2 = func2.hessian(p);
+          const SymmetricTensor<2,dim> dh = h1 - h2;
+          AssertThrow ( dh.norm() <= h1.norm() * 1e-10,
+                        ExcInternalError());
+        }
+  deallog << "OK"<< std::endl;
+}
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  check<3>();
+}
+
diff --git a/tests/base/functions_spherical_01.output b/tests/base/functions_spherical_01.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/base/functions_spherical_02.cc b/tests/base/functions_spherical_02.cc
new file mode 100644 (file)
index 0000000..91434ce
--- /dev/null
@@ -0,0 +1,223 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test value, gradient and Hessian of a Functions::Spherical
+// for the case of Px orbital
+// (sqrt(3/4pi) sin(phi) cos(theta) == sqrt(3/pi)*x/r)
+// multiplied with r*r.
+
+/* MWE in Maxima:
+
+spherical:
+f: r*r*cos(theta)*sin(phi);
+[diff(f,r), diff(f,theta), diff(f,phi)];
+[diff(f,r,2), diff(f,theta,2), diff(f,phi,2), diff(f,r,1,theta,1), diff(f,r,1,phi,1), diff (f,theta,1,phi,1)];
+
+Cartesian:
+
+assume(r>0);
+srule: [r= sqrt(x*x+y*y+z*z)];
+srule2: [sqrt(x*x+y*y+z*z) = r, x*x+y*y+z*z = r^2, expand((x*x+y*y+z*z)^2) = r^4];
+f: subst(srule,r*r*x/r);
+df: fullratsimp([diff(f,x), diff(f,y), diff(f,z)]);
+df2 :(subst(srule2,expand(df)));
+ddf : fullratsimp([diff(f,x,2), diff(f,y,2), diff(f,z,2), diff(f,x,1,y,1), diff(f,x,1,z,1), diff(f,y,1,z,1)]);
+ddf2: (subst(srule2,expand(ddf)));
+
+ */
+
+#include "../tests.h"
+#include <cmath>
+#include <deal.II/base/function_spherical.h>
+#include <deal.II/base/geometric_utilities.h>
+
+
+
+template<int dim>
+class RefFunc : public Function<dim>
+{
+public:
+  RefFunc(const Point<dim> origin = Point<dim>())
+    : Function<dim>(1),
+      origin(origin)
+  {}
+
+  virtual double value(const Point<dim> &point,
+                       const unsigned int component = 0) const
+  {
+    Tensor<1,dim> dist = point-origin;
+    const double r = dist.norm();
+    return r*dist[0];
+  }
+
+  virtual Tensor< 1, dim> gradient (const Point<dim > &p,
+                                    const unsigned int component=0) const
+  {
+    Tensor<1,dim> dist = p-origin;
+    const double r = dist.norm();
+    Assert (r>0.0, ExcMessage("r is not positive"));
+    Tensor<1,dim> res;
+    const double x = dist[0];
+    const double y = dist[1];
+    const double z = dist[2];
+    const double x2= x*x;
+    const double y2= y*y;
+    const double z2= z*z;
+
+    res[0] = z2/r+y2/r+(2.*x2)/r;
+    res[1] = (x*y)/r;
+    res[2] = (x*z)/r;
+
+    return res;
+  }
+
+  virtual SymmetricTensor<2,dim> hessian (const Point<dim> &p,
+                                          const unsigned int component=0) const
+  {
+    const Tensor<1,dim> dist = p-origin;
+    const double r = dist.norm();
+    Assert (r>0.0, ExcMessage("r is not positive"));
+    const double x = dist[0];
+    const double y = dist[1];
+    const double z = dist[2];
+    const double z2= z*z;
+    const double y2= y*y;
+    const double x2= x*x;
+    const double x3= x2*x;
+    const double z3= z2*z;
+    const double y3= y2*y;
+    const double r3= r*r*r;
+
+
+    SymmetricTensor<2,dim> res;
+    res[0][0] = (3.*x*z2)/r3+(3*x*y2)/r3+(2.*x3)/r3;
+    res[1][1] = (x*z2)/r3+x3/r3;
+    res[2][2] = (x*y2)/r3+x3/r3;
+    res[0][1] = (y*z2)/r3+y3/r3;
+    res[0][2] = z3/r3+(y2*z)/r3;
+    res[1][2] =-(x*y*z)/r3;
+
+    return res;
+  }
+
+private:
+  const Point<dim> origin;
+};
+
+
+// same as above but using Functions::Spherical
+template<int dim>
+class SphFunc : public Functions::Spherical<dim>
+{
+public:
+  SphFunc(const Point<dim> origin = Point<dim>())
+    : Functions::Spherical<dim>(origin)
+  {}
+
+private:
+  virtual double svalue(const std_cxx11::array<double, dim> &sp) const
+  {
+    return sp[0]*sp[0]*std::cos(sp[1])*std::sin(sp[2]);
+  }
+
+  virtual std_cxx11::array<double, dim> sgradient(const std_cxx11::array<double, dim> &sp) const
+  {
+    std_cxx11::array<double, dim> res;
+    const double r     = sp[0];
+    const double theta = sp[1];
+    const double phi   = sp[2];
+    res[0] = 2.*sin(phi)*r*cos(theta);
+    res[1] = -sin(phi)*r*r*sin(theta);
+    res[2] =  cos(phi)*r*r*cos(theta);
+    return res;
+  }
+
+  virtual std_cxx11::array<double, 6> shessian (const std_cxx11::array<double, dim> &sp) const
+  {
+    std_cxx11::array<double, 6> res;
+    const double r = sp[0];
+    const double theta = sp[1];
+    const double phi = sp[2];
+    const double r2 = r*r;
+    res[0] = 2.*sin(phi)*cos(theta);
+    res[1] = -sin(phi)*r2*cos(theta);
+    res[2] = -sin(phi)*r2*cos(theta);
+    res[3] = -2.*sin(phi)*r*sin(theta);
+    res[4] = 2.*cos(phi)*r*cos(theta);
+    res[5] = -cos(phi)*r2*sin(theta);
+    return res;
+  }
+};
+
+template <int dim>
+void check()
+{
+  Point<dim> center;
+  center[1] = 2.0;
+  if (dim>2)
+    center[2] = -1.5;
+
+  RefFunc<dim>  func(center);
+  SphFunc<dim>  func2(center);
+
+  for (double r = 0.1; r < 10; r+= 0.35)
+    for (double theta = 0; theta < 2*numbers::PI; theta+= numbers::PI/3.)
+      for (double phi = 0.01; phi <= numbers::PI; phi+= numbers::PI/4.)
+        {
+          std_cxx11::array<double, dim> sp;
+          sp[0] = r;
+          sp[1] = theta;
+          sp[2] = phi;
+          Point<dim> p = GeometricUtilities::Coordinates::from_spherical(sp);
+          for (unsigned int i = 0; i < dim; i++)
+            p[i]+= center[i];
+
+          // check values:
+          const double v1 = func.value(p);
+          const double v2 = func2.value(p);
+          AssertThrow( std::fabs(v1 - v2) <= std::abs(v1)*1e-10,
+                       ExcInternalError());
+
+          // check gradients:
+          const Tensor<1,dim> g1 = func.gradient(p);
+          const Tensor<1,dim> g2 = func2.gradient(p);
+          const Tensor<1,dim> gd = g1-g2;
+          AssertThrow ( gd.norm() <= g1.norm() * 1e-10,
+                        ExcInternalError());
+
+
+          // check hessian:
+          const SymmetricTensor<2,dim> h1 = func.hessian(p);
+          const SymmetricTensor<2,dim> h2 = func2.hessian(p);
+          const SymmetricTensor<2,dim> dh = h1 - h2;
+          AssertThrow ( dh.norm() <= h1.norm() * 1e-10,
+                        ExcInternalError());
+
+
+        }
+  deallog << "OK"<< std::endl;
+}
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  check<3>();
+}
+
diff --git a/tests/base/functions_spherical_02.output b/tests/base/functions_spherical_02.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.