]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Clean up some parts of tutorial 10942/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 21 Sep 2020 13:43:40 +0000 (15:43 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 22 Sep 2020 06:57:07 +0000 (08:57 +0200)
examples/step-19/doc/intro.dox
examples/step-19/step-19.cc

index 072f0fb5d193c17a5277d8908bc06103aa145bd8..faed96ac74a0f1d569076aca5f92c9c89ed8d10b 100644 (file)
@@ -27,7 +27,7 @@ time, then deal.II is clearly not your right tool. On the other hand, if
 this evolution is due to the interaction with the solution of partial differential
 equation, or if having a mesh to determine which particles interact
 with others (such as in the
-[smoothed particle hydrodynamics (SMH)](https://en.wikipedia.org/wiki/Smoothed-particle_hydrodynamics)
+[smoothed particle hydrodynamics (SPH)](https://en.wikipedia.org/wiki/Smoothed-particle_hydrodynamics)
 method), then deal.II has support for you.
 
 The case we will consider here is how electrically charged particles move through
@@ -60,7 +60,7 @@ density. This is augmented by boundary conditions that we will choose as follows
    && \text{on}\; \partial\Omega\setminus\Gamma_\text{cathode}\setminus\Gamma_\text{anode}.
 @f}
 In other words, we prescribe voltages $+V_0$ and $-V_0$ at the two electrodes
-and isolating (Neumann) boundary conditions elsewhere. Since the dynamics of the
+and insulating (Neumann) boundary conditions elsewhere. Since the dynamics of the
 particles are purely due to the electric field $\mathbf E=\nabla V$, we could
 as well have prescribed $2V_0$ and $0$ at the two electrodes -- all that matters
 is the voltage difference at the two electrodes.
@@ -426,7 +426,7 @@ electrons are then accelerated towards the (green) anode which has a
 hole in the middle through which the electrons can escape the device and
 fly on to hit the screen, where they excite the "phosphor" to then emit
 the light that we see from these old-style TV screens. The non-heated
-part of the cathode is not heated, and consequently not subject
+part of the cathode is not subject
 to the emission of electrons -- in the code, we will mark this as the
 "focussing element" of the tube, because its negative electric voltage
 repels the electrons and makes sure that they do not just fly
index 8177bd47e877141ac3047e468984ff1b83b6299f..c88266275001636654ad91dda6880cbaf5c918c6 100644 (file)
@@ -299,7 +299,7 @@ namespace Step19
       }
 
     triangulation.create_triangulation(
-      {std::begin(vertices), std::end(vertices)},
+      vertices,
       cells,
       SubCellData()); // No boundary information
 
@@ -378,9 +378,8 @@ namespace Step19
 
   // @sect4{The <code>CathodeRaySimulator::assemble_system</code> function}
 
-  // The same is true for the function that assembles the linear system to be
-  // solved in each time step. At least that is true for the computation
-  // of the matrix entries, which is again in essence a copy of the
+  // The function that computes
+  // the matrix entries is again in essence a copy of the
   // corresponding function in step-6:
   template <int dim>
   void CathodeRaySimulator<dim>::assemble_system()
@@ -452,11 +451,11 @@ namespace Step19
         // 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
-        // number of particles, we can start to collect their reference
-        // locations (with respect to the reference cell), which we have stored
-        // with each particle when they were created or when it was last moved
-        // (see below).
+        // 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
@@ -468,52 +467,16 @@ namespace Step19
         //   <a href="#extensions">"possibilities for extensions" section</a>
         //   below, and use a better approach in step-70, for example.
         if (particle_handler.n_particles_in_cell(cell) > 0)
-          {
-            std::vector<Point<dim>> particle_reference_locations;
-
-            const typename Particles::ParticleHandler<
-              dim>::particle_iterator_range particles_in_cell =
-              particle_handler.particles_in_cell(cell);
-
-            const unsigned int n_particles_in_cell =
-              particle_handler.n_particles_in_cell(cell);
-
-            particle_reference_locations.resize(n_particles_in_cell);
-
+          for (const auto &particle : particle_handler.particles_in_cell(cell))
             {
-              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)
-                particle_reference_locations[particle_index] =
-                  particle->get_reference_location();
-            }
-
-            // Now that we know where the particles on the current cell are
-            // located with regard to the reference cell's coordinate system, we
-            // can create a Quadrature object with these locations and then an
-            // FEValues object that we will use to evaluate the shape functions
-            // at these locations. The contribution to the right hand side
-            // vector then immediately follows from the formula shown above.
-            // Note again the absence of the call to FEValues::JxW that would
-            // have to be present if we were evaluating an integral.
-            const Quadrature<dim> quadrature_formula(
-              particle_reference_locations);
-            FEValues<dim> fe_value(mapping,
-                                   fe,
-                                   quadrature_formula,
-                                   update_values);
-
-            for (const unsigned int q_index :
-                 fe_values.quadrature_point_indices())
+              const Point<dim> reference_location =
+                particle.get_reference_location();
               for (const unsigned int i : fe_values.dof_indices())
                 cell_rhs(i) +=
-                  (fe_values.shape_value(i, q_index) *   // phi_i(x_p)
-                   (-Constants::electrons_per_particle * // N
-                    Constants::electron_charge));        // e
-          }
-
+                  (fe.shape_value(i, reference_location) * // phi_i(x_p)
+                   (-Constants::electrons_per_particle *   // N
+                    Constants::electron_charge));          // e
+            }
 
         // Finally, we can copy the contributions of this cell into
         // the global matrix and right hand side vector:
@@ -687,8 +650,6 @@ namespace Step19
   {
     const double dt = time.get_next_step_size();
 
-    std::vector<Point<dim>>     particle_positions;
-    std::vector<Tensor<1, dim>> field_gradients;
 
     for (const auto &cell : dof_handler.active_cell_iterators())
       if (particle_handler.n_particles_in_cell(cell) > 0)
@@ -697,33 +658,25 @@ namespace Step19
             dim>::particle_iterator_range particles_in_cell =
             particle_handler.particles_in_cell(cell);
 
-          const unsigned int n_particles_in_cell =
-            particle_handler.n_particles_in_cell(cell);
-
-          particle_positions.resize(n_particles_in_cell);
-          {
-            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)
-              particle_positions[particle_index] =
-                particle->get_reference_location();
-          }
+          std::vector<Point<dim>> particle_positions;
+          for (const auto &particle : particles_in_cell)
+            particle_positions.push_back(particle.get_reference_location());
 
           const Quadrature<dim> quadrature_formula(particle_positions);
-          FEValues<dim>         fe_value(mapping,
-                                 fe,
-                                 quadrature_formula,
-                                 update_gradients);
+          FEValues<dim>         particle_position_fe_values(mapping,
+                                                    fe,
+                                                    quadrature_formula,
+                                                    update_gradients);
 
-          fe_value.reinit(cell);
+          particle_position_fe_values.reinit(cell);
 
           // Then we can ask the FEValues object for the gradients of the
           // solution (i.e., the electric field $\mathbf E$) at these locations
           // and loop over the individual particles:
-          field_gradients.resize(n_particles_in_cell);
-          fe_value.get_function_gradients(solution, field_gradients);
+          std::vector<Tensor<1, dim>> field_gradients(
+            quadrature_formula.size());
+          particle_position_fe_values.get_function_gradients(solution,
+                                                             field_gradients);
 
           {
             typename Particles::ParticleHandler<dim>::particle_iterator
@@ -760,9 +713,9 @@ namespace Step19
                 particle->set_properties(make_array_view(new_velocity));
 
                 // With the new velocity, we can then also update the location
-                // of the particle and set tell the particle about it.
+                // of the particle and tell the particle about it.
                 const Point<dim> new_location =
-                  particle->get_location() + dt * old_velocity;
+                  particle->get_location() + dt * new_velocity;
                 particle->set_location(new_location);
               }
           }

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.