16 double *diffL,
const int dim) {
29 for (
int d = 0; d !=
dim; ++d) {
30 diffL[d * (
p + 1) + 0] = 0;
43 for (
int d = 0; d !=
dim; ++d) {
44 diffL[d * (
p + 1) + 1] = diff_s[d];
54 L[
l + 1] =
A * s *
L[
l] -
B *
L[
l - 1];
56 for (
int d = 0; d !=
dim; ++d) {
57 diffL[d * (
p + 1) +
l + 1] =
58 A * (diff_s[d] *
L[
l] + s * diffL[d * (
p + 1) +
l]) -
59 B * diffL[d * (
p + 1) +
l - 1];
68 double *diff_x,
double *diff_t,
double *
L,
69 double *diffL,
const int dim) {
81 diffL[0 * (
p + 1) + 0] = 0;
83 diffL[1 * (
p + 1) + 0] = 0;
85 diffL[2 * (
p + 1) + 0] = 0;
91 L[1] = 2 * x -
t + alpha * x;
98 double d_t = (diff_t) ? diff_t[0] : 0;
99 diffL[0 * (
p + 1) + 1] = (2 + alpha) * diff_x[0] - d_t;
101 double d_t = (diff_t) ? diff_t[1] : 0;
102 diffL[1 * (
p + 1) + 1] = (2 + alpha) * diff_x[1] - d_t;
104 double d_t = (diff_t) ? diff_t[2] : 0;
105 diffL[2 * (
p + 1) + 1] = (2 + alpha) * diff_x[2] - d_t;
114 double a = 2 * lp1 * (lp1 + alpha) * (2 * lp1 + alpha - 2);
115 double b = 2 * lp1 + alpha - 1;
116 double c = (2 * lp1 + alpha) * (2 * lp1 + alpha - 2);
117 double d = 2 * (lp1 + alpha - 1) * (lp1 - 1) * (2 * lp1 + alpha);
118 double A = b * (
c * (2 * x -
t) + alpha * alpha *
t) /
a;
119 double B = d *
t *
t /
a;
120 L[lp1] =
A *
L[
l] -
B *
L[
l - 1];
122 double d_t = (diff_t) ? diff_t[0] : 0;
123 double diffA = b * (
c * (2 * diff_x[0] - d_t) + alpha * alpha * d_t) /
a;
124 double diffB = d * 2 *
t * d_t /
a;
125 diffL[0 * (
p + 1) + lp1] =
A * diffL[0 * (
p + 1) +
l] -
126 B * diffL[0 * (
p + 1) +
l - 1] + diffA *
L[
l] -
129 double d_t = (diff_t) ? diff_t[1] : 0;
131 b * (
c * (2 * diff_x[1] - d_t) + alpha * alpha * d_t) /
a;
132 double diffB = d * 2 *
t * d_t /
a;
133 diffL[1 * (
p + 1) + lp1] =
A * diffL[1 * (
p + 1) +
l] -
134 B * diffL[1 * (
p + 1) +
l - 1] +
135 diffA *
L[
l] - diffB *
L[
l - 1];
137 double d_t = (diff_t) ? diff_t[2] : 0;
139 b * (
c * (2 * diff_x[2] - d_t) + alpha * alpha * d_t) /
a;
140 double diffB = d * 2 *
t * d_t /
a;
141 diffL[2 * (
p + 1) + lp1] =
A * diffL[2 * (
p + 1) +
l] -
142 B * diffL[2 * (
p + 1) +
l - 1] +
143 diffA *
L[
l] - diffB *
L[
l - 1];
152 double t,
double *diff_x,
153 double *diff_t,
double *
L,
154 double *diffL,
const int dim) {
167 for (; d !=
dim; ++d) {
168 diffL[d *
p + 0] = diff_x[d];
173 double jacobi[(
p + 1)];
174 double diff_jacobi[(
p + 1) *
dim];
181 const double a = (
i + alpha) / ((2 *
i + alpha - 1) * (2 *
i + alpha));
182 const double b = alpha / ((2 *
i + alpha - 2) * (2 *
i + alpha));
183 const double c = (
i - 1) / ((2 *
i + alpha - 2) * (2 *
i + alpha - 1));
184 L[
l] =
a * jacobi[
i] + b *
t * jacobi[
i - 1] -
c *
t *
t * jacobi[
i - 2];
187 for (; dd !=
dim; ++dd) {
188 diffL[dd *
p +
l] =
a * diff_jacobi[dd * (
p + 1) +
i] +
189 b * (
t * diff_jacobi[dd * (
p + 1) +
i - 1] +
190 diff_t[dd] * jacobi[
i - 1]) -
191 c * (
t *
t * diff_jacobi[dd * (
p + 1) +
i - 2] +
192 2 *
t * diff_t[dd] * jacobi[
i - 2]);
200 double *diffL,
const int dim) {
217 for (
int d = 0; d !=
dim; ++d) {
218 diffL[d * (
p + 1) + 0] = 0;
224 if (diff_s == NULL) {
228 for (
int d = 0; d !=
dim; ++d) {
229 diffL[d * (
p + 1) + 1] = diff_s[d];
234 for (
int k = 2;
k <=
p;
k++) {
235 const double factor = 2 * (2 *
k - 1);
236 L[
k] = 1.0 / factor * (
l[
k] -
l[
k - 2]);
240 for (
int k = 2;
k <=
p;
k++) {
241 double a =
l[
k - 1] / 2.;
242 for (
int d = 0; d !=
dim; ++d) {
243 diffL[d * (
p + 1) +
k] =
a * diff_s[d];
346 double *
L,
double *diffL,
358 "Polynomial beyond order 9 is not implemented");
362 for (;
l !=
p + 1;
l++) {
368 if (diff_s == NULL) {
373 for (;
l !=
p + 1;
l++) {
375 diffL[0 * (
p + 1) +
l] = diff_s[0] *
a;
377 diffL[1 * (
p + 1) +
l] = diff_s[1] *
a;
379 diffL[2 * (
p + 1) +
l] = diff_s[2] *
a;
static double f_phi4(double x)
static double f_phi1(double x)
static double f_phi0(double x)
static double f_phi8x(double x)
PetscErrorCode IntegratedJacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate integrated Jacobi approximation basis.
static double(* f_phi[])(double x)
static double f_phi7(double x)
static double f_phi0x(double x)
static double f_phi4x(double x)
static PetscErrorCode ierr
static double f_phi8(double x)
static double(* f_phix[])(double x)
PetscErrorCode Jacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate Jacobi approximation basis.
static double f_phi2x(double x)
static double f_phi1x(double x)
static double f_phi9x(double x)
static double f_phi5(double x)
static double f_phi7x(double x)
static double f_phi3x(double x)
static double f_phi3(double x)
static double f_phi5x(double x)
static double f_phi6x(double x)
static double f_phi9(double x)
static double f_phi2(double x)
static double f_phi6(double x)
#define LOBATTO_PHI0X(x)
Derivatives of kernel functions for Lobbatto base.
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
useful compiler directives and definitions
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Lobatto base functions .
PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Kernel Lobatto base functions.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'k', 3 > k
constexpr double t
plate stiffness