From d058779eea80d4c71db6cd300d0a3fbbd8d530b9 Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 8 Jul 2013 14:01:00 +0000
Subject: [PATCH] Minimal form of the test. This fails with my patches to
 work_stream.h because we step on some other array's memory...

git-svn-id: https://svn.dealii.org/trunk@29953 0785d39b-7218-0410-832d-ea1e28bc413d
---
 tests/base/work_stream_03.cc | 188 +++++++++++------------------------
 1 file changed, 58 insertions(+), 130 deletions(-)

diff --git a/tests/base/work_stream_03.cc b/tests/base/work_stream_03.cc
index 61c094d4ef..cea42e56fd 100644
--- a/tests/base/work_stream_03.cc
+++ b/tests/base/work_stream_03.cc
@@ -20,7 +20,6 @@
 // verifying that it indeed works
 
 #include "../tests.h"
-#include <deal.II/base/function.h>
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/work_stream.h>
@@ -28,29 +27,12 @@
 
 #include <deal.II/grid/tria.h>
 #include <deal.II/dofs/dof_handler.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_refinement.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/fe/fe_abf.h>
-#include <deal.II/fe/fe_dgp.h>
-#include <deal.II/fe/fe_dgp_monomial.h>
-#include <deal.II/fe/fe_dgp_nonparametric.h>
-#include <deal.II/fe/fe_dgq.h>
-#include <deal.II/fe/fe_nedelec.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_q_hierarchical.h>
-#include <deal.II/fe/fe_raviart_thomas.h>
-#include <deal.II/fe/fe_system.h>
-#include <deal.II/numerics/matrix_tools.h>
-#include <deal.II/hp/fe_values.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
 
 #include <fstream>
 #include <vector>
@@ -60,39 +42,29 @@ char logname[] = "work_stream_03/output";
 
 
 template <int dim>
-class F :  public Function<dim>
+double value (const Point<dim> &p)
 {
-  public:
-    virtual double value (const Point<dim> &p,
-                          const unsigned int component) const
-      {
-        Assert ((component == 0) && (this->n_components == 1),
-                ExcInternalError());
-        double val = 0;
-        for (unsigned int d=0; d<dim; ++d)
-          for (unsigned int i=0; i<=1; ++i)
-            val += (d+1)*(i+1)*std::pow (p[d], 1.*i);
-        return val;
-      }
-};
+  double val = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    for (unsigned int i=0; i<=1; ++i)
+      val += std::pow (p[d], 1.*i);
+  return val;
+}
 
 namespace
 {
-  template <int dim,
-           int spacedim>
+  template <int dim>
   struct Scratch
   {
-    Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe,
-             const UpdateFlags         update_flags,
-             const ::dealii::hp::QCollection<dim> &quadrature)
+    Scratch (const FiniteElement<dim> &fe,
+             const Quadrature<dim> &quadrature)
       :
       fe_collection (fe),
       quadrature_collection (quadrature),
       x_fe_values (fe_collection,
                    quadrature_collection,
-                   update_flags),
-      rhs_values(quadrature_collection.max_n_quadrature_points()),
-      update_flags (update_flags)
+                   update_q_points),
+      rhs_values(quadrature_collection.size())
     {}
 
     Scratch (const Scratch &data)
@@ -101,19 +73,16 @@ namespace
       quadrature_collection (data.quadrature_collection),
       x_fe_values (fe_collection,
                    quadrature_collection,
-                   data.update_flags),
-      rhs_values (data.rhs_values),
-      update_flags (data.update_flags)
+                   update_q_points),
+      rhs_values (data.rhs_values)
     {}
 
-    const ::dealii::hp::FECollection<dim,spacedim>      &fe_collection;
-    const ::dealii::hp::QCollection<dim>                &quadrature_collection;
+    const FiniteElement<dim>      &fe_collection;
+    const Quadrature<dim>                &quadrature_collection;
 
-    ::dealii::hp::FEValues<dim,spacedim> x_fe_values;
+    FEValues<dim> x_fe_values;
 
     std::vector<double>                  rhs_values;
-
-    const UpdateFlags update_flags;
   };
 
 
@@ -123,108 +92,67 @@ namespace
   };
 }
 
-template<int dim, int spacedim, typename CellIterator>
+template<int dim, typename CellIterator>
   void
-  mass_assembler(const CellIterator &cell, Scratch<dim, spacedim> &data,
+  mass_assembler(const CellIterator &cell, Scratch<dim> &data,
       CopyData &copy_data)
   {
     data.x_fe_values.reinit(cell);
-    const FEValues<dim, spacedim> &fe_values =
-        data.x_fe_values.get_present_fe_values();
-
-    const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points =
-        fe_values.n_quadrature_points;
-    const FiniteElement<dim, spacedim> &fe = fe_values.get_fe();
-    const unsigned int n_components = fe.n_components();
-
-      {
-        copy_data.cell_rhs.reinit(dofs_per_cell);
-      }
-
-    data.rhs_values.resize(n_q_points);
-    F<dim>().value_list(fe_values.get_quadrature_points(),
-        data.rhs_values);
-
-    for (unsigned int i = 0; i < dofs_per_cell; ++i)
-      {
-        const unsigned int component_i = fe.system_to_component_index(i).first;
-        const double *phi_i = &fe_values.shape_value(i, 0);
-
-        double add_data = 0;
-        for (unsigned int point = 0; point < n_q_points; ++point)
-          add_data += phi_i[point] * fe_values.JxW(point)
-              * data.rhs_values[point];
-        copy_data.cell_rhs(i) = add_data;
-      }
+
+    // this appears to be the key: the following line overwrites some of the memory
+    // in which we store the quadrature point location. if the line is moved down,
+    // the comparison in the if() always succeeds...
+    copy_data.cell_rhs = 0;
+    if (cell->center().distance (data.x_fe_values.quadrature_point(0)) >= 1e-6 *cell->diameter())
+      std::cout << '.' << std::flush;
+    else
+      std::cout << '*' << std::flush;
+
+    copy_data.cell_rhs[0] = value(data.x_fe_values.quadrature_point(0));
   }
 
 
-template <typename VectorType>
 void copy_local_to_global (const CopyData &data,
-                           VectorType *right_hand_side)
+                           double *sum)
 {
-//  std::cout << data.cell_rhs.l2_norm() << ' ';
-  right_hand_side->push_back (data.cell_rhs.l2_norm());
+  *sum += data.cell_rhs[0];
 }
 
 
-
-template <int dim, int spacedim>
-void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
-                         const Quadrature<dim>    &q,
-                         std::vector<double>           &rhs_vector)
+void do_project ()
 {
-  hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
-  hp::QCollection<dim>                q_collection (q);
-  Scratch<dim, spacedim>
-  assembler_data (fe_collection,
-                  update_values |
-                  update_JxW_values | update_quadrature_points,
-                  q_collection);
-  CopyData copy_data;
-  copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-
-  dealii::WorkStream::run (dof.begin_active(),
-                   static_cast<typename DoFHandler<dim>::active_cell_iterator>(dof.end()),
-                   &mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>,
-                   std_cxx1x::bind(&
-                                   copy_local_to_global<std::vector<double> >,
-                                   std_cxx1x::_1, &rhs_vector),
-                   assembler_data,
-                   copy_data,
-                   2*multithread_info.n_default_threads,
-                   1);
-}
+  static const int dim = 3;
 
-
-template <int dim>
-void do_project (const FiniteElement<dim> &fe)
-{
   Triangulation<dim>     triangulation;
   GridGenerator::hyper_cube (triangulation);
   triangulation.refine_global (2);
 
+  FE_Nothing<3> fe;
   DoFHandler<dim>        dof_handler(triangulation);
   dof_handler.distribute_dofs (fe);
 
-  std::vector<double> tmp;
-
-  create_mass_matrix (dof_handler, QGauss<dim>(3),
-      tmp);
-
-  double sum=0;
-  for (unsigned int i=0; i<tmp.size(); ++i)
-    sum += std::fabs(tmp[i]);
-  printf ("\nCheck: %5.13f\n", sum);
-}
-
-
+  QMidpoint<dim> q;
 
-template <int dim>
-void test_no_hanging_nodes (const FiniteElement<dim> &fe)
-{
   for (unsigned int i=0; i<12; ++i)
-    do_project (fe);
+    {
+      std::vector<double> tmp;
+
+      double sum = 0;
+      Scratch<dim> assembler_data (dof_handler.get_fe(), q);
+      CopyData copy_data;
+      copy_data.cell_rhs.reinit (8);
+      dealii::WorkStream::run (dof_handler.begin_active(),
+                       dof_handler.end(),
+                       &mass_assembler<dim, typename DoFHandler<dim>::active_cell_iterator>,
+                       std_cxx1x::bind(&
+                                       copy_local_to_global,
+                                       std_cxx1x::_1, &sum),
+                       assembler_data,
+                       copy_data,
+                       8,
+                       1);
+      printf ("\nCheck: %5.3f\n", sum);
+    }
 }
 
 
@@ -237,5 +165,5 @@ int main ()
   deallog.depth_console(0);
   deallog.threshold_double(1.e-10);
 
-  test_no_hanging_nodes (FE_Q<3>(1));
+  do_project ();
 }
-- 
2.39.5