]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add a base class for Matrix-free operators 3252/head
authorDenis Davydov <davydden@gmail.com>
Mon, 17 Oct 2016 08:34:23 +0000 (10:34 +0200)
committerDenis Davydov <davydden@gmail.com>
Mon, 17 Oct 2016 13:33:32 +0000 (15:33 +0200)
doc/news/changes.h
include/deal.II/matrix_free/operators.h
tests/matrix_free/parallel_multigrid_adaptive_02.cc [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_adaptive_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/matrix_free/parallel_multigrid_adaptive_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=7.output [new file with mode: 0644]

index 4cada38b9d0a1d5e71a1f36a8ecf5c8c91cbc712..f14bfe495a9a0d0c64b48ab0d2536a6f505737c3 100644 (file)
@@ -392,6 +392,11 @@ inconvenience this causes.
 
 <ol>
 
+ <li> New: Add a base class for matrix-free operators MatrixFreeOperators::Base.
+ <br>
+ (Denis Davydov, 2016/10/16)
+ </li>
+
  <li> New: There is now a function FEEvaluation::JxW() to return the Jacobian
  determinant times the quadrature weight in the matrix-free evaluation
  routines similary to FEValues.
index f1c2e02428b23650a5a1435d47de303b48a9fe97..01f1c924fb6c16040d39fcd144aa5fc23d635673 100644 (file)
@@ -29,6 +29,198 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace MatrixFreeOperators
 {
+  /**
+   * Abstract base class for matrix-free operators which can be used both at
+   * the finest mesh or at a certain level in geometric multigrid.
+   *
+   * A derived class has to implement apply_add() method as well as
+   * compute_diagonal() to fill the protected member inverse_diagonal_entries.
+   * In case of a non-symmetric operator, Tapply_add() should be additionally
+   * implemented.
+   *
+   * @author Denis Davydov, 2016
+   */
+  template <int dim, typename number = double>
+  class Base : public Subscriptor
+  {
+  public:
+    /**
+     * Number typedef.
+     */
+    typedef number number_type;
+
+    /**
+     * size_type needed for preconditioner classes.
+     */
+    typedef typename LinearAlgebra::distributed::Vector<number>::size_type size_type;
+
+    /**
+     * Default constructor.
+     */
+    Base ();
+
+    /**
+     * Virtual destructor.
+     */
+    virtual ~Base();
+
+    /**
+     * Release all memory and return to a state just like after having called
+     * the default constructor.
+     */
+    void clear();
+
+    /**
+     * Initialize operator on fine scale.
+     */
+    void initialize (const MatrixFree<dim,number> &data);
+
+    /**
+     * Initialize operator on a level @p level.
+     */
+    void initialize (const MatrixFree<dim,number> &data,
+                     const MGConstrainedDoFs &mg_constrained_dofs,
+                     const unsigned int level);
+
+    /**
+     * Return the dimension of the codomain (or range) space.
+     */
+    size_type m () const;
+
+    /**
+     * Return the dimension of the domain space.
+     */
+    size_type n () const;
+
+    /**
+     * vmult operator for interface.
+     */
+    void vmult_interface_down(LinearAlgebra::distributed::Vector<number> &dst,
+                              const LinearAlgebra::distributed::Vector<number> &src) const;
+
+    /**
+     * vmult operator for interface.
+     */
+    void vmult_interface_up(LinearAlgebra::distributed::Vector<number> &dst,
+                            const LinearAlgebra::distributed::Vector<number> &src) const;
+
+    /**
+     * Matrix-vector multiplication.
+     */
+    void vmult (LinearAlgebra::distributed::Vector<number> &dst,
+                const LinearAlgebra::distributed::Vector<number> &src) const;
+
+    /**
+     * Transpose matrix-vector multiplication.
+     */
+    void Tvmult (LinearAlgebra::distributed::Vector<number> &dst,
+                 const LinearAlgebra::distributed::Vector<number> &src) const;
+
+    /**
+     * Adding Matrix-vector multiplication.
+     */
+    void vmult_add (LinearAlgebra::distributed::Vector<number> &dst,
+                    const LinearAlgebra::distributed::Vector<number> &src) const;
+
+    /**
+     * Adding transpose matrix-vector multiplication.
+     */
+    void Tvmult_add (LinearAlgebra::distributed::Vector<number> &dst,
+                     const LinearAlgebra::distributed::Vector<number> &src) const;
+
+    /**
+     * Returns the value of the matrix entry (row,col). In matrix-free context
+     * this function is valid only for row==col when diagonal is initialized.
+     */
+    number el (const unsigned int row,
+               const unsigned int col) const;
+
+    /**
+     * Determine an estimate for the memory consumption (in bytes) of this object.
+     */
+    virtual std::size_t memory_consumption () const;
+
+    /**
+     * A wrapper for initialize_dof_vector() of MatrixFree object.
+     */
+    void initialize_dof_vector (LinearAlgebra::distributed::Vector<number> &vec) const;
+
+    /**
+     * Compute diagonal of this operator.
+     *
+     * A derived class needs to implement this function and resize and fill
+     * the protected member inverse_diagonal_entries accordingly.
+     */
+    virtual void compute_diagonal () = 0;
+
+    /**
+     * Get read access to diagonal of this operator.
+     */
+    const LinearAlgebra::distributed::Vector<number> &get_matrix_diagonal_inverse() const;
+
+    /**
+     * Apply the Jacobi preconditioner, which multiplies every element of the
+     * <tt>src</tt> vector by the inverse of the respective diagonal element and
+     * multiplies the result with the relaxation factor <tt>omega</tt>.
+     */
+    void precondition_Jacobi(LinearAlgebra::distributed::Vector<number> &dst,
+                             const LinearAlgebra::distributed::Vector<number> &src,
+                             const number omega) const;
+
+  protected:
+
+    /**
+     * Apply operator to @p src and add result in @p dst.
+     */
+    virtual void apply_add(LinearAlgebra::distributed::Vector<number> &dst,
+                           const LinearAlgebra::distributed::Vector<number> &src) const = 0;
+
+    /**
+     * Apply transpose operator to @p src and add result in @p dst.
+     *
+     * Default implementation is to call apply_add().
+     */
+    virtual void Tapply_add(LinearAlgebra::distributed::Vector<number> &dst,
+                            const LinearAlgebra::distributed::Vector<number> &src) const;
+
+    /**
+     * MatrixFree object to be used with this operator.
+     */
+    SmartPointer<const MatrixFree<dim,number>, Base<dim,number> > data;
+
+    /**
+     * A vector to store inverse of diagonal elements.
+     */
+    LinearAlgebra::distributed::Vector<number> inverse_diagonal_entries;
+
+    /**
+     * Indices of DoFs on edge in case the operator is used in GMG context.
+     */
+    std::vector<unsigned int> edge_constrained_indices;
+
+  private:
+    /**
+     * Auxiliary vector.
+     */
+    mutable std::vector<std::pair<number,number> > edge_constrained_values;
+
+    /**
+     * A flag which determines whether or not this operator has interface
+     * matrices in GMG context.
+     */
+    bool have_interface_matrices;
+
+    /**
+     * Function which implements vmult_add (@p transpose = false) and
+     * Tvmult_add (@p transpose = true).
+     */
+    void mult_add (LinearAlgebra::distributed::Vector<number> &dst,
+                   const LinearAlgebra::distributed::Vector<number> &src,
+                   const bool transpose) const;
+  };
+
+
+
   /**
    * This class implements the operation of the action of the inverse of a
    * mass matrix on an element for the special case of an evaluation object
@@ -205,6 +397,322 @@ namespace MatrixFreeOperators
       }
   }
 
+  //----------------- Base operator -----------------------------
+  template <int dim, typename number>
+  Base<dim,number>::~Base ()
+  {
+  }
+
+
+
+  template <int dim, typename number>
+  Base<dim,number>::Base ()
+    :
+    Subscriptor(),
+    data(NULL),
+    have_interface_matrices(false)
+  {
+  }
+
+
+
+  template <int dim, typename number>
+  typename Base<dim,number>::size_type
+  Base<dim,number>::m () const
+  {
+    Assert(data != NULL,
+           ExcNotInitialized());
+    return data->get_vector_partitioner()->size();
+  }
+
+
+
+  template <int dim, typename number>
+  typename Base<dim,number>::size_type
+  Base<dim,number>::n () const
+  {
+    return m();
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::clear ()
+  {
+    data = NULL;
+    inverse_diagonal_entries.reinit(0);
+  }
+
+
+
+  template <int dim, typename number>
+  number
+  Base<dim,number>::el (const unsigned int row,
+                        const unsigned int col) const
+  {
+    Assert (row == col, ExcNotImplemented());
+    Assert (inverse_diagonal_entries.size() > 0, ExcNotInitialized());
+    return 1.0/inverse_diagonal_entries(row);
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::initialize_dof_vector (LinearAlgebra::distributed::Vector<number> &vec) const
+  {
+    Assert(data != NULL,
+           ExcNotInitialized());
+    if (!vec.partitioners_are_compatible(*data->get_dof_info(0).vector_partitioner))
+      data->initialize_dof_vector(vec);
+    Assert(vec.partitioners_are_globally_compatible(*data->get_dof_info(0).vector_partitioner),
+           ExcInternalError());
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::
+  initialize (const MatrixFree<dim,number> &data_)
+  {
+    data =  SmartPointer<const MatrixFree<dim,number>, Base<dim,number> >(&data_,typeid(*this).name());
+    edge_constrained_indices.clear();
+    have_interface_matrices = false;
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::
+  initialize (const MatrixFree<dim,number> &data_,
+              const MGConstrainedDoFs &mg_constrained_dofs,
+              const unsigned int level)
+  {
+    AssertThrow (level != numbers::invalid_unsigned_int,
+                 ExcMessage("level is not set"));
+
+    data =  SmartPointer<const MatrixFree<dim,number>, Base<dim,number> >(&data_,typeid(*this).name());
+
+    // setup edge_constrained indices
+    std::vector<types::global_dof_index> interface_indices;
+    mg_constrained_dofs.get_refinement_edge_indices(level).fill_index_vector(interface_indices);
+    edge_constrained_indices.clear();
+    edge_constrained_indices.reserve(interface_indices.size());
+    edge_constrained_values.resize(interface_indices.size());
+    const IndexSet &locally_owned = data->get_dof_handler().locally_owned_mg_dofs(level);
+    for (unsigned int i=0; i<interface_indices.size(); ++i)
+      if (locally_owned.is_element(interface_indices[i]))
+        edge_constrained_indices.push_back(locally_owned.index_within_set(interface_indices[i]));
+    have_interface_matrices = Utilities::MPI::max((unsigned int)edge_constrained_indices.size(),
+                                                  data_.get_vector_partitioner()->get_communicator()) > 0;
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::vmult (LinearAlgebra::distributed::Vector<number>       &dst,
+                           const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    dst = 0;
+    vmult_add (dst, src);
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::vmult_add (LinearAlgebra::distributed::Vector<number> &dst,
+                               const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    mult_add (dst, src, false);
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::Tvmult_add (LinearAlgebra::distributed::Vector<number> &dst,
+                                const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    mult_add (dst, src, true);
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::mult_add (LinearAlgebra::distributed::Vector<number> &dst,
+                              const LinearAlgebra::distributed::Vector<number> &src,
+                              const bool transpose) const
+  {
+    Assert(src.partitioners_are_globally_compatible(*data->get_dof_info(0).vector_partitioner), ExcInternalError());
+    Assert(dst.partitioners_are_globally_compatible(*data->get_dof_info(0).vector_partitioner), ExcInternalError());
+
+    // set zero Dirichlet values on the input vector (and remember the src and
+    // dst values because we need to reset them at the end)
+    for (unsigned int i=0; i<edge_constrained_indices.size(); ++i)
+      {
+        edge_constrained_values[i] =
+          std::pair<number,number>(src.local_element(edge_constrained_indices[i]),
+                                   dst.local_element(edge_constrained_indices[i]));
+        const_cast<LinearAlgebra::distributed::Vector<number>&>(src).local_element(edge_constrained_indices[i]) = 0.;
+      }
+
+    if (transpose)
+      Tapply_add(dst,src);
+    else
+      apply_add(dst,src);
+
+    const std::vector<unsigned int> &
+    constrained_dofs = data->get_constrained_dofs();
+    for (unsigned int i=0; i<constrained_dofs.size(); ++i)
+      dst.local_element(constrained_dofs[i]) += src.local_element(constrained_dofs[i]);
+
+    // reset edge constrained values, multiply by unit matrix and add into
+    // destination
+    for (unsigned int i=0; i<edge_constrained_indices.size(); ++i)
+      {
+        const_cast<LinearAlgebra::distributed::Vector<number>&>(src).local_element(edge_constrained_indices[i]) = edge_constrained_values[i].first;
+        dst.local_element(edge_constrained_indices[i]) = edge_constrained_values[i].second + edge_constrained_values[i].first;
+      }
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::
+  vmult_interface_down(LinearAlgebra::distributed::Vector<number> &dst,
+                       const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    Assert(src.partitioners_are_globally_compatible(*data->get_dof_info(0).vector_partitioner), ExcInternalError());
+    Assert(dst.partitioners_are_globally_compatible(*data->get_dof_info(0).vector_partitioner), ExcInternalError());
+
+    dst = 0;
+
+    if (!have_interface_matrices)
+      return;
+
+    // set zero Dirichlet values on the input vector (and remember the src and
+    // dst values because we need to reset them at the end)
+    for (unsigned int i=0; i<edge_constrained_indices.size(); ++i)
+      {
+        edge_constrained_values[i] =
+          std::pair<number,number>(src.local_element(edge_constrained_indices[i]),
+                                   dst.local_element(edge_constrained_indices[i]));
+        const_cast<LinearAlgebra::distributed::Vector<number>&>(src).local_element(edge_constrained_indices[i]) = 0.;
+      }
+
+    apply_add(dst,src);
+
+    unsigned int c=0;
+    for (unsigned int i=0; i<edge_constrained_indices.size(); ++i)
+      {
+        for ( ; c<edge_constrained_indices[i]; ++c)
+          dst.local_element(c) = 0.;
+        ++c;
+
+        // reset the src values
+        const_cast<LinearAlgebra::distributed::Vector<number>&>(src).local_element(edge_constrained_indices[i]) = edge_constrained_values[i].first;
+      }
+    for ( ; c<dst.local_size(); ++c)
+      dst.local_element(c) = 0.;
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::
+  vmult_interface_up(LinearAlgebra::distributed::Vector<number> &dst,
+                     const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    Assert(src.partitioners_are_globally_compatible(*data->get_dof_info(0).vector_partitioner), ExcInternalError());
+    Assert(dst.partitioners_are_globally_compatible(*data->get_dof_info(0).vector_partitioner), ExcInternalError());
+
+    dst = 0;
+
+    if (!have_interface_matrices)
+      return;
+
+    LinearAlgebra::distributed::Vector<number> src_cpy (src);
+    unsigned int c=0;
+    for (unsigned int i=0; i<edge_constrained_indices.size(); ++i)
+      {
+        for ( ; c<edge_constrained_indices[i]; ++c)
+          src_cpy.local_element(c) = 0.;
+        ++c;
+      }
+    for ( ; c<src_cpy.local_size(); ++c)
+      src_cpy.local_element(c) = 0.;
+
+    apply_add(dst,src_cpy);
+
+    for (unsigned int i=0; i<edge_constrained_indices.size(); ++i)
+      {
+        dst.local_element(edge_constrained_indices[i]) = 0.;
+      }
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::Tvmult (LinearAlgebra::distributed::Vector<number>       &dst,
+                            const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    dst = 0;
+    Tvmult_add (dst,src);
+  }
+
+
+
+  template <int dim, typename number>
+  std::size_t
+  Base<dim,number>::memory_consumption () const
+  {
+    return inverse_diagonal_entries.memory_consumption();
+  }
+
+
+
+  template <int dim, typename number>
+  const LinearAlgebra::distributed::Vector<number> &
+  Base<dim,number>::get_matrix_diagonal_inverse() const
+  {
+    Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized());
+    return inverse_diagonal_entries;
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::Tapply_add(LinearAlgebra::distributed::Vector<number> &dst,
+                               const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    apply_add(dst,src);
+  }
+
+
+
+  template <int dim, typename number>
+  void
+  Base<dim,number>::precondition_Jacobi(LinearAlgebra::distributed::Vector<number> &dst,
+                                        const LinearAlgebra::distributed::Vector<number> &src,
+                                        const number omega) const
+  {
+    Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized());
+
+    dst = src;
+    dst.scale(inverse_diagonal_entries);
+    dst*= omega;
+  }
+
 } // end of namespace MatrixFreeOperators
 
 
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_02.cc b/tests/matrix_free/parallel_multigrid_adaptive_02.cc
new file mode 100644 (file)
index 0000000..279217e
--- /dev/null
@@ -0,0 +1,453 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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 parallel_multigird_adaptive, but derive from
+// MatrixFreeOperators::Base class.
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+std::ofstream logfile("output");
+
+using namespace dealii::MatrixFreeOperators;
+
+
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1, typename number=double>
+class LaplaceOperator : public MatrixFreeOperators::Base<dim, number>
+{
+public:
+  LaplaceOperator()
+    :
+    MatrixFreeOperators::Base<dim, number>()
+  {};
+
+  void compute_diagonal ()
+  {
+    unsigned int dummy = 0;
+    LinearAlgebra::distributed::Vector<number> &inverse_diagonal_entries = Base<dim,number>::inverse_diagonal_entries;
+    this->initialize_dof_vector(inverse_diagonal_entries);
+    Base<dim,number>::
+    data->cell_loop (&LaplaceOperator::local_diagonal_cell,
+                     this, inverse_diagonal_entries, dummy);
+
+    const std::vector<unsigned int> &
+    constrained_dofs = Base<dim,number>::data->get_constrained_dofs();
+    for (unsigned int i=0; i<constrained_dofs.size(); ++i)
+      inverse_diagonal_entries.local_element(constrained_dofs[i]) = 0.;
+    for (unsigned int i=0; i<Base<dim,number>::edge_constrained_indices.size(); ++i)
+      inverse_diagonal_entries.local_element(Base<dim,number>::edge_constrained_indices[i]) = 0.;
+
+    for (unsigned int i=0; i<inverse_diagonal_entries.local_size(); ++i)
+      if (std::abs(inverse_diagonal_entries.local_element(i)) > 1e-10)
+        inverse_diagonal_entries.local_element(i) = 1./inverse_diagonal_entries.local_element(i);
+      else
+        inverse_diagonal_entries.local_element(i) = 1.;
+  }
+
+protected:
+
+
+  void apply_add (LinearAlgebra::distributed::Vector<number>       &dst,
+                  const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    Base<dim,number>::
+    data->cell_loop (&LaplaceOperator::local_apply, this, dst, src);
+  }
+
+
+private:
+
+  void
+  local_apply (const MatrixFree<dim,number>                &data,
+               LinearAlgebra::distributed::Vector<number>       &dst,
+               const LinearAlgebra::distributed::Vector<number> &src,
+               const std::pair<unsigned int,unsigned int>  &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,n_q_points_1d,1,number> phi (data);
+
+    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_gradient (phi.get_gradient(q), q);
+        phi.integrate (false,true);
+        phi.distribute_local_to_global (dst);
+      }
+  }
+
+  void
+  local_diagonal_cell (const MatrixFree<dim,number>                &data,
+                       LinearAlgebra::distributed::Vector<number>       &dst,
+                       const unsigned int &,
+                       const std::pair<unsigned int,unsigned int>  &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,n_q_points_1d,1,number> phi (data);
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        phi.reinit (cell);
+
+        VectorizedArray<number> local_diagonal_vector[phi.tensor_dofs_per_cell];
+        for (unsigned int i=0; i<phi.dofs_per_cell; ++i)
+          {
+            for (unsigned int j=0; j<phi.dofs_per_cell; ++j)
+              phi.begin_dof_values()[j] = VectorizedArray<number>();
+            phi.begin_dof_values()[i] = 1.;
+            phi.evaluate (false,true,false);
+            for (unsigned int q=0; q<phi.n_q_points; ++q)
+              phi.submit_gradient (phi.get_gradient(q), q);
+            phi.integrate (false,true);
+            local_diagonal_vector[i] = phi.begin_dof_values()[i];
+          }
+        for (unsigned int i=0; i<phi.tensor_dofs_per_cell; ++i)
+          phi.begin_dof_values()[i] = local_diagonal_vector[i];
+        phi.distribute_local_to_global (dst);
+      }
+  }
+
+};
+
+
+
+template <typename LAPLACEOPERATOR>
+class MGInterfaceMatrix : public Subscriptor
+{
+public:
+  void initialize (const LAPLACEOPERATOR &laplace)
+  {
+    this->laplace = &laplace;
+  }
+
+  void vmult (LinearAlgebra::distributed::Vector<double> &dst,
+              const LinearAlgebra::distributed::Vector<double> &src) const
+  {
+    laplace->vmult_interface_down(dst, src);
+  }
+
+  void Tvmult (LinearAlgebra::distributed::Vector<double> &dst,
+               const LinearAlgebra::distributed::Vector<double> &src) const
+  {
+    laplace->vmult_interface_up(dst, src);
+  }
+
+private:
+  SmartPointer<const LAPLACEOPERATOR> laplace;
+};
+
+
+
+template <typename LAPLACEOPERATOR>
+class MGTransferMF : public MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> >
+{
+public:
+  MGTransferMF(const MGLevelObject<LAPLACEOPERATOR> &laplace,
+               const MGConstrainedDoFs &mg_constrained_dofs)
+    :
+    MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> >(mg_constrained_dofs),
+    laplace_operator (laplace)
+  {
+  }
+
+  /**
+   * Overload copy_to_mg from MGTransferPrebuilt to get the vectors compatible
+   * with MatrixFree and bypass the crude vector initialization in
+   * MGTransferPrebuilt
+   */
+  template <int dim, class InVector, int spacedim>
+  void
+  copy_to_mg (const DoFHandler<dim,spacedim> &mg_dof_handler,
+              MGLevelObject<LinearAlgebra::distributed::Vector<double> > &dst,
+              const InVector &src) const
+  {
+    for (unsigned int level=dst.min_level();
+         level<=dst.max_level(); ++level)
+      laplace_operator[level].initialize_dof_vector(dst[level]);
+    MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> >::
+    copy_to_mg(mg_dof_handler, dst, src);
+  }
+
+private:
+  const MGLevelObject<LAPLACEOPERATOR> &laplace_operator;
+};
+
+
+
+template<typename MatrixType, typename Number>
+class MGCoarseIterative : public MGCoarseGridBase<LinearAlgebra::distributed::Vector<Number> >
+{
+public:
+  MGCoarseIterative() {}
+
+  void initialize(const MatrixType &matrix)
+  {
+    coarse_matrix = &matrix;
+  }
+
+  virtual void operator() (const unsigned int   level,
+                           LinearAlgebra::distributed::Vector<double> &dst,
+                           const LinearAlgebra::distributed::Vector<double> &src) const
+  {
+    ReductionControl solver_control (1e4, 1e-50, 1e-10);
+    SolverCG<LinearAlgebra::distributed::Vector<double> > solver_coarse (solver_control);
+    solver_coarse.solve (*coarse_matrix, dst, src, PreconditionIdentity());
+  }
+
+  const MatrixType *coarse_matrix;
+};
+
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, typename number>
+void do_test (const DoFHandler<dim>  &dof)
+{
+  if (types_are_equal<number,float>::value == true)
+    {
+      deallog.push("float");
+      deallog.threshold_double(1e-6);
+    }
+  else
+    {
+      deallog.threshold_double(5.e-11);
+    }
+
+  deallog << "Testing " << dof.get_fe().get_name();
+  deallog << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs);
+
+  // Dirichlet BC
+  ZeroFunction<dim> zero_function;
+  typename FunctionMap<dim>::type dirichlet_boundary;
+  dirichlet_boundary[0] = &zero_function;
+
+  // fine-level constraints
+  ConstraintMatrix constraints;
+  constraints.reinit(locally_relevant_dofs);
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, dirichlet_boundary,
+                                           constraints);
+  constraints.close();
+
+  // level constraints:
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof, dirichlet_boundary);
+
+  MappingQ<dim> mapping(fe_degree+1);
+
+  LaplaceOperator<dim,fe_degree,n_q_points_1d,number> fine_matrix;
+  MatrixFree<dim,number> fine_level_data;
+
+  typename MatrixFree<dim,number>::AdditionalData fine_level_additional_data;
+  fine_level_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+  fine_level_additional_data.tasks_block_size = 3;
+  fine_level_additional_data.mpi_communicator = MPI_COMM_WORLD;
+  fine_level_data.reinit (mapping, dof, constraints, QGauss<1>(n_q_points_1d),
+                          fine_level_additional_data);
+
+  fine_matrix.initialize(fine_level_data);
+  fine_matrix.compute_diagonal();
+
+
+  LinearAlgebra::distributed::Vector<number> in, sol;
+  fine_matrix.initialize_dof_vector(in);
+  fine_matrix.initialize_dof_vector(sol);
+
+  // set constant rhs vector
+  {
+    // this is to make it consistent with parallel_multigrid_adaptive.cc
+    ConstraintMatrix hanging_node_constraints;
+    hanging_node_constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof, hanging_node_constraints);
+    hanging_node_constraints.close();
+
+    for (unsigned int i=0; i<in.local_size(); ++i)
+      if (!hanging_node_constraints.is_constrained(in.get_partitioner()->local_to_global(i)))
+        in.local_element(i) = 1.;
+  }
+
+  // set up multigrid in analogy to step-37
+  typedef LaplaceOperator<dim,fe_degree,n_q_points_1d,number> LevelMatrixType;
+
+  MGLevelObject<LevelMatrixType> mg_matrices;
+  MGLevelObject<MatrixFree<dim,number> > mg_level_data;
+  mg_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
+  mg_level_data.resize(0, dof.get_triangulation().n_global_levels()-1);
+  for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
+    {
+      typename MatrixFree<dim,number>::AdditionalData mg_additional_data;
+      mg_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+      mg_additional_data.tasks_block_size = 3;
+      mg_additional_data.mpi_communicator = MPI_COMM_WORLD;
+      mg_additional_data.level_mg_handler = level;
+
+      ConstraintMatrix level_constraints;
+      IndexSet relevant_dofs;
+      DoFTools::extract_locally_relevant_level_dofs(dof, level,
+                                                    relevant_dofs);
+      level_constraints.reinit(relevant_dofs);
+      level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
+      level_constraints.close();
+
+      mg_level_data[level].reinit (mapping, dof, level_constraints, QGauss<1>(n_q_points_1d),
+                                   mg_additional_data);
+      mg_matrices[level].initialize(mg_level_data[level],
+                                    mg_constrained_dofs,
+                                    level);
+      mg_matrices[level].compute_diagonal();
+    }
+  MGLevelObject<MGInterfaceMatrix<LevelMatrixType> > mg_interface_matrices;
+  mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
+  for (unsigned int level=0; level<dof.get_triangulation().n_global_levels(); ++level)
+    mg_interface_matrices[level].initialize(mg_matrices[level]);
+
+  MGTransferMF<LevelMatrixType> mg_transfer(mg_matrices,
+                                            mg_constrained_dofs);
+  mg_transfer.build_matrices(dof);
+
+  MGCoarseIterative<LevelMatrixType,number> mg_coarse;
+  mg_coarse.initialize(mg_matrices[0]);
+
+  typedef PreconditionChebyshev<LevelMatrixType,LinearAlgebra::distributed::Vector<number> > SMOOTHER;
+  MGSmootherPrecondition<LevelMatrixType, SMOOTHER, LinearAlgebra::distributed::Vector<number> >
+  mg_smoother;
+
+  MGLevelObject<typename SMOOTHER::AdditionalData> smoother_data;
+  smoother_data.resize(0, dof.get_triangulation().n_global_levels()-1);
+  for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
+    {
+      smoother_data[level].smoothing_range = 15.;
+      smoother_data[level].degree = 5;
+      smoother_data[level].eig_cg_n_iterations = 15;
+      smoother_data[level].matrix_diagonal_inverse =
+        mg_matrices[level].get_matrix_diagonal_inverse();
+    }
+  mg_smoother.initialize(mg_matrices, smoother_data);
+
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> >
+  mg_matrix(mg_matrices);
+  mg::Matrix<LinearAlgebra::distributed::Vector<double> >
+  mg_interface(mg_interface_matrices);
+
+  Multigrid<LinearAlgebra::distributed::Vector<double> > mg(dof,
+                                                            mg_matrix,
+                                                            mg_coarse,
+                                                            mg_transfer,
+                                                            mg_smoother,
+                                                            mg_smoother);
+  mg.set_edge_matrices(mg_interface, mg_interface);
+  PreconditionMG<dim, LinearAlgebra::distributed::Vector<double>,
+                 MGTransferMF<LevelMatrixType> >
+                 preconditioner(dof, mg, mg_transfer);
+
+  {
+    ReductionControl control(30, 1e-20, 1e-7);
+    SolverCG<LinearAlgebra::distributed::Vector<double> > solver(control);
+    solver.solve(fine_matrix, sol, in, preconditioner);
+  }
+
+  if (types_are_equal<number,float>::value == true)
+    deallog.pop();
+
+  fine_matrix.clear();
+  for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
+    mg_matrices[level].clear();
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD,
+                                                 Triangulation<dim>::limit_level_difference_at_vertices,
+                                                 parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(6-dim);
+  const unsigned int n_runs = fe_degree == 1 ? 6-dim : 5-dim;
+  for (unsigned int i=0; i<n_runs; ++i)
+    {
+      for (typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+           cell != tria.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            (cell->center().norm() < 0.5 && (cell->level() < 5 ||
+                                             cell->center().norm() > 0.45)
+             ||
+             (dim == 2 && cell->center().norm() > 1.2)))
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+      FE_Q<dim> fe (fe_degree);
+      DoFHandler<dim> dof (tria);
+      dof.distribute_dofs(fe);
+      dof.distribute_mg_dofs(fe);
+
+      do_test<dim, fe_degree, fe_degree+1, double> (dof);
+    }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      deallog.attach(logfile);
+      deallog << std::setprecision (4);
+    }
+
+  {
+    deallog.threshold_double(1.e-10);
+    deallog.push("2d");
+    test<2,1>();
+    test<2,3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/parallel_multigrid_adaptive_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..e9349c5
--- /dev/null
@@ -0,0 +1,311 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9984
+DEAL:2d:cg::Failure step 15 value 0.001385
+DEAL:2d:cg::Starting value 21.93
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 5 value 7.812e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 881
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9988
+DEAL:2d:cg::Failure step 15 value 0.002581
+DEAL:2d:cg::Starting value 0.9989
+DEAL:2d:cg::Failure step 15 value 0.0001487
+DEAL:2d:cg::Starting value 27.77
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 6 value 9.245e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 2147
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9988
+DEAL:2d:cg::Failure step 15 value 0.003288
+DEAL:2d:cg::Starting value 0.9993
+DEAL:2d:cg::Failure step 15 value 0.005495
+DEAL:2d:cg::Starting value 0.9997
+DEAL:2d:cg::Failure step 15 value 0.05100
+DEAL:2d:cg::Starting value 43.26
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 6 value 1.145e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 6517
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9989
+DEAL:2d:cg::Failure step 15 value 0.005519
+DEAL:2d:cg::Starting value 0.9994
+DEAL:2d:cg::Failure step 15 value 0.01098
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 0.09572
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 0.7527
+DEAL:2d:cg::Starting value 76.92
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 6 value 1.559e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 4259
+DEAL:2d:cg::Starting value 0.9682
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.9897
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d:cg::Starting value 0.9970
+DEAL:2d:cg::Failure step 15 value 2.456e-05
+DEAL:2d:cg::Starting value 0.9992
+DEAL:2d:cg::Failure step 15 value 0.2156
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.976
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.763
+DEAL:2d:cg::Starting value 64.26
+DEAL:2d:cg:cg::Starting value 108.9
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.03777
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.0007755
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 7.094e-06
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 2.710e-07
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg::Convergence step 5 value 1.631e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 7457
+DEAL:2d:cg::Starting value 0.9682
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.9897
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d:cg::Starting value 0.9970
+DEAL:2d:cg::Failure step 15 value 2.456e-05
+DEAL:2d:cg::Starting value 0.9992
+DEAL:2d:cg::Failure step 15 value 0.2156
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.976
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.672
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 0.9896
+DEAL:2d:cg::Starting value 83.11
+DEAL:2d:cg:cg::Starting value 165.9
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.07398
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.003378
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 2.891e-05
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 1.842e-06
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg::Convergence step 5 value 7.691e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 18535
+DEAL:2d:cg::Starting value 0.9682
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.9897
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d:cg::Starting value 0.9970
+DEAL:2d:cg::Failure step 15 value 2.456e-05
+DEAL:2d:cg::Starting value 0.9992
+DEAL:2d:cg::Failure step 15 value 0.2156
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.976
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 1.886
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 0.8395
+DEAL:2d:cg::Starting value 1.000
+DEAL:2d:cg::Failure step 15 value 2.204
+DEAL:2d:cg::Starting value 131.0
+DEAL:2d:cg:cg::Starting value 384.0
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.2113
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.002203
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 3.160e-05
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 1.004e-06
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg::Convergence step 5 value 6.045e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+DEAL:3d:cg::Starting value 0.9354
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 13 value 9.706e-11
+DEAL:3d:cg::Starting value 0.9989
+DEAL:3d:cg::Convergence step 14 value 6.535e-11
+DEAL:3d:cg::Starting value 31.21
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg::Convergence step 6 value 3.615e-08
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 3694
+DEAL:3d:cg::Starting value 0.9354
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 13 value 9.706e-11
+DEAL:3d:cg::Starting value 0.9994
+DEAL:3d:cg::Failure step 15 value 1.839e-08
+DEAL:3d:cg::Starting value 0.9998
+DEAL:3d:cg::Failure step 15 value 7.469e-05
+DEAL:3d:cg::Starting value 55.00
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg::Convergence step 7 value 3.941e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 9566
+DEAL:3d:cg::Starting value 0.9354
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 13 value 9.706e-11
+DEAL:3d:cg::Starting value 0.9995
+DEAL:3d:cg::Failure step 15 value 1.479e-07
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.0004644
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 1.830e-10
+DEAL:3d:cg::Starting value 76.18
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg::Convergence step 8 value 5.376e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 7494
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 15 value 5.276e-11
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.001572
+DEAL:3d:cg::Starting value 0.9998
+DEAL:3d:cg::Failure step 15 value 0.0008241
+DEAL:3d:cg::Starting value 82.90
+DEAL:3d:cg:cg::Starting value 7.225
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.01157
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.0009262
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.0001114
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 3.295e-06
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 2.088e-07
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 2.026e-08
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg::Convergence step 7 value 8.696e-07
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 26548
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 15 value 5.276e-11
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.001572
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.008306
+DEAL:3d:cg::Starting value 1.000
+DEAL:3d:cg::Failure step 15 value 0.1026
+DEAL:3d:cg::Starting value 152.6
+DEAL:3d:cg:cg::Starting value 17.46
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.008363
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.001190
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.0001843
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 8.789e-06
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 7.697e-07
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 6.053e-08
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg::Convergence step 7 value 3.971e-06
diff --git a/tests/matrix_free/parallel_multigrid_adaptive_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=7.output b/tests/matrix_free/parallel_multigrid_adaptive_02.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=7.output
new file mode 100644 (file)
index 0000000..e9349c5
--- /dev/null
@@ -0,0 +1,311 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 507
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9984
+DEAL:2d:cg::Failure step 15 value 0.001385
+DEAL:2d:cg::Starting value 21.93
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 5 value 7.812e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 881
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9988
+DEAL:2d:cg::Failure step 15 value 0.002581
+DEAL:2d:cg::Starting value 0.9989
+DEAL:2d:cg::Failure step 15 value 0.0001487
+DEAL:2d:cg::Starting value 27.77
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 6 value 9.245e-07
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 2147
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9988
+DEAL:2d:cg::Failure step 15 value 0.003288
+DEAL:2d:cg::Starting value 0.9993
+DEAL:2d:cg::Failure step 15 value 0.005495
+DEAL:2d:cg::Starting value 0.9997
+DEAL:2d:cg::Failure step 15 value 0.05100
+DEAL:2d:cg::Starting value 43.26
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 6 value 1.145e-06
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Number of degrees of freedom: 6517
+DEAL:2d:cg::Starting value 0.8660
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9428
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d:cg::Starting value 0.9798
+DEAL:2d:cg::Convergence step 4 value 0
+DEAL:2d:cg::Starting value 0.9938
+DEAL:2d:cg::Convergence step 10 value 0
+DEAL:2d:cg::Starting value 0.9983
+DEAL:2d:cg::Failure step 15 value 2.432e-05
+DEAL:2d:cg::Starting value 0.9989
+DEAL:2d:cg::Failure step 15 value 0.005519
+DEAL:2d:cg::Starting value 0.9994
+DEAL:2d:cg::Failure step 15 value 0.01098
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 0.09572
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 0.7527
+DEAL:2d:cg::Starting value 76.92
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg:cg::Convergence step 0 value 0
+DEAL:2d:cg::Convergence step 6 value 1.559e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 4259
+DEAL:2d:cg::Starting value 0.9682
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.9897
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d:cg::Starting value 0.9970
+DEAL:2d:cg::Failure step 15 value 2.456e-05
+DEAL:2d:cg::Starting value 0.9992
+DEAL:2d:cg::Failure step 15 value 0.2156
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.976
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.763
+DEAL:2d:cg::Starting value 64.26
+DEAL:2d:cg:cg::Starting value 108.9
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.03777
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.0007755
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 7.094e-06
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 2.710e-07
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg::Convergence step 5 value 1.631e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 7457
+DEAL:2d:cg::Starting value 0.9682
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.9897
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d:cg::Starting value 0.9970
+DEAL:2d:cg::Failure step 15 value 2.456e-05
+DEAL:2d:cg::Starting value 0.9992
+DEAL:2d:cg::Failure step 15 value 0.2156
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.976
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.672
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 0.9896
+DEAL:2d:cg::Starting value 83.11
+DEAL:2d:cg:cg::Starting value 165.9
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.07398
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.003378
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 2.891e-05
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 1.842e-06
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg::Convergence step 5 value 7.691e-06
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Number of degrees of freedom: 18535
+DEAL:2d:cg::Starting value 0.9682
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.9897
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d:cg::Starting value 0.9970
+DEAL:2d:cg::Failure step 15 value 2.456e-05
+DEAL:2d:cg::Starting value 0.9992
+DEAL:2d:cg::Failure step 15 value 0.2156
+DEAL:2d:cg::Starting value 0.9998
+DEAL:2d:cg::Failure step 15 value 1.976
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 1.886
+DEAL:2d:cg::Starting value 0.9999
+DEAL:2d:cg::Failure step 15 value 0.8395
+DEAL:2d:cg::Starting value 1.000
+DEAL:2d:cg::Failure step 15 value 2.204
+DEAL:2d:cg::Starting value 131.0
+DEAL:2d:cg:cg::Starting value 384.0
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.2113
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 0.002203
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 3.160e-05
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg:cg::Starting value 1.004e-06
+DEAL:2d:cg:cg::Convergence step 3 value 0
+DEAL:2d:cg::Convergence step 5 value 6.045e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 1103
+DEAL:3d:cg::Starting value 0.9354
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 13 value 9.706e-11
+DEAL:3d:cg::Starting value 0.9989
+DEAL:3d:cg::Convergence step 14 value 6.535e-11
+DEAL:3d:cg::Starting value 31.21
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg::Convergence step 6 value 3.615e-08
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 3694
+DEAL:3d:cg::Starting value 0.9354
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 13 value 9.706e-11
+DEAL:3d:cg::Starting value 0.9994
+DEAL:3d:cg::Failure step 15 value 1.839e-08
+DEAL:3d:cg::Starting value 0.9998
+DEAL:3d:cg::Failure step 15 value 7.469e-05
+DEAL:3d:cg::Starting value 55.00
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg::Convergence step 7 value 3.941e-06
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Number of degrees of freedom: 9566
+DEAL:3d:cg::Starting value 0.9354
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 13 value 9.706e-11
+DEAL:3d:cg::Starting value 0.9995
+DEAL:3d:cg::Failure step 15 value 1.479e-07
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.0004644
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 1.830e-10
+DEAL:3d:cg::Starting value 76.18
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg:cg::Convergence step 0 value 0
+DEAL:3d:cg::Convergence step 8 value 5.376e-06
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 7494
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 15 value 5.276e-11
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.001572
+DEAL:3d:cg::Starting value 0.9998
+DEAL:3d:cg::Failure step 15 value 0.0008241
+DEAL:3d:cg::Starting value 82.90
+DEAL:3d:cg:cg::Starting value 7.225
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.01157
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.0009262
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.0001114
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 3.295e-06
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 2.088e-07
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 2.026e-08
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg::Convergence step 7 value 8.696e-07
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Number of degrees of freedom: 26548
+DEAL:3d:cg::Starting value 0.9813
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d:cg::Starting value 0.9960
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d:cg::Starting value 0.9993
+DEAL:3d:cg::Convergence step 15 value 5.276e-11
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.001572
+DEAL:3d:cg::Starting value 0.9999
+DEAL:3d:cg::Failure step 15 value 0.008306
+DEAL:3d:cg::Starting value 1.000
+DEAL:3d:cg::Failure step 15 value 0.1026
+DEAL:3d:cg::Starting value 152.6
+DEAL:3d:cg:cg::Starting value 17.46
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.008363
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.001190
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 0.0001843
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 8.789e-06
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 7.697e-07
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg:cg::Starting value 6.053e-08
+DEAL:3d:cg:cg::Convergence step 1 value 0
+DEAL:3d:cg::Convergence step 7 value 3.971e-06

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.