--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Step-20 on a simplex mesh. Following modifications had to be made:
+// - Create a hypercube with simplices
+// - Change the FE_RT elements to Simplex::FE_P (2nd degree) and FE_DGQ to
+// Simplex::FE_DGP (1st degree)
+// These spaces are an alternative to solve the mixed Laplacian
+// (https://link.springer.com/article/10.1007/s10092-009-0009-6)
+// - Change QGauss to Simplex::QGauss and use MappingFE instead of default
+// mapping.
+
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor_function.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace Step20
+{
+ using namespace dealii;
+
+ template <int dim>
+ class MixedLaplaceProblem
+ {
+ public:
+ MixedLaplaceProblem(const unsigned int degree);
+ void
+ run();
+
+ private:
+ void
+ make_grid_and_dofs();
+ void
+ assemble_system();
+ void
+ solve();
+ void
+ compute_errors() const;
+ void
+ output_results() const;
+
+ const unsigned int degree;
+
+ Triangulation<dim> triangulation;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ BlockSparsityPattern sparsity_pattern;
+ BlockSparseMatrix<double> system_matrix;
+
+ BlockVector<double> solution;
+ BlockVector<double> system_rhs;
+ };
+
+ namespace PrescribedSolution
+ {
+ constexpr double alpha = 0.3;
+ constexpr double beta = 1;
+
+
+ template <int dim>
+ class RightHandSide : public Function<dim>
+ {
+ public:
+ RightHandSide()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> & p,
+ const unsigned int component = 0) const override;
+ };
+
+
+
+ template <int dim>
+ class PressureBoundaryValues : public Function<dim>
+ {
+ public:
+ PressureBoundaryValues()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> & p,
+ const unsigned int component = 0) const override;
+ };
+
+
+ template <int dim>
+ class ExactSolution : public Function<dim>
+ {
+ public:
+ ExactSolution()
+ : Function<dim>(dim + 1)
+ {}
+
+ virtual void
+ vector_value(const Point<dim> &p, Vector<double> &value) const override;
+ };
+
+
+ template <int dim>
+ double
+ RightHandSide<dim>::value(const Point<dim> & /*p*/,
+ const unsigned int /*component*/) const
+ {
+ return 0;
+ }
+
+
+
+ template <int dim>
+ double
+ PressureBoundaryValues<dim>::value(const Point<dim> &p,
+ const unsigned int /*component*/) const
+ {
+ return -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
+ alpha * p[0] * p[0] * p[0] / 6);
+ }
+
+
+
+ template <int dim>
+ void
+ ExactSolution<dim>::vector_value(const Point<dim> &p,
+ Vector<double> & values) const
+ {
+ Assert(values.size() == dim + 1,
+ ExcDimensionMismatch(values.size(), dim + 1));
+
+ values(0) = alpha * p[1] * p[1] / 2 + beta - alpha * p[0] * p[0] / 2;
+ values(1) = alpha * p[0] * p[1];
+ values(2) = -(alpha * p[0] * p[1] * p[1] / 2 + beta * p[0] -
+ alpha * p[0] * p[0] * p[0] / 6);
+ }
+
+ template <int dim>
+ class KInverse : public TensorFunction<2, dim>
+ {
+ public:
+ KInverse()
+ : TensorFunction<2, dim>()
+ {}
+
+ virtual void
+ value_list(const std::vector<Point<dim>> &points,
+ std::vector<Tensor<2, dim>> & values) const override;
+ };
+
+ template <int dim>
+ void
+ KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<Tensor<2, dim>> & values) const
+ {
+ (void)points;
+ AssertDimension(points.size(), values.size());
+
+ for (auto &value : values)
+ value = unit_symmetric_tensor<dim>();
+ }
+ } // namespace PrescribedSolution
+
+ template <int dim>
+ MixedLaplaceProblem<dim>::MixedLaplaceProblem(const unsigned int degree)
+ : degree(degree)
+ , fe(FESystem<dim>(Simplex::FE_P<dim>(degree), dim),
+ 1,
+ Simplex::FE_DGP<dim>(degree - 1),
+ 1)
+ , dof_handler(triangulation)
+ {}
+
+ template <int dim>
+ void
+ MixedLaplaceProblem<dim>::make_grid_and_dofs()
+ {
+ GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+ 2,
+ -1,
+ 1);
+ triangulation.refine_global(3);
+
+ dof_handler.distribute_dofs(fe);
+
+ DoFRenumbering::component_wise(dof_handler);
+
+ const std::vector<types::global_dof_index> dofs_per_component =
+ DoFTools::count_dofs_per_fe_component(dof_handler);
+ const unsigned int n_u = dofs_per_component[0] + dofs_per_component[1],
+ n_p = dofs_per_component[dim];
+
+ deallog << "Number of active cells: " << triangulation.n_active_cells()
+ << std::endl
+ << "Total number of cells: " << triangulation.n_cells() << std::endl
+ << "Number of degrees of freedom: " << dof_handler.n_dofs() << " ("
+ << n_u << '+' << n_p << ')' << std::endl;
+
+ BlockDynamicSparsityPattern dsp(2, 2);
+ dsp.block(0, 0).reinit(n_u, n_u);
+ dsp.block(1, 0).reinit(n_p, n_u);
+ dsp.block(0, 1).reinit(n_u, n_p);
+ dsp.block(1, 1).reinit(n_p, n_p);
+ dsp.collect_sizes();
+ DoFTools::make_sparsity_pattern(dof_handler, dsp);
+
+ sparsity_pattern.copy_from(dsp);
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(2);
+ solution.block(0).reinit(n_u);
+ solution.block(1).reinit(n_p);
+ solution.collect_sizes();
+
+ system_rhs.reinit(2);
+ system_rhs.block(0).reinit(n_u);
+ system_rhs.block(1).reinit(n_p);
+ system_rhs.collect_sizes();
+ }
+
+ template <int dim>
+ void
+ MixedLaplaceProblem<dim>::assemble_system()
+ {
+ Simplex::QGauss<dim> quadrature_formula(degree + 1);
+ Simplex::QGauss<dim - 1> face_quadrature_formula(degree + 1);
+ MappingFE<dim> mapping(Simplex::FE_DGP<dim>(1));
+
+ FEValues<dim> fe_values(mapping,
+ fe,
+ quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+ FEFaceValues<dim> fe_face_values(mapping,
+ fe,
+ face_quadrature_formula,
+ update_values | update_normal_vectors |
+ update_quadrature_points |
+ update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+ const unsigned int n_q_points = quadrature_formula.size();
+ const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+ FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> local_rhs(dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+ const PrescribedSolution::RightHandSide<dim> right_hand_side;
+ const PrescribedSolution::PressureBoundaryValues<dim>
+ pressure_boundary_values;
+ const PrescribedSolution::KInverse<dim> k_inverse;
+
+ std::vector<double> rhs_values(n_q_points);
+ std::vector<double> boundary_values(n_face_q_points);
+ std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
+
+ const FEValuesExtractors::Vector velocities(0);
+ const FEValuesExtractors::Scalar pressure(dim);
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ fe_values.reinit(cell);
+ local_matrix = 0;
+ local_rhs = 0;
+
+ right_hand_side.value_list(fe_values.get_quadrature_points(),
+ rhs_values);
+ k_inverse.value_list(fe_values.get_quadrature_points(),
+ k_inverse_values);
+
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);
+ const double div_phi_i_u = fe_values[velocities].divergence(i, q);
+ const double phi_i_p = fe_values[pressure].value(i, q);
+
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ {
+ const Tensor<1, dim> phi_j_u =
+ fe_values[velocities].value(j, q);
+ const double div_phi_j_u =
+ fe_values[velocities].divergence(j, q);
+ const double phi_j_p = fe_values[pressure].value(j, q);
+
+ local_matrix(i, j) +=
+ (phi_i_u * k_inverse_values[q] * phi_j_u //
+ - phi_i_p * div_phi_j_u //
+ - div_phi_i_u * phi_j_p) //
+ * fe_values.JxW(q);
+ }
+
+ local_rhs(i) += -phi_i_p * rhs_values[q] * fe_values.JxW(q);
+ }
+
+ for (const auto &face : cell->face_iterators())
+ if (face->at_boundary())
+ {
+ fe_face_values.reinit(cell, face);
+
+ pressure_boundary_values.value_list(
+ fe_face_values.get_quadrature_points(), boundary_values);
+
+ for (unsigned int q = 0; q < n_face_q_points; ++q)
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ local_rhs(i) += -(fe_face_values[velocities].value(i, q) * //
+ fe_face_values.normal_vector(q) * //
+ boundary_values[q] * //
+ fe_face_values.JxW(q));
+ }
+
+ cell->get_dof_indices(local_dof_indices);
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ system_matrix.add(local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i, j));
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ system_rhs(local_dof_indices[i]) += local_rhs(i);
+ }
+ }
+
+ template <int dim>
+ void
+ MixedLaplaceProblem<dim>::solve()
+ {
+ const auto &M = system_matrix.block(0, 0);
+ const auto &B = system_matrix.block(0, 1);
+
+ const auto &F = system_rhs.block(0);
+ const auto &G = system_rhs.block(1);
+
+ auto &U = solution.block(0);
+ auto &P = solution.block(1);
+
+ const auto op_M = linear_operator(M);
+ const auto op_B = linear_operator(B);
+
+ ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10);
+ SolverCG<Vector<double>> solver_M(reduction_control_M);
+ PreconditionJacobi<SparseMatrix<double>> preconditioner_M;
+
+ preconditioner_M.initialize(M);
+
+ const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
+
+ const auto op_S = transpose_operator(op_B) * op_M_inv * op_B;
+ const auto op_aS =
+ transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B;
+
+ IterationNumberControl iteration_number_control_aS(30, 1.e-18);
+ SolverCG<Vector<double>> solver_aS(iteration_number_control_aS);
+
+ const auto preconditioner_S =
+ inverse_operator(op_aS, solver_aS, PreconditionIdentity());
+
+ const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G;
+
+ SolverControl solver_control_S(2000, 1.e-12);
+ SolverCG<Vector<double>> solver_S(solver_control_S);
+
+ const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S);
+
+ P = op_S_inv * schur_rhs;
+
+ deallog << solver_control_S.last_step()
+ << " CG Schur complement iterations to obtain convergence."
+ << std::endl;
+
+ U = op_M_inv * (F - op_B * P);
+ }
+
+ template <int dim>
+ void
+ MixedLaplaceProblem<dim>::compute_errors() const
+ {
+ const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
+ const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
+ dim + 1);
+
+ PrescribedSolution::ExactSolution<dim> exact_solution;
+ Vector<double> cellwise_errors(triangulation.n_active_cells());
+
+ QTrapezoid<1> q_trapez;
+ QIterated<dim> quadrature(q_trapez, degree + 2);
+ MappingFE<dim> mapping(Simplex::FE_DGP<dim>(1));
+
+ VectorTools::integrate_difference(mapping,
+ dof_handler,
+ solution,
+ exact_solution,
+ cellwise_errors,
+ quadrature,
+ VectorTools::L2_norm,
+ &pressure_mask);
+ const double p_l2_error =
+ VectorTools::compute_global_error(triangulation,
+ cellwise_errors,
+ VectorTools::L2_norm);
+
+ VectorTools::integrate_difference(mapping,
+ dof_handler,
+ solution,
+ exact_solution,
+ cellwise_errors,
+ quadrature,
+ VectorTools::L2_norm,
+ &velocity_mask);
+ const double u_l2_error =
+ VectorTools::compute_global_error(triangulation,
+ cellwise_errors,
+ VectorTools::L2_norm);
+
+ deallog << "Errors: ||e_p||_L2 = " << p_l2_error
+ << ", ||e_u||_L2 = " << u_l2_error << std::endl;
+ }
+
+ template <int dim>
+ void
+ MixedLaplaceProblem<dim>::output_results() const
+ {
+ MappingFE<dim> mapping(Simplex::FE_DGP<dim>(1));
+ std::vector<std::string> solution_names(dim, "u");
+ solution_names.emplace_back("p");
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ interpretation(dim,
+ DataComponentInterpretation::component_is_part_of_vector);
+ interpretation.push_back(DataComponentInterpretation::component_is_scalar);
+
+ DataOut<dim> data_out;
+ data_out.add_data_vector(dof_handler,
+ solution,
+ solution_names,
+ interpretation);
+
+ data_out.build_patches(mapping, degree);
+
+ std::ofstream output("solution.vtu");
+ data_out.write_vtu(output);
+ }
+
+
+ template <int dim>
+ void
+ MixedLaplaceProblem<dim>::run()
+ {
+ make_grid_and_dofs();
+ assemble_system();
+ solve();
+ compute_errors();
+ output_results();
+ }
+} // namespace Step20
+
+int
+main()
+{
+ initlog();
+
+ try
+ {
+ using namespace Step20;
+
+ const unsigned int fe_degree = 2;
+ MixedLaplaceProblem<2> mixed_laplace_problem(fe_degree);
+ mixed_laplace_problem.run();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
--- /dev/null
+
+DEAL::Number of active cells: 512
+DEAL::Total number of cells: 680
+DEAL::Number of degrees of freedom: 3714 (2178+1536)
+DEAL:cg::Starting value 0.624584
+DEAL:cg::Convergence step 25 value 4.62032e-11
+DEAL:cg::Starting value 37.4651
+DEAL:cg:cg::Starting value 37.4651
+DEAL:cg:cg::Convergence step 30 value 0.0848672
+DEAL:cg:cg::Starting value 0.826345
+DEAL:cg:cg::Convergence step 25 value 5.36874e-11
+DEAL:cg:cg::Starting value 8.99626
+DEAL:cg:cg::Convergence step 30 value 0.108473
+DEAL:cg:cg::Starting value 0.257746
+DEAL:cg:cg::Convergence step 25 value 1.58401e-11
+DEAL:cg:cg::Starting value 3.18381
+DEAL:cg:cg::Convergence step 30 value 0.0283153
+DEAL:cg:cg::Starting value 0.0877011
+DEAL:cg:cg::Convergence step 25 value 6.78488e-12
+DEAL:cg:cg::Starting value 1.08833
+DEAL:cg:cg::Convergence step 30 value 0.0265170
+DEAL:cg:cg::Starting value 0.0390461
+DEAL:cg:cg::Convergence step 25 value 3.22787e-12
+DEAL:cg:cg::Starting value 0.408285
+DEAL:cg:cg::Convergence step 30 value 0.00738618
+DEAL:cg:cg::Starting value 0.0169890
+DEAL:cg:cg::Convergence step 25 value 1.31488e-12
+DEAL:cg:cg::Starting value 0.205984
+DEAL:cg:cg::Convergence step 30 value 0.00300423
+DEAL:cg:cg::Starting value 0.00694787
+DEAL:cg:cg::Convergence step 25 value 5.53097e-13
+DEAL:cg:cg::Starting value 0.0555164
+DEAL:cg:cg::Convergence step 30 value 0.00172572
+DEAL:cg:cg::Starting value 0.00208394
+DEAL:cg:cg::Convergence step 25 value 1.92428e-13
+DEAL:cg:cg::Starting value 0.0202742
+DEAL:cg:cg::Convergence step 30 value 0.000727924
+DEAL:cg:cg::Starting value 0.000839806
+DEAL:cg:cg::Convergence step 25 value 7.63218e-14
+DEAL:cg:cg::Starting value 0.00899883
+DEAL:cg:cg::Convergence step 30 value 0.000386101
+DEAL:cg:cg::Starting value 0.000366998
+DEAL:cg:cg::Convergence step 25 value 2.97977e-14
+DEAL:cg:cg::Starting value 0.00371456
+DEAL:cg:cg::Convergence step 30 value 0.000143672
+DEAL:cg:cg::Starting value 0.000149522
+DEAL:cg:cg::Convergence step 25 value 1.25747e-14
+DEAL:cg:cg::Starting value 0.00151238
+DEAL:cg:cg::Convergence step 30 value 2.76485e-05
+DEAL:cg:cg::Starting value 5.80561e-05
+DEAL:cg:cg::Convergence step 25 value 4.95193e-15
+DEAL:cg:cg::Starting value 0.000568491
+DEAL:cg:cg::Convergence step 30 value 1.51566e-05
+DEAL:cg:cg::Starting value 2.14606e-05
+DEAL:cg:cg::Convergence step 25 value 1.85827e-15
+DEAL:cg:cg::Starting value 0.000231159
+DEAL:cg:cg::Convergence step 30 value 4.26986e-06
+DEAL:cg:cg::Starting value 8.83789e-06
+DEAL:cg:cg::Convergence step 25 value 7.59785e-16
+DEAL:cg:cg::Starting value 9.93855e-05
+DEAL:cg:cg::Convergence step 30 value 1.66407e-06
+DEAL:cg:cg::Starting value 3.98372e-06
+DEAL:cg:cg::Convergence step 25 value 3.25992e-16
+DEAL:cg:cg::Starting value 3.74185e-05
+DEAL:cg:cg::Convergence step 30 value 1.12542e-06
+DEAL:cg:cg::Starting value 1.44245e-06
+DEAL:cg:cg::Convergence step 25 value 1.13963e-16
+DEAL:cg:cg::Starting value 1.41827e-05
+DEAL:cg:cg::Convergence step 30 value 3.57707e-07
+DEAL:cg:cg::Starting value 5.33178e-07
+DEAL:cg:cg::Convergence step 25 value 4.28754e-17
+DEAL:cg:cg::Starting value 5.62814e-06
+DEAL:cg:cg::Convergence step 30 value 1.27950e-07
+DEAL:cg:cg::Starting value 2.23213e-07
+DEAL:cg:cg::Convergence step 25 value 1.81003e-17
+DEAL:cg:cg::Starting value 2.12806e-06
+DEAL:cg:cg::Convergence step 30 value 4.90715e-08
+DEAL:cg:cg::Starting value 8.59275e-08
+DEAL:cg:cg::Convergence step 25 value 6.89855e-18
+DEAL:cg:cg::Starting value 8.70957e-07
+DEAL:cg:cg::Convergence step 30 value 1.25220e-08
+DEAL:cg:cg::Starting value 3.33966e-08
+DEAL:cg:cg::Convergence step 25 value 2.74359e-18
+DEAL:cg:cg::Starting value 2.94918e-07
+DEAL:cg:cg::Convergence step 30 value 7.19353e-09
+DEAL:cg:cg::Starting value 1.20206e-08
+DEAL:cg:cg::Convergence step 25 value 9.72271e-19
+DEAL:cg:cg::Starting value 1.20514e-07
+DEAL:cg:cg::Convergence step 30 value 2.48469e-09
+DEAL:cg:cg::Starting value 4.68186e-09
+DEAL:cg:cg::Convergence step 25 value 3.97912e-19
+DEAL:cg:cg::Starting value 5.18972e-08
+DEAL:cg:cg::Convergence step 30 value 1.83964e-09
+DEAL:cg:cg::Starting value 2.00490e-09
+DEAL:cg:cg::Convergence step 24 value 4.00689e-19
+DEAL:cg:cg::Starting value 1.99515e-08
+DEAL:cg:cg::Convergence step 30 value 9.89302e-10
+DEAL:cg:cg::Starting value 7.90511e-10
+DEAL:cg:cg::Convergence step 23 value 4.17184e-19
+DEAL:cg:cg::Starting value 7.90659e-09
+DEAL:cg:cg::Convergence step 30 value 3.01626e-10
+DEAL:cg:cg::Starting value 3.15023e-10
+DEAL:cg:cg::Convergence step 22 value 4.65119e-19
+DEAL:cg:cg::Starting value 3.05591e-09
+DEAL:cg:cg::Convergence step 30 value 8.78930e-11
+DEAL:cg:cg::Starting value 1.23874e-10
+DEAL:cg:cg::Convergence step 21 value 4.01642e-19
+DEAL:cg:cg::Starting value 1.44274e-09
+DEAL:cg:cg::Convergence step 30 value 3.28342e-11
+DEAL:cg:cg::Starting value 5.36043e-11
+DEAL:cg:cg::Convergence step 20 value 4.90919e-19
+DEAL:cg:cg::Starting value 5.39812e-10
+DEAL:cg:cg::Convergence step 30 value 1.21691e-11
+DEAL:cg:cg::Starting value 2.06924e-11
+DEAL:cg:cg::Convergence step 19 value 5.18474e-19
+DEAL:cg:cg::Starting value 1.95959e-10
+DEAL:cg:cg::Convergence step 30 value 3.48296e-12
+DEAL:cg:cg::Starting value 7.18708e-12
+DEAL:cg:cg::Convergence step 18 value 4.56473e-19
+DEAL:cg:cg::Starting value 7.93580e-11
+DEAL:cg:cg::Convergence step 30 value 1.28153e-12
+DEAL:cg:cg::Starting value 2.99859e-12
+DEAL:cg:cg::Convergence step 17 value 4.66330e-19
+DEAL:cg:cg::Starting value 3.15137e-11
+DEAL:cg:cg::Convergence step 30 value 5.03608e-13
+DEAL:cg:cg::Starting value 1.21338e-12
+DEAL:cg:cg::Convergence step 16 value 4.39333e-19
+DEAL:cg:cg::Starting value 1.28779e-11
+DEAL:cg:cg::Convergence step 30 value 3.26557e-13
+DEAL:cg:cg::Starting value 4.68242e-13
+DEAL:cg:cg::Convergence step 15 value 4.13608e-19
+DEAL:cg:cg::Starting value 4.57873e-12
+DEAL:cg:cg::Convergence step 30 value 1.99147e-13
+DEAL:cg:cg::Starting value 1.70285e-13
+DEAL:cg:cg::Convergence step 14 value 4.69229e-19
+DEAL:cg:cg::Starting value 1.63301e-12
+DEAL:cg:cg::Convergence step 30 value 5.33467e-14
+DEAL:cg:cg::Starting value 6.71379e-14
+DEAL:cg:cg::Convergence step 13 value 4.72485e-19
+DEAL:cg::Convergence step 33 value 7.32840e-13
+DEAL::33 CG Schur complement iterations to obtain convergence.
+DEAL:cg::Starting value 0.143776
+DEAL:cg::Convergence step 19 value 1.12728e-11
+DEAL::Errors: ||e_p||_L2 = 0.00221522, ||e_u||_L2 = 2.44171e-09
\ No newline at end of file