From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Wed, 5 Nov 2014 13:02:16 +0000 (+0100)
Subject: Further actions to make MatrixFree/FEEvaluation work in 1D
X-Git-Tag: v8.2.0-rc1~86^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5ae0d4c9e45cf1a114f73f4a67babcee9b1527c7;p=dealii.git

Further actions to make MatrixFree/FEEvaluation work in 1D

- Added test for matrix-vector product
- Made compile with clang.
---

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index c143e2f5ad..9dfa0eabcf 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -1236,12 +1236,11 @@ template <typename Number>
 class FEEvaluationAccess<1,1,Number> : public FEEvaluationBase<1,1,Number>
 {
 public:
-  static const int dim                         = 1;
   typedef Number                                 number_type;
   typedef VectorizedArray<Number>                value_type;
-  typedef Tensor<1,dim,VectorizedArray<Number> > gradient_type;
-  static const unsigned int dimension          = dim;
-  typedef FEEvaluationBase<dim,1,Number>         BaseClass;
+  typedef Tensor<1,1,VectorizedArray<Number> >   gradient_type;
+  static const unsigned int dimension          = 1;
+  typedef FEEvaluationBase<1,1,Number>           BaseClass;
 
   /**
    * Returns the value stored for the local degree of freedom with index @p
@@ -4265,7 +4264,7 @@ FEEvaluationAccess<1,1,Number>
 template <typename Number>
 inline
 FEEvaluationAccess<1,1,Number>
-::FEEvaluationAccess (const FEEvaluationAccess<dim,1,Number> &other)
+::FEEvaluationAccess (const FEEvaluationAccess<1,1,Number> &other)
   :
   FEEvaluationBase <1,1,Number>(other)
 {}
@@ -4437,10 +4436,7 @@ FEEvaluationAccess<1,1,Number>
         this->cell_type == internal::MatrixFreeFunctions::general ?
         this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
 
-      VectorizedArray<Number> new_val = jac[0][0] * grad_in[0];
-      for (unsigned int e=1; e<dim; ++e)
-        new_val += jac[e][0] * grad_in[e];
-      this->gradients_quad[0][0][q_point] = new_val * JxW;
+      this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW;
     }
 }
 
diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index b0c3826589..d2f74ad086 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -191,7 +191,7 @@ namespace internal
               current_data.n_q_points.push_back (n_q_points);
 
               current_data.n_q_points_face.push_back
-              (Utilities::fixed_power<dim-1>(n_q_points_1d[q]));
+              (dim>1 ? Utilities::fixed_power<dim-1>(n_q_points_1d[q]) : 1);
               current_data.quadrature.push_back
               (Quadrature<dim>(quad[my_q][q]));
               current_data.face_quadrature.push_back
diff --git a/tests/matrix_free/matrix_vector_18.cc b/tests/matrix_free/matrix_vector_18.cc
new file mode 100644
index 0000000000..f3c5e596c0
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_18.cc
@@ -0,0 +1,261 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// this tests matrix-free matrix-vector products in 1D (otherwise similar as
+// the other matrix-vector tests)
+
+
+#include "../tests.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/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/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <iostream>
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_operation,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+
+  void local_operation (const MatrixFree<dim,Number> &data,
+                        VECTOR &out,
+                        const VECTOR &in,
+                        const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+    const unsigned int n_q_points = fe_eval.n_q_points;
+
+    Tensor<1,dim,VectorizedArray<Number> > ones;
+    for (unsigned int d=0; d<dim; ++d)
+      ones[d] = Number(1);
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        fe_eval.reinit (cell);
+        fe_eval.read_dof_values (in);
+        fe_eval.evaluate (true, true, true);
+        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) +
+                                     make_vectorized_array<Number>(3.2221)*
+                                     (fe_eval.get_hessian(q)*ones),q);
+          }
+        fe_eval.integrate (true,true);
+        fe_eval.distribute_local_to_global (out);
+      }
+  }
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+              const ConstraintMatrix &constraints,
+              const unsigned int     parallel_option = 0)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  if (parallel_option > 0)
+    deallog << "Parallel option: " << parallel_option << std::endl;
+
+  MatrixFree<dim,number> mf_data;
+  {
+    const QGauss<1> quad (fe_degree+1);
+    typename MatrixFree<dim,number>::AdditionalData data;
+    if (parallel_option == 1)
+      data.tasks_parallel_scheme =
+        MatrixFree<dim,number>::AdditionalData::partition_color;
+    else if (parallel_option == 2)
+      data.tasks_parallel_scheme =
+        MatrixFree<dim,number>::AdditionalData::color;
+    else
+      {
+        Assert (parallel_option == 0, ExcInternalError());
+        data.tasks_parallel_scheme =
+          MatrixFree<dim,number>::AdditionalData::partition_partition;
+      }
+    data.tasks_block_size = 7;
+    data.mapping_update_flags |= update_second_derivatives;
+
+    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 + 3.2221 * \nabla^2 u * ones) + (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+1);
+
+    FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_JxW_values | update_second_derivatives);
+
+    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);
+
+    Tensor<1,dim> ones;
+    for (unsigned int d=0; d<dim; ++d)
+      ones[d] = 1.;
+
+    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) +
+                                       3.2221*
+                                       (fe_values.shape_hessian(j,q_point) * ones))
+                                      +
+                                      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_cube (tria);
+  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_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction<dim>(),
+                                            constraints);
+  constraints.close();
+
+  do_test<dim, fe_degree, double> (dof, constraints);
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (3);
+
+  {
+    deallog.threshold_double(1.e-9);
+    deallog.push("1d");
+    test<1,1>();
+    test<1,2>();
+    test<1,3>();
+    deallog.pop();
+    deallog.push("2d");
+    test<2,1>();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/matrix_free/matrix_vector_18.output b/tests/matrix_free/matrix_vector_18.output
new file mode 100644
index 0000000000..35142b13b8
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_18.output
@@ -0,0 +1,13 @@
+
+DEAL:1d::Testing FE_Q<1>(1)
+DEAL:1d::Norm of difference: 0
+DEAL:1d::
+DEAL:1d::Testing FE_Q<1>(2)
+DEAL:1d::Norm of difference: 0
+DEAL:1d::
+DEAL:1d::Testing FE_Q<1>(3)
+DEAL:1d::Norm of difference: 0
+DEAL:1d::
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::