From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Fri, 28 Feb 2014 16:05:32 +0000 (+0000)
Subject: Fix bug
X-Git-Tag: v8.2.0-rc1~741
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=920f0c1f959e301b748152f29755bfbc7ee84d86;p=dealii.git

Fix bug

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

diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h
index 05f544e929..3d7aed87b9 100644
--- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h
+++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h
@@ -1753,10 +1753,10 @@ namespace internal
 template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
           int n_components_ = 1, typename Number = double >
 class FEEvaluationDGP :
-  public FEEvaluationGeneral<dim,fe_degree,fe_degree+1,n_components_,Number>
+  public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
 {
 public:
-  typedef FEEvaluationGeneral<dim,fe_degree,fe_degree+1,n_components_,Number> BaseClass;
+  typedef FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number> BaseClass;
   typedef Number                            number_type;
   typedef typename BaseClass::value_type    value_type;
   typedef typename BaseClass::gradient_type gradient_type;
@@ -2397,6 +2397,8 @@ FEEvaluationBase<dim,n_components_,Number>
                                                     mapped_geometry->get_cell()->level(),
                                                     mapped_geometry->get_cell()->index(),
                                                     dof_handler);
+      Assert (sizeof(src) == sizeof(VectorType *),
+              ExcMessage("Only one single vector allowed"));
       local_dof_indices.resize(dof_handler->get_fe().dofs_per_cell);
       cell->get_dof_indices(local_dof_indices);
 
@@ -4382,7 +4384,7 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
       // points
       unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
                    proposed_quad_comp = numbers::invalid_unsigned_int;
-      if (dofs_per_cell == this->matrix_info->get_dof_info(fe_no).dofs_per_cell[this->active_fe_index])
+      if (dofs_per_cell == this->data->dofs_per_cell)
         proposed_dof_comp = fe_no;
       else
         for (unsigned int no=0; no<this->matrix_info->n_components(); ++no)
@@ -4476,6 +4478,14 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
   BaseClass (geometry, dof_handler, first_selected_component)
 {
   set_data_pointers();
+
+  Assert ((dofs_per_cell == this->data->dofs_per_cell ||
+           internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value ==
+           this->data->dofs_per_cell
+           )
+          &&
+          n_q_points == this->data->n_q_points,
+          ExcMessage("Underlying element and template arguments do not match"));
 }
 
 
@@ -6625,7 +6635,11 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
                    const unsigned int quad_no)
   :
   BaseClass (data_in, fe_no, quad_no)
-{}
+{
+  // reset values_dofs pointer as it has wider gaps in the base class
+  for (unsigned int c=0; c<n_components_; ++c)
+    this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
+}
 
 
 
@@ -6638,7 +6652,10 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
                 const unsigned int                      first_selected_component)
   :
   BaseClass (geometry, dof_handler, first_selected_component)
-{}
+{
+  for (unsigned int c=0; c<n_components_; ++c)
+    this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
+}
 
 
 
@@ -6649,7 +6666,10 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::FEEvaluationDGP (const FEEvaluationDGP &other)
   :
   BaseClass (other)
-{}
+{
+  for (unsigned int c=0; c<n_components_; ++c)
+    this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
+}
 
 
 
@@ -6666,10 +6686,10 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
           internal::ExcAccessToUninitializedField());
 
   // expand dof_values to tensor product
-  VectorizedArray<Number> data_array[n_components*Utilities::fixed_int_power<fe_degree+1,dim>::value];
+  VectorizedArray<Number> data_array[n_components*BaseClass::dofs_per_cell];
   VectorizedArray<Number> *expanded_dof_values[n_components];
   for (unsigned int c=0; c<n_components; ++c)
-    expanded_dof_values[c] = &data_array[c*Utilities::fixed_int_power<fe_degree+1,dim>::value];
+    expanded_dof_values[c] = &data_array[c*BaseClass::dofs_per_cell];
 
   unsigned int count_p = 0, count_q = 0;
   for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i)
@@ -6677,15 +6697,13 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
       for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j)
         {
           for (unsigned int k=0; k<fe_degree+1-j-i; ++k, ++count_p, ++count_q)
-            {
-              for (unsigned int c=0; c<n_components; ++c)
-                expanded_dof_values[c][count_q] = this->values_dofs[c][count_p];
-            }
+            for (unsigned int c=0; c<n_components; ++c)
+              expanded_dof_values[c][count_q] = this->values_dofs[c][count_p];
           for (unsigned int k=fe_degree+1-j-i; k<fe_degree+1; ++k, ++count_q)
             for (unsigned int c=0; c<n_components; ++c)
               expanded_dof_values[c][count_q] = VectorizedArray<Number>();
         }
-      for (unsigned int j=(dim>1?fe_degree+1-i:1); j<(dim>1?fe_degree+1:1); ++j)
+      for (unsigned int j=fe_degree+1-i; j<fe_degree+1; ++j)
         for (unsigned int k=0; k<fe_degree+1; ++k, ++count_q)
           for (unsigned int c=0; c<n_components; ++c)
             expanded_dof_values[c][count_q] = VectorizedArray<Number>();
@@ -6724,10 +6742,10 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
     Assert (this->gradients_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
 
-  VectorizedArray<Number> data_array[n_components*Utilities::fixed_int_power<fe_degree+1,dim>::value];
+  VectorizedArray<Number> data_array[n_components*BaseClass::dofs_per_cell];
   VectorizedArray<Number> *expanded_dof_values[n_components];
   for (unsigned int c=0; c<n_components; ++c)
-    expanded_dof_values[c] = &data_array[c*Utilities::fixed_int_power<fe_degree+1,dim>::value];
+    expanded_dof_values[c] = &data_array[c*BaseClass::dofs_per_cell];
   internal::do_integrate (*this, expanded_dof_values, this->values_quad,
                           this->gradients_quad, integrate_val, integrate_grad,
                           internal::int2type<dim>());
diff --git a/tests/matrix_free/matrix_vector_14b.cc b/tests/matrix_free/matrix_vector_14b.cc
new file mode 100644
index 0000000000..04e33dd37f
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_14b.cc
@@ -0,0 +1,245 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Same as matrix_vector_14, but using more quadrature points than the FE
+// degree would suggest
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.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/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/lac/vector.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename VECTOR>
+void
+helmholtz_operator_dgp (const MatrixFree<dim,typename VECTOR::value_type>  &data,
+                        VECTOR       &dst,
+                        const VECTOR &src,
+                        const std::pair<unsigned int,unsigned int> &cell_range)
+{
+  typedef typename VECTOR::value_type Number;
+  FEEvaluationDGP<dim,fe_degree,fe_degree+2,1,Number> fe_eval (data);
+  const unsigned int n_q_points = fe_eval.n_q_points;
+
+  for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+    {
+      fe_eval.reinit (cell);
+      fe_eval.read_dof_values (src);
+      fe_eval.evaluate (true, true, false);
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
+          fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
+        }
+      fe_eval.integrate (true,true);
+      fe_eval.distribute_local_to_global (dst);
+    }
+}
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+  typedef VectorizedArray<Number> vector_t;
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    dst = 0;
+    const std_cxx1x::function<void(const MatrixFree<dim,typename VECTOR::value_type> &,
+                                   VECTOR &,
+                                   const VECTOR &,
+                                   const std::pair<unsigned int,unsigned int> &)>
+      wrap = helmholtz_operator_dgp<dim,fe_degree,VECTOR>;
+    data.cell_loop (wrap, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+              const ConstraintMatrix &constraints)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim,number> mf_data;
+  {
+    const QGauss<1> quad (fe_degree+2);
+    typename MatrixFree<dim,number>::AdditionalData data;
+    data.tasks_parallel_scheme =
+      MatrixFree<dim,number>::AdditionalData::partition_color;
+    data.tasks_block_size = 3;
+
+    mf_data.reinit (dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+  Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
+  Vector<number> in_dist (dof.n_dofs());
+  Vector<number> out_dist (in_dist);
+
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      const double entry = Testing::rand()/(double)RAND_MAX;
+      in(i) = entry;
+      in_dist(i) = entry;
+    }
+
+  mf.vmult (out_dist, in_dist);
+
+
+  // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
+  SparsityPattern sparsity;
+  {
+    CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
+    sparsity.copy_from(csp);
+  }
+  SparseMatrix<double> sparse_matrix (sparsity);
+  {
+    QGauss<dim>  quadrature_formula(fe_degree+2);
+
+    FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_JxW_values);
+
+    const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof.begin_active(),
+    endc = dof.end();
+    for (; cell!=endc; ++cell)
+      {
+        cell_matrix = 0;
+        fe_values.reinit (cell);
+
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            {
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+                                      fe_values.shape_grad(j,q_point)
+                                      +
+                                      10. *
+                                      fe_values.shape_value(i,q_point) *
+                                      fe_values.shape_value(j,q_point)) *
+                                     fe_values.JxW(q_point));
+            }
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global (cell_matrix,
+                                                local_dof_indices,
+                                                sparse_matrix);
+      }
+  }
+
+  sparse_matrix.vmult (out, in);
+  out -= out_dist;
+  const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
+
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball (tria);
+  static const HyperBallBoundary<dim> boundary;
+  tria.set_boundary (0, boundary);
+  if (dim < 3 || fe_degree < 2)
+    tria.refine_global(1);
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  typename Triangulation<dim>::active_cell_iterator
+  cell = tria.begin_active (),
+  endc = tria.end();
+  for (; cell!=endc; ++cell)
+    if (cell->center().norm()<1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  FE_DGP<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+
+  do_test<dim, fe_degree, double> (dof, constraints);
+}
+
+
+
+int main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (3);
+
+  {
+    deallog.threshold_double(5.e-11);
+    deallog.push("2d");
+    test<2,1>();
+    test<2,2>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
+
+
diff --git a/tests/matrix_free/matrix_vector_14b.output b/tests/matrix_free/matrix_vector_14b.output
new file mode 100644
index 0000000000..2a6f032686
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_14b.output
@@ -0,0 +1,13 @@
+
+DEAL:2d::Testing FE_DGP<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_DGP<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_DGP<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_DGP<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::