to a severe bottleneck in terms of memory. Consequently, an alternative strategy is sought
in this step.
-When looking at the formulation for both the penalization approach ($\mathcal{L}^2$ or $\mathcal{H}^1$)
-and the Nitsche method, we can come to the conclusion that we only require limited information
-related to $\Gamma$ on its quadrature points that is:
-- Their reference location with regard to the embedding fluid cell;
-- The weights times the Jacobian of the transformation, i.e. `JxW`.
+<h4>Using particles to track Γ</h4>
+
+Remember that for both the penalization approach ($\mathcal{L}^2$ or $\mathcal{H}^1$)
+and the Nitsche method, we want to compute integrals that are approximated by
+the quadrature. That is, we need to compute
+\f[
+\beta (\textbf{v},\textbf{u})_{\Gamma} =
+\sum_{K\in \Gamma} \sum_{i=1}^{n_q} \big(\hat{\textbf{u}}(\hat x_i) (\textbf{v} \circ F_{K}) (\hat x_i) J_K (\hat x_i) w_i \big)
+\f]
+If you followed the discussion above, then you will recall that $\textbf{u}$
+and $\textbf{v}$ are shape functions defined on the fluid mesh.
+The only things defined on the solid mesh are:
+$F_K(\hat x_i)$, which is the location of a quadrature point on a solid cell that
+is part of $\Gamma$, $J_K$ is the determinant of its Jacobian, and $w_i$ the corresponding
+quadrature weight.
+
+The important part to realize is now this: $w_i$ is a property of
+the quadrature formula and does not change with time. Furthermore,
+the Jacobian matrix of $F_K$ itself changes as the solid obstacle
+moves around in the fluid, but because the solid is considered
+non-deforming (it only translates and rotates, but doesn't dilate),
+the determinant of the Jacobian remains constant. As a consequence,
+the product $J_K(\hat x_i) w_i$ (which we typically denote by `JxW`)
+remains constant for each quadrature point. So the only thing we need
+keep track of are the positions $x_i=F_K(\hat x_i)$ -- but these
+move with the velocity of the solid domain.
+
+In other words, we don't actually need to keep the solid mesh at all.
+All we need is the positions $x_i(t)$ and corresponding `JxW` values.
Since both of these properties are point-properties (or point-vectors) that are
attached to the solid material, they can be idealized as a set of disconnected
-infinitesimally small particles, which carry the required information with the
-movement of the solid. deal.II already possesses the ability to distribute and
-store such a set of particles in large-scale parallel computations in form of
+infinitesimally small "particles", which carry the required `JxW` information with the
+movement of the solid. deal.II has the ability to distribute and
+store such a set of particles in large-scale parallel computations in the form of
the ParticleHandler class (for details on the implementation see Gassmöller et
al., 2018), and we will make use of this functionality in this tutorial.
Thus, the approach taken in this step is as follows:
-- Create a parallel::distributed::Triangulation for the domain $\Gamma$
-- Create Particles::Particle at the position of the quadrature points on $\Gamma$
+- Create a parallel::distributed::Triangulation for the domain $\Gamma$;
+- Create Particles::Particle at the positions of the quadrature points on $\Gamma$;
- Call the Particles::ParticleHandler::insert_global_particles() function,
to distribute the particles across processors, *following the solid
- triangulation*
-- Attach the quadrature weight as a property to each Particles::Particle object.
+ triangulation*;
+- Attach the `JxW` values as a "property" to each Particles::Particle object.
This structure is relatively expensive to generate, but must only be generated
once per simulation. Once the Particles::ParticleHandler is generated and the