https://mooseframework.inl.gov
MathUtils.h
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #pragma once
11 
12 #include "Moose.h"
13 #include "MooseError.h"
14 #include "MooseTypes.h"
15 #include "libmesh/libmesh.h"
16 #include "libmesh/utility.h"
17 #include "libmesh/numeric_vector.h"
18 #include "libmesh/compare_types.h"
19 #include "libmesh/point.h"
20 
21 namespace MathUtils
22 {
23 
25 static constexpr Real sqrt2 = 1.4142135623730951;
26 
27 Real poly1Log(Real x, Real tol, unsigned int derivative_order);
28 Real poly2Log(Real x, Real tol, unsigned int derivative_order);
29 Real poly3Log(Real x, Real tol, unsigned int derivative_order);
30 Real poly4Log(Real x, Real tol, unsigned int derivative_order);
31 Real taylorLog(Real x);
40 Point barycentricToCartesian2D(const Point & p0,
41  const Point & p1,
42  const Point & p2,
43  const Real b0,
44  const Real b1,
45  const Real b2);
54 Point barycentricToCartesian3D(const Point & p0,
55  const Point & p1,
56  const Point & p2,
57  const Point & p3,
58  const Real b0,
59  const Real b1,
60  const Real b2,
61  const Real b3);
67 Point circumcenter2D(const Point & p0, const Point & p1, const Point & p2);
73 Point circumcenter3D(const Point & p0, const Point & p1, const Point & p2, const Point & p3);
74 
75 template <typename T>
76 T
77 round(T x)
78 {
79  return ::round(x); // use round from math.h
80 }
81 
82 template <typename T>
83 T
84 sign(T x)
85 {
86  return x >= 0.0 ? 1.0 : -1.0;
87 }
88 
89 template <typename T>
90 T
91 pow(T x, int e)
92 {
93  bool neg = false;
94  T result = 1.0;
95 
96  if (e < 0)
97  {
98  neg = true;
99  e = -e;
100  }
101 
102  while (e)
103  {
104  // if bit 0 is set multiply the current power of two factor of the exponent
105  if (e & 1)
106  result *= x;
107 
108  // x is incrementally set to consecutive powers of powers of two
109  x *= x;
110 
111  // bit shift the exponent down
112  e >>= 1;
113  }
114 
115  return neg ? 1.0 / result : result;
116 }
117 
118 template <typename T>
119 T
121 {
122  return x < 0.0 ? 0.0 : 1.0;
123 }
124 
125 template <typename T>
126 T
127 regularizedHeavyside(T x, Real smoothing_length)
128 {
129  if (x <= -smoothing_length)
130  return 0.0;
131  else if (x < smoothing_length)
132  return 0.5 * (1 + std::sin(libMesh::pi * x / 2 / smoothing_length));
133  else
134  return 1.0;
135 }
136 
137 template <typename T>
138 T
139 regularizedHeavysideDerivative(T x, Real smoothing_length)
140 {
141  if (x < smoothing_length && x > -smoothing_length)
142  return 0.25 * libMesh::pi / smoothing_length *
143  (std::cos(libMesh::pi * x / 2 / smoothing_length));
144  else
145  return 0.0;
146 }
147 
148 template <typename T>
149 T
151 {
152  return x > 0.0 ? x : 0.0;
153 }
154 
155 template <typename T>
156 T
158 {
159  return x < 0.0 ? x : 0.0;
160 }
161 
162 template <
163  typename T,
164  typename T2,
165  typename T3,
166  typename std::enable_if<libMesh::ScalarTraits<T>::value && libMesh::ScalarTraits<T2>::value &&
168  int>::type = 0>
169 void
170 addScaled(const T & a, const T2 & b, T3 & result)
171 {
172  result += a * b;
173 }
174 
175 template <typename T,
176  typename T2,
177  typename T3,
178  typename std::enable_if<libMesh::ScalarTraits<T>::value, int>::type = 0>
179 void
180 addScaled(const T & scalar,
181  const libMesh::NumericVector<T2> & numeric_vector,
183 {
184  result.add(scalar, numeric_vector);
185 }
186 
187 template <
188  typename T,
189  typename T2,
190  template <typename>
191  class W,
192  template <typename>
193  class W2,
194  typename std::enable_if<std::is_same<typename W<T>::index_type, unsigned int>::value &&
195  std::is_same<typename W2<T2>::index_type, unsigned int>::value,
196  int>::type = 0>
198 dotProduct(const W<T> & a, const W2<T2> & b)
199 {
200  return a * b;
201 }
202 
203 template <typename T,
204  typename T2,
205  template <typename>
206  class W,
207  template <typename>
208  class W2,
209  typename std::enable_if<std::is_same<typename W<T>::index_type,
210  std::tuple<unsigned int, unsigned int>>::value &&
211  std::is_same<typename W2<T2>::index_type,
212  std::tuple<unsigned int, unsigned int>>::value,
213  int>::type = 0>
215 dotProduct(const W<T> & a, const W2<T2> & b)
216 {
217  return a.contract(b);
218 }
219 
233 template <typename C,
234  typename T,
236 R
237 poly(const C & c, const T x, const bool derivative = false)
238 {
239  const auto size = c.size();
240  if (size == 0)
241  return 0.0;
242 
243  R value = c[0];
244  if (derivative)
245  {
246  value *= size - 1;
247  for (std::size_t i = 1; i < size - 1; ++i)
248  value = value * x + c[i] * (size - i - 1);
249  }
250  else
251  {
252  for (std::size_t i = 1; i < size; ++i)
253  value = value * x + c[i];
254  }
255 
256  return value;
257 }
258 
268 template <typename C,
269  typename T,
271 R
272 polynomial(const C & c, const T x)
273 {
274  auto size = c.size();
275  if (size == 0)
276  return 0.0;
277 
278  size--;
279  R value = c[size];
280  for (std::size_t i = 1; i <= size; ++i)
281  value = value * x + c[size - i];
282 
283  return value;
284 }
285 
289 template <typename C,
290  typename T,
292 R
293 polynomialDerivative(const C & c, const T x)
294 {
295  auto size = c.size();
296  if (size <= 1)
297  return 0.0;
298 
299  size--;
300  R value = c[size] * size;
301  for (std::size_t i = 1; i < size; ++i)
302  value = value * x + c[size - i] * (size - i);
303 
304  return value;
305 }
306 
307 template <typename T, typename T2>
308 T
309 clamp(const T & x, T2 lowerlimit, T2 upperlimit)
310 {
311  if (x < lowerlimit)
312  return lowerlimit;
313  if (x > upperlimit)
314  return upperlimit;
315  return x;
316 }
317 
318 template <typename T, typename T2>
319 T
320 smootherStep(T x, T2 start, T2 end, bool derivative = false)
321 {
322  mooseAssert("start < end", "Start value must be lower than end value for smootherStep");
323  if (x <= start)
324  return 0.0;
325  else if (x >= end)
326  {
327  if (derivative)
328  return 0.0;
329  else
330  return 1.0;
331  }
332  x = (x - start) / (end - start);
333  if (derivative)
334  return 30.0 * libMesh::Utility::pow<2>(x) * (x * (x - 2.0) + 1.0) / (end - start);
335  return libMesh::Utility::pow<3>(x) * (x * (x * 6.0 - 15.0) + 10.0);
336 }
337 
338 enum class ComputeType
339 {
340  value,
341  derivative
342 };
343 
344 template <ComputeType compute_type, typename X, typename S, typename E>
345 auto
346 smootherStep(const X & x, const S & start, const E & end)
347 {
348  mooseAssert("start < end", "Start value must be lower than end value for smootherStep");
349  if (x <= start)
350  return 0.0;
351  else if (x >= end)
352  {
353  if constexpr (compute_type == ComputeType::derivative)
354  return 0.0;
355  if constexpr (compute_type == ComputeType::value)
356  return 1.0;
357  }
358  const auto u = (x - start) / (end - start);
359  if constexpr (compute_type == ComputeType::derivative)
360  return 30.0 * libMesh::Utility::pow<2>(u) * (u * (u - 2.0) + 1.0) / (end - start);
361  if constexpr (compute_type == ComputeType::value)
362  return libMesh::Utility::pow<3>(u) * (u * (u * 6.0 - 15.0) + 10.0);
363 }
364 
370 template <typename T>
371 inline void
373 {
379  v = 0;
380 }
381 template <typename T>
382 inline void
384 {
385  mooseError("mooseSetToZero does not accept pointers");
386 }
387 
388 template <>
389 inline void
390 mooseSetToZero(std::vector<Real> & vec)
391 {
392  for (auto & v : vec)
393  v = 0.;
394 }
395 
413 std::vector<std::vector<unsigned int>> multiIndex(unsigned int dim, unsigned int order);
414 
415 template <ComputeType compute_type, typename X, typename X1, typename X2, typename Y1, typename Y2>
416 auto
417 linearInterpolation(const X & x, const X1 & x1, const X2 & x2, const Y1 & y1, const Y2 & y2)
418 {
419  const auto m = (y2 - y1) / (x2 - x1);
420  if constexpr (compute_type == ComputeType::derivative)
421  return m;
422  if constexpr (compute_type == ComputeType::value)
423  return m * (x - x1) + y1;
424 }
425 
432 template <typename T1, typename T2>
433 std::size_t
434 euclideanMod(T1 dividend, T2 divisor)
435 {
436  return (dividend % divisor + divisor) % divisor;
437 }
438 
443 template <typename T>
444 T
445 gradName(const T & base_prop_name)
446 {
447  return "grad_" + base_prop_name;
448 }
449 
454 template <typename T>
455 T
456 timeDerivName(const T & base_prop_name)
457 {
458  return "d" + base_prop_name + "_dt";
459 }
460 
467 void kron(RealEigenMatrix & product, const RealEigenMatrix & mat_A, const RealEigenMatrix & mat_B);
468 
469 } // namespace MathUtils
470 
472 std::vector<std::vector<unsigned int>> multiIndexHelper(unsigned int N, unsigned int K);
void addScaled(const T &a, const T2 &b, T3 &result)
Definition: MathUtils.h:170
void mooseSetToZero(T &v)
Helper function templates to set a variable to zero.
Definition: MathUtils.h:372
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sin(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tan
R poly(const C &c, const T x, const bool derivative=false)
Evaluate a polynomial with the coefficients c at x.
Definition: MathUtils.h:237
Point circumcenter2D(const Point &p0, const Point &p1, const Point &p2)
Evaluate circumcenter of a triangle given three arbitrary points.
Definition: MathUtils.C:260
Point barycentricToCartesian2D(const Point &p0, const Point &p1, const Point &p2, const Real b0, const Real b1, const Real b2)
Evaluate Cartesian coordinates of any center point of a triangle given Barycentric coordinates of cen...
Definition: MathUtils.C:217
void kron(RealEigenMatrix &product, const RealEigenMatrix &mat_A, const RealEigenMatrix &mat_B)
Computes the Kronecker product of two matrices.
Definition: MathUtils.C:17
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
T regularizedHeavysideDerivative(T x, Real smoothing_length)
Definition: MathUtils.h:139
Real poly1Log(Real x, Real tol, unsigned int derivative_order)
Definition: MathUtils.C:49
Real poly2Log(Real x, Real tol, unsigned int derivative_order)
Definition: MathUtils.C:76
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:154
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)
Evaluate Cartesian coordinates of any center point of a tetrahedron given Barycentric coordinates of ...
Definition: MathUtils.C:237
std::size_t euclideanMod(T1 dividend, T2 divisor)
perform modulo operator for Euclidean division that ensures a non-negative result ...
Definition: MathUtils.h:434
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
T round(T x)
Definition: MathUtils.h:77
libMesh::CompareTypes< T, T2 >::supertype dotProduct(const W< T > &a, const W2< T2 > &b)
Definition: MathUtils.h:198
T sign(T x)
Definition: MathUtils.h:84
Real taylorLog(Real x)
Definition: MathUtils.C:172
R polynomial(const C &c, const T x)
Evaluate a polynomial with the coefficients c at x.
Definition: MathUtils.h:272
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > RealEigenMatrix
Definition: MooseTypes.h:149
Real poly3Log(Real x, Real tol, unsigned int derivative_order)
Definition: MathUtils.C:104
T negativePart(T x)
Definition: MathUtils.h:157
T smootherStep(T x, T2 start, T2 end, bool derivative=false)
Definition: MathUtils.h:320
Point circumcenter3D(const Point &p0, const Point &p1, const Point &p2, const Point &p3)
Evaluate circumcenter of a tetrahedrom given four arbitrary points.
Definition: MathUtils.C:289
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T gradName(const T &base_prop_name)
automatic prefixing for naming material properties based on gradients of coupled variables/functors ...
Definition: MathUtils.h:445
T timeDerivName(const T &base_prop_name)
automatic prefixing for naming material properties based on time derivatives of coupled variables/fun...
Definition: MathUtils.h:456
T positivePart(T x)
Definition: MathUtils.h:150
std::vector< std::vector< unsigned int > > multiIndex(unsigned int dim, unsigned int order)
generate a complete multi index table for given dimension and order i.e.
Definition: MathUtils.C:186
T heavyside(T x)
Definition: MathUtils.h:120
Real poly4Log(Real x, Real tol, unsigned int derivative_order)
Definition: MathUtils.C:133
T regularizedHeavyside(T x, Real smoothing_length)
Definition: MathUtils.h:127
T clamp(const T &x, T2 lowerlimit, T2 upperlimit)
Definition: MathUtils.h:309
T pow(T x, int e)
Definition: MathUtils.h:91
virtual void add(const numeric_index_type i, const T value)=0
std::vector< std::vector< unsigned int > > multiIndexHelper(unsigned int N, unsigned int K)
A helper function for MathUtils::multiIndex.
Definition: MathUtils.C:333
static constexpr Real sqrt2
std::sqrt is not constexpr, so we add sqrt(2) as a constant (used in Mandel notation) ...
Definition: MathUtils.h:25
auto linearInterpolation(const X &x, const X1 &x1, const X2 &x2, const Y1 &y1, const Y2 &y2)
Definition: MathUtils.h:417
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template pow< 2 >(tan(_arg))+1.0) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(sqrt
R polynomialDerivative(const C &c, const T x)
Returns the derivative of polynomial(c, x) with respect to x.
Definition: MathUtils.h:293
const Real pi