From 21723a95531bd80b3a400f3a6b7f2d96ba92179b Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 8 Sep 2015 08:50:48 +0200 Subject: [PATCH] added Function::hessian() fixed CosineFunction in function_lib to comply. --- doc/news/changes.h | 7 ++++ include/deal.II/base/function.h | 29 +++++++++++++ include/deal.II/base/function.templates.h | 51 +++++++++++++++++++++++ include/deal.II/base/function_lib.h | 6 +-- source/base/function_lib.cc | 20 ++++----- 5 files changed, 98 insertions(+), 15 deletions(-) diff --git a/doc/news/changes.h b/doc/news/changes.h index 3c7e850128..8d6842ed75 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -160,6 +160,7 @@ inconvenience this causes. be sliced back to a LinearOperator.
(Matthias Maier, 2015/08/27) +
  • Improved: Support for complex number types throughout the library. Several parts of the library have been reorganized to support complex @@ -271,12 +272,18 @@ inconvenience this causes.
      +
    1. Introduced Hessian-related functions to the Function class. +
      + (Denis Davydov, 2015/09/08) +
      2. +
    2. Memory consumption during compilation has been reduced by splitting instantiation files. For this make_instantiations now supports additional logic to split the the instantiations in .inst files into groups. This is used in fe_values.cc, error_estimator.cc, and others.
      (Timo Heister, 2015/09/05) +
    3. New: There is now a function SparsityPattern::print_svg() which prints the sparsity of the matrix in a .svg file which can be opened in a web browser. diff --git a/include/deal.II/base/function.h b/include/deal.II/base/function.h index d0e99b4671..2f30365480 100644 --- a/include/deal.II/base/function.h +++ b/include/deal.II/base/function.h @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -309,6 +310,34 @@ public: virtual void vector_laplacian_list (const std::vector > &points, std::vector > &values) const; + /** + * Compute the Hessian of a given component at point p, + * that is the gradient of the gradient of the function. + */ + virtual SymmetricTensor<2,dim,Number> hessian (const Point &p, + const unsigned int component = 0) const; + + /** + * Compute the Hessian of all components at point p and store + * them in values. + */ + virtual void vector_hessian (const Point &p, + std::vector> &values) const; + + /** + * Compute the Hessian of one component at a set of points. + */ + virtual void hessian_list (const std::vector > &points, + std::vector > &values, + const unsigned int component = 0) const; + + /** + * Compute the Hessians of all components at a set of points. + */ + virtual void vector_hessian_list (const std::vector > &points, + std::vector > > &values) const; + + /** * Return an estimate for the memory consumption, in bytes, of this object. * This is not exact (but will usually be close) because calculating the diff --git a/include/deal.II/base/function.templates.h b/include/deal.II/base/function.templates.h index 0cef312b70..6afc6c8c19 100644 --- a/include/deal.II/base/function.templates.h +++ b/include/deal.II/base/function.templates.h @@ -233,6 +233,57 @@ void Function::vector_laplacian_list ( } +template +SymmetricTensor<2,dim,Number> Function::hessian (const Point &, + const unsigned int) const +{ + Assert (false, ExcPureFunctionCalled()); + return SymmetricTensor<2,dim,Number>(); +} + + +template +void Function::vector_hessian ( + const Point &p, + std::vector > &v) const +{ + AssertDimension(v.size(), this->n_components); + for (unsigned int i=0; in_components; ++i) + v[i] = hessian(p, i); +} + + +template +void Function::hessian_list ( + const std::vector > &points, + std::vector > &hessians, + const unsigned int component) const +{ + Assert (hessians.size() == points.size(), + ExcDimensionMismatch(hessians.size(), points.size())); + + for (unsigned int i=0; i +void Function::vector_hessian_list ( + const std::vector > &points, + std::vector > > &hessians) const +{ + Assert (hessians.size() == points.size(), + ExcDimensionMismatch(hessians.size(), points.size())); + + for (unsigned int i=0; i std::size_t diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index 8ce8ee203c..d05b71e1b8 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -236,14 +236,14 @@ namespace Functions /** * Second derivatives at a single point. */ - virtual Tensor<2,dim> hessian (const Point &p, - const unsigned int component = 0) const; + virtual SymmetricTensor<2,dim> hessian (const Point &p, + const unsigned int component = 0) const; /** * Second derivatives at multiple points. */ virtual void hessian_list (const std::vector > &points, - std::vector > &hessians, + std::vector > &hessians, const unsigned int component = 0) const; }; diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc index 20fcd3f331..e86263e0f6 100644 --- a/source/base/function_lib.cc +++ b/source/base/function_lib.cc @@ -586,13 +586,13 @@ namespace Functions } template - Tensor<2,dim> + SymmetricTensor<2,dim> CosineFunction::hessian (const Point &p, const unsigned int) const { const double pi2 = M_PI_2*M_PI_2; - Tensor<2,dim> result; + SymmetricTensor<2,dim> result; switch (dim) { case 1: @@ -605,8 +605,8 @@ namespace Functions const double sisi = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)); result[0][0] = coco; result[1][1] = coco; + // for SymmetricTensor we assign [ij] and [ji] simultaneously: result[0][1] = sisi; - result[1][0] = sisi; } break; case 3: @@ -620,12 +620,10 @@ namespace Functions result[0][0] = cococo; result[1][1] = cococo; result[2][2] = cococo; + // for SymmetricTensor we assign [ij] and [ji] simultaneously: 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: @@ -636,8 +634,8 @@ namespace Functions template void - CosineFunction::hessian_list (const std::vector > &points, - std::vector > &hessians, + CosineFunction::hessian_list (const std::vector > &points, + std::vector > &hessians, const unsigned int) const { Assert (hessians.size() == points.size(), @@ -660,8 +658,8 @@ namespace Functions const double sisi = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)); hessians[i][0][0] = coco; hessians[i][1][1] = coco; + // for SymmetricTensor we assign [ij] and [ji] simultaneously: hessians[i][0][1] = sisi; - hessians[i][1][0] = sisi; } break; case 3: @@ -675,12 +673,10 @@ namespace Functions hessians[i][0][0] = cococo; hessians[i][1][1] = cococo; hessians[i][2][2] = cococo; + // for SymmetricTensor we assign [ij] and [ji] simultaneously: 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: -- 2.39.5