DUBeat  1.0.1
High-order discontinuous Galerkin methods and applications to cardiac electrophysiology
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 
48 template <unsigned int dim>
49 class DUBValues
50 {
51 protected:
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 
74 public:
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 
119 template <>
120 std::array<unsigned int, 2>
121 DUBValues<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 
141 template <>
142 std::array<unsigned int, 3>
143 DUBValues<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 
174 template <unsigned int dim>
175 double
176 DUBValues<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 
221 template <unsigned int dim>
222 unsigned int
223 DUBValues<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 
241 template <unsigned int dim>
242 unsigned int
243 DUBValues<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 
262 template <>
263 double
264 DUBValues<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 
297 template <>
298 double
299 DUBValues<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 
339 template <>
340 dealii::Tensor<1, 2>
341 DUBValues<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 
409 template <>
410 dealii::Tensor<1, 3>
411 DUBValues<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::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.
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::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::~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
generate_convergence_plots.k
k
Definition: generate_convergence_plots.py:49
generate_convergence_plots.n
n
Definition: generate_convergence_plots.py:41