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 16 : round(T x)
78 : {
79 16 : return ::round(x); // use round from math.h
80 : }
81 :
82 : template <typename T>
83 : T
84 41877 : sign(T x)
85 : {
86 41877 : return x >= 0.0 ? 1.0 : -1.0;
87 : }
88 :
89 : template <typename T>
90 : T
91 615961759 : pow(T x, int e)
92 : {
93 615961759 : bool neg = false;
94 615961759 : T result = 1.0;
95 :
96 615961759 : if (e < 0)
97 : {
98 4 : neg = true;
99 4 : e = -e;
100 : }
101 :
102 1231442207 : while (e)
103 : {
104 : // if bit 0 is set multiply the current power of two factor of the exponent
105 615480448 : if (e & 1)
106 410401020 : result *= x;
107 :
108 : // x is incrementally set to consecutive powers of powers of two
109 615480448 : x *= x;
110 :
111 : // bit shift the exponent down
112 615480448 : e >>= 1;
113 : }
114 :
115 615962207 : return neg ? 1.0 / result : result;
116 448 : }
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(const T & x, Real smoothing_length)
128 : {
129 : if (x <= -smoothing_length)
130 : return 0.0;
131 : else if (x < smoothing_length)
132 : {
133 : using std::sin;
134 : return 0.5 * (1 + sin(libMesh::pi * x / 2 / smoothing_length));
135 : }
136 : else
137 : return 1.0;
138 : }
139 :
140 : template <typename T>
141 : T
142 : regularizedHeavysideDerivative(const T & x, Real smoothing_length)
143 : {
144 : if (x < smoothing_length && x > -smoothing_length)
145 : {
146 : using std::cos;
147 : return 0.25 * libMesh::pi / smoothing_length * (cos(libMesh::pi * x / 2 / smoothing_length));
148 : }
149 : else
150 : return 0.0;
151 : }
152 :
153 : template <typename T>
154 : T
155 : positivePart(T x)
156 : {
157 : return x > 0.0 ? x : 0.0;
158 : }
159 :
160 : template <typename T>
161 : T
162 : negativePart(T x)
163 : {
164 : return x < 0.0 ? x : 0.0;
165 : }
166 :
167 : template <
168 : typename T,
169 : typename T2,
170 : typename T3,
171 : typename std::enable_if<libMesh::ScalarTraits<T>::value && libMesh::ScalarTraits<T2>::value &&
172 : libMesh::ScalarTraits<T3>::value,
173 : int>::type = 0>
174 : void
175 5646488 : addScaled(const T & a, const T2 & b, T3 & result)
176 : {
177 5646488 : result += a * b;
178 5646488 : }
179 :
180 : template <typename T,
181 : typename T2,
182 : typename T3,
183 : typename std::enable_if<libMesh::ScalarTraits<T>::value, int>::type = 0>
184 : void
185 62316 : addScaled(const T & scalar,
186 : const libMesh::NumericVector<T2> & numeric_vector,
187 : libMesh::NumericVector<T3> & result)
188 : {
189 62316 : result.add(scalar, numeric_vector);
190 62316 : }
191 :
192 : template <
193 : typename T,
194 : typename T2,
195 : template <typename>
196 : class W,
197 : template <typename>
198 : class W2,
199 : typename std::enable_if<std::is_same<typename W<T>::index_type, unsigned int>::value &&
200 : std::is_same<typename W2<T2>::index_type, unsigned int>::value,
201 : int>::type = 0>
202 : typename libMesh::CompareTypes<T, T2>::supertype
203 133584148 : dotProduct(const W<T> & a, const W2<T2> & b)
204 : {
205 133584148 : return a * b;
206 : }
207 :
208 : template <typename T,
209 : typename T2,
210 : template <typename>
211 : class W,
212 : template <typename>
213 : class W2,
214 : typename std::enable_if<std::is_same<typename W<T>::index_type,
215 : std::tuple<unsigned int, unsigned int>>::value &&
216 : std::is_same<typename W2<T2>::index_type,
217 : std::tuple<unsigned int, unsigned int>>::value,
218 : int>::type = 0>
219 : typename libMesh::CompareTypes<T, T2>::supertype
220 0 : dotProduct(const W<T> & a, const W2<T2> & b)
221 : {
222 0 : return a.contract(b);
223 : }
224 :
225 : /**
226 : * Evaluate a polynomial with the coefficients c at x. Note that the Polynomial
227 : * form is
228 : * 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]
229 : * where s = c.size()-1 , which is counter intuitive!
230 : *
231 : * This function will be DEPRECATED soon (10/22/2020)
232 : *
233 : * The coefficient container type can be any container that provides an index
234 : * operator [] and a .size() method (e.g. std::vector, std::array). The return
235 : * type is the supertype of the container value type and the argument x.
236 : * The supertype is the type that can represent both number types.
237 : */
238 : template <typename C,
239 : typename T,
240 : typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
241 : R
242 368 : poly(const C & c, const T x, const bool derivative = false)
243 : {
244 368 : const auto size = c.size();
245 368 : if (size == 0)
246 0 : return 0.0;
247 :
248 368 : R value = c[0];
249 368 : if (derivative)
250 : {
251 94 : value *= size - 1;
252 916 : for (std::size_t i = 1; i < size - 1; ++i)
253 822 : value = value * x + c[i] * (size - i - 1);
254 : }
255 : else
256 : {
257 1370 : for (std::size_t i = 1; i < size; ++i)
258 1096 : value = value * x + c[i];
259 : }
260 :
261 368 : return value;
262 : }
263 :
264 : /**
265 : * Evaluate a polynomial with the coefficients c at x. Note that the Polynomial
266 : * form is
267 : * c[0] + c[1] * x + c[2] * x^2 + ...
268 : * The coefficient container type can be any container that provides an index
269 : * operator [] and a .size() method (e.g. std::vector, std::array). The return
270 : * type is the supertype of the container value type and the argument x.
271 : * The supertype is the type that can represent both number types.
272 : */
273 : template <typename C,
274 : typename T,
275 : typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
276 : R
277 10 : polynomial(const C & c, const T x)
278 : {
279 10 : auto size = c.size();
280 10 : if (size == 0)
281 0 : return 0.0;
282 :
283 10 : size--;
284 10 : R value = c[size];
285 54 : for (std::size_t i = 1; i <= size; ++i)
286 44 : value = value * x + c[size - i];
287 :
288 10 : return value;
289 : }
290 :
291 : /**
292 : * Returns the derivative of polynomial(c, x) with respect to x
293 : */
294 : template <typename C,
295 : typename T,
296 : typename R = typename libMesh::CompareTypes<typename C::value_type, T>::supertype>
297 : R
298 10 : polynomialDerivative(const C & c, const T x)
299 : {
300 10 : auto size = c.size();
301 10 : if (size <= 1)
302 0 : return 0.0;
303 :
304 10 : size--;
305 10 : R value = c[size] * size;
306 44 : for (std::size_t i = 1; i < size; ++i)
307 34 : value = value * x + c[size - i] * (size - i);
308 :
309 10 : return value;
310 : }
311 :
312 : template <typename T, typename T2>
313 : T
314 90 : clamp(const T & x, T2 lowerlimit, T2 upperlimit)
315 : {
316 90 : if (x < lowerlimit)
317 18 : return lowerlimit;
318 72 : if (x > upperlimit)
319 18 : return upperlimit;
320 54 : return x;
321 : }
322 :
323 : template <typename T, typename T2>
324 : T
325 180 : smootherStep(T x, T2 start, T2 end, bool derivative = false)
326 : {
327 : mooseAssert("start < end", "Start value must be lower than end value for smootherStep");
328 180 : if (x <= start)
329 36 : return 0.0;
330 144 : else if (x >= end)
331 : {
332 36 : if (derivative)
333 18 : return 0.0;
334 : else
335 18 : return 1.0;
336 : }
337 108 : x = (x - start) / (end - start);
338 108 : if (derivative)
339 54 : return 30.0 * libMesh::Utility::pow<2>(x) * (x * (x - 2.0) + 1.0) / (end - start);
340 54 : return libMesh::Utility::pow<3>(x) * (x * (x * 6.0 - 15.0) + 10.0);
341 : }
342 :
343 : enum class ComputeType
344 : {
345 : value,
346 : derivative
347 : };
348 :
349 : template <ComputeType compute_type, typename X, typename S, typename E>
350 : auto
351 12 : smootherStep(const X & x, const S & start, const E & end)
352 : {
353 : mooseAssert("start < end", "Start value must be lower than end value for smootherStep");
354 12 : if (x <= start)
355 4 : return 0.0;
356 8 : else if (x >= end)
357 : {
358 : if constexpr (compute_type == ComputeType::derivative)
359 2 : return 0.0;
360 : if constexpr (compute_type == ComputeType::value)
361 2 : return 1.0;
362 : }
363 4 : const auto u = (x - start) / (end - start);
364 : if constexpr (compute_type == ComputeType::derivative)
365 2 : return 30.0 * libMesh::Utility::pow<2>(u) * (u * (u - 2.0) + 1.0) / (end - start);
366 : if constexpr (compute_type == ComputeType::value)
367 2 : return libMesh::Utility::pow<3>(u) * (u * (u * 6.0 - 15.0) + 10.0);
368 : }
369 :
370 : /**
371 : * Helper function templates to set a variable to zero.
372 : * Specializations may have to be implemented (for examples see
373 : * RankTwoTensor, RankFourTensor).
374 : */
375 : template <typename T>
376 : inline void
377 56993 : mooseSetToZero(T & v)
378 : {
379 : /**
380 : * The default for non-pointer types is to assign zero.
381 : * This should either do something sensible, or throw a compiler error.
382 : * Otherwise the T type is designed badly.
383 : */
384 56993 : v = 0;
385 56993 : }
386 : template <typename T>
387 : inline void
388 : mooseSetToZero(T *&)
389 : {
390 : mooseError("mooseSetToZero does not accept pointers");
391 : }
392 :
393 : template <>
394 : inline void
395 : mooseSetToZero(std::vector<Real> & vec)
396 : {
397 : for (auto & v : vec)
398 : v = 0.;
399 : }
400 :
401 : /**
402 : * generate a complete multi index table for given dimension and order
403 : * i.e. given dim = 2, order = 2, generated table will have the following content
404 : * 0 0
405 : * 1 0
406 : * 0 1
407 : * 2 0
408 : * 1 1
409 : * 0 2
410 : * The first number in each entry represents the order of the first variable, i.e. x;
411 : * The second number in each entry represents the order of the second variable, i.e. y.
412 : * Multiplication is implied between numbers in each entry, i.e. 1 1 represents x^1 * y^1
413 : *
414 : * @param dim dimension of the multi-index, here dim = mesh dimension
415 : * @param order generate the multi-index up to certain order
416 : * @return a data structure holding entries representing the complete multi index
417 : */
418 : std::vector<std::vector<unsigned int>> multiIndex(unsigned int dim, unsigned int order);
419 :
420 : template <ComputeType compute_type, typename X, typename X1, typename X2, typename Y1, typename Y2>
421 : auto
422 4 : linearInterpolation(const X & x, const X1 & x1, const X2 & x2, const Y1 & y1, const Y2 & y2)
423 : {
424 4 : const auto m = (y2 - y1) / (x2 - x1);
425 : if constexpr (compute_type == ComputeType::derivative)
426 2 : return m;
427 : if constexpr (compute_type == ComputeType::value)
428 2 : return m * (x - x1) + y1;
429 : }
430 :
431 : /**
432 : * perform modulo operator for Euclidean division that ensures a non-negative result
433 : * @param dividend dividend of the modulo operation
434 : * @param divisor divisor of the modulo operation
435 : * @return the non-negative remainder when the dividend is divided by the divisor
436 : */
437 : template <typename T1, typename T2>
438 : std::size_t
439 10560 : euclideanMod(T1 dividend, T2 divisor)
440 : {
441 10560 : return (dividend % divisor + divisor) % divisor;
442 : }
443 :
444 : /**
445 : * automatic prefixing for naming material properties based on gradients of coupled
446 : * variables/functors
447 : */
448 : template <typename T>
449 : T
450 2083 : gradName(const T & base_prop_name)
451 : {
452 2083 : return "grad_" + base_prop_name;
453 : }
454 :
455 : /**
456 : * automatic prefixing for naming material properties based on time derivatives of coupled
457 : * variables/functors
458 : */
459 : template <typename T>
460 : T
461 4429 : timeDerivName(const T & base_prop_name)
462 : {
463 4429 : return "d" + base_prop_name + "_dt";
464 : }
465 :
466 : /**
467 : * Computes the Kronecker product of two matrices.
468 : * @param product Reference to the product matrix
469 : * @param mat_A Reference to the first matrix
470 : * @param mat_B Reference to the other matrix
471 : */
472 : void kron(RealEigenMatrix & product, const RealEigenMatrix & mat_A, const RealEigenMatrix & mat_B);
473 :
474 : } // namespace MathUtils
475 :
476 : /// A helper function for MathUtils::multiIndex
477 : std::vector<std::vector<unsigned int>> multiIndexHelper(unsigned int N, unsigned int K);
|
