From: kormann <kormann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 24 Jun 2012 14:36:54 +0000 (+0000)
Subject: Fix bug pointed out by Saul Teukolsky and add submit_divergence function.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=22a562162ac6ed0501430c2c22fdb91c59b629d5;p=dealii-svn.git

Fix bug pointed out by Saul Teukolsky and add submit_divergence function.

git-svn-id: https://svn.dealii.org/trunk@25639 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 580042ef36..f68e429ea4 100644
--- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h
+++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h
@@ -429,14 +429,16 @@ public:
   gradient_type get_gradient (const unsigned int q_point) const;
 
                                 /**
-                                 * Write a gradient to the field containing
-                                 * the values on quadrature points with
-                                 * component @p q_point. Access to the same
-                                 * field as through @p get_gradient. If
+                                 * Write a contribution that is tested
+                                 * by the gradient to the field
+                                 * containing the values on quadrature
+                                 * points with component @p
+                                 * q_point. Access to the same field
+                                 * as through @p get_gradient. If
                                  * applied before the function @p
-                                 * integrate(...,true) is called,
-                                 * this specifies the gradient which is tested
-                                 * by all basis function gradients on the
+                                 * integrate(...,true) is called, this
+                                 * specifies what is tested by all
+                                 * basis function gradients on the
                                  * current cell and integrated over.
                                  *
                                  * Note that the derived class
@@ -1081,17 +1083,17 @@ class FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,1,Number> :
   gradient_type get_gradient (const unsigned int q_point) const;
 
                                 /**
-                                 * Write a gradient to the field
+                                 * Write a contribution that is tested
+                                 * by the gradient to the field
                                  * containing the values on quadrature
                                  * points with component @p
                                  * q_point. Access to the same field
                                  * as through @p get_gradient. If
                                  * applied before the function @p
                                  * integrate(...,true) is called, this
-                                 * specifies the gradient which is
-                                 * tested by all basis function
-                                 * gradients on the current cell and
-                                 * integrated over.
+                                 * specifies what is tested by all
+                                 * basis function gradients on the
+                                 * current cell and integrated over.
                                  */
   void submit_gradient(const gradient_type grad_in,
                        const unsigned int  q_point);
@@ -1241,45 +1243,70 @@ class FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number> :
   gradient_type get_hessian_diagonal (const unsigned int q_point) const;
 
                                 /**
-                                 * Write a gradient to the field containing
-                                 * the values on quadrature points with
-                                 * component @p q_point. Access to the same
-                                 * field as through @p get_gradient. If
+                                 * Write a contribution that is tested
+                                 * by the gradient to the field
+                                 * containing the values on quadrature
+                                 * points with component @p
+                                 * q_point. Access to the same field
+                                 * as through @p get_gradient. If
                                  * applied before the function @p
-                                 * integrate(...,true) is called,
-                                 * this specifies the gradient which is tested
-                                 * by all basis function gradients on the
+                                 * integrate(...,true) is called, this
+                                 * specifies what is tested by all
+                                 * basis function gradients on the
                                  * current cell and integrated over.
                                  */
   void submit_gradient(const gradient_type grad_in,
                        const unsigned int  q_point);
 
                                 /**
-                                 * Write a gradient to the field containing
-                                 * the values on quadrature points with
-                                 * component @p q_point. This function is an
-                                 * alternative to the other submit_gradient
-                                 * function when using a system of fixed
-                                 * number of equations which happens to
-                                 * coincide with the dimension for some
-                                 * dimensions, but not all. To allow for
-                                 * dimension-independent programming, this
-                                 * function can be used instead.
+                                 * Write a contribution that is tested
+                                 * by the gradient to the field
+                                 * containing the values on quadrature
+                                 * points with component @p
+                                 * q_point. This function is an
+                                 * alternative to the other
+                                 * submit_gradient function when using
+                                 * a system of fixed number of
+                                 * equations which happens to coincide
+                                 * with the dimension for some
+                                 * dimensions, but not all. To allow
+                                 * for dimension-independent
+                                 * programming, this function can be
+                                 * used instead.
                                  */
   void submit_gradient(const Tensor<1,dim,Tensor<1,dim,VectorizedArray<Number> > > grad_in,
                        const unsigned int q_point);
 
                                 /**
-                                 * Write a gradient to the field containing
-                                 * the values on quadrature points with
-                                 * component @p q_point. Access to the same
-                                 * field as through @p get_gradient. If
+                                 * Write a constribution that is
+                                 * tested by the divergence to the field
+                                 * containing the values on quadrature
+                                 * points with component @p
+                                 * q_point. Access to the same field
+                                 * as through @p get_gradient. If
                                  * applied before the function @p
-                                 * integrate(...,true) is called,
-                                 * this specifies the gradient which is tested
-                                 * by all basis function gradients on the
+                                 * integrate(...,true) is called, this
+                                 * specifies what is tested by all
+                                 * basis function gradients on the
                                  * current cell and integrated over.
                                  */
+  void submit_divergence (const VectorizedArray<Number> div_in,
+                          const unsigned int q_point);
+
+                                /**
+                                 * Write a contribution that is tested
+                                 * by the gradient to the field
+                                 * containing the values on quadrature
+                                 * points with component @p
+                                 * q_point. Access to the same field
+                                 * as through @p get_gradient. If
+                                 * applied before the function @p
+                                 * integrate(...,true) is called, this
+                                 * specifies the gradient which is
+                                 * tested by all basis function
+                                 * gradients on the current cell and
+                                 * integrated over.
+                                 */
   void submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray<Number> > grad_in,
                                  const unsigned int      q_point);
 
@@ -3878,6 +3905,51 @@ FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
   BaseClass::submit_gradient(grad_in, q_point);
 }
 
+template <int dim, int dofs_per_cell_, int n_q_points_,
+          typename Number>
+inline
+void
+FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
+::submit_divergence (const VectorizedArray<Number> div_in,
+                     const unsigned int q_point)
+{
+#ifdef DEBUG
+  Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  AssertIndexRange (q_point, n_q_points);
+  this->gradients_quad_submitted = true;
+#endif
+  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
+    {
+      const VectorizedArray<Number> fac = this->J_value[0] *
+        this->quadrature_weights[q_point] * div_in;
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          this->gradients_quad[d][d][q_point] = (fac *
+                                                 this->cartesian_data[0][d]);
+          for(unsigned int e=d+1;e<dim;++e)
+            {
+              this->gradients_quad[d][e][q_point] = VectorizedArray<Number>();
+              this->gradients_quad[e][d][q_point] = VectorizedArray<Number>();
+            }
+        }
+    }
+  else
+    {
+      const Tensor<2,dim,VectorizedArray<Number> > &jac =
+        this->cell_type == internal::MatrixFreeFunctions::general ?
+        this->jacobian[q_point] : this->jacobian[0];
+      const VectorizedArray<Number> fac =
+        (this->cell_type == internal::MatrixFreeFunctions::general ?
+         this->J_value[q_point] : this->J_value[0] * 
+         this->quadrature_weights[q_point]) * div_in;
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          for (unsigned e=0; e<dim; ++e)
+            this->gradients_quad[d][e][q_point] = jac[d][e] * fac;
+        }
+    }
+}
+
 
 
 template <int dim, int dofs_per_cell_,  int n_q_points_, typename Number>
diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.h b/deal.II/include/deal.II/matrix_free/matrix_free.h
index bacc960834..864307f03f 100644
--- a/deal.II/include/deal.II/matrix_free/matrix_free.h
+++ b/deal.II/include/deal.II/matrix_free/matrix_free.h
@@ -1218,7 +1218,7 @@ MatrixFree<dim,Number>::initialize_dof_vector(VectorType &vec,
                                                   const unsigned int comp) const
 {
   AssertIndexRange (comp, n_components());
-  vec.reinit(dof_info[comp].vector_partitioner->global_size);
+  vec.reinit(dof_info[comp].vector_partitioner->size());
 }
 
 
diff --git a/tests/matrix_free/matrix_vector_div.cc b/tests/matrix_free/matrix_vector_div.cc
new file mode 100644
index 0000000000..4ecf10ca7d
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_div.cc
@@ -0,0 +1,289 @@
+//------------------  matrix_vector_div.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//------------------  matrix_vector_div.cc  ------------------------
+
+
+// this function tests the correctness of the implementation of matrix free
+// matrix-vector products by comparing with the result of deal.II sparse
+// matrix. The mesh uses a hyperball mesh with hanging nodes for a
+// vector-valued problem (div-div operator which does not really make a lot
+// of sense from a problem point of view, though).
+
+#include "../tests.h"
+
+std::ofstream logfile("matrix_vector_div/output");
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_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/dofs/dof_renumbering.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.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/numerics/vectors.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+#include "create_mesh.h"
+
+const double global_coefficient = 0.1;
+
+
+template <int dim, int degree, typename VectorType>
+class MatrixFreeTest
+{
+ public:
+  typedef typename DoFHandler<dim>::active_cell_iterator CellIterator;
+  typedef double Number;
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void
+  local_apply (const MatrixFree<dim,Number> &data,
+               VectorType          &dst,
+               const VectorType    &src,
+               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    typedef VectorizedArray<Number> vector_t;
+    FEEvaluation<dim,degree,degree+1,dim,Number> phi (data);
+    vector_t coeff = make_vectorized_array(global_coefficient);
+
+    for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
+      {
+        phi.reinit (cell);
+        phi.read_dof_values (src);
+        phi.evaluate (false,true,false);
+
+        for (unsigned int q=0; q<phi.n_q_points; ++q)
+          phi.submit_divergence (coeff * phi.get_divergence(q), q);
+
+        phi.integrate (false,true);
+        phi.distribute_local_to_global (dst);
+      }
+  }
+
+  void vmult (VectorType &dst,
+              const VectorType &src) const
+  {
+    AssertDimension (dst.size(), dim);
+    for (unsigned int d=0; d<dim; ++d)
+      dst[d] = 0;
+    data.cell_loop (&MatrixFreeTest<dim,degree,VectorType>::local_apply,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim>   tria;
+  create_mesh (tria);
+  tria.refine_global(4-dim);
+
+                                // refine a few cells
+  for (unsigned int i=0; i<10-3*dim; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator
+        cell = tria.begin_active (),
+        endc = tria.end();
+      unsigned int counter = 0;
+      for (; cell!=endc; ++cell, ++counter)
+        if (counter % (7-i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+
+  FE_Q<dim>            fe_sca (QGaussLobatto<1>(fe_degree+1));
+  FESystem<dim>        fe (fe_sca, dim);
+  DoFHandler<dim>      dof_handler_sca (tria);
+  DoFHandler<dim>      dof_handler (tria);
+
+  MatrixFree<dim,double> mf_data;
+
+  ConstraintMatrix     constraints;
+
+  BlockSparsityPattern      sparsity_pattern;
+  BlockSparseMatrix<double> system_matrix;
+
+  BlockVector<double> solution;
+  BlockVector<double> system_rhs;
+  std::vector<Vector<double> > vec1, vec2;
+
+  dof_handler.distribute_dofs (fe);
+  dof_handler_sca.distribute_dofs (fe_sca);
+  DoFRenumbering::component_wise (dof_handler);
+
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+  constraints.close ();
+
+  const unsigned int dofs_per_block = dof_handler_sca.n_dofs();
+  {
+    BlockCompressedSimpleSparsityPattern csp (dim,dim);
+    for (unsigned int d=0; d<dim; ++d)
+      for (unsigned int e=0; e<dim; ++e)
+        csp.block(d,e).reinit (dofs_per_block, dofs_per_block);
+
+    csp.collect_sizes();
+
+    DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+    sparsity_pattern.copy_from (csp);
+  }
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dim);
+  for (unsigned int i=0; i<dim; ++i)
+    solution.block(i).reinit (dofs_per_block);
+  solution.collect_sizes ();
+
+  system_rhs.reinit (solution);
+
+  vec1.resize (dim);
+  vec2.resize (dim);
+  vec1[0].reinit (dofs_per_block);
+  vec2[0].reinit (vec1[0]);
+  for (unsigned int i=1; i<dim; ++i)
+    {
+      vec1[i].reinit (vec1[0]);
+      vec2[i].reinit (vec1[0]);
+    }
+
+                                // assemble curl-curl operator
+  {
+    QGauss<dim>   quadrature_formula(fe_degree+1);
+
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values    |
+                             update_JxW_values |
+                             update_gradients);
+
+    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+    const unsigned int   n_q_points      = quadrature_formula.size();
+
+    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+    const FEValuesExtractors::Vector sc (0);
+
+    std::vector<double> phi_div (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+        local_matrix = 0;
+
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            for (unsigned int k=0; k<dofs_per_cell; ++k)
+              {
+                const Tensor<2,dim> phi_grad = fe_values[sc].gradient(k,q);
+                phi_div[k] = trace(phi_grad);
+              }
+
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                for (unsigned int j=0; j<=i; ++j)
+                  {
+                    local_matrix(i,j) += (phi_div[i] * phi_div[j] *
+                                          global_coefficient)
+                      * fe_values.JxW(q);
+                  }
+              }
+          }
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+            local_matrix(i,j) = local_matrix(j,i);
+
+        cell->get_dof_indices (local_dof_indices);
+        constraints.distribute_local_to_global (local_matrix,
+                                                local_dof_indices,
+                                                system_matrix);
+      }
+  }
+
+                                // first system_rhs with random numbers
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<system_rhs.block(i).size(); ++j)
+      {
+        const double val = -1. + 2.*(double)rand()/double(RAND_MAX);
+        system_rhs.block(i)(j) = val;
+      }
+  constraints.condense(system_rhs);
+  for (unsigned int i=0; i<dim; ++i)
+    vec1[i] = system_rhs.block(i);
+
+                                // setup matrix-free structure
+  {
+    QGauss<1> quad(fe_degree+1);
+    mf_data.reinit (dof_handler_sca, constraints, quad,
+                    typename MatrixFree<dim>::AdditionalData
+                    (MPI_COMM_WORLD,
+                     MatrixFree<dim>::AdditionalData::none));
+  }
+
+  system_matrix.vmult (solution, system_rhs);
+
+  typedef std::vector<Vector<double> > VectorType;
+  MatrixFreeTest<dim,fe_degree,VectorType> mf (mf_data);
+  mf.vmult (vec2, vec1);
+
+                                // Verification
+  double error = 0.;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<system_rhs.block(i).size(); ++j)
+      error += std::fabs (solution.block(i)(j)-vec2[i](j));
+  double relative = solution.block(0).l1_norm();
+  deallog << "  Verification fe degree " << fe_degree  <<  ": "
+          << error/relative << std::endl << std::endl;
+}
+
+
+
+int main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (3);
+
+  {
+    deallog << std::endl << "Test with doubles" << std::endl << std::endl;
+    deallog.threshold_double(1.e-12);
+    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_div/cmp/generic b/tests/matrix_free/matrix_vector_div/cmp/generic
new file mode 100644
index 0000000000..619c272047
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_div/cmp/generic
@@ -0,0 +1,12 @@
+
+DEAL::
+DEAL::Test with doubles
+DEAL::
+DEAL:2d::  Verification fe degree 1: 0
+DEAL:2d::
+DEAL:2d::  Verification fe degree 2: 0
+DEAL:2d::
+DEAL:3d::  Verification fe degree 1: 0
+DEAL:3d::
+DEAL:3d::  Verification fe degree 2: 0
+DEAL:3d::