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