Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
FindPtFuncs.h
Go to the documentation of this file.
00001 #ifndef FINDPTFUNCS_H
00002 #define FINDPTFUNCS_H
00003 
00004 #include "float.h"
00005 #include "math.h"
00006 
00007 /*======================================================
00008 / from types.h
00009 /======================================================
00010 */
00011 
00012 #define MOAB_POLY_EPS ( 128 * DBL_EPSILON )
00013 #define MOAB_POLY_PI  3.1415926535897932384626433832795028841971693993751058209749445923
00014 
00015 #define mbsqrt  sqrt
00016 #define mbabs   fabs
00017 #define mbcos   cos
00018 #define mbsin   sin
00019 #define mbfloor floor
00020 #define mbceil  ceil
00021 
00022 #ifdef INTEGER
00023 #undef INTEGER
00024 #endif
00025 
00026 #ifdef real
00027 #undef real
00028 #endif
00029 
00030 /* integer type to use for everything */
00031 #if defined( USE_LONG )
00032 #define INTEGER long
00033 #elif defined( USE_LONG_LONG )
00034 #define INTEGER long long
00035 #elif defined( USE_SHORT )
00036 #define INTEGER short
00037 #else
00038 #define INTEGER int
00039 #endif
00040 
00041 /* when defined, use the given type for global indices instead of INTEGER */
00042 #if defined( USE_GLOBAL_LONG_LONG )
00043 #define GLOBAL_INT long long
00044 #elif defined( USE_GLOBAL_LONG )
00045 #define GLOBAL_INT long
00046 #else
00047 #define GLOBAL_INT long
00048 #endif
00049 
00050 /* floating point type to use for everything */
00051 #if defined( USE_FLOAT )
00052 typedef float realType;
00053 #elif defined( USE_LONG_DOUBLE )
00054 typedef long double realType;
00055 #else
00056 typedef double realType;
00057 #endif
00058 
00059 /* apparently uint and ulong can be defined already in standard headers */
00060 #ifndef uint
00061 typedef signed INTEGER sint;
00062 #endif
00063 
00064 #ifndef uint
00065 typedef unsigned INTEGER uint;
00066 #endif
00067 #undef INTEGER
00068 
00069 #ifndef slong
00070 #ifdef GLOBAL_INT
00071 typedef signed GLOBAL_INT slong;
00072 #else
00073 typedef sint slong;
00074 #endif
00075 #endif
00076 
00077 #ifndef ulong
00078 #ifdef GLOBAL_INT
00079 typedef unsigned GLOBAL_INT ulong;
00080 #else
00081 typedef uint ulong;
00082 #endif
00083 #endif
00084 
00085 /*======================================================
00086 / from poly.h
00087 /======================================================
00088 */
00089 
00090 /* 
00091   For brevity's sake, some names have been shortened
00092   Quadrature rules
00093     Gauss   -> Gauss-Legendre quadrature (open)
00094     Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends)
00095   Polynomial bases
00096     Legendre -> Legendre basis
00097     Gauss    -> Lagrangian basis using Gauss   quadrature nodes
00098     Lobatto  -> Lagrangian basis using Lobatto quadrature nodes
00099 */
00100 
00101 /*--------------------------------------------------------------------------
00102    Legendre Polynomial Matrix Computation
00103    (compute P_i(x_j) for i = 0, ..., n and a given set of x)
00104   --------------------------------------------------------------------------*/
00105 
00106 /* precondition: n >= 1
00107    inner index is x index (0 ... m-1);
00108    outer index is Legendre polynomial number (0 ... n)
00109  */
00110 void legendre_matrix( const realType* x, int m, realType* P, int n );
00111 
00112 /* precondition: n >= 1
00113    inner index is Legendre polynomial number (0 ... n)
00114    outer index is x index (0 ... m-1);
00115  */
00116 void legendre_matrix_t( const realType* x, int m, realType* P, int n );
00117 
00118 /* precondition: n >= 1
00119    compute P_i(x) with i = 0 ... n
00120  */
00121 void legendre_row( realType x, realType* P, int n );
00122 
00123 /*--------------------------------------------------------------------------
00124    Quadrature Nodes and Weights Calculation
00125    
00126    call the _nodes function before calling the _weights function
00127   --------------------------------------------------------------------------*/
00128 
00129 void gauss_nodes( realType* z, int n );   /* n nodes (order = 2n-1) */
00130 void lobatto_nodes( realType* z, int n ); /* n nodes (order = 2n-3) */
00131 
00132 void gauss_weights( const realType* z, realType* w, int n );
00133 void lobatto_weights( const realType* z, realType* w, int n );
00134 
00135 /*--------------------------------------------------------------------------
00136    Lagrangian to Legendre Change-of-basis Matrix
00137   --------------------------------------------------------------------------*/
00138 
00139 /* precondition: n >= 2
00140    given the Gauss quadrature rule (z,w,n), compute the square matrix J
00141    for transforming from the Gauss basis to the Legendre basis:
00142    
00143       u_legendre(i) = sum_j J(i,j) u_gauss(j)
00144 
00145    computes J   = .5 (2i+1) w  P (z )
00146              ij              j  i  j
00147              
00148    in column major format (inner index is i, the Legendre index)
00149  */
00150 void gauss_to_legendre( const realType* z, const realType* w, int n, realType* J );
00151 
00152 /* precondition: n >= 2
00153    same as above, but
00154    in row major format (inner index is j, the Gauss index)
00155  */
00156 void gauss_to_legendre_t( const realType* z, const realType* w, int n, realType* J );
00157 
00158 /* precondition: n >= 3
00159    given the Lobatto quadrature rule (z,w,n), compute the square matrix J
00160    for transforming from the Lobatto basis to the Legendre basis:
00161    
00162       u_legendre(i) = sum_j J(i,j) u_lobatto(j)
00163 
00164    in column major format (inner index is i, the Legendre index)
00165  */
00166 void lobatto_to_legendre( const realType* z, const realType* w, int n, realType* J );
00167 
00168 /*--------------------------------------------------------------------------
00169    Lagrangian basis function evaluation
00170   --------------------------------------------------------------------------*/
00171 
00172 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
00173    evaluate all Lagrangian basis functions at all points x
00174    
00175    inner index of output J is the basis function index (row-major format)
00176    provide work array with space for 4*n doubles
00177  */
00178 void lagrange_weights( const realType* z, unsigned n, const realType* x, unsigned m, realType* J, realType* work );
00179 
00180 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
00181    evaluate all Lagrangian basis functions and their derivatives
00182    
00183    inner index of outputs J,D is the basis function index (row-major format)
00184    provide work array with space for 6*n doubles
00185  */
00186 void lagrange_weights_deriv( const realType* z,
00187                              unsigned n,
00188                              const realType* x,
00189                              unsigned m,
00190                              realType* J,
00191                              realType* D,
00192                              realType* work );
00193 
00194 /*--------------------------------------------------------------------------
00195    Speedy Lagrangian Interpolation
00196    
00197    Usage:
00198    
00199      lagrange_data p;
00200      lagrange_setup(&p,z,n);    *  setup for nodes z[0 ... n-1] *
00201      
00202      the weights
00203        p->J [0 ... n-1]     interpolation weights
00204        p->D [0 ... n-1]     1st derivative weights
00205        p->D2[0 ... n-1]     2nd derivative weights
00206      are computed for a given x with:
00207        lagrange_0(p,x);  *  compute p->J *
00208        lagrange_1(p,x);  *  compute p->J, p->D *
00209        lagrange_2(p,x);  *  compute p->J, p->D, p->D2 *
00210        lagrange_2u(p);   *  compute p->D2 after call of lagrange_1(p,x); *
00211      These functions use the z array supplied to setup
00212        (that pointer should not be freed between calls)
00213      Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
00214        p->J_z0, etc. and p->J_zn, etc.
00215 
00216      lagrange_free(&p);  *  deallocate memory allocated by setup *
00217   --------------------------------------------------------------------------*/
00218 
00219 typedef struct
00220 {
00221     unsigned n;                                    /* number of Lagrange nodes            */
00222     const realType* z;                             /* Lagrange nodes (user-supplied)      */
00223     realType *J, *D, *D2;                          /* weights for 0th,1st,2nd derivatives */
00224     realType *J_z0, *D_z0, *D2_z0;                 /* ditto at z[0]   (computed at setup) */
00225     realType *J_zn, *D_zn, *D2_zn;                 /* ditto at z[n-1] (computed at setup) */
00226     realType *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data           */
00227 } lagrange_data;
00228 
00229 void lagrange_setup( lagrange_data* p, const realType* z, unsigned n );
00230 void lagrange_free( lagrange_data* p );
00231 
00232 void lagrange_0( lagrange_data* p, realType x );
00233 void lagrange_1( lagrange_data* p, realType x );
00234 void lagrange_2( lagrange_data* p, realType x );
00235 void lagrange_2u( lagrange_data* p );
00236 
00237 /*======================================================
00238 / from tensor.h
00239 /======================================================
00240 */
00241 
00242 /*--------------------------------------------------------------------------
00243    1-,2-,3-d Tensor Application
00244    
00245    the 3d case:
00246    tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2)
00247      gives v = [ R (x) S (x) T ] u
00248      where R is mr x nr, S is ms x ns, T is mt x nt,
00249        each in row- or column-major format according to f := r | c
00250      u is nr x ns x nt in column-major format (inner index is r)
00251      v is mr x ms x mt in column-major format (inner index is r)
00252   --------------------------------------------------------------------------*/
00253 
00254 void tensor_c1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v );
00255 void tensor_r1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v );
00256 
00257 /* work holds mr*ns realTypes */
00258 void tensor_c2( const realType* R,
00259                 unsigned mr,
00260                 unsigned nr,
00261                 const realType* S,
00262                 unsigned ms,
00263                 unsigned ns,
00264                 const realType* u,
00265                 realType* v,
00266                 realType* work );
00267 void tensor_r2( const realType* R,
00268                 unsigned mr,
00269                 unsigned nr,
00270                 const realType* S,
00271                 unsigned ms,
00272                 unsigned ns,
00273                 const realType* u,
00274                 realType* v,
00275                 realType* work );
00276 
00277 /* work1 holds mr*ns*nt realTypes,
00278    work2 holds mr*ms*nt realTypes */
00279 void tensor_c3( const realType* R,
00280                 unsigned mr,
00281                 unsigned nr,
00282                 const realType* S,
00283                 unsigned ms,
00284                 unsigned ns,
00285                 const realType* T,
00286                 unsigned mt,
00287                 unsigned nt,
00288                 const realType* u,
00289                 realType* v,
00290                 realType* work1,
00291                 realType* work2 );
00292 void tensor_r3( const realType* R,
00293                 unsigned mr,
00294                 unsigned nr,
00295                 const realType* S,
00296                 unsigned ms,
00297                 unsigned ns,
00298                 const realType* T,
00299                 unsigned mt,
00300                 unsigned nt,
00301                 const realType* u,
00302                 realType* v,
00303                 realType* work1,
00304                 realType* work2 );
00305 
00306 /*--------------------------------------------------------------------------
00307    1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
00308    
00309    the 3d case:
00310    v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
00311    same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
00312      gives v = [ Jr (x) Js (x) Jt ] u
00313      where Jr, Js, Jt are row vectors (interpolation weights)
00314      u is nr x ns x nt in column-major format (inner index is r)
00315      v is a scalar
00316   --------------------------------------------------------------------------*/
00317 
00318 realType tensor_i1( const realType* Jr, unsigned nr, const realType* u );
00319 
00320 /* work holds ns realTypes */
00321 realType tensor_i2( const realType* Jr,
00322                     unsigned nr,
00323                     const realType* Js,
00324                     unsigned ns,
00325                     const realType* u,
00326                     realType* work );
00327 
00328 /* work holds ns*nt + nt realTypes */
00329 realType tensor_i3( const realType* Jr,
00330                     unsigned nr,
00331                     const realType* Js,
00332                     unsigned ns,
00333                     const realType* Jt,
00334                     unsigned nt,
00335                     const realType* u,
00336                     realType* work );
00337 
00338 /*--------------------------------------------------------------------------
00339    1-,2-,3-d Tensor Application of Row Vectors
00340              for simultaneous Interpolation and Gradient computation
00341    
00342    the 3d case:
00343    v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
00344      gives v   = [ Jr (x) Js (x) Jt ] u
00345            g_0 = [ Dr (x) Js (x) Jt ] u
00346            g_1 = [ Jr (x) Ds (x) Jt ] u
00347            g_2 = [ Jr (x) Js (x) Dt ] u
00348      where Jr,Dr,Js,Ds,Jt,Dt are row vectors
00349        (interpolation & derivative weights)
00350      u is nr x ns x nt in column-major format (inner index is r)
00351      v is a scalar, g is an array of 3 realTypes
00352   --------------------------------------------------------------------------*/
00353 
00354 realType tensor_ig1( const realType* Jr, const realType* Dr, unsigned nr, const realType* u, realType* g );
00355 
00356 /* work holds 2*ns realTypes */
00357 realType tensor_ig2( const realType* Jr,
00358                      const realType* Dr,
00359                      unsigned nr,
00360                      const realType* Js,
00361                      const realType* Ds,
00362                      unsigned ns,
00363                      const realType* u,
00364                      realType* g,
00365                      realType* work );
00366 
00367 /* work holds 2*ns*nt + 3*ns realTypes */
00368 realType tensor_ig3( const realType* Jr,
00369                      const realType* Dr,
00370                      unsigned nr,
00371                      const realType* Js,
00372                      const realType* Ds,
00373                      unsigned ns,
00374                      const realType* Jt,
00375                      const realType* Dt,
00376                      unsigned nt,
00377                      const realType* u,
00378                      realType* g,
00379                      realType* work );
00380 
00381 /*======================================================
00382 / from findpt.h
00383 /======================================================
00384 */
00385 
00386 typedef struct
00387 {
00388     realType x[2], A[4], axis_bnd[4];
00389 } obbox_2;
00390 
00391 typedef struct
00392 {
00393     realType x[3], A[9], axis_bnd[6];
00394 } obbox_3;
00395 
00396 typedef struct
00397 {
00398     unsigned hash_n;
00399     realType bnd[4]; /* bounds for all elements */
00400     realType fac[2]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
00401     obbox_2* obb;    /* obb[nel] -- bounding box info for each element */
00402     uint* offset;    /* hash table -- for cell i,j:
00403                          uint index = j*hash_n+i,
00404                                   b = offset[index  ],
00405                                   e = offset[index+1];
00406                          elements in cell are
00407                            offset[b], offset[b+1], ..., offset[e-1] */
00408     unsigned max;    /* maximum # of elements in any cell */
00409 } hash_data_2;
00410 
00411 typedef struct
00412 {
00413     unsigned hash_n;
00414     realType bnd[6]; /* bounds for all elements */
00415     realType fac[3]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
00416     obbox_3* obb;    /* obb[nel] -- bounding box info for each element */
00417     uint* offset;    /* hash table -- for cell i,j,k:
00418                          uint index = (k*hash_n+j)*hash_n+i,
00419                                   b = offset[index  ],
00420                                   e = offset[index+1];
00421                          elements in cell are
00422                            offset[b], offset[b+1], ..., offset[e-1] */
00423     unsigned max;    /* maximum # of elements in any cell */
00424 } hash_data_3;
00425 
00426 typedef struct
00427 {
00428     uint el;
00429     realType r[3];
00430     realType dist;
00431 } findpt_listel;
00432 
00433 typedef struct
00434 {
00435     unsigned constraints;
00436     unsigned de, d1;
00437     realType *x[2], *fd1[2];
00438 } opt_edge_data_2;
00439 
00440 typedef struct
00441 {
00442     unsigned constraints;
00443     realType x[2], jac[4];
00444 } opt_point_data_2;
00445 
00446 typedef struct
00447 {
00448     lagrange_data* ld;
00449     unsigned size[3];
00450     const realType* elx[2];
00451     opt_edge_data_2 ed;
00452     opt_point_data_2 pd;
00453     realType* work;
00454     realType x[2], jac[4];
00455 } opt_data_2;
00456 
00457 typedef struct
00458 {
00459     const realType* xw[2]; /* geometry data */
00460     realType* z[2];        /* lobatto nodes */
00461     lagrange_data ld[2];   /* interpolation, derivative weights & data */
00462     unsigned nptel;        /* nodes per element */
00463     hash_data_2* hash;     /* geometric hashing data */
00464     findpt_listel *list, **sorted, **end;
00465     /* pre-allocated list of elements to
00466                                            check (found by hashing), and
00467                                            pre-allocated list of pointers into
00468                                            the first list for sorting */
00469     opt_data_2* od; /* data for the optimization algorithm */
00470     realType* od_work;
00471 } findpt_data_2;
00472 
00473 typedef struct
00474 {
00475     unsigned constraints;
00476     unsigned dn, d1, d2;
00477     realType *x[3], *fdn[3];
00478 } opt_face_data_3;
00479 
00480 typedef struct
00481 {
00482     unsigned constraints;
00483     unsigned de, d1, d2;
00484     realType *x[3], *fd1[3], *fd2[3];
00485 } opt_edge_data_3;
00486 
00487 typedef struct
00488 {
00489     unsigned constraints;
00490     realType x[3], jac[9];
00491 } opt_point_data_3;
00492 
00493 typedef struct
00494 {
00495     lagrange_data* ld;
00496     unsigned size[4];
00497     const realType* elx[3];
00498     opt_face_data_3 fd;
00499     opt_edge_data_3 ed;
00500     opt_point_data_3 pd;
00501     realType* work;
00502     realType x[3], jac[9];
00503 } opt_data_3;
00504 
00505 typedef struct
00506 {
00507     const realType* xw[3]; /* geometry data */
00508     realType* z[3];        /* lobatto nodes */
00509     lagrange_data ld[3];   /* interpolation, derivative weights & data */
00510     unsigned nptel;        /* nodes per element */
00511     hash_data_3* hash;     /* geometric hashing data */
00512     findpt_listel *list, **sorted, **end;
00513     /* pre-allocated list of elements to
00514                                            check (found by hashing), and
00515                                            pre-allocated list of pointers into
00516                                            the first list for sorting */
00517     opt_data_3* od; /* data for the optimization algorithm */
00518     realType* od_work;
00519 } findpt_data_3;
00520 
00521 findpt_data_2* findpt_setup_2( const realType* const xw[2],
00522                                const unsigned n[2],
00523                                uint nel,
00524                                uint max_hash_size,
00525                                realType bbox_tol );
00526 findpt_data_3* findpt_setup_3( const realType* const xw[3],
00527                                const unsigned n[3],
00528                                uint nel,
00529                                uint max_hash_size,
00530                                realType bbox_tol );
00531 
00532 void findpt_free_2( findpt_data_2* p );
00533 void findpt_free_3( findpt_data_3* p );
00534 
00535 const realType* findpt_allbnd_2( const findpt_data_2* p );
00536 const realType* findpt_allbnd_3( const findpt_data_3* p );
00537 
00538 typedef int ( *findpt_func )( void*, const realType*, int, uint*, realType*, realType* );
00539 int findpt_2( findpt_data_2* p, const realType x[2], int guess, uint* el, realType r[2], realType* dist );
00540 int findpt_3( findpt_data_3* p, const realType x[3], int guess, uint* el, realType r[3], realType* dist );
00541 
00542 void findpt_weights_2( findpt_data_2* p, const realType r[2] );
00543 
00544 void findpt_weights_3( findpt_data_3* p, const realType r[3] );
00545 
00546 double findpt_eval_2( findpt_data_2* p, const realType* u );
00547 
00548 double findpt_eval_3( findpt_data_3* p, const realType* u );
00549 
00550 /*======================================================
00551 / from extrafindpt.h
00552 /======================================================
00553 */
00554 
00555 void opt_alloc_3( opt_data_3* p, lagrange_data* ld );
00556 void opt_free_3( opt_data_3* p );
00557 double opt_findpt_3( opt_data_3* p,
00558                      const realType* const elx[3],
00559                      const realType xstar[3],
00560                      realType r[3],
00561                      unsigned* constr );
00562 void opt_vol_set_intp_3( opt_data_3* p, const realType r[3] );
00563 
00564 /* for 2d spectralQuad */
00565 /*--------------------------------------------------------------------------
00566 
00567    2 - D
00568 
00569   --------------------------------------------------------------------------*/
00570 
00571 void opt_alloc_2( opt_data_2* p, lagrange_data* ld );
00572 void opt_free_2( opt_data_2* p );
00573 double opt_findpt_2( opt_data_2* p,
00574                      const realType* const elx[2],
00575                      const realType xstar[2],
00576                      realType r[2],
00577                      unsigned* constr );
00578 
00579 extern const unsigned opt_no_constraints_2;
00580 extern const unsigned opt_no_constraints_3;
00581 
00582 /*======================================================
00583 / from errmem.h
00584 /======================================================
00585 */
00586 
00587 /* requires:
00588      <stdlib.h> for malloc, calloc, realloc, free
00589 */
00590 
00591 /*--------------------------------------------------------------------------
00592    Error Reporting
00593    Memory Allocation Wrappers to Catch Out-of-memory
00594   --------------------------------------------------------------------------*/
00595 
00596 /* #include "malloc.h" */
00597 #include <stdlib.h>
00598 
00599 #ifdef __GNUC__
00600 void fail( const char* fmt, ... ) __attribute__( ( noreturn ) );
00601 #define MAYBE_UNUSED __attribute__( ( unused ) )
00602 #else
00603 void fail( const char* fmt, ... );
00604 #define MAYBE_UNUSED
00605 #endif
00606 
00607 #if 0
00608 {}
00609 #endif
00610 
00611 static void* smalloc( size_t size, const char* file ) MAYBE_UNUSED;
00612 static void* smalloc( size_t size, const char* file )
00613 {
00614     void* res = malloc( size );
00615     if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
00616     return res;
00617 }
00618 
00619 static void* scalloc( size_t nmemb, size_t size, const char* file ) MAYBE_UNUSED;
00620 static void* scalloc( size_t nmemb, size_t size, const char* file )
00621 {
00622     void* res = calloc( nmemb, size );
00623     if( !res && nmemb ) fail( "%s: allocation of %d bytes failed\n", file, (int)size * nmemb );
00624     return res;
00625 }
00626 
00627 static void* srealloc( void* ptr, size_t size, const char* file ) MAYBE_UNUSED;
00628 static void* srealloc( void* ptr, size_t size, const char* file )
00629 {
00630     void* res = realloc( ptr, size );
00631     if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
00632     return res;
00633 }
00634 
00635 #define tmalloc( type, count )       ( (type*)smalloc( ( count ) * sizeof( type ), __FILE__ ) )
00636 #define tcalloc( type, count )       ( (type*)scalloc( ( count ), sizeof( type ), __FILE__ ) )
00637 #define trealloc( type, ptr, count ) ( (type*)srealloc( ( ptr ), ( count ) * sizeof( type ), __FILE__ ) )
00638 /*
00639 typedef struct { size_t size; void *ptr; } buffer;
00640 
00641 static void buffer_init_(buffer *b, size_t size, const char *file) MAYBE_UNUSED;
00642 static void buffer_init_(buffer *b, size_t size, const char *file)
00643 {
00644   b->size=size, b->ptr=smalloc(size,file);
00645 }
00646 static void buffer_reserve_(buffer *b, size_t min, const char *file)
00647   MAYBE_UNUSED;
00648 static void buffer_reserve_(buffer *b, size_t min, const char *file)
00649 {
00650   size_t size = b->size;
00651   if(size<min) {
00652     size+=size/2+1;
00653     if(size<min) size=min;
00654     b->ptr=srealloc(b->ptr,size,file);
00655   }
00656 }
00657 static void buffer_free(buffer *b) MAYBE_UNUSED;
00658 static void buffer_free(buffer *b) { free(b->ptr); }
00659 
00660 #define buffer_init(b,size) buffer_init_(b,size,__FILE__)
00661 #define buffer_reserve(b,min) buffer_reserve_(b,min,__FILE__)
00662 */
00663 
00664 /*======================================================
00665 / from minmax.h
00666 /======================================================
00667 */
00668 
00669 /*--------------------------------------------------------------------------
00670    Min, Max, Norm
00671   --------------------------------------------------------------------------*/
00672 
00673 #ifdef __GNUC__
00674 #define MAYBE_UNUSED __attribute__( ( unused ) )
00675 #else
00676 #define MAYBE_UNUSED
00677 #endif
00678 
00679 #define DECLMINMAX( type, prefix )                                                             \
00680     static type prefix##min_2( type a, type b ) MAYBE_UNUSED;                                  \
00681     static type prefix##min_2( type a, type b )                                                \
00682     {                                                                                          \
00683         return b < a ? b : a;                                                                  \
00684     }                                                                                          \
00685     static type prefix##max_2( type a, type b ) MAYBE_UNUSED;                                  \
00686     static type prefix##max_2( type a, type b )                                                \
00687     {                                                                                          \
00688         return a > b ? a : b;                                                                  \
00689     }                                                                                          \
00690     static void prefix##minmax_2( type* min, type* max, type a, type b ) MAYBE_UNUSED;         \
00691     static void prefix##minmax_2( type* min, type* max, type a, type b )                       \
00692     {                                                                                          \
00693         if( b < a )                                                                            \
00694             *min = b, *max = a;                                                                \
00695         else                                                                                   \
00696             *min = a, *max = b;                                                                \
00697     }                                                                                          \
00698     static type prefix##min_3( type a, type b, type c ) MAYBE_UNUSED;                          \
00699     static type prefix##min_3( type a, type b, type c )                                        \
00700     {                                                                                          \
00701         return b < a ? ( c < b ? c : b ) : ( c < a ? c : a );                                  \
00702     }                                                                                          \
00703     static type prefix##max_3( type a, type b, type c ) MAYBE_UNUSED;                          \
00704     static type prefix##max_3( type a, type b, type c )                                        \
00705     {                                                                                          \
00706         return a > b ? ( a > c ? a : c ) : ( b > c ? b : c );                                  \
00707     }                                                                                          \
00708     static void prefix##minmax_3( type* min, type* max, type a, type b, type c ) MAYBE_UNUSED; \
00709     static void prefix##minmax_3( type* min, type* max, type a, type b, type c )               \
00710     {                                                                                          \
00711         if( b < a )                                                                            \
00712             *min = prefix##min_2( b, c ), *max = prefix##max_2( a, c );                        \
00713         else                                                                                   \
00714             *min = prefix##min_2( a, c ), *max = prefix##max_2( b, c );                        \
00715     }
00716 
00717 DECLMINMAX( int, i )
00718 DECLMINMAX( unsigned, u )
00719 DECLMINMAX( realType, r )
00720 #undef DECLMINMAX
00721 
00722 static realType r1norm_1( realType a ) MAYBE_UNUSED;
00723 static realType r1norm_1( realType a )
00724 {
00725     return mbabs( a );
00726 }
00727 static realType r1norm_2( realType a, realType b ) MAYBE_UNUSED;
00728 static realType r1norm_2( realType a, realType b )
00729 {
00730     return mbabs( a ) + mbabs( b );
00731 }
00732 static realType r1norm_3( realType a, realType b, realType c ) MAYBE_UNUSED;
00733 static realType r1norm_3( realType a, realType b, realType c )
00734 {
00735     return mbabs( a ) + mbabs( b ) + mbabs( c );
00736 }
00737 static realType r2norm_1( realType a ) MAYBE_UNUSED;
00738 static realType r2norm_1( realType a )
00739 {
00740     return mbsqrt( a * a );
00741 }
00742 static realType r2norm_2( realType a, realType b ) MAYBE_UNUSED;
00743 static realType r2norm_2( realType a, realType b )
00744 {
00745     return mbsqrt( a * a + b * b );
00746 }
00747 static realType r2norm_3( realType a, realType b, realType c ) MAYBE_UNUSED;
00748 static realType r2norm_3( realType a, realType b, realType c )
00749 {
00750     return mbsqrt( a * a + b * b + c * c );
00751 }
00752 
00753 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines