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>

namespace aspect
{
  namespace Particle
  {
    template <int dim>
    class ParticlePDF
    {
      public:

        enum class KernelFunction
        {
          gaussian,
          triangular,
          uniform,
          cutoff_function_w1_dealii
        };

        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 unsigned int n_particles_in_cell);

        void
        insert_kernel_sum_from_particle_range(const Point<dim> reference_point,
                                              std::array<unsigned int,dim> table_index,
                                              const unsigned int n_particles_in_cell,
                                              const typename Particles::ParticleHandler<dim>::particle_iterator_range particle_range);

        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;

        double
        get_standard_deviation() const;

      private:
        Table<dim,double> function_output_table;

        double bandwidth;

        KernelFunction kernel_function;

        unsigned int granularity;

        double max;

        double min;

        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;
    };
  }
}

#endif