utilities.h

Go to the documentation of this file.

/**
 *  @file utilities.h
 *  Various templated functions that carry out common vector and polynomial operations
 *  (see @ref mathTemplates).
 */
 
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
 
/**
 * @defgroup mathTemplates Templated Array and Polynomial Operations
 *
 * These are templates to perform various simple operations on arrays. Note that
 * the compiler will inline these, so using them carries no performance penalty.
 *
 * @ingroup numerics
 */
 
#ifndef CT_UTILITIES_H
#define CT_UTILITIES_H
 
#include "ct_defs.h"
#include "cantera/base/ctexceptions.h"
#include <numeric>
#include <iterator>
 
namespace Cantera
{
//! @addtogroup mathTemplates
//! @{
 
//! Templated Inner product of two vectors of length 4.
/*!
 * If either @e x or @e y has length greater than 4, only the first 4 elements
 * will be used.
 *
 * @param x   first reference to the templated class V
 * @param y   second reference to the templated class V
 * @return This class returns a hard-coded type, double.
 */
template<class X, class Y>
inline auto dot4(const X& x, const Y& y)
{
    auto x_ = span{std::data(x), std::size(x)};
    auto y_ = span{std::data(y), std::size(y)};
#ifndef NDEBUG
    checkArraySize("dot4: x", x_.size(), 4);
    checkArraySize("dot4: y", y_.size(), 4);
#endif
    return x_[0]*y_[0] + x_[1]*y_[1] + x_[2]*y_[2] + x_[3]*y_[3];
}
 
//! Templated Inner product of two vectors of length 5
/*!
 * If either @e x or @e y has length greater than 4, only the first 4 elements
 * will be used.
 *
 * @param x   first reference to the templated class V
 * @param y   second reference to the templated class V
 * @return This class returns a hard-coded type, double.
 */
template<class X, class Y>
inline auto dot5(const X& x, const Y& y)
{
    auto x_ = span{std::data(x), std::size(x)};
    auto y_ = span{std::data(y), std::size(y)};
#ifndef NDEBUG
    checkArraySize("dot5: x", x_.size(), 5);
    checkArraySize("dot5: y", y_.size(), 5);
#endif
    return x_[0]*y_[0] + x_[1]*y_[1] + x_[2]*y_[2] + x_[3]*y_[3] +
           x_[4]*y_[4];
}
 
//! Function that calculates a templated inner product.
/*!
 * This inner product is templated twice. The output variable is hard coded
 * to return a double.
 *
 * template<class InputIter, class InputIter2>
 *
 * @code
 *     double x[8], y[8];
 *     double dsum = dot<double *,double *>(x, &x+7, y);
 * @endcode
 *
 * @param x_begin  Iterator pointing to the beginning, belonging to the
 *                 iterator class InputIter.
 * @param x_end    Iterator pointing to the end, belonging to the
 *                 iterator class InputIter.
 * @param y_begin Iterator pointing to the beginning of y, belonging to the
 *               iterator class InputIter2.
 * @return The return is hard-coded to return a double.
 */
template<class InputIter, class InputIter2>
inline double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
{
    return std::inner_product(x_begin, x_end, y_begin, 0.0);
}
 
//! Multiply elements of an array by a scale factor.
/*!
 * @code
 * vector<double> in(8, 1.0), out(8);
 * scale(in.begin(), in.end(), out.begin(), factor);
 * @endcode
 *
 * @param begin  Iterator pointing to the beginning, belonging to the
 *               iterator class InputIter.
 * @param end    Iterator pointing to the end, belonging to the
 *               iterator class InputIter.
 * @param out    Iterator pointing to the beginning of out, belonging to the
 *               iterator class OutputIter. This is the output variable
 *               for this routine.
 * @param scale_factor  input scale factor belonging to the class S.
 */
template<class InputIter, class OutputIter, class S>
inline void scale(InputIter begin, InputIter end,
                  OutputIter out, S scale_factor)
{
    std::transform(begin, end, out,
        [scale_factor](double x) { return x * scale_factor; });
}
 
//! Templated evaluation of a polynomial of order 6
/*!
 * @param x   Value of the independent variable - First template parameter
 * @param c   Pointer to the polynomial - Second template parameter
 */
template<class D, class C>
auto poly6(D x, const C& c)
{
    auto c_ = span{std::data(c), std::size(c)};
#ifndef NDEBUG
    checkArraySize("poly6: c", c_.size(), 7);
#endif
    return ((((((c_[6]*x + c_[5])*x + c_[4])*x + c_[3])*x +
              c_[2])*x + c_[1])*x + c_[0]);
}
 
//! Templated evaluation of a polynomial of order 8
/*!
 * @param x   Value of the independent variable - First template parameter
 * @param c   Pointer to the polynomial - Second template parameter
 */
template<class D, class C>
auto poly8(D x, const C& c)
{
    auto c_ = span{std::data(c), std::size(c)};
#ifndef NDEBUG
    checkArraySize("poly8: c", c_.size(), 9);
#endif
    return ((((((((c_[8]*x + c_[7])*x + c_[6])*x + c_[5])*x + c_[4])*x + c_[3])*x +
              c_[2])*x + c_[1])*x + c_[0]);
}
 
//! Templated evaluation of a polynomial of order 5
/*!
 * @param x   Value of the independent variable - First template parameter
 * @param c   Pointer to the polynomial - Second template parameter
 */
template<class D, class C>
auto poly5(D x, const C& c)
{
    auto c_ = span{std::data(c), std::size(c)};
#ifndef NDEBUG
    checkArraySize("poly5: c", c_.size(), 6);
#endif
    return (((((c_[5]*x + c_[4])*x + c_[3])*x +
              c_[2])*x + c_[1])*x + c_[0]);
}
 
//! Evaluates a polynomial of order 4.
/*!
 * @param x   Value of the independent variable.
 * @param c   Pointer to the polynomial coefficient array.
 */
template<class D, class C>
auto poly4(D x, const C& c)
{
    auto c_ = span{std::data(c), std::size(c)};
#ifndef NDEBUG
    checkArraySize("poly4: c", c_.size(), 5);
#endif
    return ((((c_[4]*x + c_[3])*x +
              c_[2])*x + c_[1])*x + c_[0]);
}
 
//! Templated evaluation of a polynomial of order 3
/*!
 * @param x   Value of the independent variable - First template parameter
 * @param c   Pointer to the polynomial - Second template parameter
 */
template<class D, class C>
auto poly3(D x, const C& c)
{
    auto c_ = span{std::data(c), std::size(c)};
#ifndef NDEBUG
    checkArraySize("poly3: c", c_.size(), 4);
#endif
    return (((c_[3]*x + c_[2])*x + c_[1])*x + c_[0]);
}
 
//! @}
 
//! Check to see that a number is finite (not NaN, +Inf or -Inf)
void checkFinite(const double tmp);
 
//! Check to see that all elements in an array are finite
/*!
 * Throws an exception if any element is NaN, +Inf, or -Inf
 * @param name    Name to be used in the exception message if the check fails
 * @param values  Array of values to be checked
 */
void checkFinite(const string& name, span<const double> values);
 
//! Const accessor for a value in a map.
/*
 * Similar to map.at(key), but returns *default_val* if the key is not
 * found instead of throwing an exception.
 */
template <class T, class U>
const U& getValue(const map<T, U>& m, const T& key, const U& default_val) {
    typename map<T,U>::const_iterator iter = m.find(key);
    return (iter == m.end()) ? default_val : iter->second;
}
 
//! Get the size of a container, cast to a signed integer type
//! @ingroup mathTemplates
template <class T, class U=int>
U len(const T& container) {
    return static_cast<U>(container.size());
}
 
//! A macro for generating member function detectors, which can then be used in
//! combination with `if constexpr` to condition behavior on the availability of that
//! member function. See MultiRate for examples of use.
#define CT_DEFINE_HAS_MEMBER(detector_name, func_name) \
    template<class T, class=void> \
    struct detector_name : std::false_type {}; \
    template<class T> \
    struct detector_name<T, std::void_t<decltype(&T::func_name)>> : std::true_type {};
 
}
 
#endif