v0.14.0
Loading...
Searching...
No Matches
base_functions.c
Go to the documentation of this file.
1/** \file base_functions.c
2
3*/
4
5#include <cblas.h>
6#include <petscsys.h>
7#include <phg-quadrule/quad.h>
8
9#include <definitions.h>
10
11#include <base_functions.h>
12
13static PetscErrorCode ierr;
14
15PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L,
16 double *diffL, const int dim) {
17 MoFEMFunctionBeginHot;
18#ifndef NDEBUG
19 if (dim < 1)
20 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
21 if (dim > 3)
22 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
23 if (p < 0)
24 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
25#endif // NDEBUG
26
27 L[0] = 1;
28 if (diffL != NULL) {
29 for (int d = 0; d != dim; ++d) {
30 diffL[d * (p + 1) + 0] = 0;
31 }
32 }
33 if (p == 0)
34 MoFEMFunctionReturnHot(0);
35
36 L[1] = s;
37 if (diffL != NULL) {
38#ifndef NDEBUG
39 if (diff_s == NULL) {
40 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
41 }
42#endif // NDEBUG
43 for (int d = 0; d != dim; ++d) {
44 diffL[d * (p + 1) + 1] = diff_s[d];
45 }
46 }
47 if (p == 1)
48 MoFEMFunctionReturnHot(0);
49
50 int l = 1;
51 for (; l < p; l++) {
52 double A = ((2 * (double)l + 1) / ((double)l + 1));
53 double B = ((double)l / ((double)l + 1));
54 L[l + 1] = A * s * L[l] - B * L[l - 1];
55 if (diffL != NULL) {
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];
60 }
61 }
62 }
63
64 MoFEMFunctionReturnHot(0);
65}
66
67PetscErrorCode Jacobi_polynomials(int p, double alpha, double x, double t,
68 double *diff_x, double *diff_t, double *L,
69 double *diffL, const int dim) {
70 MoFEMFunctionBeginHot;
71#ifndef NDEBUG
72 if (dim < 1)
73 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
74 if (dim > 3)
75 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
76 if (p < 0)
77 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
78
79 if (diffL != NULL) {
80 if (diff_x == NULL) {
81 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
82 }
83 }
84
85#endif // NDEBUG
86 L[0] = 1;
87 if (diffL != NULL) {
88 diffL[0 * (p + 1) + 0] = 0;
89 if (dim >= 2) {
90 diffL[1 * (p + 1) + 0] = 0;
91 if (dim == 3) {
92 diffL[2 * (p + 1) + 0] = 0;
93 }
94 }
95 }
96 if (p == 0)
97 MoFEMFunctionReturnHot(0);
98 L[1] = 2 * x - t + alpha * x;
99 if (diffL != NULL) {
100 int d = 0;
101 for (; d < dim; ++d) {
102 double d_t = (diff_t) ? diff_t[d] : 0;
103 diffL[d * (p + 1) + 1] = (2 + alpha) * diff_x[d] - d_t;
104 }
105 }
106 if (p == 1)
107 MoFEMFunctionReturnHot(0);
108 int l = 1;
109 for (; l < p; l++) {
110 int lp1 = l + 1;
111 double a = 2 * lp1 * (lp1 + alpha) * (2 * lp1 + alpha - 2);
112 double b = 2 * lp1 + alpha - 1;
113 double c = (2 * lp1 + alpha) * (2 * lp1 + alpha - 2);
114 double d = 2 * (lp1 + alpha - 1) * (lp1 - 1) * (2 * lp1 + alpha);
115 double A = b * (c * (2 * x - t) + alpha * alpha * t) / a;
116 double B = d * t * t / a;
117 L[lp1] = A * L[l] - B * L[l - 1];
118 if (diffL != NULL) {
119 int z = 0;
120 for (; z < dim; ++z) {
121 double d_t = (diff_t) ? diff_t[z] : 0;
122 double diffA =
123 b * (c * (2 * diff_x[z] - d_t) + alpha * alpha * d_t) / a;
124 double diffB = d * 2 * t * d_t / a;
125 diffL[z * (p + 1) + lp1] = A * diffL[z * (p + 1) + l] -
126 B * diffL[z * (p + 1) + l - 1] +
127 diffA * L[l] - diffB * L[l - 1];
128 }
129 }
130 }
131 MoFEMFunctionReturnHot(0);
132}
133
134PetscErrorCode IntegratedJacobi_polynomials(int p, double alpha, double x,
135 double t, double *diff_x,
136 double *diff_t, double *L,
137 double *diffL, const int dim) {
138 MoFEMFunctionBeginHot;
139#ifndef NDEBUG
140 if (dim < 1)
141 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
142 if (dim > 3)
143 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
144 if (p < 1)
145 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 1");
146#endif // NDEBUG
147 L[0] = x;
148 if (diffL != NULL) {
149 int d = 0;
150 for (; d != dim; ++d) {
151 diffL[d * p + 0] = diff_x[d];
152 }
153 }
154 if (p == 0)
155 MoFEMFunctionReturnHot(0);
156 double jacobi[(p + 1)];
157 double diff_jacobi[(p + 1) * dim];
158 ierr = Jacobi_polynomials(p, alpha, x, t, diff_x, diff_t, jacobi, diff_jacobi,
159 dim);
160 CHKERRQ(ierr);
161 int l = 1;
162 for (; l < p; l++) {
163 int i = l + 1;
164 const double a = (i + alpha) / ((2 * i + alpha - 1) * (2 * i + alpha));
165 const double b = alpha / ((2 * i + alpha - 2) * (2 * i + alpha));
166 const double c = (i - 1) / ((2 * i + alpha - 2) * (2 * i + alpha - 1));
167 L[l] = a * jacobi[i] + b * t * jacobi[i - 1] - c * t * t * jacobi[i - 2];
168 if (diffL != NULL) {
169 int d = 0;
170 for (; d != dim; ++d) {
171 diffL[d * p + l] = a * diff_jacobi[d * (p + 1) + i] +
172 b * (t * diff_jacobi[d * (p + 1) + i - 1] +
173 diff_t[d] * jacobi[i - 1]) -
174 c * (t * t * diff_jacobi[d * (p + 1) + i - 2] +
175 2 * t * diff_t[d] * jacobi[i - 2]);
176 }
177 }
178 }
179 MoFEMFunctionReturnHot(0);
180}
181
182PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L,
183 double *diffL, const int dim) {
184
185 MoFEMFunctionBeginHot;
186#ifndef NDEBUG
187 if (dim < 1)
188 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
189 if (dim > 3)
190 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
191 if (p < 2)
192 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 2");
193#endif // NDEBUG
194 double l[p + 1];
195 ierr = Legendre_polynomials(p, s, NULL, l, NULL, 1);
196 CHKERRQ(ierr);
197
198 L[0] = 1;
199 if (diffL != NULL) {
200 for (int d = 0; d != dim; ++d) {
201 diffL[d * (p + 1) + 0] = 0;
202 }
203 }
204 L[1] = s;
205 if (diffL != NULL) {
206#ifndef NDEBUG
207 if (diff_s == NULL) {
208 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
209 }
210#endif // NDEBUG
211 for (int d = 0; d != dim; ++d) {
212 diffL[d * (p + 1) + 1] = diff_s[d];
213 }
214 }
215
216 // Integrated Legendre
217 for (int k = 2; k <= p; k++) {
218 const double factor = 2 * (2 * k - 1);
219 L[k] = 1.0 / factor * (l[k] - l[k - 2]);
220 }
221
222 if (diffL != NULL) {
223 for (int k = 2; k <= p; k++) {
224 double a = l[k - 1] / 2.;
225 for (int d = 0; d != dim; ++d) {
226 diffL[d * (p + 1) + k] = a * diff_s[d];
227 }
228 }
229 }
230 MoFEMFunctionReturnHot(0);
231}
232
233static double f_phi0(double x) { return LOBATTO_PHI0(x); }
234static double f_phi1(double x) { return LOBATTO_PHI1(x); }
235static double f_phi2(double x) { return LOBATTO_PHI2(x); }
236static double f_phi3(double x) { return LOBATTO_PHI3(x); }
237static double f_phi4(double x) { return LOBATTO_PHI4(x); }
238static double f_phi5(double x) { return LOBATTO_PHI5(x); }
239static double f_phi6(double x) { return LOBATTO_PHI6(x); }
240static double f_phi7(double x) { return LOBATTO_PHI7(x); }
241static double f_phi8(double x) { return LOBATTO_PHI8(x); }
242static double f_phi9(double x) { return LOBATTO_PHI9(x); }
243
244static double (*f_phi[])(double x) = {f_phi0, f_phi1, f_phi2, f_phi3, f_phi4,
245 f_phi5, f_phi6, f_phi7, f_phi8, f_phi9};
246
247static double f_phi0x(double x) { return LOBATTO_PHI0X(x); }
248static double f_phi1x(double x) { return LOBATTO_PHI1X(x); }
249static double f_phi2x(double x) { return LOBATTO_PHI2X(x); }
250static double f_phi3x(double x) { return LOBATTO_PHI3X(x); }
251static double f_phi4x(double x) { return LOBATTO_PHI4X(x); }
252static double f_phi5x(double x) { return LOBATTO_PHI5X(x); }
253static double f_phi6x(double x) { return LOBATTO_PHI6X(x); }
254static double f_phi7x(double x) { return LOBATTO_PHI7X(x); }
255static double f_phi8x(double x) { return LOBATTO_PHI8X(x); }
256static double f_phi9x(double x) { return LOBATTO_PHI9X(x); }
257
261
262// /// Legendre polynomials
263// #define Legendre0(x) (1.0)
264// #define Legendre1(x) (x)
265// #define Legendre2(x) (1.0 / 2.0 * (3 * (x) * (x) - 1))
266// #define Legendre3(x) (1.0 / 2.0 * (5 * (x) * (x) - 3) * (x))
267// #define Legendre4(x) (1.0 / 8.0 * ((35 * (x) * (x) - 30) * (x) * (x) + 3))
268// #define Legendre5(x) (1.0 / 8.0 * ((63 * (x) * (x) - 70) * (x) * (x) + 15) *
269// (x)) #define Legendre6(x) (1.0 / 16.0 * (((231 * (x) * (x) - 315) * (x) * (x)
270// + 105) * (x) * (x) - 5)) #define Legendre7(x) (1.0 / 16.0 * (((429 * (x) *
271// (x) - 693) * (x) * (x) + 315) * (x) * (x) - 35) * (x)) #define Legendre8(x)
272// (1.0 / 128.0 * ((((6435 * (x) * (x) - 12012) * (x) * (x) + 6930) * (x) * (x)
273// - 1260) * (x) * (x) + 35)) #define Legendre9(x) (1.0 / 128.0 * ((((12155 *
274// (x) * (x) - 25740) * (x) * (x) + 18018) * (x) * (x) - 4620) * (x) * (x) +
275// 315) * (x)) #define Legendre10(x) (1.0 / 256.0 * (((((46189 * (x) * (x) -
276// 109395) * (x) * (x) + 90090) * (x) * (x) - 30030) * (x) * (x) + 3465) * (x) *
277// (x) - 63))
278//
279// /// derivatives of Legendre polynomials
280// #define Legendre0x(x) (0.0)
281// #define Legendre1x(x) (1.0)
282// #define Legendre2x(x) (3.0 * (x))
283// #define Legendre3x(x) (15.0 / 2.0 * (x) * (x) - 3.0 / 2.0)
284// #define Legendre4x(x) (5.0 / 2.0 * (x) * (7.0 * (x) * (x) - 3.0))
285// #define Legendre5x(x) ((315.0 / 8.0 * (x) * (x) - 105.0 / 4.0) * (x) * (x)
286// + 15.0 / 8.0) #define Legendre6x(x) (21.0 / 8.0 * (x) * ((33.0 * (x) * (x)
287// - 30.0) * (x) * (x) + 5.0)) #define Legendre7x(x) (((3003.0 / 16.0 * (x) *
288// (x) - 3465.0 / 16.0) * (x) * (x) + 945.0 / 16.0) * (x) * (x) - 35.0 / 16.0)
289// #define Legendre8x(x) (9.0 / 16.0 * (x) * (((715.0 * (x) * (x) - 1001.0) *
290// (x) * (x) + 385.0) * (x) * (x) - 35.0)) #define Legendre9x(x) ((((109395.0 /
291// 128.0 * (x) * (x) - 45045.0 / 32.0) * (x) * (x) + 45045.0 / 64.0) * (x) * (x)
292// - 3465.0 / 32.0) * (x) * (x) + 315.0 / 128.0) #define Legendre10x(x) (2.0 /
293// 256.0 * (x) * ((((230945.0 * (x) * (x) - 437580.0) * (x) * (x) + 270270.0) *
294// (x) * (x) - 60060.0) * (x) * (x) + 3465.0))
295//
296// /// first two Lobatto shape functions
297// #define l0(x) ((1.0 - (x)) * 0.5)
298// #define l1(x) ((1.0 + (x)) * 0.5)
299//
300// #define l0l1(x) ((1.0 - (x)*(x)) * 0.25)
301//
302// /// other Lobatto shape functions
303// #define l2(x) (phi0(x) * l0l1(x))
304// #define l3(x) (phi1(x) * l0l1(x))
305// #define l4(x) (phi2(x) * l0l1(x))
306// #define l5(x) (phi3(x) * l0l1(x))
307// #define l6(x) (phi4(x) * l0l1(x))
308// #define l7(x) (phi5(x) * l0l1(x))
309// #define l8(x) (phi6(x) * l0l1(x))
310// #define l9(x) (phi7(x) * l0l1(x))
311// #define l10(x) (phi8(x) * l0l1(x))
312// #define l11(x) (phi9(x) * l0l1(x))
313//
314// /// derivatives of Lobatto functions
315// #define dl0(x) (-0.5)
316// #define dl1(x) (0.5)
317// #define dl2(x) (sqrt(3.0/2.0) * Legendre1(x))
318// #define dl3(x) (sqrt(5.0/2.0) * Legendre2(x))
319// #define dl4(x) (sqrt(7.0/2.0) * Legendre3(x))
320// #define dl5(x) (sqrt(9.0/2.0) * Legendre4(x))
321// #define dl6(x) (sqrt(11.0/2.0) * Legendre5(x))
322// #define dl7(x) (sqrt(13.0/2.0) * Legendre6(x))
323// #define dl8(x) (sqrt(15.0/2.0) * Legendre7(x))
324// #define dl9(x) (sqrt(17.0/2.0) * Legendre8(x))
325// #define dl10(x) (sqrt(19.0/2.0) * Legendre9(x))
326// #define dl11(x) (sqrt(21.0/2.0) * Legendre10(x))
327
328PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s,
329 double *L, double *diffL,
330 const int dim) {
331 MoFEMFunctionBeginHot;
332#ifndef NDEBUG
333 if (dim < 1)
334 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
335 if (dim > 3)
336 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
337 if (p < 0)
338 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
339 if (p > 9)
340 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
341 "Polynomial beyond order 9 is not implemented");
342 if (diffL != NULL) {
343 if (diff_s == NULL) {
344 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
345 }
346 }
347
348#endif // NDEBUG
349 if (L) {
350 int l = 0;
351 for (; l != p + 1; l++) {
352 L[l] = f_phi[l](s);
353 }
354 }
355 if (diffL != NULL) {
356 int l = 0;
357 for (; l != p + 1; l++) {
358 double a = f_phix[l](s);
359 int d = 0;
360 for (; d < dim; ++d) {
361 diffL[d * (p + 1) + l] = diff_s[d] * a;
362 }
363 }
364 }
365 MoFEMFunctionReturnHot(0);
366}
a
constexpr double a
Definition approx_sphere.cpp:30
f_phi4
static double f_phi4(double x)
Definition base_functions.c:237
f_phi1
static double f_phi1(double x)
Definition base_functions.c:234
f_phi0
static double f_phi0(double x)
Definition base_functions.c:233
f_phi8x
static double f_phi8x(double x)
Definition base_functions.c:255
IntegratedJacobi_polynomials
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.
Definition base_functions.c:134
f_phi
static double(* f_phi[])(double x)
Definition base_functions.c:244
f_phi7
static double f_phi7(double x)
Definition base_functions.c:240
f_phi0x
static double f_phi0x(double x)
Definition base_functions.c:247
f_phi4x
static double f_phi4x(double x)
Definition base_functions.c:251
ierr
static PetscErrorCode ierr
Definition base_functions.c:13
f_phi8
static double f_phi8(double x)
Definition base_functions.c:241
f_phix
static double(* f_phix[])(double x)
Definition base_functions.c:258
Jacobi_polynomials
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.
Definition base_functions.c:67
f_phi2x
static double f_phi2x(double x)
Definition base_functions.c:249
f_phi1x
static double f_phi1x(double x)
Definition base_functions.c:248
f_phi9x
static double f_phi9x(double x)
Definition base_functions.c:256
f_phi5
static double f_phi5(double x)
Definition base_functions.c:238
f_phi7x
static double f_phi7x(double x)
Definition base_functions.c:254
f_phi3x
static double f_phi3x(double x)
Definition base_functions.c:250
f_phi3
static double f_phi3(double x)
Definition base_functions.c:236
f_phi5x
static double f_phi5x(double x)
Definition base_functions.c:252
f_phi6x
static double f_phi6x(double x)
Definition base_functions.c:253
f_phi9
static double f_phi9(double x)
Definition base_functions.c:242
f_phi2
static double f_phi2(double x)
Definition base_functions.c:235
f_phi6
static double f_phi6(double x)
Definition base_functions.c:239
base_functions.h
LOBATTO_PHI2
#define LOBATTO_PHI2(x)
Definition base_functions.h:128
LOBATTO_PHI5
#define LOBATTO_PHI5(x)
Definition base_functions.h:135
LOBATTO_PHI8
#define LOBATTO_PHI8(x)
Definition base_functions.h:144
LOBATTO_PHI3
#define LOBATTO_PHI3(x)
Definition base_functions.h:130
LOBATTO_PHI7X
#define LOBATTO_PHI7X(x)
Definition base_functions.h:172
LOBATTO_PHI0X
#define LOBATTO_PHI0X(x)
Derivatives of kernel functions for Lobbatto base.
Definition base_functions.h:157
LOBATTO_PHI2X
#define LOBATTO_PHI2X(x)
Definition base_functions.h:159
LOBATTO_PHI5X
#define LOBATTO_PHI5X(x)
Definition base_functions.h:166
LOBATTO_PHI0
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
Definition base_functions.h:126
LOBATTO_PHI1
#define LOBATTO_PHI1(x)
Definition base_functions.h:127
LOBATTO_PHI9X
#define LOBATTO_PHI9X(x)
Definition base_functions.h:179
LOBATTO_PHI4
#define LOBATTO_PHI4(x)
Definition base_functions.h:132
LOBATTO_PHI8X
#define LOBATTO_PHI8X(x)
Definition base_functions.h:175
LOBATTO_PHI7
#define LOBATTO_PHI7(x)
Definition base_functions.h:141
LOBATTO_PHI6
#define LOBATTO_PHI6(x)
Definition base_functions.h:138
LOBATTO_PHI1X
#define LOBATTO_PHI1X(x)
Definition base_functions.h:158
LOBATTO_PHI9
#define LOBATTO_PHI9(x)
Definition base_functions.h:149
LOBATTO_PHI4X
#define LOBATTO_PHI4X(x)
Definition base_functions.h:163
LOBATTO_PHI6X
#define LOBATTO_PHI6X(x)
Definition base_functions.h:169
LOBATTO_PHI3X
#define LOBATTO_PHI3X(x)
Definition base_functions.h:161
definitions.h
useful compiler derivatives and definitions
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:460
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:453
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition base_functions.c:15
Lobatto_polynomials
PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Lobatto base functions fuentes2015353.
Definition base_functions.c:182
LobattoKernel_polynomials
PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Kernel Lobatto base functions.
Definition base_functions.c:328
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
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
t
constexpr double t
plate stiffness
Definition plate.cpp:58
Generated by Doxygen 1.12.0 and hosted at University of Glasgow