Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Functions
poly.c File Reference
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include "moab/FindPtFuncs.h"
+ Include dependency graph for poly.c:

Go to the source code of this file.

Functions

void legendre_matrix (const realType *x, int m, realType *P, int n)
 
static void legendre_row_even (realType x, realType *P, int n)
 
static void legendre_row_odd (realType x, realType *P, int n)
 
void legendre_row (realType x, realType *P, int n)
 
void legendre_matrix_t (const realType *x, int m, realType *P, int n)
 
static realType legendre (int n, realType x)
 
static realType legendre_d1 (int n, realType x)
 
static realType legendre_d2 (int n, realType x)
 
void gauss_nodes (realType *z, int n)
 
static void lobatto_nodes_aux (realType *z, int n)
 
void lobatto_nodes (realType *z, int n)
 
void gauss_weights (const realType *z, realType *w, int n)
 
void lobatto_weights (const realType *z, realType *w, int n)
 
void gauss_to_legendre (const realType *z, const realType *w, int n, realType *J)
 
void gauss_to_legendre_t (const realType *z, const realType *w, int n, realType *J)
 
void lobatto_to_legendre (const realType *z, const realType *w, int n, realType *J)
 
void lagrange_weights (const realType *z, unsigned n, const realType *x, unsigned m, realType *J, realType *work)
 
void lagrange_weights_deriv (const realType *z, unsigned n, const realType *x, unsigned m, realType *J, realType *D, realType *work)
 
void lagrange_0 (lagrange_data *p, realType x)
 
void lagrange_1 (lagrange_data *p, realType x)
 
void lagrange_2 (lagrange_data *p, realType x)
 
void lagrange_2u (lagrange_data *p)
 
void lagrange_setup (lagrange_data *p, const realType *z, unsigned n)
 
void lagrange_free (lagrange_data *p)
 

Function Documentation

◆ gauss_nodes()

void gauss_nodes ( realType *  z,
int  n 
)

Definition at line 155 of file poly.c.

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 }

References legendre(), legendre_d1(), mbabs, mbcos, MOAB_POLY_EPS, and MOAB_POLY_PI.

◆ gauss_to_legendre()

void gauss_to_legendre ( const realType *  z,
const realType *  w,
int  n,
realType *  J 
)

Definition at line 239 of file poly.c.

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 }

References legendre_matrix_t().

◆ gauss_to_legendre_t()

void gauss_to_legendre_t ( const realType *  z,
const realType *  w,
int  n,
realType *  J 
)

Definition at line 255 of file poly.c.

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 }

References legendre_matrix().

◆ gauss_weights()

void gauss_weights ( const realType *  z,
realType *  w,
int  n 
)

Definition at line 199 of file poly.c.

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 }

References legendre().

◆ lagrange_0()

void lagrange_0 ( lagrange_data *  p,
realType  x 
)

Definition at line 414 of file poly.c.

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 }

References lagrange_data::d, lagrange_data::J, lagrange_data::n, lagrange_data::u0, lagrange_data::v0, lagrange_data::w, and lagrange_data::z.

Referenced by moab::SpectralQuad::evalFcn(), moab::Element::SpectralHex::evaluate(), moab::Element::SpectralQuad::evaluate(), moab::element_utility::Spectral_hex_map< _Matrix >::evaluate(), moab::Element::SpectralHex::evaluate_scalar_field(), moab::Element::SpectralQuad::evaluate_scalar_field(), moab::element_utility::Spectral_hex_map< _Matrix >::evaluate_scalar_field(), findpt_weights_2(), findpt_weights_3(), and moab::ElemUtil::hex_eval().

◆ lagrange_1()

void lagrange_1 ( lagrange_data *  p,
realType  x 
)

Definition at line 427 of file poly.c.

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 }

References lagrange_data::D, lagrange_data::d, lagrange_data::J, lagrange_data::n, lagrange_data::u0, lagrange_data::u1, lagrange_data::v0, lagrange_data::v1, lagrange_data::w, and lagrange_data::z.

Referenced by opt_area_set_2(), opt_edge_set_2(), opt_edge_set_3(), opt_face_set_3(), and opt_vol_set_3().

◆ lagrange_2()

void lagrange_2 ( lagrange_data *  p,
realType  x 
)

Definition at line 440 of file poly.c.

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 }

References lagrange_data::D, lagrange_data::d, lagrange_data::D2, lagrange_data::J, lagrange_data::n, lagrange_data::u0, lagrange_data::u1, lagrange_data::u2, lagrange_data::v0, lagrange_data::v1, lagrange_data::v2, lagrange_data::w, and lagrange_data::z.

Referenced by lagrange_setup().

◆ lagrange_2u()

void lagrange_2u ( lagrange_data *  p)

Definition at line 456 of file poly.c.

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 }

References lagrange_data::d, lagrange_data::D2, lagrange_data::n, lagrange_data::u0, lagrange_data::u1, lagrange_data::u2, lagrange_data::v0, lagrange_data::v1, lagrange_data::v2, and lagrange_data::w.

Referenced by opt_edge_hess_2(), opt_edge_hess_3(), and opt_face_hess_3().

◆ lagrange_free()

void lagrange_free ( lagrange_data *  p)

Definition at line 497 of file poly.c.

498 {
499  free( p->w );
500 }

References lagrange_data::w.

Referenced by moab::element_utility::Spectral_hex_map< _Matrix >::free_data(), moab::Element::SpectralHex::freedata(), moab::Element::SpectralQuad::freedata(), moab::ElemUtil::hex_eval(), and moab::ElemUtil::hex_findpt().

◆ lagrange_setup()

void lagrange_setup ( lagrange_data *  p,
const realType *  z,
unsigned  n 
)

Definition at line 467 of file poly.c.

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 }

References lagrange_data::D, lagrange_data::d, lagrange_data::D2, lagrange_data::D2_z0, lagrange_data::D2_zn, lagrange_data::D_z0, lagrange_data::D_zn, lagrange_data::J, lagrange_data::J_z0, lagrange_data::J_zn, lagrange_2(), lagrange_data::n, tmalloc, lagrange_data::u0, lagrange_data::u1, lagrange_data::u2, lagrange_data::v0, lagrange_data::v1, lagrange_data::v2, lagrange_data::w, and lagrange_data::z.

Referenced by findpt_setup_2(), findpt_setup_3(), moab::ElemUtil::hex_eval(), moab::ElemUtil::hex_findpt(), moab::Element::SpectralHex::Init(), moab::Element::SpectralQuad::Init(), moab::SpectralHex::initFcn(), and moab::element_utility::Spectral_hex_map< _Matrix >::initialize_spectral_hex().

◆ lagrange_weights()

void lagrange_weights ( const realType *  z,
unsigned  n,
const realType *  x,
unsigned  m,
realType *  J,
realType *  work 
)

Definition at line 320 of file poly.c.

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 }

◆ lagrange_weights_deriv()

void lagrange_weights_deriv ( const realType *  z,
unsigned  n,
const realType *  x,
unsigned  m,
realType *  J,
realType *  D,
realType *  work 
)

Definition at line 354 of file poly.c.

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 }

Referenced by lob_bnd_base_setup(), obbox_setup_2(), and obbox_setup_3().

◆ legendre()

static realType legendre ( int  n,
realType  x 
)
static

Definition at line 104 of file poly.c.

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 }

Referenced by gauss_nodes(), gauss_weights(), and lobatto_weights().

◆ legendre_d1()

static realType legendre_d1 ( int  n,
realType  x 
)
static

Definition at line 119 of file poly.c.

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 }

Referenced by gauss_nodes(), and lobatto_nodes_aux().

◆ legendre_d2()

static realType legendre_d2 ( int  n,
realType  x 
)
static

Definition at line 134 of file poly.c.

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 }

Referenced by lobatto_nodes_aux().

◆ legendre_matrix()

void legendre_matrix ( const realType *  x,
int  m,
realType *  P,
int  n 
)

Definition at line 30 of file poly.c.

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 }

Referenced by gauss_to_legendre_t().

◆ legendre_matrix_t()

void legendre_matrix_t ( const realType *  x,
int  m,
realType *  P,
int  n 
)

Definition at line 87 of file poly.c.

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 }

References legendre_row_even(), and legendre_row_odd().

Referenced by gauss_to_legendre().

◆ legendre_row()

void legendre_row ( realType  x,
realType *  P,
int  n 
)

Definition at line 75 of file poly.c.

76 {
77  if( n & 1 )
78  legendre_row_odd( x, P, n );
79  else
80  legendre_row_even( x, P, n );
81 }

References legendre_row_even(), and legendre_row_odd().

◆ legendre_row_even()

static void legendre_row_even ( realType  x,
realType *  P,
int  n 
)
static

Definition at line 48 of file poly.c.

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 }

Referenced by legendre_matrix_t(), legendre_row(), and lobatto_to_legendre().

◆ legendre_row_odd()

static void legendre_row_odd ( realType  x,
realType *  P,
int  n 
)
static

Definition at line 61 of file poly.c.

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 }

Referenced by legendre_matrix_t(), legendre_row(), and lobatto_to_legendre().

◆ lobatto_nodes()

void lobatto_nodes ( realType *  z,
int  n 
)

Definition at line 193 of file poly.c.

194 {
195  z[0] = -1, z[n - 1] = 1;
196  lobatto_nodes_aux( &z[1], n - 2 );
197 }

References lobatto_nodes_aux().

Referenced by findpt_setup_2(), findpt_setup_3(), hash_getbb_2(), hash_getbb_3(), moab::ElemUtil::hex_eval(), moab::ElemUtil::hex_findpt(), moab::Element::SpectralHex::Init(), moab::Element::SpectralQuad::Init(), moab::SpectralHex::initFcn(), and moab::element_utility::Spectral_hex_map< _Matrix >::initialize_spectral_hex().

◆ lobatto_nodes_aux()

static void lobatto_nodes_aux ( realType *  z,
int  n 
)
static

Definition at line 174 of file poly.c.

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 }

References legendre_d1(), legendre_d2(), mbabs, mbcos, MOAB_POLY_EPS, and MOAB_POLY_PI.

Referenced by lobatto_nodes().

◆ lobatto_to_legendre()

void lobatto_to_legendre ( const realType *  z,
const realType *  w,
int  n,
realType *  J 
)

Definition at line 275 of file poly.c.

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 }

References legendre_row_even(), legendre_row_odd(), and moab::sum().

◆ lobatto_weights()

void lobatto_weights ( const realType *  z,
realType *  w,
int  n 
)

Definition at line 211 of file poly.c.

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 }

References legendre().

Referenced by hash_getbb_2(), and hash_getbb_3().