MathUtils.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include "Moose.h"
#include "MooseError.h"
#include "MooseTypes.h"
#include "MooseUtils.h"
#include "libmesh/libmesh.h"
#include "libmesh/utility.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/compare_types.h"
#include "libmesh/point.h"
 
namespace MathUtils
{
 
static constexpr Real sqrt2 = 1.4142135623730951;
 
Real poly1Log(Real x, Real tol, unsigned int derivative_order);
Real poly2Log(Real x, Real tol, unsigned int derivative_order);
Real poly3Log(Real x, Real tol, unsigned int derivative_order);
Real poly4Log(Real x, Real tol, unsigned int derivative_order);
Real taylorLog(Real x);
Point barycentricToCartesian2D(const Point & p0,
                               const Point & p1,
                               const Point & p2,
                               const Real b0,
                               const Real b1,
                               const Real b2);
Point barycentricToCartesian3D(const Point & p0,
                               const Point & p1,
                               const Point & p2,
                               const Point & p3,
                               const Real b0,
                               const Real b1,
                               const Real b2,
                               const Real b3);
Point circumcenter2D(const Point & p0, const Point & p1, const Point & p2);
Point circumcenter3D(const Point & p0, const Point & p1, const Point & p2, const Point & p3);
 
template <typename T>
T
round(T x)
{
  return ::round(x); // use round from math.h
}
 
template <typename T>
T
sign(T x)
{
  return x >= 0.0 ? 1.0 : -1.0;
}
 
template <typename T>
T
pow(T x, int e)
{
  bool neg = false;
  T result = 1.0;
 
  if (e < 0)
  {
    neg = true;
    e = -e;
  }
 
  while (e)
  {
    // if bit 0 is set multiply the current power of two factor of the exponent
    if (e & 1)
      result *= x;
 
    // x is incrementally set to consecutive powers of powers of two
    x *= x;
 
    // bit shift the exponent down
    e >>= 1;
  }
 
  return neg ? 1.0 / result : result;
}
 
template <typename T>
T
heavyside(T x)
{
  return x < 0.0 ? 0.0 : 1.0;
}
 
template <typename T>
T
regularizedHeavyside(const T & x, Real smoothing_length)
{
  if (x <= -smoothing_length)
    return 0.0;
  else if (x < smoothing_length)
  {
    using std::sin;
    return 0.5 * (1 + sin(libMesh::pi * x / 2 / smoothing_length));
  }
  else
    return 1.0;
}
 
template <typename T>
T
regularizedHeavysideDerivative(const T & x, Real smoothing_length)
{
  if (x < smoothing_length && x > -smoothing_length)
  {
    using std::cos;
    return 0.25 * libMesh::pi / smoothing_length * (cos(libMesh::pi * x / 2 / smoothing_length));
  }
  else
    return 0.0;
}
 
template <typename T>
T
positivePart(T x)
{
  return x > 0.0 ? x : 0.0;
}
 
template <typename T>
T
negativePart(T x)
{
  return x < 0.0 ? x : 0.0;
}
 
template <
    typename T,
    typename T2,
    typename T3,
    typename std::enable_if<libMesh::ScalarTraits<T>::value && libMesh::ScalarTraits<T2>::value &&
                                libMesh::ScalarTraits<T3>::value,
                            int>::type = 0>
void
addScaled(const T & a, const T2 & b, T3 & result)
{
  result += a * b;
}
 
template <typename T,
          typename T2,
          typename T3,
          typename std::enable_if<libMesh::ScalarTraits<T>::value, int>::type = 0>
void
addScaled(const T & scalar,
          const libMesh::NumericVector<T2> & numeric_vector,
          libMesh::NumericVector<T3> & result)
{
  result.add(scalar, numeric_vector);
}
 
template <
    typename T,
    typename T2,
    template <typename> class W,
    template <typename> class W2,
    typename std::enable_if<std::is_same<typename W<T>::index_type, unsigned int>::value &&
                                std::is_same<typename W2<T2>::index_type, unsigned int>::value,
                            int>::type = 0>
typename libMesh::CompareTypes<T, T2>::supertype
dotProduct(const W<T> & a, const W2<T2> & b)
{
  return a * b;
}
 
template <typename T,
          typename T2,
          template <typename> class W,
          template <typename> class W2,
          typename std::enable_if<std::is_same<typename W<T>::index_type,
                                               std::tuple<unsigned int, unsigned int>>::value &&
                                      std::is_same<typename W2<T2>::index_type,
                                                   std::tuple<unsigned int, unsigned int>>::value,
                                  int>::type = 0>
typename libMesh::CompareTypes<T, T2>::supertype
dotProduct(const W<T> & a, const W2<T2> & b)
{
  return a.contract(b);
}
 
template <typename T>
auto
normSquared(const T & value)
{
  if constexpr (MooseUtils::Has_size<T>::value)
  {
    decltype(normSquared(*value.begin())) sum{};
    for (const auto & component : value)
      sum += normSquared(component);
    return sum;
  }
  else if constexpr (libMesh::TensorTools::TensorTraits<T>::rank == 0)
    return value * value;
  else if constexpr (libMesh::TensorTools::TensorTraits<T>::rank == 1)
  {
    auto sum = value(0) * value(0);
    for (const auto i : make_range(std::size_t(1), Moose::dim))
      sum += value(i) * value(i);
    return sum;
  }
  else
  {
    auto sum = value(0, 0) * value(0, 0);
    for (const auto j : make_range(std::size_t(1), Moose::dim))
      sum += value(0, j) * value(0, j);
    for (const auto i : make_range(std::size_t(1), Moose::dim))
      for (const auto j : make_range(Moose::dim))
        sum += value(i, j) * value(i, j);
    return sum;
  }
}
 
template <typename T>
auto
norm(const T & value)
{
  using std::sqrt;
  if (MooseUtils::isZero(value))
    return decltype(normSquared(value)){};
  return sqrt(normSquared(value));
}
 
template <typename C,
          typename T,
          typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
R
poly(const C & c, const T x, const bool derivative = false)
{
  const auto size = c.size();
  if (size == 0)
    return 0.0;
 
  R value = c[0];
  if (derivative)
  {
    value *= size - 1;
    for (std::size_t i = 1; i < size - 1; ++i)
      value = value * x + c[i] * (size - i - 1);
  }
  else
  {
    for (std::size_t i = 1; i < size; ++i)
      value = value * x + c[i];
  }
 
  return value;
}
 
template <typename C,
          typename T,
          typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
R
polynomial(const C & c, const T x)
{
  auto size = c.size();
  if (size == 0)
    return 0.0;
 
  size--;
  R value = c[size];
  for (std::size_t i = 1; i <= size; ++i)
    value = value * x + c[size - i];
 
  return value;
}
 
template <typename C,
          typename T,
          typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
R
polynomialDerivative(const C & c, const T x)
{
  auto size = c.size();
  if (size <= 1)
    return 0.0;
 
  size--;
  R value = c[size] * size;
  for (std::size_t i = 1; i < size; ++i)
    value = value * x + c[size - i] * (size - i);
 
  return value;
}
 
template <typename T, typename T2>
T
clamp(const T & x, T2 lowerlimit, T2 upperlimit)
{
  if (x < lowerlimit)
    return lowerlimit;
  if (x > upperlimit)
    return upperlimit;
  return x;
}
 
template <typename T, typename T2>
T
smootherStep(T x, T2 start, T2 end, bool derivative = false)
{
  mooseAssert("start < end", "Start value must be lower than end value for smootherStep");
  if (x <= start)
    return 0.0;
  else if (x >= end)
  {
    if (derivative)
      return 0.0;
    else
      return 1.0;
  }
  x = (x - start) / (end - start);
  if (derivative)
    return 30.0 * libMesh::Utility::pow<2>(x) * (x * (x - 2.0) + 1.0) / (end - start);
  return libMesh::Utility::pow<3>(x) * (x * (x * 6.0 - 15.0) + 10.0);
}
 
enum class ComputeType
{
  value,
  derivative
};
 
template <ComputeType compute_type, typename X, typename S, typename E>
auto
smootherStep(const X & x, const S & start, const E & end)
{
  mooseAssert("start < end", "Start value must be lower than end value for smootherStep");
  if (x <= start)
    return 0.0;
  else if (x >= end)
  {
    if constexpr (compute_type == ComputeType::derivative)
      return 0.0;
    if constexpr (compute_type == ComputeType::value)
      return 1.0;
  }
  const auto u = (x - start) / (end - start);
  if constexpr (compute_type == ComputeType::derivative)
    return 30.0 * libMesh::Utility::pow<2>(u) * (u * (u - 2.0) + 1.0) / (end - start);
  if constexpr (compute_type == ComputeType::value)
    return libMesh::Utility::pow<3>(u) * (u * (u * 6.0 - 15.0) + 10.0);
}
 
template <typename T>
inline void
mooseSetToZero(T & v)
{
  v = 0;
}
template <typename T>
inline void
mooseSetToZero(T *&)
{
  mooseError("mooseSetToZero does not accept pointers");
}
 
template <>
inline void
mooseSetToZero(std::vector<Real> & vec)
{
  for (auto & v : vec)
    v = 0.;
}
 
std::vector<std::vector<unsigned int>> multiIndex(unsigned int dim, unsigned int order);
 
template <ComputeType compute_type, typename X, typename X1, typename X2, typename Y1, typename Y2>
auto
linearInterpolation(const X & x, const X1 & x1, const X2 & x2, const Y1 & y1, const Y2 & y2)
{
  const auto m = (y2 - y1) / (x2 - x1);
  if constexpr (compute_type == ComputeType::derivative)
    return m;
  if constexpr (compute_type == ComputeType::value)
    return m * (x - x1) + y1;
}
 
template <typename T1, typename T2>
std::size_t
euclideanMod(T1 dividend, T2 divisor)
{
  return (dividend % divisor + divisor) % divisor;
}
 
template <typename T>
T
gradName(const T & base_prop_name)
{
  return "grad_" + base_prop_name;
}
 
template <typename T>
T
timeDerivName(const T & base_prop_name)
{
  return "d" + base_prop_name + "_dt";
}
 
void kron(RealEigenMatrix & product, const RealEigenMatrix & mat_A, const RealEigenMatrix & mat_B);
 
} // namespace MathUtils
 
std::vector<std::vector<unsigned int>> multiIndexHelper(unsigned int N, unsigned int K);