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 : #include "MathUtils.h"
11 : #include "MooseUtils.h"
12 :
13 : namespace MathUtils
14 : {
15 :
16 : void
17 0 : kron(RealEigenMatrix & product, const RealEigenMatrix & mat_A, const RealEigenMatrix & mat_B)
18 : {
19 0 : product.resize(mat_A.rows() * mat_B.rows(), mat_A.cols() * mat_B.cols());
20 0 : for (unsigned int i = 0; i < mat_A.rows(); i++)
21 0 : for (unsigned int j = 0; j < mat_A.cols(); j++)
22 0 : for (unsigned int k = 0; k < mat_B.rows(); k++)
23 0 : for (unsigned int l = 0; l < mat_B.cols(); l++)
24 0 : product(((i * mat_B.rows()) + k), ((j * mat_B.cols()) + l)) = mat_A(i, j) * mat_B(k, l);
25 0 : }
26 :
27 : Real
28 16 : plainLog(Real x, unsigned int derivative_order)
29 : {
30 16 : switch (derivative_order)
31 : {
32 4 : case 0:
33 4 : return std::log(x);
34 :
35 4 : case 1:
36 4 : return 1.0 / x;
37 :
38 4 : case 2:
39 4 : return -1.0 / (x * x);
40 :
41 4 : case 3:
42 4 : return 2.0 / (x * x * x);
43 : }
44 :
45 0 : mooseError("Unsupported derivative order ", derivative_order);
46 : }
47 :
48 : Real
49 8 : poly1Log(Real x, Real tol, unsigned int derivative_order)
50 : {
51 8 : if (x >= tol)
52 4 : return plainLog(x, derivative_order);
53 :
54 2 : const auto c1 = [&]() { return 1.0 / tol; };
55 1 : const auto c2 = [&]() { return std::log(tol) - 1.0; };
56 :
57 4 : switch (derivative_order)
58 : {
59 1 : case 0:
60 1 : return c1() * x + c2();
61 :
62 1 : case 1:
63 1 : return c1();
64 :
65 1 : case 2:
66 1 : return 0.0;
67 :
68 1 : case 3:
69 1 : return 0.0;
70 : }
71 :
72 0 : mooseError("Unsupported derivative order ", derivative_order);
73 : }
74 :
75 : Real
76 8 : poly2Log(Real x, Real tol, unsigned int derivative_order)
77 : {
78 8 : if (x >= tol)
79 4 : return plainLog(x, derivative_order);
80 :
81 3 : const auto c1 = [&]() { return -0.5 / (tol * tol); };
82 2 : const auto c2 = [&]() { return 2.0 / tol; };
83 1 : const auto c3 = [&]() { return std::log(tol) - 3.0 / 2.0; };
84 :
85 4 : switch (derivative_order)
86 : {
87 1 : case 0:
88 1 : return c1() * x * x + c2() * x + c3();
89 :
90 1 : case 1:
91 1 : return 2.0 * c1() * x + c2();
92 :
93 1 : case 2:
94 1 : return 2.0 * c1();
95 :
96 1 : case 3:
97 1 : return 0.0;
98 : }
99 :
100 0 : mooseError("Unsupported derivative order ", derivative_order);
101 : }
102 :
103 : Real
104 8 : poly3Log(Real x, Real tol, unsigned int derivative_order)
105 : {
106 8 : if (x >= tol)
107 4 : return plainLog(x, derivative_order);
108 :
109 4 : const auto c1 = [&]() { return 1.0 / (3.0 * tol * tol * tol); };
110 3 : const auto c2 = [&]() { return -3.0 / (2.0 * tol * tol); };
111 2 : const auto c3 = [&]() { return 3.0 / tol; };
112 1 : const auto c4 = [&]() { return std::log(tol) - 11.0 / 6.0; };
113 :
114 4 : switch (derivative_order)
115 : {
116 1 : case 0:
117 1 : return c1() * x * x * x + c2() * x * x + c3() * x + c4();
118 :
119 1 : case 1:
120 1 : return 3.0 * c1() * x * x + 2.0 * c2() * x + c3();
121 :
122 1 : case 2:
123 1 : return 6.0 * c1() * x + 2.0 * c2();
124 :
125 1 : case 3:
126 1 : return 6.0 * c1();
127 : }
128 :
129 0 : mooseError("Unsupported derivative order ", derivative_order);
130 : }
131 :
132 : Real
133 8 : poly4Log(Real x, Real tol, unsigned int derivative_order)
134 : {
135 8 : if (x >= tol)
136 4 : return plainLog(x, derivative_order);
137 :
138 4 : switch (derivative_order)
139 : {
140 1 : case 0:
141 2 : return std::log(tol) + (x - tol) / tol -
142 1 : Utility::pow<2>(x - tol) / (2.0 * Utility::pow<2>(tol)) +
143 1 : Utility::pow<3>(x - tol) / (3.0 * Utility::pow<3>(tol)) -
144 1 : Utility::pow<4>(x - tol) / (4.0 * Utility::pow<4>(tol)) +
145 1 : Utility::pow<5>(x - tol) / (5.0 * Utility::pow<5>(tol)) -
146 1 : Utility::pow<6>(x - tol) / (6.0 * Utility::pow<6>(tol));
147 :
148 1 : case 1:
149 1 : return 1.0 / tol - (x - tol) / Utility::pow<2>(tol) +
150 1 : Utility::pow<2>(x - tol) / Utility::pow<3>(tol) -
151 1 : Utility::pow<3>(x - tol) / Utility::pow<4>(tol) +
152 1 : Utility::pow<4>(x - tol) / Utility::pow<5>(tol) -
153 1 : Utility::pow<5>(x - tol) / Utility::pow<6>(tol);
154 :
155 1 : case 2:
156 1 : return -1.0 / Utility::pow<2>(tol) + 2.0 * (x - tol) / Utility::pow<3>(tol) -
157 1 : 3.0 * Utility::pow<2>(x - tol) / Utility::pow<4>(tol) +
158 1 : 4.0 * Utility::pow<3>(x - tol) / Utility::pow<5>(tol) -
159 1 : 5.0 * Utility::pow<4>(x - tol) / Utility::pow<6>(tol);
160 :
161 1 : case 3:
162 1 : return 2.0 / Utility::pow<3>(tol) - 6.0 * (x - tol) / Utility::pow<4>(tol) +
163 1 : 12.0 * Utility::pow<2>(x - tol) / Utility::pow<5>(tol) -
164 1 : 20.0 * Utility::pow<3>(x - tol) / Utility::pow<6>(tol);
165 : }
166 :
167 0 : mooseError("Unsupported derivative order ", derivative_order);
168 : }
169 :
170 : /// \todo This can be done without std::pow!
171 : Real
172 8 : taylorLog(Real x)
173 : {
174 8 : Real y = (x - 1.0) / (x + 1.0);
175 8 : Real val = 1.0;
176 48 : for (unsigned int i = 0; i < 5; ++i)
177 : {
178 40 : Real exponent = i + 2.0;
179 40 : val += 1.0 / (2.0 * (i + 1.0) + 1.0) * std::pow(y, exponent);
180 : }
181 :
182 8 : return val * 2.0 * y;
183 : }
184 :
185 : std::vector<std::vector<unsigned int>>
186 68 : multiIndex(unsigned int dim, unsigned int order)
187 : {
188 : // first row all zero
189 68 : std::vector<std::vector<unsigned int>> multi_index;
190 68 : std::vector<std::vector<unsigned int>> n_choose_k;
191 68 : std::vector<unsigned int> row(dim, 0);
192 68 : multi_index.push_back(row);
193 :
194 68 : if (dim == 1)
195 4 : for (unsigned int q = 1; q <= order; q++)
196 : {
197 3 : row[0] = q;
198 3 : multi_index.push_back(row);
199 : }
200 : else
201 137 : for (unsigned int q = 1; q <= order; q++)
202 : {
203 70 : n_choose_k = multiIndexHelper(dim + q - 1, dim - 1);
204 218 : for (unsigned int r = 0; r < n_choose_k.size(); r++)
205 : {
206 148 : row.clear();
207 453 : for (unsigned int c = 1; c < n_choose_k[0].size(); c++)
208 305 : row.push_back(n_choose_k[r][c] - n_choose_k[r][c - 1] - 1);
209 148 : multi_index.push_back(row);
210 : }
211 : }
212 :
213 136 : return multi_index;
214 68 : }
215 :
216 : Point
217 0 : barycentricToCartesian2D(const Point & p0,
218 : const Point & p1,
219 : const Point & p2,
220 : const Real b0,
221 : const Real b1,
222 : const Real b2)
223 : {
224 : mooseAssert(!MooseUtils::isZero(b0 + b1 + b2 - 1.0), "Barycentric coordinates must sum to one!");
225 :
226 0 : Point center;
227 :
228 0 : for (unsigned int d = 0; d < 2; ++d)
229 0 : center(d) = p0(d) * b0 + p1(d) * b1 + p2(d) * b2;
230 : // p0, p1, p2 are vertices of triangle
231 : // b0, b1, b2 are Barycentric coordinates of the triangle center
232 :
233 0 : return center;
234 : }
235 :
236 : Point
237 0 : barycentricToCartesian3D(const Point & p0,
238 : const Point & p1,
239 : const Point & p2,
240 : const Point & p3,
241 : const Real b0,
242 : const Real b1,
243 : const Real b2,
244 : const Real b3)
245 : {
246 : mooseAssert(!MooseUtils::isZero(b0 + b1 + b2 + b3 - 1.0),
247 : "Barycentric coordinates must sum to one!");
248 :
249 0 : Point center;
250 :
251 0 : for (unsigned int d = 0; d < 3; ++d)
252 0 : center(d) = p0(d) * b0 + p1(d) * b1 + p2(d) * b2 + p3(d) * b3;
253 : // p0, p1, p2, p3 are vertices of tetrahedron
254 : // b0, b1, b2, b3 are Barycentric coordinates of the tetrahedron center
255 :
256 0 : return center;
257 : }
258 :
259 : Point
260 0 : circumcenter2D(const Point & p0, const Point & p1, const Point & p2)
261 : {
262 : // Square of triangle edge lengths
263 0 : Real edge01 = (p0 - p1).norm_sq();
264 0 : Real edge02 = (p0 - p2).norm_sq();
265 0 : Real edge12 = (p1 - p2).norm_sq();
266 :
267 : // Barycentric weights for circumcenter
268 0 : Real weight0 = edge12 * (edge01 + edge02 - edge12);
269 0 : Real weight1 = edge02 * (edge01 + edge12 - edge02);
270 0 : Real weight2 = edge01 * (edge02 + edge12 - edge01);
271 :
272 0 : Real sum_weights = weight0 + weight1 + weight2;
273 :
274 : // Check to make sure vertices are not collinear
275 0 : if (MooseUtils::isZero(sum_weights))
276 0 : mooseError("Cannot evaluate circumcenter. Points should be non-collinear.");
277 :
278 0 : Real inv_sum_weights = 1.0 / sum_weights;
279 :
280 : // Barycentric coordinates
281 0 : Real b0 = weight0 * inv_sum_weights;
282 0 : Real b1 = weight1 * inv_sum_weights;
283 0 : Real b2 = weight2 * inv_sum_weights;
284 :
285 0 : return MathUtils::barycentricToCartesian2D(p0, p1, p2, b0, b1, b2);
286 : }
287 :
288 : Point
289 0 : circumcenter3D(const Point & p0, const Point & p1, const Point & p2, const Point & p3)
290 : {
291 : // Square of tetrahedron edge lengths
292 0 : Real edge01 = (p0 - p1).norm_sq();
293 0 : Real edge02 = (p0 - p2).norm_sq();
294 0 : Real edge03 = (p0 - p3).norm_sq();
295 0 : Real edge12 = (p1 - p2).norm_sq();
296 0 : Real edge13 = (p1 - p3).norm_sq();
297 0 : Real edge23 = (p2 - p3).norm_sq();
298 :
299 : // Barycentric weights for circumcenter
300 0 : Real weight0 = -2 * edge12 * edge13 * edge23 + edge01 * edge23 * (edge13 + edge12 - edge23) +
301 0 : edge02 * edge13 * (edge12 + edge23 - edge13) +
302 0 : edge03 * edge12 * (edge13 + edge23 - edge12);
303 0 : Real weight1 = -2 * edge02 * edge03 * edge23 + edge01 * edge23 * (edge02 + edge03 - edge23) +
304 0 : edge13 * edge02 * (edge03 + edge23 - edge02) +
305 0 : edge12 * edge03 * (edge02 + edge23 - edge03);
306 0 : Real weight2 = -2 * edge01 * edge03 * edge13 + edge02 * edge13 * (edge01 + edge03 - edge13) +
307 0 : edge23 * edge01 * (edge03 + edge13 - edge01) +
308 0 : edge12 * edge03 * (edge01 + edge13 - edge03);
309 0 : Real weight3 = -2 * edge01 * edge02 * edge12 + edge03 * edge12 * (edge01 + edge02 - edge12) +
310 0 : edge23 * edge01 * (edge02 + edge12 - edge01) +
311 0 : edge13 * edge02 * (edge01 + edge12 - edge02);
312 :
313 0 : Real sum_weights = weight0 + weight1 + weight2 + weight3;
314 :
315 : // Check to make sure vertices are not coplanar
316 0 : if (MooseUtils::isZero(sum_weights))
317 0 : mooseError("Cannot evaluate circumcenter. Points should be non-coplanar.");
318 :
319 0 : Real inv_sum_weights = 1.0 / sum_weights;
320 :
321 : // Barycentric coordinates
322 0 : Real b0 = weight0 * inv_sum_weights;
323 0 : Real b1 = weight1 * inv_sum_weights;
324 0 : Real b2 = weight2 * inv_sum_weights;
325 0 : Real b3 = weight3 * inv_sum_weights;
326 :
327 0 : return MathUtils::barycentricToCartesian3D(p0, p1, p2, p3, b0, b1, b2, b3);
328 : }
329 :
330 : } // namespace MathUtils
331 :
332 : std::vector<std::vector<unsigned int>>
333 70 : multiIndexHelper(unsigned int N, unsigned int K)
334 : {
335 70 : std::vector<std::vector<unsigned int>> n_choose_k;
336 70 : std::vector<unsigned int> row;
337 70 : std::string bitmask(K, 1); // K leading 1's
338 70 : bitmask.resize(N, 0); // N-K trailing 0's
339 :
340 : do
341 : {
342 148 : row.clear();
343 148 : row.push_back(0);
344 470 : for (unsigned int i = 0; i < N; ++i) // [0..N-1] integers
345 322 : if (bitmask[i])
346 157 : row.push_back(i + 1);
347 148 : row.push_back(N + 1);
348 148 : n_choose_k.push_back(row);
349 148 : } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
350 :
351 140 : return n_choose_k;
352 70 : }
|
