DUBeat 1.0.1
High-order discontinuous Galerkin methods and applications to cardiac electrophysiology
Loading...
Searching...
No Matches
DUBValues.hpp
Go to the documentation of this file.
1/********************************************************************************
2 Copyright (C) 2024 by the DUBeat authors.
3
4 This file is part of DUBeat.
5
6 DUBeat is free software; you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 DUBeat is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public License
17 along with DUBeat. If not, see <http://www.gnu.org/licenses/>.
18********************************************************************************/
19
27#ifndef DUBValues_HPP_
28#define DUBValues_HPP_
29
30#include <deal.II/base/quadrature_lib.h>
31#include <deal.II/base/tensor.h>
32
33#include <deal.II/fe/fe.h>
34#include <deal.II/fe/fe_q.h>
35
36#include <math.h>
37
38#include <cmath>
39#include <utility>
40
48template <unsigned int dim>
49class DUBValues
50{
51protected:
53 const unsigned int poly_degree;
54
56 const double tol = 1e-10;
57
60 double
61 eval_jacobi_polynomial(const unsigned int n,
62 const unsigned int alpha,
63 const unsigned int beta,
64 const double eval_point) const;
65
71 std::array<unsigned int, dim>
72 fun_coeff_conversion(const unsigned int n) const;
73
74public:
76 DUBValues<dim>(const unsigned int degree)
77 : poly_degree(degree)
78 {
79 dofs_per_cell = get_dofs_per_cell();
80 }
81
83 unsigned int dofs_per_cell;
84
86 DUBValues<dim>(DUBValues<dim> &DUBValues) = default;
87
89 DUBValues<dim>(const DUBValues<dim> &DUBValues) = default;
90
92 DUBValues<dim>(DUBValues<dim> &&DUBValues) = default;
93
95 unsigned int
96 get_dofs_per_cell() const;
97
99 unsigned int
100 get_dofs_per_cell(const unsigned int degree) const;
101
103 double
104 shape_value(const unsigned int function_no,
105 const dealii::Point<dim> &quadrature_point) const;
106
108 dealii::Tensor<1, dim>
109 shape_grad(const unsigned int function_no,
110 const dealii::Point<dim> &quadrature_point) const;
111
113 virtual ~DUBValues() = default;
114};
115
119template <>
120std::array<unsigned int, 2>
121DUBValues<2>::fun_coeff_conversion(const unsigned int n) const
122{
123 unsigned int i, j;
124 unsigned int c = poly_degree + 1;
125 unsigned int c_seq = c;
126 while (n > c_seq - 1)
127 {
128 c -= 1;
129 c_seq += c;
130 }
131
132 i = n - (c_seq - c);
133 j = poly_degree + 1 - c;
134
135 return {{i, j}};
136}
137
141template <>
142std::array<unsigned int, 3>
143DUBValues<3>::fun_coeff_conversion(const unsigned int n) const
144{
145 unsigned int i, j, k;
146 unsigned int p = poly_degree;
147 unsigned int c1 = (p + 1) * (p + 2) / 2;
148 unsigned int c1_seq = c1;
149
150 while (n > c1_seq - 1)
151 {
152 p -= 1;
153 c1 = (p + 1) * (p + 2) / 2;
154 c1_seq += c1;
155 }
156
157 i = poly_degree - p;
158 unsigned int n2 = n - (c1_seq - c1) + 1;
159 unsigned int c2 = p + 1;
160 unsigned int c2_seq = c2;
161
162 while (n2 > c2_seq)
163 {
164 c2 -= 1;
165 c2_seq += c2;
166 }
167
168 j = poly_degree - i + 1 - c2;
169 k = n2 - (c2_seq - c2) - 1;
170
171 return {{i, j, k}};
172}
173
174template <unsigned int dim>
175double
176DUBValues<dim>::eval_jacobi_polynomial(const unsigned int n,
177 const unsigned int alpha,
178 const unsigned int beta,
179 const double eval_point) const
180{
181 if (n == 0)
182 return 1.0;
183 else if (n == 1)
184 return static_cast<double>(
185 (alpha - beta + (alpha + beta + 2.0) * eval_point) / 2.0);
186 else
187 {
188 const double apb = alpha + beta;
189 double poly = 0.0;
190 double polyn2 = 1.0;
191 double polyn1 = static_cast<double>(
192 (alpha - beta + (alpha + beta + 2.0) * eval_point) / 2.0);
193 double a1, a2, a3, a4;
194
195 for (unsigned int k = 2; k <= n; ++k)
196 {
197 a1 = 2.0 * static_cast<double>(k) * static_cast<double>(k + apb) *
198 static_cast<double>(2.0 * k + apb - 2.0);
199 a2 = static_cast<double>(2.0 * k + apb - 1.0) *
200 static_cast<double>(alpha * alpha - beta * beta);
201 a3 = static_cast<double>(2.0 * k + apb - 2.0) *
202 static_cast<double>(2.0 * k + apb - 1.0) *
203 static_cast<double>(2.0 * k + apb);
204 a4 = 2.0 * static_cast<double>(k + alpha - 1.0) *
205 static_cast<double>(k + beta - 1.0) *
206 static_cast<double>(2.0 * k + apb);
207
208 a2 = a2 / a1;
209 a3 = a3 / a1;
210 a4 = a4 / a1;
211
212 poly = (a2 + a3 * eval_point) * polyn1 - a4 * polyn2;
213 polyn2 = polyn1;
214 polyn1 = poly;
215 }
216
217 return poly;
218 }
219}
220
221template <unsigned int dim>
222unsigned int
223DUBValues<dim>::get_dofs_per_cell() const
224{
225 // The analytical formula is:
226 // n_dof_per_cell = (p+1)*(p+2)*...(p+d) / d!,
227 // where p is the space order and d the space dimension..
228
229 unsigned int denominator = 1;
230 unsigned int nominator = 1;
231
232 for (unsigned int i = 1; i <= dim; i++)
233 {
234 denominator *= i;
235 nominator *= poly_degree + i;
236 }
237
238 return (int)(nominator / denominator);
239}
240
241template <unsigned int dim>
242unsigned int
243DUBValues<dim>::get_dofs_per_cell(const unsigned int degree) const
244{
245 // The analytical formula is:
246 // n_dof_per_cell = (p+1)*(p+2)*...(p+d) / d!,
247 // where p is the space order and d the space dimension..
248
249 unsigned int denominator = 1;
250 unsigned int nominator = 1;
251
252 for (unsigned int i = 1; i <= dim; i++)
253 {
254 denominator *= i;
255 nominator *= degree + i;
256 }
257
258 return (int)(nominator / denominator);
259}
260
262template <>
263double
264DUBValues<2>::shape_value(const unsigned int function_no,
265 const dealii::Point<2> &quadrature_point) const
266{
267 AssertThrow(function_no < dofs_per_cell,
268 dealii::StandardExceptions::ExcMessage(
269 "function_no outside the limit."));
270
271 const auto fun_coeff = fun_coeff_conversion(function_no);
272
273 const int i = fun_coeff[0];
274 const int j = fun_coeff[1];
275 const double csi = quadrature_point[0];
276 const double eta = quadrature_point[1];
277
278 // If the point is outside the reference cell, the evaluation is zero.
279 if (csi + eta > 1.0 + tol || csi < -tol || eta < -tol)
280 return 0.0;
281
282 double a, b;
283 a = b = 0;
284
285 if (std::abs(1 - eta) > this->tol)
286 a = 2.0 * csi / (1 - eta) - 1.0;
287 b = 2.0 * eta - 1.0;
288
289 const double r = sqrt((2 * i + 1) * 2 * (i + j + 1) / pow(4, i));
290 const double pi = this->eval_jacobi_polynomial(i, 0, 0, a);
291 const double pj = this->eval_jacobi_polynomial(j, 2 * i + 1, 0, b);
292
293 return r * pow(2, i) * pow(1 - eta, i) * pi * pj;
294}
295
297template <>
298double
299DUBValues<3>::shape_value(const unsigned int function_no,
300 const dealii::Point<3> &quadrature_point) const
301{
302 AssertThrow(function_no < dofs_per_cell,
303 dealii::StandardExceptions::ExcMessage(
304 "function_no outside the limit."));
305
306 const auto fun_coeff = fun_coeff_conversion(function_no);
307
308 const int i = fun_coeff[0];
309 const int j = fun_coeff[1];
310 const int k = fun_coeff[2];
311 const double csi = quadrature_point[0];
312 const double eta = quadrature_point[1];
313 const double ni = quadrature_point[2];
314
315 // If the point is outside the reference cell, the evaluation is zero.
316 if (csi + eta + ni > 1.0 + tol || csi < -tol || eta < -tol || ni < -tol)
317 return 0.0;
318
319 double a, b, c;
320 a = b = c = 0;
321
322 if (std::abs(-eta - ni + 1.) > this->tol)
323 a = 2. * csi / (-eta - ni + 1.) - 1.;
324 if (std::abs(1 - ni) > this->tol)
325 b = 2. * eta / (1. - ni) - 1.;
326 c = 2. * ni - 1.;
327
328 const double r =
329 sqrt(8 * pow(2, 4 * i + 2 * j + 3) * (2 * i + 2 * j + 2 * k + 3) *
330 (2 * i + 2 * j + 2) * (2 * i + 1) / pow(2, 4 * i + 2 * j + 6));
331 const double pi = this->eval_jacobi_polynomial(i, 0, 0, a);
332 const double pj = this->eval_jacobi_polynomial(j, 2 * i + 1, 0, b);
333 const double pk = this->eval_jacobi_polynomial(k, 2 * i + 2 * j + 2, 0, c);
334
335 return r * pi * pow((1. - b) / 2., i) * pj * pow((1. - c) / 2., i + j) * pk;
336}
337
339template <>
340dealii::Tensor<1, 2>
341DUBValues<2>::shape_grad(const unsigned int function_no,
342 const dealii::Point<2> &quadrature_point) const
343{
344 AssertThrow(function_no < dofs_per_cell,
345 dealii::StandardExceptions::ExcMessage(
346 "function_no outside the limit."));
347
348 const auto fun_coeff = fun_coeff_conversion(function_no);
349
350 const int i = fun_coeff[0];
351 const int j = fun_coeff[1];
352 const double csi = quadrature_point[0];
353 const double eta = quadrature_point[1];
354
355 dealii::Tensor<1, 2> grad;
356
357 // If the point is outside the reference cell, the evaluation is zero.
358 if (csi + eta > 1.0 + tol || csi < -tol || eta < -tol)
359 return grad;
360
361 double a, b;
362 a = b = 0;
363
364 if (std::abs(1 - eta) > this->tol)
365 a = 2.0 * csi / (1 - eta) - 1.0;
366 b = 2.0 * eta - 1.0;
367
368 const double r = sqrt((2 * i + 1) * 2 * (i + j + 1) / pow(4, i));
369 const double pi = this->eval_jacobi_polynomial(i, 0, 0, a);
370 const double pj = this->eval_jacobi_polynomial(j, 2 * i + 1, 0, b);
371
372
373 if (i == 0 && j == 0)
374 {
375 grad[0] = 0.0;
376 grad[1] = 0.0;
377 }
378 else if (i == 0 && j != 0)
379 {
380 grad[0] = 0.0;
381 grad[1] = r * (j + 2) * eval_jacobi_polynomial(j - 1, 2, 1, b);
382 }
383 else if (i != 0 && j == 0)
384 {
385 grad[0] = r * pow(2, i) * pow(1 - eta, i - 1) * (i + 1) *
386 eval_jacobi_polynomial(i - 1, 1, 1, a);
387 grad[1] = r * pow(2, i) *
388 (-i * pow(1 - eta, i - 1) * pi +
389 csi * pow(1 - eta, i - 2) * (i + 1) *
390 eval_jacobi_polynomial(i - 1, 1, 1, a));
391 }
392 else
393 {
394 grad[0] = r * pow(2, i) * pow(1 - eta, i - 1) * (i + 1) *
395 eval_jacobi_polynomial(i - 1, 1, 1, a) *
396 eval_jacobi_polynomial(j, 2 * i + 1, 0, b);
397 grad[1] = r * pow(2, i) *
398 (-i * pow(1 - eta, i - 1) * pi * pj +
399 csi * pow(1 - eta, i - 2) * (i + 1) *
400 eval_jacobi_polynomial(i - 1, 1, 1, a) * pj +
401 pow(1 - eta, i) * (2 * i + j + 2) * pi *
402 eval_jacobi_polynomial(j - 1, 2 * i + 2, 1, b));
403 }
404
405 return grad;
406}
407
409template <>
410dealii::Tensor<1, 3>
411DUBValues<3>::shape_grad(const unsigned int function_no,
412 const dealii::Point<3> &quadrature_point) const
413{
414 AssertThrow(function_no < dofs_per_cell,
415 dealii::StandardExceptions::ExcMessage(
416 "function_no outside the limit."));
417
418 const auto fun_coeff = fun_coeff_conversion(function_no);
419
420 const int i = fun_coeff[0];
421 const int j = fun_coeff[1];
422 const int k = fun_coeff[2];
423 const double csi = quadrature_point[0];
424 const double eta = quadrature_point[1];
425 const double ni = quadrature_point[2];
426
427 dealii::Tensor<1, 3> grad;
428
429 // If the point is outside the reference cell, the evaluation is zero.
430 if (csi + eta + ni > 1.0 + tol || csi < -tol || eta < -tol || ni < -tol)
431 return grad;
432
433 double a, b, c;
434 a = b = c = 0.0;
435 if (std::abs(-eta - ni + 1.) > this->tol)
436 a = 2. * csi / (-eta - ni + 1.) - 1;
437 if (std::abs(1. - ni) > this->tol)
438 b = 2. * eta / (1. - ni) - 1.;
439 c = 2. * ni - 1.;
440
441 const double r =
442 sqrt(8 * pow(2, 4 * i + 2 * j + 3) * (2 * i + 2 * j + 2 * k + 3) *
443 (2 * i + 2 * j + 2) * (2 * i + 1) / pow(2, 4 * i + 2 * j + 6));
444 const double pi = this->eval_jacobi_polynomial(i, 0, 0, a);
445 const double pj = this->eval_jacobi_polynomial(j, 2 * i + 1, 0, b);
446 const double pk = this->eval_jacobi_polynomial(k, 2 * i + 2 * j + 2, 0, c);
447
448
449 double dPi_csi = 0.0, dPi_eta = 0.0, dPi_ni = 0.0;
450 double db_csi = 0.0, db_eta = 0.0, db_ni = 0.0;
451 double dPj_csi = 0.0, dPj_eta = 0.0, dPj_ni = 0.0;
452 double dc_csi = 0.0, dc_eta = 0.0, dc_ni = 0.0;
453 double dPk_csi = 0.0, dPk_eta = 0.0, dPk_ni = 0.0;
454
455 if (i != 0)
456 {
457 dPi_csi = static_cast<double>(i + 1.) * 1. / (-eta - ni + 1.) *
458 eval_jacobi_polynomial(i - 1, 1, 1, a);
459 dPi_eta = static_cast<double>(i + 1.) * csi / pow(-eta - ni + 1., 2.) *
460 eval_jacobi_polynomial(i - 1, 1, 1, a);
461 dPi_ni = static_cast<double>(i + 1.) * csi / pow(-eta - ni + 1., 2.) *
462 eval_jacobi_polynomial(i - 1, 1, 1, a);
463
464 db_csi = 0.0;
465 db_eta = -static_cast<double>(i) / (1. - ni) * pow((1. - b) / 2., i - 1.);
466 db_ni = -static_cast<double>(i) * eta / pow(1. - ni, 2) *
467 pow((1. - b) / 2, i - 1.);
468 }
469
470 if (j != 0)
471 {
472 dPj_csi = 0.0;
473 dPj_eta = static_cast<double>(2. * i + j + 2.) / (1. - ni) *
474 eval_jacobi_polynomial(j - 1, 2 * i + 2, 1, b);
475 dPj_ni = eta * static_cast<double>(2. * i + j + 2.) / pow(1. - ni, 2.) *
476 eval_jacobi_polynomial(j - 1, 2 * i + 2, 1, b);
477 }
478
479 if ((i + j) != 0)
480 {
481 dc_csi = 0.0;
482 dc_eta = 0.0;
483 dc_ni = -static_cast<double>(i + j) * pow((1. - c) / 2., i + j - 1.);
484 }
485
486 if (k != 0)
487 {
488 dPk_csi = 0.0;
489 dPk_eta = 0.0;
490 dPk_ni = static_cast<double>(2. * i + 2. * j + k + 3.) *
491 eval_jacobi_polynomial(k - 1, 2 * i + 2 * j + 3, 1, c);
492 }
493
494 grad[0] =
495 r * (dPi_csi * pow((1. - b) / 2., i) * pj * pow((1. - c) / 2., i + j) * pk +
496 pi * db_csi * pj * pow((1. - c) / 2., i + j) * pk +
497 pi * pow((1. - b) / 2., i) * dPj_csi * pow((1. - c) / 2., i + j) * pk +
498 pi * pow((1. - b) / 2., i) * pj * dc_csi * pk +
499 pi * pow((1. - b) / 2., i) * pj * pow((1. - c) / 2., i + j) * dPk_csi);
500
501 grad[1] =
502 r * (dPi_eta * pow((1. - b) / 2., i) * pj * pow((1. - c) / 2., i + j) * pk +
503 pi * db_eta * pj * pow((1. - c) / 2., i + j) * pk +
504 pi * pow((1. - b) / 2., i) * dPj_eta * pow((1. - c) / 2., i + j) * pk +
505 pi * pow((1. - b) / 2., i) * pj * dc_eta * pk +
506 pi * pow((1. - b) / 2., i) * pj * pow((1. - c) / 2., i + j) * dPk_eta);
507
508 grad[2] =
509 r * (dPi_ni * pow((1. - b) / 2., i) * pj * pow((1. - c) / 2., i + j) * pk +
510 pi * db_ni * pj * pow((1. - c) / 2., i + j) * pk +
511 pi * pow((1. - b) / 2., i) * dPj_ni * pow((1. - c) / 2., i + j) * pk +
512 pi * pow((1. - b) / 2., i) * pj * dc_ni * pk +
513 pi * pow((1. - b) / 2., i) * pj * pow((1. - c) / 2., i + j) * dPk_ni);
514
515 return grad;
516}
517
518#endif /* DUBValues_HPP_*/
DUBValues
Class representing the Dubiner basis functions definitions.
Definition DUBValues.hpp:50
DUBValues::fun_coeff_conversion
std::array< unsigned int, dim > fun_coeff_conversion(const unsigned int n) const
Conversion from the FEM basis taxonomy ( basis function, basis function ) to the Dubiner basis taxo...
DUBValues::eval_jacobi_polynomial
double eval_jacobi_polynomial(const unsigned int n, const unsigned int alpha, const unsigned int beta, const double eval_point) const
Evaluation of the -Jacobi polynomial (used to evaluate Dubiner basis).
Definition DUBValues.hpp:176
DUBValues::~DUBValues
virtual ~DUBValues()=default
Default destructor.
DUBValues::tol
const double tol
Default tolerance.
Definition DUBValues.hpp:56
DUBValues::poly_degree
const unsigned int poly_degree
Polynomial degree used.
Definition DUBValues.hpp:53
DUBValues::shape_value
double shape_value(const unsigned int function_no, const dealii::Point< dim > &quadrature_point) const
Evaluation of the Dubiner basis functions.
DUBValues::get_dofs_per_cell
unsigned int get_dofs_per_cell(const unsigned int degree) const
Same as get_dofs_per_cell() but with variable space degree.
Definition DUBValues.hpp:243
DUBValues::dofs_per_cell
unsigned int dofs_per_cell
Number of degrees of freedom.
Definition DUBValues.hpp:83
DUBValues::get_dofs_per_cell
unsigned int get_dofs_per_cell() const
Return the number of degrees of freedom per element.
Definition DUBValues.hpp:223
DUBValues::shape_grad
dealii::Tensor< 1, dim > shape_grad(const unsigned int function_no, const dealii::Point< dim > &quadrature_point) const
Evaluation of the Dubiner basis functions gradients.