From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 6 Aug 2007 13:23:36 +0000 (+0000)
Subject: working on flow functions
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=35651943d3e35106e12cbb01081a9811c8035f48;p=dealii-svn.git

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 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 <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);
+}