MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #ifndef SPECTRALFUNCS_HPP 00002 #define SPECTRALFUNCS_HPP 00003 00004 #include "float.h" 00005 00006 //====================================================== 00007 // from types.h 00008 //====================================================== 00009 00010 /* integer type to use for everything */ 00011 #if defined( USE_LONG ) 00012 #define INTEGER long 00013 #elif defined( USE_LONG_LONG ) 00014 #define INTEGER long long 00015 #elif defined( USE_SHORT ) 00016 #define INTEGER short 00017 #else 00018 #define INTEGER int 00019 #endif 00020 00021 /* when defined, use the given type for global indices instead of INTEGER */ 00022 #if defined( USE_GLOBAL_LONG_LONG ) 00023 #define GLOBAL_INT long long 00024 #elif defined( USE_GLOBAL_LONG ) 00025 #define GLOBAL_INT long 00026 #else 00027 #define GLOBAL_INT long 00028 #endif 00029 00030 /* floating point type to use for everything */ 00031 #if defined( USE_FLOAT ) 00032 typedef float real; 00033 #define floorr floorf 00034 #define ceilr ceilf 00035 #define sqrtr sqrtf 00036 #define fabsr fabsf 00037 #define cosr cosf 00038 #define sinr sinf 00039 #define EPS ( 128 * FLT_EPSILON ) 00040 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923F 00041 #elif defined( USE_LONG_DOUBLE ) 00042 typedef long double real; 00043 #define floorr floorl 00044 #define ceilr ceill 00045 #define sqrtr sqrtl 00046 #define fabsr fabsl 00047 #define cosr cosl 00048 #define sinr sinl 00049 #define EPS ( 128 * LDBL_EPSILON ) 00050 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923L 00051 #else 00052 typedef double real; 00053 #define floorr floor 00054 #define ceilr ceil 00055 #define sqrtr sqrt 00056 #define fabsr fabs 00057 #define cosr cos 00058 #define sinr sin 00059 #define EPS ( 128 * DBL_EPSILON ) 00060 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923 00061 #endif 00062 00063 /* apparently uint and ulong can be defined already in standard headers */ 00064 #define uint uint_ 00065 #define ulong ulong_ 00066 #define sint sint_ 00067 #define slong slong_ 00068 00069 typedef signed INTEGER sint; 00070 typedef unsigned INTEGER uint; 00071 #undef INTEGER 00072 00073 #ifdef GLOBAL_INT 00074 typedef signed GLOBAL_INT slong; 00075 typedef unsigned GLOBAL_INT ulong; 00076 #else 00077 typedef sint slong; 00078 typedef uint ulong; 00079 #endif 00080 00081 //====================================================== 00082 // from poly.h 00083 //====================================================== 00084 00085 /* 00086 For brevity's sake, some names have been shortened 00087 Quadrature rules 00088 Gauss -> Gauss-Legendre quadrature (open) 00089 Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends) 00090 Polynomial bases 00091 Legendre -> Legendre basis 00092 Gauss -> Lagrangian basis using Gauss quadrature nodes 00093 Lobatto -> Lagrangian basis using Lobatto quadrature nodes 00094 */ 00095 00096 /*-------------------------------------------------------------------------- 00097 Legendre Polynomial Matrix Computation 00098 (compute P_i(x_j) for i = 0, ..., n and a given set of x) 00099 --------------------------------------------------------------------------*/ 00100 00101 /* precondition: n >= 1 00102 inner index is x index (0 ... m-1); 00103 outer index is Legendre polynomial number (0 ... n) 00104 */ 00105 void legendre_matrix( const real* x, int m, real* P, int n ); 00106 00107 /* precondition: n >= 1 00108 inner index is Legendre polynomial number (0 ... n) 00109 outer index is x index (0 ... m-1); 00110 */ 00111 void legendre_matrix_t( const real* x, int m, real* P, int n ); 00112 00113 /* precondition: n >= 1 00114 compute P_i(x) with i = 0 ... n 00115 */ 00116 void legendre_row( real x, real* P, int n ); 00117 00118 /*-------------------------------------------------------------------------- 00119 Quadrature Nodes and Weights Calculation 00120 00121 call the _nodes function before calling the _weights function 00122 --------------------------------------------------------------------------*/ 00123 00124 void gauss_nodes( real* z, int n ); /* n nodes (order = 2n-1) */ 00125 void lobatto_nodes( real* z, int n ); /* n nodes (order = 2n-3) */ 00126 00127 void gauss_weights( const real* z, real* w, int n ); 00128 void lobatto_weights( const real* z, real* w, int n ); 00129 00130 /*-------------------------------------------------------------------------- 00131 Lagrangian to Legendre Change-of-basis Matrix 00132 --------------------------------------------------------------------------*/ 00133 00134 /* precondition: n >= 2 00135 given the Gauss quadrature rule (z,w,n), compute the square matrix J 00136 for transforming from the Gauss basis to the Legendre basis: 00137 00138 u_legendre(i) = sum_j J(i,j) u_gauss(j) 00139 00140 computes J = .5 (2i+1) w P (z ) 00141 ij j i j 00142 00143 in column major format (inner index is i, the Legendre index) 00144 */ 00145 void gauss_to_legendre( const real* z, const real* w, int n, real* J ); 00146 00147 /* precondition: n >= 2 00148 same as above, but 00149 in row major format (inner index is j, the Gauss index) 00150 */ 00151 void gauss_to_legendre_t( const real* z, const real* w, int n, real* J ); 00152 00153 /* precondition: n >= 3 00154 given the Lobatto quadrature rule (z,w,n), compute the square matrix J 00155 for transforming from the Lobatto basis to the Legendre basis: 00156 00157 u_legendre(i) = sum_j J(i,j) u_lobatto(j) 00158 00159 in column major format (inner index is i, the Legendre index) 00160 */ 00161 void lobatto_to_legendre( const real* z, const real* w, int n, real* J ); 00162 00163 /*-------------------------------------------------------------------------- 00164 Lagrangian basis function evaluation 00165 --------------------------------------------------------------------------*/ 00166 00167 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00168 evaluate all Lagrangian basis functions at all points x 00169 00170 inner index of output J is the basis function index (row-major format) 00171 provide work array with space for 4*n doubles 00172 */ 00173 void lagrange_weights( const real* z, unsigned n, const real* x, unsigned m, real* J, real* work ); 00174 00175 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00176 evaluate all Lagrangian basis functions and their derivatives 00177 00178 inner index of outputs J,D is the basis function index (row-major format) 00179 provide work array with space for 6*n doubles 00180 */ 00181 void lagrange_weights_deriv( const real* z, unsigned n, const real* x, unsigned m, real* J, real* D, real* work ); 00182 00183 /*-------------------------------------------------------------------------- 00184 Speedy Lagrangian Interpolation 00185 00186 Usage: 00187 00188 lagrange_data p; 00189 lagrange_setup(&p,z,n); * setup for nodes z[0 ... n-1] * 00190 00191 the weights 00192 p->J [0 ... n-1] interpolation weights 00193 p->D [0 ... n-1] 1st derivative weights 00194 p->D2[0 ... n-1] 2nd derivative weights 00195 are computed for a given x with: 00196 lagrange_0(p,x); * compute p->J * 00197 lagrange_1(p,x); * compute p->J, p->D * 00198 lagrange_2(p,x); * compute p->J, p->D, p->D2 * 00199 lagrange_2u(p); * compute p->D2 after call of lagrange_1(p,x); * 00200 These functions use the z array supplied to setup 00201 (that pointer should not be freed between calls) 00202 Weights for x=z[0] and x=z[n-1] are computed during setup; access as: 00203 p->J_z0, etc. and p->J_zn, etc. 00204 00205 lagrange_free(&p); * deallocate memory allocated by setup * 00206 --------------------------------------------------------------------------*/ 00207 00208 typedef struct 00209 { 00210 unsigned n; /* number of Lagrange nodes */ 00211 const real* z; /* Lagrange nodes (user-supplied) */ 00212 real *J, *D, *D2; /* weights for 0th,1st,2nd derivatives */ 00213 real *J_z0, *D_z0, *D2_z0; /* ditto at z[0] (computed at setup) */ 00214 real *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */ 00215 real *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data */ 00216 } lagrange_data; 00217 00218 void lagrange_setup( lagrange_data* p, const real* z, unsigned n ); 00219 void lagrange_free( lagrange_data* p ); 00220 00221 void lagrange_0( lagrange_data* p, real x ); 00222 void lagrange_1( lagrange_data* p, real x ); 00223 void lagrange_2( lagrange_data* p, real x ); 00224 void lagrange_2u( lagrange_data* p ); 00225 00226 //====================================================== 00227 // from tensor.h 00228 //====================================================== 00229 00230 /*-------------------------------------------------------------------------- 00231 1-,2-,3-d Tensor Application 00232 00233 the 3d case: 00234 tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2) 00235 gives v = [ R (x) S (x) T ] u 00236 where R is mr x nr, S is ms x ns, T is mt x nt, 00237 each in row- or column-major format according to f := r | c 00238 u is nr x ns x nt in column-major format (inner index is r) 00239 v is mr x ms x mt in column-major format (inner index is r) 00240 --------------------------------------------------------------------------*/ 00241 00242 void tensor_c1( const real* R, unsigned mr, unsigned nr, const real* u, real* v ); 00243 void tensor_r1( const real* R, unsigned mr, unsigned nr, const real* u, real* v ); 00244 00245 /* work holds mr*ns reals */ 00246 void tensor_c2( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* u, 00247 real* v, real* work ); 00248 void tensor_r2( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* u, 00249 real* v, real* work ); 00250 00251 /* work1 holds mr*ns*nt reals, 00252 work2 holds mr*ms*nt reals */ 00253 void tensor_c3( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* T, 00254 unsigned mt, unsigned nt, const real* u, real* v, real* work1, real* work2 ); 00255 void tensor_r3( const real* R, unsigned mr, unsigned nr, const real* S, unsigned ms, unsigned ns, const real* T, 00256 unsigned mt, unsigned nt, const real* u, real* v, real* work1, real* work2 ); 00257 00258 /*-------------------------------------------------------------------------- 00259 1-,2-,3-d Tensor Application of Row Vectors (for Interpolation) 00260 00261 the 3d case: 00262 v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work) 00263 same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2): 00264 gives v = [ Jr (x) Js (x) Jt ] u 00265 where Jr, Js, Jt are row vectors (interpolation weights) 00266 u is nr x ns x nt in column-major format (inner index is r) 00267 v is a scalar 00268 --------------------------------------------------------------------------*/ 00269 00270 real tensor_i1( const real* Jr, unsigned nr, const real* u ); 00271 00272 /* work holds ns reals */ 00273 real tensor_i2( const real* Jr, unsigned nr, const real* Js, unsigned ns, const real* u, real* work ); 00274 00275 /* work holds ns*nt + nt reals */ 00276 real tensor_i3( const real* Jr, unsigned nr, const real* Js, unsigned ns, const real* Jt, unsigned nt, const real* u, 00277 real* work ); 00278 00279 /*-------------------------------------------------------------------------- 00280 1-,2-,3-d Tensor Application of Row Vectors 00281 for simultaneous Interpolation and Gradient computation 00282 00283 the 3d case: 00284 v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work) 00285 gives v = [ Jr (x) Js (x) Jt ] u 00286 g_0 = [ Dr (x) Js (x) Jt ] u 00287 g_1 = [ Jr (x) Ds (x) Jt ] u 00288 g_2 = [ Jr (x) Js (x) Dt ] u 00289 where Jr,Dr,Js,Ds,Jt,Dt are row vectors 00290 (interpolation & derivative weights) 00291 u is nr x ns x nt in column-major format (inner index is r) 00292 v is a scalar, g is an array of 3 reals 00293 --------------------------------------------------------------------------*/ 00294 00295 real tensor_ig1( const real* Jr, const real* Dr, unsigned nr, const real* u, real* g ); 00296 00297 /* work holds 2*ns reals */ 00298 real tensor_ig2( const real* Jr, const real* Dr, unsigned nr, const real* Js, const real* Ds, unsigned ns, 00299 const real* u, real* g, real* work ); 00300 00301 /* work holds 2*ns*nt + 3*ns reals */ 00302 real tensor_ig3( const real* Jr, const real* Dr, unsigned nr, const real* Js, const real* Ds, unsigned ns, 00303 const real* Jt, const real* Dt, unsigned nt, const real* u, real* g, real* work ); 00304 00305 //====================================================== 00306 // from findpt.h 00307 //====================================================== 00308 00309 typedef struct 00310 { 00311 const real* xw[2]; /* geometry data */ 00312 real* z[2]; /* lobatto nodes */ 00313 lagrange_data ld[2]; /* interpolation, derivative weights & data */ 00314 unsigned nptel; /* nodes per element */ 00315 struct findpt_hash_data_2* hash; /* geometric hashing data */ 00316 struct findpt_listel *list, **sorted, **end; 00317 /* pre-allocated list of elements to 00318 check (found by hashing), and 00319 pre-allocated list of pointers into 00320 the first list for sorting */ 00321 struct findpt_opt_data_2* od; /* data for the optimization algorithm */ 00322 real* od_work; 00323 } findpt_data_2; 00324 00325 typedef struct 00326 { 00327 const real* xw[3]; /* geometry data */ 00328 real* z[3]; /* lobatto nodes */ 00329 lagrange_data ld[3]; /* interpolation, derivative weights & data */ 00330 unsigned nptel; /* nodes per element */ 00331 struct findpt_hash_data_3* hash; /* geometric hashing data */ 00332 struct findpt_listel *list, **sorted, **end; 00333 /* pre-allocated list of elements to 00334 check (found by hashing), and 00335 pre-allocated list of pointers into 00336 the first list for sorting */ 00337 struct findpt_opt_data_3* od; /* data for the optimization algorithm */ 00338 real* od_work; 00339 } findpt_data_3; 00340 00341 findpt_data_2* findpt_setup_2( const real* const xw[2], const unsigned n[2], uint nel, uint max_hash_size, 00342 real bbox_tol ); 00343 findpt_data_3* findpt_setup_3( const real* const xw[3], const unsigned n[3], uint nel, uint max_hash_size, 00344 real bbox_tol ); 00345 00346 void findpt_free_2( findpt_data_2* p ); 00347 void findpt_free_3( findpt_data_3* p ); 00348 00349 const real* findpt_allbnd_2( const findpt_data_2* p ); 00350 const real* findpt_allbnd_3( const findpt_data_3* p ); 00351 00352 typedef int ( *findpt_func )( void*, const real*, int, uint*, real*, real* ); 00353 int findpt_2( findpt_data_2* p, const real x[2], int guess, uint* el, real r[2], real* dist ); 00354 int findpt_3( findpt_data_3* p, const real x[3], int guess, uint* el, real r[3], real* dist ); 00355 00356 inline void findpt_weights_2( findpt_data_2* p, const real r[2] ) 00357 { 00358 lagrange_0( &p->ld[0], r[0] ); 00359 lagrange_0( &p->ld[1], r[1] ); 00360 } 00361 00362 inline void findpt_weights_3( findpt_data_3* p, const real r[3] ) 00363 { 00364 lagrange_0( &p->ld[0], r[0] ); 00365 lagrange_0( &p->ld[1], r[1] ); 00366 lagrange_0( &p->ld[2], r[2] ); 00367 } 00368 00369 inline double findpt_eval_2( findpt_data_2* p, const real* u ) 00370 { 00371 return tensor_i2( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, u, p->od_work ); 00372 } 00373 00374 inline double findpt_eval_3( findpt_data_3* p, const real* u ) 00375 { 00376 return tensor_i3( p->ld[0].J, p->ld[0].n, p->ld[1].J, p->ld[1].n, p->ld[2].J, p->ld[2].n, u, p->od_work ); 00377 } 00378 00379 //====================================================== 00380 // from extrafindpt.h 00381 //====================================================== 00382 00383 typedef struct 00384 { 00385 unsigned constraints; 00386 unsigned dn, d1, d2; 00387 real *x[3], *fdn[3]; 00388 } opt_face_data_3; 00389 00390 typedef struct 00391 { 00392 unsigned constraints; 00393 unsigned de, d1, d2; 00394 real *x[3], *fd1[3], *fd2[3]; 00395 } opt_edge_data_3; 00396 00397 typedef struct 00398 { 00399 unsigned constraints; 00400 real x[3], jac[9]; 00401 } opt_point_data_3; 00402 00403 typedef struct 00404 { 00405 lagrange_data* ld; 00406 unsigned size[4]; 00407 const real* elx[3]; 00408 opt_face_data_3 fd; 00409 opt_edge_data_3 ed; 00410 opt_point_data_3 pd; 00411 real* work; 00412 real x[3], jac[9]; 00413 } opt_data_3; 00414 00415 void opt_alloc_3( opt_data_3* p, lagrange_data* ld ); 00416 void opt_free_3( opt_data_3* p ); 00417 double opt_findpt_3( opt_data_3* p, const real* const elx[3], const real xstar[3], real r[3], unsigned* constr ); 00418 void opt_vol_set_intp_3( opt_data_3* p, const real r[3] ); 00419 00420 const unsigned opt_no_constraints_2 = 3 + 1; 00421 const unsigned opt_no_constraints_3 = 9 + 3 + 1; 00422 00423 /* for 2d spectralQuad */ 00424 /*-------------------------------------------------------------------------- 00425 00426 2 - D 00427 00428 --------------------------------------------------------------------------*/ 00429 00430 typedef struct 00431 { 00432 unsigned constraints; 00433 unsigned de, d1; 00434 real *x[2], *fd1[2]; 00435 } opt_edge_data_2; 00436 00437 typedef struct 00438 { 00439 unsigned constraints; 00440 real x[2], jac[4]; 00441 } opt_point_data_2; 00442 00443 typedef struct 00444 { 00445 lagrange_data* ld; 00446 unsigned size[3]; 00447 const real* elx[2]; 00448 opt_edge_data_2 ed; 00449 opt_point_data_2 pd; 00450 real* work; 00451 real x[2], jac[4]; 00452 } opt_data_2; 00453 void opt_alloc_2( opt_data_2* p, lagrange_data* ld ); 00454 void opt_free_2( opt_data_2* p ); 00455 double opt_findpt_2( opt_data_2* p, const real* const elx[2], const real xstar[2], real r[2], unsigned* constr ); 00456 00457 #endif