]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert step-20 to a simplex test 11545/head
authorLaura Prieto Saavedra <lauraprietosaavedra@gmail.com>
Thu, 14 Jan 2021 23:33:10 +0000 (18:33 -0500)
committerLaura Prieto Saavedra <lauraprietosaavedra@gmail.com>
Thu, 14 Jan 2021 23:33:10 +0000 (18:33 -0500)
tests/simplex/step-20.cc [new file with mode: 0644]
tests/simplex/step-20.with_simplex_support=on.output [new file with mode: 0644]

diff --git a/tests/simplex/step-20.cc b/tests/simplex/step-20.cc
new file mode 100644 (file)
index 0000000..34e482f
--- /dev/null
@@ -0,0 +1,549 @@
+// ---------------------------------------------------------------------
+//
+// 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;
+}
diff --git a/tests/simplex/step-20.with_simplex_support=on.output b/tests/simplex/step-20.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..157e9c1
--- /dev/null
@@ -0,0 +1,144 @@
+
+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

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.