]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Some documentation.
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 28 Apr 2020 21:36:19 +0000 (23:36 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 14 May 2020 22:28:43 +0000 (00:28 +0200)
examples/step-70/step-70.cc

index b1507ebef464fb21183ace78dd3bd2be3f63ab3c..1fdf21aa8f66681c31bc42a6120d7a5e815307d6 100644 (file)
  * Authors: Luca Heltai, Bruno Blais, 2019
  */
 
+// @sect3{Include files}
+// Most of these have been introduced elsewhere, we'll comment only on the new
+// ones.
+
 #include <deal.II/base/function.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/timer.h>
@@ -26,8 +30,6 @@
 #include <deal.II/lac/linear_operator.h>
 #include <deal.II/lac/linear_operator_tools.h>
 
-#include <deal.II/particles/data_out.h>
-
 #define FORCE_USE_OF_TRILINOS
 
 namespace LA
@@ -88,6 +90,39 @@ namespace LA
 #include <deal.II/numerics/error_estimator.h>
 #include <deal.II/numerics/vector_tools.h>
 
+// These are the only new include files w.r.t. step-60. In this tutorial,
+// the non-matching coupling between the solid and the fluid is computed using
+// an intermediate data structure that keeps track of how the quadrature points
+// of the solid evolve w.r.t. the fluid mesh. This data structure needs to keep
+// track of the position of the quadrature points on each cell describing the
+// solid domain, of the quadrature weights, and possibly of the normal vector
+// to each point, if the solid domain is of co-dimension one.
+//
+// Deal.II offers these facilities on the Particles namespace, through the
+// ParticleHandler class. ParticleHandler is a class that allows you to manage
+// a collection of particles (objects of type Particles::Particle), representing
+// a collection of points with some attached properties floating on a
+// parallel::distributed::Triangulation. The methods and classes on the
+// namespace Particles allows one to easily implement Particle In Cell methods
+// and particle tracing on distributed triangulations.
+//
+// We "abuse" this data structure to store information about the location of
+// solid quadrature points w.r.t. to the surrounding fluid grid, including
+// integration weights, and possibly surface normals. The reason why we use this
+// additional data structure is related to the fact that the solid and the fluid
+// grids are non-overlapping, and distributed independently among processes.
+//
+// In order to couple the two problems, we rely on the ParticleHandler class,
+// storing in each particle the position of a solid quadrature point (which is
+// in general not aligned to any of the fluid quadrature points), its weight,
+// and any other information that may be required to couple the two problems.
+//
+// Ownership of the solid quadrature points is inherited by the MPI partitioning
+// on the solid mesh itslef. The Particles so generated are later distributed to
+// the fluid mesh using the methods of the ParticleHandler class. This allows
+// transparent exchange of information between MPI processes about the
+// overlapping pattern between fluid cells and solid quadrature points.
+#include <deal.II/particles/data_out.h>
 #include <deal.II/particles/generators.h>
 #include <deal.II/particles/particle_handler.h>
 
@@ -100,6 +135,7 @@ namespace Step70
 {
   using namespace dealii;
 
+  // REMOVE THIS FUNCTION ONCE #9891 is merged.
   template <int dim,
             int spacedim,
             typename InputVectorType,
@@ -170,141 +206,163 @@ namespace Step70
     interpolated_field.compress(VectorOperation::add);
   }
 
+  // Similiarly to what we have done in step-60, we set up a class that holds
+  // all the parameters of our problem and derive it from the ParameterAcceptor
+  // class to simplify the management and creation of parameter files.
+  //
+  // The ParameterAcceptor paradigm requires all parameters to be writeable by
+  // the ParameterAcceptor methods. In order to avoid bugs that would be very
+  // difficult to trace down (such as witing things like `time = 0` instead of
+  // `time == 0`), we declare all the parameters in an external class, which is
+  // initialized before the actual StokesImmersedProblem class, and pass it to
+  // the main class as a const reference.
   template <int dim, int spacedim = dim>
   class StokesImmersedProblemParameters : public ParameterAcceptor
   {
   public:
-    StokesImmersedProblemParameters()
-      : ParameterAcceptor("Stokes Immersed Problem/")
-      , rhs("Right hand side", spacedim + 1)
-      , angular_velocity("Angular velocity", spacedim == 3 ? spacedim : 1)
-    {
-      add_parameter("Velocity degree",
-                    velocity_degree,
-                    "",
-                    this->prm,
-                    Patterns::Integer(1));
-
-      add_parameter("Number of time steps", number_of_time_steps);
-      add_parameter("Output frequency", mod_output);
-
-      add_parameter("Final time", final_time);
-
-      add_parameter("Viscosity", viscosity);
-
-      add_parameter("Nitsche penalty term", penalty_term);
-
-      add_parameter("Initial fluid refinement",
-                    initial_fluid_refinement,
-                    "Initial mesh refinement used for the fluid domain Omega");
-
-      add_parameter("Initial solid refinement",
-                    initial_solid_refinement,
-                    "Initial mesh refinement used for the solid domain Gamma");
-
-      add_parameter(
-        "Particle insertion refinement",
-        particle_insertion_refinement,
-        "Refinement of the volumetric mesh used to insert the particles");
-
-      add_parameter(
-        "Homogeneous Dirichlet boundary ids",
-        homogeneous_dirichlet_ids,
-        "Boundary Ids over which homogeneous Dirichlet boundary conditions are applied");
-
-      enter_my_subsection(this->prm);
-      this->prm.enter_subsection("Grid generation");
-      this->prm.add_parameter("Grid one generator", name_of_grid1);
-      this->prm.add_parameter("Grid one generator arguments",
-                              arguments_for_grid1);
-
-      this->prm.add_parameter("Grid two generator", name_of_grid2);
-      this->prm.add_parameter("Grid two generator arguments",
-                              arguments_for_grid2);
-
-      this->prm.add_parameter("Particle grid generator", name_of_particle_grid);
-      this->prm.add_parameter("Particle grid generator arguments",
-                              arguments_for_particle_grid);
-      this->prm.leave_subsection();
-
-      leave_my_subsection(this->prm);
-
-
-
-      enter_my_subsection(this->prm);
-      this->prm.enter_subsection("Refinement and remeshing");
-      this->prm.add_parameter("Refinement step frequency", mod_refinement);
-      this->prm.add_parameter("Refinement maximal level", max_level_refinement);
-      this->prm.add_parameter("Refinement strategy",
-                              refinement_strategy,
-                              "",
-                              Patterns::Selection(
-                                "fixed_fraction|fixed_number"));
-      this->prm.add_parameter("Refinement coarsening fraction",
-                              coarsening_fraction);
-      this->prm.add_parameter("Refinement fraction", refinement_fraction);
-      this->prm.add_parameter("Maximum number of cells", max_cells);
-
-      this->prm.leave_subsection();
-      leave_my_subsection(this->prm);
-
-      // correct the default dimension for the functions
-      rhs.declare_parameters_call_back.connect([&]() {
-        Functions::ParsedFunction<spacedim>::declare_parameters(this->prm,
-                                                                spacedim + 1);
-      });
-      angular_velocity.declare_parameters_call_back.connect([&]() {
-        Functions::ParsedFunction<spacedim>::declare_parameters(
-          this->prm, spacedim == 3 ? spacedim : 1);
-      });
-    }
-
+    // The constructor is responsible for the connection between the members of
+    // this class and the corresponding entries in the ParameterHandler. Thanks
+    // to the use of the ParameterHandler::add_parameter() method, this
+    // connection is trivial, but requires all members of this class to be
+    // writeable
+    StokesImmersedProblemParameters();
+
+    // however, since this class will be passed as a const reference to the
+    // StokesImmersedProblem class, we have to make sure we can still set the
+    // time correctly in the objects derived by the Function class defined
+    // here. In order to do so, we declare both the
+    // StokesImmersedProblemParameters::rhs and
+    // StokesImmersedProblemParameters::angular_velocity members to be mutable,
+    // and define this little helper method that sets their time to the correct
+    // value.
     void set_time(const double &time) const
     {
       rhs.set_time(time);
       angular_velocity.set_time(time);
     }
 
-    unsigned int                  velocity_degree               = 2;
-    unsigned int                  number_of_time_steps          = 1;
-    double                        viscosity                     = 1.0;
-    double                        final_time                    = 1.0;
-    unsigned int                  initial_fluid_refinement      = 3;
-    unsigned int                  initial_solid_refinement      = 3;
-    unsigned int                  particle_insertion_refinement = 1;
-    double                        penalty_term                  = 1e3;
-    std::list<types::boundary_id> homogeneous_dirichlet_ids{0, 1, 2, 3};
-    std::string                   name_of_grid1       = "hyper_cube";
-    std::string                   arguments_for_grid1 = "-1: 1: false";
-    std::string                   name_of_grid2       = "hyper_rectangle";
-    std::string                   arguments_for_grid2 =
+    // We will use a Taylor-Hood function space of arbitrary order. This
+    // parameter is used to initialize the FiniteElement space with the corret
+    // FESystem object
+    unsigned int velocity_degree = 2;
+
+    // Instead of defining a time step increment, in this tutorial we prefer to
+    // let the user choose a final simulation time, and the number of steps in
+    // which we want to reach the final time
+    unsigned int number_of_time_steps = 1;
+    double       final_time           = 1.0;
+
+    // Instead of producing an output at every time step, we allow the user to
+    // select the frequency at which output is produced:
+    unsigned int output_frequency = 1;
+
+    // We allow every grid to be refined independently. In this tutorial, no
+    // physics is resolved on the solid grid, and its velocity is given as a
+    // datum. However it relatively straight forward to incorporate some
+    // elasticity model in this tutorial, and transform it in a fully fledged
+    // FSI solver.
+    unsigned int initial_fluid_refinement      = 3;
+    unsigned int initial_solid_refinement      = 3;
+    unsigned int particle_insertion_refinement = 1;
+
+    // The only two parameters used in the equations are the viscosity of the
+    // fluid, and the penalty term used in the Nitsche formulation:
+    double viscosity    = 1.0;
+    double penalty_term = 1e3;
+
+    // By default, we create a hyper_cube without colorisation, and we use
+    // homogenous Dirichlet boundary conditions. In this set we store the
+    // boundary ids to use when setting the boundary conditions:
+    std::list<types::boundary_id> homogeneous_dirichlet_ids{0};
+
+    // We illustrate here another way to create a Triangulation from a parameter
+    // file, using the method GridGenerator::generate_from_name_and_arguments(),
+    // that takes the name of a function in the GridGenerator namespace, and its
+    // arguments as a single string representing the arguments as a tuple.
+    //
+    // The mechanism with which the arguments are parsed from and to a string is
+    // explained in detail in the Patterns::Tools::Convert class, which is
+    // used to translate from strings to most of the basic STL types (vectors,
+    // maps, tuples) and basic dealii types (Point, Tensor, BoundingBox, etc.).
+    //
+    // In general objects that can be represented by rank 1 uniform elements
+    // (i.e., std::vector<double>, Point<dim>, std::set<int>, etc.) are comma
+    // separated. Additional ranks take a semicolon, allowing you to parse
+    // strings into objects of type `std::vector<std::vector<double>>`, or,
+    // for example, `std::vector<Point<dim>>`, as `0.0, 0.1; 0.1, 0.2`. This
+    // string could be interpreted as a vector of two Point objects, or a vector
+    // of vector of doubles.
+    //
+    // When the entries are not uniform, as in the tuple case, we use a colon
+    // to separate the various entries. For example, a string like `5: 0.1, 0.2`
+    // could be used to parse an object of type `std::pair<int, Point<2>>` or a
+    // `std::tuple<int, std::vector<double>>`.
+    //
+    // In our case most of the arguments are Point objects (representing
+    // centers, corners, subdivision elements, etc.), integer values (number of
+    // subdivisions), double values (radius, lengths, etc.), or boolean options
+    // (such as the `colorize` option that many GridGenerator functions take).
+    //
+    // In the example below, we set reasonable default values, but these can be
+    // changed at run time by selecting any other supported function of the
+    // GridGenerator namespace.
+    //
+    // We do this for each of the generated grids, to be as generic as possible:
+    std::string name_of_grid1       = "hyper_cube";
+    std::string arguments_for_grid1 = "-1: 1: false";
+    std::string name_of_grid2       = "hyper_rectangle";
+    std::string arguments_for_grid2 =
       dim == 2 ? "-.5, -.1: .5, .1: false" : "-.5, -.1, -.1: .5, .1, .1: false";
     std::string name_of_particle_grid = "hyper_ball";
     std::string arguments_for_particle_grid =
       dim == 2 ? "0.3, 0.3: 0.1: false" : "0.3, 0.3, 0.3 : 0.1: false";
 
-    // Refinement parameters
+    // Similarly, we allow for different local refinement strategies. In
+    // particular, we limit the maximum number of refinement levels, in order
+    // to control the minimum size of the fluid grid, and guarantee that it is
+    // compatible with the solid grid, and we perform local refinement based
+    // on standard error estimators on the fluid velocity field.
+    //
+    // We permit the user to choose between the
+    // two most common refinement strategies, namely `fixed_number` or
+    // `fixed_fraction`, that refer to the methods
+    // GridRefinement::refine_and_coarsen_fixed_fraction() and
+    // GridRefinement::refine_and_coarsen_fixed_number().
+    //
+    // Refinement may be done every few time steps, instead of continuosly, and
+    // we control this value by the `refinement_frequency` parameter:
     int          max_level_refinement = 5;
     std::string  refinement_strategy  = "fixed_fraction";
     double       coarsening_fraction  = 0.3;
     double       refinement_fraction  = 0.3;
     unsigned int max_cells            = 1000;
-    int          mod_refinement       = 5;
-    int          mod_output           = 1;
-
+    int          refinement_frequency = 5;
+
+    // These two functions are used to control the source term of Stokes flow
+    // and the angular velocity at which we move solid. In a more realistic
+    // simulation, the solid velocity or its deformation would come from the
+    // solution of an auxiliary problem on the solid domain. In this example
+    // step we leave this part aside, and simply impose a fixed rotational
+    // velocity field on the immersed solid, governed by function that can be
+    // specified in the parameter file:
     mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>> rhs;
     mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
       angular_velocity;
   }; // namespace Step70
 
 
+  // Once the angular velocity is provided as a Function object, we reconstruct
+  // the pointwise solid velocity thrugh the following class.
   template <int spacedim>
   class SolidVelocity : public Function<spacedim>
   {
   public:
     SolidVelocity(const Functions::ParsedFunction<spacedim> &angular_velocity)
       : angular_velocity(angular_velocity)
-    {}
+    {
+      static_assert(spacedim > 1,
+                    "Cannot instatiate SolidVelocity for spacedim == 1");
+    }
 
     virtual double value(const Point<spacedim> &p,
                          unsigned int           component = 0) const
@@ -318,8 +376,7 @@ namespace Step70
 
           velocity = cross_product_3d(p, omega);
         }
-
-      if (spacedim == 2)
+      else if (spacedim == 2)
         {
           double omega = angular_velocity.value(p, 0);
 
@@ -334,7 +391,10 @@ namespace Step70
     const Functions::ParsedFunction<spacedim> &angular_velocity;
   };
 
-
+  // Similarly, we assume that the incremental solid displacement can be
+  // computed simply by a one step time integration process (here using a
+  // trivial forward Euler method), so that at each time step, the solid simply
+  // displaces by `v*dt`.
   template <int spacedim>
   class SolidDisplacement : public Function<spacedim>
   {
@@ -345,7 +405,10 @@ namespace Step70
       : Function<spacedim>(spacedim)
       , angular_velocity(angular_velocity)
       , time_step(time_step)
-    {}
+    {
+      static_assert(spacedim > 1,
+                    "Cannot instatiate SolidDisplacement for spacedim == 1");
+    }
 
     virtual double value(const Point<spacedim> &p,
                          unsigned int           component = 0) const
@@ -370,6 +433,7 @@ namespace Step70
     double                                     time_step;
   };
 
+  // We are now ready to introduce the main class of our tutorial program.
   template <int dim, int spacedim = dim>
   class StokesImmersedProblem
   {
@@ -377,18 +441,44 @@ namespace Step70
     StokesImmersedProblem(
       const StokesImmersedProblemParameters<dim, spacedim> &par);
 
+    // As usual, we leave a single public entry point to the user: the run
+    // method. Everything else is left private, and accessed through the run
+    // method itself.
     void run();
 
   private:
     void make_grid();
+
+    // These two methods are new w.r.t. previous examples, and initiliaze the
+    // ParticleHandler objects used in this class. We have two such objects: one
+    // is a passive tracer, used to plot the trajectories of fluid particles,
+    // while the the other is composed of the actual solid quadrature points,
+    // and represent material particles of the solid.
     void setup_tracer_particles();
     void setup_solid_particles();
+
+    // The setup is split in two parts: create all objects that are needed once
+    // per simulation,
     void initial_setup();
+    // followed by all objects that need to be reinitialized at every refinement
+    // step.
     void setup_dofs();
+
+    // The assembly rutine is identical to other Stokes assembly rutines,
     void assemble_stokes_system();
+    // with the exception of the Nistche restriction part, which exploits one of
+    // the particle handlers to integrate on a non-matching part of the fluid
+    // domain, corresponding to the position of the solid.
     void assemble_nitche_restriction();
+
     void solve();
-    void refine_grid();
+
+    // The refine_and_transfer() method is called only every
+    // `refinement_frequency` steps, and makes sure that all the fields
+    // that were computed on the time step before refinement are transfered
+    // correctly to the new grid. This includes vector fields, as well as
+    // particle information.
+    void refine_and_transfer();
     void output_results(const unsigned int cycle, const double time) const;
 
     void
@@ -397,29 +487,52 @@ namespace Step70
                      const unsigned int                               iter,
                      const double time) const;
 
+    // As noted before, make sure we cannot modify this object from within this
+    // class, by making it a const reference.
     const StokesImmersedProblemParameters<dim, spacedim> &par;
 
     MPI_Comm mpi_communicator;
 
-    std::unique_ptr<FiniteElement<spacedim>>      fe1;
-    std::unique_ptr<FiniteElement<dim, spacedim>> fe2;
-
-    parallel::distributed::Triangulation<spacedim>      tria1;
-    parallel::distributed::Triangulation<dim, spacedim> tria2;
-
-    DoFHandler<spacedim>      dh1;
-    DoFHandler<dim, spacedim> dh2;
-
-    std::unique_ptr<MappingFEField<dim, spacedim>> mapping2;
-
-    std::vector<IndexSet> owned1;
-    std::vector<IndexSet> owned2;
-
-    std::vector<IndexSet> relevant1;
-    std::vector<IndexSet> relevant2;
-
-    IndexSet owned_tracer_particles;
-    IndexSet relevant_tracer_particles;
+    // For the current implemenation, only `fluid_fe` would be really necessary.
+    // For completeness, and to allow easy extension, we keep also the
+    // `solid_fe` around, which is however initialized to a FE_Nothing finite
+    // element space, i.e., one that has no degrees of freedom.
+    //
+    // We declare both finite element spaces as unique pointers, to allow their
+    // generation after StokesImmersedProblemParameters has been initialized. In
+    // particular, we assume that they are filled only after initial_setup() has
+    // been called.
+    std::unique_ptr<FiniteElement<spacedim>>      fluid_fe;
+    std::unique_ptr<FiniteElement<dim, spacedim>> solid_fe;
+
+    // This is one of the main novelty w.r.t. the tutorial step-60. Here we
+    // assume that both the solid and the fluid are fully distributed
+    // triangulations. This allows the problem to scale to a very large number
+    // of degrees of freedom, at the cost of communicating all the overlapping
+    // regions between non matching triangulations. This is especially tricky,
+    // since we make no assumptions on the relative position or distribution of
+    // the various subdomains. In particular, we assume that ever process owns
+    // only a part of the solid_tria, and only a part of the fluid_tria, not
+    // necessarily in the same physical region, and not necessarily overlapping.
+    //
+    // In order to couple the overlapping regions, we exploit the facilities
+    // implemented in the ParticleHandler class.
+    parallel::distributed::Triangulation<spacedim>      fluid_tria;
+    parallel::distributed::Triangulation<dim, spacedim> solid_tria;
+
+    DoFHandler<spacedim>      fluid_dh;
+    DoFHandler<dim, spacedim> solid_dh;
+
+    std::unique_ptr<MappingFEField<dim, spacedim>> solid_mapping;
+
+    // Similarly to how things are done in step-32, we use a block system to
+    // treat the Stokes part of the problem, and follow very closely what was
+    // done there.
+    std::vector<IndexSet> fluid_owned_dofs;
+    std::vector<IndexSet> solid_owned_dofs;
+
+    std::vector<IndexSet> fluid_relevant_dofs;
+    std::vector<IndexSet> solid_relevant_dofs;
 
     AffineConstraints<double> constraints;
 
@@ -431,9 +544,40 @@ namespace Step70
     LA::MPI::BlockVector       locally_relevant_solution;
     LA::MPI::BlockVector       system_rhs;
 
+    // For every tracer particle, we need to compute the velocity field in its
+    // current position, and update its position using a discrete time stepping
+    // scheme. We do this using distributed linear algebra objects, where the
+    // owner of a particle is set to be equal to the process that generated that
+    // particle at time t=0. This information is stored for every process in the
+    // `owned_tracer_particles` IndexSet, that indicates which particles the
+    // current process owns.
+    //
+    // Once the particles have been distributed around to match the process that
+    // owns the region where the particle lives, we will need read access from
+    // that process on the corresponding velocity field. We achieve this by
+    // filling a read only velocity vector field, that contains the relevant
+    // information in ghost entries. This is achieved using the
+    // `relevant_tracer_particles` IndexSet, that keeps track of how things
+    // change during the simulation, i.e., it keeps track of where particles
+    // that I own have ended up being, and who owns the particles that ended up
+    // in my subdomain.
+    //
+    // While this is not the most efficient strategy, we keep it this way to
+    // illustrate how things would work in a real FSI problem. If a particle
+    // is linked to a specific solid degree of freedom, we are not free to
+    // choose who owns it, and we have to communicate this information around.
+    // We illustrate this here, and show that the communication pattern is
+    // point-to-point, and negligible in terms of total cost of the algorithm.
+    IndexSet owned_tracer_particles;
+    IndexSet relevant_tracer_particles;
+
+    // These vectors are used to store the particles velocities (read-only, with
+    // ghost entries) and their displacement (read/write, no ghost entries).
     LA::MPI::Vector tracer_particle_velocities;
     LA::MPI::Vector relevant_tracer_particle_displacements;
 
+    // We fix once the quadrature formula that is used to integrate the solid
+    // domain.
     std::unique_ptr<Quadrature<dim>> quadrature_formula;
 
     Particles::ParticleHandler<dim, spacedim> tracer_particle_handler;
@@ -450,16 +594,16 @@ namespace Step70
     const StokesImmersedProblemParameters<dim, spacedim> &par)
     : par(par)
     , mpi_communicator(MPI_COMM_WORLD)
-    , tria1(mpi_communicator,
-            typename Triangulation<spacedim>::MeshSmoothing(
-              Triangulation<spacedim>::smoothing_on_refinement |
-              Triangulation<spacedim>::smoothing_on_coarsening))
-    , tria2(mpi_communicator,
-            typename Triangulation<dim, spacedim>::MeshSmoothing(
-              Triangulation<dim, spacedim>::smoothing_on_refinement |
-              Triangulation<dim, spacedim>::smoothing_on_coarsening))
-    , dh1(tria1)
-    , dh2(tria2)
+    , fluid_tria(mpi_communicator,
+                 typename Triangulation<spacedim>::MeshSmoothing(
+                   Triangulation<spacedim>::smoothing_on_refinement |
+                   Triangulation<spacedim>::smoothing_on_coarsening))
+    , solid_tria(mpi_communicator,
+                 typename Triangulation<dim, spacedim>::MeshSmoothing(
+                   Triangulation<dim, spacedim>::smoothing_on_refinement |
+                   Triangulation<dim, spacedim>::smoothing_on_coarsening))
+    , fluid_dh(fluid_tria)
+    , solid_dh(solid_tria)
     , pcout(std::cout,
             (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
     , computing_timer(mpi_communicator,
@@ -472,15 +616,15 @@ namespace Step70
   template <int dim, int spacedim>
   void StokesImmersedProblem<dim, spacedim>::make_grid()
   {
-    GridGenerator::generate_from_name_and_arguments(tria1,
+    GridGenerator::generate_from_name_and_arguments(fluid_tria,
                                                     par.name_of_grid1,
                                                     par.arguments_for_grid1);
-    tria1.refine_global(par.initial_fluid_refinement);
+    fluid_tria.refine_global(par.initial_fluid_refinement);
 
-    GridGenerator::generate_from_name_and_arguments(tria2,
+    GridGenerator::generate_from_name_and_arguments(solid_tria,
                                                     par.name_of_grid2,
                                                     par.arguments_for_grid2);
-    tria2.refine_global(par.initial_solid_refinement);
+    solid_tria.refine_global(par.initial_solid_refinement);
   }
 
   template <int dim, int spacedim>
@@ -503,7 +647,7 @@ namespace Step70
     particles_dof_handler.distribute_dofs(particles_fe);
 
     // Create the particle handler associated with the fluid triangulation
-    tracer_particle_handler.initialize(tria1,
+    tracer_particle_handler.initialize(fluid_tria,
                                        StaticMappingQ1<spacedim>::mapping);
 
 
@@ -511,7 +655,7 @@ namespace Step70
     // The generation of the global bounding boxes requires an all-to-all
     // communication
     auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
-      tria1, IteratorFilters::LocallyOwnedCell());
+      fluid_tria, IteratorFilters::LocallyOwnedCell());
     auto global_bounding_boxes =
       Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
 
@@ -529,11 +673,11 @@ namespace Step70
     relevant_tracer_particles = owned_tracer_particles;
 
     // Now make sure that upon refinement, particles are correctly transferred
-    tria1.signals.pre_distributed_refinement.connect(std::bind(
+    fluid_tria.signals.pre_distributed_refinement.connect(std::bind(
       &Particles::ParticleHandler<spacedim>::register_store_callback_function,
       &tracer_particle_handler));
 
-    tria1.signals.post_distributed_refinement.connect(std::bind(
+    fluid_tria.signals.post_distributed_refinement.connect(std::bind(
       &Particles::ParticleHandler<dim,
                                   spacedim>::register_load_callback_function,
       &tracer_particle_handler,
@@ -546,28 +690,28 @@ namespace Step70
   template <int dim, int spacedim>
   void StokesImmersedProblem<dim, spacedim>::setup_solid_particles()
   {
-    QGauss<dim> quadrature(fe1->degree + 1);
+    QGauss<dim> quadrature(fluid_fe->degree + 1);
     // In codimension one case, we store also the normal, else only the
     // quadrature weight.
     const unsigned int n_properties = (dim == spacedim) ? 1 : spacedim + 1;
-    solid_particle_handler.initialize(tria1,
+    solid_particle_handler.initialize(fluid_tria,
                                       StaticMappingQ1<dim>::mapping,
                                       n_properties);
 
     std::vector<Point<spacedim>> quadrature_points_vec(
-      quadrature.size() * tria2.n_locally_owned_active_cells());
+      quadrature.size() * solid_tria.n_locally_owned_active_cells());
 
     std::vector<std::vector<double>> properties(
-      quadrature.size() * tria2.n_locally_owned_active_cells(),
+      quadrature.size() * solid_tria.n_locally_owned_active_cells(),
       std::vector<double>(n_properties));
 
     UpdateFlags flags = update_JxW_values | update_quadrature_points;
     if (spacedim > dim)
       flags |= update_normal_vectors;
-    FEValues<dim, spacedim> fe_v(*fe2, quadrature, flags);
+    FEValues<dim, spacedim> fe_v(*solid_fe, quadrature, flags);
 
     unsigned int cell_index = 0;
-    for (const auto &cell : dh2.active_cell_iterators())
+    for (const auto &cell : solid_dh.active_cell_iterators())
       if (cell->is_locally_owned())
         {
           fe_v.reinit(cell);
@@ -591,7 +735,7 @@ namespace Step70
     // Distribute the local points to the processor that owns
     // them on the triangulation
     auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
-      tria1, IteratorFilters::LocallyOwnedCell());
+      fluid_tria, IteratorFilters::LocallyOwnedCell());
 
     auto global_bounding_boxes =
       Utilities::MPI::all_gather(mpi_communicator, my_bounding_box);
@@ -603,11 +747,11 @@ namespace Step70
 
 
     // Now make sure that upon refinement, particles are correctly transferred
-    tria1.signals.pre_distributed_refinement.connect(std::bind(
+    fluid_tria.signals.pre_distributed_refinement.connect(std::bind(
       &Particles::ParticleHandler<spacedim>::register_store_callback_function,
       &solid_particle_handler));
 
-    tria1.signals.post_distributed_refinement.connect(std::bind(
+    fluid_tria.signals.post_distributed_refinement.connect(std::bind(
       &Particles::ParticleHandler<dim,
                                   spacedim>::register_load_callback_function,
       &solid_particle_handler,
@@ -622,7 +766,7 @@ namespace Step70
   {
     TimerOutput::Scope t(computing_timer, "initial setup");
 
-    fe1 =
+    fluid_fe =
       std::make_unique<FESystem<spacedim>>(FE_Q<spacedim>(par.velocity_degree),
                                            spacedim,
                                            FE_Q<spacedim>(par.velocity_degree -
@@ -630,8 +774,8 @@ namespace Step70
                                            1);
 
 
-    fe2 = std::make_unique<FE_Nothing<dim, spacedim>>();
-    dh2.distribute_dofs(*fe2);
+    solid_fe = std::make_unique<FE_Nothing<dim, spacedim>>();
+    solid_dh.distribute_dofs(*solid_fe);
     quadrature_formula = std::make_unique<QGauss<dim>>(par.velocity_degree + 1);
   }
 
@@ -642,43 +786,44 @@ namespace Step70
   {
     TimerOutput::Scope t(computing_timer, "setup dofs");
 
-    dh1.distribute_dofs(*fe1);
+    fluid_dh.distribute_dofs(*fluid_fe);
 
     std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
     stokes_sub_blocks[dim] = 1;
-    DoFRenumbering::component_wise(dh1, stokes_sub_blocks);
+    DoFRenumbering::component_wise(fluid_dh, stokes_sub_blocks);
 
     auto dofs_per_block =
-      DoFTools::count_dofs_per_fe_block(dh1, stokes_sub_blocks);
+      DoFTools::count_dofs_per_fe_block(fluid_dh, stokes_sub_blocks);
 
     const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
 
-    pcout << "   Number of degrees of freedom: " << dh1.n_dofs() << " (" << n_u
-          << '+' << n_p << " -- " << solid_particle_handler.n_global_particles()
-          << '+' << tracer_particle_handler.n_global_particles() << ')'
-          << std::endl;
+    pcout << "   Number of degrees of freedom: " << fluid_dh.n_dofs() << " ("
+          << n_u << '+' << n_p << " -- "
+          << solid_particle_handler.n_global_particles() << '+'
+          << tracer_particle_handler.n_global_particles() << ')' << std::endl;
 
-    owned1.resize(2);
-    owned1[0] = dh1.locally_owned_dofs().get_view(0, n_u);
-    owned1[1] = dh1.locally_owned_dofs().get_view(n_u, n_u + n_p);
+    fluid_owned_dofs.resize(2);
+    fluid_owned_dofs[0] = fluid_dh.locally_owned_dofs().get_view(0, n_u);
+    fluid_owned_dofs[1] =
+      fluid_dh.locally_owned_dofs().get_view(n_u, n_u + n_p);
 
     IndexSet locally_relevant_dofs;
-    DoFTools::extract_locally_relevant_dofs(dh1, locally_relevant_dofs);
-    relevant1.resize(2);
-    relevant1[0] = locally_relevant_dofs.get_view(0, n_u);
-    relevant1[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p);
+    DoFTools::extract_locally_relevant_dofs(fluid_dh, locally_relevant_dofs);
+    fluid_relevant_dofs.resize(2);
+    fluid_relevant_dofs[0] = locally_relevant_dofs.get_view(0, n_u);
+    fluid_relevant_dofs[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p);
 
     {
       constraints.reinit(locally_relevant_dofs);
 
       FEValuesExtractors::Vector velocities(0);
-      DoFTools::make_hanging_node_constraints(dh1, constraints);
-      VectorTools::interpolate_boundary_values(dh1,
-                                               0,
-                                               ZeroFunction<spacedim>(spacedim +
-                                                                      1),
-                                               constraints,
-                                               fe1->component_mask(velocities));
+      DoFTools::make_hanging_node_constraints(fluid_dh, constraints);
+      VectorTools::interpolate_boundary_values(
+        fluid_dh,
+        0,
+        ZeroFunction<spacedim>(spacedim + 1),
+        constraints,
+        fluid_fe->component_mask(velocities));
       constraints.close();
     }
 
@@ -697,15 +842,16 @@ namespace Step70
 
       BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
 
-      DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false);
+      DoFTools::make_sparsity_pattern(
+        fluid_dh, coupling, dsp, constraints, false);
 
       SparsityTools::distribute_sparsity_pattern(
         dsp,
-        dh1.compute_locally_owned_dofs_per_processor(),
+        fluid_dh.compute_locally_owned_dofs_per_processor(),
         mpi_communicator,
         locally_relevant_dofs);
 
-      system_matrix.reinit(owned1, dsp, mpi_communicator);
+      system_matrix.reinit(fluid_owned_dofs, dsp, mpi_communicator);
     }
 
     {
@@ -721,18 +867,21 @@ namespace Step70
 
       BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
 
-      DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false);
+      DoFTools::make_sparsity_pattern(
+        fluid_dh, coupling, dsp, constraints, false);
       SparsityTools::distribute_sparsity_pattern(
         dsp,
-        dh1.compute_locally_owned_dofs_per_processor(),
+        fluid_dh.compute_locally_owned_dofs_per_processor(),
         mpi_communicator,
         locally_relevant_dofs);
-      preconditioner_matrix.reinit(owned1, dsp, mpi_communicator);
+      preconditioner_matrix.reinit(fluid_owned_dofs, dsp, mpi_communicator);
     }
 
-    locally_relevant_solution.reinit(owned1, relevant1, mpi_communicator);
-    system_rhs.reinit(owned1, mpi_communicator);
-    solution.reinit(owned1, mpi_communicator);
+    locally_relevant_solution.reinit(fluid_owned_dofs,
+                                     fluid_relevant_dofs,
+                                     mpi_communicator);
+    system_rhs.reinit(fluid_owned_dofs, mpi_communicator);
+    solution.reinit(fluid_owned_dofs, mpi_communicator);
   }
 
 
@@ -747,13 +896,13 @@ namespace Step70
     TimerOutput::Scope t(computing_timer, "Stokes_assembly");
 
 
-    FEValues<spacedim> fe_values(*fe1,
+    FEValues<spacedim> fe_values(*fluid_fe,
                                  *quadrature_formula,
                                  update_values | update_gradients |
                                    update_quadrature_points |
                                    update_JxW_values);
 
-    const unsigned int dofs_per_cell = fe1->dofs_per_cell;
+    const unsigned int dofs_per_cell = fluid_fe->dofs_per_cell;
     const unsigned int n_q_points    = quadrature_formula->size();
 
     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
@@ -771,7 +920,7 @@ namespace Step70
     const FEValuesExtractors::Vector     velocities(0);
     const FEValuesExtractors::Scalar     pressure(spacedim);
 
-    for (const auto &cell : dh1.active_cell_iterators())
+    for (const auto &cell : fluid_dh.active_cell_iterators())
       if (cell->is_locally_owned())
         {
           cell_matrix  = 0;
@@ -805,7 +954,7 @@ namespace Step70
                     }
 
                   const unsigned int component_i =
-                    fe1->system_to_component_index(i).first;
+                    fluid_fe->system_to_component_index(i).first;
                   cell_rhs(i) += fe_values.shape_value(i, q) *
                                  rhs_values[q](component_i) * fe_values.JxW(q);
                 }
@@ -838,19 +987,20 @@ namespace Step70
 
     SolidVelocity<spacedim> solid_velocity(par.angular_velocity);
 
-    std::vector<types::global_dof_index> dof_indices1(fe1->dofs_per_cell);
+    std::vector<types::global_dof_index> dof_indices1(fluid_fe->dofs_per_cell);
 
-    FullMatrix<double>     local_matrix(fe1->dofs_per_cell, fe1->dofs_per_cell);
-    dealii::Vector<double> local_rhs(fe1->dofs_per_cell);
+    FullMatrix<double>     local_matrix(fluid_fe->dofs_per_cell,
+                                    fluid_fe->dofs_per_cell);
+    dealii::Vector<double> local_rhs(fluid_fe->dofs_per_cell);
 
     auto particle = solid_particle_handler.begin();
     while (particle != solid_particle_handler.end())
       {
         local_matrix     = 0;
         local_rhs        = 0;
-        const auto &cell = particle->get_surrounding_cell(tria1);
+        const auto &cell = particle->get_surrounding_cell(fluid_tria);
         const auto &dh_cell =
-          typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &dh1);
+          typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &fluid_dh);
         dh_cell->get_dof_indices(dof_indices1);
 
         const auto pic = solid_particle_handler.particles_in_cell(cell);
@@ -861,23 +1011,24 @@ namespace Step70
             const auto  real_q     = p.get_location();
             const auto  properties = p.get_properties();
             const auto &JxW        = properties[0];
-            for (unsigned int i = 0; i < fe1->dofs_per_cell; ++i)
+            for (unsigned int i = 0; i < fluid_fe->dofs_per_cell; ++i)
               {
-                const auto comp_i = fe1->system_to_component_index(i).first;
+                const auto comp_i =
+                  fluid_fe->system_to_component_index(i).first;
                 if (comp_i < spacedim)
                   {
-                    for (unsigned int j = 0; j < fe1->dofs_per_cell; ++j)
+                    for (unsigned int j = 0; j < fluid_fe->dofs_per_cell; ++j)
                       {
                         const auto comp_j =
-                          fe1->system_to_component_index(j).first;
+                          fluid_fe->system_to_component_index(j).first;
                         if (comp_i == comp_j)
                           local_matrix(i, j) +=
-                            par.penalty_term * fe1->shape_value(i, ref_q) *
-                            fe1->shape_value(j, ref_q) * JxW;
+                            par.penalty_term * fluid_fe->shape_value(i, ref_q) *
+                            fluid_fe->shape_value(j, ref_q) * JxW;
                       }
                     local_rhs(i) += par.penalty_term *
                                     solid_velocity.value(real_q, comp_i) *
-                                    fe1->shape_value(i, ref_q) * JxW;
+                                    fluid_fe->shape_value(i, ref_q) * JxW;
                   }
               }
           }
@@ -948,7 +1099,7 @@ namespace Step70
 
     locally_relevant_solution = solution;
     const double mean_pressure =
-      VectorTools::compute_mean_value(dh1,
+      VectorTools::compute_mean_value(fluid_dh,
                                       QGauss<spacedim>(par.velocity_degree + 2),
                                       locally_relevant_solution,
                                       spacedim);
@@ -959,23 +1110,23 @@ namespace Step70
 
 
   template <int dim, int spacedim>
-  void StokesImmersedProblem<dim, spacedim>::refine_grid()
+  void StokesImmersedProblem<dim, spacedim>::refine_and_transfer()
   {
     TimerOutput::Scope               t(computing_timer, "refine");
     const FEValuesExtractors::Vector velocity(0);
 
-    Vector<float> error_per_cell(tria1.n_active_cells());
-    KellyErrorEstimator<dim>::estimate(dh1,
+    Vector<float> error_per_cell(fluid_tria.n_active_cells());
+    KellyErrorEstimator<dim>::estimate(fluid_dh,
                                        QGauss<dim - 1>(par.velocity_degree + 1),
                                        {},
                                        locally_relevant_solution,
                                        error_per_cell,
-                                       fe1->component_mask(velocity));
+                                       fluid_fe->component_mask(velocity));
 
     if (par.refinement_strategy == "fixed_fraction")
       {
         parallel::distributed::GridRefinement::
-          refine_and_coarsen_fixed_fraction(tria1,
+          refine_and_coarsen_fixed_fraction(fluid_tria,
                                             error_per_cell,
                                             par.refinement_fraction,
                                             par.coarsening_fraction);
@@ -983,22 +1134,22 @@ namespace Step70
     else if (par.refinement_strategy == "fixed_number")
       {
         parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
-          tria1,
+          fluid_tria,
           error_per_cell,
           par.refinement_fraction,
           par.coarsening_fraction,
           par.max_cells);
       }
 
-    for (const auto &cell : tria1.active_cell_iterators())
+    for (const auto &cell : fluid_tria.active_cell_iterators())
       if (cell->refine_flag_set() && cell->level() == par.max_level_refinement)
         cell->clear_refine_flag();
 
     parallel::distributed::SolutionTransfer<dim, LA::MPI::BlockVector> transfer(
-      dh1);
-    tria1.prepare_coarsening_and_refinement();
+      fluid_dh);
+    fluid_tria.prepare_coarsening_and_refinement();
     transfer.prepare_for_coarsening_and_refinement(locally_relevant_solution);
-    tria1.execute_coarsening_and_refinement();
+    fluid_tria.execute_coarsening_and_refinement();
     setup_dofs();
     transfer.interpolate(solution);
     constraints.distribute(solution);
@@ -1023,20 +1174,20 @@ namespace Step70
       DataComponentInterpretation::component_is_scalar);
 
     DataOut<spacedim> data_out;
-    data_out.attach_dof_handler(dh1);
+    data_out.attach_dof_handler(fluid_dh);
     data_out.add_data_vector(locally_relevant_solution,
                              solution_names,
                              DataOut<spacedim>::type_dof_data,
                              data_component_interpretation);
 
     LA::MPI::BlockVector interpolated;
-    interpolated.reinit(owned1, MPI_COMM_WORLD);
-    VectorTools::interpolate(dh1,
+    interpolated.reinit(fluid_owned_dofs, MPI_COMM_WORLD);
+    VectorTools::interpolate(fluid_dh,
                              ConstantFunction<spacedim>(1.0, spacedim + 1),
                              interpolated);
 
-    LA::MPI::BlockVector interpolated_relevant(owned1,
-                                               relevant1,
+    LA::MPI::BlockVector interpolated_relevant(fluid_owned_dofs,
+                                               fluid_relevant_dofs,
                                                MPI_COMM_WORLD);
     interpolated_relevant = interpolated;
     {
@@ -1049,9 +1200,9 @@ namespace Step70
     }
 
 
-    Vector<float> subdomain(tria1.n_active_cells());
+    Vector<float> subdomain(fluid_tria.n_active_cells());
     for (unsigned int i = 0; i < subdomain.size(); ++i)
-      subdomain(i) = tria1.locally_owned_subdomain();
+      subdomain(i) = fluid_tria.locally_owned_subdomain();
     data_out.add_data_vector(subdomain, "subdomain");
 
     data_out.build_patches();
@@ -1134,7 +1285,7 @@ namespace Step70
           }
         {
           TimerOutput::Scope t(computing_timer, "Set tracer particle motion");
-          interpolate_field_on_particles(dh1,
+          interpolate_field_on_particles(fluid_dh,
                                          tracer_particle_handler,
                                          locally_relevant_solution,
                                          tracer_particle_velocities,
@@ -1160,7 +1311,7 @@ namespace Step70
         assemble_nitche_restriction();
         solve();
 
-        if (cycle % par.mod_output == 0)
+        if (cycle % par.output_frequency == 0)
           {
             static unsigned int output_cycle = 0;
             output_results(output_cycle, time);
@@ -1180,11 +1331,104 @@ namespace Step70
             }
             ++output_cycle;
           }
-        if (cycle % par.mod_refinement == 0 &&
+        if (cycle % par.refinement_frequency == 0 &&
             cycle != par.number_of_time_steps - 1)
-          refine_grid();
+          refine_and_transfer();
       }
   }
+
+  template <int dim, int spacedim>
+  StokesImmersedProblemParameters<dim,
+                                  spacedim>::StokesImmersedProblemParameters()
+    : ParameterAcceptor("Stokes Immersed Problem/")
+    , rhs("Right hand side", spacedim + 1)
+    , angular_velocity("Angular velocity", spacedim == 3 ? spacedim : 1)
+  {
+    // We split the parameters in various cathegories, by putting them in
+    // different sections of the ParameterHandler class. We begin by declaring
+    // all the global parameters used by StokesImmersedProblem in the global
+    // scope:
+    add_parameter(
+      "Velocity degree", velocity_degree, "", this->prm, Patterns::Integer(1));
+
+    add_parameter("Number of time steps", number_of_time_steps);
+    add_parameter("Output frequency", output_frequency);
+
+    add_parameter("Final time", final_time);
+
+    add_parameter("Viscosity", viscosity);
+
+    add_parameter("Nitsche penalty term", penalty_term);
+
+    add_parameter("Initial fluid refinement",
+                  initial_fluid_refinement,
+                  "Initial mesh refinement used for the fluid domain Omega");
+
+    add_parameter("Initial solid refinement",
+                  initial_solid_refinement,
+                  "Initial mesh refinement used for the solid domain Gamma");
+
+    add_parameter(
+      "Particle insertion refinement",
+      particle_insertion_refinement,
+      "Refinement of the volumetric mesh used to insert the particles");
+
+    add_parameter(
+      "Homogeneous Dirichlet boundary ids",
+      homogeneous_dirichlet_ids,
+      "Boundary Ids over which homogeneous Dirichlet boundary conditions are applied");
+
+    // Next section is dedicated to the parameters used to create the various
+    // grids. We will need three different triangulations: `Grid one` is used
+    // to define the fluid domain, `Grid two` defines the solid domain, and
+    // `Particle grid` is used to distribute some tracer particles, that are
+    // advected with the velocity and only used as passive tracers.
+    enter_my_subsection(this->prm);
+    this->prm.enter_subsection("Grid generation");
+    this->prm.add_parameter("Grid one generator", name_of_grid1);
+    this->prm.add_parameter("Grid one generator arguments",
+                            arguments_for_grid1);
+
+    this->prm.add_parameter("Grid two generator", name_of_grid2);
+    this->prm.add_parameter("Grid two generator arguments",
+                            arguments_for_grid2);
+
+    this->prm.add_parameter("Particle grid generator", name_of_particle_grid);
+    this->prm.add_parameter("Particle grid generator arguments",
+                            arguments_for_particle_grid);
+    this->prm.leave_subsection();
+
+    leave_my_subsection(this->prm);
+
+
+
+    enter_my_subsection(this->prm);
+    this->prm.enter_subsection("Refinement and remeshing");
+    this->prm.add_parameter("Refinement step frequency", refinement_frequency);
+    this->prm.add_parameter("Refinement maximal level", max_level_refinement);
+    this->prm.add_parameter("Refinement strategy",
+                            refinement_strategy,
+                            "",
+                            Patterns::Selection("fixed_fraction|fixed_number"));
+    this->prm.add_parameter("Refinement coarsening fraction",
+                            coarsening_fraction);
+    this->prm.add_parameter("Refinement fraction", refinement_fraction);
+    this->prm.add_parameter("Maximum number of cells", max_cells);
+
+    this->prm.leave_subsection();
+    leave_my_subsection(this->prm);
+
+    // correct the default dimension for the functions
+    rhs.declare_parameters_call_back.connect([&]() {
+      Functions::ParsedFunction<spacedim>::declare_parameters(this->prm,
+                                                              spacedim + 1);
+    });
+    angular_velocity.declare_parameters_call_back.connect([&]() {
+      Functions::ParsedFunction<spacedim>::declare_parameters(
+        this->prm, spacedim == 3 ? spacedim : 1);
+    });
+  }
+
 } // namespace Step70
 
 
@@ -1199,9 +1443,6 @@ int main(int argc, char *argv[])
       Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
 
       StokesImmersedProblemParameters<2> par;
-      par.declare_all_parameters();
-      std::ofstream out("default.prm");
-      par.prm.print_parameters(out, ParameterHandler::ShortText);
       ParameterAcceptor::initialize("parameters.prm", "used_parameters.prm");
 
       StokesImmersedProblem<2> problem(par);

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.