v0.13.2
Loading...
Searching...
No Matches
BernsteinBezier.cpp
Go to the documentation of this file.
1/** \file BernsteinBezier.cpp
2\brief Bernstein-Bezier polynomials for H1 space
3
4Implementation based on \cite ainsworth2011bernstein
5
6*/
7
8template <int D, int Side>
9MoFEMErrorCode BernsteinBezier::generateIndicesVertex(const int N, int *alpha) {
10 MoFEMFunctionBeginHot;
11 alpha[(D + 1) * Side + Side] = N;
12 MoFEMFunctionReturnHot(0);
13}
14
15template <int D, int Side>
16MoFEMErrorCode BernsteinBezier::generateIndicesEdgeOnSimplex(const int N,
17 int *alpha) {
18 constexpr int edge_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
19 {0, 3}, {1, 3}, {2, 3}};
20 MoFEMFunctionBeginHot;
21 std::fill(alpha, &alpha[(D + 1) * NBEDGE_H1(N)], 0);
22 for (int n = N - 1; n != 0; --n) {
23
24 // + o o o o +
25
26 alpha[edge_nodes[Side][0]] = n;
27 alpha[edge_nodes[Side][1]] = N - n;
28 alpha += D + 1;
29 }
30 MoFEMFunctionReturnHot(0);
31}
32
33template <int D, int Side>
34MoFEMErrorCode BernsteinBezier::generateIndicesTriOnSimplex(const int N,
35 int *alpha) {
36 constexpr int tri_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
37 MoFEMFunctionBeginHot;
38 std::fill(alpha, &alpha[(D + 1) * NBFACETRI_H1(N)], 0);
39 for (int n0 = 1; n0 != N - 1; ++n0) {
40 for (int n1 = 1; n1 != N - n0; ++n1) {
41
42 // +
43 // + +
44 // + o +
45 // + o o +
46 // + o o o +
47 // + + + + + +
48
49 alpha[tri_nodes[Side][0]] = n0;
50 alpha[tri_nodes[Side][1]] = n1;
51 alpha[tri_nodes[Side][2]] = N - n0 - n1;
52
53 alpha += D + 1;
54 }
55 }
56 MoFEMFunctionReturnHot(0);
57}
58
59MoFEMErrorCode BernsteinBezier::generateIndicesTetOnSimplex(const int N,
60 int *alpha) {
61 MoFEMFunctionBeginHot;
62 std::fill(alpha, &alpha[4 * NBVOLUMETET_H1(N)], 0);
63 for (int n0 = 1; n0 != N - 1; ++n0) {
64 for (int n1 = 1; n1 != N - n0; ++n1) {
65 for (int n2 = 1; n2 != N - n0 - n1; ++n2) {
66 alpha[0] = n0;
67 alpha[1] = n1;
68 alpha[2] = n2;
69 alpha[3] = N - n0 - n1 - n2;
70 alpha += 4;
71 }
72 }
73 }
74 MoFEMFunctionReturnHot(0);
75}
76
77MoFEMErrorCode BernsteinBezier::generateIndicesVertexEdge(const int N,
78 int *alpha) {
79 MoFEMFunctionBeginHot;
80 CHKERR generateIndicesVertex<1, 0>(N, alpha);
81 CHKERR generateIndicesVertex<1, 1>(N, alpha);
82 MoFEMFunctionReturnHot(0);
83}
84
85MoFEMErrorCode BernsteinBezier::generateIndicesEdgeEdge(const int N,
86 int *alpha) {
87 return generateIndicesEdgeOnSimplex<1, 0>(N, alpha);
88}
89
90MoFEMErrorCode BernsteinBezier::generateIndicesVertexTri(const int N,
91 int *alpha) {
92 MoFEMFunctionBeginHot;
93 CHKERR generateIndicesVertex<2, 0>(N, alpha);
94 CHKERR generateIndicesVertex<2, 1>(N, alpha);
95 CHKERR generateIndicesVertex<2, 2>(N, alpha);
96 MoFEMFunctionReturnHot(0);
97}
98
99MoFEMErrorCode BernsteinBezier::generateIndicesEdgeTri(const int N[],
100 int *alpha[]) {
101 MoFEMFunctionBeginHot;
102 CHKERR generateIndicesEdgeOnSimplex<2, 0>(N[0], alpha[0]);
103 CHKERR generateIndicesEdgeOnSimplex<2, 1>(N[1], alpha[1]);
104 CHKERR generateIndicesEdgeOnSimplex<2, 2>(N[2], alpha[2]);
105 MoFEMFunctionReturnHot(0);
106}
107
108MoFEMErrorCode BernsteinBezier::generateIndicesEdgeTri(const int side,
109 const int N,
110 int *alpha) {
111 switch (side) {
112 case 0:
113 return generateIndicesEdgeOnSimplex<2, 0>(N, alpha);
114 case 1:
115 return generateIndicesEdgeOnSimplex<2, 1>(N, alpha);
116 case 2:
117 return generateIndicesEdgeOnSimplex<2, 2>(N, alpha);
118 default:
119 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
120 "Wrong side number, on triangle are edges from 0 to 2");
121 }
122}
123
124MoFEMErrorCode BernsteinBezier::generateIndicesTriTri(const int N, int *alpha) {
125 return generateIndicesTriOnSimplex<2, 3>(N, alpha);
126}
127
128MoFEMErrorCode BernsteinBezier::generateIndicesVertexTet(const int N,
129 int *alpha) {
130 MoFEMFunctionBeginHot;
131 CHKERR generateIndicesVertex<3, 0>(N, alpha);
132 CHKERR generateIndicesVertex<3, 1>(N, alpha);
133 CHKERR generateIndicesVertex<3, 2>(N, alpha);
134 CHKERR generateIndicesVertex<3, 3>(N, alpha);
135 MoFEMFunctionReturnHot(0);
136}
137
138MoFEMErrorCode BernsteinBezier::generateIndicesEdgeTet(const int N[],
139 int *alpha[]) {
140 MoFEMFunctionBeginHot;
141 CHKERR generateIndicesEdgeOnSimplex<3, 0>(N[0], alpha[0]);
142 CHKERR generateIndicesEdgeOnSimplex<3, 1>(N[1], alpha[1]);
143 CHKERR generateIndicesEdgeOnSimplex<3, 2>(N[2], alpha[2]);
144 CHKERR generateIndicesEdgeOnSimplex<3, 3>(N[3], alpha[3]);
145 CHKERR generateIndicesEdgeOnSimplex<3, 4>(N[4], alpha[4]);
146 CHKERR generateIndicesEdgeOnSimplex<3, 5>(N[5], alpha[5]);
147 MoFEMFunctionReturnHot(0);
148}
149
150MoFEMErrorCode BernsteinBezier::generateIndicesEdgeTet(const int side,
151 const int N,
152 int *alpha) {
153 MoFEMFunctionBeginHot;
154 switch (side) {
155 case 0:
156 return generateIndicesEdgeOnSimplex<3, 0>(N, alpha);
157 case 1:
158 return generateIndicesEdgeOnSimplex<3, 1>(N, alpha);
159 case 2:
160 return generateIndicesEdgeOnSimplex<3, 2>(N, alpha);
161 case 3:
162 return generateIndicesEdgeOnSimplex<3, 3>(N, alpha);
163 case 4:
164 return generateIndicesEdgeOnSimplex<3, 4>(N, alpha);
165 case 5:
166 return generateIndicesEdgeOnSimplex<3, 5>(N, alpha);
167 default:
168 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
169 "Wrong side number, on tetrahedron are edges from 0 to 5");
170 }
171 MoFEMFunctionReturnHot(0);
172}
173
174MoFEMErrorCode BernsteinBezier::generateIndicesTriTet(const int N[],
175 int *alpha[]) {
176 MoFEMFunctionBeginHot;
177 CHKERR generateIndicesTriOnSimplex<3, 0>(N[0], alpha[0]);
178 CHKERR generateIndicesTriOnSimplex<3, 1>(N[1], alpha[1]);
179 CHKERR generateIndicesTriOnSimplex<3, 2>(N[2], alpha[2]);
180 CHKERR generateIndicesTriOnSimplex<3, 3>(N[3], alpha[3]);
181 MoFEMFunctionReturnHot(0);
182}
183
184MoFEMErrorCode BernsteinBezier::generateIndicesTriTet(const int side,
185 const int N, int *alpha) {
186 MoFEMFunctionBeginHot;
187 switch (side) {
188 case 0:
189 return generateIndicesTriOnSimplex<3, 0>(N, alpha);
190 case 1:
191 return generateIndicesTriOnSimplex<3, 1>(N, alpha);
192 case 2:
193 return generateIndicesTriOnSimplex<3, 2>(N, alpha);
194 case 3:
195 return generateIndicesTriOnSimplex<3, 3>(N, alpha);
196 default:
197 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
198 "Wrong side number, on tetrahedron are from from 0 to 3");
199 }
200 MoFEMFunctionReturnHot(0);
201}
202
203MoFEMErrorCode BernsteinBezier::generateIndicesTetTet(const int N, int *alpha) {
204 return generateIndicesTetOnSimplex(N, alpha);
205}
206
207template <int D>
208MoFEMErrorCode BernsteinBezier::genrateDerivativeIndices(
209 const int N, const int n_alpha, const int *alpha, const int *diff,
210 const int n_alpha_diff, const int *alpha_diff, double *c) {
211 MoFEMFunctionBeginHot;
212
213 std::fill(c, &c[n_alpha * n_alpha_diff], 0);
214 int abs_diff = diff[0];
215 for (int d = 1; d != D + 1; ++d)
216 abs_diff += diff[d];
217
218 const double b = boost::math::factorial<double>(N) /
219 boost::math::factorial<double>(N - abs_diff);
220
221 for (int i = 0; i != n_alpha; ++i) {
222 for (int j = 0; j != n_alpha_diff; ++j) {
223 int d = 0;
224 for (; d != D + 1; ++d)
225 if (alpha_diff[(D + 1) * j + d] + diff[d] != alpha[(D + 1) * i + d])
226 break;
227 if (d == D + 1) {
228 c[i * n_alpha_diff + j] = b;
229 break;
230 }
231 }
232 }
233
234 MoFEMFunctionReturnHot(0);
235}
236
237MoFEMErrorCode BernsteinBezier::genrateDerivativeIndicesEdges(
238 const int N, const int n_alpha, const int *alpha, const int *diff,
239 const int n_alpha_diff, const int *alpha_diff, double *c) {
240 return genrateDerivativeIndices<1>(N, n_alpha, alpha, diff, n_alpha_diff,
241 alpha_diff, c);
242}
243
244MoFEMErrorCode BernsteinBezier::genrateDerivativeIndicesTri(
245 const int N, const int n_alpha, const int *alpha, const int *diff,
246 const int n_alpha_diff, const int *alpha_diff, double *c) {
247 return genrateDerivativeIndices<2>(N, n_alpha, alpha, diff, n_alpha_diff,
248 alpha_diff, c);
249}
250
251MoFEMErrorCode BernsteinBezier::genrateDerivativeIndicesTet(
252 const int N, const int n_alpha, const int *alpha, const int *diff,
253 const int n_alpha_diff, const int *alpha_diff, double *c) {
254 return genrateDerivativeIndices<3>(N, n_alpha, alpha, diff, n_alpha_diff,
255 alpha_diff, c);
256}
257
258template <int D>
259MoFEMErrorCode
260BernsteinBezier::domainPoints(const int N, const int n_x, const int n_alpha,
261 const int *alpha, const double *x_k,
262 double *x_alpha) {
263 FTensor::Index<'i', D> i;
264 MoFEMFunctionBeginHot;
265
266 auto t_x_alpha = FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
267 x_alpha, &x_alpha[1], &x_alpha[2]);
268 for (int n0 = 0; n0 != n_alpha; ++n0) {
269 t_x_alpha(i) = 0;
270 auto t_x_k = FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3>(
271 x_k, &x_k[1], &x_k[2]);
272 for (int n1 = 0; n1 != n_x; ++n1) {
273 t_x_alpha(i) += static_cast<double>(*alpha) * t_x_k(i);
274 ++t_x_k;
275 ++alpha;
276 }
277 t_x_alpha(i) /= static_cast<double>(N);
278 ++t_x_alpha;
279 }
280
281 MoFEMFunctionReturnHot(0);
282}
283
284MoFEMErrorCode BernsteinBezier::domainPoints3d(const int N, const int n_x,
285 const int n_alpha,
286 const int *alpha,
287 const double *x_k,
288 double *x_alpha) {
289 return domainPoints<3>(N, n_x, n_alpha, alpha, x_k, x_alpha);
290}
291
292template <int D, bool GRAD_BASE>
293MoFEMErrorCode
294BernsteinBezier::baseFunctions(const int N, const int gdim, const int n_alpha,
295 const int *alpha, const double *lambda,
296 const double *grad_lambda, double *base,
297 double *grad_base) {
298 MoFEMFunctionBeginHot;
299
300 const int *const alpha0 = alpha;
301 const double fN = boost::math::factorial<double>(N);
302 constexpr int MAX_ALPHA = 12;
303 int max_alpha = *std::max_element(alpha, &alpha[(D + 1) * n_alpha]);
304 if(max_alpha > MAX_ALPHA)
305 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
306 "Is assumed maximal order not to be bigger than %d", MAX_ALPHA);
307 std::array<double, (D + 1) * (MAX_ALPHA + 1)> pow_alpha;
308 std::array<double, (MAX_ALPHA + 1)> factorial_alpha;
309 std::array<double, D + 1> terms, diff_terms;
310 const double *const grad_lambda0 = grad_lambda;
311
312 factorial_alpha[0] = 1;
313 if (max_alpha >= 1)
314 factorial_alpha[1] = 1;
315 if (max_alpha >= 2)
316 for (int a = 2; a <= max_alpha; ++a)
317 factorial_alpha[a] = factorial_alpha[a - 1] * a;
318
319 for (int g = 0; g != gdim; ++g) {
320
321 for (int n = 0; n != D + 1; ++n, ++lambda) {
322 const size_t shift = (MAX_ALPHA + 1) * n;
323 double *pow_alpha_ptr = &pow_alpha[shift];
324 *pow_alpha_ptr = 1;
325
326 if (max_alpha >= 1) {
327 ++pow_alpha_ptr;
328 *pow_alpha_ptr = *lambda;
329 }
330
331 if (max_alpha >= 2)
332 for (int a = 2; a <= max_alpha; ++a) {
333 const double p = (*pow_alpha_ptr) * (*lambda);
334 ++pow_alpha_ptr;
335 *pow_alpha_ptr = p;
336 }
337 }
338
339 for (int n0 = 0; n0 != n_alpha; ++n0) {
340
341 grad_lambda = grad_lambda0;
342
343 double f = factorial_alpha[(*alpha)];
344 double *terms_ptr = terms.data();
345 double *diff_terms_ptr = diff_terms.data();
346 *terms_ptr = pow_alpha[(*alpha)];
347 if (GRAD_BASE) {
348 if (*alpha > 0)
349 *diff_terms_ptr = (*alpha) * pow_alpha[(*alpha) - 1];
350 else
351 *diff_terms_ptr = 0;
352 }
353 *base = *terms_ptr;
354 ++alpha;
355 ++terms_ptr;
356 ++diff_terms_ptr;
357
358 for (int n1 = 1; n1 < D + 1;
359 ++n1, ++alpha, ++terms_ptr, ++diff_terms_ptr) {
360 f *= factorial_alpha[(*alpha)];
361 const size_t shift = (MAX_ALPHA + 1) * n1;
362 *terms_ptr = pow_alpha[shift + (*alpha)];
363 if (GRAD_BASE) {
364 if (*alpha > 0)
365 *diff_terms_ptr = (*alpha) * pow_alpha[shift + (*alpha) - 1];
366 else
367 *diff_terms_ptr = 0;
368 }
369 *base *= *terms_ptr;
370 }
371
372 const double b = fN / f;
373 *base *= b;
374 ++base;
375
376 if (GRAD_BASE) {
377 double *terms_ptr = terms.data();
378 double *diff_terms_ptr = diff_terms.data();
379 double z = *diff_terms_ptr;
380 ++terms_ptr;
381 ++diff_terms_ptr;
382 for (int n2 = 1; n2 != D + 1; ++n2, ++terms_ptr)
383 z *= *terms_ptr;
384
385 for (int d = 0; d != D; ++d, ++grad_lambda)
386 grad_base[d] = z * (*grad_lambda);
387
388 for (int n1 = 1; n1 < D + 1; ++n1) {
389 z = *diff_terms_ptr;
390
391 int n2 = 0;
392 for (terms_ptr = terms.data(); n2 != n1; ++n2, ++terms_ptr)
393 z *= *terms_ptr;
394
395 ++n2;
396 ++terms_ptr;
397
398 for (; n2 < D + 1; ++n2, ++terms_ptr)
399 z *= *terms_ptr;
400
401 for (int d = 0; d != D; ++d, ++grad_lambda)
402 grad_base[d] += z * (*grad_lambda);
403
404 ++diff_terms_ptr;
405 }
406 for (int d = 0; d != D; ++d)
407 grad_base[d] *= b;
408 grad_base += D;
409 }
410
411 }
412
413 alpha = alpha0;
414 }
415
416 MoFEMFunctionReturnHot(0);
417}
418
419MoFEMErrorCode BernsteinBezier::baseFunctionsEdge(
420 const int N, const int gdim, const int n_alpha, const int *alpha,
421 const double *lambda, const double *grad_lambda, double *base,
422 double *grad_base) {
423 if (grad_base)
424 return baseFunctions<1, true>(N, gdim, n_alpha, alpha, lambda, grad_lambda,
425 base, grad_base);
426 else
427 return baseFunctions<1, false>(N, gdim, n_alpha, alpha, lambda, grad_lambda,
428 base, grad_base);
429}
430
431MoFEMErrorCode BernsteinBezier::baseFunctionsTri(
432 const int N, const int gdim, const int n_alpha, const int *alpha,
433 const double *lambda, const double *grad_lambda, double *base,
434 double *grad_base) {
435 if (grad_base)
436 return baseFunctions<2, true>(N, gdim, n_alpha, alpha, lambda, grad_lambda,
437 base, grad_base);
438 else
439 return baseFunctions<2, false>(N, gdim, n_alpha, alpha, lambda, grad_lambda,
440 base, grad_base);
441}
442
443MoFEMErrorCode BernsteinBezier::baseFunctionsTet(
444 const int N, const int gdim, const int n_alpha, const int *alpha,
445 const double *lambda, const double *grad_lambda, double *base,
446 double *grad_base) {
447 if (grad_base)
448 return baseFunctions<3, true>(N, gdim, n_alpha, alpha, lambda, grad_lambda,
449 base, grad_base);
450 else
451 return baseFunctions<3, false>(N, gdim, n_alpha, alpha, lambda, grad_lambda,
452 base, grad_base);
453}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
lambda
constexpr double lambda
surface tension
Definition: free_surface.cpp:165
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
D
double D
Definition: initial_diffusion.cpp:44
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
g
constexpr double g
Definition: shallow_wave.cpp:63
N
const int N
Definition: speed_test.cpp:3
Generated by Doxygen 1.9.5 and hosted at University of Glasgow