solution_evaluator.h

Go to the documentation of this file.

// ------------------------------------------------------
//
// SPDX-License-Identifier: GPL-2.0-or-later
//
// SPDX-FileCopyrightText: Copyright (C) 2024-2025 by the ASPECT authors.
//
// This file is part of ASPECT.
//
// Detailed license information governing the source code and
// contributions can be found in the folder LICENSES and in
// LICENSE and CONTRIBUTING.md at the top level directory.
//
// ------------------------------------------------------

#ifndef _aspect_solution_evaluator_h
#define _aspect_solution_evaluator_h

#include <aspect/global.h>


#include <aspect/simulator_access.h>

#include <deal.II/base/array_view.h>
#include <deal.II/matrix_free/fe_point_evaluation.h>
#include <deal.II/non_matching/mapping_info.h>
#include <deal.II/dofs/dof_handler.h>

namespace aspect
{
  namespace internal
  {
    template <int dim>
    class DynamicFEPointEvaluation
    {
      public:
        DynamicFEPointEvaluation(const unsigned int first_component, const unsigned int n_components)
          : first_component (first_component),
            n_components (n_components)
        {}

        virtual ~DynamicFEPointEvaluation() = default;

        virtual void evaluate(const ArrayView<double> &solution_values,
                              const EvaluationFlags::EvaluationFlags flags) = 0;

        virtual
        small_vector<Tensor<1,dim>>
        get_gradient(const unsigned int evaluation_point) const = 0;

        virtual
        void
        get_gradient(const unsigned int evaluation_point,
                     const ArrayView<Tensor<1,dim>> &gradients) const = 0;

        virtual
        small_vector<double>
        get_value(const unsigned int evaluation_point) const = 0;

        virtual
        void
        get_value(const unsigned int evaluation_point,
                  const ArrayView<double> &solution) const = 0;

        virtual
        unsigned int
        get_first_component() const final
        {
          return first_component;
        }

        virtual
        unsigned int
        get_n_components() const final
        {
          return n_components;
        }

      private:
        unsigned int first_component;

        unsigned int n_components;
    };
  }

  template <int dim>
  class SolutionEvaluator
  {
    public:
      SolutionEvaluator(const SimulatorAccess<dim> &simulator,
                        const UpdateFlags update_flags);

      void
      reinit(const typename DoFHandler<dim>::active_cell_iterator &cell,
             const ArrayView<Point<dim>> &positions);

      void
      evaluate(const ArrayView<double> &solution_values,
               const std::vector<EvaluationFlags::EvaluationFlags> &evaluation_flags);

      void get_solution(const unsigned int evaluation_point,
                        const ArrayView<double> &solution,
                        const std::vector<EvaluationFlags::EvaluationFlags> &evaluation_flags) const;

      void get_gradients(const unsigned int evaluation_point,
                         const ArrayView<Tensor<1, dim>> &gradients,
                         const std::vector<EvaluationFlags::EvaluationFlags> &evaluation_flags) const;

      FEPointEvaluation<dim, dim> &
      get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity);

      NonMatching::MappingInfo<dim> &
      get_mapping_info();

      unsigned int
      n_components() const;

    private:
      NonMatching::MappingInfo<dim> mapping_info;

      FEPointEvaluation<dim, dim> velocity;
      std::unique_ptr<FEPointEvaluation<1, dim>> pressure;
      FEPointEvaluation<1, dim> temperature;

      std::vector<std::unique_ptr<internal::DynamicFEPointEvaluation<dim>>> compositions;

      std::unique_ptr<FEPointEvaluation<dim, dim>> fluid_velocity;
      std::unique_ptr<FEPointEvaluation<1, dim>> compaction_pressure;
      std::unique_ptr<FEPointEvaluation<1, dim>> fluid_pressure;

      std::array<unsigned int, 3> melt_component_indices;

      const SimulatorAccess<dim> &simulator_access;
  };

  template <int dim>
  std::unique_ptr<SolutionEvaluator<dim>>
  construct_solution_evaluator(const SimulatorAccess<dim> &simulator_access,
                               const UpdateFlags update_flags);
}

#endif