]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add some new tests.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Apr 2003 16:21:57 +0000 (16:21 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Apr 2003 16:21:57 +0000 (16:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@7423 0785d39b-7218-0410-832d-ea1e28bc413d

tests/base/Makefile
tests/base/anisotropic_1.cc [new file with mode: 0644]
tests/base/anisotropic_2.cc [new file with mode: 0644]
tests/base/hierarchical.cc [new file with mode: 0644]

index 2747de2d9afce883e87b0c73218571e9f91f3ef3..dc044d762737ad62017436d958e6bfb9bdc5ede4 100644 (file)
@@ -34,12 +34,15 @@ tensor.exe         : tensor.g.$(OBJEXT)                             $(libraries)
 timer.exe          : timer.g.$(OBJEXT)                              $(libraries)
 threads.exe        : threads.g.$(OBJEXT)                            $(libraries)
 auto_derivative_function.exe : auto_derivative_function.g.$(OBJEXT) $(libraries)
-
+anisotropic_1.exe  : anisotropic_1.g.$(OBJEXT)                      $(lib-base.g)
+anisotropic_2.exe  : anisotropic_2.g.$(OBJEXT)                      $(lib-base.g)
+hierarchical.exe   : hierarchical.g.$(OBJEXT)                       $(lib-base.g)
 
 tests = $(sort \
           logtest reference quadrature_test table tensor \
           timer threads polynomial1d polynomial_test \
-         auto_derivative_function)
+         auto_derivative_function anisotropic_1 anisotropic_2 \
+          hierarchical)
 
 ############################################################
 
diff --git a/tests/base/anisotropic_1.cc b/tests/base/anisotropic_1.cc
new file mode 100644 (file)
index 0000000..9c3f727
--- /dev/null
@@ -0,0 +1,174 @@
+//----------------------------  anisotropic_1.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  anisotropic_1.cc  ---------------------------
+
+
+// check that TensorProducPolynomials and AnisotropicPolynomials
+// compute the same thing if the basis polynomials for the anisotropic
+// class happen to be the same for all coordinate directions
+
+
+#include <fstream>
+#include <cmath>
+
+#include <base/logstream.h>
+#include <base/tensor_product_polynomials.h>
+
+
+using namespace Polynomials;
+
+
+template<int dim, class POLY1, class POLY2>
+void check_poly(const Point<dim> &x,
+               const POLY1      &p,
+                const POLY2      &q)
+{
+  const unsigned int n = p.n();
+  
+  std::vector<double> values1 (n), values2 (n);
+  std::vector<Tensor<1,dim> > gradients1(n), gradients2(n);
+  std::vector<Tensor<2,dim> > second1(n), second2(n);
+  
+  p.compute (x, values1, gradients1, second1);
+  q.compute (x, values2, gradients2, second2);
+  
+  for (unsigned int k=0; k<n; ++k)
+    {
+                                      // Check if compute_value is ok
+      double val1 = p.compute_value(k,x);
+      if (fabs(val1 - values1[k]) > 5.0E-16)
+       deallog << 'P' << k << ": values differ " << val1 << " != "
+               << values1[k] << std::endl;
+      double val2 = q.compute_value(k,x);
+      if (fabs(val2 - values2[k]) > 5.0E-16)
+       deallog << 'Q' << k << ": values differ " << val2 << " != "
+               << values2[k] << std::endl;
+      if (fabs(val2 - val1) > 5.0E-16)
+       deallog << "PQ" << k << ": values differ " << val1 << " != "
+               << val2 << std::endl;
+
+                                       // Check if compute_grad is ok
+      Tensor<1,dim> grad1 = p.compute_grad(k,x);
+      if (grad1 != gradients1[k])
+       deallog << 'P' << k << ": gradients differ " << grad1 << " != "
+               << gradients1[k] << std::endl;
+      Tensor<1,dim> grad2 = q.compute_grad(k,x);
+      if (grad2 != gradients2[k])
+       deallog << 'Q' << k << ": gradients differ " << grad1 << " != "
+               << gradients2[k] << std::endl;
+      if (grad2 != grad1)
+       deallog << "PQ" << k << ": gradients differ " << grad1 << " != "
+               << grad2 << std::endl;
+      
+                                      // Check if compute_grad_grad is ok
+      Tensor<2,dim> grad_grad1 = p.compute_grad_grad(k,x);
+      if (grad_grad1 != second1[k])
+       deallog << 'P' << k << ": second derivatives differ " << grad_grad1 << " != "
+               << second1[k] << std::endl;
+      Tensor<2,dim> grad_grad2 = q.compute_grad_grad(k,x);
+      if (grad_grad2 != second2[k])
+       deallog << 'Q' << k << ": second derivatives differ " << grad_grad2 << " != "
+               << second2[k] << std::endl;
+      if (grad_grad2 != grad_grad1)
+       deallog << "PQ" << k << ": second derivatives differ " << grad_grad1 << " != "
+               << grad_grad2 << std::endl;
+
+
+                                       // finally output values,
+                                       // gradients, etc, to make sure
+                                       // that they are not only
+                                       // consistent, but also
+                                       // correct. Multiply them
+                                       // somewhat to make them
+                                       // significant despite our
+                                       // two-post-dot-digits limit
+      values1[k] *= pow(10, dim);
+      gradients1[k] *= pow(10, dim);
+      
+      deallog << 'P' << k << "\t= " << values1[k]
+             << "\tgradient\t";
+      for (unsigned int d=0;d<dim;++d)
+       deallog << gradients1[k][d] << '\t';
+      deallog << "\t2nd\t";
+      for (unsigned int d1=0;d1<dim;++d1)
+       for (unsigned int d2=0;d2<dim;++d2)
+         deallog << second1[k][d1][d2] << '\t';
+      deallog << std::endl;
+    }
+  deallog << std::endl;
+}
+
+
+template <int dim>
+void
+check_tensor (const std::vector<Polynomial<double> >& v,
+             const Point<dim>& x)
+{
+  TensorProductPolynomials<dim> p(v);
+
+  std::vector<std::vector<Polynomial<double> > > pols (dim, v);
+  AnisotropicPolynomials<dim> q(pols);
+  
+  check_poly (x, p, q);
+}
+
+
+void
+check_dimensions (const std::vector<Polynomial<double> >& p)
+{
+  deallog.push("1d");
+  check_tensor(p, Point<1>(.5));
+  deallog.pop();
+
+  deallog.push("2d");
+  check_tensor(p, Point<2>(.5, .2));
+  deallog.pop();
+
+  deallog.push("3d");
+  check_tensor(p, Point<3>(.5, .2, .3));
+  deallog.pop();
+}
+
+int main()
+{
+  std::ofstream logfile("anisotropic_1.output");
+  logfile.precision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog.push("Lagrange");
+  std::vector<Polynomial<double> > p;
+  for (unsigned int i=0;i<3;++i)
+    p.push_back (LagrangeEquidistant(3, i));
+
+  check_dimensions(p);
+
+  deallog.pop();
+  deallog.push("Legendre");
+  
+  p.clear ();
+  for (unsigned int i=0;i<3;++i)
+    p.push_back (Legendre<double>(i));
+
+  check_dimensions(p);
+
+  deallog.pop();
+  deallog.push("Hierarchical");
+
+  p.clear ();
+  for (unsigned int i=0;i<3;++i)
+    p.push_back (Hierarchical<double>(i));
+
+  check_dimensions(p);
+
+  deallog.pop();
+}
diff --git a/tests/base/anisotropic_2.cc b/tests/base/anisotropic_2.cc
new file mode 100644 (file)
index 0000000..6fd4d51
--- /dev/null
@@ -0,0 +1,166 @@
+//----------------------------  anisotropic_2.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  anisotropic_2.cc  ---------------------------
+
+
+// check AnisotropicPolynomials
+
+
+#include <fstream>
+#include <cmath>
+
+#include <base/logstream.h>
+#include <base/tensor_product_polynomials.h>
+
+
+using namespace Polynomials;
+
+typedef std::vector<Polynomial<double> > PolVector;
+
+
+void print_2d (const AnisotropicPolynomials<2> &aniso)
+{
+  const unsigned int N=10, M=13;
+  for (unsigned int i=0; i<=N; ++i)
+    {
+      deallog << std::endl;
+      for (unsigned int j=0; j<=M; ++j)
+        {
+          deallog << 1.*i/N << " "
+                  << 1.*j/M << " ";
+          for (unsigned int k=0; k<aniso.n(); ++k)
+            deallog << aniso.compute_value (k, Point<2>(1.*i/N, 1.*j/M))
+                    << " ";
+          deallog << std::endl;
+          for (unsigned int k=0; k<aniso.n(); ++k)
+            deallog << aniso.compute_grad (k, Point<2>(1.*i/N, 1.*j/M))
+                    << " ";
+          deallog << std::endl;
+        }
+    }
+}
+
+
+
+template <class Pol>
+void check_2d () 
+{
+                                   // two checks with higher degree in
+                                   // x or y direction
+  {
+    PolVector pols[2] = { Pol::generate_complete_basis (3),
+                          Pol::generate_complete_basis (1) };
+    AnisotropicPolynomials<2> aniso (std::vector<PolVector>(&pols[0],
+                                                            &pols[2]));
+    print_2d (aniso);
+  }
+  {
+    PolVector pols[2] = { Pol::generate_complete_basis (2),
+                          Pol::generate_complete_basis (3) };
+    AnisotropicPolynomials<2> aniso (std::vector<PolVector>(&pols[0],
+                                                            &pols[2]));
+
+    print_2d (aniso);
+  }
+}
+
+
+void print_3d (const AnisotropicPolynomials<3> &aniso)
+{
+  const unsigned int N=4, M=3, P=5;
+  for (unsigned int i=0; i<=N; ++i)
+    {
+      deallog << std::endl;
+      for (unsigned int j=0; j<=M; ++j)
+        {
+          deallog << std::endl;
+          for (unsigned int k=0; k<=P; ++k)
+            {
+              deallog << 1.*i/N << " "
+                      << 1.*j/M << " "
+                      << 1.*k/P;
+              for (unsigned int k=0; k<aniso.n(); ++k)
+                deallog << aniso.compute_value (k, Point<3>(1.*i/N, 1.*j/M, 1.*k/P))
+                        << " ";
+              deallog << std::endl;
+              for (unsigned int k=0; k<aniso.n(); ++k)
+                deallog << aniso.compute_grad (k, Point<3>(1.*i/N, 1.*j/M, 1.*k/P))
+                        << " ";
+              deallog << std::endl;
+            }
+        }
+    }
+}
+
+
+
+template <class Pol>
+void check_3d () 
+{
+                                   // three checks with higher degree
+                                   // in x, y or z direction
+  {
+    PolVector pols[3] = { Pol::generate_complete_basis (3),
+                          Pol::generate_complete_basis (1),
+                          Pol::generate_complete_basis (1) };
+    AnisotropicPolynomials<3> aniso (std::vector<PolVector>(&pols[0],
+                                                            &pols[3]));
+    print_3d (aniso);
+  }
+  {
+    PolVector pols[3] = { Pol::generate_complete_basis (1),
+                          Pol::generate_complete_basis (3),
+                          Pol::generate_complete_basis (1) };
+    AnisotropicPolynomials<3> aniso (std::vector<PolVector>(&pols[0],
+                                                            &pols[3]));
+    print_3d (aniso);
+  }
+  {
+    PolVector pols[3] = { Pol::generate_complete_basis (1),
+                          Pol::generate_complete_basis (2),
+                          Pol::generate_complete_basis (3) };
+    AnisotropicPolynomials<3> aniso (std::vector<PolVector>(&pols[0],
+                                                            &pols[3]));
+    print_3d (aniso);
+  }
+}
+
+
+
+template <class Pol>
+void check () 
+{
+  check_2d<Pol> ();
+  check_3d<Pol> ();
+}
+
+
+
+int main()
+{
+  std::ofstream logfile("anisotropic_2.output");
+  logfile.precision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+//   deallog.push("Lagrange");
+//   check<LagrangeEquidistant> ();
+//   deallog.pop();
+
+//   deallog.push("Legendre");
+//   check<Legendre<double> > ();
+//   deallog.pop();
+
+  deallog.push("Hierarchical");
+  check<Hierarchical<double> > ();
+  deallog.pop();
+}
diff --git a/tests/base/hierarchical.cc b/tests/base/hierarchical.cc
new file mode 100644 (file)
index 0000000..fdab924
--- /dev/null
@@ -0,0 +1,48 @@
+//-----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------------------------
+
+// generate a hierarchical basis and display the values of the shape
+// functions at equidistant points. (I needed this output at one point
+// in time, so why not make it a testcase -- WB)
+
+#include <fstream>
+#include <cmath>
+
+#include <base/logstream.h>
+#include <base/polynomial.h>
+
+
+using namespace Polynomials;
+
+
+int main ()
+{
+  std::ofstream logfile("hierarchical.output");
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  const std::vector<Polynomial<double> >
+    p = Hierarchical<double>::generate_complete_basis (10);
+
+  const unsigned int div=30;
+  for (unsigned int i=0; i<=div; ++i)
+    {
+      const double x = 1.*i/div;
+      deallog << x << " ";
+      for (unsigned int j=0; j<p.size(); ++j)
+        deallog << p[j].value(x) << " ";
+      deallog << std::endl;
+    }
+}
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.