]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement point_gradient similarly to point_value
authorDaniel Arndt <d.arndt@math.uni-goettingen.de>
Sun, 26 Apr 2015 21:19:54 +0000 (23:19 +0200)
committerDaniel Arndt <d.arndt@math.uni-goettingen.de>
Sun, 3 May 2015 22:42:34 +0000 (00:42 +0200)
13 files changed:
include/deal.II/numerics/vector_tools.h
include/deal.II/numerics/vector_tools.templates.h
source/numerics/CMakeLists.txt
source/numerics/vector_tools_point_gradient.cc [new file with mode: 0644]
source/numerics/vector_tools_point_gradient.inst.in [new file with mode: 0644]
tests/bits/point_gradient_01.cc [new file with mode: 0644]
tests/bits/point_gradient_01.output [new file with mode: 0644]
tests/bits/point_gradient_02.cc [new file with mode: 0644]
tests/bits/point_gradient_02.output [new file with mode: 0644]
tests/bits/point_gradient_hp_01.cc [new file with mode: 0644]
tests/bits/point_gradient_hp_01.output [new file with mode: 0644]
tests/bits/point_gradient_hp_02.cc [new file with mode: 0644]
tests/bits/point_gradient_hp_02.output [new file with mode: 0644]

index 78c3c1953afbab4fb5201f7e4ad9ad568beaaba8..7ce3c469e65b46cfb1a8d431959fa45d4391e93c 100644 (file)
@@ -2106,6 +2106,130 @@ namespace VectorTools
                const InVector                            &fe_function,
                const Point<spacedim>                     &point);
 
+  /**
+    * Evaluate a possibly vector-valued finite element function defined by the
+    * given DoFHandler and nodal vector at the given point, and return the
+    * (vector) gradient of this function through the last argument.
+    *
+    * This is a wrapper function using a Q1-mapping for cell boundaries to call
+    * the other point_difference() function.
+    *
+    * @note If the cell in which the point is found is not locally owned, an
+    * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+    */
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const DoFHandler<dim,spacedim>    &dof,
+                  const InVector                    &fe_function,
+                  const Point<spacedim>             &point,
+                  std::vector<Tensor<1, spacedim> > &value);
+
+  /**
+   * Same as above for hp.
+   *
+   * @note If the cell in which the point is found is not locally owned, an
+   * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+   */
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
+                  const InVector                     &fe_function,
+                  const Point<spacedim>              &point,
+                  std::vector<Tensor<1, spacedim> >  &value);
+
+  /**
+   * Evaluate a scalar finite element function defined by the given DoFHandler
+   * and nodal vector at the given point, and return the gradient of this
+   * function.
+   *
+   * Compared with the other function of the same name, this is a wrapper
+   * function using a Q1-mapping for cells.
+   *
+   * @note If the cell in which the point is found is not locally owned, an
+   * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+   */
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const DoFHandler<dim,spacedim> &dof,
+                  const InVector                 &fe_function,
+                  const Point<spacedim>          &point);
+
+  /**
+   * Same as above for hp.
+   *
+   * @note If the cell in which the point is found is not locally owned, an
+   * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+   */
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
+                  const InVector                     &fe_function,
+                  const Point<spacedim>              &point);
+
+  /**
+   * Evaluate a possibly vector-valued finite element function defined by the
+   * given DoFHandler and nodal vector at the given point, and return the
+   * gradients of this function through the last argument.
+   *
+   * Compared with the other function of the same name, this function uses an
+   * arbitrary mapping for evaluation.
+   *
+   * @note If the cell in which the point is found is not locally owned, an
+   * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+   */
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const Mapping<dim, spacedim>      &mapping,
+                  const DoFHandler<dim,spacedim>    &dof,
+                  const InVector                    &fe_function,
+                  const Point<spacedim>             &point,
+                  std::vector<Tensor<1, spacedim> > &value);
+
+  /**
+   * Same as above for hp.
+   *
+   * @note If the cell in which the point is found is not locally owned, an
+   * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+   */
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const hp::MappingCollection<dim, spacedim> &mapping,
+                  const hp::DoFHandler<dim,spacedim>         &dof,
+                  const InVector                             &fe_function,
+                  const Point<spacedim>                      &point,
+                  std::vector<Tensor<1, spacedim> >          &value);
+
+  /**
+   * Evaluate a scalar finite element function defined by the given DoFHandler
+   * and nodal vector at the given point, and return the gradient of this
+   * function.
+   *
+   * Compared with the other function of the same name, this function uses an
+   * arbitrary mapping for evaluation.
+   *
+   * @note If the cell in which the point is found is not locally owned, an
+   * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+   */
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const Mapping<dim,spacedim>    &mapping,
+                  const DoFHandler<dim,spacedim> &dof,
+                  const InVector                 &fe_function,
+                  const Point<spacedim>          &point);
+
+  /**
+   * Same as above for hp.
+   *
+   * @note If the cell in which the point is found is not locally owned, an
+   * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
+   */
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const hp::MappingCollection<dim,spacedim> &mapping,
+                  const hp::DoFHandler<dim,spacedim>        &dof,
+                  const InVector                            &fe_function,
+                  const Point<spacedim>                     &point);
+
   //@}
   /**
    * Mean value operations
index e576dcd2eaed73b80442bffcd3e2adc3761a5920..08ecab5eff741c2219113b1ebb192d90778c8341 100644 (file)
@@ -6793,6 +6793,218 @@ namespace VectorTools
 
 
 
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const DoFHandler<dim,spacedim>    &dof,
+                  const InVector                    &fe_function,
+                  const Point<spacedim>             &point,
+                  std::vector<Tensor<1, spacedim> > &gradients)
+  {
+
+    point_gradient (StaticMappingQ1<dim,spacedim>::mapping,
+                    dof,
+                    fe_function,
+                    point,
+                    gradients);
+  }
+
+
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
+                  const InVector        &fe_function,
+                  const Point<spacedim>      &point,
+                  std::vector<Tensor<1, spacedim> > &gradients)
+  {
+    point_gradient(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
+                   dof,
+                   fe_function,
+                   point,
+                   gradients);
+  }
+
+
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const DoFHandler<dim,spacedim> &dof,
+                  const InVector        &fe_function,
+                  const Point<spacedim>      &point)
+  {
+    return point_gradient (StaticMappingQ1<dim,spacedim>::mapping,
+                           dof,
+                           fe_function,
+                           point);
+  }
+
+
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
+                  const InVector        &fe_function,
+                  const Point<spacedim>      &point)
+  {
+    return point_gradient(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
+                          dof,
+                          fe_function,
+                          point);
+  }
+
+
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const Mapping<dim, spacedim>    &mapping,
+                  const DoFHandler<dim,spacedim> &dof,
+                  const InVector        &fe_function,
+                  const Point<spacedim>      &point,
+                  std::vector<Tensor<1, spacedim> > &gradient)
+  {
+    const FiniteElement<dim> &fe = dof.get_fe();
+
+    Assert(gradient.size() == fe.n_components(),
+           ExcDimensionMismatch(gradient.size(), fe.n_components()));
+
+    // first find the cell in which this point
+    // is, initialize a quadrature rule with
+    // it, and then a FEValues object
+    const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+    cell_point
+      = GridTools::find_active_cell_around_point (mapping, dof, point);
+
+    AssertThrow(cell_point.first->is_locally_owned(),
+                ExcPointNotAvailableHere());
+    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
+           ExcInternalError());
+
+    const Quadrature<dim>
+    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+
+    FEValues<dim> fe_values(mapping, fe, quadrature, update_gradients);
+    fe_values.reinit(cell_point.first);
+
+    // then use this to get the gradients of
+    // the given fe_function at this point
+    std::vector<std::vector<Tensor<1, dim> > > u_gradient(1, std::vector<Tensor<1, dim> > (fe.n_components()));
+    fe_values.get_function_gradients(fe_function, u_gradient);
+
+    gradient = u_gradient[0];
+  }
+
+
+  template <int dim, class InVector, int spacedim>
+  void
+  point_gradient (const hp::MappingCollection<dim, spacedim>    &mapping,
+                  const hp::DoFHandler<dim,spacedim> &dof,
+                  const InVector        &fe_function,
+                  const Point<spacedim>      &point,
+                  std::vector<Tensor<1, spacedim> > &gradient)
+  {
+    const hp::FECollection<dim, spacedim> &fe = dof.get_fe();
+
+    Assert(gradient.size() == fe.n_components(),
+           ExcDimensionMismatch(gradient.size(), fe.n_components()));
+
+    // first find the cell in which this point
+    // is, initialize a quadrature rule with
+    // it, and then a FEValues object
+    const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+    cell_point =  GridTools::find_active_cell_around_point (mapping, dof, point);
+
+    AssertThrow(cell_point.first->is_locally_owned(),
+                ExcPointNotAvailableHere());
+    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
+           ExcInternalError());
+
+    const Quadrature<dim>
+    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+    hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe, hp::QCollection<dim>(quadrature), update_gradients);
+    hp_fe_values.reinit(cell_point.first);
+    const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
+
+    // then use this to get the gradients of
+    // the given fe_function at this point
+    std::vector<std::vector<Tensor<1, dim> > > u_gradient(1, std::vector<Tensor<1, dim> > (fe.n_components()));
+    fe_values.get_function_gradients(fe_function, u_gradient);
+
+    gradient = u_gradient[0];
+  }
+
+
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const Mapping<dim, spacedim>    &mapping,
+                  const DoFHandler<dim,spacedim> &dof,
+                  const InVector        &fe_function,
+                  const Point<spacedim>      &point)
+  {
+    const FiniteElement<dim> &fe = dof.get_fe();
+
+    Assert(fe.n_components() == 1,
+           ExcMessage ("Finite element is not scalar as is necessary for this function"));
+
+    // first find the cell in which this point
+    // is, initialize a quadrature rule with
+    // it, and then a FEValues object
+    const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+    cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
+
+    AssertThrow(cell_point.first->is_locally_owned(),
+                ExcPointNotAvailableHere());
+    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
+           ExcInternalError());
+
+    const Quadrature<dim>
+    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+    FEValues<dim> fe_values(mapping, fe, quadrature, update_gradients);
+    fe_values.reinit(cell_point.first);
+
+    // then use this to get the gradients of
+    // the given fe_function at this point
+    std::vector<Tensor<1, dim> > u_gradient(1);
+    fe_values.get_function_gradients(fe_function, u_gradient);
+
+    return u_gradient[0];
+  }
+
+
+  template <int dim, class InVector, int spacedim>
+  Tensor<1, spacedim>
+  point_gradient (const hp::MappingCollection<dim, spacedim>    &mapping,
+                  const hp::DoFHandler<dim,spacedim> &dof,
+                  const InVector        &fe_function,
+                  const Point<spacedim>      &point)
+  {
+    const hp::FECollection<dim, spacedim> &fe = dof.get_fe();
+
+    Assert(fe.n_components() == 1,
+           ExcMessage ("Finite element is not scalar as is necessary for this function"));
+
+    // first find the cell in which this point
+    // is, initialize a quadrature rule with
+    // it, and then a FEValues object
+    const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+    cell_point =  GridTools::find_active_cell_around_point (mapping, dof, point);
+
+    AssertThrow(cell_point.first->is_locally_owned(),
+                ExcPointNotAvailableHere());
+    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
+           ExcInternalError());
+
+    const Quadrature<dim>
+    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+    hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe, hp::QCollection<dim>(quadrature), update_gradients);
+    hp_fe_values.reinit(cell_point.first);
+    const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
+
+    // then use this to get the gradients of
+    // the given fe_function at this point
+    std::vector<Tensor<1, dim > > u_gradient(1);
+    fe_values.get_function_gradients(fe_function, u_gradient);
+
+    return u_gradient[0];
+  }
+
+
+
   template <class VECTOR>
   void
   subtract_mean_value(VECTOR                  &v,
index c1f0cc2868b803d6260ffb11f5ef91ecb365f75e..04e28cdd9fa18338c0bb446ed48d6476edf55ea1 100644 (file)
@@ -41,6 +41,7 @@ SET(_src
   vector_tools_interpolate.cc
   vector_tools_mean_value.cc
   vector_tools_point_value.cc
+  vector_tools_point_gradient.cc
   vector_tools_project.cc
   vector_tools_rhs.cc
   )
@@ -67,6 +68,7 @@ SET(_inst
   vector_tools_interpolate.inst.in
   vector_tools_mean_value.inst.in
   vector_tools_point_value.inst.in
+  vector_tools_point_gradient.inst.in
   vector_tools_project.inst.in
   vector_tools_rhs.inst.in
   )
diff --git a/source/numerics/vector_tools_point_gradient.cc b/source/numerics/vector_tools_point_gradient.cc
new file mode 100644 (file)
index 0000000..dc6867f
--- /dev/null
@@ -0,0 +1,24 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2014 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/numerics/vector_tools.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// ---------------------------- explicit instantiations --------------------
+#include "vector_tools_point_gradient.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/numerics/vector_tools_point_gradient.inst.in b/source/numerics/vector_tools_point_gradient.inst.in
new file mode 100644 (file)
index 0000000..0844b7c
--- /dev/null
@@ -0,0 +1,80 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension == deal_II_space_dimension
+
+    namespace VectorTools \{
+
+  template
+    void point_gradient (
+      const hp::DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&,
+      std::vector<Tensor<1,deal_II_space_dimension> >&);
+
+  template
+    Tensor<1,deal_II_space_dimension> point_gradient (
+      const hp::DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&);
+
+  template
+     void point_gradient (
+      const hp::MappingCollection<deal_II_dimension>&,
+      const hp::DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&,
+      std::vector<Tensor<1,deal_II_space_dimension> >&);
+
+  template
+    Tensor<1,deal_II_space_dimension> point_gradient (
+     const hp::MappingCollection<deal_II_dimension>&,
+      const hp::DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&);
+
+  template
+    void point_gradient (
+      const DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&,
+      std::vector<Tensor<1,deal_II_space_dimension> >&);
+
+  template
+    Tensor<1,deal_II_space_dimension> point_gradient (
+      const DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&);
+
+  template
+    void point_gradient (
+      const Mapping<deal_II_dimension>&,
+      const DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&,
+      std::vector<Tensor<1,deal_II_space_dimension> >&);
+
+  template
+    Tensor<1,deal_II_space_dimension> point_gradient (
+      const Mapping<deal_II_dimension>&,
+      const DoFHandler<deal_II_dimension>&,
+      const VEC&,
+      const Point<deal_II_dimension>&);
+
+      \}
+#endif
+  }
diff --git a/tests/bits/point_gradient_01.cc b/tests/bits/point_gradient_01.cc
new file mode 100644 (file)
index 0000000..7023844
--- /dev/null
@@ -0,0 +1,212 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::point_gradient
+
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <cmath>
+#include <iomanip>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+public:
+  MySquareFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return (component+1)*p.square()+1;
+  }
+
+  virtual void vector_value (const Point<dim>   &p,
+                             Vector<double> &values) const
+  {
+    values[0] = value(p,0);
+  }
+
+
+  virtual Tensor<1, dim> gradient (const Point<dim>   &p,
+                                   const unsigned int  component) const
+  {
+    Tensor<1, dim> return_value;
+    for (unsigned int i=0; i<dim; ++i)
+      return_value[i] = 2*(component+1)*p(component);
+    return return_value;
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim> > &gradients) const
+  {
+    gradients[0] = gradient(p,0);
+  }
+};
+
+
+template<int dim>
+class MyExpFunction : public Function<dim>
+{
+public:
+  MyExpFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_value (const Point<dim>   &p,
+                               Vector<double>     &values) const
+  {
+    values(0) = value(p,0);
+  }
+
+  virtual double gradient (const Point<dim>   &p,
+                           const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim > > &gradients) const
+  {
+    gradients(0) = gradient(p,0);
+  }
+
+};
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+  // refine the mesh in a random way so as to
+  // generate as many cells with
+  // hanging nodes as possible
+  tria.refine_global (4-dim);
+  const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+  for (unsigned int i=0; i<steps[dim]; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+      cell = tria.begin_active();
+      for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+        if (index % (3*dim) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement ();
+    }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tria;
+  make_mesh (tria);
+
+  FE_Q<dim> element(3);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(element);
+
+  // test with two different functions: one
+  // that is exactly representable on the
+  // chosen finite element space, and one
+  // that isn't
+  for (unsigned int i=0; i<2; ++i)
+    {
+      static const MySquareFunction<dim>          function_1;
+      static const Functions::CosineFunction<dim> function_2;
+
+      const Function<dim> &
+      function = (i==0 ?
+                  static_cast<const Function<dim>&>(function_1) :
+                  static_cast<const Function<dim>&>(function_2));
+
+      Vector<double> v (dof.n_dofs());
+      VectorTools::interpolate (dof, function, v);
+
+      // for the following points, check the
+      // function gradient, output it, and
+      // verify that the gradient retrieved from
+      // the interpolated function is close
+      // enough to that of the real function
+      //
+      // also verify that the actual value is
+      // roughly correct
+      Point<dim> p[3];
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          p[0][d] = 0;
+          p[1][d] = 0.5;
+          p[2][d] = 1./3.;
+        }
+      std::vector<Tensor<1, dim> > gradient(1);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          VectorTools::point_gradient (dof, v, p[i], gradient);
+          deallog << -gradient[0] << std::endl;
+
+          Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4,
+                  ExcInternalError());
+
+          const Tensor<1, dim> scalar_gradient = VectorTools::point_gradient (dof, v, p[i]);
+          Assert (std::abs((gradient[0] - scalar_gradient).norm_square()) < 1e-4,
+                  ExcInternalError());
+        }
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (4);
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/bits/point_gradient_01.output b/tests/bits/point_gradient_01.output
new file mode 100644 (file)
index 0000000..da86705
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL:1d::0.000
+DEAL:1d::-1.000
+DEAL:1d::-0.6667
+DEAL:1d::2.150e-07
+DEAL:1d::1.111
+DEAL:1d::0.7854
+DEAL:1d::OK
+DEAL:2d::0.000 0.000
+DEAL:2d::-1.000 -1.000
+DEAL:2d::-0.6667 -0.6667
+DEAL:2d::-0.0008683 0.0008683
+DEAL:2d::0.7855 0.7855
+DEAL:2d::0.6802 0.6802
+DEAL:2d::OK
+DEAL:3d::2.463e-16 -2.463e-16 0.000
+DEAL:3d::-1.000 -1.000 -1.000
+DEAL:3d::-0.6667 -0.6667 -0.6667
+DEAL:3d::-0.0008683 0.0008683 0.0008683
+DEAL:3d::0.5585 0.5585 0.5585
+DEAL:3d::0.5874 0.5874 0.5874
+DEAL:3d::OK
diff --git a/tests/bits/point_gradient_02.cc b/tests/bits/point_gradient_02.cc
new file mode 100644 (file)
index 0000000..0a0e654
--- /dev/null
@@ -0,0 +1,214 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::point_gradient, alternative algorithm with mapping
+
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <cmath>
+#include <iomanip>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+public:
+  MySquareFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return (component+1)*p.square()+1;
+  }
+
+  virtual void vector_value (const Point<dim>   &p,
+                             Vector<double> &values) const
+  {
+    values[0] = value(p,0);
+  }
+
+
+  virtual Tensor<1, dim> gradient (const Point<dim>   &p,
+                                   const unsigned int  component) const
+  {
+    Tensor<1, dim> return_value;
+    for (unsigned int i=0; i<dim; ++i)
+      return_value[i] = 2*(component+1)*p(component);
+    return return_value;
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim> > &gradients) const
+  {
+    gradients[0] = gradient(p,0);
+  }
+};
+
+
+template<int dim>
+class MyExpFunction : public Function<dim>
+{
+public:
+  MyExpFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_value (const Point<dim>   &p,
+                               Vector<double>     &values) const
+  {
+    values(0) = value(p,0);
+  }
+
+  virtual double gradient (const Point<dim>   &p,
+                           const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim > > &gradients) const
+  {
+    gradients(0) = gradient(p,0);
+  }
+
+};
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+  // refine the mesh in a random way so as to
+  // generate as many cells with
+  // hanging nodes as possible
+  tria.refine_global (4-dim);
+  const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+  for (unsigned int i=0; i<steps[dim]; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+      cell = tria.begin_active();
+      for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+        if (index % (3*dim) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement ();
+    }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tria;
+  make_mesh (tria);
+
+  FE_Q<dim> element(3);
+  DoFHandler<dim> dof(tria);
+  MappingQ1<dim> mapping;
+  dof.distribute_dofs(element);
+
+  // test with two different functions: one
+  // that is exactly representable on the
+  // chosen finite element space, and one
+  // that isn't
+  for (unsigned int i=0; i<2; ++i)
+    {
+      static const MySquareFunction<dim>          function_1;
+      static const Functions::CosineFunction<dim> function_2;
+
+      const Function<dim> &
+      function = (i==0 ?
+                  static_cast<const Function<dim>&>(function_1) :
+                  static_cast<const Function<dim>&>(function_2));
+
+      Vector<double> v (dof.n_dofs());
+      VectorTools::interpolate (dof, function, v);
+
+      // for the following points, check the
+      // function gradient, output it, and
+      // verify that the gradient retrieved from
+      // the interpolated function is close
+      // enough to that of the real function
+      //
+      // also verify that the actual value is
+      // roughly correct
+      Point<dim> p[3];
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          p[0][d] = 0;
+          p[1][d] = 0.5;
+          p[2][d] = 1./3.;
+        }
+      std::vector<Tensor<1, dim> > gradient(1);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          VectorTools::point_gradient (mapping, dof, v, p[i], gradient);
+          deallog << -gradient[0] << std::endl;
+
+          Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4,
+                  ExcInternalError());
+
+          const Tensor<1, dim> scalar_gradient = VectorTools::point_gradient (dof, v, p[i]);
+          Assert (std::abs((gradient[0] - scalar_gradient).norm_square()) < 1e-4,
+                  ExcInternalError());
+        }
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (4);
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/bits/point_gradient_02.output b/tests/bits/point_gradient_02.output
new file mode 100644 (file)
index 0000000..da86705
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL:1d::0.000
+DEAL:1d::-1.000
+DEAL:1d::-0.6667
+DEAL:1d::2.150e-07
+DEAL:1d::1.111
+DEAL:1d::0.7854
+DEAL:1d::OK
+DEAL:2d::0.000 0.000
+DEAL:2d::-1.000 -1.000
+DEAL:2d::-0.6667 -0.6667
+DEAL:2d::-0.0008683 0.0008683
+DEAL:2d::0.7855 0.7855
+DEAL:2d::0.6802 0.6802
+DEAL:2d::OK
+DEAL:3d::2.463e-16 -2.463e-16 0.000
+DEAL:3d::-1.000 -1.000 -1.000
+DEAL:3d::-0.6667 -0.6667 -0.6667
+DEAL:3d::-0.0008683 0.0008683 0.0008683
+DEAL:3d::0.5585 0.5585 0.5585
+DEAL:3d::0.5874 0.5874 0.5874
+DEAL:3d::OK
diff --git a/tests/bits/point_gradient_hp_01.cc b/tests/bits/point_gradient_hp_01.cc
new file mode 100644 (file)
index 0000000..4319a3e
--- /dev/null
@@ -0,0 +1,226 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::point_gradient for a hp dofhandler.
+
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <cmath>
+#include <iomanip>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+public:
+  MySquareFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return (component+1)*p.square()+1;
+  }
+
+  virtual void vector_value (const Point<dim>   &p,
+                             Vector<double> &values) const
+  {
+    values[0] = value(p,0);
+  }
+
+
+  virtual Tensor<1, dim> gradient (const Point<dim>   &p,
+                                   const unsigned int  component) const
+  {
+    Tensor<1, dim> return_value;
+    for (unsigned int i=0; i<dim; ++i)
+      return_value[i] = 2*(component+1)*p(component);
+    return return_value;
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim> > &gradients) const
+  {
+    gradients[0] = gradient(p,0);
+  }
+};
+
+
+template<int dim>
+class MyExpFunction : public Function<dim>
+{
+public:
+  MyExpFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_value (const Point<dim>   &p,
+                               Vector<double>     &values) const
+  {
+    values(0) = value(p,0);
+  }
+
+  virtual double gradient (const Point<dim>   &p,
+                           const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim > > &gradients) const
+  {
+    gradients(0) = gradient(p,0);
+  }
+
+};
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+  // refine the mesh in a random way so as to
+  // generate as many cells with
+  // hanging nodes as possible
+  tria.refine_global (4-dim);
+  const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+  for (unsigned int i=0; i<steps[dim]; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+      cell = tria.begin_active();
+      for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+        if (index % (3*dim) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement ();
+    }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tria;
+  make_mesh (tria);
+
+  hp::FECollection<dim> fe;
+  fe.push_back (FE_Q<dim>(3));
+  fe.push_back (FE_Q<dim>(4));
+  fe.push_back (FE_Q<dim>(5));
+
+  hp::DoFHandler<dim> dof_handler (tria);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell->set_active_fe_index (Testing::rand() % fe.size());
+    }
+
+  dof_handler.distribute_dofs (fe);
+
+  // test with two different functions: one
+  // that is exactly representable on the
+  // chosen finite element space, and one
+  // that isn't
+  for (unsigned int i=0; i<2; ++i)
+    {
+      static const MySquareFunction<dim>          function_1;
+      static const Functions::CosineFunction<dim> function_2;
+
+      const Function<dim> &
+      function = (i==0 ?
+                  static_cast<const Function<dim>&>(function_1) :
+                  static_cast<const Function<dim>&>(function_2));
+
+      Vector<double> v (dof_handler.n_dofs());
+      VectorTools::interpolate (dof_handler, function, v);
+
+      // for the following points, check the
+      // function gradient, output it, and
+      // verify that the gradient retrieved from
+      // the interpolated function is close
+      // enough to that of the real function
+      //
+      // also verify that the actual value is
+      // roughly correct
+      Point<dim> p[3];
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          p[0][d] = 0;
+          p[1][d] = 0.5;
+          p[2][d] = 1./3.;
+        }
+      std::vector<Tensor<1, dim> > gradient(1);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          VectorTools::point_gradient (dof_handler, v, p[i], gradient);
+          deallog << -gradient[0] << std::endl;
+
+          Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4,
+                  ExcInternalError());
+
+          const Tensor<1, dim> scalar_gradient = VectorTools::point_gradient (dof_handler, v, p[i]);
+          Assert (std::abs((gradient[0] - scalar_gradient).norm_square()) < 1e-4,
+                  ExcInternalError());
+        }
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (4);
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/bits/point_gradient_hp_01.output b/tests/bits/point_gradient_hp_01.output
new file mode 100644 (file)
index 0000000..d6adfac
--- /dev/null
@@ -0,0 +1,22 @@
+
+DEAL:1d::0.000
+DEAL:1d::-1.000
+DEAL:1d::-0.6667
+DEAL:1d::2.150e-07
+DEAL:1d::1.111
+DEAL:1d::0.7854
+DEAL:1d::OK
+DEAL:2d::-3.553e-15 0.000
+DEAL:2d::-1.000 -1.000
+DEAL:2d::-0.6667 -0.6667
+DEAL:2d::-4.747e-06 4.747e-06
+DEAL:2d::0.7854 0.7854
+DEAL:2d::0.6802 0.6802
+DEAL:2d::OK
+DEAL:3d::-3.031e-15 -5.216e-16 0.000
+DEAL:3d::-1.000 -1.000 -1.000
+DEAL:3d::-0.6667 -0.6667 -0.6667
+DEAL:3d::-4.747e-06 4.747e-06 4.747e-06
+DEAL:3d::0.5553 0.5553 0.5553
+DEAL:3d::0.5890 0.5890 0.5890
+DEAL:3d::OK
diff --git a/tests/bits/point_gradient_hp_02.cc b/tests/bits/point_gradient_hp_02.cc
new file mode 100644 (file)
index 0000000..4bcbe1c
--- /dev/null
@@ -0,0 +1,244 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check VectorTools::point_gradient for hp, alternative algorithm with mapping
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <cmath>
+#include <iomanip>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+public:
+  MySquareFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return (component+1)*p.square()+1;
+  }
+
+  virtual void vector_value (const Point<dim>   &p,
+                             Vector<double> &values) const
+  {
+    values[0] = value(p,0);
+  }
+
+
+  virtual Tensor<1, dim> gradient (const Point<dim>   &p,
+                                   const unsigned int  component) const
+  {
+    Tensor<1, dim> return_value;
+    for (unsigned int i=0; i<dim; ++i)
+      return_value[i] = 2*(component+1)*p(component);
+    return return_value;
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim> > &gradients) const
+  {
+    gradients[0] = gradient(p,0);
+  }
+};
+
+
+template<int dim>
+class MyExpFunction : public Function<dim>
+{
+public:
+  MyExpFunction () : Function<dim> () {}
+
+  virtual double value (const Point<dim>   &p,
+                        const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_value (const Point<dim>   &p,
+                               Vector<double>     &values) const
+  {
+    values(0) = value(p,0);
+  }
+
+  virtual double gradient (const Point<dim>   &p,
+                           const unsigned int  component) const
+  {
+    return std::exp (p(0));
+  }
+
+  virtual void   vector_gradient (const Point<dim>   &p,
+                                  std::vector<Tensor<1, dim > > &gradients) const
+  {
+    gradients(0) = gradient(p,0);
+  }
+
+};
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+
+  GridGenerator::hyper_cube(tria, -1, 1);
+
+  // refine the mesh in a random way so as to
+  // generate as many cells with
+  // hanging nodes as possible
+  tria.refine_global (4-dim);
+  const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+  for (unsigned int i=0; i<steps[dim]; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+      cell = tria.begin_active();
+      for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+        if (index % (3*dim) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement ();
+    }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tria;
+  make_mesh (tria);
+
+  hp::FECollection<dim> fe;
+  fe.push_back (FE_Q<dim>(3));
+  fe.push_back (FE_Q<dim>(4));
+  fe.push_back (FE_Q<dim>(5));
+
+  hp::MappingCollection<dim> mapping_1;
+  mapping_1.push_back(MappingQ1<dim>());
+  mapping_1.push_back(MappingQ1<dim>());
+  mapping_1.push_back(MappingQ1<dim>());
+
+  hp::MappingCollection<dim> mapping_2;
+  mapping_2.push_back(MappingQ1<dim>());
+
+  hp::DoFHandler<dim> dof_handler (tria);
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell->set_active_fe_index (Testing::rand() % fe.size());
+    }
+
+  dof_handler.distribute_dofs (fe);
+
+  // test with two different functions: one
+  // that is exactly representable on the
+  // chosen finite element space, and one
+  // that isn't
+  for (unsigned int i=0; i<2; ++i)
+    {
+      static const MySquareFunction<dim>          function_1;
+      static const Functions::CosineFunction<dim> function_2;
+
+      const Function<dim> &
+      function = (i==0 ?
+                  static_cast<const Function<dim>&>(function_1) :
+                  static_cast<const Function<dim>&>(function_2));
+
+      Vector<double> v (dof_handler.n_dofs());
+      VectorTools::interpolate (dof_handler, function, v);
+
+      // for the following points, check the
+      // function gradient, output it, and
+      // verify that the gradient retrieved from
+      // the interpolated function is close
+      // enough to that of the real function
+      //
+      // also verify that the actual value is
+      // roughly correct
+      Point<dim> p[3];
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          p[0][d] = 0;
+          p[1][d] = 0.5;
+          p[2][d] = 1./3.;
+        }
+      std::vector<Tensor<1, dim> > gradient(1);
+      for (unsigned int i=0; i<3; ++i)
+        {
+          VectorTools::point_gradient (mapping_1, dof_handler, v, p[i], gradient);
+          deallog << -gradient[0] << std::endl;
+
+          Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4,
+                  ExcInternalError());
+          VectorTools::point_gradient (mapping_2, dof_handler, v, p[i], gradient);
+          deallog << -gradient[0] << std::endl;
+
+          Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4,
+                  ExcInternalError());
+
+          const Tensor<1, dim> scalar_gradient_1 = VectorTools::point_gradient (mapping_1, dof_handler, v, p[i]);
+          Assert (std::abs((gradient[0] - scalar_gradient_1).norm_square()) < 1e-4,
+                  ExcInternalError());
+
+          const Tensor<1, dim> scalar_gradient_2 = VectorTools::point_gradient (mapping_2, dof_handler, v, p[i]);
+          Assert (std::abs((gradient[0] - scalar_gradient_2).norm_square()) < 1e-4,
+                  ExcInternalError());
+        }
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (4);
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/bits/point_gradient_hp_02.output b/tests/bits/point_gradient_hp_02.output
new file mode 100644 (file)
index 0000000..600d18b
--- /dev/null
@@ -0,0 +1,40 @@
+
+DEAL:1d::0.000
+DEAL:1d::0.000
+DEAL:1d::-1.000
+DEAL:1d::-1.000
+DEAL:1d::-0.6667
+DEAL:1d::-0.6667
+DEAL:1d::2.150e-07
+DEAL:1d::2.150e-07
+DEAL:1d::1.111
+DEAL:1d::1.111
+DEAL:1d::0.7854
+DEAL:1d::0.7854
+DEAL:1d::OK
+DEAL:2d::-3.553e-15 0.000
+DEAL:2d::-3.553e-15 0.000
+DEAL:2d::-1.000 -1.000
+DEAL:2d::-1.000 -1.000
+DEAL:2d::-0.6667 -0.6667
+DEAL:2d::-0.6667 -0.6667
+DEAL:2d::-4.747e-06 4.747e-06
+DEAL:2d::-4.747e-06 4.747e-06
+DEAL:2d::0.7854 0.7854
+DEAL:2d::0.7854 0.7854
+DEAL:2d::0.6802 0.6802
+DEAL:2d::0.6802 0.6802
+DEAL:2d::OK
+DEAL:3d::-3.031e-15 -5.216e-16 0.000
+DEAL:3d::-3.031e-15 -5.216e-16 0.000
+DEAL:3d::-1.000 -1.000 -1.000
+DEAL:3d::-1.000 -1.000 -1.000
+DEAL:3d::-0.6667 -0.6667 -0.6667
+DEAL:3d::-0.6667 -0.6667 -0.6667
+DEAL:3d::-4.747e-06 4.747e-06 4.747e-06
+DEAL:3d::-4.747e-06 4.747e-06 4.747e-06
+DEAL:3d::0.5553 0.5553 0.5553
+DEAL:3d::0.5553 0.5553 0.5553
+DEAL:3d::0.5890 0.5890 0.5890
+DEAL:3d::0.5890 0.5890 0.5890
+DEAL:3d::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.