v0.13.1
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) {
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 L[0] = 1;
27 if (diffL != NULL) {
28 diffL[0 * (p + 1) + 0] = 0;
29 if (dim >= 2) {
30 diffL[1 * (p + 1) + 0] = 0;
31 if (dim == 3) {
32 diffL[2 * (p + 1) + 0] = 0;
33 }
34 }
35 }
36 if (p == 0)
38 L[1] = s;
39 if (diffL != NULL) {
40#ifndef NDEBUG
41 if (diff_s == NULL) {
42 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
43 }
44#endif // NDEBUG
45 diffL[0 * (p + 1) + 1] = diff_s[0];
46 if (dim >= 2) {
47 diffL[1 * (p + 1) + 1] = diff_s[1];
48 if (dim == 3) {
49 diffL[2 * (p + 1) + 1] = diff_s[2];
50 }
51 }
52 }
53 if (p == 1)
55 int l = 1;
56 for (; l < p; l++) {
57 double A = ((2 * (double)l + 1) / ((double)l + 1));
58 double B = ((double)l / ((double)l + 1));
59 L[l + 1] = A * s * L[l] - B * L[l - 1];
60 if (diffL != NULL) {
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];
64 if (dim >= 2) {
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];
68 if (dim == 3) {
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];
72 }
73 }
74 }
75 }
77}
78
79PetscErrorCode Jacobi_polynomials(int p, double alpha, double x, double t,
80 double *diff_x, double *diff_t, double *L,
81 double *diffL, const int dim) {
83#ifndef NDEBUG
84 if (dim < 1)
85 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
86 if (dim > 3)
87 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
88 if (p < 0)
89 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
90#endif // NDEBUG
91 L[0] = 1;
92 if (diffL != NULL) {
93 diffL[0 * (p + 1) + 0] = 0;
94 if (dim >= 2) {
95 diffL[1 * (p + 1) + 0] = 0;
96 if (dim == 3) {
97 diffL[2 * (p + 1) + 0] = 0;
98 }
99 }
100 }
101 if (p == 0)
103 L[1] = 2 * x - t + alpha * x;
104 if (diffL != NULL) {
105#ifndef NDEBUG
106 if (diff_x == NULL) {
107 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
108 }
109#endif // NDEBUG
110 double d_t = (diff_t) ? diff_t[0] : 0;
111 diffL[0 * (p + 1) + 1] = (2 + alpha) * diff_x[0] - d_t;
112 if (dim >= 2) {
113 double d_t = (diff_t) ? diff_t[1] : 0;
114 diffL[1 * (p + 1) + 1] = (2 + alpha) * diff_x[1] - d_t;
115 if (dim == 3) {
116 double d_t = (diff_t) ? diff_t[2] : 0;
117 diffL[2 * (p + 1) + 1] = (2 + alpha) * diff_x[2] - d_t;
118 }
119 }
120 }
121 if (p == 1)
123 int l = 1;
124 for (; l < p; l++) {
125 int lp1 = l + 1;
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];
133 if (diffL != NULL) {
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] -
139 diffB * L[l - 1];
140 if (dim >= 2) {
141 double d_t = (diff_t) ? diff_t[1] : 0;
142 double diffA =
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];
148 if (dim == 3) {
149 double d_t = (diff_t) ? diff_t[2] : 0;
150 double diffA =
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];
156 }
157 }
158 }
159 }
161}
162
163PetscErrorCode IntegratedJacobi_polynomials(int p, double alpha, double x,
164 double t, double *diff_x,
165 double *diff_t, double *L,
166 double *diffL, const int dim) {
168#ifndef NDEBUG
169 if (dim < 1)
170 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
171 if (dim > 3)
172 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
173 if (p < 1)
174 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 1");
175#endif // NDEBUG
176 L[0] = x;
177 if (diffL != NULL) {
178 int d = 0;
179 for (; d != dim; ++d) {
180 diffL[d * p + 0] = diff_x[d];
181 }
182 }
183 if (p == 0)
185 double jacobi[(p + 1)];
186 double diff_jacobi[(p + 1) * dim];
187 ierr = Jacobi_polynomials(p, alpha, x, t, diff_x, diff_t, jacobi, diff_jacobi,
188 dim);
189 CHKERRQ(ierr);
190 int l = 1;
191 for (; l < p; l++) {
192 int i = l + 1;
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];
197 if (diffL != NULL) {
198 int dd = 0;
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]);
205 }
206 }
207 }
209}
210
211PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L,
212 double *diffL, const int dim) {
213
215#ifndef NDEBUG
216 if (dim < 1)
217 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
218 if (dim > 3)
219 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
220 if (p < 2)
221 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 2");
222#endif // NDEBUG
223 double l[p + 5];
224 ierr = Legendre_polynomials(p + 4, s, NULL, l, NULL, 1);
225 CHKERRQ(ierr);
226
227 L[0] = 1;
228 if (diffL != NULL) {
229 diffL[0 * (p + 1) + 0] = 0;
230 if (dim >= 2) {
231 diffL[1 * (p + 1) + 0] = 0;
232 if (dim == 3) {
233 diffL[2 * (p + 1) + 0] = 0;
234 }
235 }
236 }
237 L[1] = s;
238 if (diffL != NULL) {
239#ifndef NDEBUG
240 if (diff_s == NULL) {
241 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
242 }
243#endif // NDEBUG
244 diffL[0 * (p + 1) + 1] = diff_s[0];
245 if (dim >= 2) {
246 diffL[1 * (p + 1) + 1] = diff_s[1];
247 if (dim == 3) {
248 diffL[2 * (p + 1) + 1] = diff_s[2];
249 }
250 }
251 }
252
253 // Integrated Legendre
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]);
257
258 if (diffL != NULL) {
259 if (diff_s == NULL) {
260 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
261 }
262 double a = l[k - 1] / 2.;
263 diffL[0 * (p + 1) + k] = a * diff_s[0];
264 if (dim >= 2) {
265 diffL[1 * (p + 1) + k] = a * diff_s[1];
266 if (dim == 3) {
267 diffL[2 * (p + 1) + k] = a * diff_s[2];
268 }
269 }
270 }
271 }
273}
274
275static double f_phi0(double x) { return LOBATTO_PHI0(x); }
276static double f_phi1(double x) { return LOBATTO_PHI1(x); }
277static double f_phi2(double x) { return LOBATTO_PHI2(x); }
278static double f_phi3(double x) { return LOBATTO_PHI3(x); }
279static double f_phi4(double x) { return LOBATTO_PHI4(x); }
280static double f_phi5(double x) { return LOBATTO_PHI5(x); }
281static double f_phi6(double x) { return LOBATTO_PHI6(x); }
282static double f_phi7(double x) { return LOBATTO_PHI7(x); }
283static double f_phi8(double x) { return LOBATTO_PHI8(x); }
284static double f_phi9(double x) { return LOBATTO_PHI9(x); }
285
286static double (*f_phi[])(double x) = {f_phi0, f_phi1, f_phi2, f_phi3, f_phi4,
288
289static double f_phi0x(double x) { return LOBATTO_PHI0X(x); }
290static double f_phi1x(double x) { return LOBATTO_PHI1X(x); }
291static double f_phi2x(double x) { return LOBATTO_PHI2X(x); }
292static double f_phi3x(double x) { return LOBATTO_PHI3X(x); }
293static double f_phi4x(double x) { return LOBATTO_PHI4X(x); }
294static double f_phi5x(double x) { return LOBATTO_PHI5X(x); }
295static double f_phi6x(double x) { return LOBATTO_PHI6X(x); }
296static double f_phi7x(double x) { return LOBATTO_PHI7X(x); }
297static double f_phi8x(double x) { return LOBATTO_PHI8X(x); }
298static double f_phi9x(double x) { return LOBATTO_PHI9X(x); }
299
300static double (*f_phix[])(double x) = {f_phi0x, f_phi1x, f_phi2x, f_phi3x,
303
304// /// Legendre polynomials
305// #define Legendre0(x) (1.0)
306// #define Legendre1(x) (x)
307// #define Legendre2(x) (1.0 / 2.0 * (3 * (x) * (x) - 1))
308// #define Legendre3(x) (1.0 / 2.0 * (5 * (x) * (x) - 3) * (x))
309// #define Legendre4(x) (1.0 / 8.0 * ((35 * (x) * (x) - 30) * (x) * (x) + 3))
310// #define Legendre5(x) (1.0 / 8.0 * ((63 * (x) * (x) - 70) * (x) * (x) + 15) *
311// (x)) #define Legendre6(x) (1.0 / 16.0 * (((231 * (x) * (x) - 315) * (x) * (x)
312// + 105) * (x) * (x) - 5)) #define Legendre7(x) (1.0 / 16.0 * (((429 * (x) *
313// (x) - 693) * (x) * (x) + 315) * (x) * (x) - 35) * (x)) #define Legendre8(x)
314// (1.0 / 128.0 * ((((6435 * (x) * (x) - 12012) * (x) * (x) + 6930) * (x) * (x)
315// - 1260) * (x) * (x) + 35)) #define Legendre9(x) (1.0 / 128.0 * ((((12155 *
316// (x) * (x) - 25740) * (x) * (x) + 18018) * (x) * (x) - 4620) * (x) * (x) +
317// 315) * (x)) #define Legendre10(x) (1.0 / 256.0 * (((((46189 * (x) * (x) -
318// 109395) * (x) * (x) + 90090) * (x) * (x) - 30030) * (x) * (x) + 3465) * (x) *
319// (x) - 63))
320//
321// /// derivatives of Legendre polynomials
322// #define Legendre0x(x) (0.0)
323// #define Legendre1x(x) (1.0)
324// #define Legendre2x(x) (3.0 * (x))
325// #define Legendre3x(x) (15.0 / 2.0 * (x) * (x) - 3.0 / 2.0)
326// #define Legendre4x(x) (5.0 / 2.0 * (x) * (7.0 * (x) * (x) - 3.0))
327// #define Legendre5x(x) ((315.0 / 8.0 * (x) * (x) - 105.0 / 4.0) * (x) * (x)
328// + 15.0 / 8.0) #define Legendre6x(x) (21.0 / 8.0 * (x) * ((33.0 * (x) * (x)
329// - 30.0) * (x) * (x) + 5.0)) #define Legendre7x(x) (((3003.0 / 16.0 * (x) *
330// (x) - 3465.0 / 16.0) * (x) * (x) + 945.0 / 16.0) * (x) * (x) - 35.0 / 16.0)
331// #define Legendre8x(x) (9.0 / 16.0 * (x) * (((715.0 * (x) * (x) - 1001.0) *
332// (x) * (x) + 385.0) * (x) * (x) - 35.0)) #define Legendre9x(x) ((((109395.0 /
333// 128.0 * (x) * (x) - 45045.0 / 32.0) * (x) * (x) + 45045.0 / 64.0) * (x) * (x)
334// - 3465.0 / 32.0) * (x) * (x) + 315.0 / 128.0) #define Legendre10x(x) (2.0 /
335// 256.0 * (x) * ((((230945.0 * (x) * (x) - 437580.0) * (x) * (x) + 270270.0) *
336// (x) * (x) - 60060.0) * (x) * (x) + 3465.0))
337//
338// /// first two Lobatto shape functions
339// #define l0(x) ((1.0 - (x)) * 0.5)
340// #define l1(x) ((1.0 + (x)) * 0.5)
341//
342// #define l0l1(x) ((1.0 - (x)*(x)) * 0.25)
343//
344// /// other Lobatto shape functions
345// #define l2(x) (phi0(x) * l0l1(x))
346// #define l3(x) (phi1(x) * l0l1(x))
347// #define l4(x) (phi2(x) * l0l1(x))
348// #define l5(x) (phi3(x) * l0l1(x))
349// #define l6(x) (phi4(x) * l0l1(x))
350// #define l7(x) (phi5(x) * l0l1(x))
351// #define l8(x) (phi6(x) * l0l1(x))
352// #define l9(x) (phi7(x) * l0l1(x))
353// #define l10(x) (phi8(x) * l0l1(x))
354// #define l11(x) (phi9(x) * l0l1(x))
355//
356// /// derivatives of Lobatto functions
357// #define dl0(x) (-0.5)
358// #define dl1(x) (0.5)
359// #define dl2(x) (sqrt(3.0/2.0) * Legendre1(x))
360// #define dl3(x) (sqrt(5.0/2.0) * Legendre2(x))
361// #define dl4(x) (sqrt(7.0/2.0) * Legendre3(x))
362// #define dl5(x) (sqrt(9.0/2.0) * Legendre4(x))
363// #define dl6(x) (sqrt(11.0/2.0) * Legendre5(x))
364// #define dl7(x) (sqrt(13.0/2.0) * Legendre6(x))
365// #define dl8(x) (sqrt(15.0/2.0) * Legendre7(x))
366// #define dl9(x) (sqrt(17.0/2.0) * Legendre8(x))
367// #define dl10(x) (sqrt(19.0/2.0) * Legendre9(x))
368// #define dl11(x) (sqrt(21.0/2.0) * Legendre10(x))
369
370PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s,
371 double *L, double *diffL,
372 const int dim) {
374#ifndef NDEBUG
375 if (dim < 1)
376 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
377 if (dim > 3)
378 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
379 if (p < 0)
380 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
381 if (p > 9)
382 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
383 "Polynomial beyond order 9 is not implemented");
384#endif // NDEBUG
385 if (L) {
386 int l = 0;
387 for (; l != p + 1; l++) {
388 L[l] = f_phi[l](s);
389 }
390 }
391 if (diffL != NULL) {
392#ifndef NDEBUG
393 if (diff_s == NULL) {
394 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
395 }
396#endif // NDEBUG
397 int l = 0;
398 for (; l != p + 1; l++) {
399 double a = f_phix[l](s);
400 diffL[0 * (p + 1) + l] = diff_s[0] * a;
401 if (dim >= 2) {
402 diffL[1 * (p + 1) + l] = diff_s[1] * a;
403 if (dim == 3) {
404 diffL[2 * (p + 1) + l] = diff_s[2] * a;
405 }
406 }
407 }
408 }
410}
static Index< 'L', 3 > L
static Index< 'p', 3 > p
constexpr double 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_PHI2(x)
#define LOBATTO_PHI5(x)
#define LOBATTO_PHI8(x)
#define LOBATTO_PHI3(x)
#define LOBATTO_PHI7X(x)
#define LOBATTO_PHI0X(x)
Derivatives of kernel functions for Lobbatto base.
#define LOBATTO_PHI2X(x)
#define LOBATTO_PHI5X(x)
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
#define LOBATTO_PHI1(x)
#define LOBATTO_PHI9X(x)
#define LOBATTO_PHI4(x)
#define LOBATTO_PHI8X(x)
#define LOBATTO_PHI7(x)
#define LOBATTO_PHI6(x)
#define LOBATTO_PHI1X(x)
#define LOBATTO_PHI9(x)
#define LOBATTO_PHI4X(x)
#define LOBATTO_PHI6X(x)
#define LOBATTO_PHI3X(x)
useful compiler directives and definitions
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
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 AssemblyType A
Definition: plastic.cpp:35
constexpr double t
plate stiffness
Definition: plate.cpp:59