From ec6c08b586db851477e59e3e8abf164d9c7bd8c2 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 30 Jan 2017 18:22:55 +0100 Subject: [PATCH] add spherical function --- doc/news/changes/minor/20170128DenisDavydov | 4 + include/deal.II/base/function_spherical.h | 118 +++++++ source/base/CMakeLists.txt | 1 + source/base/function_spherical.cc | 348 ++++++++++++++++++++ tests/base/functions_spherical_01.cc | 171 ++++++++++ tests/base/functions_spherical_01.output | 2 + tests/base/functions_spherical_02.cc | 223 +++++++++++++ tests/base/functions_spherical_02.output | 2 + 8 files changed, 869 insertions(+) create mode 100644 doc/news/changes/minor/20170128DenisDavydov create mode 100644 include/deal.II/base/function_spherical.h create mode 100644 source/base/function_spherical.cc create mode 100644 tests/base/functions_spherical_01.cc create mode 100644 tests/base/functions_spherical_01.output create mode 100644 tests/base/functions_spherical_02.cc create mode 100644 tests/base/functions_spherical_02.output diff --git a/doc/news/changes/minor/20170128DenisDavydov b/doc/news/changes/minor/20170128DenisDavydov new file mode 100644 index 0000000000..129c8d12bc --- /dev/null +++ b/doc/news/changes/minor/20170128DenisDavydov @@ -0,0 +1,4 @@ +New: Implement Functions::Spherical to represent a function given in spherical +coordinates. +
+(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 index 0000000000..1ed772627e --- /dev/null +++ b/include/deal.II/base/function_spherical.h @@ -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 + +#include +#include +#include + +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 + class Spherical : public Function + { + public: + /** + * Constructor which should be provided with @p center defining the origin + * of the coordinate system. + */ + Spherical(const Point ¢er = Point()); + + /** + * 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 &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 &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 &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 &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 sgradient(const std_cxx11::array &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 shessian (const std_cxx11::array &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 + diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 66ece18044..0bff1fa583 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -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 index 0000000000..56b8117e21 --- /dev/null +++ b/source/base/function_spherical.cc @@ -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 +#include +#include + +#include +#include + +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 + 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 + 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 + 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 + Spherical::Spherical(const Point &p) + : + Function(1), + coordinate_system_offset(p) + { + AssertThrow(dim==3, + ExcNotImplemented()); + } + + + + template + double + Spherical::value (const Point &p_, + const unsigned int) const + { + const Point p = p_ - coordinate_system_offset; + const std_cxx11::array sp = GeometricUtilities::Coordinates::to_spherical(p); + return svalue(sp); + } + + + + template + Tensor<1,dim> + Spherical::gradient (const Point &p_, + const unsigned int) const + { + const Point p = p_ - coordinate_system_offset; + const std_cxx11::array sp = GeometricUtilities::Coordinates::to_spherical(p); + const std_cxx11::array 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 + SymmetricTensor<2,dim> + Spherical::hessian (const Point &p_, + const unsigned int component) const + { + const Point p = p_ - coordinate_system_offset; + const std_cxx11::array sp = GeometricUtilities::Coordinates::to_spherical(p); + const std_cxx11::array sg = sgradient(sp); + const std_cxx11::array 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 + std::size_t + Spherical::memory_consumption () const + { + return sizeof(Spherical); + } + + + + template + double + Spherical::svalue(const std_cxx11::array &sp) const + { + AssertThrow(false, + ExcNotImplemented()); + return 0.; + } + + + + template + std_cxx11::array + Spherical::sgradient(const std_cxx11::array &sp) const + { + AssertThrow(false, + ExcNotImplemented()); + return std_cxx11::array(); + } + + + + template + std_cxx11::array + Spherical::shessian (const std_cxx11::array &sp) const + { + AssertThrow(false, + ExcNotImplemented()); + return std_cxx11::array(); + } + + + +// 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 index 0000000000..8ce3bcf600 --- /dev/null +++ b/tests/base/functions_spherical_01.cc @@ -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 +#include +#include + + +// reference function f(x) = exp (-Z*r) +template +class ExpFunc : public Function +{ +public: + ExpFunc(const Point &origin, + const double &Z) + : Function(1), + origin(origin), + Z(Z) + {} + + virtual double value(const Point &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 &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 &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()/r ); + } + +private: + const Point origin; + const double Z; +}; + + +// same as above but using Functions::Spherical +template +class ExpFunc2 : public Functions::Spherical +{ +public: + ExpFunc2(const Point &origin, + const double &Z) + : Functions::Spherical(origin), + Z(Z) + {} + +private: + virtual double svalue(const std_cxx11::array &sp) const + { + return std::exp(-Z*sp[0]); + } + + virtual std_cxx11::array sgradient(const std_cxx11::array &sp) const + { + std_cxx11::array 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 shessian (const std_cxx11::array &sp) const + { + std_cxx11::array 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 +void check() +{ + Point center; + const double Z = 2.5; + center[1] = 2.0; + if (dim>2) + center[2] = -1.5; + + ExpFunc func(center, Z); + ExpFunc2 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 sp; + sp[0] = r; + sp[1] = theta; + sp[2] = phi; + Point 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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/base/functions_spherical_01.output @@ -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 index 0000000000..91434ce200 --- /dev/null +++ b/tests/base/functions_spherical_02.cc @@ -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 +#include +#include + + + +template +class RefFunc : public Function +{ +public: + RefFunc(const Point origin = Point()) + : Function(1), + origin(origin) + {} + + virtual double value(const Point &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 &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 &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 origin; +}; + + +// same as above but using Functions::Spherical +template +class SphFunc : public Functions::Spherical +{ +public: + SphFunc(const Point origin = Point()) + : Functions::Spherical(origin) + {} + +private: + virtual double svalue(const std_cxx11::array &sp) const + { + return sp[0]*sp[0]*std::cos(sp[1])*std::sin(sp[2]); + } + + virtual std_cxx11::array sgradient(const std_cxx11::array &sp) const + { + std_cxx11::array 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 shessian (const std_cxx11::array &sp) const + { + std_cxx11::array 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 +void check() +{ + Point center; + center[1] = 2.0; + if (dim>2) + center[2] = -1.5; + + RefFunc func(center); + SphFunc 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 sp; + sp[0] = r; + sp[1] = theta; + sp[2] = phi; + Point 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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/base/functions_spherical_02.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5