moab
|
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 /*-------------------------------------------------------------------------- 00120 Quadrature Nodes and Weights Calculation 00121 00122 call the _nodes function before calling the _weights function 00123 --------------------------------------------------------------------------*/ 00124 00125 void gauss_nodes(real *z, int n); /* n nodes (order = 2n-1) */ 00126 void lobatto_nodes(real *z, int n); /* n nodes (order = 2n-3) */ 00127 00128 void gauss_weights(const real *z, real *w, int n); 00129 void lobatto_weights(const real *z, real *w, int n); 00130 00131 /*-------------------------------------------------------------------------- 00132 Lagrangian to Legendre Change-of-basis Matrix 00133 --------------------------------------------------------------------------*/ 00134 00135 /* precondition: n >= 2 00136 given the Gauss quadrature rule (z,w,n), compute the square matrix J 00137 for transforming from the Gauss basis to the Legendre basis: 00138 00139 u_legendre(i) = sum_j J(i,j) u_gauss(j) 00140 00141 computes J = .5 (2i+1) w P (z ) 00142 ij j i j 00143 00144 in column major format (inner index is i, the Legendre index) 00145 */ 00146 void gauss_to_legendre(const real *z, const real *w, int n, real *J); 00147 00148 /* precondition: n >= 2 00149 same as above, but 00150 in row major format (inner index is j, the Gauss index) 00151 */ 00152 void gauss_to_legendre_t(const real *z, const real *w, int n, real *J); 00153 00154 /* precondition: n >= 3 00155 given the Lobatto quadrature rule (z,w,n), compute the square matrix J 00156 for transforming from the Lobatto basis to the Legendre basis: 00157 00158 u_legendre(i) = sum_j J(i,j) u_lobatto(j) 00159 00160 in column major format (inner index is i, the Legendre index) 00161 */ 00162 void lobatto_to_legendre(const real *z, const real *w, int n, real *J); 00163 00164 /*-------------------------------------------------------------------------- 00165 Lagrangian basis function evaluation 00166 --------------------------------------------------------------------------*/ 00167 00168 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00169 evaluate all Lagrangian basis functions at all points x 00170 00171 inner index of output J is the basis function index (row-major format) 00172 provide work array with space for 4*n doubles 00173 */ 00174 void lagrange_weights(const real *z, unsigned n, 00175 const real *x, unsigned m, 00176 real *J, real *work); 00177 00178 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00179 evaluate all Lagrangian basis functions and their derivatives 00180 00181 inner index of outputs J,D is the basis function index (row-major format) 00182 provide work array with space for 6*n doubles 00183 */ 00184 void lagrange_weights_deriv(const real *z, unsigned n, 00185 const real *x, unsigned m, 00186 real *J, real *D, real *work); 00187 00188 /*-------------------------------------------------------------------------- 00189 Speedy Lagrangian Interpolation 00190 00191 Usage: 00192 00193 lagrange_data p; 00194 lagrange_setup(&p,z,n); * setup for nodes z[0 ... n-1] * 00195 00196 the weights 00197 p->J [0 ... n-1] interpolation weights 00198 p->D [0 ... n-1] 1st derivative weights 00199 p->D2[0 ... n-1] 2nd derivative weights 00200 are computed for a given x with: 00201 lagrange_0(p,x); * compute p->J * 00202 lagrange_1(p,x); * compute p->J, p->D * 00203 lagrange_2(p,x); * compute p->J, p->D, p->D2 * 00204 lagrange_2u(p); * compute p->D2 after call of lagrange_1(p,x); * 00205 These functions use the z array supplied to setup 00206 (that pointer should not be freed between calls) 00207 Weights for x=z[0] and x=z[n-1] are computed during setup; access as: 00208 p->J_z0, etc. and p->J_zn, etc. 00209 00210 lagrange_free(&p); * deallocate memory allocated by setup * 00211 --------------------------------------------------------------------------*/ 00212 00213 typedef struct { 00214 unsigned n; /* number of Lagrange nodes */ 00215 const real *z; /* Lagrange nodes (user-supplied) */ 00216 real *J, *D, *D2; /* weights for 0th,1st,2nd derivatives */ 00217 real *J_z0, *D_z0, *D2_z0; /* ditto at z[0] (computed at setup) */ 00218 real *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */ 00219 real *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data */ 00220 } lagrange_data; 00221 00222 void lagrange_setup(lagrange_data *p, const real *z, unsigned n); 00223 void lagrange_free(lagrange_data *p); 00224 00225 void lagrange_0(lagrange_data *p, real x) ; 00226 void lagrange_1(lagrange_data *p, real x) ; 00227 void lagrange_2(lagrange_data *p, real x) ; 00228 void lagrange_2u(lagrange_data *p) ; 00229 00230 //====================================================== 00231 // from tensor.h 00232 //====================================================== 00233 00234 /*-------------------------------------------------------------------------- 00235 1-,2-,3-d Tensor Application 00236 00237 the 3d case: 00238 tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2) 00239 gives v = [ R (x) S (x) T ] u 00240 where R is mr x nr, S is ms x ns, T is mt x nt, 00241 each in row- or column-major format according to f := r | c 00242 u is nr x ns x nt in column-major format (inner index is r) 00243 v is mr x ms x mt in column-major format (inner index is r) 00244 --------------------------------------------------------------------------*/ 00245 00246 void tensor_c1(const real *R, unsigned mr, unsigned nr, 00247 const real *u, real *v); 00248 void tensor_r1(const real *R, unsigned mr, unsigned nr, 00249 const real *u, real *v); 00250 00251 /* work holds mr*ns reals */ 00252 void tensor_c2(const real *R, unsigned mr, unsigned nr, 00253 const real *S, unsigned ms, unsigned ns, 00254 const real *u, real *v, real *work); 00255 void tensor_r2(const real *R, unsigned mr, unsigned nr, 00256 const real *S, unsigned ms, unsigned ns, 00257 const real *u, real *v, real *work); 00258 00259 /* work1 holds mr*ns*nt reals, 00260 work2 holds mr*ms*nt reals */ 00261 void tensor_c3(const real *R, unsigned mr, unsigned nr, 00262 const real *S, unsigned ms, unsigned ns, 00263 const real *T, unsigned mt, unsigned nt, 00264 const real *u, real *v, real *work1, real *work2); 00265 void tensor_r3(const real *R, unsigned mr, unsigned nr, 00266 const real *S, unsigned ms, unsigned ns, 00267 const real *T, unsigned mt, unsigned nt, 00268 const real *u, real *v, real *work1, real *work2); 00269 00270 /*-------------------------------------------------------------------------- 00271 1-,2-,3-d Tensor Application of Row Vectors (for Interpolation) 00272 00273 the 3d case: 00274 v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work) 00275 same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2): 00276 gives v = [ Jr (x) Js (x) Jt ] u 00277 where Jr, Js, Jt are row vectors (interpolation weights) 00278 u is nr x ns x nt in column-major format (inner index is r) 00279 v is a scalar 00280 --------------------------------------------------------------------------*/ 00281 00282 real tensor_i1(const real *Jr, unsigned nr, const real *u); 00283 00284 /* work holds ns reals */ 00285 real tensor_i2(const real *Jr, unsigned nr, 00286 const real *Js, unsigned ns, 00287 const real *u, real *work); 00288 00289 /* work holds ns*nt + nt reals */ 00290 real tensor_i3(const real *Jr, unsigned nr, 00291 const real *Js, unsigned ns, 00292 const real *Jt, unsigned nt, 00293 const real *u, real *work); 00294 00295 /*-------------------------------------------------------------------------- 00296 1-,2-,3-d Tensor Application of Row Vectors 00297 for simultaneous Interpolation and Gradient computation 00298 00299 the 3d case: 00300 v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work) 00301 gives v = [ Jr (x) Js (x) Jt ] u 00302 g_0 = [ Dr (x) Js (x) Jt ] u 00303 g_1 = [ Jr (x) Ds (x) Jt ] u 00304 g_2 = [ Jr (x) Js (x) Dt ] u 00305 where Jr,Dr,Js,Ds,Jt,Dt are row vectors 00306 (interpolation & derivative weights) 00307 u is nr x ns x nt in column-major format (inner index is r) 00308 v is a scalar, g is an array of 3 reals 00309 --------------------------------------------------------------------------*/ 00310 00311 real tensor_ig1(const real *Jr, const real *Dr, unsigned nr, 00312 const real *u, real *g); 00313 00314 /* work holds 2*ns reals */ 00315 real tensor_ig2(const real *Jr, const real *Dr, unsigned nr, 00316 const real *Js, const real *Ds, unsigned ns, 00317 const real *u, real *g, real *work); 00318 00319 /* work holds 2*ns*nt + 3*ns reals */ 00320 real tensor_ig3(const real *Jr, const real *Dr, unsigned nr, 00321 const real *Js, const real *Ds, unsigned ns, 00322 const real *Jt, const real *Dt, unsigned nt, 00323 const real *u, real *g, real *work); 00324 00325 //====================================================== 00326 // from findpt.h 00327 //====================================================== 00328 00329 typedef struct { 00330 const real *xw[2]; /* geometry data */ 00331 real *z[2]; /* lobatto nodes */ 00332 lagrange_data ld[2]; /* interpolation, derivative weights & data */ 00333 unsigned nptel; /* nodes per element */ 00334 struct findpt_hash_data_2 *hash; /* geometric hashing data */ 00335 struct findpt_listel *list, **sorted, **end; 00336 /* pre-allocated list of elements to 00337 check (found by hashing), and 00338 pre-allocated list of pointers into 00339 the first list for sorting */ 00340 struct findpt_opt_data_2 *od; /* data for the optimization algorithm */ 00341 real *od_work; 00342 } findpt_data_2; 00343 00344 typedef struct { 00345 const real *xw[3]; /* geometry data */ 00346 real *z[3]; /* lobatto nodes */ 00347 lagrange_data ld[3]; /* interpolation, derivative weights & data */ 00348 unsigned nptel; /* nodes per element */ 00349 struct findpt_hash_data_3 *hash; /* geometric hashing data */ 00350 struct findpt_listel *list, **sorted, **end; 00351 /* pre-allocated list of elements to 00352 check (found by hashing), and 00353 pre-allocated list of pointers into 00354 the first list for sorting */ 00355 struct findpt_opt_data_3 *od; /* data for the optimization algorithm */ 00356 real *od_work; 00357 } findpt_data_3; 00358 00359 findpt_data_2 *findpt_setup_2( 00360 const real *const xw[2], const unsigned n[2], uint nel, 00361 uint max_hash_size, real bbox_tol); 00362 findpt_data_3 *findpt_setup_3( 00363 const real *const xw[3], const unsigned n[3], uint nel, 00364 uint max_hash_size, real bbox_tol); 00365 00366 void findpt_free_2(findpt_data_2 *p); 00367 void findpt_free_3(findpt_data_3 *p); 00368 00369 const real *findpt_allbnd_2(const findpt_data_2 *p); 00370 const real *findpt_allbnd_3(const findpt_data_3 *p); 00371 00372 typedef int (*findpt_func)(void *, const real *, int, uint *, real *, real *); 00373 int findpt_2(findpt_data_2 *p, const real x[2], int guess, 00374 uint *el, real r[2], real *dist); 00375 int findpt_3(findpt_data_3 *p, const real x[3], int guess, 00376 uint *el, real r[3], real *dist); 00377 00378 inline void findpt_weights_2(findpt_data_2 *p, const real r[2]) 00379 { 00380 lagrange_0(&p->ld[0],r[0]); 00381 lagrange_0(&p->ld[1],r[1]); 00382 } 00383 00384 inline void findpt_weights_3(findpt_data_3 *p, const real r[3]) 00385 { 00386 lagrange_0(&p->ld[0],r[0]); 00387 lagrange_0(&p->ld[1],r[1]); 00388 lagrange_0(&p->ld[2],r[2]); 00389 } 00390 00391 inline double findpt_eval_2(findpt_data_2 *p, const real *u) 00392 { 00393 return tensor_i2(p->ld[0].J,p->ld[0].n, 00394 p->ld[1].J,p->ld[1].n, 00395 u, p->od_work); 00396 } 00397 00398 inline double findpt_eval_3(findpt_data_3 *p, const real *u) 00399 { 00400 return tensor_i3(p->ld[0].J,p->ld[0].n, 00401 p->ld[1].J,p->ld[1].n, 00402 p->ld[2].J,p->ld[2].n, 00403 u, p->od_work); 00404 } 00405 00406 //====================================================== 00407 // from extrafindpt.h 00408 //====================================================== 00409 00410 typedef struct { 00411 unsigned constraints; 00412 unsigned dn, d1, d2; 00413 real *x[3], *fdn[3]; 00414 } opt_face_data_3; 00415 00416 typedef struct { 00417 unsigned constraints; 00418 unsigned de, d1, d2; 00419 real *x[3], *fd1[3], *fd2[3]; 00420 } opt_edge_data_3; 00421 00422 typedef struct { 00423 unsigned constraints; 00424 real x[3], jac[9]; 00425 } opt_point_data_3; 00426 00427 typedef struct { 00428 lagrange_data *ld; 00429 unsigned size[4]; 00430 const real *elx[3]; 00431 opt_face_data_3 fd; 00432 opt_edge_data_3 ed; 00433 opt_point_data_3 pd; 00434 real *work; 00435 real x[3], jac[9]; 00436 } opt_data_3; 00437 00438 00439 void opt_alloc_3(opt_data_3 *p, lagrange_data *ld); 00440 void opt_free_3(opt_data_3 *p); 00441 double opt_findpt_3(opt_data_3 *p, const real *const elx[3], 00442 const real xstar[3], real r[3], unsigned *constr); 00443 void opt_vol_set_intp_3(opt_data_3 *p, const real r[3]); 00444 00445 const unsigned opt_no_constraints_2 = 3+1; 00446 const unsigned opt_no_constraints_3 = 9+3+1; 00447 00448 /* for 2d spectralQuad */ 00449 /*-------------------------------------------------------------------------- 00450 00451 2 - D 00452 00453 --------------------------------------------------------------------------*/ 00454 00455 typedef struct { 00456 unsigned constraints; 00457 unsigned de, d1; 00458 real *x[2], *fd1[2]; 00459 } opt_edge_data_2; 00460 00461 typedef struct { 00462 unsigned constraints; 00463 real x[2], jac[4]; 00464 } opt_point_data_2; 00465 00466 typedef struct { 00467 lagrange_data *ld; 00468 unsigned size[3]; 00469 const real *elx[2]; 00470 opt_edge_data_2 ed; 00471 opt_point_data_2 pd; 00472 real *work; 00473 real x[2], jac[4]; 00474 } opt_data_2; 00475 void opt_alloc_2(opt_data_2 *p, lagrange_data *ld); 00476 void opt_free_2(opt_data_2 *p); 00477 double opt_findpt_2(opt_data_2 *p, const real *const elx[2], 00478 const real xstar[2], real r[2], unsigned *constr); 00479 00480 #endif 00481