Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
poly.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <stdarg.h>
4 #include <math.h> /* for cos, fabs */
5 #include <string.h> /* for memcpy */
6 #include <float.h>
7 
8 #include "moab/FindPtFuncs.h"
9 
10 /*
11  For brevity's sake, some names have been shortened
12  Quadrature rules
13  Gauss -> Gauss-Legendre quadrature (open)
14  Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends)
15  Polynomial bases
16  Legendre -> Legendre basis
17  Gauss -> Lagrangian basis using Gauss quadrature nodes
18  Lobatto -> Lagrangian basis using Lobatto quadrature nodes
19 */
20 
21 /*--------------------------------------------------------------------------
22  Legendre Polynomial Matrix Computation
23  (compute P_i(x_j) for i = 0, ..., n and a given set of x)
24  --------------------------------------------------------------------------*/
25 
26 /* precondition: n >= 1
27  inner index is x index (0 ... m-1);
28  outer index is Legendre polynomial number (0 ... n)
29  */
30 void legendre_matrix( const realType* x, int m, realType* P, int n )
31 {
32  int i, j;
33  realType *Pjm1 = P, *Pj = Pjm1 + m, *Pjp1 = Pj + m;
34  for( i = 0; i < m; ++i )
35  Pjm1[i] = 1;
36  for( i = 0; i < m; ++i )
37  Pj[i] = x[i];
38  for( j = 1; j < n; ++j )
39  {
40  realType c = 1 / (realType)( j + 1 ), a = c * ( 2 * j + 1 ), b = c * j;
41  for( i = 0; i < m; ++i )
42  Pjp1[i] = a * x[i] * Pj[i] - b * Pjm1[i];
43  Pjp1 += m, Pj += m, Pjm1 += m;
44  }
45 }
46 
47 /* precondition: n >= 1, n even */
48 static void legendre_row_even( realType x, realType* P, int n )
49 {
50  int i;
51  P[0] = 1, P[1] = x;
52  for( i = 1; i <= n - 2; i += 2 )
53  {
54  P[i + 1] = ( ( 2 * i + 1 ) * x * P[i] - i * P[i - 1] ) / ( i + 1 );
55  P[i + 2] = ( ( 2 * i + 3 ) * x * P[i - 1] - ( i + 1 ) * P[i] ) / ( i + 2 );
56  }
57  P[n] = ( ( 2 * n - 1 ) * x * P[n - 1] - ( n - 1 ) * P[n - 2] ) / n;
58 }
59 
60 /* precondition: n >= 1, n odd */
61 static void legendre_row_odd( realType x, realType* P, int n )
62 {
63  int i;
64  P[0] = 1, P[1] = x;
65  for( i = 1; i <= n - 2; i += 2 )
66  {
67  P[i + 1] = ( ( 2 * i + 1 ) * x * P[i] - i * P[i - 1] ) / ( i + 1 );
68  P[i + 2] = ( ( 2 * i + 3 ) * x * P[i - 1] - ( i + 1 ) * P[i] ) / ( i + 2 );
69  }
70 }
71 
72 /* precondition: n >= 1
73  compute P_i(x) with i = 0 ... n
74  */
75 void legendre_row( realType x, realType* P, int n )
76 {
77  if( n & 1 )
78  legendre_row_odd( x, P, n );
79  else
80  legendre_row_even( x, P, n );
81 }
82 
83 /* precondition: n >= 1
84  inner index is Legendre polynomial number (0 ... n)
85  outer index is x index (0 ... m-1);
86  */
87 void legendre_matrix_t( const realType* x, int m, realType* P, int n )
88 {
89  int i;
90  if( n & 1 )
91  for( i = 0; i < m; ++i, P += n + 1 )
92  legendre_row_odd( x[i], P, n );
93  else
94  for( i = 0; i < m; ++i, P += n + 1 )
95  legendre_row_even( x[i], P, n );
96 }
97 
98 /*--------------------------------------------------------------------------
99  Legendre Polynomial Computation
100  compute P_n(x) or P_n'(x) or P_n''(x)
101  --------------------------------------------------------------------------*/
102 
103 /* precondition: n >= 0 */
104 static realType legendre( int n, realType x )
105 {
106  realType p[2];
107  int i;
108  p[0] = 1;
109  p[1] = x;
110  for( i = 1; i < n; i += 2 )
111  {
112  p[0] = ( ( 2 * i + 1 ) * x * p[1] - i * p[0] ) / ( i + 1 );
113  p[1] = ( ( 2 * i + 3 ) * x * p[0] - ( i + 1 ) * p[1] ) / ( i + 2 );
114  }
115  return p[n & 1];
116 }
117 
118 /* precondition: n > 0 */
119 static realType legendre_d1( int n, realType x )
120 {
121  realType p[2];
122  int i;
123  p[0] = 3 * x;
124  p[1] = 1;
125  for( i = 2; i < n; i += 2 )
126  {
127  p[1] = ( ( 2 * i + 1 ) * x * p[0] - ( i + 1 ) * p[1] ) / i;
128  p[0] = ( ( 2 * i + 3 ) * x * p[1] - ( i + 2 ) * p[0] ) / ( i + 1 );
129  }
130  return p[n & 1];
131 }
132 
133 /* precondition: n > 1 */
134 static realType legendre_d2( int n, realType x )
135 {
136  realType p[2];
137  int i;
138  p[0] = 3;
139  p[1] = 15 * x;
140  for( i = 3; i < n; i += 2 )
141  {
142  p[0] = ( ( 2 * i + 1 ) * x * p[1] - ( i + 2 ) * p[0] ) / ( i - 1 );
143  p[1] = ( ( 2 * i + 3 ) * x * p[0] - ( i + 3 ) * p[1] ) / i;
144  }
145  return p[n & 1];
146 }
147 
148 /*--------------------------------------------------------------------------
149  Quadrature Nodes and Weights Calculation
150  compute the n Gauss-Legendre nodes and weights or
151  the n Gauss-Lobatto-Legendre nodes and weights
152  --------------------------------------------------------------------------*/
153 
154 /* n nodes */
155 void gauss_nodes( realType* z, int n )
156 {
157  int i, j;
158  for( i = 0; i <= n / 2 - 1; ++i )
159  {
160  realType ox, x = mbcos( ( 2 * n - 2 * i - 1 ) * ( MOAB_POLY_PI / 2 ) / n );
161  do
162  {
163  ox = x;
164  x -= legendre( n, x ) / legendre_d1( n, x );
165  } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS );
166  z[i] = x - legendre( n, x ) / legendre_d1( n, x );
167  }
168  if( n & 1 ) z[n / 2] = 0;
169  for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
170  z[j] = -z[i];
171 }
172 
173 /* n inner lobatto nodes (excluding -1,1) */
174 static void lobatto_nodes_aux( realType* z, int n )
175 {
176  int i, j, np = n + 1;
177  for( i = 0; i <= n / 2 - 1; ++i )
178  {
179  realType ox, x = mbcos( ( n - i ) * MOAB_POLY_PI / np );
180  do
181  {
182  ox = x;
183  x -= legendre_d1( np, x ) / legendre_d2( np, x );
184  } while( mbabs( x - ox ) > -x * MOAB_POLY_EPS );
185  z[i] = x - legendre_d1( np, x ) / legendre_d2( np, x );
186  }
187  if( n & 1 ) z[n / 2] = 0;
188  for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
189  z[j] = -z[i];
190 }
191 
192 /* n lobatto nodes */
193 void lobatto_nodes( realType* z, int n )
194 {
195  z[0] = -1, z[n - 1] = 1;
196  lobatto_nodes_aux( &z[1], n - 2 );
197 }
198 
199 void gauss_weights( const realType* z, realType* w, int n )
200 {
201  int i, j;
202  for( i = 0; i <= ( n - 1 ) / 2; ++i )
203  {
204  realType d = ( n + 1 ) * legendre( n + 1, z[i] );
205  w[i] = 2 * ( 1 - z[i] * z[i] ) / ( d * d );
206  }
207  for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
208  w[j] = w[i];
209 }
210 
211 void lobatto_weights( const realType* z, realType* w, int n )
212 {
213  int i, j;
214  for( i = 0; i <= ( n - 1 ) / 2; ++i )
215  {
216  realType d = legendre( n - 1, z[i] );
217  w[i] = 2 / ( ( n - 1 ) * n * d * d );
218  }
219  for( j = ( n + 1 ) / 2, i = n / 2 - 1; j < n; ++j, --i )
220  w[j] = w[i];
221 }
222 
223 /*--------------------------------------------------------------------------
224  Lagrangian to Legendre Change-of-basis Matrix
225  where the nodes of the Lagrangian basis are the GL or GLL quadrature nodes
226  --------------------------------------------------------------------------*/
227 
228 /* precondition: n >= 2
229  given the Gauss quadrature rule (z,w,n), compute the square matrix J
230  for transforming from the Gauss basis to the Legendre basis:
231 
232  u_legendre(i) = sum_j J(i,j) u_gauss(j)
233 
234  computes J = .5 (2i+1) w P (z )
235  ij j i j
236 
237  in column major format (inner index is i, the Legendre index)
238  */
239 void gauss_to_legendre( const realType* z, const realType* w, int n, realType* J )
240 {
241  int i, j;
242  legendre_matrix_t( z, n, J, n - 1 );
243  for( j = 0; j < n; ++j )
244  {
245  realType ww = w[j];
246  for( i = 0; i < n; ++i )
247  *J++ *= ( 2 * i + 1 ) * ww / 2;
248  }
249 }
250 
251 /* precondition: n >= 2
252  same as above, but
253  in row major format (inner index is j, the Gauss index)
254  */
255 void gauss_to_legendre_t( const realType* z, const realType* w, int n, realType* J )
256 {
257  int i, j;
258  legendre_matrix( z, n, J, n - 1 );
259  for( i = 0; i < n; ++i )
260  {
261  realType ii = (realType)( 2 * i + 1 ) / 2;
262  for( j = 0; j < n; ++j )
263  *J++ *= ii * w[j];
264  }
265 }
266 
267 /* precondition: n >= 3
268  given the Lobatto quadrature rule (z,w,n), compute the square matrix J
269  for transforming from the Gauss basis to the Legendre basis:
270 
271  u_legendre(i) = sum_j J(i,j) u_lobatto(j)
272 
273  in column major format (inner index is i, the Legendre index)
274  */
275 void lobatto_to_legendre( const realType* z, const realType* w, int n, realType* J )
276 {
277  int i, j, m = ( n + 1 ) / 2;
278  realType *p = J, *q;
279  realType ww, sum;
280  if( n & 1 )
281  for( j = 0; j < m; ++j, p += n )
282  legendre_row_odd( z[j], p, n - 2 );
283  else
284  for( j = 0; j < m; ++j, p += n )
285  legendre_row_even( z[j], p, n - 2 );
286  p = J;
287  for( j = 0; j < m; ++j )
288  {
289  ww = w[j], sum = 0;
290  for( i = 0; i < n - 1; ++i )
291  *p *= ( 2 * i + 1 ) * ww / 2, sum += *p++;
292  *p++ = -sum;
293  }
294  q = J + ( n / 2 - 1 ) * n;
295  if( n & 1 )
296  for( ; j < n; ++j, p += n, q -= n )
297  {
298  for( i = 0; i < n - 1; i += 2 )
299  p[i] = q[i], p[i + 1] = -q[i + 1];
300  p[i] = q[i];
301  }
302  else
303  for( ; j < n; ++j, p += n, q -= n )
304  {
305  for( i = 0; i < n - 1; i += 2 )
306  p[i] = q[i], p[i + 1] = -q[i + 1];
307  }
308 }
309 
310 /*--------------------------------------------------------------------------
311  Lagrangian to Lagrangian change-of-basis matrix, and derivative matrix
312  --------------------------------------------------------------------------*/
313 
314 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
315  evaluate all Lagrangian basis functions at all points x
316 
317  inner index of output J is the basis function index (row-major format)
318  provide work array with space for 4*n reals
319  */
320 void lagrange_weights( const realType* z, unsigned n, const realType* x, unsigned m, realType* J, realType* work )
321 {
322  unsigned i, j;
323  realType *w = work, *d = w + n, *u = d + n, *v = u + n;
324  for( i = 0; i < n; ++i )
325  {
326  realType ww = 1, zi = z[i];
327  for( j = 0; j < i; ++j )
328  ww *= zi - z[j];
329  for( ++j; j < n; ++j )
330  ww *= zi - z[j];
331  w[i] = 1 / ww;
332  }
333  u[0] = v[n - 1] = 1;
334  for( i = 0; i < m; ++i )
335  {
336  realType xi = x[i];
337  for( j = 0; j < n; ++j )
338  d[j] = xi - z[j];
339  for( j = 0; j < n - 1; ++j )
340  u[j + 1] = d[j] * u[j];
341  for( j = n - 1; j; --j )
342  v[j - 1] = d[j] * v[j];
343  for( j = 0; j < n; ++j )
344  *J++ = w[j] * u[j] * v[j];
345  }
346 }
347 
348 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
349  evaluate all Lagrangian basis functions and their derivatives
350 
351  inner index of outputs J,D is the basis function index (row-major format)
352  provide work array with space for 6*n reals
353  */
354 void lagrange_weights_deriv( const realType* z,
355  unsigned n,
356  const realType* x,
357  unsigned m,
358  realType* J,
359  realType* D,
360  realType* work )
361 {
362  unsigned i, j;
363  realType *w = work, *d = w + n, *u = d + n, *v = u + n, *up = v + n, *vp = up + n;
364  for( i = 0; i < n; ++i )
365  {
366  realType ww = 1, zi = z[i];
367  for( j = 0; j < i; ++j )
368  ww *= zi - z[j];
369  for( ++j; j < n; ++j )
370  ww *= zi - z[j];
371  w[i] = 1 / ww;
372  }
373  u[0] = v[n - 1] = 1;
374  up[0] = vp[n - 1] = 0;
375  for( i = 0; i < m; ++i )
376  {
377  realType xi = x[i];
378  for( j = 0; j < n; ++j )
379  d[j] = xi - z[j];
380  for( j = 0; j < n - 1; ++j )
381  u[j + 1] = d[j] * u[j], up[j + 1] = d[j] * up[j] + u[j];
382  for( j = n - 1; j; --j )
383  v[j - 1] = d[j] * v[j], vp[j - 1] = d[j] * vp[j] + v[j];
384  for( j = 0; j < n; ++j )
385  *J++ = w[j] * u[j] * v[j], *D++ = w[j] * ( up[j] * v[j] + u[j] * vp[j] );
386  }
387 }
388 
389 /*--------------------------------------------------------------------------
390  Speedy Lagrangian Interpolation
391 
392  Usage:
393 
394  lagrange_data p;
395  lagrange_setup(&p,z,n); // setup for nodes z[0 ... n-1]
396 
397  the weights
398  p->J [0 ... n-1] interpolation weights
399  p->D [0 ... n-1] 1st derivative weights
400  p->D2[0 ... n-1] 2nd derivative weights
401  are computed for a given x with:
402  lagrange_0(p,x); // compute p->J
403  lagrange_1(p,x); // compute p->J, p->D
404  lagrange_2(p,x); // compute p->J, p->D, p->D2
405  lagrange_2u(p); // compute p->D2 after call of lagrange_1(p,x);
406  These functions use the z array supplied to setup
407  (that pointer should not be freed between calls)
408  Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
409  p->J_z0, etc. and p->J_zn, etc.
410 
411  lagrange_free(&p); // deallocate memory allocated by setup
412  --------------------------------------------------------------------------*/
413 
414 void lagrange_0( lagrange_data* p, realType x )
415 {
416  unsigned i, n = p->n;
417  for( i = 0; i < n; ++i )
418  p->d[i] = x - p->z[i];
419  for( i = 0; i < n - 1; ++i )
420  p->u0[i + 1] = p->d[i] * p->u0[i];
421  for( i = n - 1; i; --i )
422  p->v0[i - 1] = p->d[i] * p->v0[i];
423  for( i = 0; i < n; ++i )
424  p->J[i] = p->w[i] * p->u0[i] * p->v0[i];
425 }
426 
427 void lagrange_1( lagrange_data* p, realType x )
428 {
429  unsigned i, n = p->n;
430  for( i = 0; i < n; ++i )
431  p->d[i] = x - p->z[i];
432  for( i = 0; i < n - 1; ++i )
433  p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i];
434  for( i = n - 1; i; --i )
435  p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i];
436  for( i = 0; i < n; ++i )
437  p->J[i] = p->w[i] * p->u0[i] * p->v0[i], p->D[i] = p->w[i] * ( p->u1[i] * p->v0[i] + p->u0[i] * p->v1[i] );
438 }
439 
440 void lagrange_2( lagrange_data* p, realType x )
441 {
442  unsigned i, n = p->n;
443  for( i = 0; i < n; ++i )
444  p->d[i] = x - p->z[i];
445  for( i = 0; i < n - 1; ++i )
446  p->u0[i + 1] = p->d[i] * p->u0[i], p->u1[i + 1] = p->d[i] * p->u1[i] + p->u0[i],
447  p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
448  for( i = n - 1; i; --i )
449  p->v0[i - 1] = p->d[i] * p->v0[i], p->v1[i - 1] = p->d[i] * p->v1[i] + p->v0[i],
450  p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
451  for( i = 0; i < n; ++i )
452  p->J[i] = p->w[i] * p->u0[i] * p->v0[i], p->D[i] = p->w[i] * ( p->u1[i] * p->v0[i] + p->u0[i] * p->v1[i] ),
453  p->D2[i] = p->w[i] * ( p->u2[i] * p->v0[i] + 2 * p->u1[i] * p->v1[i] + p->u0[i] * p->v2[i] );
454 }
455 
456 void lagrange_2u( lagrange_data* p )
457 {
458  unsigned i, n = p->n;
459  for( i = 0; i < n - 1; ++i )
460  p->u2[i + 1] = p->d[i] * p->u2[i] + 2 * p->u1[i];
461  for( i = n - 1; i; --i )
462  p->v2[i - 1] = p->d[i] * p->v2[i] + 2 * p->v1[i];
463  for( i = 0; i < n; ++i )
464  p->D2[i] = p->w[i] * ( p->u2[i] * p->v0[i] + 2 * p->u1[i] * p->v1[i] + p->u0[i] * p->v2[i] );
465 }
466 
467 void lagrange_setup( lagrange_data* p, const realType* z, unsigned n )
468 {
469  unsigned i, j;
470  p->n = n, p->z = z;
471  p->w = tmalloc( realType, 17 * n );
472  p->d = p->w + n;
473  p->J = p->d + n, p->D = p->J + n, p->D2 = p->D + n;
474  p->u0 = p->D2 + n, p->v0 = p->u0 + n;
475  p->u1 = p->v0 + n, p->v1 = p->u1 + n;
476  p->u2 = p->v1 + n, p->v2 = p->u2 + n;
477  p->J_z0 = p->v2 + n, p->D_z0 = p->J_z0 + n, p->D2_z0 = p->D_z0 + n;
478  p->J_zn = p->D2_z0 + n, p->D_zn = p->J_zn + n, p->D2_zn = p->D_zn + n;
479  for( i = 0; i < n; ++i )
480  {
481  realType ww = 1, zi = z[i];
482  for( j = 0; j < i; ++j )
483  ww *= zi - z[j];
484  for( ++j; j < n; ++j )
485  ww *= zi - z[j];
486  p->w[i] = 1 / ww;
487  }
488  p->u0[0] = p->v0[n - 1] = 1;
489  p->u1[0] = p->v1[n - 1] = 0;
490  p->u2[0] = p->v2[n - 1] = 0;
491  lagrange_2( p, z[0] );
492  memcpy( p->J_z0, p->J, 3 * n * sizeof( realType ) );
493  lagrange_2( p, z[n - 1] );
494  memcpy( p->J_zn, p->J, 3 * n * sizeof( realType ) );
495 }
496 
497 void lagrange_free( lagrange_data* p )
498 {
499  free( p->w );
500 }