From 35651943d3e35106e12cbb01081a9811c8035f48 Mon Sep 17 00:00:00 2001 From: kanschat Date: Mon, 6 Aug 2007 13:23:36 +0000 Subject: [PATCH] working on flow functions git-svn-id: https://svn.dealii.org/trunk@14907 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/base/functions_01.cc | 177 +++++++++++++++++++++++++++++++++++-- 1 file changed, 170 insertions(+), 7 deletions(-) diff --git a/tests/base/functions_01.cc b/tests/base/functions_01.cc index 544c95c45c..4cf49787d0 100644 --- a/tests/base/functions_01.cc +++ b/tests/base/functions_01.cc @@ -11,13 +11,14 @@ // //----------------------------------------------------------------------------- -// Plot functions in library and check their consistency +// Plot flow functions in library and check their consistency #include "../tests.h" #include #include #include #include +#include #include #include @@ -27,13 +28,37 @@ #include +/** + * A class replacing the implemented derivatives of a function with + * difference quotients. This way, the correctness of the + * implementation can be tested. + */ +template +class DerivativeTestFunction : + public AutoDerivativeFunction +{ + public: + DerivativeTestFunction(const Function&, const double h); + ~DerivativeTestFunction(); + + virtual void vector_value (const Point& points, Vector& value) const; + virtual double value (const Point& 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& func; +}; + template void -check_function(const Function& f, +check_function(const Functions::FlowFunction& f, unsigned int sub, std::ostream& out) { + DerivativeTestFunction dtest1(f, 1.e-2); + DerivativeTestFunction dtest2(f, 2.e-2); // Prepare a vector with a single // patch stretching over the cube // [-1,1]^dim @@ -74,13 +99,88 @@ check_function(const Function& f, } std::vector > values(points.size(), Vector(f.n_components)); + std::vector > values2(f.n_components, + std::vector(points.size())); f.vector_value_list(points, values); + f.vector_values(points, values2); for (unsigned int i=0;i > > + gradients(points.size(), std::vector >(f.n_components)); + std::vector > > + gradients1(points.size(), std::vector >(f.n_components)); + std::vector > > + gradients2(points.size(), std::vector >(f.n_components)); + + f.vector_gradient_list(points, gradients); + dtest1.vector_gradient_list(points, gradients1); + dtest2.vector_gradient_list(points, gradients2); + + // Compare gradients and difference quotients + for (unsigned int k=0;k 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 (d2.norm() < 12.* d1.norm()) + { + deallog << "Gradient error " << k << " " << i + << " " << d1.norm() << " " << d2.norm() + << std::endl; +// for (unsigned int i=0;i 1.e-13) + deallog << "Divergence " << k << " " + << div << std::endl; + } + deallog << "tested" << std::endl; + + f.vector_laplacian_list(points, values); + f.vector_laplacians(points, values2); + double sum = 0.; + for (unsigned int i=0;i names(f.n_components); for (unsigned int i=0;i& f, names[i] = std::string("comp"); } - DataOutBase dout; - DataOutBase::DXFlags flags; - dout.write_dx(patches, names, flags, out); + DataOutBase::DXFlags dxflags; + DataOutBase::GnuplotFlags gflags; + if (dim==2) + dout.write_gnuplot(patches, names, gflags, out); + else + dout.write_dx(patches, names, dxflags, out); } int main() { std::ofstream logfile("functions_01/output"); - + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + if (true) { + deallog << " Functions::StokesCosine<2>" << std::endl; + Functions::StokesCosine<2> f(1.); + check_function(f, 10, logfile); + } + + if (false) + { + deallog << "Functions::StokesLSingularity" << std::endl; Functions::StokesLSingularity f; check_function(f, 4, logfile); } + if (true) { + deallog << "Functions::PoisseuilleFlow<2>" << std::endl; Functions::PoisseuilleFlow<2> f(.8, 10.); check_function(f, 4, logfile); } if (true) { + deallog << "Functions::PoisseuilleFlow<3>" << std::endl; Functions::PoisseuilleFlow<3> f(.8, 10.); check_function(f, 4, logfile); } } + + +template +DerivativeTestFunction::DerivativeTestFunction(const Function& f, + const double h) + : + AutoDerivativeFunction(h, f.n_components), + func(f) +{ + this->set_formula(AutoDerivativeFunction::FourthOrder); +} + + +template +DerivativeTestFunction::~DerivativeTestFunction() +{} + + +template +void +DerivativeTestFunction::vector_value_list ( + const std::vector< Point< dim > > &points, + std::vector< Vector< double > > &values) const +{ + func.vector_value_list(points, values); +} + + +template +void DerivativeTestFunction::vector_value ( + const Point& point, + Vector& value) const +{ + func.vector_value(point, value); +} + + +template +double DerivativeTestFunction::value ( + const Point& point, + const unsigned int comp) const +{ +// std::cerr << '[' << point << '!' << func.value(point, comp) << ']'; + + return func.value(point, comp); +} -- 2.39.5