From: Bruno Blais <blais.bruno@gmail.com>
Date: Tue, 4 Jul 2023 16:41:19 +0000 (-0400)
Subject: Update step-68 to use FEPointEvaluation
X-Git-Tag: relicensing~756^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9538b39c873ae274cafcf401079c723c212593f5;p=dealii.git

Update step-68 to use FEPointEvaluation
---

diff --git a/examples/step-68/step-68.cc b/examples/step-68/step-68.cc
index 47680ef395..9e23985f3c 100644
--- a/examples/step-68/step-68.cc
+++ b/examples/step-68/step-68.cc
@@ -536,12 +536,16 @@ namespace Step68
   // In contrast to the previous function in this function we
   // integrate the particle trajectories by interpolating the value of
   // the velocity field at the degrees of freedom to the position of
-  // the particles.
+  // the particles. This is achieved using the FEPointEvaluation object.
   template <int dim>
   void ParticleTracking<dim>::euler_step_interpolated(const double dt)
   {
     Vector<double> local_dof_values(fluid_fe.dofs_per_cell);
 
+    FEPointEvaluation<dim, dim> evaluator(mapping, fluid_fe, update_values);
+    const FEValuesExtractors::Vector velocity(0);
+
+
     // We loop over all the local particles. Although this could be achieved
     // directly by looping over all the cells, this would force us
     // to loop over numerous cells which do not contain particles.
@@ -561,43 +565,36 @@ namespace Step68
 
         // Next, compute the velocity at the particle locations by evaluating
         // the finite element solution at the position of the particles.
-        // This is essentially an optimized version of the particle advection
-        // functionality in step 19, but instead of creating quadrature
-        // objects and FEValues objects for each cell, we do the
-        // evaluation by hand, which is somewhat more efficient and only
-        // matters for this tutorial, because the particle work is the
-        // dominant cost of the whole program.
+        // This is achieved using FEPointEvaluation object.
         const auto pic = particle_handler.particles_in_cell(cell);
         Assert(pic.begin() == particle, ExcInternalError());
+        std::vector<Point<dim>> particle_positions;
         for (auto &p : pic)
-          {
-            const Point<dim> reference_location = p.get_reference_location();
-            Tensor<1, dim>   particle_velocity;
-            for (unsigned int j = 0; j < fluid_fe.dofs_per_cell; ++j)
-              {
-                const auto comp_j = fluid_fe.system_to_component_index(j);
-
-                particle_velocity[comp_j.first] +=
-                  fluid_fe.shape_value(j, reference_location) *
-                  local_dof_values[j];
-              }
-
-            Point<dim> particle_location = particle->get_location();
-            for (int d = 0; d < dim; ++d)
-              particle_location[d] += particle_velocity[d] * dt;
-            p.set_location(particle_location);
-
-            // Again, we store the particle velocity and the processor id in the
-            // particle properties for visualization purposes.
-            ArrayView<double> properties = p.get_properties();
-            for (int d = 0; d < dim; ++d)
-              properties[d] = particle_velocity[d];
-
-            properties[dim] =
-              Utilities::MPI::this_mpi_process(mpi_communicator);
-
-            ++particle;
-          }
+            particle_positions.push_back(p.get_reference_location());
+
+        evaluator.reinit(cell, particle_positions);
+        evaluator.evaluate(make_array_view(local_dof_values),
+                             EvaluationFlags::values);
+
+        // We move the particles using the interpolated velocity field
+        for (unsigned int particle_index = 0;
+            particle != pic.end();
+            ++particle, ++particle_index)
+        {
+          Point<dim> particle_location = particle->get_location();
+          const Tensor<1, dim> &particle_velocity = evaluator.get_value(particle_index) ;                
+          particle_location += particle_velocity*dt;
+          particle->set_location(particle_location);
+
+          // Again, we store the particle velocity and the processor id in the
+          // particle properties for visualization purposes.
+          ArrayView<double> properties = particle->get_properties();
+          for (int d = 0; d < dim; ++d)
+            properties[d] = particle_velocity[d];
+
+          properties[dim] =
+            Utilities::MPI::this_mpi_process(mpi_communicator);
+        }
       }
   }