]> https://gitweb.dealii.org/ - dealii.git/commitdiff
test examples
authorTimo Heister <timo.heister@gmail.com>
Sat, 1 Jul 2023 11:04:37 +0000 (07:04 -0400)
committerMatthias Maier <tamiko@43-1.org>
Mon, 3 Jul 2023 04:03:19 +0000 (23:03 -0500)
.gitignore
tests/examples/CMakeLists.txt
tests/examples/create_diffs.sh [new file with mode: 0755]
tests/examples/step-1.diff [new file with mode: 0644]
tests/examples/step-1.output [new file with mode: 0644]
tests/examples/step-2.diff [new file with mode: 0644]
tests/examples/step-2.output [new file with mode: 0644]
tests/examples/step-56.cc [deleted file]
tests/examples/step-56.diff [new file with mode: 0644]
tests/examples/step-56.output
tests/examples/update.sh [new file with mode: 0755]

index 4d186c07f01a1985044d55b0b3e02d32d402290b..d5f9ed90111cbdf6eb2e7119375296122f529574 100644 (file)
@@ -55,3 +55,4 @@ examples/step-*/CMakeCache.txt
 examples/step-*/cmake_install.cmake
 examples/step-*/step-?
 examples/step-*/step-??
+/tests/examples/*.cc
index db9c1adcf2d84cdcde0005dfbfccf95fe0d7b872..1ff57689661719f6ede52583b8354e6412304e87 100644 (file)
@@ -1,4 +1,11 @@
 cmake_minimum_required(VERSION 3.13.4)
 include(../scripts/setup_testsubproject.cmake)
 project(testsuite CXX)
+
+# create tests from diffs:
+execute_process(
+    COMMAND "./update.sh"
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
 deal_ii_pickup_tests()
diff --git a/tests/examples/create_diffs.sh b/tests/examples/create_diffs.sh
new file mode 100755 (executable)
index 0000000..515ad3f
--- /dev/null
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+# create .diff files from current .cc files
+
+# Find all *.diff files in the current directory
+diff_files=$(find . -maxdepth 1 -type f -name "*.diff")
+
+for file in $diff_files; do
+  # Extract the filename without the extension
+  filename=$(basename "$file" .diff)
+  src="../../examples/$filename/$filename.cc"
+
+  # Check if the corresponding .cc file exists
+  if [[ -f "$src" ]]; then
+      echo "creating $file from $src and $filename.cc"
+    diff "$src" "$filename.cc" > $file
+  else
+      echo "No matching .cc file found for $diff_file"
+      exit -1
+  fi
+done
diff --git a/tests/examples/step-1.diff b/tests/examples/step-1.diff
new file mode 100644 (file)
index 0000000..e69de29
diff --git a/tests/examples/step-1.output b/tests/examples/step-1.output
new file mode 100644 (file)
index 0000000..6aa9888
--- /dev/null
@@ -0,0 +1,2 @@
+Grid written to grid-1.svg
+Grid written to grid-2.svg
diff --git a/tests/examples/step-2.diff b/tests/examples/step-2.diff
new file mode 100644 (file)
index 0000000..e69de29
diff --git a/tests/examples/step-2.output b/tests/examples/step-2.output
new file mode 100644 (file)
index 0000000..e69de29
diff --git a/tests/examples/step-56.cc b/tests/examples/step-56.cc
deleted file mode 100644 (file)
index 8f1ae48..0000000
+++ /dev/null
@@ -1,1171 +0,0 @@
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2016 - 2022 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.
- *
- * ---------------------------------------------------------------------
-
- * Author: Ryan Grove, Clemson University
- *         Timo Heister, Clemson University
- */
-
-// @sect3{Include files}
-
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/utilities.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_q.h>
-#include <deal.II/fe/fe_system.h>
-#include <deal.II/fe/fe_values.h>
-
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/grid_tools.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/affine_constraints.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/block_sparsity_pattern.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/lac/sparse_direct.h>
-#include <deal.II/lac/sparse_ilu.h>
-
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/error_estimator.h>
-#include <deal.II/numerics/matrix_tools.h>
-#include <deal.II/numerics/vector_tools.h>
-
-// We need to include the following file to do timings:
-#include <deal.II/base/timer.h>
-
-// This includes the files necessary for us to use geometric Multigrid
-#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.h>
-#include <deal.II/multigrid/multigrid.h>
-
-#include <fstream>
-#include <iostream>
-
-#include "../tests.h"
-
-namespace Step56
-{
-  // In order to make it easy to switch between the different solvers that are
-  // being used, we declare an enum that can be passed as an argument to the
-  // constructor of the main class.
-  struct SolverType
-  {
-    enum type
-    {
-      FGMRES_ILU,
-      FGMRES_GMG,
-      UMFPACK
-    };
-  };
-
-  // @sect3{Functions for Solution and Righthand side}
-  //
-  // The class Solution is used to define the boundary conditions and to
-  // compute errors of the numerical solution. Note that we need to define the
-  // values and gradients in order to compute L2 and H1 errors. Here we
-  // decided to separate the implementations for 2d and 3d using template
-  // specialization.
-  //
-  // Note that the first dim components are the velocity components
-  // and the last is the pressure.
-  template <int dim>
-  class Solution : public Function<dim>
-  {
-  public:
-    Solution()
-      : Function<dim>(dim + 1)
-    {}
-    virtual double
-    value(const Point<dim> &p, const unsigned int component = 0) const override;
-    virtual Tensor<1, dim>
-    gradient(const Point<dim> & p,
-             const unsigned int component = 0) const override;
-  };
-
-  template <>
-  double
-  Solution<2>::value(const Point<2> &p, const unsigned int component) const
-  {
-    Assert(component <= 2 + 1, ExcIndexRange(component, 0, 2 + 1));
-
-    using numbers::PI;
-    const double x = p(0);
-    const double y = p(1);
-
-    if (component == 0)
-      return sin(PI * x);
-    if (component == 1)
-      return -PI * y * cos(PI * x);
-    if (component == 2)
-      return sin(PI * x) * cos(PI * y);
-
-    return 0;
-  }
-
-  template <>
-  double
-  Solution<3>::value(const Point<3> &p, const unsigned int component) const
-  {
-    Assert(component <= 3 + 1, ExcIndexRange(component, 0, 3 + 1));
-
-    using numbers::PI;
-    const double x = p(0);
-    const double y = p(1);
-    const double z = p(2);
-
-    if (component == 0)
-      return 2.0 * sin(PI * x);
-    if (component == 1)
-      return -PI * y * cos(PI * x);
-    if (component == 2)
-      return -PI * z * cos(PI * x);
-    if (component == 3)
-      return sin(PI * x) * cos(PI * y) * sin(PI * z);
-
-    return 0;
-  }
-
-  // Note that for the gradient we need to return a Tensor<1,dim>
-  template <>
-  Tensor<1, 2>
-  Solution<2>::gradient(const Point<2> &p, const unsigned int component) const
-  {
-    Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1));
-
-    using numbers::PI;
-    const double x = p(0);
-    const double y = p(1);
-
-    Tensor<1, 2> return_value;
-    if (component == 0)
-      {
-        return_value[0] = PI * cos(PI * x);
-        return_value[1] = 0.0;
-      }
-    else if (component == 1)
-      {
-        return_value[0] = y * PI * PI * sin(PI * x);
-        return_value[1] = -PI * cos(PI * x);
-      }
-    else if (component == 2)
-      {
-        return_value[0] = PI * cos(PI * x) * cos(PI * y);
-        return_value[1] = -PI * sin(PI * x) * sin(PI * y);
-      }
-
-    return return_value;
-  }
-
-  template <>
-  Tensor<1, 3>
-  Solution<3>::gradient(const Point<3> &p, const unsigned int component) const
-  {
-    Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1));
-
-    using numbers::PI;
-    const double x = p(0);
-    const double y = p(1);
-    const double z = p(2);
-
-    Tensor<1, 3> return_value;
-    if (component == 0)
-      {
-        return_value[0] = 2 * PI * cos(PI * x);
-        return_value[1] = 0.0;
-        return_value[2] = 0.0;
-      }
-    else if (component == 1)
-      {
-        return_value[0] = y * PI * PI * sin(PI * x);
-        return_value[1] = -PI * cos(PI * x);
-        return_value[2] = 0.0;
-      }
-    else if (component == 2)
-      {
-        return_value[0] = z * PI * PI * sin(PI * x);
-        return_value[1] = 0.0;
-        return_value[2] = -PI * cos(PI * x);
-      }
-    else if (component == 3)
-      {
-        return_value[0] = PI * cos(PI * x) * cos(PI * y) * sin(PI * z);
-        return_value[1] = -PI * sin(PI * x) * sin(PI * y) * sin(PI * z);
-        return_value[2] = PI * sin(PI * x) * cos(PI * y) * cos(PI * z);
-      }
-
-    return return_value;
-  }
-
-  // Implementation of $f$. See the introduction for more information.
-  template <int dim>
-  class RightHandSide : public Function<dim>
-  {
-  public:
-    RightHandSide()
-      : Function<dim>(dim + 1)
-    {}
-
-    virtual double
-    value(const Point<dim> &p, const unsigned int component = 0) const override;
-  };
-
-  template <>
-  double
-  RightHandSide<2>::value(const Point<2> &p, const unsigned int component) const
-  {
-    Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1));
-
-    using numbers::PI;
-    double x = p(0);
-    double y = p(1);
-    if (component == 0)
-      return PI * PI * sin(PI * x) + PI * cos(PI * x) * cos(PI * y);
-    if (component == 1)
-      return -PI * PI * PI * y * cos(PI * x) - PI * sin(PI * y) * sin(PI * x);
-    if (component == 2)
-      return 0;
-
-    return 0;
-  }
-
-  template <>
-  double
-  RightHandSide<3>::value(const Point<3> &p, const unsigned int component) const
-  {
-    Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1));
-
-    using numbers::PI;
-    double x = p(0);
-    double y = p(1);
-    double z = p(2);
-    if (component == 0)
-      return 2 * PI * PI * sin(PI * x) +
-             PI * cos(PI * x) * cos(PI * y) * sin(PI * z);
-    if (component == 1)
-      return -PI * PI * PI * y * cos(PI * x) +
-             PI * (-1) * sin(PI * y) * sin(PI * x) * sin(PI * z);
-    if (component == 2)
-      return -PI * PI * PI * z * cos(PI * x) +
-             PI * cos(PI * z) * sin(PI * x) * cos(PI * y);
-    if (component == 3)
-      return 0;
-
-    return 0;
-  }
-
-
-
-  // @sect3{ASPECT BlockSchurPreconditioner}
-
-  // In the following, we will implement a preconditioner that expands
-  // on the ideas discussed in the Results section of step-22.
-  // Specifically, we
-  // 1. use an upper block-triangular preconditioner because we want to
-  // use right preconditioning.
-  // 2. optionally allow using an inner solver for the velocity block instead
-  // of a single preconditioner application.
-  // 3. do not use InverseMatrix but explicitly call SolverCG.
-  // This approach is also used in the ASPECT code
-  // (see http://aspect.dealii.org) that solves the Stokes equations in
-  // the context of simulating convection in the earth mantle, and which
-  // has been used to solve problems on many thousands of processors.
-  //
-  // The bool flag @p do_solve_A in the constructor allows us to either
-  // apply the preconditioner for the velocity block once or use an inner
-  // iterative solver for a more accurate approximation instead.
-  //
-  // Notice how we keep track of the sum of the inner iterations
-  // (preconditioner applications).
-  template <class PreconditionerAType, class PreconditionerSType>
-  class BlockSchurPreconditioner : public Subscriptor
-  {
-  public:
-    BlockSchurPreconditioner(
-      const BlockSparseMatrix<double> &system_matrix,
-      const SparseMatrix<double> &     schur_complement_matrix,
-      const PreconditionerAType &      preconditioner_A,
-      const PreconditionerSType &      preconditioner_S,
-      const bool                       do_solve_A);
-
-    void
-    vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
-
-    mutable unsigned int n_iterations_A;
-    mutable unsigned int n_iterations_S;
-
-  private:
-    const BlockSparseMatrix<double> &system_matrix;
-    const SparseMatrix<double> &     schur_complement_matrix;
-    const PreconditionerAType &      preconditioner_A;
-    const PreconditionerSType &      preconditioner_S;
-
-    const bool do_solve_A;
-  };
-
-  template <class PreconditionerAType, class PreconditionerSType>
-  BlockSchurPreconditioner<PreconditionerAType, PreconditionerSType>::
-    BlockSchurPreconditioner(
-      const BlockSparseMatrix<double> &system_matrix,
-      const SparseMatrix<double> &     schur_complement_matrix,
-      const PreconditionerAType &      preconditioner_A,
-      const PreconditionerSType &      preconditioner_S,
-      const bool                       do_solve_A)
-    : n_iterations_A(0)
-    , n_iterations_S(0)
-    , system_matrix(system_matrix)
-    , schur_complement_matrix(schur_complement_matrix)
-    , preconditioner_A(preconditioner_A)
-    , preconditioner_S(preconditioner_S)
-    , do_solve_A(do_solve_A)
-  {}
-
-
-
-  template <class PreconditionerAType, class PreconditionerSType>
-  void
-  BlockSchurPreconditioner<PreconditionerAType, PreconditionerSType>::vmult(
-    BlockVector<double> &      dst,
-    const BlockVector<double> &src) const
-  {
-    Vector<double> utmp(src.block(0));
-
-    // First solve with the approximation for S
-    {
-      SolverControl solver_control(1000, 1e-6 * src.block(1).l2_norm());
-      SolverCG<>    cg(solver_control);
-
-      dst.block(1) = 0.0;
-      cg.solve(schur_complement_matrix,
-               dst.block(1),
-               src.block(1),
-               preconditioner_S);
-
-      n_iterations_S += solver_control.last_step();
-      dst.block(1) *= -1.0;
-    }
-
-    // Second, apply the top right block (B^T)
-    {
-      system_matrix.block(0, 1).vmult(utmp, dst.block(1));
-      utmp *= -1.0;
-      utmp += src.block(0);
-    }
-
-    // Finally, either solve with the top left block
-    // or just apply one preconditioner sweep
-    if (do_solve_A == true)
-      {
-        SolverControl solver_control(10000, utmp.l2_norm() * 1e-4);
-        SolverCG<>    cg(solver_control);
-
-        dst.block(0) = 0.0;
-        cg.solve(system_matrix.block(0, 0),
-                 dst.block(0),
-                 utmp,
-                 preconditioner_A);
-
-        n_iterations_A += solver_control.last_step();
-      }
-    else
-      {
-        preconditioner_A.vmult(dst.block(0), utmp);
-        n_iterations_A += 1;
-      }
-  }
-
-  // @sect3{The StokesProblem class}
-  //
-  // This is the main class of the problem.
-  template <int dim>
-  class StokesProblem
-  {
-  public:
-    StokesProblem(const unsigned int pressure_degree,
-                  SolverType::type   solver_type);
-    void
-    run();
-
-  private:
-    void
-    setup_dofs();
-    void
-    assemble_system();
-    void
-    assemble_multigrid();
-    void
-    solve();
-    void
-    compute_errors();
-    void
-    output_results(const unsigned int refinement_cycle) const;
-
-    const unsigned int pressure_degree;
-    SolverType::type   solver_type;
-
-    Triangulation<dim> triangulation;
-    FESystem<dim>      velocity_fe;
-    FESystem<dim>      fe;
-    DoFHandler<dim>    dof_handler;
-    DoFHandler<dim>    velocity_dof_handler;
-
-    AffineConstraints<double> constraints;
-
-    BlockSparsityPattern      sparsity_pattern;
-    BlockSparseMatrix<double> system_matrix;
-    SparseMatrix<double>      pressure_mass_matrix;
-
-    BlockVector<double> solution;
-    BlockVector<double> system_rhs;
-
-    MGLevelObject<SparsityPattern>      mg_sparsity_patterns;
-    MGLevelObject<SparseMatrix<double>> mg_matrices;
-    MGLevelObject<SparseMatrix<double>> mg_interface_matrices;
-    MGConstrainedDoFs                   mg_constrained_dofs;
-
-    TimerOutput computing_timer;
-  };
-
-
-
-  template <int dim>
-  StokesProblem<dim>::StokesProblem(const unsigned int pressure_degree,
-                                    SolverType::type   solver_type)
-    : pressure_degree(pressure_degree)
-    , solver_type(solver_type)
-    , triangulation(Triangulation<dim>::maximum_smoothing)
-    ,
-    // Finite element for the velocity only:
-    velocity_fe(FE_Q<dim>(pressure_degree + 1), dim)
-    ,
-    // Finite element for the whole system:
-    fe(velocity_fe, 1, FE_Q<dim>(pressure_degree), 1)
-    , dof_handler(triangulation)
-    , velocity_dof_handler(triangulation)
-    , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
-  {}
-
-
-
-  // @sect4{StokesProblem::setup_dofs}
-
-  // This function sets up the DoFHandler, matrices, vectors, and Multigrid
-  // structures (if needed).
-  template <int dim>
-  void
-  StokesProblem<dim>::setup_dofs()
-  {
-    TimerOutput::Scope scope(computing_timer, "Setup");
-
-    system_matrix.clear();
-    pressure_mass_matrix.clear();
-
-    // The main DoFHandler only needs active DoFs, so we are not calling
-    // distribute_mg_dofs();
-    dof_handler.distribute_dofs(fe);
-
-    // This block structure separates the dim velocity components from
-    // the pressure component (used for reordering). Note that we have
-    // 2 instead of dim+1 blocks like in step-22, because our FESystem
-    // is nested and the dim velocity components appear as one block.
-    std::vector<unsigned int> block_component(2);
-    block_component[0] = 0;
-    block_component[1] = 1;
-
-    // Velocities start at component 0:
-    const FEValuesExtractors::Vector velocities(0);
-
-    // ILU behaves better if we apply a reordering to reduce fillin. There
-    // is no advantage in doing this for the other solvers.
-    if (solver_type == SolverType::FGMRES_ILU)
-      {
-        TimerOutput::Scope ilu_specific(computing_timer, "(ILU specific)");
-        DoFRenumbering::Cuthill_McKee(dof_handler);
-      }
-
-    // This ensures that all velocities DoFs are enumerated before the
-    // pressure unknowns. This allows us to use blocks for vectors and
-    // matrices and allows us to get the same DoF numbering for
-    // dof_handler and velocity_dof_handler.
-    DoFRenumbering::block_wise(dof_handler);
-
-    if (solver_type == SolverType::FGMRES_GMG)
-      {
-        TimerOutput::Scope multigrid_specific(computing_timer,
-                                              "(Multigrid specific)");
-        TimerOutput::Scope setup_multigrid(computing_timer,
-                                           "Setup - Multigrid");
-
-        // This distributes the active dofs and multigrid dofs for the
-        // velocity space in a separate DoFHandler as described in the
-        // introduction.
-        velocity_dof_handler.distribute_dofs(velocity_fe);
-        velocity_dof_handler.distribute_mg_dofs();
-
-        // The following block of code initializes the MGConstrainedDofs
-        // (using the boundary conditions for the velocity), and the
-        // sparsity patterns and matrices for each level. The resize()
-        // function of MGLevelObject<T> will destroy all existing contained
-        // objects.
-        const std::set<types::boundary_id> zero_boundary_ids = {0};
-
-        mg_constrained_dofs.clear();
-        mg_constrained_dofs.initialize(velocity_dof_handler);
-        mg_constrained_dofs.make_zero_boundary_constraints(velocity_dof_handler,
-                                                           zero_boundary_ids);
-        const unsigned int n_levels = triangulation.n_levels();
-
-        mg_interface_matrices.resize(0, n_levels - 1);
-        mg_matrices.resize(0, n_levels - 1);
-        mg_sparsity_patterns.resize(0, n_levels - 1);
-
-        for (unsigned int level = 0; level < n_levels; ++level)
-          {
-            DynamicSparsityPattern csp(velocity_dof_handler.n_dofs(level),
-                                       velocity_dof_handler.n_dofs(level));
-            MGTools::make_sparsity_pattern(velocity_dof_handler, csp, level);
-            mg_sparsity_patterns[level].copy_from(csp);
-
-            mg_matrices[level].reinit(mg_sparsity_patterns[level]);
-            mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]);
-          }
-      }
-
-    const std::vector<types::global_dof_index> dofs_per_block =
-      DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
-    const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
-
-    {
-      constraints.clear();
-      // The following makes use of a component mask for interpolation of the
-      // boundary values for the velocity only, which is further explained in
-      // the vector valued dealii step-20 tutorial.
-      DoFTools::make_hanging_node_constraints(dof_handler, constraints);
-      VectorTools::interpolate_boundary_values(dof_handler,
-                                               0,
-                                               Solution<dim>(),
-                                               constraints,
-                                               fe.component_mask(velocities));
-
-      // As discussed in the introduction, we need to fix one degree of freedom
-      // of the pressure variable to ensure solvability of the problem. We do
-      // this here by marking the first pressure dof, which has index n_u as a
-      // constrained dof.
-      if (solver_type == SolverType::UMFPACK)
-        constraints.add_line(n_u);
-
-      constraints.close();
-    }
-
-    std::cout << "\tNumber of active cells: " << triangulation.n_active_cells()
-              << std::endl
-              << "\tNumber of degrees of freedom: " << dof_handler.n_dofs()
-              << " (" << n_u << '+' << n_p << ')' << std::endl;
-
-    {
-      BlockDynamicSparsityPattern csp(dofs_per_block, dofs_per_block);
-      DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false);
-      sparsity_pattern.copy_from(csp);
-    }
-    system_matrix.reinit(sparsity_pattern);
-
-    solution.reinit(dofs_per_block);
-    system_rhs.reinit(dofs_per_block);
-  }
-
-
-  // @sect4{StokesProblem::assemble_system}
-
-  // In this function, the system matrix is assembled. We assemble the pressure
-  // mass matrix in the (1,1) block (if needed) and move it out of this location
-  // at the end of this function.
-  template <int dim>
-  void
-  StokesProblem<dim>::assemble_system()
-  {
-    TimerOutput::Scope assemble(computing_timer, "Assemble");
-    system_matrix = 0;
-    system_rhs    = 0;
-
-    // If true, we will assemble the pressure mass matrix in the (1,1) block:
-    const bool assemble_pressure_mass_matrix =
-      (solver_type == SolverType::UMFPACK) ? false : true;
-
-    QGauss<dim> quadrature_formula(pressure_degree + 2);
-
-    FEValues<dim> fe_values(fe,
-                            quadrature_formula,
-                            update_values | update_quadrature_points |
-                              update_JxW_values | update_gradients);
-
-    const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
-    const unsigned int n_q_points = quadrature_formula.size();
-
-    FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
-    Vector<double>     local_rhs(dofs_per_cell);
-
-    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
-
-    const RightHandSide<dim>    right_hand_side;
-    std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
-
-    const FEValuesExtractors::Vector velocities(0);
-    const FEValuesExtractors::Scalar pressure(dim);
-
-    std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
-    std::vector<double>                  div_phi_u(dofs_per_cell);
-    std::vector<double>                  phi_p(dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator cell =
-                                                     dof_handler.begin_active(),
-                                                   endc = dof_handler.end();
-    for (; cell != endc; ++cell)
-      {
-        fe_values.reinit(cell);
-        local_matrix = 0;
-        local_rhs    = 0;
-
-        right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
-                                          rhs_values);
-
-        for (unsigned int q = 0; q < n_q_points; ++q)
-          {
-            for (unsigned int k = 0; k < dofs_per_cell; ++k)
-              {
-                symgrad_phi_u[k] =
-                  fe_values[velocities].symmetric_gradient(k, q);
-                div_phi_u[k] = fe_values[velocities].divergence(k, q);
-                phi_p[k]     = fe_values[pressure].value(k, q);
-              }
-
-            for (unsigned int i = 0; i < dofs_per_cell; ++i)
-              {
-                for (unsigned int j = 0; j <= i; ++j)
-                  {
-                    local_matrix(i, j) +=
-                      (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) -
-                       div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] +
-                       (assemble_pressure_mass_matrix ? phi_p[i] * phi_p[j] :
-                                                        0)) *
-                      fe_values.JxW(q);
-                  }
-
-                const unsigned int component_i =
-                  fe.system_to_component_index(i).first;
-                local_rhs(i) += fe_values.shape_value(i, q) *
-                                rhs_values[q](component_i) * fe_values.JxW(q);
-              }
-          }
-
-        for (unsigned int i = 0; i < dofs_per_cell; ++i)
-          for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
-            local_matrix(i, j) = local_matrix(j, i);
-
-        cell->get_dof_indices(local_dof_indices);
-        constraints.distribute_local_to_global(local_matrix,
-                                               local_rhs,
-                                               local_dof_indices,
-                                               system_matrix,
-                                               system_rhs);
-      }
-
-    if (solver_type != SolverType::UMFPACK)
-      {
-        pressure_mass_matrix.reinit(sparsity_pattern.block(1, 1));
-        pressure_mass_matrix.copy_from(system_matrix.block(1, 1));
-        system_matrix.block(1, 1) = 0;
-      }
-  }
-
-  // @sect4{StokesProblem::assemble_multigrid}
-
-  // Here, like in step-16, we have a function that assembles the level
-  // and interface matrices necessary for the multigrid preconditioner.
-  template <int dim>
-  void
-  StokesProblem<dim>::assemble_multigrid()
-  {
-    TimerOutput::Scope multigrid_specific(computing_timer,
-                                          "(Multigrid specific)");
-    TimerOutput::Scope assemble_multigrid(computing_timer,
-                                          "Assemble Multigrid");
-
-    mg_matrices = 0.;
-
-    QGauss<dim> quadrature_formula(pressure_degree + 2);
-
-    FEValues<dim> fe_values(velocity_fe,
-                            quadrature_formula,
-                            update_values | update_quadrature_points |
-                              update_JxW_values | update_gradients);
-
-    const unsigned int dofs_per_cell = velocity_fe.dofs_per_cell;
-    const unsigned int n_q_points    = quadrature_formula.size();
-
-    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
-
-    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
-
-    const FEValuesExtractors::Vector velocities(0);
-
-    std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
-
-    std::vector<AffineConstraints<double>> boundary_constraints(
-      triangulation.n_levels());
-    std::vector<AffineConstraints<double>> boundary_interface_constraints(
-      triangulation.n_levels());
-    for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
-      {
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_refinement_edge_indices(level));
-        boundary_constraints[level].add_lines(
-          mg_constrained_dofs.get_boundary_indices(level));
-        boundary_constraints[level].close();
-
-        IndexSet idx = mg_constrained_dofs.get_refinement_edge_indices(level) &
-                       mg_constrained_dofs.get_boundary_indices(level);
-
-        boundary_interface_constraints[level].add_lines(idx);
-        boundary_interface_constraints[level].close();
-      }
-
-    // This iterator goes over all cells (not just active)
-    typename DoFHandler<dim>::cell_iterator cell = velocity_dof_handler.begin(),
-                                            endc = velocity_dof_handler.end();
-
-    for (; cell != endc; ++cell)
-      {
-        fe_values.reinit(cell);
-        cell_matrix = 0;
-
-        for (unsigned int q = 0; q < n_q_points; ++q)
-          {
-            for (unsigned int k = 0; k < dofs_per_cell; ++k)
-              symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient(k, q);
-
-            for (unsigned int i = 0; i < dofs_per_cell; ++i)
-              for (unsigned int j = 0; j <= i; ++j)
-                {
-                  cell_matrix(i, j) +=
-                    (symgrad_phi_u[i] * symgrad_phi_u[j]) * fe_values.JxW(q);
-                }
-          }
-
-        for (unsigned int i = 0; i < dofs_per_cell; ++i)
-          for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
-            cell_matrix(i, j) = cell_matrix(j, i);
-
-        cell->get_mg_dof_indices(local_dof_indices);
-
-        boundary_constraints[cell->level()].distribute_local_to_global(
-          cell_matrix, local_dof_indices, mg_matrices[cell->level()]);
-
-        for (unsigned int i = 0; i < dofs_per_cell; ++i)
-          for (unsigned int j = 0; j < dofs_per_cell; ++j)
-            if (!mg_constrained_dofs.at_refinement_edge(cell->level(),
-                                                        local_dof_indices[i]) ||
-                mg_constrained_dofs.at_refinement_edge(cell->level(),
-                                                       local_dof_indices[j]))
-              cell_matrix(i, j) = 0;
-
-        boundary_interface_constraints[cell->level()]
-          .distribute_local_to_global(cell_matrix,
-                                      local_dof_indices,
-                                      mg_interface_matrices[cell->level()]);
-      }
-  }
-
-  // @sect4{StokesProblem::solve}
-
-  // This function sets up things differently based on if you want to use ILU
-  // or GMG as a preconditioner.  Both methods share the same solver (FGMRES)
-  // but require a different preconditioner to be initialized. Here we time not
-  // only the entire solve function, but we separately time the setup of the
-  // preconditioner as well as the solve itself.
-  template <int dim>
-  void
-  StokesProblem<dim>::solve()
-  {
-    TimerOutput::Scope solve(computing_timer, "Solve");
-    constraints.set_zero(solution);
-
-    if (solver_type == SolverType::UMFPACK)
-      {
-        computing_timer.enter_subsection("(UMFPACK specific)");
-        computing_timer.enter_subsection("Solve - Initialize");
-
-        SparseDirectUMFPACK A_direct;
-        A_direct.initialize(system_matrix);
-
-        computing_timer.leave_subsection();
-        computing_timer.leave_subsection();
-
-        {
-          TimerOutput::Scope solve_backslash(computing_timer,
-                                             "Solve - Backslash");
-          A_direct.vmult(solution, system_rhs);
-        }
-
-        constraints.distribute(solution);
-        return;
-      }
-
-    // Here we must make sure to solve for the residual with "good enough"
-    // accuracy
-    SolverControl solver_control(system_matrix.m(),
-                                 1e-10 * system_rhs.l2_norm());
-    unsigned int  n_iterations_A;
-    unsigned int  n_iterations_S;
-
-    // This is used to pass whether or not we want to solve for A inside
-    // the preconditioner.  One could change this to false to see if
-    // there is still convergence and if so does the program then run
-    // faster or slower
-    const bool use_expensive = true;
-
-    SolverFGMRES<BlockVector<double>> solver(solver_control);
-
-    if (solver_type == SolverType::FGMRES_ILU)
-      {
-        computing_timer.enter_subsection("(ILU specific)");
-        computing_timer.enter_subsection("Solve - Set-up Preconditioner");
-
-        std::cout << "   Computing preconditioner..." << std::endl
-                  << std::flush;
-
-        SparseILU<double> A_preconditioner;
-        A_preconditioner.initialize(system_matrix.block(0, 0));
-
-        SparseILU<double> S_preconditioner;
-        S_preconditioner.initialize(pressure_mass_matrix);
-
-        const BlockSchurPreconditioner<SparseILU<double>, SparseILU<double>>
-          preconditioner(system_matrix,
-                         pressure_mass_matrix,
-                         A_preconditioner,
-                         S_preconditioner,
-                         use_expensive);
-
-        computing_timer.leave_subsection();
-        computing_timer.leave_subsection();
-
-        {
-          TimerOutput::Scope solve_fmgres(computing_timer, "Solve - FGMRES");
-
-          solver.solve(system_matrix, solution, system_rhs, preconditioner);
-          n_iterations_A = preconditioner.n_iterations_A;
-          n_iterations_S = preconditioner.n_iterations_S;
-        }
-      }
-    else
-      {
-        computing_timer.enter_subsection("(Multigrid specific)");
-        computing_timer.enter_subsection("Solve - Set-up Preconditioner");
-
-        // Transfer operators between levels
-        MGTransferPrebuilt<Vector<double>> mg_transfer(mg_constrained_dofs);
-        mg_transfer.build(velocity_dof_handler);
-
-        // Setup coarse grid solver
-        FullMatrix<double> coarse_matrix;
-        coarse_matrix.copy_from(mg_matrices[0]);
-        MGCoarseGridHouseholder<> coarse_grid_solver;
-        coarse_grid_solver.initialize(coarse_matrix);
-
-        using Smoother = PreconditionSOR<SparseMatrix<double>>;
-        mg::SmootherRelaxation<Smoother, Vector<double>> mg_smoother;
-        mg_smoother.initialize(mg_matrices);
-        mg_smoother.set_steps(2);
-
-        // Multigrid, when used as a preconditioner for CG, needs to be a
-        // symmetric operator, so the smoother must be symmetric
-        mg_smoother.set_symmetric(true);
-
-        mg::Matrix<Vector<double>> mg_matrix(mg_matrices);
-        mg::Matrix<Vector<double>> mg_interface_up(mg_interface_matrices);
-        mg::Matrix<Vector<double>> mg_interface_down(mg_interface_matrices);
-
-        // Now, we are ready to set up the V-cycle operator and the multilevel
-        // preconditioner.
-        Multigrid<Vector<double>> mg(
-          mg_matrix, coarse_grid_solver, mg_transfer, mg_smoother, mg_smoother);
-        mg.set_edge_matrices(mg_interface_down, mg_interface_up);
-
-        PreconditionMG<dim, Vector<double>, MGTransferPrebuilt<Vector<double>>>
-          A_Multigrid(velocity_dof_handler, mg, mg_transfer);
-
-        SparseILU<double> S_preconditioner;
-        S_preconditioner.initialize(pressure_mass_matrix,
-                                    SparseILU<double>::AdditionalData());
-
-        const BlockSchurPreconditioner<
-          PreconditionMG<dim,
-                         Vector<double>,
-                         MGTransferPrebuilt<Vector<double>>>,
-          SparseILU<double>>
-          preconditioner(system_matrix,
-                         pressure_mass_matrix,
-                         A_Multigrid,
-                         S_preconditioner,
-                         use_expensive);
-
-        computing_timer.leave_subsection();
-        computing_timer.leave_subsection();
-
-        {
-          TimerOutput::Scope solve_fmgres(computing_timer, "Solve - FGMRES");
-          solver.solve(system_matrix, solution, system_rhs, preconditioner);
-          n_iterations_A = preconditioner.n_iterations_A;
-          n_iterations_S = preconditioner.n_iterations_S;
-        }
-      }
-
-    constraints.distribute(solution);
-
-    std::cout
-      << std::endl
-      << "\tNumber of FGMRES iterations: " << solver_control.last_step()
-      << std::endl
-      << "\tTotal number of iterations used for approximation of A inverse: "
-      << n_iterations_A << std::endl
-      << "\tTotal number of iterations used for approximation of S inverse: "
-      << n_iterations_S << std::endl
-      << std::endl;
-  }
-
-
-  // @sect4{StokesProblem::process_solution}
-
-  // This function computes the L2 and H1 errors of the solution. For this,
-  // we need to make sure the pressure has mean zero.
-  template <int dim>
-  void
-  StokesProblem<dim>::compute_errors()
-  {
-    // Compute the mean pressure $\frac{1}{\Omega} \int_{\Omega} p(x) dx $
-    // and then subtract it from each pressure coefficient. This will result
-    // in a pressure with mean value zero. Here we make use of the fact that
-    // the pressure is component $dim$ and that the finite element space
-    // is nodal.
-    const double mean_pressure = VectorTools::compute_mean_value(
-      dof_handler, QGauss<dim>(pressure_degree + 2), solution, dim);
-    solution.block(1).add(-mean_pressure);
-    std::cout << "   Note: The mean value was adjusted by " << -mean_pressure
-              << std::endl;
-
-    const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
-    const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
-                                                     dim + 1);
-
-    Vector<float> difference_per_cell(triangulation.n_active_cells());
-    VectorTools::integrate_difference(dof_handler,
-                                      solution,
-                                      Solution<dim>(),
-                                      difference_per_cell,
-                                      QGauss<dim>(pressure_degree + 2),
-                                      VectorTools::L2_norm,
-                                      &velocity_mask);
-
-    const double Velocity_L2_error =
-      VectorTools::compute_global_error(triangulation,
-                                        difference_per_cell,
-                                        VectorTools::L2_norm);
-
-    VectorTools::integrate_difference(dof_handler,
-                                      solution,
-                                      Solution<dim>(),
-                                      difference_per_cell,
-                                      QGauss<dim>(pressure_degree + 2),
-                                      VectorTools::L2_norm,
-                                      &pressure_mask);
-
-    const double Pressure_L2_error =
-      VectorTools::compute_global_error(triangulation,
-                                        difference_per_cell,
-                                        VectorTools::L2_norm);
-
-    VectorTools::integrate_difference(dof_handler,
-                                      solution,
-                                      Solution<dim>(),
-                                      difference_per_cell,
-                                      QGauss<dim>(pressure_degree + 2),
-                                      VectorTools::H1_norm,
-                                      &velocity_mask);
-
-    const double Velocity_H1_error =
-      VectorTools::compute_global_error(triangulation,
-                                        difference_per_cell,
-                                        VectorTools::H1_norm);
-
-    std::cout << std::endl
-              << "   Velocity L2 Error: " << Velocity_L2_error << std::endl
-              << "   Pressure L2 Error: " << Pressure_L2_error << std::endl
-              << "   Velocity H1 Error: " << Velocity_H1_error << std::endl;
-  }
-
-
-  // @sect4{StokesProblem::output_results}
-
-  // This function generates graphical output like it is done in step-22.
-  template <int dim>
-  void
-  StokesProblem<dim>::output_results(const unsigned int refinement_cycle) const
-  {
-    std::vector<std::string> solution_names(dim, "velocity");
-    solution_names.emplace_back("pressure");
-
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-      data_component_interpretation(
-        dim, DataComponentInterpretation::component_is_part_of_vector);
-    data_component_interpretation.push_back(
-      DataComponentInterpretation::component_is_scalar);
-
-    DataOut<dim> data_out;
-    data_out.attach_dof_handler(dof_handler);
-    data_out.add_data_vector(solution,
-                             solution_names,
-                             DataOut<dim>::type_dof_data,
-                             data_component_interpretation);
-    data_out.build_patches();
-
-    std::ofstream output(
-      "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtk");
-    data_out.write_vtk(output);
-  }
-
-
-
-  // @sect4{StokesProblem::run}
-
-  // The last step in the Stokes class is, as usual, the function that
-  // generates the initial grid and calls the other functions in the
-  // respective order.
-  template <int dim>
-  void
-  StokesProblem<dim>::run()
-  {
-    GridGenerator::hyper_cube(triangulation);
-    triangulation.refine_global(2);
-
-    if (solver_type == SolverType::FGMRES_ILU)
-      std::cout << "Now running with ILU" << std::endl;
-    else if (solver_type == SolverType::FGMRES_GMG)
-      std::cout << "Now running with Multigrid" << std::endl;
-    else
-      std::cout << "Now running with UMFPACK" << std::endl;
-
-
-    for (unsigned int refinement_cycle = 0; refinement_cycle < 2;
-         ++refinement_cycle)
-      {
-        std::cout << "Refinement cycle " << refinement_cycle << std::endl;
-
-        if (refinement_cycle > 0)
-          triangulation.refine_global(1);
-
-        std::cout << "   Set-up..." << std::endl;
-        setup_dofs();
-
-        std::cout << "   Assembling..." << std::endl;
-        assemble_system();
-
-        if (solver_type == SolverType::FGMRES_GMG)
-          {
-            std::cout << "   Assembling Multigrid..." << std::endl;
-
-            assemble_multigrid();
-          }
-
-        std::cout << "   Solving..." << std::flush;
-        solve();
-
-        compute_errors();
-
-        output_results(refinement_cycle);
-
-        computing_timer.reset();
-      }
-  }
-} // namespace Step56
-
-// @sect3{The main function}
-int
-main()
-{
-  try
-    {
-      using namespace Step56;
-
-      const int degree = 1;
-      const int dim    = 3;
-      // options for SolverType: UMFPACK FGMRES_ILU FGMRES_GMG
-      {
-        StokesProblem<dim> flow_problem(degree, SolverType::FGMRES_GMG);
-        flow_problem.run();
-      }
-      {
-        StokesProblem<dim> flow_problem(degree, SolverType::FGMRES_ILU);
-        flow_problem.run();
-      }
-    }
-  catch (const 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/examples/step-56.diff b/tests/examples/step-56.diff
new file mode 100644 (file)
index 0000000..cb1bf1d
--- /dev/null
@@ -0,0 +1,29 @@
+1057c1057
+<     triangulation.refine_global(6 - dim);
+---
+>     triangulation.refine_global(2);
+1097c1097
+<         std::cout << "   VM Peak: " << mem.VmPeak << std::endl;
+---
+>         // std::cout << "   VM Peak: " << mem.VmPeak << std::endl;
+1099c1099
+<         computing_timer.print_summary();
+---
+>         // computing_timer.print_summary();
+1113c1113
+<       const int dim    = 3;
+---
+>       const int dim    = 2;
+1115,1117c1115,1122
+<       StokesProblem<dim> flow_problem(degree, SolverType::FGMRES_GMG);
+< 
+<       flow_problem.run();
+---
+>       {
+>         StokesProblem<dim> flow_problem(degree, SolverType::FGMRES_GMG);
+>         flow_problem.run();
+>       }
+>       {
+>         StokesProblem<dim> flow_problem(degree, SolverType::FGMRES_ILU);
+>         flow_problem.run();
+>       }
index b4f4b729b1eaec2c54719e6ddb7d0e1314123ced..5edf0cac0872e6f07afeadcf98f9cc77f0c1b508 100644 (file)
@@ -1,66 +1,98 @@
 Now running with Multigrid
 Refinement cycle 0
    Set-up...
-       Number of active cells: 64
-       Number of degrees of freedom: 2312 (2187+125)
+       Number of active cells: 16
+       Number of degrees of freedom: 187 (162+25)
    Assembling...
    Assembling Multigrid...
    Solving...
-       Number of FGMRES iterations: 22
-       Total number of iterations used for approximation of A inverse: 69
-       Total number of iterations used for approximation of S inverse: 23
+       Number of FGMRES iterations: 14
+       Total number of iterations used for approximation of A inverse: 45
+       Total number of iterations used for approximation of S inverse: 15
 
-   Note: The mean value was adjusted by -3.36485e-15
+   Note: The mean value was adjusted by 8.38229e-15
 
-   Velocity L2 Error: 0.00539774
-   Pressure L2 Error: 0.019341
-   Velocity H1 Error: 0.165265
+   Velocity L2 Error: 0.00340778
+   Pressure L2 Error: 0.017624
+   Velocity H1 Error: 0.105243
 Refinement cycle 1
    Set-up...
-       Number of active cells: 512
-       Number of degrees of freedom: 15468 (14739+729)
+       Number of active cells: 64
+       Number of degrees of freedom: 659 (578+81)
+   Assembling...
+   Assembling Multigrid...
+   Solving...
+       Number of FGMRES iterations: 16
+       Total number of iterations used for approximation of A inverse: 52
+       Total number of iterations used for approximation of S inverse: 17
+
+   Note: The mean value was adjusted by -1.75663e-14
+
+   Velocity L2 Error: 0.000426317
+   Pressure L2 Error: 0.00414007
+   Velocity H1 Error: 0.0264094
+Refinement cycle 2
+   Set-up...
+       Number of active cells: 256
+       Number of degrees of freedom: 2467 (2178+289)
    Assembling...
    Assembling Multigrid...
    Solving...
-       Number of FGMRES iterations: 21
-       Total number of iterations used for approximation of A inverse: 67
-       Total number of iterations used for approximation of S inverse: 22
+       Number of FGMRES iterations: 16
+       Total number of iterations used for approximation of A inverse: 52
+       Total number of iterations used for approximation of S inverse: 17
 
-   Note: The mean value was adjusted by -1.51964e-13
+   Note: The mean value was adjusted by -2.40243e-14
 
-   Velocity L2 Error: 0.000670888
-   Pressure L2 Error: 0.0036533
-   Velocity H1 Error: 0.0414704
+   Velocity L2 Error: 5.33246e-05
+   Pressure L2 Error: 0.00102063
+   Velocity H1 Error: 0.00660864
 Now running with ILU
 Refinement cycle 0
    Set-up...
-       Number of active cells: 64
-       Number of degrees of freedom: 2312 (2187+125)
+       Number of active cells: 16
+       Number of degrees of freedom: 187 (162+25)
    Assembling...
    Solving...   Computing preconditioner...
 
-       Number of FGMRES iterations: 22
-       Total number of iterations used for approximation of A inverse: 113
-       Total number of iterations used for approximation of S inverse: 23
+       Number of FGMRES iterations: 14
+       Total number of iterations used for approximation of A inverse: 75
+       Total number of iterations used for approximation of S inverse: 15
 
-   Note: The mean value was adjusted by 1.55695e-14
+   Note: The mean value was adjusted by 3.10407e-15
 
-   Velocity L2 Error: 0.00539774
-   Pressure L2 Error: 0.019341
-   Velocity H1 Error: 0.165265
+   Velocity L2 Error: 0.00340778
+   Pressure L2 Error: 0.017624
+   Velocity H1 Error: 0.105243
 Refinement cycle 1
    Set-up...
-       Number of active cells: 512
-       Number of degrees of freedom: 15468 (14739+729)
+       Number of active cells: 64
+       Number of degrees of freedom: 659 (578+81)
+   Assembling...
+   Solving...   Computing preconditioner...
+
+       Number of FGMRES iterations: 16
+       Total number of iterations used for approximation of A inverse: 134
+       Total number of iterations used for approximation of S inverse: 17
+
+   Note: The mean value was adjusted by -3.36177e-15
+
+   Velocity L2 Error: 0.000426317
+   Pressure L2 Error: 0.00414007
+   Velocity H1 Error: 0.0264094
+Refinement cycle 2
+   Set-up...
+       Number of active cells: 256
+       Number of degrees of freedom: 2467 (2178+289)
    Assembling...
    Solving...   Computing preconditioner...
 
-       Number of FGMRES iterations: 21
-       Total number of iterations used for approximation of A inverse: 180
-       Total number of iterations used for approximation of S inverse: 22
+       Number of FGMRES iterations: 16
+       Total number of iterations used for approximation of A inverse: 229
+       Total number of iterations used for approximation of S inverse: 17
 
-   Note: The mean value was adjusted by -1.52194e-13
+   Note: The mean value was adjusted by -1.16485e-15
 
-   Velocity L2 Error: 0.000670888
-   Pressure L2 Error: 0.0036533
-   Velocity H1 Error: 0.0414704
+   Velocity L2 Error: 5.33247e-05
+   Pressure L2 Error: 0.00102063
+   Velocity H1 Error: 0.00660864
diff --git a/tests/examples/update.sh b/tests/examples/update.sh
new file mode 100755 (executable)
index 0000000..cede5e4
--- /dev/null
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+# use .diff files to create .cc files
+
+# Find all *.diff files in the current directory
+diff_files=$(find . -maxdepth 1 -type f -name "*.diff")
+
+for file in $diff_files; do
+  # Extract the filename without the extension
+  filename=$(basename "$file" .diff)
+  src="../../examples/$filename/$filename.cc"
+  
+  # Check if the corresponding .cc file exists
+  if [[ -f "$src" ]]; then
+    # Apply the diff and save it to a new file
+    patch "$src" < "$file" -o "$filename.cc"
+  else
+      echo "No matching .cc file found for $diff_file"
+      exit -1
+  fi
+done

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.