v0.9.0
base_functions.c
Go to the documentation of this file.
1 /** \file base_functions.c
2 
3 */
4 
5 /**
6  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
7  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
8  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
9  * License for more details.
10  *
11  * You should have received a copy of the GNU Lesser General Public
12  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
13  */
14 
15 #include <cblas.h>
16 #include <petscsys.h>
17 #include <phg-quadrule/quad.h>
18 
19 #include <definitions.h>
20 
21 #include <base_functions.h>
22 
23 static PetscErrorCode ierr;
24 
25 PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L,
26  double *diffL, const int dim) {
28  if (dim < 1)
29  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
30  if (dim > 3)
31  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
32  if (p < 0)
33  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
34  L[0] = 1;
35  if (diffL != NULL) {
36  diffL[0 * (p + 1) + 0] = 0;
37  if (dim >= 2) {
38  diffL[1 * (p + 1) + 0] = 0;
39  if (dim == 3) {
40  diffL[2 * (p + 1) + 0] = 0;
41  }
42  }
43  }
44  if (p == 0)
46  L[1] = s;
47  if (diffL != NULL) {
48  if (diff_s == NULL) {
49  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
50  }
51  diffL[0 * (p + 1) + 1] = diff_s[0];
52  if (dim >= 2) {
53  diffL[1 * (p + 1) + 1] = diff_s[1];
54  if (dim == 3) {
55  diffL[2 * (p + 1) + 1] = diff_s[2];
56  }
57  }
58  }
59  if (p == 1)
61  int l = 1;
62  for (; l < p; l++) {
63  double A = ((2 * (double)l + 1) / ((double)l + 1));
64  double B = ((double)l / ((double)l + 1));
65  L[l + 1] = A * s * L[l] - B * L[l - 1];
66  if (diffL != NULL) {
67  diffL[0 * (p + 1) + l + 1] =
68  A * (s * diffL[0 * (p + 1) + l] + diff_s[0] * L[l]) -
69  B * diffL[0 * (p + 1) + l - 1];
70  if (dim >= 2) {
71  diffL[1 * (p + 1) + l + 1] =
72  A * (s * diffL[1 * (p + 1) + l] + diff_s[1] * L[l]) -
73  B * diffL[1 * (p + 1) + l - 1];
74  if (dim == 3) {
75  diffL[2 * (p + 1) + l + 1] =
76  A * (s * diffL[2 * (p + 1) + l] + diff_s[2] * L[l]) -
77  B * diffL[2 * (p + 1) + l - 1];
78  }
79  }
80  }
81  }
83 }
84 
85 PetscErrorCode Jacobi_polynomials(int p, double alpha, double x, double t,
86  double *diff_x, double *diff_t, double *L,
87  double *diffL, const int dim) {
89  if (dim < 1)
90  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
91  if (dim > 3)
92  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
93  if (p < 0)
94  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
95  L[0] = 1;
96  if (diffL != NULL) {
97  diffL[0 * (p + 1) + 0] = 0;
98  if (dim >= 2) {
99  diffL[1 * (p + 1) + 0] = 0;
100  if (dim == 3) {
101  diffL[2 * (p + 1) + 0] = 0;
102  }
103  }
104  }
105  if (p == 0)
107  L[1] = 2 * x - t + alpha * x;
108  if (diffL != NULL) {
109  if (diff_x == NULL) {
110  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
111  }
112  double d_t = (diff_t) ? diff_t[0] : 0;
113  diffL[0 * (p + 1) + 1] = (2 + alpha) * diff_x[0] - d_t;
114  if (dim >= 2) {
115  double d_t = (diff_t) ? diff_t[1] : 0;
116  diffL[1 * (p + 1) + 1] = (2 + alpha) * diff_x[1] - d_t;
117  if (dim == 3) {
118  double d_t = (diff_t) ? diff_t[2] : 0;
119  diffL[2 * (p + 1) + 1] = (2 + alpha) * diff_x[2] - d_t;
120  }
121  }
122  }
123  if (p == 1)
125  int l = 1;
126  for (; l < p; l++) {
127  int lp1 = l + 1;
128  double a = 2 * lp1 * (lp1 + alpha) * (2 * lp1 + alpha - 2);
129  double b = 2 * lp1 + alpha - 1;
130  double c = (2 * lp1 + alpha) * (2 * lp1 + alpha - 2);
131  double d = 2 * (lp1 + alpha - 1) * (lp1 - 1) * (2 * lp1 + alpha);
132  double A = b * (c * (2 * x - t) + alpha * alpha * t) / a;
133  double B = d * t * t / a;
134  L[lp1] = A * L[l] - B * L[l - 1];
135  if (diffL != NULL) {
136  double d_t = (diff_t) ? diff_t[0] : 0;
137  double diffA = b * (c * (2 * diff_x[0] - d_t) + alpha * alpha * d_t) / a;
138  double diffB = d * 2 * t * d_t / a;
139  diffL[0 * (p + 1) + lp1] = A * diffL[0 * (p + 1) + l] -
140  B * diffL[0 * (p + 1) + l - 1] + diffA * L[l] -
141  diffB * L[l - 1];
142  if (dim >= 2) {
143  double d_t = (diff_t) ? diff_t[1] : 0;
144  double diffA =
145  b * (c * (2 * diff_x[1] - d_t) + alpha * alpha * d_t) / a;
146  double diffB = d * 2 * t * d_t / a;
147  diffL[1 * (p + 1) + lp1] = A * diffL[1 * (p + 1) + l] -
148  B * diffL[1 * (p + 1) + l - 1] +
149  diffA * L[l] - diffB * L[l - 1];
150  if (dim == 3) {
151  double d_t = (diff_t) ? diff_t[2] : 0;
152  double diffA =
153  b * (c * (2 * diff_x[2] - d_t) + alpha * alpha * d_t) / a;
154  double diffB = d * 2 * t * d_t / a;
155  diffL[2 * (p + 1) + lp1] = A * diffL[2 * (p + 1) + l] -
156  B * diffL[2 * (p + 1) + l - 1] +
157  diffA * L[l] - diffB * L[l - 1];
158  }
159  }
160  }
161  }
163 }
164 
165 PetscErrorCode IntegratedJacobi_polynomials(int p, double alpha, double x,
166  double t, double *diff_x,
167  double *diff_t, double *L,
168  double *diffL, const int dim) {
170  if (dim < 1)
171  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
172  if (dim > 3)
173  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
174  if (p < 1)
175  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 1");
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 
211 PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L,
212  double *diffL, const int dim) {
213 
215  if (dim < 1)
216  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
217  if (dim > 3)
218  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
219  if (p < 2)
220  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 2");
221  p -= 2; // polynomial order starts from 2
222  double l[p + 2];
223  ierr = Legendre_polynomials(p + 1, s, NULL, l, NULL, 1);
224  CHKERRQ(ierr);
225  {
226  // Derivatives
227  int k = 0;
228  for (; k <= p; k++) {
229  if (diffL != NULL) {
230  if (diff_s == NULL) {
231  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
232  }
233  double a = l[k + 1];
234  diffL[0 * (p + 1) + k] = a * diff_s[0];
235  if (dim >= 2) {
236  diffL[1 * (p + 1) + k] = a * diff_s[1];
237  if (dim == 3) {
238  diffL[2 * (p + 1) + k] = a * diff_s[2];
239  }
240  }
241  }
242  }
243  }
244  {
245  // Functions
246  bzero(L, (p + 1) * sizeof(double));
247  int nb_gauss_pts = QUAD_1D_TABLE[p + 2]->npoints;
248  double *points = QUAD_1D_TABLE[p + 2]->points;
249  double *weights = QUAD_1D_TABLE[p + 2]->weights;
250  s = s + 1;
251  int gg = 0;
252  for (; gg != nb_gauss_pts; gg++) {
253  double ksi = points[2 * gg + 1];
254  double zeta = s * ksi - 1;
255  ierr = Legendre_polynomials(p + 1, zeta, NULL, l, NULL, 1);
256  CHKERRQ(ierr);
257  double w = s * weights[gg];
258  cblas_daxpy(p + 1, w, &l[1], 1, &L[0], 1);
259  }
260  }
261  {
262  int k = 0;
263  for (; k <= p; k++) {
264  double a = 4 * sqrt(k + 2 - 0.5);
265  if (L != NULL)
266  L[k] *= a;
267  if (diffL != NULL)
268  diffL[k] *= a;
269  }
270  }
272 }
273 
274 static double f_phi0(double x) { return LOBATTO_PHI0(x); }
275 static double f_phi1(double x) { return LOBATTO_PHI1(x); }
276 static double f_phi2(double x) { return LOBATTO_PHI2(x); }
277 static double f_phi3(double x) { return LOBATTO_PHI3(x); }
278 static double f_phi4(double x) { return LOBATTO_PHI4(x); }
279 static double f_phi5(double x) { return LOBATTO_PHI5(x); }
280 static double f_phi6(double x) { return LOBATTO_PHI6(x); }
281 static double f_phi7(double x) { return LOBATTO_PHI7(x); }
282 static double f_phi8(double x) { return LOBATTO_PHI8(x); }
283 static double f_phi9(double x) { return LOBATTO_PHI9(x); }
284 
285 static double (*f_phi[])(double x) = {f_phi0, f_phi1, f_phi2, f_phi3, f_phi4,
287 
288 static double f_phi0x(double x) { return LOBATTO_PHI0X(x); }
289 static double f_phi1x(double x) { return LOBATTO_PHI1X(x); }
290 static double f_phi2x(double x) { return LOBATTO_PHI2X(x); }
291 static double f_phi3x(double x) { return LOBATTO_PHI3X(x); }
292 static double f_phi4x(double x) { return LOBATTO_PHI4X(x); }
293 static double f_phi5x(double x) { return LOBATTO_PHI5X(x); }
294 static double f_phi6x(double x) { return LOBATTO_PHI6X(x); }
295 static double f_phi7x(double x) { return LOBATTO_PHI7X(x); }
296 static double f_phi8x(double x) { return LOBATTO_PHI8X(x); }
297 static double f_phi9x(double x) { return LOBATTO_PHI9X(x); }
298 
299 static double (*f_phix[])(double x) = {f_phi0x, f_phi1x, f_phi2x, f_phi3x,
301  f_phi8x, f_phi9x};
302 
303 // /// Legendre polynomials
304 // #define Legendre0(x) (1.0)
305 // #define Legendre1(x) (x)
306 // #define Legendre2(x) (1.0 / 2.0 * (3 * (x) * (x) - 1))
307 // #define Legendre3(x) (1.0 / 2.0 * (5 * (x) * (x) - 3) * (x))
308 // #define Legendre4(x) (1.0 / 8.0 * ((35 * (x) * (x) - 30) * (x) * (x) + 3))
309 // #define Legendre5(x) (1.0 / 8.0 * ((63 * (x) * (x) - 70) * (x) * (x) + 15) *
310 // (x)) #define Legendre6(x) (1.0 / 16.0 * (((231 * (x) * (x) - 315) * (x) * (x)
311 // + 105) * (x) * (x) - 5)) #define Legendre7(x) (1.0 / 16.0 * (((429 * (x) *
312 // (x) - 693) * (x) * (x) + 315) * (x) * (x) - 35) * (x)) #define Legendre8(x)
313 // (1.0 / 128.0 * ((((6435 * (x) * (x) - 12012) * (x) * (x) + 6930) * (x) * (x)
314 // - 1260) * (x) * (x) + 35)) #define Legendre9(x) (1.0 / 128.0 * ((((12155 *
315 // (x) * (x) - 25740) * (x) * (x) + 18018) * (x) * (x) - 4620) * (x) * (x) +
316 // 315) * (x)) #define Legendre10(x) (1.0 / 256.0 * (((((46189 * (x) * (x) -
317 // 109395) * (x) * (x) + 90090) * (x) * (x) - 30030) * (x) * (x) + 3465) * (x) *
318 // (x) - 63))
319 //
320 // /// derivatives of Legendre polynomials
321 // #define Legendre0x(x) (0.0)
322 // #define Legendre1x(x) (1.0)
323 // #define Legendre2x(x) (3.0 * (x))
324 // #define Legendre3x(x) (15.0 / 2.0 * (x) * (x) - 3.0 / 2.0)
325 // #define Legendre4x(x) (5.0 / 2.0 * (x) * (7.0 * (x) * (x) - 3.0))
326 // #define Legendre5x(x) ((315.0 / 8.0 * (x) * (x) - 105.0 / 4.0) * (x) * (x)
327 // + 15.0 / 8.0) #define Legendre6x(x) (21.0 / 8.0 * (x) * ((33.0 * (x) * (x)
328 // - 30.0) * (x) * (x) + 5.0)) #define Legendre7x(x) (((3003.0 / 16.0 * (x) *
329 // (x) - 3465.0 / 16.0) * (x) * (x) + 945.0 / 16.0) * (x) * (x) - 35.0 / 16.0)
330 // #define Legendre8x(x) (9.0 / 16.0 * (x) * (((715.0 * (x) * (x) - 1001.0) *
331 // (x) * (x) + 385.0) * (x) * (x) - 35.0)) #define Legendre9x(x) ((((109395.0 /
332 // 128.0 * (x) * (x) - 45045.0 / 32.0) * (x) * (x) + 45045.0 / 64.0) * (x) * (x)
333 // - 3465.0 / 32.0) * (x) * (x) + 315.0 / 128.0) #define Legendre10x(x) (2.0 /
334 // 256.0 * (x) * ((((230945.0 * (x) * (x) - 437580.0) * (x) * (x) + 270270.0) *
335 // (x) * (x) - 60060.0) * (x) * (x) + 3465.0))
336 //
337 // /// first two Lobatto shape functions
338 // #define l0(x) ((1.0 - (x)) * 0.5)
339 // #define l1(x) ((1.0 + (x)) * 0.5)
340 //
341 // #define l0l1(x) ((1.0 - (x)*(x)) * 0.25)
342 //
343 // /// other Lobatto shape functions
344 // #define l2(x) (phi0(x) * l0l1(x))
345 // #define l3(x) (phi1(x) * l0l1(x))
346 // #define l4(x) (phi2(x) * l0l1(x))
347 // #define l5(x) (phi3(x) * l0l1(x))
348 // #define l6(x) (phi4(x) * l0l1(x))
349 // #define l7(x) (phi5(x) * l0l1(x))
350 // #define l8(x) (phi6(x) * l0l1(x))
351 // #define l9(x) (phi7(x) * l0l1(x))
352 // #define l10(x) (phi8(x) * l0l1(x))
353 // #define l11(x) (phi9(x) * l0l1(x))
354 //
355 // /// derivatives of Lobatto functions
356 // #define dl0(x) (-0.5)
357 // #define dl1(x) (0.5)
358 // #define dl2(x) (sqrt(3.0/2.0) * Legendre1(x))
359 // #define dl3(x) (sqrt(5.0/2.0) * Legendre2(x))
360 // #define dl4(x) (sqrt(7.0/2.0) * Legendre3(x))
361 // #define dl5(x) (sqrt(9.0/2.0) * Legendre4(x))
362 // #define dl6(x) (sqrt(11.0/2.0) * Legendre5(x))
363 // #define dl7(x) (sqrt(13.0/2.0) * Legendre6(x))
364 // #define dl8(x) (sqrt(15.0/2.0) * Legendre7(x))
365 // #define dl9(x) (sqrt(17.0/2.0) * Legendre8(x))
366 // #define dl10(x) (sqrt(19.0/2.0) * Legendre9(x))
367 // #define dl11(x) (sqrt(21.0/2.0) * Legendre10(x))
368 
369 PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s,
370  double *L, double *diffL,
371  const int dim) {
373  if (dim < 1)
374  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim < 1");
375  if (dim > 3)
376  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "dim > 3");
377  if (p < 0)
378  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "p < 0");
379  if (p > 9)
380  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
381  "Polynomial beyond order 9 is not implemented");
382  if (L) {
383  int l = 0;
384  for (; l != p + 1; l++) {
385  L[l] = f_phi[l](s);
386  }
387  }
388  if (diffL != NULL) {
389  if (diff_s == NULL) {
390  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "diff_s == NULL");
391  }
392  int l = 0;
393  for (; l != p + 1; l++) {
394  double a = f_phix[l](s);
395  diffL[0 * (p + 1) + l] = diff_s[0] * a;
396  if (dim >= 2) {
397  diffL[1 * (p + 1) + l] = diff_s[1] * a;
398  if (dim == 3) {
399  diffL[2 * (p + 1) + l] = diff_s[2] * a;
400  }
401  }
402  }
403  }
405 }
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
static double f_phi6(double x)
static double f_phi8x(double x)
#define LOBATTO_PHI1X(x)
#define LOBATTO_PHI6X(x)
#define LOBATTO_PHI2X(x)
double * points
Definition: quad.h:30
static double f_phi3(double x)
static double f_phi9x(double x)
static double f_phi1x(double x)
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
static double f_phi5(double x)
#define LOBATTO_PHI1(x)
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
static double f_phi2x(double x)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
int npoints
Definition: quad.h:29
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_phi1(double x)
#define LOBATTO_PHI5X(x)
static double f_phi0(double x)
static PetscErrorCode ierr
#define LOBATTO_PHI3X(x)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static double f_phi6x(double x)
static double f_phi0x(double x)
static double f_phi3x(double x)
#define LOBATTO_PHI9(x)
static double f_phi2(double x)
#define LOBATTO_PHI3(x)
static double f_phi7x(double x)
#define LOBATTO_PHI9X(x)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
#define LOBATTO_PHI7(x)
#define LOBATTO_PHI8X(x)
#define LOBATTO_PHI7X(x)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Kernel Lobatto base functions.
#define LOBATTO_PHI8(x)
CHKERRQ(ierr)
useful compiler directives and definitions
static double f_phi5x(double x)
static double f_phi8(double x)
static double f_phi7(double x)
double * weights
Definition: quad.h:31
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
static double f_phi9(double x)
#define LOBATTO_PHI0X(x)
Derivatives of kernel functions for Lobbatto base.
#define LOBATTO_PHI4X(x)
static double f_phi4x(double x)
#define LOBATTO_PHI5(x)
#define LOBATTO_PHI4(x)
PetscErrorCode Lobatto_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Lobatto base functions.
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_phi4(double x)
#define LOBATTO_PHI6(x)
#define LOBATTO_PHI2(x)