added Function::hessian()

fixed CosineFunction in function_lib to comply.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 3c7e850128b3b7c6a1aee4a280e6c19eafec815d..8d6842ed75ddfb746ef7d9f38fd3e7c7690f87c8 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -160,6 +160,7 @@ inconvenience this causes.
   be sliced back to a LinearOperator.
   <br>
   (Matthias Maier, 2015/08/27)
+  </li>
 
   <li> 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.
 
 
 <ol>
+  <li> Introduced Hessian-related functions to the Function class.
+  <br>
+  (Denis Davydov, 2015/09/08)
+  </li>
+
   <li> 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.
   <br>
   (Timo Heister, 2015/09/05)
+  </li>
 
   <li> 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 d0e99b4671939019a10cc1f03d5be949adf49692..2f30365480498dad2d9c05f98f8f5592fa50e8e4 100644 (file)
--- a/include/deal.II/base/function.h
+++ b/include/deal.II/base/function.h
@@ -22,6 +22,7 @@
 #include <deal.II/base/function_time.h>
 #include <deal.II/base/subscriptor.h>
 #include <deal.II/base/tensor.h>
+#include <deal.II/base/symmetric_tensor.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/std_cxx11/function.h>
 
@@ -309,6 +310,34 @@ public:
   virtual void vector_laplacian_list (const std::vector<Point<dim> > &points,
                                       std::vector<Vector<Number> >   &values) const;
 
+  /**
+   * Compute the Hessian of a given component at point <tt>p</tt>,
+   * that is the gradient of the gradient of the function.
+   */
+  virtual SymmetricTensor<2,dim,Number> hessian (const Point<dim>   &p,
+                                                 const unsigned int          component = 0) const;
+
+  /**
+   * Compute the Hessian of all components at point <tt>p</tt> and store
+   * them in <tt>values</tt>.
+   */
+  virtual void vector_hessian (const Point<dim>                           &p,
+                               std::vector<SymmetricTensor<2,dim,Number>> &values) const;
+
+  /**
+   * Compute the Hessian of one component at a set of points.
+   */
+  virtual void hessian_list (const std::vector<Point<dim> >              &points,
+                             std::vector<SymmetricTensor<2,dim,Number> > &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<Point<dim> >                            &points,
+                                    std::vector<std::vector<SymmetricTensor<2,dim,Number> > > &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 0cef312b70a4d0046ddfc00558dacdd0dfdfcdec..6afc6c8c193f36bc9c5da9e722ffc83227ca5a8f 100644 (file)
--- a/include/deal.II/base/function.templates.h
+++ b/include/deal.II/base/function.templates.h
@@ -233,6 +233,57 @@ void Function<dim, Number>::vector_laplacian_list (
 }
 
 
+template <int dim, typename Number>
+SymmetricTensor<2,dim,Number> Function<dim, Number>::hessian (const Point<dim> &,
+    const unsigned int) const
+{
+  Assert (false, ExcPureFunctionCalled());
+  return SymmetricTensor<2,dim,Number>();
+}
+
+
+template <int dim, typename Number>
+void Function<dim, Number>::vector_hessian (
+  const Point<dim> &p,
+  std::vector<SymmetricTensor<2,dim,Number> > &v) const
+{
+  AssertDimension(v.size(), this->n_components);
+  for (unsigned int i=0; i<this->n_components; ++i)
+    v[i] = hessian(p, i);
+}
+
+
+template <int dim, typename Number>
+void Function<dim, Number>::hessian_list (
+  const std::vector<Point<dim> >     &points,
+  std::vector<SymmetricTensor<2,dim,Number> > &hessians,
+  const unsigned int                  component) const
+{
+  Assert (hessians.size() == points.size(),
+          ExcDimensionMismatch(hessians.size(), points.size()));
+
+  for (unsigned int i=0; i<points.size(); ++i)
+    hessians[i] = hessian(points[i], component);
+}
+
+
+template <int dim, typename Number>
+void Function<dim, Number>::vector_hessian_list (
+  const std::vector<Point<dim> >                   &points,
+  std::vector<std::vector<SymmetricTensor<2,dim,Number> > > &hessians) const
+{
+  Assert (hessians.size() == points.size(),
+          ExcDimensionMismatch(hessians.size(), points.size()));
+
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      Assert (hessians[i].size() == n_components,
+              ExcDimensionMismatch(hessians[i].size(), n_components));
+      vector_hessian (points[i], hessians[i]);
+    }
+}
+
+
 
 template <int dim, typename Number>
 std::size_t

diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h
index 8ce8ee203ce16ee2f8648772b9abebeb5d539e27..d05b71e1b84fe53c4ebbc096d31e2d8275be5958 100644 (file)
--- 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<dim>   &p,
-                                   const unsigned int  component = 0) const;
+    virtual SymmetricTensor<2,dim> hessian (const Point<dim>   &p,
+                                            const unsigned int  component = 0) const;
 
     /**
      * Second derivatives at multiple points.
      */
     virtual void hessian_list (const std::vector<Point<dim> > &points,
-                               std::vector<Tensor<2,dim> >    &hessians,
+                               std::vector<SymmetricTensor<2,dim> >    &hessians,
                                const unsigned int              component = 0) const;
   };
 

diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc
index 20fcd3f331dbdd268b1c2ab62862d99a55a5a03c..e86263e0f6d71ee14340e77587e6f77f3085dc03 100644 (file)
--- a/source/base/function_lib.cc
+++ b/source/base/function_lib.cc
@@ -586,13 +586,13 @@ namespace Functions
   }
 
   template<int dim>
-  Tensor<2,dim>
+  SymmetricTensor<2,dim>
   CosineFunction<dim>::hessian (const Point<dim>   &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<int dim>
   void
-  CosineFunction<dim>::hessian_list (const std::vector<Point<dim> > &points,
-                                     std::vector<Tensor<2,dim> >    &hessians,
+  CosineFunction<dim>::hessian_list (const std::vector<Point<dim> >       &points,
+                                     std::vector<SymmetricTensor<2,dim> > &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: