Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
FindPtFuncs.h
Go to the documentation of this file.
1 #ifndef FINDPTFUNCS_H
2 #define FINDPTFUNCS_H
3 
4 #include "float.h"
5 #include "math.h"
6 
7 /*======================================================
8 / from types.h
9 /======================================================
10 */
11 
12 #define MOAB_POLY_EPS ( 128 * DBL_EPSILON )
13 #define MOAB_POLY_PI 3.1415926535897932384626433832795028841971693993751058209749445923
14 
15 #define mbsqrt sqrt
16 #define mbabs fabs
17 #define mbcos cos
18 #define mbsin sin
19 #define mbfloor floor
20 #define mbceil ceil
21 
22 #ifdef INTEGER
23 #undef INTEGER
24 #endif
25 
26 #ifdef real
27 #undef real
28 #endif
29 
30 /* integer type to use for everything */
31 #if defined( USE_LONG )
32 #define INTEGER long
33 #elif defined( USE_LONG_LONG )
34 #define INTEGER long long
35 #elif defined( USE_SHORT )
36 #define INTEGER short
37 #else
38 #define INTEGER int
39 #endif
40 
41 /* when defined, use the given type for global indices instead of INTEGER */
42 #if defined( USE_GLOBAL_LONG_LONG )
43 #define GLOBAL_INT long long
44 #elif defined( USE_GLOBAL_LONG )
45 #define GLOBAL_INT long
46 #else
47 #define GLOBAL_INT long
48 #endif
49 
50 /* floating point type to use for everything */
51 #if defined( USE_FLOAT )
52 typedef float realType;
53 #elif defined( USE_LONG_DOUBLE )
54 typedef long double realType;
55 #else
56 typedef double realType;
57 #endif
58 
59 /* apparently uint, sint, and ulong can be defined already in standard headers;
60  they may be typedefs (not macros), so #ifndef may not catch them — suppress
61  the resulting duplicate_typedef warning for compilers that emit it (e.g. NVHPC). */
62 #ifdef __NVCOMPILER
63 #pragma diag_suppress duplicate_typedef
64 #endif
65 
66 #ifndef sint
67 typedef signed INTEGER sint;
68 #endif
69 
70 #ifndef uint
71 typedef unsigned INTEGER uint;
72 #endif
73 #undef INTEGER
74 
75 #ifndef slong
76 #ifdef GLOBAL_INT
77 typedef signed GLOBAL_INT slong;
78 #else
79 typedef sint slong;
80 #endif
81 #endif
82 
83 #ifndef ulong
84 #ifdef GLOBAL_INT
85 typedef unsigned GLOBAL_INT ulong;
86 #else
87 typedef uint ulong;
88 #endif
89 #endif
90 
91 #ifdef __NVCOMPILER
92 #pragma diag_default duplicate_typedef
93 #endif
94 
95 /*======================================================
96 / from poly.h
97 /======================================================
98 */
99 
100 /*
101  For brevity's sake, some names have been shortened
102  Quadrature rules
103  Gauss -> Gauss-Legendre quadrature (open)
104  Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends)
105  Polynomial bases
106  Legendre -> Legendre basis
107  Gauss -> Lagrangian basis using Gauss quadrature nodes
108  Lobatto -> Lagrangian basis using Lobatto quadrature nodes
109 */
110 
111 /*--------------------------------------------------------------------------
112  Legendre Polynomial Matrix Computation
113  (compute P_i(x_j) for i = 0, ..., n and a given set of x)
114  --------------------------------------------------------------------------*/
115 
116 /* precondition: n >= 1
117  inner index is x index (0 ... m-1);
118  outer index is Legendre polynomial number (0 ... n)
119  */
120 void legendre_matrix( const realType* x, int m, realType* P, int n );
121 
122 /* precondition: n >= 1
123  inner index is Legendre polynomial number (0 ... n)
124  outer index is x index (0 ... m-1);
125  */
126 void legendre_matrix_t( const realType* x, int m, realType* P, int n );
127 
128 /* precondition: n >= 1
129  compute P_i(x) with i = 0 ... n
130  */
131 void legendre_row( realType x, realType* P, int n );
132 
133 /*--------------------------------------------------------------------------
134  Quadrature Nodes and Weights Calculation
135 
136  call the _nodes function before calling the _weights function
137  --------------------------------------------------------------------------*/
138 
139 void gauss_nodes( realType* z, int n ); /* n nodes (order = 2n-1) */
140 void lobatto_nodes( realType* z, int n ); /* n nodes (order = 2n-3) */
141 
142 void gauss_weights( const realType* z, realType* w, int n );
143 void lobatto_weights( const realType* z, realType* w, int n );
144 
145 /*--------------------------------------------------------------------------
146  Lagrangian to Legendre Change-of-basis Matrix
147  --------------------------------------------------------------------------*/
148 
149 /* precondition: n >= 2
150  given the Gauss quadrature rule (z,w,n), compute the square matrix J
151  for transforming from the Gauss basis to the Legendre basis:
152 
153  u_legendre(i) = sum_j J(i,j) u_gauss(j)
154 
155  computes J = .5 (2i+1) w P (z )
156  ij j i j
157 
158  in column major format (inner index is i, the Legendre index)
159  */
160 void gauss_to_legendre( const realType* z, const realType* w, int n, realType* J );
161 
162 /* precondition: n >= 2
163  same as above, but
164  in row major format (inner index is j, the Gauss index)
165  */
166 void gauss_to_legendre_t( const realType* z, const realType* w, int n, realType* J );
167 
168 /* precondition: n >= 3
169  given the Lobatto quadrature rule (z,w,n), compute the square matrix J
170  for transforming from the Lobatto basis to the Legendre basis:
171 
172  u_legendre(i) = sum_j J(i,j) u_lobatto(j)
173 
174  in column major format (inner index is i, the Legendre index)
175  */
176 void lobatto_to_legendre( const realType* z, const realType* w, int n, realType* J );
177 
178 /*--------------------------------------------------------------------------
179  Lagrangian basis function evaluation
180  --------------------------------------------------------------------------*/
181 
182 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
183  evaluate all Lagrangian basis functions at all points x
184 
185  inner index of output J is the basis function index (row-major format)
186  provide work array with space for 4*n doubles
187  */
188 void lagrange_weights( const realType* z, unsigned n, const realType* x, unsigned m, realType* J, realType* work );
189 
190 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
191  evaluate all Lagrangian basis functions and their derivatives
192 
193  inner index of outputs J,D is the basis function index (row-major format)
194  provide work array with space for 6*n doubles
195  */
196 void lagrange_weights_deriv( const realType* z,
197  unsigned n,
198  const realType* x,
199  unsigned m,
200  realType* J,
201  realType* D,
202  realType* work );
203 
204 /*--------------------------------------------------------------------------
205  Speedy Lagrangian Interpolation
206 
207  Usage:
208 
209  lagrange_data p;
210  lagrange_setup(&p,z,n); * setup for nodes z[0 ... n-1] *
211 
212  the weights
213  p->J [0 ... n-1] interpolation weights
214  p->D [0 ... n-1] 1st derivative weights
215  p->D2[0 ... n-1] 2nd derivative weights
216  are computed for a given x with:
217  lagrange_0(p,x); * compute p->J *
218  lagrange_1(p,x); * compute p->J, p->D *
219  lagrange_2(p,x); * compute p->J, p->D, p->D2 *
220  lagrange_2u(p); * compute p->D2 after call of lagrange_1(p,x); *
221  These functions use the z array supplied to setup
222  (that pointer should not be freed between calls)
223  Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
224  p->J_z0, etc. and p->J_zn, etc.
225 
226  lagrange_free(&p); * deallocate memory allocated by setup *
227  --------------------------------------------------------------------------*/
228 
229 typedef struct
230 {
231  unsigned n; /* number of Lagrange nodes */
232  const realType* z; /* Lagrange nodes (user-supplied) */
233  realType *J, *D, *D2; /* weights for 0th,1st,2nd derivatives */
234  realType *J_z0, *D_z0, *D2_z0; /* ditto at z[0] (computed at setup) */
235  realType *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */
236  realType *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data */
237 } lagrange_data;
238 
239 void lagrange_setup( lagrange_data* p, const realType* z, unsigned n );
240 void lagrange_free( lagrange_data* p );
241 
242 void lagrange_0( lagrange_data* p, realType x );
243 void lagrange_1( lagrange_data* p, realType x );
244 void lagrange_2( lagrange_data* p, realType x );
245 void lagrange_2u( lagrange_data* p );
246 
247 /*======================================================
248 / from tensor.h
249 /======================================================
250 */
251 
252 /*--------------------------------------------------------------------------
253  1-,2-,3-d Tensor Application
254 
255  the 3d case:
256  tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2)
257  gives v = [ R (x) S (x) T ] u
258  where R is mr x nr, S is ms x ns, T is mt x nt,
259  each in row- or column-major format according to f := r | c
260  u is nr x ns x nt in column-major format (inner index is r)
261  v is mr x ms x mt in column-major format (inner index is r)
262  --------------------------------------------------------------------------*/
263 
264 void tensor_c1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v );
265 void tensor_r1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v );
266 
267 /* work holds mr*ns realTypes */
268 void tensor_c2( const realType* R,
269  unsigned mr,
270  unsigned nr,
271  const realType* S,
272  unsigned ms,
273  unsigned ns,
274  const realType* u,
275  realType* v,
276  realType* work );
277 void tensor_r2( const realType* R,
278  unsigned mr,
279  unsigned nr,
280  const realType* S,
281  unsigned ms,
282  unsigned ns,
283  const realType* u,
284  realType* v,
285  realType* work );
286 
287 /* work1 holds mr*ns*nt realTypes,
288  work2 holds mr*ms*nt realTypes */
289 void tensor_c3( const realType* R,
290  unsigned mr,
291  unsigned nr,
292  const realType* S,
293  unsigned ms,
294  unsigned ns,
295  const realType* T,
296  unsigned mt,
297  unsigned nt,
298  const realType* u,
299  realType* v,
300  realType* work1,
301  realType* work2 );
302 void tensor_r3( const realType* R,
303  unsigned mr,
304  unsigned nr,
305  const realType* S,
306  unsigned ms,
307  unsigned ns,
308  const realType* T,
309  unsigned mt,
310  unsigned nt,
311  const realType* u,
312  realType* v,
313  realType* work1,
314  realType* work2 );
315 
316 /*--------------------------------------------------------------------------
317  1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
318 
319  the 3d case:
320  v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
321  same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
322  gives v = [ Jr (x) Js (x) Jt ] u
323  where Jr, Js, Jt are row vectors (interpolation weights)
324  u is nr x ns x nt in column-major format (inner index is r)
325  v is a scalar
326  --------------------------------------------------------------------------*/
327 
328 realType tensor_i1( const realType* Jr, unsigned nr, const realType* u );
329 
330 /* work holds ns realTypes */
331 realType tensor_i2( const realType* Jr,
332  unsigned nr,
333  const realType* Js,
334  unsigned ns,
335  const realType* u,
336  realType* work );
337 
338 /* work holds ns*nt + nt realTypes */
339 realType tensor_i3( const realType* Jr,
340  unsigned nr,
341  const realType* Js,
342  unsigned ns,
343  const realType* Jt,
344  unsigned nt,
345  const realType* u,
346  realType* work );
347 
348 /*--------------------------------------------------------------------------
349  1-,2-,3-d Tensor Application of Row Vectors
350  for simultaneous Interpolation and Gradient computation
351 
352  the 3d case:
353  v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
354  gives v = [ Jr (x) Js (x) Jt ] u
355  g_0 = [ Dr (x) Js (x) Jt ] u
356  g_1 = [ Jr (x) Ds (x) Jt ] u
357  g_2 = [ Jr (x) Js (x) Dt ] u
358  where Jr,Dr,Js,Ds,Jt,Dt are row vectors
359  (interpolation & derivative weights)
360  u is nr x ns x nt in column-major format (inner index is r)
361  v is a scalar, g is an array of 3 realTypes
362  --------------------------------------------------------------------------*/
363 
364 realType tensor_ig1( const realType* Jr, const realType* Dr, unsigned nr, const realType* u, realType* g );
365 
366 /* work holds 2*ns realTypes */
367 realType tensor_ig2( const realType* Jr,
368  const realType* Dr,
369  unsigned nr,
370  const realType* Js,
371  const realType* Ds,
372  unsigned ns,
373  const realType* u,
374  realType* g,
375  realType* work );
376 
377 /* work holds 2*ns*nt + 3*ns realTypes */
378 realType tensor_ig3( const realType* Jr,
379  const realType* Dr,
380  unsigned nr,
381  const realType* Js,
382  const realType* Ds,
383  unsigned ns,
384  const realType* Jt,
385  const realType* Dt,
386  unsigned nt,
387  const realType* u,
388  realType* g,
389  realType* work );
390 
391 /*======================================================
392 / from findpt.h
393 /======================================================
394 */
395 
396 typedef struct
397 {
398  realType x[2], A[4], axis_bnd[4];
399 } obbox_2;
400 
401 typedef struct
402 {
403  realType x[3], A[9], axis_bnd[6];
404 } obbox_3;
405 
406 typedef struct
407 {
408  unsigned hash_n;
409  realType bnd[4]; /* bounds for all elements */
410  realType fac[2]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
411  obbox_2* obb; /* obb[nel] -- bounding box info for each element */
412  uint* offset; /* hash table -- for cell i,j:
413  uint index = j*hash_n+i,
414  b = offset[index ],
415  e = offset[index+1];
416  elements in cell are
417  offset[b], offset[b+1], ..., offset[e-1] */
418  unsigned max; /* maximum # of elements in any cell */
419 } hash_data_2;
420 
421 typedef struct
422 {
423  unsigned hash_n;
424  realType bnd[6]; /* bounds for all elements */
425  realType fac[3]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
426  obbox_3* obb; /* obb[nel] -- bounding box info for each element */
427  uint* offset; /* hash table -- for cell i,j,k:
428  uint index = (k*hash_n+j)*hash_n+i,
429  b = offset[index ],
430  e = offset[index+1];
431  elements in cell are
432  offset[b], offset[b+1], ..., offset[e-1] */
433  unsigned max; /* maximum # of elements in any cell */
434 } hash_data_3;
435 
436 typedef struct
437 {
438  uint el;
439  realType r[3];
440  realType dist;
441 } findpt_listel;
442 
443 typedef struct
444 {
445  unsigned constraints;
446  unsigned de, d1;
447  realType *x[2], *fd1[2];
448 } opt_edge_data_2;
449 
450 typedef struct
451 {
452  unsigned constraints;
453  realType x[2], jac[4];
454 } opt_point_data_2;
455 
456 typedef struct
457 {
458  lagrange_data* ld;
459  unsigned size[3];
460  const realType* elx[2];
461  opt_edge_data_2 ed;
462  opt_point_data_2 pd;
463  realType* work;
464  realType x[2], jac[4];
465 } opt_data_2;
466 
467 typedef struct
468 {
469  const realType* xw[2]; /* geometry data */
470  realType* z[2]; /* lobatto nodes */
471  lagrange_data ld[2]; /* interpolation, derivative weights & data */
472  unsigned nptel; /* nodes per element */
473  hash_data_2* hash; /* geometric hashing data */
474  findpt_listel *list, **sorted, **end;
475  /* pre-allocated list of elements to
476  check (found by hashing), and
477  pre-allocated list of pointers into
478  the first list for sorting */
479  opt_data_2* od; /* data for the optimization algorithm */
480  realType* od_work;
481 } findpt_data_2;
482 
483 typedef struct
484 {
485  unsigned constraints;
486  unsigned dn, d1, d2;
487  realType *x[3], *fdn[3];
488 } opt_face_data_3;
489 
490 typedef struct
491 {
492  unsigned constraints;
493  unsigned de, d1, d2;
494  realType *x[3], *fd1[3], *fd2[3];
495 } opt_edge_data_3;
496 
497 typedef struct
498 {
499  unsigned constraints;
500  realType x[3], jac[9];
501 } opt_point_data_3;
502 
503 typedef struct
504 {
505  lagrange_data* ld;
506  unsigned size[4];
507  const realType* elx[3];
508  opt_face_data_3 fd;
509  opt_edge_data_3 ed;
510  opt_point_data_3 pd;
511  realType* work;
512  realType x[3], jac[9];
513 } opt_data_3;
514 
515 typedef struct
516 {
517  const realType* xw[3]; /* geometry data */
518  realType* z[3]; /* lobatto nodes */
519  lagrange_data ld[3]; /* interpolation, derivative weights & data */
520  unsigned nptel; /* nodes per element */
521  hash_data_3* hash; /* geometric hashing data */
522  findpt_listel *list, **sorted, **end;
523  /* pre-allocated list of elements to
524  check (found by hashing), and
525  pre-allocated list of pointers into
526  the first list for sorting */
527  opt_data_3* od; /* data for the optimization algorithm */
528  realType* od_work;
529 } findpt_data_3;
530 
531 findpt_data_2* findpt_setup_2( const realType* const xw[2],
532  const unsigned n[2],
533  uint nel,
534  uint max_hash_size,
535  realType bbox_tol );
536 findpt_data_3* findpt_setup_3( const realType* const xw[3],
537  const unsigned n[3],
538  uint nel,
539  uint max_hash_size,
540  realType bbox_tol );
541 
542 void findpt_free_2( findpt_data_2* p );
543 void findpt_free_3( findpt_data_3* p );
544 
545 const realType* findpt_allbnd_2( const findpt_data_2* p );
546 const realType* findpt_allbnd_3( const findpt_data_3* p );
547 
548 typedef int ( *findpt_func )( void*, const realType*, int, uint*, realType*, realType* );
549 int findpt_2( findpt_data_2* p, const realType x[2], int guess, uint* el, realType r[2], realType* dist );
550 int findpt_3( findpt_data_3* p, const realType x[3], int guess, uint* el, realType r[3], realType* dist );
551 
552 void findpt_weights_2( findpt_data_2* p, const realType r[2] );
553 
554 void findpt_weights_3( findpt_data_3* p, const realType r[3] );
555 
556 double findpt_eval_2( findpt_data_2* p, const realType* u );
557 
558 double findpt_eval_3( findpt_data_3* p, const realType* u );
559 
560 /*======================================================
561 / from extrafindpt.h
562 /======================================================
563 */
564 
565 void opt_alloc_3( opt_data_3* p, lagrange_data* ld );
566 void opt_free_3( opt_data_3* p );
567 double opt_findpt_3( opt_data_3* p,
568  const realType* const elx[3],
569  const realType xstar[3],
570  realType r[3],
571  unsigned* constr );
572 void opt_vol_set_intp_3( opt_data_3* p, const realType r[3] );
573 
574 /* for 2d spectralQuad */
575 /*--------------------------------------------------------------------------
576 
577  2 - D
578 
579  --------------------------------------------------------------------------*/
580 
581 void opt_alloc_2( opt_data_2* p, lagrange_data* ld );
582 void opt_free_2( opt_data_2* p );
583 double opt_findpt_2( opt_data_2* p,
584  const realType* const elx[2],
585  const realType xstar[2],
586  realType r[2],
587  unsigned* constr );
588 
589 extern const unsigned opt_no_constraints_2;
590 extern const unsigned opt_no_constraints_3;
591 
592 /*======================================================
593 / from errmem.h
594 /======================================================
595 */
596 
597 /* requires:
598  <stdlib.h> for malloc, calloc, realloc, free
599 */
600 
601 /*--------------------------------------------------------------------------
602  Error Reporting
603  Memory Allocation Wrappers to Catch Out-of-memory
604  --------------------------------------------------------------------------*/
605 
606 /* #include "malloc.h" */
607 #include <stdlib.h>
608 
609 #ifdef __GNUC__
610 void fail( const char* fmt, ... ) __attribute__( ( noreturn ) );
611 #define MAYBE_UNUSED __attribute__( ( unused ) )
612 #else
613 void fail( const char* fmt, ... );
614 #define MAYBE_UNUSED
615 #endif
616 
617 #if 0
618 {}
619 #endif
620 
621 static void* smalloc( size_t size, const char* file ) MAYBE_UNUSED;
622 static void* smalloc( size_t size, const char* file )
623 {
624  void* res = malloc( size );
625  if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
626  return res;
627 }
628 
629 static void* scalloc( size_t nmemb, size_t size, const char* file ) MAYBE_UNUSED;
630 static void* scalloc( size_t nmemb, size_t size, const char* file )
631 {
632  void* res = calloc( nmemb, size );
633  if( !res && nmemb ) fail( "%s: allocation of %d bytes failed\n", file, (int)size * nmemb );
634  return res;
635 }
636 
637 static void* srealloc( void* ptr, size_t size, const char* file ) MAYBE_UNUSED;
638 static void* srealloc( void* ptr, size_t size, const char* file )
639 {
640  void* res = realloc( ptr, size );
641  if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
642  return res;
643 }
644 
645 #define tmalloc( type, count ) ( (type*)smalloc( ( count ) * sizeof( type ), __FILE__ ) )
646 #define tcalloc( type, count ) ( (type*)scalloc( ( count ), sizeof( type ), __FILE__ ) )
647 #define trealloc( type, ptr, count ) ( (type*)srealloc( ( ptr ), ( count ) * sizeof( type ), __FILE__ ) )
648 /*
649 typedef struct { size_t size; void *ptr; } buffer;
650 
651 static void buffer_init_(buffer *b, size_t size, const char *file) MAYBE_UNUSED;
652 static void buffer_init_(buffer *b, size_t size, const char *file)
653 {
654  b->size=size, b->ptr=smalloc(size,file);
655 }
656 static void buffer_reserve_(buffer *b, size_t min, const char *file)
657  MAYBE_UNUSED;
658 static void buffer_reserve_(buffer *b, size_t min, const char *file)
659 {
660  size_t size = b->size;
661  if(size<min) {
662  size+=size/2+1;
663  if(size<min) size=min;
664  b->ptr=srealloc(b->ptr,size,file);
665  }
666 }
667 static void buffer_free(buffer *b) MAYBE_UNUSED;
668 static void buffer_free(buffer *b) { free(b->ptr); }
669 
670 #define buffer_init(b,size) buffer_init_(b,size,__FILE__)
671 #define buffer_reserve(b,min) buffer_reserve_(b,min,__FILE__)
672 */
673 
674 /*======================================================
675 / from minmax.h
676 /======================================================
677 */
678 
679 /*--------------------------------------------------------------------------
680  Min, Max, Norm
681  --------------------------------------------------------------------------*/
682 
683 #ifdef __GNUC__
684 #define MAYBE_UNUSED __attribute__( ( unused ) )
685 #else
686 #define MAYBE_UNUSED
687 #endif
688 
689 #define DECLMINMAX( type, prefix ) \
690  static type prefix##min_2( type a, type b ) MAYBE_UNUSED; \
691  static type prefix##min_2( type a, type b ) \
692  { \
693  return b < a ? b : a; \
694  } \
695  static type prefix##max_2( type a, type b ) MAYBE_UNUSED; \
696  static type prefix##max_2( type a, type b ) \
697  { \
698  return a > b ? a : b; \
699  } \
700  static void prefix##minmax_2( type* min, type* max, type a, type b ) MAYBE_UNUSED; \
701  static void prefix##minmax_2( type* min, type* max, type a, type b ) \
702  { \
703  if( b < a ) \
704  *min = b, *max = a; \
705  else \
706  *min = a, *max = b; \
707  } \
708  static type prefix##min_3( type a, type b, type c ) MAYBE_UNUSED; \
709  static type prefix##min_3( type a, type b, type c ) \
710  { \
711  return b < a ? ( c < b ? c : b ) : ( c < a ? c : a ); \
712  } \
713  static type prefix##max_3( type a, type b, type c ) MAYBE_UNUSED; \
714  static type prefix##max_3( type a, type b, type c ) \
715  { \
716  return a > b ? ( a > c ? a : c ) : ( b > c ? b : c ); \
717  } \
718  static void prefix##minmax_3( type* min, type* max, type a, type b, type c ) MAYBE_UNUSED; \
719  static void prefix##minmax_3( type* min, type* max, type a, type b, type c ) \
720  { \
721  if( b < a ) \
722  *min = prefix##min_2( b, c ), *max = prefix##max_2( a, c ); \
723  else \
724  *min = prefix##min_2( a, c ), *max = prefix##max_2( b, c ); \
725  }
726 
727 DECLMINMAX( int, i )
728 DECLMINMAX( unsigned, u )
729 DECLMINMAX( realType, r )
730 #undef DECLMINMAX
731 
732 static realType r1norm_1( realType a ) MAYBE_UNUSED;
733 static realType r1norm_1( realType a )
734 {
735  return mbabs( a );
736 }
737 static realType r1norm_2( realType a, realType b ) MAYBE_UNUSED;
738 static realType r1norm_2( realType a, realType b )
739 {
740  return mbabs( a ) + mbabs( b );
741 }
742 static realType r1norm_3( realType a, realType b, realType c ) MAYBE_UNUSED;
743 static realType r1norm_3( realType a, realType b, realType c )
744 {
745  return mbabs( a ) + mbabs( b ) + mbabs( c );
746 }
747 static realType r2norm_1( realType a ) MAYBE_UNUSED;
748 static realType r2norm_1( realType a )
749 {
750  return mbsqrt( a * a );
751 }
752 static realType r2norm_2( realType a, realType b ) MAYBE_UNUSED;
753 static realType r2norm_2( realType a, realType b )
754 {
755  return mbsqrt( a * a + b * b );
756 }
757 static realType r2norm_3( realType a, realType b, realType c ) MAYBE_UNUSED;
758 static realType r2norm_3( realType a, realType b, realType c )
759 {
760  return mbsqrt( a * a + b * b + c * c );
761 }
762 
763 #endif