From: Bruno Blais 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 void ParticleTracking::euler_step_interpolated(const double dt) { Vector local_dof_values(fluid_fe.dofs_per_cell); + FEPointEvaluation 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> particle_positions; for (auto &p : pic) - { - const Point 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 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 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 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 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); + } } }