working on flow functions

git-svn-id: https://svn.dealii.org/trunk@14907 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/tests/base/functions_01.cc b/tests/base/functions_01.cc
index 544c95c45caabf12deeec7324bb71adfd589132c..4cf49787d0f8dd5fd5009252031fbff9230bd69d 100644 (file)
--- 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 <base/data_out_base.h>
 #include <base/logstream.h>
 #include <base/quadrature_lib.h>
 #include <base/function_lib.h>
+#include <base/auto_derivative_function.h>
 #include <base/flow_function.h>
 #include <lac/vector.h>
 
@@ -27,13 +28,37 @@
 #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(const Function<dim>& f,
+check_function(const Functions::FlowFunction<dim>& f,
               unsigned int sub,
               std::ostream& out)
 {
+  DerivativeTestFunction<dim> dtest1(f, 1.e-2);
+  DerivativeTestFunction<dim> 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<dim>& f,
        }
   
   std::vector<Vector<double> > values(points.size(), Vector<double>(f.n_components));
+  std::vector<std::vector<double> > values2(f.n_components,
+                                           std::vector<double>(points.size()));
   f.vector_value_list(points, values);
+  f.vector_values(points, values2);
   for (unsigned int i=0;i<values.size();++i)
     for (unsigned int j=0;j<values[i].size();++j)
       {
-//     deallog << points[i] << '\t' << values[i](j) << std::endl;
        patches[0].data(j, i) = values[i](j);
+       Assert (values[i](j) == values2[j][i], ExcInternalError());
+      }
+  
+  deallog << "Gradients ";
+                                  // Compute gradients and difference approximations
+  std::vector<std::vector<Tensor<1,dim> > >
+    gradients(points.size(), std::vector<Tensor<1,dim> >(f.n_components));
+  std::vector<std::vector<Tensor<1,dim> > >
+    gradients1(points.size(), std::vector<Tensor<1,dim> >(f.n_components));
+  std::vector<std::vector<Tensor<1,dim> > >
+    gradients2(points.size(), std::vector<Tensor<1,dim> >(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<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 (d2.norm() < 12.* d1.norm())
+             {
+               deallog << "Gradient error " << k << " " << i
+                       << " " << d1.norm() << " " << d2.norm()
+                       << std::endl;
+//             for (unsigned int i=0;i<f.n_components;++i)
+                 for (unsigned int d=0;d<dim;++d)
+                   deallog
+                     << " " << gradients[k][i][d]
+                     << " " << gradients1[k][i][d]
+                     << std::endl;
+             }
+         }
       }
+  deallog << "tested" << std::endl;
+
+                                  // Check if divergence is zero
+  deallog << "Divergence ";
+  
+  for (unsigned int k=0;k<gradients.size();++k)
+    {
+      double div = 0.;
+      for (unsigned int d=0;d<dim;++d)
+       div += gradients[k][d][d];
+      if (std::fabs(div)> 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<values.size();++i)
+    for (unsigned int j=0;j<values[i].size();++j)
+      {
+       sum += values[i](j)*values[i](j);
+       Assert (values[i](j) == values2[j][i], ExcInternalError());
+      }
+  deallog << "Laplacians " << std::sqrt(sum)/points.size() << std::endl;
   
   std::vector<std::string> names(f.n_components);
   for (unsigned int i=0;i<names.size();++i)
@@ -88,30 +188,93 @@ check_function(const Function<dim>& 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 <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);
+}