From: Guido Kanschat Date: Mon, 5 Mar 2001 17:49:45 +0000 (+0000) Subject: second derivatives of cosine function X-Git-Tag: v8.0.0~19657 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9ee2674ff6190b6590a2f6081e4aed9a4b3fe1a8;p=dealii.git second derivatives of cosine function git-svn-id: https://svn.dealii.org/trunk@4108 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/function_lib.h b/deal.II/base/include/base/function_lib.h index c7aab5d4eb..6d5b489bfe 100644 --- a/deal.II/base/include/base/function_lib.h +++ b/deal.II/base/include/base/function_lib.h @@ -242,6 +242,19 @@ class CosineFunction : public Function virtual void laplacian_list (const typename std::vector > &points, std::vector &values, const unsigned int component = 0) const; + + /** + * Gradient at a single point. + */ + virtual Tensor<2,dim> hessian (const Point &p, + const unsigned int component = 0) const; + + /** + * Gradients at multiple points. + */ + virtual void hessian_list (const typename std::vector > &points, + typename std::vector > &hessians, + const unsigned int component = 0) const; }; diff --git a/deal.II/base/source/function_lib.cc b/deal.II/base/source/function_lib.cc index afe6e5ca2d..a9d4a1ee72 100644 --- a/deal.II/base/source/function_lib.cc +++ b/deal.II/base/source/function_lib.cc @@ -12,6 +12,7 @@ //---------------------------- function_lib.cc --------------------------- +#include #include #include @@ -519,6 +520,110 @@ CosineFunction::gradient_list (const typename std::vector > &poi } } +template +Tensor<2,dim> +CosineFunction::hessian (const Point &p, + const unsigned int) const +{ + const double pi2 = M_PI_2*M_PI_2; + + Tensor<2,dim> result; + switch(dim) + { + case 1: + result[0][0] = -pi2* cos(M_PI_2*p(0)); + break; + case 2: + if (true) + { + const double coco = -pi2*cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)); + const double sisi = pi2*sin(M_PI_2*p(0)) * sin(M_PI_2*p(1)); + result[0][0] = coco; + result[1][1] = coco; + result[0][1] = sisi; + result[1][0] = sisi; + } + break; + case 3: + if (true) + { + const double cococo = -pi2*cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * cos(M_PI_2*p(2)); + const double sisico = pi2*sin(M_PI_2*p(0)) * sin(M_PI_2*p(1)) * cos(M_PI_2*p(2)); + const double sicosi = pi2*sin(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * sin(M_PI_2*p(2)); + const double cosisi = pi2*cos(M_PI_2*p(0)) * sin(M_PI_2*p(1)) * sin(M_PI_2*p(2)); + + result[0][0] = cococo; + result[1][1] = cococo; + result[2][2] = cococo; + result[0][1] = sisico; + result[1][0] = sisico; + result[0][2] = sicosi; + result[2][0] = sicosi; + result[1][2] = cosisi; + result[2][1] = cosisi; + } + break; + default: + Assert(false, ExcNotImplemented()); + } + return result; +} + +template +void +CosineFunction::hessian_list (const typename std::vector > &points, + typename std::vector > &hessians, + const unsigned int) const +{ + Assert (hessians.size() == points.size(), + ExcDimensionMismatch(hessians.size(), points.size())); + + const double pi2 = M_PI_2*M_PI_2; + + for (unsigned int i=0;i& p = points[i]; + switch(dim) + { + case 1: + hessians[i][0][0] = -pi2* cos(M_PI_2*p(0)); + break; + case 2: + if (true) + { + const double coco = -pi2*cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)); + const double sisi = pi2*sin(M_PI_2*p(0)) * sin(M_PI_2*p(1)); + hessians[i][0][0] = coco; + hessians[i][1][1] = coco; + hessians[i][0][1] = sisi; + hessians[i][1][0] = sisi; + } + break; + case 3: + if (true) + { + const double cococo = -pi2*cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * cos(M_PI_2*p(2)); + const double sisico = pi2*sin(M_PI_2*p(0)) * sin(M_PI_2*p(1)) * cos(M_PI_2*p(2)); + const double sicosi = pi2*sin(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * sin(M_PI_2*p(2)); + const double cosisi = pi2*cos(M_PI_2*p(0)) * sin(M_PI_2*p(1)) * sin(M_PI_2*p(2)); + + hessians[i][0][0] = cococo; + hessians[i][1][1] = cococo; + hessians[i][2][2] = cococo; + hessians[i][0][1] = sisico; + hessians[i][1][0] = sisico; + hessians[i][0][2] = sicosi; + hessians[i][2][0] = sicosi; + hessians[i][1][2] = cosisi; + hessians[i][2][1] = cosisi; + } + break; + default: + Assert(false, ExcNotImplemented()); + } + } +} + ////////////////////////////////////////////////////////////////////// template