Add test case for large polynomial degree

diff --git a/tests/matrix_free/matrix_vector_large_degree.cc b/tests/matrix_free/matrix_vector_large_degree.cc
new file mode 100644 (file)
index 0000000..b327043
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_large_degree.cc
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test that FEEvaluation can correctly handle large polynomial degrees of 20
+// in 3D. The old implementation of FEEvaluation that allocated all data
+// needed by the kernel on the stack ran into stack overflows in that case
+// since degree 20 needs around 3.5 MB of data. Since we cannot compare with
+// matrix-based results in this case, this test simply checks that the
+// matrix-vector product runs without error
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_mf.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(1);
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim,double> mf_data;
+  {
+    const QGauss<1> quad (fe_degree+2);
+    typename MatrixFree<dim,double>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::none;
+    mf_data.reinit (dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim,fe_degree,double,Vector<double>,fe_degree+2> mf (mf_data);
+  Vector<double> in (dof.n_dofs()), out (dof.n_dofs());
+
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      in(i) = 1.;
+    }
+
+  mf.vmult (out, in);
+  deallog << "Result norm: " << out.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  initlog();
+
+  {
+    deallog.push("2d");
+    test<2,5>();
+    test<2,15>();
+    test<2,35>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,10>();
+    test<3,20>();
+    deallog.pop();
+  }
+}

diff --git a/tests/matrix_free/matrix_vector_large_degree.output b/tests/matrix_free/matrix_vector_large_degree.output
new file mode 100644 (file)
index 0000000..dc5e78d
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_large_degree.output
@@ -0,0 +1,11 @@
+
+DEAL:2d::Testing FE_Q<2>(5)
+DEAL:2d::Result norm: 1.14444
+DEAL:2d::Testing FE_Q<2>(15)
+DEAL:2d::Result norm: 0.399128
+DEAL:2d::Testing FE_Q<2>(35)
+DEAL:2d::Result norm: 0.173861
+DEAL:3d::Testing FE_Q<3>(10)
+DEAL:3d::Result norm: 0.143708
+DEAL:3d::Testing FE_Q<3>(20)
+DEAL:3d::Result norm: 0.0523289