 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
#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.
| 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.
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().
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().
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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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 }
| 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().
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().
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().
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().
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().
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().
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().
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().
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().
| 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().
|
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().
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().
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().
1.9.1.