From: David Wells <drwells@email.unc.edu>
Date: Sat, 1 Sep 2018 17:40:25 +0000 (-0400)
Subject: Add a test based on step-37.
X-Git-Tag: v9.1.0-rc1~758^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9ebdcd4464c8370e3d7dcb36ae6fc75997665c6b;p=dealii.git

Add a test based on step-37.

The discussion of step-37 modernizations lead to the conclusion that each
proposed inhomogeneous boundary enforcement strategy should correspond to a
distinct test. This commit includes a test based on the first such enforcement
strategy.
---

diff --git a/tests/matrix_free/step-37-inhomogeneous-1.cc b/tests/matrix_free/step-37-inhomogeneous-1.cc
new file mode 100644
index 0000000000..838c1900e8
--- /dev/null
+++ b/tests/matrix_free/step-37-inhomogeneous-1.cc
@@ -0,0 +1,657 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2009 - 2018 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.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * Authors: Katharina Kormann, Martin Kronbichler, Uppsala University,
+ * 2009-2012, updated to MPI version with parallel vectors in 2016
+ */
+
+// This test verifies that the first strategy for enforcing inhomogeneous
+// boundary conditions (i.e., using read_dof_values_plain to compute the
+// contribution of the constraints to the right hand side in assemble_rhs)
+// given in step-37 works.
+
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_matrix.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/multigrid.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step37
+{
+  using namespace dealii;
+
+  const unsigned int degree_finite_element = 2;
+
+  template <int dim>
+  class ManufacturedSolution : public Function<dim>
+  {
+  public:
+    virtual double
+    value(const Point<dim> &p,
+          const unsigned int /*component*/ = 0) const override
+    {
+      const double pi = numbers::PI;
+      return std::sin(pi * p[0]) + std::cos(pi * p[1]);
+    }
+  };
+
+
+
+  template <int dim>
+  class ManufacturedForcing : public Function<dim>
+  {
+  public:
+    virtual double
+    value(const Point<dim> &p,
+          const unsigned int /*compononent*/ = 0) const override
+    {
+      return value<double>(p);
+    }
+
+    template <typename number>
+    number
+    value(const Point<dim, number> &p) const
+    {
+      const double pi = numbers::PI;
+      const number d  = 40.0 * p.square() + 1.0;
+      return 20 * pi * pi * std::cos(pi * p[1]) / d +
+             20 * pi * pi * std::sin(pi * p[0]) / d +
+             1600 * pi * p[0] * std::cos(pi * p[0]) / (d * d) -
+             1600 * pi * p[1] * std::sin(pi * p[1]) / (d * d);
+    }
+  };
+
+
+  template <int dim>
+  class Coefficient : public Function<dim>
+  {
+  public:
+    Coefficient()
+      : Function<dim>()
+    {}
+
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+    template <typename number>
+    number
+    value(const Point<dim, number> &p, const unsigned int component = 0) const;
+  };
+
+
+
+  template <int dim>
+  template <typename number>
+  number
+  Coefficient<dim>::value(const Point<dim, number> &p,
+                          const unsigned int /*component*/) const
+  {
+    return 1. / (0.05 + 2. * p.square());
+  }
+
+
+
+  template <int dim>
+  double
+  Coefficient<dim>::value(const Point<dim> & p,
+                          const unsigned int component) const
+  {
+    return value<double>(p, component);
+  }
+
+
+  template <int dim, int fe_degree, typename number>
+  class LaplaceOperator
+    : public MatrixFreeOperators::
+        Base<dim, LinearAlgebra::distributed::Vector<number>>
+  {
+  public:
+    using value_type = number;
+
+    LaplaceOperator();
+
+    void
+    clear() override;
+
+    void
+    evaluate_coefficient(const Coefficient<dim> &coefficient_function);
+
+    virtual void
+    compute_diagonal() override;
+
+    const Table<2, VectorizedArray<double>> &
+    get_coefficient() const
+    {
+      return coefficient;
+    }
+
+  private:
+    virtual void
+    apply_add(
+      LinearAlgebra::distributed::Vector<number> &      dst,
+      const LinearAlgebra::distributed::Vector<number> &src) const override;
+
+    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;
+
+    void
+    local_compute_diagonal(
+      const MatrixFree<dim, number> &              data,
+      LinearAlgebra::distributed::Vector<number> & dst,
+      const unsigned int &                         dummy,
+      const std::pair<unsigned int, unsigned int> &cell_range) const;
+
+    Table<2, VectorizedArray<number>> coefficient;
+  };
+
+
+
+  template <int dim, int fe_degree, typename number>
+  LaplaceOperator<dim, fe_degree, number>::LaplaceOperator()
+    : MatrixFreeOperators::Base<dim,
+                                LinearAlgebra::distributed::Vector<number>>()
+  {}
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim, fe_degree, number>::clear()
+  {
+    coefficient.reinit(0, 0);
+    MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number>>::
+      clear();
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim, fe_degree, number>::evaluate_coefficient(
+    const Coefficient<dim> &coefficient_function)
+  {
+    const unsigned int n_cells = this->data->n_macro_cells();
+    FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*this->data);
+
+    coefficient.reinit(n_cells, phi.n_q_points);
+    for (unsigned int cell = 0; cell < n_cells; ++cell)
+      {
+        phi.reinit(cell);
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          coefficient(cell, q) =
+            coefficient_function.value(phi.quadrature_point(q));
+      }
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim, fe_degree, number>::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, fe_degree + 1, 1, number> phi(data);
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        AssertDimension(coefficient.size(0), data.n_macro_cells());
+        AssertDimension(coefficient.size(1), phi.n_q_points);
+
+        phi.reinit(cell);
+        phi.read_dof_values(src);
+        phi.evaluate(false, true);
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          phi.submit_gradient(coefficient(cell, q) * phi.get_gradient(q), q);
+        phi.integrate(false, true);
+        phi.distribute_local_to_global(dst);
+      }
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim, fe_degree, number>::apply_add(
+    LinearAlgebra::distributed::Vector<number> &      dst,
+    const LinearAlgebra::distributed::Vector<number> &src) const
+  {
+    this->data->cell_loop(&LaplaceOperator::local_apply, this, dst, src);
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim, fe_degree, number>::compute_diagonal()
+  {
+    this->inverse_diagonal_entries.reset(
+      new DiagonalMatrix<LinearAlgebra::distributed::Vector<number>>());
+    LinearAlgebra::distributed::Vector<number> &inverse_diagonal =
+      this->inverse_diagonal_entries->get_vector();
+    this->data->initialize_dof_vector(inverse_diagonal);
+    unsigned int dummy = 0;
+    this->data->cell_loop(&LaplaceOperator::local_compute_diagonal,
+                          this,
+                          inverse_diagonal,
+                          dummy);
+
+    this->set_constrained_entries_to_one(inverse_diagonal);
+
+    for (unsigned int i = 0; i < inverse_diagonal.local_size(); ++i)
+      {
+        Assert(inverse_diagonal.local_element(i) > 0.,
+               ExcMessage("No diagonal entry in a positive definite operator "
+                          "should be zero"));
+        inverse_diagonal.local_element(i) =
+          1. / inverse_diagonal.local_element(i);
+      }
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim, fe_degree, number>::local_compute_diagonal(
+    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, fe_degree + 1, 1, number> phi(data);
+
+    AlignedVector<VectorizedArray<number>> diagonal(phi.dofs_per_cell);
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        AssertDimension(coefficient.size(0), data.n_macro_cells());
+        AssertDimension(coefficient.size(1), phi.n_q_points);
+
+        phi.reinit(cell);
+        for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < phi.dofs_per_cell; ++j)
+              phi.submit_dof_value(VectorizedArray<number>(), j);
+            phi.submit_dof_value(make_vectorized_array<number>(1.), i);
+
+            phi.evaluate(false, true);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_gradient(coefficient(cell, q) * phi.get_gradient(q),
+                                  q);
+            phi.integrate(false, true);
+            diagonal[i] = phi.get_dof_value(i);
+          }
+        for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+          phi.submit_dof_value(diagonal[i], i);
+        phi.distribute_local_to_global(dst);
+      }
+  }
+
+
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem();
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_rhs();
+    void
+    solve();
+    void
+    output_results(const unsigned int cycle) const;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    FE_Q<dim>       fe;
+    DoFHandler<dim> dof_handler;
+
+    AffineConstraints<double> constraints;
+    using SystemMatrixType =
+      LaplaceOperator<dim, degree_finite_element, double>;
+    SystemMatrixType system_matrix;
+
+    MGConstrainedDoFs mg_constrained_dofs;
+    using LevelMatrixType = LaplaceOperator<dim, degree_finite_element, float>;
+    MGLevelObject<LevelMatrixType> mg_matrices;
+
+    LinearAlgebra::distributed::Vector<double> solution;
+    LinearAlgebra::distributed::Vector<double> system_rhs;
+
+    ConditionalOStream pcout;
+  };
+
+
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem()
+    : triangulation(MPI_COMM_WORLD,
+                    Triangulation<dim>::limit_level_difference_at_vertices,
+                    parallel::distributed::Triangulation<
+                      dim>::construct_multigrid_hierarchy)
+    , fe(degree_finite_element)
+    , dof_handler(triangulation)
+    , pcout(deallog.get_file_stream(),
+            Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+  {}
+
+
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::setup_system()
+  {
+    system_matrix.clear();
+    mg_matrices.clear_elements();
+
+    dof_handler.distribute_dofs(fe);
+    dof_handler.distribute_mg_dofs();
+
+    pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+          << std::endl;
+
+    IndexSet locally_relevant_dofs;
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    constraints.clear();
+    constraints.reinit(locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             ManufacturedSolution<dim>(),
+                                             constraints);
+    constraints.close();
+
+    {
+      typename MatrixFree<dim, double>::AdditionalData additional_data;
+      additional_data.tasks_parallel_scheme =
+        MatrixFree<dim, double>::AdditionalData::none;
+      additional_data.mapping_update_flags =
+        (update_gradients | update_JxW_values | update_quadrature_points);
+      std::shared_ptr<MatrixFree<dim, double>> system_mf_storage(
+        new MatrixFree<dim, double>());
+      system_mf_storage->reinit(dof_handler,
+                                constraints,
+                                QGauss<1>(fe.degree + 1),
+                                additional_data);
+      system_matrix.initialize(system_mf_storage);
+    }
+
+    system_matrix.evaluate_coefficient(Coefficient<dim>());
+
+    system_matrix.initialize_dof_vector(solution);
+    system_matrix.initialize_dof_vector(system_rhs);
+
+    const unsigned int nlevels = triangulation.n_global_levels();
+    mg_matrices.resize(0, nlevels - 1);
+
+    std::set<types::boundary_id> dirichlet_boundary;
+    dirichlet_boundary.insert(0);
+    mg_constrained_dofs.initialize(dof_handler);
+    mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
+                                                       dirichlet_boundary);
+
+    for (unsigned int level = 0; level < nlevels; ++level)
+      {
+        IndexSet relevant_dofs;
+        DoFTools::extract_locally_relevant_level_dofs(dof_handler,
+                                                      level,
+                                                      relevant_dofs);
+        AffineConstraints<double> level_constraints;
+        level_constraints.reinit(relevant_dofs);
+        level_constraints.add_lines(
+          mg_constrained_dofs.get_boundary_indices(level));
+        level_constraints.close();
+
+        typename MatrixFree<dim, float>::AdditionalData additional_data;
+        additional_data.tasks_parallel_scheme =
+          MatrixFree<dim, float>::AdditionalData::none;
+        additional_data.mapping_update_flags =
+          (update_gradients | update_JxW_values | update_quadrature_points);
+        additional_data.level_mg_handler = level;
+        std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(
+          new MatrixFree<dim, float>());
+        mg_mf_storage_level->reinit(dof_handler,
+                                    level_constraints,
+                                    QGauss<1>(fe.degree + 1),
+                                    additional_data);
+
+        mg_matrices[level].initialize(mg_mf_storage_level,
+                                      mg_constrained_dofs,
+                                      level);
+        mg_matrices[level].evaluate_coefficient(Coefficient<dim>());
+      }
+  }
+
+
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::assemble_rhs()
+  {
+    solution = 0.0;
+    constraints.distribute(solution);
+    solution.update_ghost_values();
+    system_rhs = 0.0;
+
+    const Table<2, VectorizedArray<double>> &coefficient =
+      system_matrix.get_coefficient();
+    ManufacturedForcing<dim> forcing;
+
+    FEEvaluation<dim, degree_finite_element> phi(
+      *system_matrix.get_matrix_free());
+    for (unsigned int cell = 0;
+         cell < system_matrix.get_matrix_free()->n_macro_cells();
+         ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values_plain(solution);
+        phi.evaluate(false, true);
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          {
+            phi.submit_gradient(-coefficient(cell, q) * phi.get_gradient(q), q);
+            phi.submit_value(forcing.value(phi.quadrature_point(q)), q);
+          }
+        phi.integrate(true, true);
+        phi.distribute_local_to_global(system_rhs);
+      }
+    system_rhs.compress(VectorOperation::add);
+  }
+
+
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::solve()
+  {
+    MGTransferMatrixFree<dim, float> mg_transfer(mg_constrained_dofs);
+    mg_transfer.build(dof_handler);
+
+    using SmootherType =
+      PreconditionChebyshev<LevelMatrixType,
+                            LinearAlgebra::distributed::Vector<float>>;
+    mg::SmootherRelaxation<SmootherType,
+                           LinearAlgebra::distributed::Vector<float>>
+                                                         mg_smoother;
+    MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
+    smoother_data.resize(0, triangulation.n_global_levels() - 1);
+    for (unsigned int level = 0; level < triangulation.n_global_levels();
+         ++level)
+      {
+        if (level > 0)
+          {
+            smoother_data[level].smoothing_range     = 15.;
+            smoother_data[level].degree              = 4;
+            smoother_data[level].eig_cg_n_iterations = 10;
+          }
+        else
+          {
+            smoother_data[0].smoothing_range = 1e-3;
+            smoother_data[0].degree          = numbers::invalid_unsigned_int;
+            smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
+          }
+        mg_matrices[level].compute_diagonal();
+        smoother_data[level].preconditioner =
+          mg_matrices[level].get_matrix_diagonal_inverse();
+      }
+    mg_smoother.initialize(mg_matrices, smoother_data);
+
+    MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<float>>
+      mg_coarse;
+    mg_coarse.initialize(mg_smoother);
+
+    mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(
+      mg_matrices);
+
+    MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMatrixType>>
+      mg_interface_matrices;
+    mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1);
+    for (unsigned int level = 0; level < triangulation.n_global_levels();
+         ++level)
+      mg_interface_matrices[level].initialize(mg_matrices[level]);
+    mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_interface(
+      mg_interface_matrices);
+
+    Multigrid<LinearAlgebra::distributed::Vector<float>> mg(
+      mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
+    mg.set_edge_matrices(mg_interface, mg_interface);
+
+    PreconditionMG<dim,
+                   LinearAlgebra::distributed::Vector<float>,
+                   MGTransferMatrixFree<dim, float>>
+      preconditioner(dof_handler, mg, mg_transfer);
+
+    SolverControl solver_control(100, 1e-12 * system_rhs.l2_norm());
+    SolverCG<LinearAlgebra::distributed::Vector<double>> cg(solver_control);
+
+    constraints.set_zero(solution);
+    cg.solve(system_matrix, solution, system_rhs, preconditioner);
+
+    constraints.distribute(solution);
+  }
+
+
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::output_results(const unsigned int cycle) const
+  {
+    if (triangulation.n_global_active_cells() > 1000000)
+      return;
+
+    Vector<double> errors;
+    errors.reinit(triangulation.n_active_cells());
+    solution.update_ghost_values();
+    VectorTools::integrate_difference(dof_handler,
+                                      solution,
+                                      ManufacturedSolution<dim>(),
+                                      errors,
+                                      QIterated<dim>(QTrapez<1>(), 4),
+                                      VectorTools::NormType::Linfty_norm);
+    double max_cell_error = 1.0;
+    if (errors.begin() != errors.end())
+      max_cell_error = *std::max_element(errors.begin(), errors.end());
+    max_cell_error = Utilities::MPI::max(max_cell_error, MPI_COMM_WORLD);
+    Assert(max_cell_error != 0.0, ExcInternalError());
+    pcout << "max error: " << max_cell_error << '\n';
+    static double error = max_cell_error;
+    pcout << "error ratio: "
+          << Utilities::MPI::max(error, MPI_COMM_WORLD) / max_cell_error
+          << '\n';
+    error = max_cell_error;
+  }
+
+
+
+  template <int dim>
+  void
+  LaplaceProblem<dim>::run()
+  {
+    for (unsigned int cycle = 0; cycle < 8 - dim; ++cycle)
+      {
+        pcout << "Cycle " << cycle << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube(triangulation, 0., 1.);
+            triangulation.refine_global(3 - dim);
+          }
+        triangulation.refine_global(1);
+        setup_system();
+        assemble_rhs();
+        solve();
+        output_results(cycle);
+        pcout << std::endl;
+      }
+  }
+} // namespace Step37
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+  MPILogInitAll                    mpi_log;
+
+  {
+    Step37::LaplaceProblem<2> laplace_problem;
+    laplace_problem.run();
+  }
+
+  {
+    Step37::LaplaceProblem<3> laplace_problem;
+    laplace_problem.run();
+  }
+}
diff --git a/tests/matrix_free/step-37-inhomogeneous-1.mpirun=4.with_p4est=true.output b/tests/matrix_free/step-37-inhomogeneous-1.mpirun=4.with_p4est=true.output
new file mode 100644
index 0000000000..6121bdd29e
--- /dev/null
+++ b/tests/matrix_free/step-37-inhomogeneous-1.mpirun=4.with_p4est=true.output
@@ -0,0 +1,150 @@
+
+Cycle 0
+Number of degrees of freedom: 81
+DEAL:0:cg::Starting value 45.0807
+DEAL:0:cg::Convergence step 6 value 2.23481e-13
+max error: 0.00772526
+error ratio: 1.00000
+
+Cycle 1
+Number of degrees of freedom: 289
+DEAL:0:cg::Starting value 68.5863
+DEAL:0:cg::Convergence step 6 value 4.06929e-13
+max error: 0.000971239
+error ratio: 7.95402
+
+Cycle 2
+Number of degrees of freedom: 1089
+DEAL:0:cg::Starting value 95.4381
+DEAL:0:cg::Convergence step 6 value 5.78305e-13
+max error: 0.000120278
+error ratio: 8.07497
+
+Cycle 3
+Number of degrees of freedom: 4225
+DEAL:0:cg::Starting value 131.732
+DEAL:0:cg::Convergence step 6 value 1.31602e-12
+max error: 1.49228e-05
+error ratio: 8.05998
+
+Cycle 4
+Number of degrees of freedom: 16641
+DEAL:0:cg::Starting value 183.473
+DEAL:0:cg::Convergence step 6 value 2.60867e-12
+max error: 1.85713e-06
+error ratio: 8.03544
+
+Cycle 5
+Number of degrees of freedom: 66049
+DEAL:0:cg::Starting value 257.379
+DEAL:0:cg::Convergence step 6 value 4.17156e-12
+max error: 2.31590e-07
+error ratio: 8.01904
+
+Cycle 0
+Number of degrees of freedom: 125
+DEAL:0:cg::Starting value 6.44073
+DEAL:0:cg::Convergence step 6 value 2.97002e-14
+max error: 0.0499110
+error ratio: 1.00000
+
+Cycle 1
+Number of degrees of freedom: 729
+DEAL:0:cg::Starting value 8.68642
+DEAL:0:cg::Convergence step 6 value 6.80139e-14
+max error: 0.00756448
+error ratio: 6.59807
+
+Cycle 2
+Number of degrees of freedom: 4913
+DEAL:0:cg::Starting value 9.47362
+DEAL:0:cg::Convergence step 6 value 6.56927e-14
+max error: 0.000965192
+error ratio: 7.83728
+
+Cycle 3
+Number of degrees of freedom: 35937
+DEAL:0:cg::Starting value 9.30402
+DEAL:0:cg::Convergence step 6 value 7.06917e-14
+max error: 0.000120048
+error ratio: 8.04008
+
+Cycle 4
+Number of degrees of freedom: 274625
+DEAL:0:cg::Starting value 8.99023
+DEAL:0:cg::Convergence step 6 value 8.15983e-14
+max error: 1.49140e-05
+error ratio: 8.04933
+
+
+DEAL:1:cg::Starting value 45.0807
+DEAL:1:cg::Convergence step 6 value 2.23481e-13
+DEAL:1:cg::Starting value 68.5863
+DEAL:1:cg::Convergence step 6 value 4.06929e-13
+DEAL:1:cg::Starting value 95.4381
+DEAL:1:cg::Convergence step 6 value 5.78305e-13
+DEAL:1:cg::Starting value 131.732
+DEAL:1:cg::Convergence step 6 value 1.31602e-12
+DEAL:1:cg::Starting value 183.473
+DEAL:1:cg::Convergence step 6 value 2.60867e-12
+DEAL:1:cg::Starting value 257.379
+DEAL:1:cg::Convergence step 6 value 4.17156e-12
+DEAL:1:cg::Starting value 6.44073
+DEAL:1:cg::Convergence step 6 value 2.97002e-14
+DEAL:1:cg::Starting value 8.68642
+DEAL:1:cg::Convergence step 6 value 6.80139e-14
+DEAL:1:cg::Starting value 9.47362
+DEAL:1:cg::Convergence step 6 value 6.56927e-14
+DEAL:1:cg::Starting value 9.30402
+DEAL:1:cg::Convergence step 6 value 7.06917e-14
+DEAL:1:cg::Starting value 8.99023
+DEAL:1:cg::Convergence step 6 value 8.15983e-14
+
+
+DEAL:2:cg::Starting value 45.0807
+DEAL:2:cg::Convergence step 6 value 2.23481e-13
+DEAL:2:cg::Starting value 68.5863
+DEAL:2:cg::Convergence step 6 value 4.06929e-13
+DEAL:2:cg::Starting value 95.4381
+DEAL:2:cg::Convergence step 6 value 5.78305e-13
+DEAL:2:cg::Starting value 131.732
+DEAL:2:cg::Convergence step 6 value 1.31602e-12
+DEAL:2:cg::Starting value 183.473
+DEAL:2:cg::Convergence step 6 value 2.60867e-12
+DEAL:2:cg::Starting value 257.379
+DEAL:2:cg::Convergence step 6 value 4.17156e-12
+DEAL:2:cg::Starting value 6.44073
+DEAL:2:cg::Convergence step 6 value 2.97002e-14
+DEAL:2:cg::Starting value 8.68642
+DEAL:2:cg::Convergence step 6 value 6.80139e-14
+DEAL:2:cg::Starting value 9.47362
+DEAL:2:cg::Convergence step 6 value 6.56927e-14
+DEAL:2:cg::Starting value 9.30402
+DEAL:2:cg::Convergence step 6 value 7.06917e-14
+DEAL:2:cg::Starting value 8.99023
+DEAL:2:cg::Convergence step 6 value 8.15983e-14
+
+
+DEAL:3:cg::Starting value 45.0807
+DEAL:3:cg::Convergence step 6 value 2.23481e-13
+DEAL:3:cg::Starting value 68.5863
+DEAL:3:cg::Convergence step 6 value 4.06929e-13
+DEAL:3:cg::Starting value 95.4381
+DEAL:3:cg::Convergence step 6 value 5.78305e-13
+DEAL:3:cg::Starting value 131.732
+DEAL:3:cg::Convergence step 6 value 1.31602e-12
+DEAL:3:cg::Starting value 183.473
+DEAL:3:cg::Convergence step 6 value 2.60867e-12
+DEAL:3:cg::Starting value 257.379
+DEAL:3:cg::Convergence step 6 value 4.17156e-12
+DEAL:3:cg::Starting value 6.44073
+DEAL:3:cg::Convergence step 6 value 2.97002e-14
+DEAL:3:cg::Starting value 8.68642
+DEAL:3:cg::Convergence step 6 value 6.80139e-14
+DEAL:3:cg::Starting value 9.47362
+DEAL:3:cg::Convergence step 6 value 6.56927e-14
+DEAL:3:cg::Starting value 9.30402
+DEAL:3:cg::Convergence step 6 value 7.06917e-14
+DEAL:3:cg::Starting value 8.99023
+DEAL:3:cg::Convergence step 6 value 8.15983e-14
+