parsed_function.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2007 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#include <deal.II/base/parsed_function.h>
#include <deal.II/base/utilities.h>

#include <cstdio>

DEAL_II_NAMESPACE_OPEN

namespace Functions
{
  template <int dim>
  ParsedFunction<dim>::ParsedFunction(const unsigned int n_components,
                                      const double       h)
    : AutoDerivativeFunction<dim>(h, n_components)
    , function_object(n_components)
  {}



  template <int dim>
  void
  ParsedFunction<dim>::declare_parameters(ParameterHandler & prm,
                                          const unsigned int n_components)
  {
    Assert(n_components > 0, ExcZero());

    std::string vnames;
    switch (dim)
      {
        case 1:
          vnames = "x,t";
          break;
        case 2:
          vnames = "x,y,t";
          break;
        case 3:
          vnames = "x,y,z,t";
          break;
        default:
          AssertThrow(false, ExcNotImplemented());
          break;
      }
    prm.declare_entry(
      "Variable names",
      vnames,
      Patterns::Anything(),
      "The names of the variables as they will be used in the "
      "function, separated by commas. By default, the names of variables "
      "at which the function will be evaluated are `x' (in 1d), `x,y' (in 2d) or "
      "`x,y,z' (in 3d) for spatial coordinates and `t' for time. You can then "
      "use these variable names in your function expression and they will be "
      "replaced by the values of these variables at which the function is "
      "currently evaluated. However, you can also choose a different set "
      "of names for the independent variables at which to evaluate your function "
      "expression. For example, if you work in spherical coordinates, you may "
      "wish to set this input parameter to `r,phi,theta,t' and then use these "
      "variable names in your function expression.");

    // The expression of the function
    std::string expr = "0";
    for (unsigned int i = 1; i < n_components; ++i)
      expr += "; 0";

    prm.declare_entry(
      "Function expression",
      expr,
      Patterns::Anything(),
      "The formula that denotes the function you want to evaluate for "
      "particular values of the independent variables. This expression "
      "may contain any of the usual operations such as addition or "
      "multiplication, as well as all of the common functions such as "
      "`sin' or `cos'. In addition, it may contain expressions like "
      "`if(x>0, 1, -1)' where the expression evaluates to the second "
      "argument if the first argument is true, and to the third argument "
      "otherwise. For a full overview of possible expressions accepted "
      "see the documentation of the muparser library at http://muparser.beltoforion.de/."
      "\n\n"
      "If the function you are describing represents a vector-valued "
      "function with multiple components, then separate the expressions "
      "for individual components by a semicolon.");
    prm.declare_entry(
      "Function constants",
      "",
      Patterns::Anything(),
      "Sometimes it is convenient to use symbolic constants in the "
      "expression that describes the function, rather than having to "
      "use its numeric value everywhere the constant appears. These "
      "values can be defined using this parameter, in the form "
      "`var1=value1, var2=value2, ...'."
      "\n\n"
      "A typical example would be to set this runtime parameter to "
      "`pi=3.1415926536' and then use `pi' in the expression of the "
      "actual formula. (That said, for convenience this class actually "
      "defines both `pi' and `Pi' by default, but you get the idea.)");
  }



  template <int dim>
  void
  ParsedFunction<dim>::parse_parameters(ParameterHandler &prm)
  {
    std::string vnames         = prm.get("Variable names");
    std::string expression     = prm.get("Function expression");
    std::string constants_list = prm.get("Function constants");

    std::vector<std::string> const_list =
      Utilities::split_string_list(constants_list, ',');
    std::map<std::string, double> constants;
    for (unsigned int i = 0; i < const_list.size(); ++i)
      {
        std::vector<std::string> this_c =
          Utilities::split_string_list(const_list[i], '=');
        AssertThrow(this_c.size() == 2, ExcMessage("Invalid format"));
        double tmp;
        AssertThrow(std::sscanf(this_c[1].c_str(), "%lf", &tmp),
                    ExcMessage("Double number?"));
        constants[this_c[0]] = tmp;
      }

    // set pi and Pi as synonyms for the corresponding value. note that
    // this overrides any value a user may have given
    constants["pi"] = numbers::PI;
    constants["Pi"] = numbers::PI;

    const unsigned int nn = (Utilities::split_string_list(vnames)).size();
    switch (nn)
      {
        case dim:
          // Time independent function
          function_object.initialize(vnames, expression, constants);
          break;
        case dim + 1:
          // Time dependent function
          function_object.initialize(vnames, expression, constants, true);
          break;
        default:
          AssertThrow(false,
                      ExcMessage(
                        "The list of variables specified is <" + vnames +
                        "> which is a list of length " +
                        Utilities::int_to_string(nn) +
                        " but it has to be a list of length equal to" +
                        " either dim (for a time-independent function)" +
                        " or dim+1 (for a time-dependent function)."));
      }
  }



  template <int dim>
  void
  ParsedFunction<dim>::vector_value(const Point<dim> &p,
                                    Vector<double> &  values) const
  {
    function_object.vector_value(p, values);
  }



  template <int dim>
  double
  ParsedFunction<dim>::value(const Point<dim> &p, unsigned int comp) const
  {
    return function_object.value(p, comp);
  }



  template <int dim>
  void
  ParsedFunction<dim>::set_time(const double newtime)
  {
    function_object.set_time(newtime);
    AutoDerivativeFunction<dim>::set_time(newtime);
  }


  // Explicit instantiations
  template class ParsedFunction<1>;
  template class ParsedFunction<2>;
  template class ParsedFunction<3>;
} // namespace Functions
DEAL_II_NAMESPACE_CLOSE