16 double *diffL,
const int dim) {
28 diffL[0 * (
p + 1) + 0] = 0;
30 diffL[1 * (
p + 1) + 0] = 0;
32 diffL[2 * (
p + 1) + 0] = 0;
45 diffL[0 * (
p + 1) + 1] = diff_s[0];
47 diffL[1 * (
p + 1) + 1] = diff_s[1];
49 diffL[2 * (
p + 1) + 1] = diff_s[2];
59 L[
l + 1] =
A * s *
L[
l] -
B *
L[
l - 1];
61 diffL[0 * (
p + 1) +
l + 1] =
62 A * (s * diffL[0 * (
p + 1) +
l] + diff_s[0] *
L[
l]) -
63 B * diffL[0 * (
p + 1) +
l - 1];
65 diffL[1 * (
p + 1) +
l + 1] =
66 A * (s * diffL[1 * (
p + 1) +
l] + diff_s[1] *
L[
l]) -
67 B * diffL[1 * (
p + 1) +
l - 1];
69 diffL[2 * (
p + 1) +
l + 1] =
70 A * (s * diffL[2 * (
p + 1) +
l] + diff_s[2] *
L[
l]) -
71 B * diffL[2 * (
p + 1) +
l - 1];
80 double *diff_x,
double *diff_t,
double *
L,
81 double *diffL,
const int dim) {
93 diffL[0 * (
p + 1) + 0] = 0;
95 diffL[1 * (
p + 1) + 0] = 0;
97 diffL[2 * (
p + 1) + 0] = 0;
103 L[1] = 2 * x -
t + alpha * x;
106 if (diff_x == NULL) {
110 double d_t = (diff_t) ? diff_t[0] : 0;
111 diffL[0 * (
p + 1) + 1] = (2 + alpha) * diff_x[0] - d_t;
113 double d_t = (diff_t) ? diff_t[1] : 0;
114 diffL[1 * (
p + 1) + 1] = (2 + alpha) * diff_x[1] - d_t;
116 double d_t = (diff_t) ? diff_t[2] : 0;
117 diffL[2 * (
p + 1) + 1] = (2 + alpha) * diff_x[2] - d_t;
126 double a = 2 * lp1 * (lp1 + alpha) * (2 * lp1 + alpha - 2);
127 double b = 2 * lp1 + alpha - 1;
128 double c = (2 * lp1 + alpha) * (2 * lp1 + alpha - 2);
129 double d = 2 * (lp1 + alpha - 1) * (lp1 - 1) * (2 * lp1 + alpha);
130 double A = b * (
c * (2 * x -
t) + alpha * alpha *
t) /
a;
131 double B = d *
t *
t /
a;
132 L[lp1] =
A *
L[
l] -
B *
L[
l - 1];
134 double d_t = (diff_t) ? diff_t[0] : 0;
135 double diffA = b * (
c * (2 * diff_x[0] - d_t) + alpha * alpha * d_t) /
a;
136 double diffB = d * 2 *
t * d_t /
a;
137 diffL[0 * (
p + 1) + lp1] =
A * diffL[0 * (
p + 1) +
l] -
138 B * diffL[0 * (
p + 1) +
l - 1] + diffA *
L[
l] -
141 double d_t = (diff_t) ? diff_t[1] : 0;
143 b * (
c * (2 * diff_x[1] - d_t) + alpha * alpha * d_t) /
a;
144 double diffB = d * 2 *
t * d_t /
a;
145 diffL[1 * (
p + 1) + lp1] =
A * diffL[1 * (
p + 1) +
l] -
146 B * diffL[1 * (
p + 1) +
l - 1] +
147 diffA *
L[
l] - diffB *
L[
l - 1];
149 double d_t = (diff_t) ? diff_t[2] : 0;
151 b * (
c * (2 * diff_x[2] - d_t) + alpha * alpha * d_t) /
a;
152 double diffB = d * 2 *
t * d_t /
a;
153 diffL[2 * (
p + 1) + lp1] =
A * diffL[2 * (
p + 1) +
l] -
154 B * diffL[2 * (
p + 1) +
l - 1] +
155 diffA *
L[
l] - diffB *
L[
l - 1];
164 double t,
double *diff_x,
165 double *diff_t,
double *
L,
166 double *diffL,
const int dim) {
179 for (; d !=
dim; ++d) {
180 diffL[d *
p + 0] = diff_x[d];
185 double jacobi[(
p + 1)];
186 double diff_jacobi[(
p + 1) *
dim];
193 const double a = (
i + alpha) / ((2 *
i + alpha - 1) * (2 *
i + alpha));
194 const double b = alpha / ((2 *
i + alpha - 2) * (2 *
i + alpha));
195 const double c = (
i - 1) / ((2 *
i + alpha - 2) * (2 *
i + alpha - 1));
196 L[
l] =
a * jacobi[
i] + b *
t * jacobi[
i - 1] -
c *
t *
t * jacobi[
i - 2];
199 for (; dd !=
dim; ++dd) {
200 diffL[dd *
p +
l] =
a * diff_jacobi[dd * (
p + 1) +
i] +
201 b * (
t * diff_jacobi[dd * (
p + 1) +
i - 1] +
202 diff_t[dd] * jacobi[
i - 1]) -
203 c * (
t *
t * diff_jacobi[dd * (
p + 1) +
i - 2] +
204 2 *
t * diff_t[dd] * jacobi[
i - 2]);
212 double *diffL,
const int dim) {
229 diffL[0 * (
p + 1) + 0] = 0;
231 diffL[1 * (
p + 1) + 0] = 0;
233 diffL[2 * (
p + 1) + 0] = 0;
240 if (diff_s == NULL) {
244 diffL[0 * (
p + 1) + 1] = diff_s[0];
246 diffL[1 * (
p + 1) + 1] = diff_s[1];
248 diffL[2 * (
p + 1) + 1] = diff_s[2];
254 for (
int k = 2;
k <=
p;
k++) {
255 const double factor = 2 * (2 *
k - 1);
256 L[
k] = 1.0 / factor * (
l[
k] -
l[
k - 2]);
259 if (diff_s == NULL) {
262 double a =
l[
k - 1] / 2.;
263 diffL[0 * (
p + 1) +
k] =
a * diff_s[0];
265 diffL[1 * (
p + 1) +
k] =
a * diff_s[1];
267 diffL[2 * (
p + 1) +
k] =
a * diff_s[2];
371 double *
L,
double *diffL,
383 "Polynomial beyond order 9 is not implemented");
387 for (;
l !=
p + 1;
l++) {
393 if (diff_s == NULL) {
398 for (;
l !=
p + 1;
l++) {
400 diffL[0 * (
p + 1) +
l] = diff_s[0] *
a;
402 diffL[1 * (
p + 1) +
l] = diff_s[1] *
a;
404 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