distribution.h

Go to the documentation of this file.

/*
  Copyright (C) 2025 - by the authors of the ASPECT code.
 
  This file is part of ASPECT.
 
  ASPECT is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2, or (at your option)
  any later version.
 
  ASPECT is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with ASPECT; see the file LICENSE.  If not see
  <http://www.gnu.org/licenses/>.
*/
 
#ifndef _aspect_particle_distribution_h
#define _aspect_particle_distribution_h
 
#include <aspect/postprocess/interface.h>
#include <aspect/simulator_access.h>
#include <algorithm>
#include <limits>
#include <deal.II/base/table.h>
#include <deal.II/particles/property_pool.h>
#include <deal.II/particles/particle_handler.h>
#include <deal.II/base/function_lib.h>
#include <vector>
 
namespace aspect
{
  namespace Particle
  {
    template <int dim>
    class ParticlePDF
    {
      public:
 
        enum class KernelFunction
        {
          gaussian,
          triangular,
          uniform,
          cutoff_function_w1_dealii,
          cutoff_function_c1_dealii,
          epanechnikov
        };
 
        ParticlePDF(const unsigned int granularity,
                    const double bandwidth,
                    const KernelFunction kernel_function);
 
        ParticlePDF(const double bandwidth,
                    const KernelFunction kernel_function);
 
        void
        fill_from_particle_range(const typename Particles::ParticleHandler<dim>::particle_iterator_range &particle_range,
                                 const std::vector<typename Particles::ParticleHandler<dim>::particle_iterator_range>
                                 &particle_ranges_to_sum_over,
                                 const unsigned int n_particles_in_cell,
                                 const typename ::Mapping<dim> &mapping,
                                 const typename Triangulation<dim>::active_cell_iterator &cell);
 
        void
        insert_kernel_sum_from_particle_range(const Point<dim> &reference_point,
                                              const std::array<unsigned int,dim> &table_index,
                                              const typename Triangulation<dim>::active_cell_iterator &cell,
                                              const std::vector<typename Particles::ParticleHandler<dim>::particle_iterator_range>
                                              &particle_ranges_to_sum_over,
                                              const typename ::Mapping<dim> &mapping);
 
        void
        add_value_to_function_table(const unsigned int x_index,
                                    const unsigned int y_index,
                                    const unsigned int z_index,
                                    const double input_value);
 
        void
        add_value_to_function_table(const std::array<unsigned int,dim> &index_point,
                                    const double input_value);
 
        void
        add_value_to_function_table(const double input_value,
                                    const types::particle_index reference_particle_id);
 
        double
        evaluate_function_at_index(const unsigned int x_index,
                                   const unsigned int y_index,
                                   const unsigned int z_index) const;
 
        void
        compute_statistical_values();
 
        double
        get_max() const;
 
        double
        get_min() const;
 
        Point<dim>
        get_max_position() const;
 
        std::vector<Point<dim>>
        get_min_positions() const;
 
        double
        get_standard_deviation() const;
 
        types::particle_index
        get_max_particle() const;
 
        types::particle_index
        get_min_particle() const;
 
      private:
        Table<dim,double> function_output_table;
 
        double bandwidth;
 
        KernelFunction kernel_function;
 
        unsigned int granularity;
 
        double max;
 
        double min;
 
        Point<dim> max_position;
 
        std::vector<Point<dim>> min_positions;
 
        double standard_deviation;
 
        double mean;
 
        small_vector<double> function_output_vector;
 
        bool is_defined_per_particle;
 
        types::particle_index max_particle_index;
 
        types::particle_index min_particle_index;
 
        double apply_selected_kernel_function(const double distance) const;
 
        double kernelfunction_uniform(const double distance) const;
 
        double kernelfunction_triangular(const double distance) const;
 
        double kernelfunction_gaussian(const double distance) const;
 
        double kernelfunction_epanechnikov(const double distance) const;
    };
  }
}
 
#endif