introspection.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_introspection_h
#define _aspect_introspection_h

#include <deal.II/base/index_set.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/fe.h>

#include <aspect/fe_variable_collection.h>
#include <aspect/parameters.h>

#include <map>

namespace aspect
{
  template <int dim>
  std::vector<VariableDeclaration<dim>>
  construct_default_variables (const Parameters<dim> &parameters);


  struct CompositionalFieldDescription
  {
    enum Type
    {
      chemical_composition = 0,
      stress = 1,
      strain = 2,
      grain_size = 3,
      porosity = 4,
      density = 5,
      entropy = 6,
      generic = 7,
      unspecified = 8
    } type;

    constexpr static unsigned int n_types = 9;

    static
    Type
    parse_type(const std::string &input)
    {
      if (input == "chemical composition")
        return CompositionalFieldDescription::chemical_composition;
      else if (input == "stress")
        return CompositionalFieldDescription::stress;
      else if (input == "strain")
        return CompositionalFieldDescription::strain;
      else if (input == "grain size")
        return CompositionalFieldDescription::grain_size;
      else if (input == "porosity")
        return CompositionalFieldDescription::porosity;
      else if (input == "density")
        return CompositionalFieldDescription::density;
      else if (input == "entropy")
        return CompositionalFieldDescription::entropy;
      else if (input == "generic")
        return CompositionalFieldDescription::generic;
      else if (input == "unspecified")
        return CompositionalFieldDescription::unspecified;
      else
        AssertThrow(false, ExcNotImplemented());

      return CompositionalFieldDescription::Type();
    }
  };

  template <int dim>
  struct Introspection: public FEVariableCollection<dim>
  {
    public:
      Introspection (const std::vector<VariableDeclaration<dim>> &variables,
                     const Parameters<dim> &parameters);

      ~Introspection ();



      const unsigned int n_components;

      const unsigned int n_compositional_fields;

      const bool use_discontinuous_temperature_discretization;

      const std::vector<bool> use_discontinuous_composition_discretization;

      struct ComponentIndices
      {
        std::array<unsigned int, dim> velocities;
        unsigned int                  pressure;
        unsigned int                  temperature;
        std::vector<unsigned int>     compositional_fields;
      };
      const ComponentIndices component_indices;

      const unsigned int n_blocks;

      struct BlockIndices
      {
        unsigned int       velocities;
        unsigned int       pressure;
        unsigned int       temperature;
        std::vector<unsigned int> compositional_fields;

        std::vector<unsigned int> compositional_field_sparsity_pattern;
      };

      const BlockIndices block_indices;

      struct Extractors
      {
        Extractors (const ComponentIndices &component_indices);

        const FEValuesExtractors::Vector              velocities;
        const FEValuesExtractors::Scalar              pressure;
        const FEValuesExtractors::Scalar              temperature;
        const std::vector<FEValuesExtractors::Scalar> compositional_fields;
      };

      const Extractors extractors;

      struct BaseElements
      {
        unsigned int              velocities;
        unsigned int              pressure;
        unsigned int              temperature;
        std::vector<unsigned int> compositional_fields;
      };

      const BaseElements base_elements;

      struct PolynomialDegree
      {
        unsigned int              max_degree;
        unsigned int              velocities;
        unsigned int              temperature;
        std::vector<unsigned int> compositional_fields;
        unsigned int              max_compositional_field;
      };

      const PolynomialDegree polynomial_degree;

      struct Quadratures
      {
        Quadrature<dim>       velocities;
        Quadrature<dim>       pressure;
        Quadrature<dim>       temperature;
        Quadrature<dim>       compositional_field_max;
        std::vector<Quadrature<dim>> compositional_fields;
        Quadrature<dim>       system;
      };

      const Quadratures quadratures;

      struct FaceQuadratures
      {
        Quadrature<dim-1>       velocities;
        Quadrature<dim-1>       pressure;
        Quadrature<dim-1>       temperature;
        Quadrature<dim-1>       compositional_fields;
        Quadrature<dim-1>       system;
      };

      const FaceQuadratures face_quadratures;

      struct ComponentMasks
      {
        ComponentMasks (const FEVariableCollection<dim> &fevs, const Introspection<dim>::ComponentIndices &indices);

        ComponentMask              velocities;

        ComponentMask              pressure;

        ComponentMask              temperature;

        std::vector<ComponentMask> compositional_fields;

        ComponentMask              compositions;
      };

      const ComponentMasks component_masks;

      std::vector<types::global_dof_index> system_dofs_per_block;

      struct IndexSets
      {
        IndexSet system_relevant_set;

        std::vector<IndexSet> system_partitioning;

        std::vector<IndexSet> system_relevant_partitioning;

        std::vector<IndexSet> stokes_partitioning;

        IndexSet locally_owned_pressure_dofs;

        IndexSet locally_owned_melt_pressure_dofs;

        IndexSet locally_owned_fluid_pressure_dofs;
      };

      IndexSets index_sets;

      typename Parameters<dim>::AdvectionFieldMethod::Kind temperature_method;

      std::vector<typename Parameters<dim>::AdvectionFieldMethod::Kind> compositional_field_methods;

      const std::vector<unsigned int> &
      get_composition_base_element_indices() const;

      const std::vector<unsigned int> &
      get_compositional_field_indices_with_base_element(const unsigned int base_element_index) const;

      unsigned int
      compositional_index_for_name (const std::string &name) const;

      std::string
      name_for_compositional_index (const unsigned int index) const;

      const std::vector<std::string> &
      get_composition_names () const;

      const std::vector<CompositionalFieldDescription> &
      get_composition_descriptions () const;

      const std::vector<std::string> &
      chemical_composition_field_names () const;

      const std::vector<unsigned int> &
      chemical_composition_field_indices () const;

      unsigned int
      n_chemical_composition_fields () const;

      bool
      composition_type_exists (const CompositionalFieldDescription::Type &type) const;

      unsigned int
      find_composition_type (const CompositionalFieldDescription::Type &type) const;

      bool
      compositional_name_exists (const std::string &name) const;

      const std::vector<unsigned int> &
      get_indices_for_fields_of_type (const CompositionalFieldDescription::Type &type) const;

      const std::vector<std::string> &
      get_names_for_fields_of_type (const CompositionalFieldDescription::Type &type) const;

      unsigned int
      get_number_of_fields_of_type (const CompositionalFieldDescription::Type &type) const;

      bool
      is_stokes_component (const unsigned int component_index) const;

      bool
      is_composition_component (const unsigned int component_index) const;

    private:
      std::vector<std::string> composition_names;

      std::vector<CompositionalFieldDescription> composition_descriptions;

      std::vector<std::vector<std::string>> composition_names_for_type;

      std::vector<std::vector<unsigned int>> composition_indices_for_type;

      std::vector<unsigned int> composition_base_element_indices;

      std::map<unsigned int, std::vector<unsigned int>> compositional_field_indices_with_base_element;
  };
}


#endif