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
+ 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,
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
)
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
)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff