interface.h

Go to the documentation of this file.

/*
  Copyright (C) 2011 - 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_material_model_interface_h
#define _aspect_material_model_interface_h
 
#include <aspect/global.h>
#include <aspect/plugins.h>
#include <aspect/material_model/utilities.h>
 
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/numerics/data_postprocessor.h>
#include <deal.II/base/signaling_nan.h>
 
namespace aspect
{
  template <int dim>
  struct Introspection;
 
 
  namespace MaterialModel
  {
    namespace NonlinearDependence
    {
      enum Dependence
      {
        uninitialized        = 0,
 
        none                 = 1,
        temperature          = 2,
        pressure             = 4,
        strain_rate          = 8,
        compositional_fields = 16,
 
        any_variable         = temperature | pressure | strain_rate | compositional_fields
      };
 
 
      inline Dependence operator | (const Dependence d1,
                                    const Dependence d2)
      {
        return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
      }
 
      inline Dependence operator |= (Dependence &d1,
                                     const Dependence d2)
      {
        d1 = (d1 | d2);
        return d1;
      }
 
      struct ModelDependence
      {
        Dependence viscosity;
 
        Dependence density;
 
        Dependence compressibility;
 
        Dependence specific_heat;
 
        Dependence thermal_conductivity;
 
        ModelDependence ();
      };
 
      bool
      identifies_single_variable(const Dependence dependence);
    }
 
    namespace MaterialProperties
    {
      enum Property
      {
        uninitialized                  = 0,
 
        none                           = 1,
        viscosity                      = 2,
        density                        = 4,
        thermal_expansion_coefficient  = 8,
        specific_heat                  = 16,
        thermal_conductivity           = 32,
        compressibility                = 64,
        entropy_derivative_pressure    = 128,
        entropy_derivative_temperature = 256,
        reaction_terms                 = 512,
        reaction_rates                 = 1024,
        additional_outputs             = 2048,
 
        equation_of_state_properties   = density |
                                         thermal_expansion_coefficient |
                                         specific_heat |
                                         compressibility |
                                         entropy_derivative_pressure |
                                         entropy_derivative_temperature,
        all_properties                 = equation_of_state_properties |
                                         viscosity |
                                         thermal_conductivity |
                                         reaction_terms |
                                         reaction_rates |
                                         additional_outputs
      };
 
      inline Property operator | (const Property d1,
                                  const Property d2)
      {
        return Property(static_cast<int>(d1) | static_cast<int>(d2));
      }
 
      inline Property operator |= (Property &d1,
                                   const Property d2)
      {
        return (d1 | d2);
      }
    }
 
 
    // Forward declaration:
    template <int dim>
    class AdditionalMaterialInputs;
 
    template <int dim>
    class MaterialModelInputs
    {
      public:
        MaterialModelInputs(const unsigned int n_points,
                            const unsigned int n_comp);
 
        MaterialModelInputs(const DataPostprocessorInputs::Vector<dim> &input_data,
                            const Introspection<dim> &introspection,
                            const bool compute_strain_rate = true);
 
 
        MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
                            const typename DoFHandler<dim>::active_cell_iterator &cell,
                            const Introspection<dim> &introspection,
                            const LinearAlgebra::BlockVector &solution_vector,
                            const bool compute_strain_rate = true);
 
        MaterialModelInputs (const MaterialModelInputs &source);
 
        MaterialModelInputs (MaterialModelInputs &&)  noexcept = default;
 
        MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
 
        MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
 
        void reinit(const FEValuesBase<dim,dim>                          &fe_values,
                    const typename DoFHandler<dim>::active_cell_iterator &cell,
                    const Introspection<dim>                             &introspection,
                    const LinearAlgebra::BlockVector                     &solution_vector,
                    const bool                                            compute_strain_rate = true);
 
        unsigned int n_evaluation_points() const;
 
        bool requests_property(const MaterialProperties::Property &property) const;
 
        std::vector<Point<dim>> position;
 
        std::vector<double> temperature;
 
        std::vector<double> pressure;
 
        std::vector<Tensor<1,dim>> pressure_gradient;
 
        std::vector<Tensor<1,dim>> velocity;
 
        std::vector<std::vector<double>> composition;
 
        std::vector<SymmetricTensor<2,dim>> strain_rate;
 
        typename DoFHandler<dim>::active_cell_iterator current_cell;
 
        MaterialProperties::Property requested_properties;
 
        template <class AdditionalInputType>
        std::shared_ptr<AdditionalInputType>
        get_additional_input_object();
 
        template <class AdditionalInputType>
        std::shared_ptr<const AdditionalInputType>
        get_additional_input_object() const;
 
        template <class AdditionalInputType>
        DEAL_II_DEPRECATED
        AdditionalInputType *
        get_additional_input();
 
        template <class AdditionalInputType>
        DEAL_II_DEPRECATED
        const AdditionalInputType *
        get_additional_input() const;
 
        template <class AdditionalInputType>
        bool
        has_additional_input_object() const;
 
        std::vector<std::shared_ptr<AdditionalMaterialInputs<dim>>> additional_inputs;
    };
 
 
    // Forward declaration:
    template <int dim>
    class AdditionalMaterialOutputs;
 
 
    template <int dim>
    class MaterialModelOutputs
    {
      public:
        MaterialModelOutputs (const unsigned int n_points,
                              const unsigned int n_comp);
 
        MaterialModelOutputs (const MaterialModelOutputs &source);
 
        MaterialModelOutputs (MaterialModelOutputs &&)  noexcept = default;
 
        MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
 
        MaterialModelOutputs &operator= (MaterialModelOutputs &&)  noexcept = default;
 
        unsigned int n_evaluation_points() const;
 
        std::vector<double> viscosities;
 
        std::vector<double> densities;
 
        std::vector<double> thermal_expansion_coefficients;
 
        std::vector<double> specific_heat;
 
        std::vector<double> thermal_conductivities;
 
        std::vector<double> compressibilities;
 
        std::vector<double> entropy_derivative_pressure;
 
        std::vector<double> entropy_derivative_temperature;
 
        std::vector<std::vector<double>> reaction_terms;
 
        std::vector<std::shared_ptr<AdditionalMaterialOutputs<dim>>> additional_outputs;
 
        template <class AdditionalOutputType>
        std::shared_ptr<AdditionalOutputType>
        get_additional_output_object();
 
        template <class AdditionalOutputType>
        std::shared_ptr<const AdditionalOutputType>
        get_additional_output_object() const;
 
        template <class AdditionalOutputType>
        DEAL_II_DEPRECATED
        AdditionalOutputType *
        get_additional_output();
 
        template <class AdditionalOutputType>
        DEAL_II_DEPRECATED
        const AdditionalOutputType *
        get_additional_output() const;
 
        template <class AdditionalInputType>
        bool
        has_additional_output_object() const;
 
        void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
    };
 
 
 
    namespace MaterialAveraging
    {
      enum AveragingOperation
      {
        none,
        arithmetic_average,
        harmonic_average,
        geometric_average,
        pick_largest,
        project_to_Q1,
        log_average,
        harmonic_average_only_viscosity,
        geometric_average_only_viscosity,
        project_to_Q1_only_viscosity,
        default_averaging
      };
 
 
      std::string get_averaging_operation_names ();
 
      AveragingOperation parse_averaging_operation_name (const std::string &s);
 
      template <int dim>
      void average (const AveragingOperation operation,
                    const typename DoFHandler<dim>::active_cell_iterator &cell,
                    const Quadrature<dim>         &quadrature_formula,
                    const Mapping<dim>            &mapping,
                    const MaterialProperties::Property &requested_properties,
                    MaterialModelOutputs<dim>     &values_out);
 
      void average_property (const AveragingOperation  operation,
                             const FullMatrix<double>      &projection_matrix,
                             const FullMatrix<double>      &expansion_matrix,
                             std::vector<double>           &values_out);
 
      AveragingOperation get_averaging_operation_for_viscosity(const AveragingOperation operation);
    }
 
    template <int dim>
    class AdditionalMaterialInputs
    {
      public:
        virtual ~AdditionalMaterialInputs() = default;
 
        virtual void
        fill (const LinearAlgebra::BlockVector &solution,
              const FEValuesBase<dim>          &fe_values,
              const Introspection<dim>         &introspection) = 0;
    };
 
 
    template <int dim>
    class AdditionalMaterialOutputs
    {
      public:
        virtual ~AdditionalMaterialOutputs() = default;
 
        virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
                              const FullMatrix<double>  &/*projection_matrix*/,
                              const FullMatrix<double>  &/*expansion_matrix*/)
        {}
    };
 
 
    template <int dim>
    class NamedAdditionalMaterialOutputs : public AdditionalMaterialOutputs<dim>
    {
      public:
        NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
 
        NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
                                       const unsigned int n_points);
 
        ~NamedAdditionalMaterialOutputs() override;
 
        const std::vector<std::string> &get_names() const;
 
        virtual std::vector<double> get_nth_output(const unsigned int idx) const;
 
        void average (const MaterialAveraging::AveragingOperation /*operation*/,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {}
 
 
        std::vector<std::vector<double>> output_values;
 
      private:
        const std::vector<std::string> names;
    };
 
 
    template <int dim>
    class SeismicAdditionalOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        SeismicAdditionalOutputs(const unsigned int n_points);
 
        std::vector<double> get_nth_output(const unsigned int idx) const override;
 
        std::vector<double> vs;
 
        std::vector<double> vp;
    };
 
 
    template <int dim>
    class ReactionRateOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        ReactionRateOutputs (const unsigned int n_points,
                             const unsigned int n_comp);
 
        std::vector<double> get_nth_output(const unsigned int idx) const override;
 
        std::vector<std::vector<double>> reaction_rates;
    };
 
 
 
    template <int dim>
    class PhaseOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        PhaseOutputs(const unsigned int n_points);
    };
 
 
 
    template <int dim>
    class PrescribedFieldOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        PrescribedFieldOutputs (const unsigned int n_points,
                                const unsigned int n_comp);
 
        std::vector<double> get_nth_output(const unsigned int idx) const override;
 
        std::vector<std::vector<double>> prescribed_field_outputs;
    };
 
    template <int dim>
    class PrescribedTemperatureOutputs : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        PrescribedTemperatureOutputs (const unsigned int n_points);
 
        std::vector<double> get_nth_output(const unsigned int idx) const override;
 
        std::vector<double> prescribed_temperature_outputs;
    };
 
    template <int dim>
    class AdditionalMaterialOutputsStokesRHS: public AdditionalMaterialOutputs<dim>
    {
      public:
        AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
          : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
        {}
 
        ~AdditionalMaterialOutputsStokesRHS() override
          = default;
 
        void average (const MaterialAveraging::AveragingOperation /*operation*/,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {
          // TODO: not implemented
        }
 
        std::vector<Tensor<1,dim>> rhs_u;
 
        std::vector<double> rhs_p;
 
        std::vector<double> rhs_melt_pc;
    };
 
 
 
    template <int dim>
    class PrescribedPlasticDilation : public NamedAdditionalMaterialOutputs<dim>
    {
      public:
        explicit PrescribedPlasticDilation (const unsigned int n_points);
 
        std::vector<double> get_nth_output(const unsigned int idx) const override;
 
        std::vector<double> dilation_lhs_term;
 
        std::vector<double> dilation_rhs_term;
    };
 
 
 
    template <int dim>
    class ElasticOutputs: public AdditionalMaterialOutputs<dim>
    {
      public:
        ElasticOutputs(const unsigned int n_points)
          : elastic_force(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>()),
            viscoelastic_strain_rate(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
        {}
 
        ~ElasticOutputs() override
          = default;
 
        void average (const MaterialAveraging::AveragingOperation operation,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {
          AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
          return;
        }
 
        std::vector<SymmetricTensor<2,dim>> elastic_force;
 
        std::vector<SymmetricTensor<2,dim>> viscoelastic_strain_rate;
    };
 
 
 
    template <int dim>
    class EnthalpyOutputs : public AdditionalMaterialOutputs<dim>
    {
      public:
        EnthalpyOutputs(const unsigned int n_points)
          : enthalpies_of_fusion(n_points, numbers::signaling_nan<double>())
        {}
 
        virtual ~EnthalpyOutputs()
          = default;
 
        void average (const MaterialAveraging::AveragingOperation operation,
                      const FullMatrix<double>  &/*projection_matrix*/,
                      const FullMatrix<double>  &/*expansion_matrix*/) override
        {
          AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
          return;
        }
 
        std::vector<double> enthalpies_of_fusion;
    };
 
 
 
    template <int dim>
    class Interface : public Plugins::InterfaceBase
    {
      public:
        using MaterialModelInputs = MaterialModel::MaterialModelInputs<dim>;
 
        using MaterialModelOutputs = MaterialModel::MaterialModelOutputs<dim>;
 
        const NonlinearDependence::ModelDependence &
        get_model_dependence () const;
 
        virtual bool is_compressible () const = 0;
        virtual
        void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
                       MaterialModel::MaterialModelOutputs<dim> &out) const = 0;
 
        virtual
        void
        create_additional_named_outputs (MaterialModelOutputs &outputs) const;
 
 
        virtual
        void
        fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
                                              const LinearAlgebra::BlockVector        &solution,
                                              const FEValuesBase<dim>                 &fe_values,
                                              const Introspection<dim>                &introspection) const;
 
      protected:
        NonlinearDependence::ModelDependence model_dependence;
    };
 
    template <int dim>
    void
    register_material_model (const std::string &name,
                             const std::string &description,
                             void (*declare_parameters_function) (ParameterHandler &),
                             std::unique_ptr<Interface<dim>> (*factory_function) ());
 
    template <int dim>
    std::unique_ptr<Interface<dim>>
    create_material_model (const std::string &model_name);
 
 
    template <int dim>
    std::unique_ptr<Interface<dim>>
    create_material_model (ParameterHandler &prm);
 
 
    template <int dim>
    std::string
    get_valid_model_names_pattern ();
 
 
    template <int dim>
    void
    declare_parameters (ParameterHandler &prm);
 
 
 
    template <int dim>
    void
    write_plugin_graph (std::ostream &output_stream);
 
 
 
// --------------------- template function definitions ----------------------------------
 
    template <int dim>
    template <class AdditionalInputType>
    std::shared_ptr<AdditionalInputType>
    MaterialModelInputs<dim>::get_additional_input_object()
    {
      for (unsigned int i=0; i<additional_inputs.size(); ++i)
        if (dynamic_cast<AdditionalInputType *> (additional_inputs[i].get()))
          return std::dynamic_pointer_cast<AdditionalInputType>(additional_inputs[i]);
 
      return nullptr;
    }
 
 
    template <int dim>
    template <class AdditionalInputType>
    std::shared_ptr<const AdditionalInputType>
    MaterialModelInputs<dim>::get_additional_input_object() const
    {
      for (unsigned int i=0; i<additional_inputs.size(); ++i)
        if (dynamic_cast<AdditionalInputType *> (additional_inputs[i].get()))
          return std::dynamic_pointer_cast<const AdditionalInputType>(additional_inputs[i]);
 
      return nullptr;
    }
 
 
    template <int dim>
    template <class AdditionalOutputType>
    std::shared_ptr<AdditionalOutputType>
    MaterialModelOutputs<dim>::get_additional_output_object()
    {
      for (unsigned int i=0; i<additional_outputs.size(); ++i)
        if (dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get()))
          return std::dynamic_pointer_cast<AdditionalOutputType>(additional_outputs[i]);
 
      return nullptr;
    }
 
 
    template <int dim>
    template <class AdditionalOutputType>
    std::shared_ptr<const AdditionalOutputType>
    MaterialModelOutputs<dim>::get_additional_output_object() const
    {
      for (unsigned int i=0; i<additional_outputs.size(); ++i)
        if (dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get()))
          return std::dynamic_pointer_cast<const AdditionalOutputType>(additional_outputs[i]);
 
      return nullptr;
    }
 
 
    // The following four functions are deprecated:
    template <int dim>
    template <class AdditionalInputType>
    AdditionalInputType *
    MaterialModelInputs<dim>::get_additional_input()
    {
      return get_additional_input_object<AdditionalInputType>().get();
    }
 
 
    template <int dim>
    template <class AdditionalInputType>
    const AdditionalInputType *
    MaterialModelInputs<dim>::get_additional_input() const
    {
      return get_additional_input_object<AdditionalInputType>().get();
    }
 
 
    template <int dim>
    template <class AdditionalOutputType>
    AdditionalOutputType *
    MaterialModelOutputs<dim>::get_additional_output()
    {
      return get_additional_output_object<AdditionalOutputType>().get();
    }
 
 
    template <int dim>
    template <class AdditionalOutputType>
    const AdditionalOutputType *
    MaterialModelOutputs<dim>::get_additional_output() const
    {
      return get_additional_output_object<AdditionalOutputType>().get();
    }
 
 
    template <int dim>
    template <class AdditionalInputType>
    bool
    MaterialModelInputs<dim>::has_additional_input_object() const
    {
      for (unsigned int i=0; i<additional_inputs.size(); ++i)
        if (dynamic_cast<AdditionalInputType *> (additional_inputs[i].get()))
          return true;
 
      return false;
    }
 
 
    template <int dim>
    template <class AdditionalOutputType>
    bool
    MaterialModelOutputs<dim>::has_additional_output_object() const
    {
      for (unsigned int i=0; i<additional_outputs.size(); ++i)
        if (dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get()))
          return true;
 
      return false;
    }
 
 
 
    template <int dim>
    void MaterialModelOutputs<dim>::move_additional_outputs_from(MaterialModelOutputs<dim> &other)
    {
      Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
      this->additional_outputs = std::move(other.additional_outputs);
    }
 
 
#define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
  template class classname<2>; \
  template class classname<3>; \
  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
  { \
    aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2>> \
    dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
                                name, description); \
    aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3>> \
    dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
                                name, description); \
  }
  }
}
 
 
#endif