--- /dev/null
+New: Implement Functions::Spherical to represent a function given in spherical
+coordinates.
+<br>
+(Denis Davydov, 2017/01/28)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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> ¢er = 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
+
function_lib.cc
function_lib_cutoff.cc
function_parser.cc
+ function_spherical.cc
function_time.cc
geometry_info.cc
geometric_utilities.cc
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
+
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
+
--- /dev/null
+
+DEAL::OK
projects / dealii.git / commitdiff