From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 31 Jan 2014 03:34:47 +0000 (+0000)
Subject: Add tests for bug report #182. At least the _28 test doesn't currently work and I... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=76a6fa7396757be95f2432b92bd48bd291cb05ea;p=dealii-svn.git

Add tests for bug report #182. At least the _28 test doesn't currently work and I will get to it over the weekend at the latest.

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

diff --git a/tests/deal.II/fe_values_view_26.cc b/tests/deal.II/fe_values_view_26.cc
new file mode 100644
index 0000000000..ec1aedf8f8
--- /dev/null
+++ b/tests/deal.II/fe_values_view_26.cc
@@ -0,0 +1,117 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2007 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// a bit like _25, but test for the curl of a function. there was a
+// bug in get_function_curls
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <fstream>
+
+
+
+Tensor<1,1> curl (const Tensor<2,2> &grads)
+{
+  return Point<1>(grads[1][0] - grads[0][1]);
+}
+
+
+Tensor<1,3> curl (const Tensor<2,3> &grads)
+{
+  return Point<3>(grads[2][1] - grads[1][2],
+		  grads[0][2] - grads[2][0],
+		  grads[1][0] - grads[0][1]);
+}
+
+
+
+template<int dim>
+void test (const Triangulation<dim> &tr,
+           const FiniteElement<dim> &fe)
+{
+  deallog << "FE=" << fe.get_name()
+          << std::endl;
+
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  Vector<double> fe_function(dof.n_dofs());
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    fe_function(i) = i+1;
+
+  const QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values (fe, quadrature,
+                           update_values | update_gradients | update_q_points);
+  fe_values.reinit (dof.begin_active());
+
+  // let the FEValues object compute the
+  // divergences at quadrature points
+  std::vector<typename dealii::internal::CurlType<dim>::type> curls (quadrature.size());
+  std::vector<Tensor<2,dim> > grads (quadrature.size());
+  FEValuesExtractors::Vector extractor(0);
+  fe_values[extractor].get_function_curls (fe_function, curls);
+  fe_values[extractor].get_function_gradients (fe_function, grads);
+
+  // now compare
+  for (unsigned int q=0; q<quadrature.size(); ++q)
+    {
+      deallog << "  curls[q]= " << curls[q] << std::endl
+	      << "  grads[q]= " << grads[q] << std::endl;
+      Assert ((curl(grads[q]) - curls[q]).norm()
+	      <= 1e-10,
+	      ExcInternalError());
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_cube()
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+
+  FESystem<dim> fe (FE_Q<dim>(1),
+                    dim);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (3);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-7);
+
+  test_hyper_cube<2>();
+  test_hyper_cube<3>();
+}
diff --git a/tests/deal.II/fe_values_view_26.output b/tests/deal.II/fe_values_view_26.output
new file mode 100644
index 0000000000..9341da7943
--- /dev/null
+++ b/tests/deal.II/fe_values_view_26.output
@@ -0,0 +1,27 @@
+
+DEAL::FE=FESystem<2>[FE_Q<2>(1)^2]
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::FE=FESystem<3>[FE_Q<3>(1)^3]
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
diff --git a/tests/deal.II/fe_values_view_27.cc b/tests/deal.II/fe_values_view_27.cc
new file mode 100644
index 0000000000..54dbb8534c
--- /dev/null
+++ b/tests/deal.II/fe_values_view_27.cc
@@ -0,0 +1,117 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2007 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// like _26, but for non-primitive elements. the bug was in code only
+// relevant to non-primitive elements
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <fstream>
+
+
+
+Tensor<1,1> curl (const Tensor<2,2> &grads)
+{
+  return Point<1>(grads[1][0] - grads[0][1]);
+}
+
+
+Tensor<1,3> curl (const Tensor<2,3> &grads)
+{
+  return Point<3>(grads[2][1] - grads[1][2],
+		  grads[0][2] - grads[2][0],
+		  grads[1][0] - grads[0][1]);
+}
+
+
+
+template<int dim>
+void test (const Triangulation<dim> &tr,
+           const FiniteElement<dim> &fe)
+{
+  deallog << "FE=" << fe.get_name()
+          << std::endl;
+
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  Vector<double> fe_function(dof.n_dofs());
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    fe_function(i) = i+1;
+
+  const QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values (fe, quadrature,
+                           update_values | update_gradients | update_q_points);
+  fe_values.reinit (dof.begin_active());
+
+  // let the FEValues object compute the
+  // divergences at quadrature points
+  std::vector<typename dealii::internal::CurlType<dim>::type> curls (quadrature.size());
+  std::vector<Tensor<2,dim> > grads (quadrature.size());
+  FEValuesExtractors::Vector extractor(0);
+  fe_values[extractor].get_function_curls (fe_function, curls);
+  fe_values[extractor].get_function_gradients (fe_function, grads);
+
+  // now compare
+  for (unsigned int q=0; q<quadrature.size(); ++q)
+    {
+      deallog << "  curls[q]= " << curls[q] << std::endl
+	      << "  grads[q]= " << grads[q] << std::endl;
+      Assert ((curl(grads[q]) - curls[q]).norm()
+	      <= 1e-10,
+	      ExcInternalError());
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_cube()
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+
+  FESystem<dim> fe (FE_Nedelec<dim>(1),
+                    2);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (3);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-7);
+
+  test_hyper_cube<2>();
+  test_hyper_cube<3>();
+}
diff --git a/tests/deal.II/fe_values_view_27.output b/tests/deal.II/fe_values_view_27.output
new file mode 100644
index 0000000000..0fffb35a02
--- /dev/null
+++ b/tests/deal.II/fe_values_view_27.output
@@ -0,0 +1,27 @@
+
+DEAL::FE=FESystem<2>[FE_Nedelec<2>(1)^2]
+DEAL::  curls[q]= 0.00
+DEAL::  grads[q]= 19.6 1.00 1.00 -10.1
+DEAL::  curls[q]= 26.0
+DEAL::  grads[q]= 19.6 -27.0 -1.00 -2.14
+DEAL::  curls[q]= -30.0
+DEAL::  grads[q]= 27.6 -1.00 -31.0 -10.1
+DEAL::  curls[q]= 4.00
+DEAL::  grads[q]= 27.6 43.0 47.0 -2.14
+DEAL::FE=FESystem<3>[FE_Nedelec<3>(1)^2]
+DEAL::  curls[q]= 0.00 0.00 1.42e-14
+DEAL::  grads[q]= -75.8 0.711 0.711 0.711 -86.9 0.711 0.711 0.711 34.6
+DEAL::  curls[q]= 0.00 -54.0 106.
+DEAL::  grads[q]= -75.8 -107. -54.7 -0.711 -83.6 0.711 -0.711 0.711 38.0
+DEAL::  curls[q]= 16.8 -1.07e-14 -109.
+DEAL::  grads[q]= -72.4 -0.711 0.711 -110. -86.9 -17.6 0.711 -0.711 46.0
+DEAL::  curls[q]= 26.1 -63.2 2.85
+DEAL::  grads[q]= -72.4 113. -63.9 116. -83.6 -26.8 -0.711 -0.711 49.4
+DEAL::  curls[q]= -35.7 71.7 1.10e-14
+DEAL::  grads[q]= -48.4 0.711 -0.711 0.711 -59.6 -0.711 -72.4 -36.4 34.6
+DEAL::  curls[q]= -44.9 30.3 115.
+DEAL::  grads[q]= -48.4 -116. 109. -0.711 -56.2 -0.711 79.2 -45.6 38.0
+DEAL::  curls[q]= -13.2 80.9 -118.
+DEAL::  grads[q]= -45.0 -0.711 -0.711 -119. -59.6 72.3 -81.6 59.2 46.0
+DEAL::  curls[q]= -13.2 30.3 2.85
+DEAL::  grads[q]= -45.0 123. 119. 126. -56.2 81.6 88.4 68.4 49.4
diff --git a/tests/deal.II/fe_values_view_28.cc b/tests/deal.II/fe_values_view_28.cc
new file mode 100644
index 0000000000..74b1736bcd
--- /dev/null
+++ b/tests/deal.II/fe_values_view_28.cc
@@ -0,0 +1,142 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2007 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// like _27, but even simpler. this is based on a testcase proposed by
+// Christoph Heiniger
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+class F : public Function<2>
+{
+public:
+  F() : Function<2>(2) {}
+
+  virtual void vector_value (const Point<2> &p,
+			     Vector<double> &v) const
+    {
+      // make the function equal to (0,x^2)
+      v[0] = 0;
+      v[1] = p[0]*p[0];
+    }
+};
+
+
+
+Tensor<1,1> curl (const Tensor<2,2> &grads)
+{
+  return Point<1>(grads[1][0] - grads[0][1]);
+}
+
+
+Tensor<1,3> curl (const Tensor<2,3> &grads)
+{
+  return Point<3>(grads[2][1] - grads[1][2],
+		  grads[0][2] - grads[2][0],
+		  grads[1][0] - grads[0][1]);
+}
+
+
+
+template<int dim>
+void test (const Triangulation<dim> &tr,
+           const FiniteElement<dim> &fe)
+{
+  deallog << "FE=" << fe.get_name()
+          << std::endl;
+
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+
+  Vector<double> fe_function(dof.n_dofs());
+  // set the elements of the vector in such a way that the function
+  // equals the vector function (0,x^2)
+  ConstraintMatrix cm;
+  cm.close ();
+  VectorTools::project (dof, cm, QGauss<2>(4), F(), fe_function);
+
+  const QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values (fe, quadrature,
+                           update_values | update_gradients | update_q_points);
+  fe_values.reinit (dof.begin_active());
+
+  // let the FEValues object compute the
+  // divergences at quadrature points
+  std::vector<typename dealii::internal::CurlType<dim>::type> curls (quadrature.size());
+  std::vector<Tensor<2,dim> > grads (quadrature.size());
+  FEValuesExtractors::Vector extractor(0);
+  fe_values[extractor].get_function_curls (fe_function, curls);
+  fe_values[extractor].get_function_gradients (fe_function, grads);
+
+  // now compare
+  for (unsigned int q=0; q<quadrature.size(); ++q)
+    {
+      deallog << "  curls[q]= " << curls[q] << std::endl
+	      << "  grads[q]= " << grads[q] << std::endl;
+      Assert ((curl(grads[q]) - curls[q]).norm()
+	      <= 1e-10,
+	      ExcInternalError());
+
+      // we know the function F=(0,x^2) and we chose an element with
+      // high enough degree to exactly represent this function, so the
+      // curl of this function should be
+      //   curl F = d_x F_y - d_y F_x = 2x
+      // verify that this is true
+      Assert (std::fabs(curls[q][0] - 2*fe_values.quadrature_point(q)[0])
+	      <= 1e-10,
+	      ExcInternalError());
+    }
+}
+
+
+
+template<int dim>
+void test_hyper_cube()
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+
+  FE_Nedelec<dim> fe(3);
+  test(tr, fe);
+}
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (3);
+
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+  deallog.threshold_double(1.e-7);
+
+  test_hyper_cube<2>();
+}
diff --git a/tests/deal.II/fe_values_view_28.output b/tests/deal.II/fe_values_view_28.output
new file mode 100644
index 0000000000..9341da7943
--- /dev/null
+++ b/tests/deal.II/fe_values_view_28.output
@@ -0,0 +1,27 @@
+
+DEAL::FE=FESystem<2>[FE_Q<2>(1)^2]
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::  curls[q]= -2.00
+DEAL::  grads[q]= 2.00 4.00 2.00 4.00
+DEAL::FE=FESystem<3>[FE_Q<3>(1)^3]
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0
+DEAL::  curls[q]= -6.00 9.00 -3.00
+DEAL::  grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0