user_constants.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
 
#include <boost/algorithm/string/predicate.hpp>
#include <boost/range/algorithm_ext/erase.hpp>
#include <boost/variant.hpp>
 
#include <prismspf/core/conditional_ostreams.h>
#include <prismspf/core/exceptions.h>
#include <prismspf/core/type_enums.h>
 
#include <prismspf/config.h>
 
#include <fstream>
#include <map>
 
PRISMS_PF_BEGIN_NAMESPACE

template <unsigned int dim>
class UserConstants
{
public:
  using InputVariant = boost::variant<double,
                                      int,
                                      bool,
                                      std::string,
                                      dealii::Tensor<1, dim>,
                                      dealii::Tensor<2, dim>,
                                      dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>>;

  InputVariant
  construct_user_constant(std::vector<std::string> &model_constants_strings);

  [[nodiscard]] double
  get_double(const std::string &constant_name) const;

  [[nodiscard]] int
  get_int(const std::string &constant_name) const;

  [[nodiscard]] bool
  get_bool(const std::string &constant_name) const;

  [[nodiscard]] std::string
  get_string(const std::string &constant_name) const;

  [[nodiscard]] dealii::Tensor<1, dim>
  get_rank_1_tensor(const std::string &constant_name) const;

  [[nodiscard]] dealii::Tensor<2, dim>
  get_rank_2_tensor(const std::string &constant_name) const;

  [[nodiscard]] dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>
  get_elasticity_tensor(const std::string &constant_name) const;

  void
  add_user_constant(const std::string        &constant_name,
                    std::vector<std::string> &model_constants_strings)
  {
    model_constants[constant_name] = construct_user_constant(model_constants_strings);
  };

  void
  print() const;
 
  static void
  strip_spaces(std::string &line);
 
  static bool
  check_keyword_match(const std::string &line, const std::string &keyword);

  static bool
  parse_line(std::string        line,
             const std::string &keyword,
             const std::string &entry_name,
             std::string       &out_string,
             bool               expect_equals_sign);

  static std::set<std::string>
  get_names(const std::string &input_file_name);

  void
  declare_parameters(dealii::ParameterHandler &parameter_handler) const;

  void
  assign_parameters(dealii::ParameterHandler &parameter_handler);

  std::string file_name;
 
private:
  unsigned int
  compute_tensor_parentheses(const unsigned int             &n_elements,
                             const std::vector<std::string> &tensor_elements);

  void
  remove_parentheses(std::vector<std::string> &tensor_elements);

  dealii::Tensor<1, dim>
  compute_rank_1_tensor_constant(const unsigned int             &n_elements,
                                 const std::vector<std::string> &tensor_elements);

  dealii::Tensor<2, dim>
  compute_rank_2_tensor_constant(const unsigned int             &n_elements,
                                 const std::vector<std::string> &tensor_elements);

  InputVariant
  primitive_model_constant(std::vector<std::string> &model_constants_strings);
 
  [[nodiscard]] dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>
  get_cij_tensor(std::vector<double> elastic_constants,
                 const std::string  &elastic_const_symmetry,
                 const StressState  &stress_state) const;
 
  [[nodiscard]] dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>
  get_cij_matrix(const ElasticityModel     &model,
                 const std::vector<double> &constants,
                 const StressState         &stress_state) const;

  std::map<std::string, InputVariant> model_constants;

  class VariantPrinter : public boost::static_visitor<>
  {
  public:
    void
    operator()(double value) const
    {
      ConditionalOStreams::pout_summary() << value;
    }
 
    void
    operator()(int value) const
    {
      ConditionalOStreams::pout_summary() << value;
    }
 
    void
    operator()(bool value) const
    {
      ConditionalOStreams::pout_summary() << std::boolalpha << value;
    }
 
    void
    operator()(const std::string &value) const
    {
      ConditionalOStreams::pout_summary() << value;
    }
 
    void
    operator()(const dealii::Tensor<1, dim> &value) const
    {
      ConditionalOStreams::pout_summary() << "Tensor<1, " << dim << ">: ";
      for (unsigned int i = 0; i < dim; ++i)
        {
          ConditionalOStreams::pout_summary() << value[i] << ' ';
        }
    }
 
    void
    operator()(const dealii::Tensor<2, dim> &value) const
    {
      ConditionalOStreams::pout_summary() << "Tensor<2, " << dim << ">: ";
      for (unsigned int i = 0; i < dim; ++i)
        {
          for (unsigned int j = 0; j < dim; ++j)
            {
              ConditionalOStreams::pout_summary() << value[i][j] << ' ';
            }
        }
    }
 
    template <unsigned int D = dim>
    void
    operator()(const dealii::Tensor<2, (2 * D) - 1 + (D / 3)> &value) const
    requires((D != ((2 * D) - 1 + (D / 3))))
    {
      constexpr unsigned int dimension = (2 * D) - 1 + (D / 3);
      ConditionalOStreams::pout_summary() << "Tensor<2, " << dimension << ">: ";
      for (unsigned int i = 0; i < dimension; ++i)
        {
          for (unsigned int j = 0; j < dimension; ++j)
            {
              ConditionalOStreams::pout_summary() << value[i][j] << ' ';
            }
        }
    }
  };
};
 
template <unsigned int dim>
inline double
UserConstants<dim>::get_double(const std::string &constant_name) const
{
  Assert(model_constants.find(constant_name) != model_constants.end(),
         dealii::ExcMessage(
           "Mismatch between constants in parameters.prm and CustomPDE.h. The constant "
           "that you attempted to access was " +
           constant_name + "."));
 
  return boost::get<double>(model_constants.at(constant_name));
}
 
template <unsigned int dim>
inline int
UserConstants<dim>::get_int(const std::string &constant_name) const
{
  Assert(model_constants.find(constant_name) != model_constants.end(),
         dealii::ExcMessage(
           "Mismatch between constants in parameters.prm and CustomPDE.h. The constant "
           "that you attempted to access was " +
           constant_name + "."));
 
  return boost::get<int>(model_constants.at(constant_name));
}
 
template <unsigned int dim>
inline bool
UserConstants<dim>::get_bool(const std::string &constant_name) const
{
  Assert(model_constants.find(constant_name) != model_constants.end(),
         dealii::ExcMessage(
           "Mismatch between constants in parameters.prm and CustomPDE.h. The constant "
           "that you attempted to access was " +
           constant_name + "."));
 
  return boost::get<bool>(model_constants.at(constant_name));
}
 
template <unsigned int dim>
inline std::string
UserConstants<dim>::get_string(const std::string &constant_name) const
{
  Assert(model_constants.find(constant_name) != model_constants.end(),
         dealii::ExcMessage(
           "Mismatch between constants in parameters.prm and CustomPDE.h. The constant "
           "that you attempted to access was " +
           constant_name + "."));
 
  return boost::get<std::string>(model_constants.at(constant_name));
}
 
template <unsigned int dim>
inline dealii::Tensor<1, dim>
UserConstants<dim>::get_rank_1_tensor(const std::string &constant_name) const
{
  Assert(model_constants.find(constant_name) != model_constants.end(),
         dealii::ExcMessage(
           " Mismatch between constants in parameters.prm and "
           "CustomPDE.h. The constant that you attempted to access was " +
           constant_name + "."));
 
  return boost::get<dealii::Tensor<1, dim>>(model_constants.at(constant_name));
}
 
template <unsigned int dim>
inline dealii::Tensor<2, dim>
UserConstants<dim>::get_rank_2_tensor(const std::string &constant_name) const
{
  Assert(model_constants.find(constant_name) != model_constants.end(),
         dealii::ExcMessage(
           "Mismatch between constants in parameters.prm and CustomPDE.h. The constant "
           "that you attempted to access was " +
           constant_name + "."));
 
  return boost::get<dealii::Tensor<2, dim>>(model_constants.at(constant_name));
}
 
template <unsigned int dim>
inline dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>
UserConstants<dim>::get_elasticity_tensor(const std::string &constant_name) const
{
  Assert(model_constants.find(constant_name) != model_constants.end(),
         dealii::ExcMessage(
           "Mismatch between constants in parameters.prm and CustomPDE.h. The constant "
           "that you attempted to access was " +
           constant_name + "."));
 
  return boost::get<dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>>(
    model_constants.at(constant_name));
}
 
template <unsigned int dim>
inline unsigned int
UserConstants<dim>::compute_tensor_parentheses(
  const unsigned int             &n_elements,
  const std::vector<std::string> &tensor_elements)
{
  unsigned int open_parentheses  = 0;
  unsigned int close_parentheses = 0;
 
  for (unsigned int element = 0; element < n_elements; element++)
    {
      for (const char character : tensor_elements.at(element))
        {
          if (character == '(')
            {
              ++open_parentheses;
            }
          else if (character == ')')
            {
              ++close_parentheses;
            }
        }
    }
 
  if (open_parentheses != close_parentheses)
    {
      AssertThrow(false,
                  dealii::ExcMessage("User-defined elastic constant list does not have "
                                     "the same number of open and close parentheses."));
    }
 
  return open_parentheses;
}
 
template <unsigned int dim>
inline void
UserConstants<dim>::remove_parentheses(std::vector<std::string> &tensor_elements)
{
  for (std::string &element : tensor_elements)
    {
      boost::range::remove_erase(element, '(');
      boost::range::remove_erase(element, ')');
    }
}
 
template <unsigned int dim>
inline dealii::Tensor<1, dim>
UserConstants<dim>::compute_rank_1_tensor_constant(
  const unsigned int             &n_elements,
  const std::vector<std::string> &tensor_elements)
{
  AssertThrow(n_elements == 3,
              dealii::ExcMessage("The columns in user-defined constant tensors must be "
                                 "equal to the maximum number of dimensions."));
 
  dealii::Tensor<1, dim> temp;
  for (unsigned int i = 0; i < dim; i++)
    {
      temp[i] = dealii::Utilities::string_to_double(tensor_elements.at(i));
    }
 
  return temp;
}
 
template <unsigned int dim>
inline dealii::Tensor<2, dim>
UserConstants<dim>::compute_rank_2_tensor_constant(
  const unsigned int             &n_elements,
  const std::vector<std::string> &tensor_elements)
{
  AssertThrow(n_elements == 9,
              dealii::ExcMessage("User-defined constant tensor does not have the "
                                 "correct number of elements, matrices must be 3x3."));
 
  const unsigned int row_length = 3;
 
  dealii::Tensor<2, dim> temp;
  for (unsigned int i = 0; i < dim; i++)
    {
      for (unsigned int j = 0; j < dim; j++)
        {
          temp[i][j] =
            dealii::Utilities::string_to_double(tensor_elements.at((i * row_length) + j));
        }
    }
 
  return temp;
}
 
template <unsigned int dim>
inline typename UserConstants<dim>::InputVariant
UserConstants<dim>::construct_user_constant(
  std::vector<std::string> &model_constants_strings)
{
  // Ensure that the input includes a value and a type
  AssertThrow(
    model_constants_strings.size() > 1,
    dealii::ExcMessage(
      "At least two fields are required for user-defined variables (value and type)."));
 
  std::vector<std::string> model_constants_type_strings =
    dealii::Utilities::split_string_list(model_constants_strings.at(
                                           model_constants_strings.size() - 1),
                                         ' ');
 
  if (model_constants_strings.size() == 2)
    {
      return primitive_model_constant(model_constants_strings);
    }
 
  if (boost::iequals(model_constants_type_strings.at(0), "tensor"))
    {
      const unsigned int n_elements = model_constants_strings.size() - 1;
 
      const unsigned int open_parentheses =
        compute_tensor_parentheses(n_elements, model_constants_strings);
 
      AssertThrow(open_parentheses <= 4,
                  FeatureNotImplemented("3rd rank tensors and above"));
 
      remove_parentheses(model_constants_strings);
 
      // Rank 1 tensor
      if (open_parentheses == 1)
        {
          return compute_rank_1_tensor_constant(n_elements, model_constants_strings);
        }
      // Rank 2 tensor
      return compute_rank_2_tensor_constant(n_elements, model_constants_strings);
    }
  if (boost::iequals(model_constants_type_strings.at(1), "elastic") &&
      boost::iequals(model_constants_type_strings.at(2), "constants"))
    {
      const unsigned int n_elements = model_constants_strings.size() - 1;
 
      remove_parentheses(model_constants_strings);
 
      // Load in the elastic constants as a vector
      std::vector<double> temp_elastic_constants(n_elements);
      for (unsigned int i = 0; i < n_elements; i++)
        {
          temp_elastic_constants[i] =
            dealii::Utilities::string_to_double(model_constants_strings.at(i));
        }
      const std::string &elastic_const_symmetry = model_constants_type_strings.at(0);
 
      // get the stress state for the 2D case
      StressState stress_state = StressState::ThreeDimension;
      if (dim == 2)
        {
          if (model_constants_type_strings.size() == 4)
            {
              if (boost::iequals(model_constants_type_strings.at(3), "planestress"))
                {
                  stress_state = StressState::PlaneStress;
                }
              else if (boost::iequals(model_constants_type_strings.at(3), "planestrain"))
                {
                  stress_state = StressState::PlaneStrain;
                }
              else
                {
                  AssertThrow(
                    false,
                    dealii::ExcMessage(
                      "Wrong 2D stress state. Use planestress or planestrain."));
                }
            }
          else if (model_constants_type_strings.size() == 3)
            {
              // Default to plane strain if nothing is specified
              stress_state = StressState::PlaneStrain;
            }
          else
            {
              AssertThrow(false,
                          dealii::ExcMessage(
                            "Wrong format of user-defined elastic constants."));
            }
        }
 
      dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)> temp =
        get_cij_tensor(temp_elastic_constants, elastic_const_symmetry, stress_state);
      return temp;
    }
 
  AssertThrow(false,
              dealii::ExcMessage(
                "Only user-defined constant tensors may have multiple elements."));
  return 0;
}
 
template <unsigned int dim>
inline typename UserConstants<dim>::InputVariant
UserConstants<dim>::primitive_model_constant(
  std::vector<std::string> &model_constants_strings)
{
  std::vector<std::string> model_constants_type_strings =
    dealii::Utilities::split_string_list(model_constants_strings.at(
                                           model_constants_strings.size() - 1),
                                         ' ');
 
  if (boost::iequals(model_constants_type_strings.at(0), "double"))
    {
      return dealii::Utilities::string_to_double(model_constants_strings.at(0));
    }
  if (boost::iequals(model_constants_type_strings.at(0), "int"))
    {
      return dealii::Utilities::string_to_int(model_constants_strings.at(0));
    }
  if (boost::iequals(model_constants_type_strings.at(0), "bool"))
    {
      return boost::iequals(model_constants_strings.at(0), "true");
    }
  if (boost::iequals(model_constants_type_strings.at(0), "string"))
    {
      return model_constants_strings.at(0);
    }
 
  AssertThrow(false,
              dealii::ExcMessage(
                "The type for user-defined variables must be `double`, `int`, "
                "`bool`, `string`, `tensor`, or `elastic constants`."));
  return 0;
}
 
template <unsigned int dim>
inline dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>
UserConstants<dim>::get_cij_tensor(std::vector<double> elastic_constants,
                                   const std::string  &elastic_const_symmetry,
                                   const StressState  &stress_state) const
{
  // First set the material model
  ElasticityModel mat_model = Isotropic;
  if (elastic_const_symmetry == "isotropic")
    {
      mat_model = ElasticityModel::Isotropic;
    }
  else if (elastic_const_symmetry == "transverse")
    {
      mat_model = ElasticityModel::Transverse;
    }
  else if (elastic_const_symmetry == "orthotropic")
    {
      mat_model = ElasticityModel::Orthotropic;
    }
  else if (elastic_const_symmetry == "anisotropic")
    {
      mat_model = ElasticityModel::Anisotropic;
    }
  else
    {
      AssertThrow(false, dealii::ExcMessage("Invalid elasticity tensor type"));
    }
 
  // If the material model is anisotropic for a 2D calculation but the elastic
  // constants are given for a 3D calculation, change the elastic constant
  // vector to the 2D form
  constexpr unsigned int max_number = 21;
  if ((mat_model == Anisotropic) && (dim == 2) && elastic_constants.size() == max_number)
    {
      std::vector<double> elastic_constants_temp = elastic_constants;
      elastic_constants.clear();
      const std::vector<unsigned int> indices_2d = {0, 1, 5, 6, 10, 14};
      std::transform(indices_2d.begin(),
                     indices_2d.end(),
                     std::back_inserter(elastic_constants),
                     [&elastic_constants_temp](unsigned int index)
                     {
                       return elastic_constants_temp.at(index);
                     });
    }
 
  return get_cij_matrix(mat_model, elastic_constants, stress_state);
}
 
template <unsigned int dim>
inline dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>
UserConstants<dim>::get_cij_matrix(const ElasticityModel     &model,
                                   const std::vector<double> &constants,
                                   const StressState         &stress_state) const
{
  // Initialize stiffness tensor
  dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)> stiffness;
 
  switch (dim)
    {
      case 1:
        {
          const int xx_dir = 0;
 
          switch (model)
            {
              // For the 1D case, it make little sense to accept anything besides an
              // isotropic elasticity tensor. One hiccup, is that if a user is debugging
              // an application and switches to 1D, they will have to modify the
              // elasticity constants to align with that. While more burdensome, there's
              // less of a chance of producing spurious behavior.
              case Isotropic:
                {
                  const double modulus = constants.at(0);
 
                  stiffness[xx_dir][xx_dir] = modulus;
                  break;
                }
              default:
                AssertThrow(false,
                            dealii::ExcMessage(
                              "Invalid elasticity model type for 1D. We only accept "
                              "isotropic elasticity tensors."));
            }
          break;
        }
      case 2:
        {
          // The voigt indexing scheme for 2 dimensions
          const int xx_dir = 0;
          const int yy_dir = 1;
          const int xy_dir = 2;
 
          switch (model)
            {
              // Like the 1D case, it is nonsensical to have transverse or orthotropic
              // stiffness tensors, so we throw an error.
              case Isotropic:
                {
                  // For isotropic stiffness tensors, we can simplify the computation to
                  // two parameters: $\lambda$ and $\mu$, where $\mu$ is the shear
                  // modulus. In cartesian coordinates,
                  // $$$
                  // c_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik}
                  // \delta_{jl} + \delta_{il} \delta_{kj})
                  // $$$
                  // In 2D, we distinguish between plane stress and plane strain
                  const double modulus = constants.at(0);
                  const double poisson = constants.at(1);
 
                  const double shear_modulus = modulus / (2 * (1 + poisson));
                  double       lambda        = 0.0;
                  switch (stress_state)
                    {
                      case PlaneStress:
                        {
                          lambda = poisson * modulus / ((1 - poisson * poisson));
                          break;
                        }
                      case PlaneStrain:
                        {
                          lambda =
                            poisson * modulus / ((1 + poisson) * (1 - 2 * poisson));
                          break;
                        }
                      default:
                        AssertThrow(false,
                                    dealii::ExcMessage("Invalid stress state type"));
                    }
 
                  stiffness[xx_dir][xx_dir] = stiffness[yy_dir][yy_dir] =
                    lambda + 2 * shear_modulus;
                  stiffness[xy_dir][xy_dir] = shear_modulus;
                  stiffness[xx_dir][yy_dir] = stiffness[yy_dir][xx_dir] = lambda;
                  break;
                }
              case Orthotropic:
                {
                  // In 2D, we distinguish between plane stress and plane strain
                  switch (stress_state)
                    {
                      case PlaneStress:
                        {
                          const double E1   = constants.at(0);
                          const double E2   = constants.at(1);
                          const double nu12 = constants.at(2);
                          const double G12  = constants.at(3);
 
                          const double nu21  = nu12 * (E2 / E1);
                          const double denom = 1.0 - (nu12 * nu21);
 
                          stiffness[xx_dir][xx_dir] = E1 / denom;
                          stiffness[yy_dir][yy_dir] = E2 / denom;
                          stiffness[xx_dir][yy_dir] = stiffness[yy_dir][xx_dir] =
                            (nu12 * E2) / denom;
                          stiffness[xy_dir][xy_dir] = G12;
                          break;
                        }
                      case PlaneStrain:
                        {
                          const double E1   = constants.at(0);
                          const double E2   = constants.at(1);
                          const double E3   = constants.at(2);
                          const double nu12 = constants.at(3);
                          const double nu13 = constants.at(4);
                          const double nu23 = constants.at(5);
                          const double G12  = constants.at(6);
 
                          const double nu21 = nu12 * (E2 / E1);
                          const double nu31 = nu13 * (E3 / E1);
                          const double nu32 = nu23 * (E3 / E2);
 
                          const double delta = 1.0 - (nu12 * nu21) - (nu23 * nu32) -
                                               (nu13 * nu31) - (2.0 * nu12 * nu23 * nu31);
 
                          stiffness[xx_dir][xx_dir] =
                            (E1 * (1.0 - (nu23 * nu32))) / delta;
                          stiffness[yy_dir][yy_dir] =
                            (E2 * (1.0 - (nu13 * nu31))) / delta;
                          stiffness[xx_dir][yy_dir] = stiffness[yy_dir][xx_dir] =
                            (E1 * (nu21 + (nu31 * nu23))) / delta;
                          stiffness[xy_dir][xy_dir] = G12;
                          break;
                        }
                      default:
                        AssertThrow(false,
                                    dealii::ExcMessage("Invalid stress state type"));
                    }
                  break;
                }
              case Anisotropic:
                {
                  // In the anisotropic case, every entry is specified (given the symmetry
                  // constraints). Also, ignore magic numbers because it is simpler to
                  // hardcode this.
 
                  // NOLINTBEGIN(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
 
                  stiffness[xx_dir][xx_dir] = constants.at(0);
                  stiffness[yy_dir][yy_dir] = constants.at(1);
                  stiffness[xy_dir][xy_dir] = constants.at(2);
                  stiffness[xx_dir][yy_dir] = stiffness[yy_dir][xx_dir] = constants.at(3);
                  stiffness[xx_dir][xy_dir] = stiffness[xy_dir][xx_dir] = constants.at(4);
                  stiffness[yy_dir][xy_dir] = stiffness[xy_dir][yy_dir] = constants.at(5);
 
                  // NOLINTEND(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
                  break;
                }
              default:
                AssertThrow(false, dealii::ExcMessage("Invalid elasticity model type"));
            }
          break;
        }
      case 3:
        {
          const int xx_dir = 0;
          const int yy_dir = 1;
          const int zz_dir = 2;
          const int yz_dir = 3;
          const int xz_dir = 4;
          const int xy_dir = 5;
 
          switch (model)
            {
              case Isotropic:
                {
                  // For isotropic stiffness tensors, we can simplify the computation to
                  // two parameters: $\lambda$ and $\mu$, where $\mu$ is the shear
                  // modulus. In cartesian coordinates,
                  // $$$
                  // c_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik}
                  // \delta_{jl} + \delta_{il} \delta_{kj})
                  // $$$
                  const double modulus = constants.at(0);
                  const double poisson = constants.at(1);
 
                  const double shear_modulus = modulus / (2 * (1 + poisson));
                  const double lambda =
                    poisson * modulus / ((1 + poisson) * (1 - 2 * poisson));
 
                  stiffness[xx_dir][xx_dir]     = stiffness[yy_dir][yy_dir] =
                    stiffness[zz_dir][zz_dir]   = lambda + 2 * shear_modulus;
                  stiffness[yz_dir][yz_dir]     = stiffness[xz_dir][xz_dir] =
                    stiffness[xy_dir][xy_dir]   = shear_modulus;
                  stiffness[xx_dir][yy_dir]     = stiffness[yy_dir][xx_dir] =
                    stiffness[xx_dir][zz_dir]   = stiffness[zz_dir][xx_dir] =
                      stiffness[yy_dir][zz_dir] = stiffness[zz_dir][yy_dir] = lambda;
                  break;
                }
              case Anisotropic:
                {
                  // In the anisotropic case, every entry is specified (given the symmetry
                  // constraints). Also, ignore magic numbers because it is simpler to
                  // hardcode this.
 
                  // NOLINTBEGIN(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
 
                  stiffness[xx_dir][xx_dir] = constants[0];
                  stiffness[yy_dir][yy_dir] = constants[1];
                  stiffness[zz_dir][zz_dir] = constants[2];
                  stiffness[yz_dir][yz_dir] = constants[3];
                  stiffness[xz_dir][xz_dir] = constants[4];
                  stiffness[xy_dir][xy_dir] = constants[5];
                  stiffness[xx_dir][yy_dir] = stiffness[yy_dir][xx_dir] = constants[6];
                  stiffness[xx_dir][zz_dir] = stiffness[zz_dir][xx_dir] = constants[7];
                  stiffness[xx_dir][yz_dir] = stiffness[yz_dir][xx_dir] = constants[8];
                  stiffness[xx_dir][xz_dir] = stiffness[xz_dir][xx_dir] = constants[9];
                  stiffness[xx_dir][xy_dir] = stiffness[xy_dir][xx_dir] = constants[10];
                  stiffness[yy_dir][zz_dir] = stiffness[zz_dir][yy_dir] = constants[11];
                  stiffness[yy_dir][yz_dir] = stiffness[yz_dir][yy_dir] = constants[12];
                  stiffness[yy_dir][xz_dir] = stiffness[xz_dir][yy_dir] = constants[13];
                  stiffness[yy_dir][xy_dir] = stiffness[xy_dir][yy_dir] = constants[14];
                  stiffness[zz_dir][yz_dir] = stiffness[yz_dir][zz_dir] = constants[15];
                  stiffness[zz_dir][xz_dir] = stiffness[xz_dir][zz_dir] = constants[16];
                  stiffness[zz_dir][xy_dir] = stiffness[xy_dir][zz_dir] = constants[17];
                  stiffness[yz_dir][xz_dir] = stiffness[xz_dir][yz_dir] = constants[18];
                  stiffness[yz_dir][xy_dir] = stiffness[xy_dir][yz_dir] = constants[19];
                  stiffness[xz_dir][xy_dir] = stiffness[xy_dir][xz_dir] = constants[20];
 
                  // NOLINTEND(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
                  break;
                }
              case Transverse:
                {
                  // TODO (landinjm): implement
                  AssertThrow(false, FeatureNotImplemented("Transverse material"));
                  break;
                }
              case Orthotropic:
                {
                  const double E1   = constants.at(0);
                  const double E2   = constants.at(1);
                  const double E3   = constants.at(2);
                  const double nu12 = constants.at(3);
                  const double nu13 = constants.at(4);
                  const double nu23 = constants.at(5);
                  const double G12  = constants.at(6);
                  const double G13  = constants.at(7);
                  const double G23  = constants.at(8);
 
                  const double nu21 = nu12 * (E2 / E1);
                  const double nu31 = nu13 * (E3 / E1);
                  const double nu32 = nu23 * (E3 / E2);
 
                  const double delta = 1.0 - (nu12 * nu21) - (nu23 * nu32) -
                                       (nu13 * nu31) - (2.0 * nu12 * nu23 * nu31);
 
                  stiffness[xx_dir][xx_dir] = (E1 * (1.0 - nu23 * nu32)) / delta;
                  stiffness[yy_dir][yy_dir] = (E2 * (1.0 - nu13 * nu31)) / delta;
                  stiffness[zz_dir][zz_dir] = (E3 * (1.0 - nu12 * nu21)) / delta;
                  stiffness[xx_dir][yy_dir] = stiffness[yy_dir][xx_dir] =
                    (E1 * (nu21 + nu31 * nu23)) / delta;
                  stiffness[xx_dir][zz_dir] = stiffness[zz_dir][xx_dir] =
                    (E1 * (nu31 + nu21 * nu32)) / delta;
                  stiffness[yy_dir][zz_dir] = stiffness[zz_dir][yy_dir] =
                    (E2 * (nu32 + nu12 * nu31)) / delta;
                  stiffness[yz_dir][yz_dir] = G23;
                  stiffness[xz_dir][xz_dir] = G13;
                  stiffness[xy_dir][xy_dir] = G12;
                  break;
                }
              default:
                AssertThrow(false, dealii::ExcMessage("Invalid elasticity model type"));
            }
          break;
        }
      default:
        {
          Assert(false, UnreachableCode());
        }
    }
 
  return stiffness;
}
 
template <unsigned int dim>
void
UserConstants<dim>::print() const
{
  if (!model_constants.empty())
    {
      ConditionalOStreams::pout_summary()
        << "================================================\n"
        << "  User Constants\n"
        << "================================================\n";
 
      for (const auto &[constant_name, variant] : model_constants)
        {
          ConditionalOStreams::pout_summary() << constant_name << ": ";
          boost::apply_visitor(VariantPrinter(), variant);
          ConditionalOStreams::pout_summary() << "\n";
        }
      ConditionalOStreams::pout_summary() << "\n" << std::flush;
    }
}
 
template <unsigned int dim>
void
UserConstants<dim>::strip_spaces(std::string &line)
{
  while ((!line.empty()) && (line[0] == ' ' || line[0] == '\t'))
    {
      line.erase(0, 1);
    }
  while ((!line.empty()) &&
         (line[line.size() - 1] == ' ' || line[line.size() - 1] == '\t'))
    {
      line.erase(line.size() - 1, std::string::npos);
    }
}
 
template <unsigned int dim>
bool
UserConstants<dim>::check_keyword_match(const std::string &line,
                                        const std::string &keyword)
{
  // Early return if the line is less than the keyword size
  if (line.size() < keyword.size())
    {
      return false;
    }
 
  // Check that the line begins with the keyword
  for (unsigned int i = 0; i < keyword.size(); i++)
    {
      if (line[i] != keyword[i])
        {
          return false;
        }
    }
 
  return true;
}
 
template <unsigned int dim>
bool
UserConstants<dim>::parse_line(std::string        line,
                               const std::string &keyword,
                               const std::string &entry_name,
                               std::string       &out_string,
                               bool               expect_equals_sign)
{
  // Remove spaces from the front and back
  strip_spaces(line);
 
  // Check whether the line starts with 'keyword'. If not, try next line (if the entry is
  // "", then zero spaces after the keyword is ok)
  if (!check_keyword_match(line, keyword))
    {
      return false;
    }
 
  if (!entry_name.empty())
    {
      if (line[keyword.size()] != ' ' && line[keyword.size()] != '\t')
        {
          return false;
        }
    }
 
  // Delete the "keyword" and any more spaces, if present
  line.erase(0, keyword.size());
  strip_spaces(line);
 
  // now see whether the next word is the word we look for
  if (!line.starts_with(entry_name))
    {
      return false;
    }
 
  line.erase(0, entry_name.size());
  strip_spaces(line);
 
  // we'd expect an equals size here if expect_equals_sign is true
  if (expect_equals_sign)
    {
      if ((line.empty()) || (line[0] != '='))
        {
          return false;
        }
    }
 
  // remove comment
  const std::string::size_type pos = line.find('#');
  if (pos != std::string::npos)
    {
      line.erase(pos);
    }
 
  // trim the equals sign at the beginning and possibly following spaces
  // as well as spaces at the end
  if (expect_equals_sign)
    {
      line.erase(0, 1);
    }
  strip_spaces(line);
 
  out_string = line;
  return true;
}
 
template <unsigned int dim>
std::set<std::string>
UserConstants<dim>::get_names(const std::string &input_file_name)
{
  const std::string keyword             = "set";
  const std::string entry_name_begining = "Model constant";
  std::ifstream     input_file;
  input_file.open(input_file_name);
 
  std::string line;
  std::string entry;
 
  std::set<std::string> entry_name_end_list;
 
  // Loop through each line
  while (std::getline(input_file, line))
    {
      if (parse_line(line, keyword, entry_name_begining, entry, false))
        {
          // Strip whitespace, the equals sign, and everything after the equals
          // sign
 
          // Strip whitespace at the beginning
          while ((!entry.empty()) && (entry[0] == ' ' || entry[0] == '\t'))
            {
              entry.erase(0, 1);
            }
 
          // Strip everything up to the equals sign
          while ((!entry.empty()) && (entry[entry.size() - 1] != '='))
            {
              entry.erase(entry.size() - 1, std::string::npos);
            }
 
          // Strip the equals sign
          entry.erase(entry.size() - 1, std::string::npos);
 
          // Strip whitespace between the entry name and the equals sign
          while ((!entry.empty()) &&
                 (entry[entry.size() - 1] == ' ' || entry[entry.size() - 1] == '\t'))
            {
              entry.erase(entry.size() - 1, std::string::npos);
            }
 
          // Add it to the list
          AssertThrow(entry_name_end_list.insert(entry).second,
                      dealii::ExcMessage(
                        "Non-unique constant name in parameters.prm. The constant that "
                        "you attempted to create was \"" +
                        entry + "\"."));
        }
    }
  return entry_name_end_list;
}
 
template <unsigned int dim>
void
UserConstants<dim>::declare_parameters(dealii::ParameterHandler &parameter_handler) const
{
  if (file_name.empty())
    {
      return;
    }
  const std::set<std::string> model_constant_names = get_names(file_name);
  for (const std::string &constant_name : model_constant_names)
    {
      std::string constants_text = "Model constant ";
      constants_text.append(constant_name);
      parameter_handler.declare_entry(constants_text,
                                      "0",
                                      dealii::Patterns::Anything(),
                                      "The value of a user-defined constant.");
    }
}
 
template <unsigned int dim>
void
UserConstants<dim>::assign_parameters(dealii::ParameterHandler &parameter_handler)
{
  if (file_name.empty())
    {
      return;
    }
  const std::set<std::string> model_constant_names = get_names(file_name);
  for (const std::string &constant_name : model_constant_names)
    {
      std::string constants_text = "Model constant ";
      constants_text.append(constant_name);
 
      std::vector<std::string> model_constants_strings =
        dealii::Utilities::split_string_list(parameter_handler.get(constants_text));
      add_user_constant(constant_name, model_constants_strings);
    }
}
 
PRISMS_PF_END_NAMESPACE