deallog << std::endl;
}
+/**
+ * A class replacing the implemented derivatives of a function with
+ * difference quotients. This way, the correctness of the
+ * implementation can be tested.
+ */
+template <int dim>
+class DerivativeTestFunction :
+ public AutoDerivativeFunction<dim>
+{
+ public:
+ DerivativeTestFunction(const Function<dim>&, const double h);
+ ~DerivativeTestFunction();
+
+ virtual void vector_value (const Point<dim>& points, Vector<double>& value) const;
+ virtual double value (const Point<dim>& points, const unsigned int component) const;
+ virtual void vector_value_list (const std::vector< Point< dim > > &points,
+ std::vector< Vector< double > > &values) const;
+
+ private:
+ const Function<dim>& func;
+};
+
+
+template <int dim>
+DerivativeTestFunction<dim>::DerivativeTestFunction(const Function<dim>& f,
+ const double h)
+ :
+ AutoDerivativeFunction<dim>(h, f.n_components),
+ func(f)
+{
+ this->set_formula(AutoDerivativeFunction<dim>::FourthOrder);
+}
+
+
+template <int dim>
+DerivativeTestFunction<dim>::~DerivativeTestFunction()
+{}
+
+
+template <int dim>
+void
+DerivativeTestFunction<dim>::vector_value_list (
+ const std::vector< Point< dim > > &points,
+ std::vector< Vector< double > > &values) const
+{
+ func.vector_value_list(points, values);
+}
+
+
+template<int dim>
+void DerivativeTestFunction<dim>::vector_value (
+ const Point<dim>& point,
+ Vector<double>& value) const
+{
+ func.vector_value(point, value);
+}
+
+
+template<int dim>
+double DerivativeTestFunction<dim>::value (
+ const Point<dim>& point,
+ const unsigned int comp) const
+{
+// std::cerr << '[' << point << '!' << func.value(point, comp) << ']';
+
+ return func.value(point, comp);
+}
+
+
+// Check whether the difference quotients converge to the gradient
+template<int dim>
+void
+check_gradient(
+ const Function<dim>& f,
+ unsigned int sub,
+ double threshold = 1./14.)
+{
+ DerivativeTestFunction<dim> dtest1(f, 1.e-2);
+ DerivativeTestFunction<dim> dtest2(f, 2.e-2);
+
+ QMidpoint<1> mid;
+ QIterated<dim> quadrature(mid, sub);
+ const std::vector<Point<dim> >& points = quadrature.get_points();
+
+ std::vector<std::vector<Tensor<1,dim> > >
+ gradients(f.n_components, std::vector<Tensor<1,dim> >(points.size()));
+ std::vector<std::vector<Tensor<1,dim> > >
+ gradients1(f.n_components, std::vector<Tensor<1,dim> >(points.size()));
+ std::vector<std::vector<Tensor<1,dim> > >
+ gradients2(f.n_components, std::vector<Tensor<1,dim> >(points.size()));
+
+ deallog << "gradients vs difference quotients";
+
+ f.vector_gradients(points, gradients);
+ dtest1.vector_gradients(points, gradients1);
+ dtest2.vector_gradients(points, gradients2);
+
+ // Compare gradients and difference quotients
+ for (unsigned int k=0;k<gradients.size();++k)
+ for (unsigned int i=0;i<gradients[k].size();++i)
+ {
+ // Compute difference
+ Tensor<1,dim> d1 = gradients1[k][i] - gradients[k][i];
+ Tensor<1,dim> d2 = gradients2[k][i] - gradients[k][i];
+
+ // If the difference is
+ // already small, we are fine
+ if (d1.norm() > 1.e-13)
+ {
+ // Check for
+ // convergence. For full
+ // 4th order, gradients2
+ // should be 16 times as
+ // large, so let's be a
+ // bit generous
+ if (threshold * d2.norm() < d1.norm())
+ {
+ deallog << "Gradient error: point " << i
+ << " (" << points[i] << " )"
+ << " comp " << k
+// << " norms " << d1.norm() << " " << d2.norm()
+ << std::endl;
+ for (unsigned int d=0;d<dim;++d)
+ deallog
+ << " " << gradients[k][i][d]
+ << " " << gradients1[i][k][d]
+ << std::endl;
+ }
+ }
+ }
+ deallog << std::endl;
+}
+
F f; \
check_function_value_consistency(f, 5); \
check_function_gradient_consistency(f, 5); \
+ check_gradient(f, 5); \
}
CHECK(Functions::CosineFunction<1>);
CHECK(Functions::CosineFunction<2>);
CHECK(Functions::CosineFunction<3>);
+ CHECK(Functions::CosineGradFunction<1>);
+ CHECK(Functions::CosineGradFunction<2>);
+ CHECK(Functions::CosineGradFunction<3>);
CHECK(Functions::ExpFunction<1>);
CHECK(Functions::ExpFunction<2>);
CHECK(Functions::ExpFunction<3>);
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::SquareFunction<2>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::SquareFunction<3>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::CosineFunction<1>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::CosineFunction<2>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::CosineFunction<3>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
+DEAL::Functions::CosineGradFunction<1>
+DEAL::value vs vector value list
+DEAL::value list vs vector value list
+DEAL::gradient vs vector gradient list
+DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
+DEAL::Functions::CosineGradFunction<2>
+DEAL::value vs vector value list
+DEAL::value list vs vector value list
+DEAL::gradient vs vector gradient list
+DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
+DEAL::Functions::CosineGradFunction<3>
+DEAL::value vs vector value list
+DEAL::value list vs vector value list
+DEAL::gradient vs vector gradient list
+DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::ExpFunction<1>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::ExpFunction<2>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
DEAL::Functions::ExpFunction<3>
DEAL::value vs vector value list
DEAL::value list vs vector value list
DEAL::gradient vs vector gradient list
DEAL::gradient list vs vector gradient list
+DEAL::gradients vs difference quotients
#include <fstream>
#include <string>
-
-/**
- * A class replacing the implemented derivatives of a function with
- * difference quotients. This way, the correctness of the
- * implementation can be tested.
- */
-template <int dim>
-class DerivativeTestFunction :
- public AutoDerivativeFunction<dim>
-{
- public:
- DerivativeTestFunction(const Function<dim>&, const double h);
- ~DerivativeTestFunction();
-
- virtual void vector_value (const Point<dim>& points, Vector<double>& value) const;
- virtual double value (const Point<dim>& points, const unsigned int component) const;
- virtual void vector_value_list (const std::vector< Point< dim > > &points,
- std::vector< Vector< double > > &values) const;
-
- private:
- const Function<dim>& func;
-};
-
+#include "functions.h"
template<int dim>
void
}
-template <int dim>
-DerivativeTestFunction<dim>::DerivativeTestFunction(const Function<dim>& f,
- const double h)
- :
- AutoDerivativeFunction<dim>(h, f.n_components),
- func(f)
-{
- this->set_formula(AutoDerivativeFunction<dim>::FourthOrder);
-}
-
-
-template <int dim>
-DerivativeTestFunction<dim>::~DerivativeTestFunction()
-{}
-
-
-template <int dim>
-void
-DerivativeTestFunction<dim>::vector_value_list (
- const std::vector< Point< dim > > &points,
- std::vector< Vector< double > > &values) const
-{
- func.vector_value_list(points, values);
-}
-
-
-template<int dim>
-void DerivativeTestFunction<dim>::vector_value (
- const Point<dim>& point,
- Vector<double>& value) const
-{
- func.vector_value(point, value);
-}
-
-
-template<int dim>
-double DerivativeTestFunction<dim>::value (
- const Point<dim>& point,
- const unsigned int comp) const
-{
-// std::cerr << '[' << point << '!' << func.value(point, comp) << ']';
-
- return func.value(point, comp);
-}
#include <fstream>
#include <string>
-
-/**
- * A class replacing the implemented derivatives of a function with
- * difference quotients. This way, the correctness of the
- * implementation can be tested.
- */
-template <int dim>
-class DerivativeTestFunction :
- public AutoDerivativeFunction<dim>
-{
- public:
- DerivativeTestFunction(const Function<dim>&, const double h);
- ~DerivativeTestFunction();
-
- virtual void vector_value (const Point<dim>& points, Vector<double>& value) const;
- virtual double value (const Point<dim>& points, const unsigned int component) const;
- virtual void vector_value_list (const std::vector< Point< dim > > &points,
- std::vector< Vector< double > > &values) const;
-
- private:
- const Function<dim>& func;
-};
-
-
-
template<int dim>
void
check_function_consistency(
check_function_consistency(fl, flg, 5);
}
-
-template <int dim>
-DerivativeTestFunction<dim>::DerivativeTestFunction(const Function<dim>& f,
- const double h)
- :
- AutoDerivativeFunction<dim>(h, f.n_components),
- func(f)
-{
- this->set_formula(AutoDerivativeFunction<dim>::FourthOrder);
-}
-
-
-template <int dim>
-DerivativeTestFunction<dim>::~DerivativeTestFunction()
-{}
-
-
-template <int dim>
-void
-DerivativeTestFunction<dim>::vector_value_list (
- const std::vector< Point< dim > > &points,
- std::vector< Vector< double > > &values) const
-{
- func.vector_value_list(points, values);
-}
-
-
-template<int dim>
-void DerivativeTestFunction<dim>::vector_value (
- const Point<dim>& point,
- Vector<double>& value) const
-{
- func.vector_value(point, value);
-}
-
-
-template<int dim>
-double DerivativeTestFunction<dim>::value (
- const Point<dim>& point,
- const unsigned int comp) const
-{
-// std::cerr << '[' << point << '!' << func.value(point, comp) << ']';
-
- return func.value(point, comp);
-}