step-19.h

Go to the documentation of this file.

dim)
 *       , next_unused_particle_id(0)
 *       , n_recently_lost_particles(0)
 *       , n_total_lost_particles(0)
 *       , n_particles_lost_through_anode(0)
 *       , time(0, 1e-4)
 *     {
 *       particle_handler.signals.particle_lost.connect(
 *         [this](const typename Particles::ParticleIterator<dim> &        particle,
 *                const typename Triangulation<dim>::active_cell_iterator &cell) {
 *           this->track_lost_particle(particle, cell);
 *         });
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ThecodeCathodeRaySimulatormake_gridcodefunction"></a> 
 * <h4>The <code>CathodeRaySimulator::make_grid</code> function</h4>
 * 
 
 * 
 * The next function is then responsible for generating the mesh on which
 * we want to solve. Recall how the domain looks like:
 * <p align="center">
 * <img
 * src="https://www.dealii.org/images/steps/developer/step-19.geometry.png"
 * alt="The geometry used in this program"
 * width="600">
 * </p>
 * We subdivide this geometry into a mesh of @f$4\times 2@f$ cells that looks
 * like this:
 * <div class=CodeFragmentInTutorialComment>
 * @code
 *   *---*---*---*---*
 *   \   |   |   |   |
 *    *--*---*---*---*
 *   /   |   |   |   |
 *   *---*---*---*---*
 * @endcode
 * </div>
 * The way this is done is by first defining where the @f$15=5\times 3@f$
 * vertices are located -- here, we say that they are on integer points
 * with the middle one on the left side moved to the right by a value of
 * `delta=0.5`.
 *   
 
 * 
 * In the following, we then have to say which vertices together form
 * the 8 cells. The following code is then entirely equivalent to what
 * we also do in @ref step_14 "step-14":
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::make_grid()
 *     {
 *       static_assert(dim == 2,
 *                     "This function is currently only implemented for 2d.");
 *   
 *       const double       delta = 0.5;
 *       const unsigned int nx    = 5;
 *       const unsigned int ny    = 3;
 *   
 *       const std::vector<Point<dim>> vertices 
 *         = {{0, 0},
 *            {1, 0},
 *            {2, 0},
 *            {3, 0},
 *            {4, 0},
 *            {delta, 1},
 *            {1, 1},
 *            {2, 1},
 *            {3, 1},
 *            {4, 1},
 *            {0, 2},
 *            {1, 2},
 *            {2, 2},
 *            {3, 2},
 *            {4, 2}};
 *       AssertDimension(vertices.size(), nx * ny);
 *   
 *       const std::vector<unsigned int> cell_vertices[(nx - 1) * (ny - 1)] = {
 *         {0, 1, nx + 0, nx + 1},
 *         {1, 2, nx + 1, nx + 2},
 *         {2, 3, nx + 2, nx + 3},
 *         {3, 4, nx + 3, nx + 4},
 *   
 *         {5, nx + 1, 2 * nx + 0, 2 * nx + 1},
 *         {nx + 1, nx + 2, 2 * nx + 1, 2 * nx + 2},
 *         {nx + 2, nx + 3, 2 * nx + 2, 2 * nx + 3},
 *         {nx + 3, nx + 4, 2 * nx + 3, 2 * nx + 4}};
 *   
 * @endcode
 * 
 * With these arrays out of the way, we can move to slightly higher
 * higher-level data structures. We create a vector of CellData
 * objects that store for each cell to be created the vertices in
 * question as well as the @ref GlossMaterialId "material id" (which
 * we will here simply set to zero since we don't use it in the program).
 *     
 
 * 
 * This information is then handed to the
 * Triangulation::create_triangulation() function, and the mesh is twice
 * globally refined.
 * 
 * @code
 *       std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
 *       for (unsigned int i = 0; i < cells.size(); ++i)
 *         {
 *           cells[i].vertices    = cell_vertices[i];
 *           cells[i].material_id = 0;
 *         }
 *   
 *       triangulation.create_triangulation(
 *         vertices,
 *         cells,
 *         SubCellData()); // No boundary information
 *   
 *       triangulation.refine_global(2);
 *   
 * @endcode
 * 
 * The remaining part of the function loops over all cells and their faces,
 * and if a face is at the boundary determines which boundary indicator
 * should be applied to it. The various conditions should make sense if
 * you compare the code with the picture of the geometry above.
 *     
 
 * 
 * Once done with this step, we refine the mesh once more globally.
 * 
 * @code
 *       for (auto &cell : triangulation.active_cell_iterators())
 *         for (auto &face : cell->face_iterators())
 *           if (face->at_boundary())
 *             {
 *               if ((face->center()[0] > 0) && (face->center()[0] < 0.5) &&
 *                   (face->center()[1] > 0) && (face->center()[1] < 2))
 *                 face->set_boundary_id(BoundaryIds::cathode);
 *               else if ((face->center()[0] > 0) && (face->center()[0] < 2))
 *                 face->set_boundary_id(BoundaryIds::focus_element);
 *               else if ((face->center()[0] > 4 - 1e-12) &&
 *                        ((face->center()[1] > 1.5) || (face->center()[1] < 0.5)))
 *                 face->set_boundary_id(BoundaryIds::anode);
 *               else
 *                 face->set_boundary_id(BoundaryIds::open);
 *             }
 *   
 *       triangulation.refine_global(1);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ThecodeCathodeRaySimulatorsetup_systemcodefunction"></a> 
 * <h4>The <code>CathodeRaySimulator::setup_system</code> function</h4>
 * 
 
 * 
 * The next function in this program deals with setting up the various
 * objects related to solving the partial differential equations. It is
 * in essence a copy of the corresponding function in @ref step_6 "step-6" and requires
 * no further discussion.
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::setup_system()
 *     {
 *       dof_handler.distribute_dofs(fe);
 *   
 *       solution.reinit(dof_handler.n_dofs());
 *       system_rhs.reinit(dof_handler.n_dofs());
 *   
 *       constraints.clear();
 *       DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 *   
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                BoundaryIds::cathode,
 *                                                Functions::ConstantFunction<dim>(
 *                                                  -Constants::V0),
 *                                                constraints);
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                BoundaryIds::focus_element,
 *                                                Functions::ConstantFunction<dim>(
 *                                                  -Constants::V0),
 *                                                constraints);
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                BoundaryIds::anode,
 *                                                Functions::ConstantFunction<dim>(
 *                                                  +Constants::V0),
 *                                                constraints);
 *       constraints.close();
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler,
 *                                       dsp,
 *                                       constraints,
 *                                       /*keep_constrained_dofs = */ false);
 *       sparsity_pattern.copy_from(dsp);
 *   
 *       system_matrix.reinit(sparsity_pattern);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ThecodeCathodeRaySimulatorassemble_systemcodefunction"></a> 
 * <h4>The <code>CathodeRaySimulator::assemble_system</code> function</h4>
 * 
 
 * 
 * The function that computes
 * the matrix entries is again in essence a copy of the
 * corresponding function in @ref step_6 "step-6":
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::assemble_system()
 *     {
 *       system_matrix = 0;
 *       system_rhs    = 0;
 *   
 *       const QGauss<dim> quadrature_formula(fe.degree + 1);
 *   
 *       FEValues<dim> fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.dofs_per_cell;
 *   
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           cell_matrix = 0;
 *           cell_rhs    = 0;
 *   
 *           fe_values.reinit(cell);
 *   
 *           for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *             for (const unsigned int i : fe_values.dof_indices())
 *               {
 *                 for (const unsigned int j : fe_values.dof_indices())
 *                   cell_matrix(i, j) +=
 *                     (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
 *                      fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
 *                      fe_values.JxW(q_index));           // dx
 *               }
 *   
 * @endcode
 * 
 * The only interesting part of this function is how it forms the right
 * hand side of the linear system. Recall that the right hand side
 * of the PDE is
 * @f[
 * \sum_p (N e)\delta(\mathbf x-\mathbf x_p),
 * @f]
 * where we have used @f$p@f$ to index the particles here to avoid
 * confusion with the shape function @f$\varphi_i@f$; @f$\mathbf x_p@f$
 * is the position of the @f$p@f$th particle.
 *         
 
 * 
 * When multiplied by a test function @f$\varphi_i@f$ and integrated over
 * the domain results in a right hand side vector
 * @f{align*}{
 * F_i &= \int_\Omega \varphi_i (\mathbf x)\left[
 * \sum_p (N e)\delta(\mathbf x-\mathbf x_p) \right] dx
 * \\  &=  \sum_p (N e) \varphi_i(\mathbf x_p).
 * @f}
 * Note that the final line no longer contains an integral, and
 * consequently also no occurrence of @f$dx@f$ which would require the
 * appearance of the `JxW` symbol in our code.
 *         
 
 * 
 * For a given cell @f$K@f$, this cell's contribution to the right hand
 * side is then
 * @f{align*}{
 * F_i^K &= \sum_{p, \mathbf x_p\in K} (N e) \varphi_i(\mathbf x_p),
 * @f}
 * i.e., we only have to worry about those particles that are actually
 * located on the current cell @f$K@f$.
 *         
 
 * 
 * In practice, what we do here is the following: If there are any
 * particles on the current cell, then we first obtain an iterator range
 * pointing to the first particle of that cell as well as the particle
 * past the last one on this cell (or the end iterator) -- i.e., a
 * half-open range as is common for C++ functions. Knowing now the list
 * of particles, we query their reference locations (with respect to
 * the reference cell), evaluate the shape functions in these reference
 * locations, and compute the force according to the formula above
 * (without any FEValues::JxW).
 *         
 
 * 
 * @note It is worth pointing out that calling the
 * Particles::ParticleHandler::particles_in_cell() and
 * Particles::ParticleHandler::n_particles_in_cell() functions is not
 * very efficient on problems with a large number of particles. But it
 * illustrates the easiest way to write this algorithm, and so we are
 * willing to incur this cost for the moment for expository purposes.
 * We discuss the issue in more detail in the
 * <a href="#extensions">"possibilities for extensions" section</a>
 * below, and use a better approach in @ref step_70 "step-70", for example.
 * 
 * @code
 *           if (particle_handler.n_particles_in_cell(cell) > 0)
 *             for (const auto &particle : particle_handler.particles_in_cell(cell))
 *               {
 *                 const Point<dim> &reference_location =
 *                   particle.get_reference_location();
 *                 for (const unsigned int i : fe_values.dof_indices())
 *                   cell_rhs(i) +=
 *                     (fe.shape_value(i, reference_location) * // phi_i(x_p)
 *                      (-Constants::electrons_per_particle *   // N
 *                       Constants::electron_charge));          // e
 *               }
 *   
 * @endcode
 * 
 * Finally, we can copy the contributions of this cell into
 * the global matrix and right hand side vector:
 * 
 * @code
 *           cell->get_dof_indices(local_dof_indices);
 *           constraints.distribute_local_to_global(
 *             cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
 *         }
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="CathodeRaySimulatorsolve"></a> 
 * <h4>CathodeRaySimulator::solve</h4>
 * 
 
 * 
 * The function that solves the linear system is then again exactly as in
 * @ref step_6 "step-6":
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::solve_field()
 *     {
 *       SolverControl            solver_control(1000, 1e-12);
 *       SolverCG<Vector<double>> solver(solver_control);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       solver.solve(system_matrix, solution, system_rhs, preconditioner);
 *   
 *       constraints.distribute(solution);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="CathodeRaySimulatorrefine_grid"></a> 
 * <h4>CathodeRaySimulator::refine_grid</h4>
 * 
 
 * 
 * The final field-related function is the one that refines the grid. We will
 * call it a number of times in the first time step to obtain a mesh that
 * is well-adapted to the structure of the solution and, in particular,
 * resolves the various singularities in the solution that are due to
 * re-entrant corners and places where the boundary condition type
 * changes. You might want to refer to @ref step_6 "step-6" again for more details:
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::refine_grid()
 *     {
 *       Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *   
 *       KellyErrorEstimator<dim>::estimate(dof_handler,
 *                                          QGauss<dim - 1>(fe.degree + 1),
 *                                          {},
 *                                          solution,
 *                                          estimated_error_per_cell);
 *   
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                       estimated_error_per_cell,
 *                                                       0.1,
 *                                                       0.03);
 *   
 *       triangulation.execute_coarsening_and_refinement();
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="CathodeRaySimulatorcreate_particles"></a> 
 * <h4>CathodeRaySimulator::create_particles</h4>
 * 
 
 * 
 * Let us now turn to the functions that deal with particles. The first one
 * is about the creation of particles. As mentioned in the introduction,
 * we want to create a particle at points of the cathode if the electric
 * field @f$\mathbf E=\nabla V@f$ exceeds a certain threshold, i.e., if
 * @f$|\mathbf E| \ge E_\text{threshold}@f$, and if furthermore the electric field
 * points into the domain (i.e., if @f$\mathbf E \cdot \mathbf n < 0@f$). As is
 * common in the finite element method, we evaluate fields (and their
 * derivatives) at specific evaluation points; typically, these are
 * "quadrature points", and so we create a "quadrature formula" that we will
 * use to designate the points at which we want to evaluate the solution.
 * Here, we will simply take QMidpoint implying that we will only check the
 * threshold condition at the midpoints of faces. We then use this to
 * initialize an object of type FEFaceValues to evaluate the solution at these
 * points.
 *   
 
 * 
 * All of this will then be used in a loop over all cells, their faces, and
 * specifically those faces that are at the boundary and, moreover, the
 * cathode part of the boundary.
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::create_particles()
 *     {
 *       FEFaceValues<dim> fe_face_values(fe,
 *                                        QMidpoint<dim - 1>(),
 *                                        update_quadrature_points |
 *                                          update_gradients |
 *                                          update_normal_vectors);
 *   
 *       std::vector<Tensor<1, dim>> solution_gradients(
 *         fe_face_values.n_quadrature_points);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           if (face->at_boundary() &&
 *               (face->boundary_id() == BoundaryIds::cathode))
 *             {
 *               fe_face_values.reinit(cell, face);
 *   
 * @endcode
 * 
 * So we have found a face on the cathode. Next, we let the
 * FEFaceValues object compute the gradient of the solution at each
 * "quadrature" point, and extract the electric field vector from
 * the gradient in the form of a Tensor variable through the methods
 * discussed in the
 * @ref vector_valued "vector-valued problems" documentation module.
 * 
 * @code
 *               const FEValuesExtractors::Scalar electric_potential(0);
 *               fe_face_values[electric_potential].get_function_gradients(
 *                 solution, solution_gradients);
 *               for (const unsigned int q_point :
 *                    fe_face_values.quadrature_point_indices())
 *                 {
 *                   const Tensor<1, dim> E = solution_gradients[q_point];
 *   
 * @endcode
 * 
 * Electrons can only escape the cathode if the electric field
 * strength exceeds a threshold and,
 * crucially, if the electric field points *into* the domain.
 * Once we have that checked, we create a new
 * Particles::Particle object at this location and insert it
 * into the Particles::ParticleHandler object with a unique ID.
 *                 
 
 * 
 * The only thing that may be not obvious here is that we also
 * associate with this particle the location in the reference
 * coordinates of the cell we are currently on. This is done
 * because we will in downstream functions compute quantities
 * such as the electric field at the location of the particle
 * (e.g., to compute the forces that act on it when updating its
 * position in each time step). Evaluating a finite element
 * field at arbitrary coordinates is quite an expensive
 * operation because shape functions are really only defined on
 * the reference cell, and so when asking for the electric field
 * at an arbitrary point requires us first to determine what the
 * reference coordinates of that point are. To avoid having to
 * do this over and over, we determine these coordinates once
 * and for all and then store these reference coordinates
 * directly with the particle.
 * 
 * @code
 *                   if ((E * fe_face_values.normal_vector(q_point) < 0) &&
 *                       (E.norm() > Constants::E_threshold))
 *                     {
 *                       const Point<dim> &location =
 *                         fe_face_values.quadrature_point(q_point);
 *   
 *                       Particles::Particle<dim> new_particle;
 *                       new_particle.set_location(location);
 *                       new_particle.set_reference_location(
 *                         mapping.transform_real_to_unit_cell(cell, location));
 *                       new_particle.set_id(next_unused_particle_id);
 *                       particle_handler.insert_particle(new_particle, cell);
 *   
 *                       ++next_unused_particle_id;
 *                     }
 *                 }
 *             }
 *   
 * @endcode
 * 
 * At the end of all of these insertions, we let the `particle_handler`
 * update some internal statistics about the particles it stores.
 * 
 * @code
 *       particle_handler.update_cached_numbers();
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="CathodeRaySimulatormove_particles"></a> 
 * <h4>CathodeRaySimulator::move_particles</h4>
 * 
 
 * 
 * The second particle-related function is the one that moves the particles
 * in each time step. To do this, we have to loop over all cells, the
 * particles in each cell, and evaluate the electric field at each of the
 * particles' positions.
 *   
 
 * 
 * The approach used here is conceptually the same used in the
 * `assemble_system()` function: We loop over all cells, find the particles
 * located there (with the same caveat about the inefficiency of the algorithm
 * used here to find these particles), and use FEPointEvaluation object to
 * evaluate the gradient at these positions:
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::move_particles()
 *     {
 *       const double dt = time.get_next_step_size();
 *   
 *       Vector<double>            solution_values(fe.n_dofs_per_cell());
 *       FEPointEvaluation<1, dim> evaluator(mapping, fe, update_gradients);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         if (particle_handler.n_particles_in_cell(cell) > 0)
 *           {
 *             const typename Particles::ParticleHandler<
 *               dim>::particle_iterator_range particles_in_cell =
 *               particle_handler.particles_in_cell(cell);
 *   
 *             std::vector<Point<dim>> particle_positions;
 *             for (const auto &particle : particles_in_cell)
 *               particle_positions.push_back(particle.get_reference_location());
 *   
 *             cell->get_dof_values(solution, solution_values);
 *   
 * @endcode
 * 
 * Then we can ask the FEPointEvaluation object for the gradients of
 * the solution (i.e., the electric field @f$\mathbf E@f$) at these
 * locations and loop over the individual particles:
 * 
 * @code
 *             evaluator.reinit(cell, particle_positions);
 *             evaluator.evaluate(make_array_view(solution_values),
 *                                EvaluationFlags::gradients);
 *   
 *             {
 *               typename Particles::ParticleHandler<dim>::particle_iterator
 *                 particle = particles_in_cell.begin();
 *               for (unsigned int particle_index = 0;
 *                    particle != particles_in_cell.end();
 *                    ++particle, ++particle_index)
 *                 {
 *                   const Tensor<1, dim> &E =
 *                     evaluator.get_gradient(particle_index);
 *   
 * @endcode
 * 
 * Having now obtained the electric field at the location of one
 * of the particles, we use this to update first the velocity
 * and then the position. To do so, let us first get the old
 * velocity out of the properties stored with the particle,
 * compute the acceleration, update the velocity, and store this
 * new velocity again in the properties of the particle. Recall
 * that this corresponds to the first of the following set of
 * update equations discussed in the introduction:
 * @f{align*}{
 * \frac{{\mathbf v}_i^{(n)}
 * -{\mathbf v}_i^{(n-1)}}{\Delta t}
 * &= \frac{e\nabla V^{(n)}}{m}
 * \\ \frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}
 * {\Delta t} &= {\mathbf v}_i^{(n)}.
 * @f}
 * 
 * @code
 *                   const Tensor<1, dim> old_velocity(particle->get_properties());
 *   
 *                   const Tensor<1, dim> acceleration =
 *                     Constants::electron_charge / Constants::electron_mass * E;
 *   
 *                   const Tensor<1, dim> new_velocity =
 *                     old_velocity + acceleration * dt;
 *   
 *                   particle->set_properties(new_velocity);
 *   
 * @endcode
 * 
 * With the new velocity, we can then also update the location
 * of the particle and tell the particle about it.
 * 
 * @code
 *                   const Point<dim> new_location =
 *                     particle->get_location() + dt * new_velocity;
 *                   particle->set_location(new_location);
 *                 }
 *             }
 *           }
 *   
 * @endcode
 * 
 * Having updated the locations and properties (i.e., velocities) of all
 * particles, we need to make sure that the `particle_handler` again knows
 * which cells they are in, and what their locations in the coordinate
 * system of the reference cell are. The following function does that. (It
 * also makes sure that, in parallel computations, particles are moved from
 * one processor to another processor if a particle moves from the subdomain
 * owned by the former to the subdomain owned by the latter.)
 * 
 * @code
 *       particle_handler.sort_particles_into_subdomains_and_cells();
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="CathodeRaySimulatortrack_lost_particle"></a> 
 * <h4>CathodeRaySimulator::track_lost_particle</h4>
 * 
 
 * 
 * The final particle-related function is the one that is called whenever a
 * particle is lost from the simulation. This typically happens if it leaves
 * the domain. If that happens, this function is called both the cell (which
 * we can ask for its new location) and the cell it was previously on. The
 * function then keeps track of updating the number of particles lost in this
 * time step, the total number of lost particles, and then estimates whether
 * the particle left through the hole in the middle of the anode. We do so by
 * first checking whether the cell it was in last had an @f$x@f$ coordinate to the
 * left of the right boundary (located at @f$x=4@f$) and the particle now has a
 * position to the right of the right boundary. If that is so, we compute a
 * direction vector of its motion that is normalized so that the @f$x@f$ component
 * of the direction vector is equal to @f$1@f$. With this direction vector, we can
 * compute where it would have intersected the line @f$x=4@f$. If this intersect
 * is between @f$0.5@f$ and @f$1.5@f$, then we claim that the particle left through
 * the hole and increment a counter.
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::track_lost_particle(
 *       const typename Particles::ParticleIterator<dim> &        particle,
 *       const typename Triangulation<dim>::active_cell_iterator &cell)
 *     {
 *       ++n_recently_lost_particles;
 *       ++n_total_lost_particles;
 *   
 *       const Point<dim> current_location              = particle->get_location();
 *       const Point<dim> approximate_previous_location = cell->center();
 *   
 *       if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
 *         {
 *           const Tensor<1, dim> direction =
 *             (current_location - approximate_previous_location) /
 *             (current_location[0] - approximate_previous_location[0]);
 *   
 *           const double right_boundary_intercept =
 *             approximate_previous_location[1] +
 *             (4 - approximate_previous_location[0]) * direction[1];
 *           if ((right_boundary_intercept > 0.5) &&
 *               (right_boundary_intercept < 1.5))
 *             ++n_particles_lost_through_anode;
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="CathodeRaySimulatorupdate_timestep_size"></a> 
 * <h4>CathodeRaySimulator::update_timestep_size</h4>
 * 
 
 * 
 * As discussed at length in the introduction, we need to respect a time step
 * condition whereby particles can not move further than one cell in one time
 * step. To ensure that this is the case, we again first compute the maximal
 * speed of all particles on each cell, and divide the cell size by that
 * speed. We then compute the next time step size as the minimum of this
 * quantity over all cells, using the safety factor discussed in the
 * introduction, and set this as the desired time step size using the
 * DiscreteTime::set_desired_time_step_size() function.
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::update_timestep_size()
 *     {
 *       if (time.get_step_number() > 0)
 *         {
 *           double min_cell_size_over_velocity = std::numeric_limits<double>::max();
 *   
 *           for (const auto &cell : dof_handler.active_cell_iterators())
 *             if (particle_handler.n_particles_in_cell(cell) > 0)
 *               {
 *                 const double cell_size = cell->minimum_vertex_distance();
 *   
 *                 double max_particle_velocity(0.0);
 *   
 *                 for (const auto &particle :
 *                      particle_handler.particles_in_cell(cell))
 *                   {
 *                     const Tensor<1, dim> velocity(particle.get_properties());
 *                     max_particle_velocity =
 *                       std::max(max_particle_velocity, velocity.norm());
 *                   }
 *   
 *                 if (max_particle_velocity > 0)
 *                   min_cell_size_over_velocity =
 *                     std::min(min_cell_size_over_velocity,
 *                              cell_size / max_particle_velocity);
 *               }
 *   
 *           constexpr double c_safety = 0.5;
 *           time.set_desired_next_step_size(c_safety * 0.5 *
 *                                           min_cell_size_over_velocity);
 *         }
 * @endcode
 * 
 * As mentioned in the introduction, we have to treat the very first
 * time step differently since there, particles are not available yet or
 * do not yet have the information associated that we need for the
 * computation of a reasonable step length. The formulas below follow the
 * discussion in the introduction.
 * 
 * @code
 *       else
 *         {
 *           const QTrapezoid<dim> vertex_quadrature;
 *           FEValues<dim> fe_values(fe, vertex_quadrature, update_gradients);
 *   
 *           std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
 *   
 *           double min_timestep = std::numeric_limits<double>::max();
 *   
 *           for (const auto &cell : dof_handler.active_cell_iterators())
 *             if (particle_handler.n_particles_in_cell(cell) > 0)
 *               {
 *                 const double cell_size = cell->minimum_vertex_distance();
 *   
 *                 fe_values.reinit(cell);
 *                 fe_values.get_function_gradients(solution, field_gradients);
 *   
 *                 double max_E = 0;
 *                 for (const auto q_point : fe_values.quadrature_point_indices())
 *                   max_E = std::max(max_E, field_gradients[q_point].norm());
 *   
 *                 if (max_E > 0)
 *                   min_timestep =
 *                     std::min(min_timestep,
 *                              std::sqrt(0.5 * cell_size *
 *                                        Constants::electron_mass /
 *                                        Constants::electron_charge / max_E));
 *               }
 *   
 *           time.set_desired_next_step_size(min_timestep);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ThecodeCathodeRaySimulatoroutput_resultscodefunction"></a> 
 * <h4>The <code>CathodeRaySimulator::output_results()</code> function</h4>
 * 
 
 * 
 * The final function implementing pieces of the overall algorithm is the one
 * that generates graphical output. In the current context, we want to output
 * both the electric potential field as well as the particle locations and
 * velocities. But we also want to output the electric field, i.e., the
 * gradient of the solution.
 *   
 
 * 
 * deal.II has a general way how one can compute derived quantities from
 * the solution and output those as well. Here, this is the electric
 * field, but it could also be some other quantity -- say, the norm of the
 * electric field, or in fact anything else one could want to compute from
 * the solution @f$V_h(\mathbf x)@f$ or its derivatives. This general solution
 * uses the DataPostprocessor class and, in cases like the one here where we
 * want to output a quantity that represents a vector field, the
 * DataPostprocessorVector class.
 *   
 
 * 
 * Rather than try and explain how this class works, let us simply refer to
 * the documentation of the DataPostprocessorVector class that has essentially
 * this case as a well-documented example.
 * 
 * @code
 *     template <int dim>
 *     class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
 *     {
 *     public:
 *       ElectricFieldPostprocessor()
 *         : DataPostprocessorVector<dim>("electric_field", update_gradients)
 *       {}
 *   
 *       virtual void evaluate_scalar_field(
 *         const DataPostprocessorInputs::Scalar<dim> &input_data,
 *         std::vector<Vector<double>> &computed_quantities) const override
 *       {
 *         AssertDimension(input_data.solution_gradients.size(),
 *                         computed_quantities.size());
 *   
 *         for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
 *           {
 *             AssertDimension(computed_quantities[p].size(), dim);
 *             for (unsigned int d = 0; d < dim; ++d)
 *               computed_quantities[p][d] = input_data.solution_gradients[p][d];
 *           }
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * With this, the `output_results()` function becomes relatively
 * straightforward: We use the DataOut class as we have in almost
 * every one of the previous tutorial programs to output the
 * solution (the "electric potential") and we use the postprocessor
 * defined above to also output its gradient (the "electric
 * field"). This all is then written into a file in VTU format after
 * also associating the current time and time step number with this
 * file, and providing physical units for the two fields to be
 * output.
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::output_results() const
 *     {
 *       {
 *         ElectricFieldPostprocessor<dim> electric_field;
 *         DataOut<dim>                    data_out;
 *         data_out.attach_dof_handler(dof_handler);
 *         data_out.add_data_vector(solution, "electric_potential");
 *         data_out.add_data_vector(solution, electric_field);
 *         data_out.build_patches();
 *   
 *         DataOutBase::VtkFlags output_flags;
 *         output_flags.time  = time.get_current_time();
 *         output_flags.cycle = time.get_step_number();
 *         output_flags.physical_units["electric_potential"] = "V";
 *         output_flags.physical_units["electric_field"]     = "V/m";
 *   
 *         data_out.set_flags(output_flags);
 *   
 *         std::ofstream output("solution-" +
 *                              Utilities::int_to_string(time.get_step_number(), 4) +
 *                              ".vtu");
 *         data_out.write_vtu(output);
 *       }
 *   
 * @endcode
 * 
 * Output the particle positions and properties is not more complicated. The
 * Particles::DataOut class plays the role of the DataOut class for
 * particles, and all we have to do is tell that class where to take
 * particles from and how to interpret the `dim` components of the
 * properties -- namely, as a single vector indicating the velocity, rather
 * than as `dim` scalar properties. The rest is then the same as above:
 * 
 * @code
 *       {
 *         Particles::DataOut<dim> particle_out;
 *         particle_out.build_patches(
 *           particle_handler,
 *           std::vector<std::string>(dim, "velocity"),
 *           std::vector<DataComponentInterpretation::DataComponentInterpretation>(
 *             dim, DataComponentInterpretation::component_is_part_of_vector));
 *   
 *         DataOutBase::VtkFlags output_flags;
 *         output_flags.time                       = time.get_current_time();
 *         output_flags.cycle                      = time.get_step_number();
 *         output_flags.physical_units["velocity"] = "m/s";
 *   
 *         particle_out.set_flags(output_flags);
 *   
 *         std::ofstream output("particles-" +
 *                              Utilities::int_to_string(time.get_step_number(), 4) +
 *                              ".vtu");
 *         particle_out.write_vtu(output);
 *       }
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="CathodeRaySimulatorrun"></a> 
 * <h4>CathodeRaySimulator::run</h4>
 * 
 
 * 
 * The last member function of the principal class of this program is then the
 * driver. At the top, it refines the mesh a number of times by solving the
 * problem (with not particles yet created) on a sequence of finer and finer
 * meshes.
 * 
 * @code
 *     template <int dim>
 *     void CathodeRaySimulator<dim>::run()
 *     {
 *       make_grid();
 *   
 * @endcode
 * 
 * do a few refinement cycles up front
 * 
 * @code
 *       const unsigned int n_pre_refinement_cycles = 3;
 *       for (unsigned int refinement_cycle = 0;
 *            refinement_cycle < n_pre_refinement_cycles;
 *            ++refinement_cycle)
 *         {
 *           setup_system();
 *           assemble_system();
 *           solve_field();
 *           refine_grid();
 *         }
 *   
 *   
 * @endcode
 * 
 * Now do the loop over time. The sequence of steps follows closely the
 * outline of the algorithm discussed in the introduction. As discussed in
 * great detail in the documentation of the DiscreteTime class, while we
 * move the field and particle information forward by one time step, the
 * time stored in the `time` variable is not consistent with where (some of)
 * these quantities are (in the diction of DiscreteTime, this is the "update
 * stage"). The call to `time.advance_time()` makes everything consistent
 * again by setting the `time` variable to the time at which the field and
 * particles already are, and once we are in this "consistent stage", we can
 * generate graphical output and write information about the current state
 * of the simulation to screen.
 * 
 * @code
 *       setup_system();
 *       do
 *         {
 *           std::cout << "Timestep " << time.get_step_number() + 1 << std::endl;
 *           std::cout << "  Field degrees of freedom:                 "
 *                     << dof_handler.n_dofs() << std::endl;
 *   
 *           assemble_system();
 *           solve_field();
 *   
 *           create_particles();
 *           std::cout << "  Total number of particles in simulation:  "
 *                     << particle_handler.n_global_particles() << std::endl;
 *   
 *           n_recently_lost_particles = 0;
 *           update_timestep_size();
 *           move_particles();
 *   
 *           time.advance_time();
 *   
 *           output_results();
 *   
 *           std::cout << "  Number of particles lost this time step:  "
 *                     << n_recently_lost_particles << std::endl;
 *           if (n_total_lost_particles > 0)
 *             std::cout << "  Fraction of particles lost through anode: "
 *                       << 1. * n_particles_lost_through_anode /
 *                            n_total_lost_particles
 *                       << std::endl;
 *   
 *           std::cout << std::endl
 *                     << "  Now at t=" << time.get_current_time()
 *                     << ", dt=" << time.get_previous_step_size() << '.'
 *                     << std::endl
 *                     << std::endl;
 *         }
 *       while (time.is_at_end() == false);
 *     }
 *   } // namespace Step19
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * The final function of the program is then again the `main()` function. It is
 * unchanged in all tutorial programs since @ref step_6 "step-6" and so there is nothing new
 * to discuss:
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         Step19::CathodeRaySimulator<2> cathode_ray_simulator_2d;
 *         cathode_ray_simulator_2d.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;
 *   }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
When this program is run, it produces output that looks as follows:
```
Timestep 1
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  20
  Number of particles lost this time step:  0
 
  Now at t=2.12647e-07, dt=2.12647e-07.
 
Timestep 2
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  40
  Number of particles lost this time step:  0
 
  Now at t=4.14362e-07, dt=2.01715e-07.
 
Timestep 3
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  60
  Number of particles lost this time step:  0
 
  Now at t=5.23066e-07, dt=1.08704e-07.
 
Timestep 4
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  80
  Number of particles lost this time step:  0
 
  Now at t=6.08431e-07, dt=8.53649e-08.
 
 
...
 
 
Timestep 1000
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  394
  Number of particles lost this time step:  2
  Fraction of particles lost through anode: 0.823188
 
  Now at t=3.83822e-05, dt=4.18519e-08.
 
Timestep 1001
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  392
  Number of particles lost this time step:  0
  Fraction of particles lost through anode: 0.823188
 
  Now at t=3.84244e-05, dt=4.22178e-08.
 
 
...
 
 
Timestep 2478
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  364
  Number of particles lost this time step:  0
  Fraction of particles lost through anode: 0.867745
 
  Now at t=9.99509e-05, dt=4.20455e-08.
 
Timestep 2479
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  364
  Number of particles lost this time step:  2
  Fraction of particles lost through anode: 0.867871
 
  Now at t=9.99929e-05, dt=4.20105e-08.
 
Timestep 2480
  Field degrees of freedom:                 4989
  Total number of particles in simulation:  362
  Number of particles lost this time step:  2
  Fraction of particles lost through anode: 0.867047
 
  Now at t=0.0001, dt=7.11401e-09.
```
 
Picking a random few time steps, we can visualize the solution in the
form of streamlines for the electric field and dots for the electrons
(the pictures were created with an earlier version of the program and
are no longer exact representations of the solution as currently
computed -- but they are conceptually still correct):
<div class="twocolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-19.solution.0000.png"
         alt="The solution at time step 0 (t=0 seconds)."
         width="500">
    <br>
    Solution at time step 0 (t=0 seconds).
    <br>
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-19.solution.1400.png"
         alt="The solution at time step 1400 (t=0.000068 seconds)."
         width="500">
    <br>
    Solution at time step 1400 (t=0.000068 seconds).
    <br>
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-19.solution.0700.png"
         alt="The solution at time step 700 (t=0.000035 seconds)."
         width="500">
    <br>
    Solution at time step 700 (t=0.000035 seconds).
    <br>
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step-19.solution.2092.png"
         alt="The solution at time step 2092 (t=0.0001 seconds)."
         width="500">
    <br>
    Solution at time step 2092 (t=0.0001 seconds).
    <br>
  </div>
</div>
 
That said, a more appropriate way to visualize the results of this
program are by creating a video that shows how these electrons move, and how
the electric field changes in response to their motion (again showing
results of an earlier version of the program):
 
@htmlonly
<p align="center">
  <iframe width="560" height="315" src="https://www.youtube.com/embed/HwUtE7xuteE"
   frameborder="0"
   allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
   allowfullscreen></iframe>
 </p>
@endhtmlonly
 
What you can see here is how the "focus element" of the boundary with its negative
voltage repels the electrons and makes sure that they do not just fly away
perpendicular from the cathode (as they do in the initial part of their
trajectories). It also shows how the electric field lines
move around over time, in response to the charges flying by -- in other words,
the feedback the particles have on the electric field that itself drives the
motion of the electrons.
 
The movie suggests that electrons move in "bunches" or "bursts". One element of
this appearance is an artifact of how the movie was created: Every frame of the
movie corresponds to one time step, but the time step length varies. More specifically,
the fastest particle moving through the smallest cell determines the length of the
time step (see the discussion in the introduction), and consequently time steps
are small whenever a (fast) particle moves through the small cells at the right
edge of the domain; time steps are longer again once the particle has left
the domain. This slowing-accelerating effect can easily be visualized by plotting
the time step length shown in the screen output.
 
The second part of this is real, however: The simulation creates a large group
of particles in the beginning, and fewer after about the 300th time step. This
is probably because of the negative charge of the particles in the simulation:
They reduce the magnitude of the electric field at the (also negatively charged
electrode) and consequently reduce the number of points on the cathode at which
the magnitude exceeds the threshold necessary to draw an electron out of the
electrode.
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
 
<a name="Morestatisticsaboutelectrons"></a><h4> More statistics about electrons </h4>
 
 
At the end of the day, we are rarely interested in the *solution* of an equation,
but in numbers that can be *extracted* from it -- in other words, we want to
*postprocess* the solution (see the results section of @ref step_4 "step-4" for an extensive
discussion of this concept). Here, what one would likely be most interested in
assessing are some statistics about the particles.
 
The program already computes the fraction of the electrons that leave the
domain through the hole in the anode. But there are other quantities one might be
interested in. For example, the average velocity of these particles. It would
not be very difficult to obtain each particle's velocity from its properties,
in the same way as we do in the `move_particles()` function, and compute
statistics from it.
 
 
<a name="Abettersynchronizedvisualization"></a><h4> A better-synchronized visualization </h4>
 
 
As discussed above, there is a varying time difference between different frames
of the video because we create output for every time step. A better way to
create movies would be to generate a new output file in fixed time intervals,
regardless of how many time steps lie between each such point.
 
 
<a name="Abettertimeintegrator"></a><h4> A better time integrator </h4>
 
 
The problem we are considering in this program is a coupled, multiphysics
problem. But the way we solve it is by first computing the (electric) potential
field and then updating first the particle velocity and finally the
particle locations. This is what is called an
"operator-splitting method", a concept we will investigate in more detail
in @ref step_58 "step-58".
 
While it is awkward to think of a way to solve this problem that does not involve
splitting the problem into a PDE piece and a particles piece, one *can*
(and probably should!) think of a better way to update the particle
locations. Specifically, the equations we use to update the particle location
are
@f{align*}{
  \frac{{\mathbf v}_i^{(n)}-{\mathbf v}_i^{(n-1)}}{\Delta t} &= \frac{e\nabla V^{(n)}}{m}
  \\
  \frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}{\Delta t} &= {\mathbf v}_i^{(n)}.
@f}
This corresponds to a Lie splitting where we first update one variable
(the velocity) and then, using the already-updated value of the one variable
to update the other (the position). It is well understood that the Lie
splitting incurs a first order error accuracy in the time step (i.e.,
the error introduced by the splitting is @f${\cal O}(\Delta t)@f$) that we know we should
avoid because we can do better. Independently, of course, we incur an error
of the same rate a second time because we use an Euler-like scheme for each of
the two updates (we replace the time derivative by a simple finite difference
quotient containing the old and new value, and we evaluate the right hand side
only at the end time of the interval, at @f$t_n@f$), and a better scheme would also
use a better way to do this kind of update.
 
Better strategies would replace the Lie splitting by something like a
[Strang splitting](https://en.wikipedia.org/wiki/Strang_splitting), and the
update of the particle position and velocity by a scheme such
as the
[leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration)
or more generally
[symplectic integrators](https://en.wikipedia.org/wiki/Symplectic_integrator)
such as the
[Verlet scheme](https://en.wikipedia.org/wiki/Verlet_integration).
 
 
<a name="Parallelization"></a><h4> Parallelization </h4>
 
 
In release mode, the program runs in about 3.5 minutes on one of the author's
laptops at the time of writing this. That's acceptable. But what if we wanted
to make the simulation three-dimensional? If we wanted to not use a maximum
of around 100 particles at any given time (as happens with the parameters
used here) but 100,000? If we needed a substantially finer mesh?
 
In those cases, one would want to run the program not just on a single processor,
but in fact on as many as one has available. This requires parallelization
of both the PDE solution as well as over particles. In practice, while there
are substantial challenges to making this efficient and scale well, these
challenges are all addressed in deal.II itself, and have been demonstrated
on simulations running on 10,000 or more cores. For example, @ref step_40 "step-40" shows
how to parallelize the finite element part, and @ref step_70 "step-70" shows how one can
then also parallelize the particles part.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-19.cc"
*/