Line data Source code
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 :
24 : /// std::sqrt is not constexpr, so we add sqrt(2) as a constant (used in Mandel notation)
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);
32 : /**
33 : * Evaluate Cartesian coordinates of any center point of a triangle given Barycentric
34 : * coordinates of center point and Cartesian coordinates of triangle's vertices
35 : * @param p0,p1,p2 are the three non-collinear vertices in Cartesian coordinates
36 : * @param b0,b1,b2 is the center point in barycentric coordinates with b0+b1+b2=1, e.g.
37 : * (1/3,1/3,1/3) for a centroid
38 : * @return the center point of triangle in Cartesian coordinates
39 : */
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);
46 : /**
47 : * Evaluate Cartesian coordinates of any center point of a tetrahedron given Barycentric
48 : * coordinates of center point and Cartesian coordinates of tetrahedon's vertices
49 : * @param p0,p1,p2,p3 are the three non-coplanar vertices in Cartesian coordinates
50 : * @param b0,b1,b2,b3 is the center point in barycentric coordinates with b0+b1+b2+b3=1, e.g.
51 : * (1/4,1/4,1/4,1/4) for a centroid.
52 : * @return the center point of tetrahedron in Cartesian coordinates
53 : */
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);
62 : /**
63 : * Evaluate circumcenter of a triangle given three arbitrary points
64 : * @param p0,p1,p2 are the three non-collinear vertices in Cartesian coordinates
65 : * @return the circumcenter in Cartesian coordinates
66 : */
67 : Point circumcenter2D(const Point & p0, const Point & p1, const Point & p2);
68 : /**
69 : * Evaluate circumcenter of a tetrahedrom given four arbitrary points
70 : * @param p0,p1,p2,p3 are the four non-coplanar vertices in Cartesian coordinates
71 : * @return the circumcenter in Cartesian coordinates
72 : */
73 : Point circumcenter3D(const Point & p0, const Point & p1, const Point & p2, const Point & p3);
74 :
75 : template <typename T>
76 : T
77 8 : round(T x)
78 : {
79 8 : return ::round(x); // use round from math.h
80 : }
81 :
82 : template <typename T>
83 : T
84 38080 : sign(T x)
85 : {
86 38080 : return x >= 0.0 ? 1.0 : -1.0;
87 : }
88 :
89 : template <typename T>
90 : T
91 547153991 : pow(T x, int e)
92 : {
93 547153991 : bool neg = false;
94 547153991 : T result = 1.0;
95 :
96 547153991 : if (e < 0)
97 : {
98 2 : neg = true;
99 2 : e = -e;
100 : }
101 :
102 1093889316 : while (e)
103 : {
104 : // if bit 0 is set multiply the current power of two factor of the exponent
105 546735325 : if (e & 1)
106 364560383 : result *= x;
107 :
108 : // x is incrementally set to consecutive powers of powers of two
109 546735325 : x *= x;
110 :
111 : // bit shift the exponent down
112 546735325 : e >>= 1;
113 : }
114 :
115 547154215 : return neg ? 1.0 / result : result;
116 224 : }
117 :
118 : template <typename T>
119 : T
120 : heavyside(T x)
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
150 : positivePart(T x)
151 : {
152 : return x > 0.0 ? x : 0.0;
153 : }
154 :
155 : template <typename T>
156 : T
157 : negativePart(T x)
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 &&
167 : libMesh::ScalarTraits<T3>::value,
168 : int>::type = 0>
169 : void
170 4394576 : addScaled(const T & a, const T2 & b, T3 & result)
171 : {
172 4394576 : result += a * b;
173 4394576 : }
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 70367 : addScaled(const T & scalar,
181 : const libMesh::NumericVector<T2> & numeric_vector,
182 : libMesh::NumericVector<T3> & result)
183 : {
184 70367 : result.add(scalar, numeric_vector);
185 70367 : }
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>
197 : typename libMesh::CompareTypes<T, T2>::supertype
198 124284548 : dotProduct(const W<T> & a, const W2<T2> & b)
199 : {
200 124284548 : 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>
214 : typename libMesh::CompareTypes<T, T2>::supertype
215 0 : dotProduct(const W<T> & a, const W2<T2> & b)
216 : {
217 0 : return a.contract(b);
218 : }
219 :
220 : /**
221 : * Evaluate a polynomial with the coefficients c at x. Note that the Polynomial
222 : * form is
223 : * c[0]*x^s + c[1]*x^(s-1) + c[2]*x^(s-2) + ... + c[s-2]*x^2 + c[s-1]*x + c[s]
224 : * where s = c.size()-1 , which is counter intuitive!
225 : *
226 : * This function will be DEPRECATED soon (10/22/2020)
227 : *
228 : * The coefficient container type can be any container that provides an index
229 : * operator [] and a .size() method (e.g. std::vector, std::array). The return
230 : * type is the supertype of the container value type and the argument x.
231 : * The supertype is the type that can represent both number types.
232 : */
233 : template <typename C,
234 : typename T,
235 : typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
236 : R
237 324 : poly(const C & c, const T x, const bool derivative = false)
238 : {
239 324 : const auto size = c.size();
240 324 : if (size == 0)
241 0 : return 0.0;
242 :
243 324 : R value = c[0];
244 324 : if (derivative)
245 : {
246 82 : value *= size - 1;
247 808 : for (std::size_t i = 1; i < size - 1; ++i)
248 726 : value = value * x + c[i] * (size - i - 1);
249 : }
250 : else
251 : {
252 1210 : for (std::size_t i = 1; i < size; ++i)
253 968 : value = value * x + c[i];
254 : }
255 :
256 324 : return value;
257 : }
258 :
259 : /**
260 : * Evaluate a polynomial with the coefficients c at x. Note that the Polynomial
261 : * form is
262 : * c[0] + c[1] * x + c[2] * x^2 + ...
263 : * The coefficient container type can be any container that provides an index
264 : * operator [] and a .size() method (e.g. std::vector, std::array). The return
265 : * type is the supertype of the container value type and the argument x.
266 : * The supertype is the type that can represent both number types.
267 : */
268 : template <typename C,
269 : typename T,
270 : typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
271 : R
272 5 : polynomial(const C & c, const T x)
273 : {
274 5 : auto size = c.size();
275 5 : if (size == 0)
276 0 : return 0.0;
277 :
278 5 : size--;
279 5 : R value = c[size];
280 27 : for (std::size_t i = 1; i <= size; ++i)
281 22 : value = value * x + c[size - i];
282 :
283 5 : return value;
284 : }
285 :
286 : /**
287 : * Returns the derivative of polynomial(c, x) with respect to x
288 : */
289 : template <typename C,
290 : typename T,
291 : typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
292 : R
293 5 : polynomialDerivative(const C & c, const T x)
294 : {
295 5 : auto size = c.size();
296 5 : if (size <= 1)
297 0 : return 0.0;
298 :
299 5 : size--;
300 5 : R value = c[size] * size;
301 22 : for (std::size_t i = 1; i < size; ++i)
302 17 : value = value * x + c[size - i] * (size - i);
303 :
304 5 : return value;
305 : }
306 :
307 : template <typename T, typename T2>
308 : T
309 80 : clamp(const T & x, T2 lowerlimit, T2 upperlimit)
310 : {
311 80 : if (x < lowerlimit)
312 16 : return lowerlimit;
313 64 : if (x > upperlimit)
314 16 : return upperlimit;
315 48 : return x;
316 : }
317 :
318 : template <typename T, typename T2>
319 : T
320 160 : 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 160 : if (x <= start)
324 32 : return 0.0;
325 128 : else if (x >= end)
326 : {
327 32 : if (derivative)
328 16 : return 0.0;
329 : else
330 16 : return 1.0;
331 : }
332 96 : x = (x - start) / (end - start);
333 96 : if (derivative)
334 48 : return 30.0 * libMesh::Utility::pow<2>(x) * (x * (x - 2.0) + 1.0) / (end - start);
335 48 : 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 6 : 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 6 : if (x <= start)
350 2 : return 0.0;
351 4 : else if (x >= end)
352 : {
353 : if constexpr (compute_type == ComputeType::derivative)
354 1 : return 0.0;
355 : if constexpr (compute_type == ComputeType::value)
356 1 : return 1.0;
357 : }
358 2 : const auto u = (x - start) / (end - start);
359 : if constexpr (compute_type == ComputeType::derivative)
360 1 : return 30.0 * libMesh::Utility::pow<2>(u) * (u * (u - 2.0) + 1.0) / (end - start);
361 : if constexpr (compute_type == ComputeType::value)
362 1 : return libMesh::Utility::pow<3>(u) * (u * (u * 6.0 - 15.0) + 10.0);
363 : }
364 :
365 : /**
366 : * Helper function templates to set a variable to zero.
367 : * Specializations may have to be implemented (for examples see
368 : * RankTwoTensor, RankFourTensor).
369 : */
370 : template <typename T>
371 : inline void
372 58074 : mooseSetToZero(T & v)
373 : {
374 : /**
375 : * The default for non-pointer types is to assign zero.
376 : * This should either do something sensible, or throw a compiler error.
377 : * Otherwise the T type is designed badly.
378 : */
379 58074 : v = 0;
380 58074 : }
381 : template <typename T>
382 : inline void
383 : mooseSetToZero(T *&)
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 :
396 : /**
397 : * generate a complete multi index table for given dimension and order
398 : * i.e. given dim = 2, order = 2, generated table will have the following content
399 : * 0 0
400 : * 1 0
401 : * 0 1
402 : * 2 0
403 : * 1 1
404 : * 0 2
405 : * The first number in each entry represents the order of the first variable, i.e. x;
406 : * The second number in each entry represents the order of the second variable, i.e. y.
407 : * Multiplication is implied between numbers in each entry, i.e. 1 1 represents x^1 * y^1
408 : *
409 : * @param dim dimension of the multi-index, here dim = mesh dimension
410 : * @param order generate the multi-index up to certain order
411 : * @return a data structure holding entries representing the complete multi index
412 : */
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 2 : linearInterpolation(const X & x, const X1 & x1, const X2 & x2, const Y1 & y1, const Y2 & y2)
418 : {
419 2 : const auto m = (y2 - y1) / (x2 - x1);
420 : if constexpr (compute_type == ComputeType::derivative)
421 1 : return m;
422 : if constexpr (compute_type == ComputeType::value)
423 1 : return m * (x - x1) + y1;
424 : }
425 :
426 : /**
427 : * perform modulo operator for Euclidean division that ensures a non-negative result
428 : * @param dividend dividend of the modulo operation
429 : * @param divisor divisor of the modulo operation
430 : * @return the non-negative remainder when the dividend is divided by the divisor
431 : */
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 :
439 : /**
440 : * automatic prefixing for naming material properties based on gradients of coupled
441 : * variables/functors
442 : */
443 : template <typename T>
444 : T
445 1978 : gradName(const T & base_prop_name)
446 : {
447 1978 : return "grad_" + base_prop_name;
448 : }
449 :
450 : /**
451 : * automatic prefixing for naming material properties based on time derivatives of coupled
452 : * variables/functors
453 : */
454 : template <typename T>
455 : T
456 4206 : timeDerivName(const T & base_prop_name)
457 : {
458 4206 : return "d" + base_prop_name + "_dt";
459 : }
460 :
461 : /**
462 : * Computes the Kronecker product of two matrices.
463 : * @param product Reference to the product matrix
464 : * @param mat_A Reference to the first matrix
465 : * @param mat_B Reference to the other matrix
466 : */
467 : void kron(RealEigenMatrix & product, const RealEigenMatrix & mat_A, const RealEigenMatrix & mat_B);
468 :
469 : } // namespace MathUtils
470 :
471 : /// A helper function for MathUtils::multiIndex
472 : std::vector<std::vector<unsigned int>> multiIndexHelper(unsigned int N, unsigned int K);
|
