solution_evaluator.h

Go to the documentation of this file.

/*
 Copyright (C) 2024 - 2024 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_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
{
  using namespace dealii;

  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