]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Step-70
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 21 Nov 2019 19:28:22 +0000 (20:28 +0100)
committerMatthias Maier <tamiko@43-1.org>
Fri, 15 May 2020 03:16:35 +0000 (22:16 -0500)
Work in progress for step-70.

Co-authored-by: Bruno Blais <blais.bruno@gmail.com>
examples/step-70/CMakeLists.txt [new file with mode: 0644]
examples/step-70/step-70.cc [new file with mode: 0644]

diff --git a/examples/step-70/CMakeLists.txt b/examples/step-70/CMakeLists.txt
new file mode 100644 (file)
index 0000000..760967f
--- /dev/null
@@ -0,0 +1,62 @@
+##
+#  CMake script
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-70")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+FIND_PACKAGE(deal.II 9.2.0 QUIET
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+#
+# Are all dependencies fulfilled?
+#
+IF(NOT ((DEAL_II_WITH_PETSC AND NOT DEAL_II_PETSC_WITH_COMPLEX) OR DEAL_II_WITH_TRILINOS) OR NOT DEAL_II_WITH_P4EST) # keep in one line
+  MESSAGE(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following options:
+    DEAL_II_WITH_PETSC = ON
+    DEAL_II_PETSC_WITH_COMPLEX = OFF
+    DEAL_II_WITH_P4EST = ON
+or
+    DEAL_II_WITH_TRILINOS = ON
+    DEAL_II_WITH_P4EST = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options
+    DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC}
+    DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX}
+    DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST}
+    DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS}
+which conflict with the requirements.
+One or both of the aforementioned combinations of prerequisites are not met by your installation, but at least one is required for this tutorial step."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+SET(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.vtk *.ucd *.d2 *.vtu *.pvtu)
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-70/step-70.cc b/examples/step-70/step-70.cc
new file mode 100644 (file)
index 0000000..fa0d265
--- /dev/null
@@ -0,0 +1,1137 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2020 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Authors: Luca Heltai, Bruno Blais, 2019
+ */
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/lac/block_linear_operator.h>
+#include <deal.II/lac/generic_linear_algebra.h>
+#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
+{
+#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
+  !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
+  using namespace dealii::LinearAlgebraPETSc;
+#  define USE_PETSC_LA
+#elif defined(DEAL_II_WITH_TRILINOS)
+  using namespace dealii::LinearAlgebraTrilinos;
+#else
+#  error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
+#endif
+} // namespace LA
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/parameter_acceptor.h>
+#include <deal.II/base/parsed_function.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.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/fe/mapping_fe_field.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.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/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/solver_minres.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/particles/generators.h>
+#include <deal.II/particles/particle_handler.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <memory>
+
+namespace Step70
+{
+  using namespace dealii;
+
+  template <int dim,
+            int spacedim,
+            typename InputVectorType,
+            typename OutputVectorType>
+  void interpolate_field_on_particles(
+    const DoFHandler<dim, spacedim> &                field_dh,
+    const Particles::ParticleHandler<dim, spacedim> &particle_handler,
+    const InputVectorType &                          field_vector,
+    OutputVectorType &                               interpolated_field,
+    const ComponentMask &                            field_comps)
+  {
+    if (particle_handler.n_locally_owned_particles() == 0)
+      {
+        interpolated_field.compress(VectorOperation::add);
+        return; // nothing else to do here
+      }
+
+    const auto &tria     = field_dh.get_triangulation();
+    const auto &fe       = field_dh.get_fe();
+    auto        particle = particle_handler.begin();
+
+    // Take care of components
+    const ComponentMask comps =
+      (field_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
+                                 field_comps);
+    AssertDimension(comps.size(), fe.n_components());
+    const auto n_comps = comps.n_selected_components();
+
+    AssertDimension(field_vector.size(), field_dh.n_dofs());
+    AssertDimension(interpolated_field.size(),
+                    particle_handler.get_next_free_particle_index() * n_comps);
+    // Add check on locally owned indices
+
+    // Global to local indices
+    std::vector<unsigned int> space_gtl(fe.n_components(),
+                                        numbers::invalid_unsigned_int);
+    for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
+      if (comps[i])
+        space_gtl[i] = j++;
+
+    std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
+
+    while (particle != particle_handler.end())
+      {
+        const auto &cell = particle->get_surrounding_cell(tria);
+        const auto &dh_cell =
+          typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &field_dh);
+        dh_cell->get_dof_indices(dof_indices);
+        const auto pic = particle_handler.particles_in_cell(cell);
+        Assert(pic.begin() == particle, ExcInternalError());
+        for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
+          {
+            const auto &reference_location = particle->get_reference_location();
+
+            const auto id = particle->get_id();
+
+            for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
+              {
+                const auto comp_j =
+                  space_gtl[fe.system_to_component_index(j).first];
+                if (comp_j != numbers::invalid_unsigned_int)
+                  interpolated_field[id * n_comps + comp_j] +=
+                    fe.shape_value(j, reference_location) *
+                    field_vector(dof_indices[j]);
+              }
+          }
+      }
+    interpolated_field.compress(VectorOperation::add);
+  }
+
+  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("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);
+
+      // 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);
+      });
+    }
+
+    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 =
+      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";
+
+    mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>> rhs;
+    mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
+      angular_velocity;
+  }; // namespace Step70
+
+
+  template <int spacedim>
+  class SolidVelocity : public Function<spacedim>
+  {
+  public:
+    SolidVelocity(const Functions::ParsedFunction<spacedim> &angular_velocity)
+      : angular_velocity(angular_velocity)
+    {}
+
+    virtual double value(const Point<spacedim> &p,
+                         unsigned int           component = 0) const
+    {
+      Tensor<1, spacedim> velocity;
+      if (spacedim == 3)
+        {
+          Tensor<1, spacedim> omega;
+          for (unsigned int i = 0; i < spacedim; ++i)
+            omega[i] = angular_velocity.value(p, i);
+
+          velocity = cross_product_3d(p, omega);
+        }
+
+      if (spacedim == 2)
+        {
+          double omega = angular_velocity.value(p, 0);
+
+          velocity[0] = -omega * p[1];
+          velocity[1] = omega * p[0];
+        }
+
+      return velocity[component];
+    }
+
+  private:
+    const Functions::ParsedFunction<spacedim> &angular_velocity;
+  };
+
+
+  template <int spacedim>
+  class SolidDisplacement : public Function<spacedim>
+  {
+  public:
+    SolidDisplacement(
+      const Functions::ParsedFunction<spacedim> &angular_velocity,
+      const double                               time_step)
+      : Function<spacedim>(spacedim)
+      , angular_velocity(angular_velocity)
+      , time_step(time_step)
+    {}
+
+    virtual double value(const Point<spacedim> &p,
+                         unsigned int           component = 0) const
+    {
+      Tensor<1, spacedim> displacement;
+
+      double dtheta = angular_velocity.value(p, 0) * time_step;
+
+      displacement[0] = std::cos(dtheta) * p[0] - std::sin(dtheta) * p[1];
+      displacement[1] = std::sin(dtheta) * p[0] + std::cos(dtheta) * p[1];
+
+      return displacement[component];
+    }
+
+    void set_time_step(const double new_time_step)
+    {
+      time_step = new_time_step;
+    }
+
+  private:
+    const Functions::ParsedFunction<spacedim> &angular_velocity;
+    double                                     time_step;
+  };
+
+  template <int dim, int spacedim = dim>
+  class StokesImmersedProblem
+  {
+  public:
+    StokesImmersedProblem(
+      const StokesImmersedProblemParameters<dim, spacedim> &par);
+
+    void run();
+
+  private:
+    void make_grid();
+    void setup_tracer_particles();
+    void setup_solid_particles();
+    void setup_system();
+    void assemble_stokes_system();
+    void assemble_nitche_restriction();
+    void solve();
+    void refine_grid();
+    void output_results(const unsigned int cycle) const;
+
+    void
+    output_particles(const Particles::ParticleHandler<dim, spacedim> &particles,
+                     std::string                                      fprefix,
+                     const unsigned int iter) const;
+
+    const StokesImmersedProblemParameters<dim, spacedim> &par;
+
+    MPI_Comm mpi_communicator;
+
+    std::unique_ptr<FESystem<spacedim>>      fe1;
+    std::unique_ptr<FESystem<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;
+
+    AffineConstraints<double> constraints;
+
+    LA::MPI::BlockSparseMatrix system_matrix;
+    LA::MPI::BlockSparseMatrix coupling_matrix;
+
+    LA::MPI::BlockSparseMatrix preconditioner_matrix;
+    LA::MPI::BlockVector       solution;
+    LA::MPI::BlockVector       locally_relevant_solution;
+    LA::MPI::BlockVector       system_rhs;
+
+    LA::MPI::Vector tracer_particle_velocities;
+    LA::MPI::Vector relevant_tracer_particle_displacements;
+
+    std::unique_ptr<Quadrature<dim>> quadrature_formula;
+
+    //    std::unique_ptr<NonMatching::DoFHandlerCoupling<dim, spacedim>>
+    //      dof_coupling;
+
+    Particles::ParticleHandler<dim, spacedim> tracer_particle_handler;
+    Particles::ParticleHandler<dim, spacedim> solid_particle_handler;
+
+    ConditionalOStream  pcout;
+    mutable TimerOutput computing_timer;
+  };
+
+
+
+  template <int dim, int spacedim>
+  StokesImmersedProblem<dim, spacedim>::StokesImmersedProblem(
+    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)
+    , pcout(std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+    , computing_timer(mpi_communicator,
+                      pcout,
+                      TimerOutput::summary,
+                      TimerOutput::wall_times)
+  {}
+
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::make_grid()
+  {
+    GridGenerator::generate_from_name_and_arguments(tria1,
+                                                    par.name_of_grid1,
+                                                    par.arguments_for_grid1);
+    tria1.refine_global(par.initial_fluid_refinement);
+
+    GridGenerator::generate_from_name_and_arguments(tria2,
+                                                    par.name_of_grid2,
+                                                    par.arguments_for_grid2);
+    tria2.refine_global(par.initial_solid_refinement);
+  }
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::setup_tracer_particles()
+  {
+    // Generate a triangulation that will be used to decide the position
+    // of the particles to insert. In this case we choose an hyper_ball, a
+    // circle (spacedim==2) or a sphere (spacedim==3) filled with particles are
+    // the position of the support points of the triangulation
+    parallel::distributed::Triangulation<spacedim> particle_insert_tria(
+      mpi_communicator);
+    GridGenerator::generate_from_name_and_arguments(
+      particle_insert_tria,
+      par.name_of_particle_grid,
+      par.arguments_for_particle_grid);
+    particle_insert_tria.refine_global(par.particle_insertion_refinement);
+
+    // Generate the support point on the triangulation that will be used as
+    // particle insertion point
+    DoFHandler<dim, spacedim> particles_dof_handler(particle_insert_tria);
+    FE_Q<dim, spacedim>       particles_fe(1);
+    particles_dof_handler.distribute_dofs(particles_fe);
+
+    // Create the particle handler associated with the fluid triangulation
+    tracer_particle_handler.initialize(tria1,
+                                       StaticMappingQ1<spacedim>::mapping);
+
+
+    // Generate the necessary local and global bounding boxes for the generator.
+    // 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());
+    auto global_bounding_boxes =
+      Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
+
+
+    // Finally generate the particles from the support point of the
+    // particle_insert_tria triangulation
+    Particles::Generators::dof_support_points(particles_dof_handler,
+                                              global_bounding_boxes,
+                                              tracer_particle_handler);
+
+    owned_tracer_particles =
+      tracer_particle_handler.locally_relevant_ids().tensor_product(
+        complete_index_set(spacedim));
+
+    relevant_tracer_particles = owned_tracer_particles;
+  }
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::setup_solid_particles()
+  {
+    QGauss<dim> quadrature(fe1->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,
+                                      StaticMappingQ1<dim>::mapping,
+                                      n_properties);
+
+    std::vector<Point<spacedim>> quadrature_points_vec(
+      quadrature.size() * tria2.n_locally_owned_active_cells());
+
+    std::vector<std::vector<double>> properties(
+      quadrature.size() * tria2.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);
+
+    unsigned int cell_index = 0;
+    for (const auto &cell : dh2.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          fe_v.reinit(cell);
+          const auto &points = fe_v.get_quadrature_points();
+          const auto &JxW    = fe_v.get_JxW_values();
+
+          for (unsigned int q = 0; q < points.size(); ++q)
+            {
+              const auto i             = cell_index * points.size() + q;
+              quadrature_points_vec[i] = points[q];
+              properties[i][0]         = JxW[q];
+              if (dim < spacedim)
+                for (unsigned int d = 0; d < spacedim; ++d)
+                  {
+                    properties[i][d + 1] = fe_v.normal_vector(q)[d];
+                  }
+            }
+          ++cell_index;
+        }
+
+    // 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());
+
+    auto global_bounding_boxes =
+      Utilities::MPI::all_gather(mpi_communicator, my_bounding_box);
+
+    auto cpu_to_index =
+      solid_particle_handler.insert_global_particles(quadrature_points_vec,
+                                                     global_bounding_boxes,
+                                                     properties);
+  }
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::setup_system()
+  {
+    TimerOutput::Scope t(computing_timer, "setup");
+
+    fe1 =
+      std::make_unique<FESystem<spacedim>>(FE_Q<spacedim>(par.velocity_degree),
+                                           spacedim,
+                                           FE_Q<spacedim>(par.velocity_degree -
+                                                          1),
+                                           1);
+
+    fe2 = std::make_unique<FESystem<dim, spacedim>>(
+      FE_Q<dim, spacedim>(par.velocity_degree), spacedim);
+
+    quadrature_formula = std::make_unique<QGauss<dim>>(par.velocity_degree + 1);
+
+    dh1.distribute_dofs(*fe1);
+    dh2.distribute_dofs(*fe2);
+
+    std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
+    stokes_sub_blocks[dim] = 1;
+    DoFRenumbering::component_wise(dh1, stokes_sub_blocks);
+
+    auto dofs_per_block =
+      DoFTools::count_dofs_per_fe_block(dh1, 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 << ')' << 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);
+
+    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);
+
+    {
+      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));
+      constraints.close();
+    }
+
+    {
+      system_matrix.clear();
+
+      Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
+      for (unsigned int c = 0; c < dim + 1; ++c)
+        for (unsigned int d = 0; d < dim + 1; ++d)
+          if (c == dim && d == dim)
+            coupling[c][d] = DoFTools::none;
+          else if (c == dim || d == dim || c == d)
+            coupling[c][d] = DoFTools::always;
+          else
+            coupling[c][d] = DoFTools::none;
+
+      BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
+
+      DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false);
+
+      SparsityTools::distribute_sparsity_pattern(
+        dsp,
+        dh1.compute_locally_owned_dofs_per_processor(),
+        mpi_communicator,
+        locally_relevant_dofs);
+
+      system_matrix.reinit(owned1, dsp, mpi_communicator);
+    }
+
+    {
+      preconditioner_matrix.clear();
+
+      Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
+      for (unsigned int c = 0; c < dim + 1; ++c)
+        for (unsigned int d = 0; d < dim + 1; ++d)
+          if (c == dim && d == dim)
+            coupling[c][d] = DoFTools::always;
+          else
+            coupling[c][d] = DoFTools::none;
+
+      BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
+
+      DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false);
+      SparsityTools::distribute_sparsity_pattern(
+        dsp,
+        dh1.compute_locally_owned_dofs_per_processor(),
+        mpi_communicator,
+        locally_relevant_dofs);
+      preconditioner_matrix.reinit(owned1, dsp, mpi_communicator);
+    }
+
+    locally_relevant_solution.reinit(owned1, relevant1, mpi_communicator);
+    system_rhs.reinit(owned1, mpi_communicator);
+    solution.reinit(owned1, mpi_communicator);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::assemble_stokes_system()
+  {
+    system_matrix         = 0;
+    preconditioner_matrix = 0;
+    system_rhs            = 0;
+
+    TimerOutput::Scope t(computing_timer, "Stokes_assembly");
+
+
+    FEValues<spacedim> fe_values(*fe1,
+                                 *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 n_q_points    = quadrature_formula->size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    FullMatrix<double> cell_matrix2(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<Vector<double>> rhs_values(n_q_points,
+                                           Vector<double>(spacedim + 1));
+
+    std::vector<Tensor<2, spacedim>> grad_phi_u(dofs_per_cell);
+    std::vector<double>              div_phi_u(dofs_per_cell);
+    std::vector<double>              phi_p(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+    const FEValuesExtractors::Vector     velocities(0);
+    const FEValuesExtractors::Scalar     pressure(spacedim);
+
+    for (const auto &cell : dh1.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix  = 0;
+          cell_matrix2 = 0;
+          cell_rhs     = 0;
+
+          fe_values.reinit(cell);
+          par.rhs.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)
+                {
+                  grad_phi_u[k] = fe_values[velocities].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 < dofs_per_cell; ++j)
+                    {
+                      cell_matrix(i, j) +=
+                        (par.viscosity *
+                           scalar_product(grad_phi_u[i], grad_phi_u[j]) -
+                         div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
+                        fe_values.JxW(q);
+
+                      cell_matrix2(i, j) += 1.0 / par.viscosity * phi_p[i] *
+                                            phi_p[j] * fe_values.JxW(q);
+                    }
+
+                  const unsigned int component_i =
+                    fe1->system_to_component_index(i).first;
+                  cell_rhs(i) += fe_values.shape_value(i, q) *
+                                 rhs_values[q](component_i) * fe_values.JxW(q);
+                }
+            }
+
+
+          cell->get_dof_indices(local_dof_indices);
+          constraints.distribute_local_to_global(cell_matrix,
+                                                 cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix,
+                                                 system_rhs);
+
+          constraints.distribute_local_to_global(cell_matrix2,
+                                                 local_dof_indices,
+                                                 preconditioner_matrix);
+        }
+
+    system_matrix.compress(VectorOperation::add);
+    preconditioner_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::assemble_nitche_restriction()
+  {
+    TimerOutput::Scope t(computing_timer, "Nitsche_assembly");
+
+    SolidVelocity<spacedim> solid_velocity(par.angular_velocity);
+
+    std::vector<types::global_dof_index> dof_indices1(fe1->dofs_per_cell);
+
+    FullMatrix<double>     local_matrix(fe1->dofs_per_cell, fe1->dofs_per_cell);
+    dealii::Vector<double> local_rhs(fe1->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 &dh_cell =
+          typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &dh1);
+        dh_cell->get_dof_indices(dof_indices1);
+
+        const auto pic = solid_particle_handler.particles_in_cell(cell);
+        Assert(pic.begin() == particle, ExcInternalError());
+        for (const auto &p : pic)
+          {
+            const auto  ref_q      = p.get_reference_location();
+            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)
+              {
+                const auto comp_i = fe1->system_to_component_index(i).first;
+                if (comp_i < spacedim)
+                  {
+                    for (unsigned int j = 0; j < fe1->dofs_per_cell; ++j)
+                      {
+                        const auto comp_j =
+                          fe1->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;
+                      }
+                    local_rhs(i) += par.penalty_term *
+                                    solid_velocity.value(real_q, comp_i) *
+                                    fe1->shape_value(i, ref_q) * JxW;
+                  }
+              }
+          }
+        constraints.distribute_local_to_global(
+          local_matrix, local_rhs, dof_indices1, system_matrix, system_rhs);
+        particle = pic.end();
+      }
+  }
+
+
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::solve()
+  {
+    TimerOutput::Scope t(computing_timer, "solve");
+
+    LA::MPI::PreconditionAMG prec_A;
+    {
+      LA::MPI::PreconditionAMG::AdditionalData data;
+
+#ifdef USE_PETSC_LA
+      data.symmetric_operator = true;
+#endif
+      prec_A.initialize(system_matrix.block(0, 0), data);
+    }
+
+    LA::MPI::PreconditionAMG prec_S;
+    {
+      LA::MPI::PreconditionAMG::AdditionalData data;
+
+#ifdef USE_PETSC_LA
+      data.symmetric_operator = true;
+#endif
+      prec_S.initialize(preconditioner_matrix.block(1, 1), data);
+    }
+
+    const auto A = linear_operator<LA::MPI::Vector>(system_matrix.block(0, 0));
+    const auto amgA = linear_operator(A, prec_A);
+
+    const auto S =
+      linear_operator<LA::MPI::Vector>(preconditioner_matrix.block(1, 1));
+    const auto amgS = linear_operator(S, prec_S);
+
+    ReductionControl          inner_solver_control(10,
+                                          1e-8 * system_rhs.l2_norm(),
+                                          1.e-2);
+    SolverCG<LA::MPI::Vector> cg(inner_solver_control);
+
+    const auto invS = inverse_operator(S, cg, amgS);
+
+    const auto P =
+      block_diagonal_operator<2, LA::MPI::BlockVector>({amgA, amgS});
+
+    SolverControl solver_control(system_matrix.m(),
+                                 1e-10 * system_rhs.l2_norm());
+
+    SolverMinRes<LA::MPI::BlockVector> solver(solver_control);
+
+    constraints.set_zero(solution);
+
+    solver.solve(system_matrix, solution, system_rhs, P);
+
+
+    pcout << "   Solved in " << solver_control.last_step() << " iterations."
+          << std::endl;
+
+    constraints.distribute(solution);
+
+    locally_relevant_solution = solution;
+    const double mean_pressure =
+      VectorTools::compute_mean_value(dh1,
+                                      QGauss<spacedim>(par.velocity_degree + 2),
+                                      locally_relevant_solution,
+                                      spacedim);
+    solution.block(1).add(-mean_pressure);
+    locally_relevant_solution.block(1) = solution.block(1);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::refine_grid()
+  {
+    TimerOutput::Scope t(computing_timer, "refine");
+
+    tria1.refine_global();
+  }
+
+
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::output_results(
+    const unsigned int cycle) const
+  {
+    TimerOutput::Scope       t(computing_timer, "Output fluid");
+    std::vector<std::string> solution_names(spacedim, "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<spacedim> data_out;
+    data_out.attach_dof_handler(dh1);
+    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,
+                             ConstantFunction<spacedim>(1.0, spacedim + 1),
+                             interpolated);
+
+    LA::MPI::BlockVector interpolated_relevant(owned1,
+                                               relevant1,
+                                               MPI_COMM_WORLD);
+    interpolated_relevant = interpolated;
+    {
+      std::vector<std::string> solution_names(dim, "ref_u");
+      solution_names.emplace_back("ref_p");
+      data_out.add_data_vector(interpolated_relevant,
+                               solution_names,
+                               DataOut<spacedim>::type_dof_data,
+                               data_component_interpretation);
+    }
+
+
+    Vector<float> subdomain(tria1.n_active_cells());
+    for (unsigned int i = 0; i < subdomain.size(); ++i)
+      subdomain(i) = tria1.locally_owned_subdomain();
+    data_out.add_data_vector(subdomain, "subdomain");
+
+    data_out.build_patches();
+
+    const std::string filename =
+      ("solution-" + Utilities::int_to_string(cycle, 2) + "." +
+       Utilities::int_to_string(tria1.locally_owned_subdomain(), 4));
+    std::ofstream output((filename + ".vtu"));
+    data_out.write_vtu(output);
+
+    if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+      {
+        std::vector<std::string> filenames;
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          filenames.push_back("solution-" + Utilities::int_to_string(cycle, 2) +
+                              "." + Utilities::int_to_string(i, 4) + ".vtu");
+
+        std::ofstream master_output(
+          "solution-" + Utilities::int_to_string(cycle, 2) + ".pvtu");
+        data_out.write_pvtu_record(master_output, filenames);
+      }
+  }
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::output_particles(
+    const Particles::ParticleHandler<dim, spacedim> &particles,
+    std::string                                      fprefix,
+    const unsigned int                               iter) const
+  {
+    Particles::DataOut<dim, spacedim> particles_out;
+    particles_out.build_patches(particles);
+    const std::string filename =
+      (fprefix + "-" + Utilities::int_to_string(iter, 2) + "." +
+       Utilities::int_to_string(tria1.locally_owned_subdomain(), 4));
+    std::ofstream output((filename + ".vtu"));
+    particles_out.write_vtu(output);
+
+    if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+      {
+        std::vector<std::string> filenames;
+        for (unsigned int i = 0;
+             i < Utilities::MPI::n_mpi_processes(mpi_communicator);
+             ++i)
+          filenames.push_back(fprefix + "-" +
+                              Utilities::int_to_string(iter, 2) + "." +
+                              Utilities::int_to_string(i, 4) + ".vtu");
+
+        std::ofstream master_output(
+          fprefix + "-" + Utilities::int_to_string(iter, 2) + ".pvtu");
+        particles_out.write_pvtu_record(master_output, filenames);
+      }
+  }
+
+
+  template <int dim, int spacedim>
+  void StokesImmersedProblem<dim, spacedim>::run()
+  {
+#ifdef USE_PETSC_LA
+    pcout << "Running using PETSc." << std::endl;
+#else
+    pcout << "Running using Trilinos." << std::endl;
+#endif
+
+    ComponentMask velocity_mask(spacedim + 1, true);
+    velocity_mask.set(spacedim, false);
+
+    const double time_step = par.final_time / (par.number_of_time_steps - 1);
+    double       time      = 0;
+    for (unsigned int cycle = 0; cycle < par.number_of_time_steps;
+         ++cycle, time += time_step)
+      {
+        par.set_time(time);
+        pcout << "Cycle " << cycle << ':' << std::endl
+              << "Time : " << time << ", time step: " << time_step << std::endl;
+
+        if (cycle == 0)
+          {
+            make_grid();
+            setup_system();
+            setup_tracer_particles();
+            setup_solid_particles();
+            tracer_particle_velocities.reinit(owned_tracer_particles,
+                                              mpi_communicator);
+          }
+        else
+          {
+            TimerOutput::Scope t(computing_timer,
+                                 "Set solid particle position");
+
+            SolidDisplacement<spacedim> solid_displacement(par.angular_velocity,
+                                                           time_step);
+            solid_particle_handler.set_particle_positions(solid_displacement,
+                                                          false);
+          }
+        {
+          TimerOutput::Scope t(computing_timer, "Set tracer particle motion");
+          interpolate_field_on_particles(dh1,
+                                         tracer_particle_handler,
+                                         locally_relevant_solution,
+                                         tracer_particle_velocities,
+                                         velocity_mask);
+
+          tracer_particle_velocities *= time_step;
+
+          relevant_tracer_particles =
+            tracer_particle_handler.locally_relevant_ids().tensor_product(
+              complete_index_set(spacedim));
+
+          relevant_tracer_particle_displacements.reinit(
+            owned_tracer_particles,
+            relevant_tracer_particles,
+            mpi_communicator);
+
+          relevant_tracer_particle_displacements = tracer_particle_velocities;
+
+          tracer_particle_handler.set_particle_positions(
+            relevant_tracer_particle_displacements);
+        }
+        assemble_stokes_system();
+        assemble_nitche_restriction();
+        solve();
+
+        if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32)
+          {
+            output_results(cycle);
+            {
+              TimerOutput::Scope t(computing_timer, "Output tracer particles");
+              output_particles(tracer_particle_handler, "tracer", cycle);
+            }
+            {
+              TimerOutput::Scope t(computing_timer, "Output solid particles");
+              output_particles(solid_particle_handler, "solid", cycle);
+            }
+          }
+      }
+  }
+} // namespace Step70
+
+
+
+int main(int argc, char *argv[])
+{
+  using namespace Step70;
+  using namespace dealii;
+  deallog.depth_console(1);
+  try
+    {
+      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);
+      problem.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}

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.