FreeFOAM The Cross-Platform CFD Toolkit
geompack.C
Go to the documentation of this file.
1 # include <cstdlib>
2 # include <iostream>
3 # include <iomanip>
4 # include <fstream>
5 # include <cmath>
6 # include <ctime>
7 # include <cstring>
8 
9 using namespace std;
10 
11 # include <meshTools/geompack.H>
12 
13 //******************************************************************************
14 
15 double d_epsilon ( void )
16 
17 //******************************************************************************
18 //
19 // Purpose:
20 //
21 // D_EPSILON returns the round off unit for double precision arithmetic.
22 //
23 // Discussion:
24 //
25 // D_EPSILON is a number R which is a power of 2 with the property that,
26 // to the precision of the computer's arithmetic,
27 // 1 < 1 + R
28 // but
29 // 1 = ( 1 + R / 2 )
30 //
31 // Modified:
32 //
33 // 06 May 2003
34 //
35 // Author:
36 //
37 // John Burkardt
38 //
39 // Parameters:
40 //
41 // Output, double D_EPSILON, the floating point round-off unit.
42 //
43 {
44  double r = 1.0;
45 
46  while (1.0 < 1.0 + r)
47  {
48  r = r/2.0;
49  }
50 
51  return 2.0*r;
52 }
53 
54 
55 //*********************************************************************
56 
57 double d_max ( double x, double y )
58 
59 //*********************************************************************
60 //
61 // Purpose:
62 //
63 // D_MAX returns the maximum of two real values.
64 //
65 // Modified:
66 //
67 // 10 January 2002
68 //
69 // Author:
70 //
71 // John Burkardt
72 //
73 // Parameters:
74 //
75 // Input, double X, Y, the quantities to compare.
76 //
77 // Output, double D_MAX, the maximum of X and Y.
78 //
79 {
80  if ( y < x )
81  {
82  return x;
83  }
84  else
85  {
86  return y;
87  }
88 }
89 //*********************************************************************
90 
91 double d_min ( double x, double y )
92 
93 //*********************************************************************
94 //
95 // Purpose:
96 //
97 // D_MIN returns the minimum of two real values.
98 //
99 // Modified:
100 //
101 // 09 May 2003
102 //
103 // Author:
104 //
105 // John Burkardt
106 //
107 // Parameters:
108 //
109 // Input, double X, Y, the quantities to compare.
110 //
111 // Output, double D_MIN, the minimum of X and Y.
112 //
113 {
114  if ( y < x )
115  {
116  return y;
117  }
118  else
119  {
120  return x;
121  }
122 }
123 //******************************************************************************
124 
125 void d2vec_part_quick_a ( int n, double a[], int *l, int *r )
126 
127 //******************************************************************************
128 //
129 // Purpose:
130 //
131 // D2VEC_PART_QUICK_A reorders an R2 vector as part of a quick sort.
132 //
133 // Discussion:
134 //
135 // The routine reorders the entries of A. Using A(1:2,1) as a
136 // key, all entries of A that are less than or equal to the key will
137 // precede the key, which precedes all entries that are greater than the key.
138 //
139 // Example:
140 //
141 // Input:
142 //
143 // N = 8
144 //
145 // A = ( (2,4), (8,8), (6,2), (0,2), (10,6), (10,0), (0,6), (4,8) )
146 //
147 // Output:
148 //
149 // L = 2, R = 4
150 //
151 // A = ( (0,2), (0,6), (2,4), (8,8), (6,2), (10,6), (10,0), (4,8) )
152 // ----------- ----------------------------------
153 // LEFT KEY RIGHT
154 //
155 // Modified:
156 //
157 // 01 September 2003
158 //
159 // Author:
160 //
161 // John Burkardt
162 //
163 // Parameters:
164 //
165 // Input, int N, the number of entries of A.
166 //
167 // Input/output, double A[N*2]. On input, the array to be checked.
168 // On output, A has been reordered as described above.
169 //
170 // Output, int *L, *R, the indices of A that define the three segments.
171 // Let KEY = the input value of A(1:2,1). Then
172 // I <= L A(1:2,I) < KEY;
173 // L < I < R A(1:2,I) = KEY;
174 // R <= I A(1:2,I) > KEY.
175 //
176 {
177  int i;
178  int j;
179  double key[2];
180  int ll;
181  int m;
182  int rr;
183 
184  if ( n < 1 )
185  {
186  cout << "\n";
187  cout << "D2VEC_PART_QUICK_A - Fatal error!\n";
188  cout << " N < 1.\n";
189  exit ( 1 );
190  }
191 
192  if ( n == 1 )
193  {
194  *l = 0;
195  *r = 2;
196  return;
197  }
198 
199  key[0] = a[2*0+0];
200  key[1] = a[2*0+1];
201  m = 1;
202 //
203 // The elements of unknown size have indices between L+1 and R-1.
204 //
205  ll = 1;
206  rr = n + 1;
207 
208  for ( i = 2; i <= n; i++ )
209  {
210  if ( dvec_gt ( 2, a+2*ll, key ) )
211  {
212  rr = rr - 1;
213  dvec_swap ( 2, a+2*(rr-1), a+2*ll );
214  }
215  else if ( dvec_eq ( 2, a+2*ll, key ) )
216  {
217  m = m + 1;
218  dvec_swap ( 2, a+2*(m-1), a+2*ll );
219  ll = ll + 1;
220  }
221  else if ( dvec_lt ( 2, a+2*ll, key ) )
222  {
223  ll = ll + 1;
224  }
225 
226  }
227 //
228 // Now shift small elements to the left, and KEY elements to center.
229 //
230  for ( i = 0; i < ll - m; i++ )
231  {
232  for ( j = 0; j < 2; j++ )
233  {
234  a[2*i+j] = a[2*(i+m)+j];
235  }
236  }
237 
238  ll = ll - m;
239 
240  for ( i = ll; i < ll+m; i++ )
241  {
242  for ( j = 0; j < 2; j++ )
243  {
244  a[2*i+j] = key[j];
245  }
246  }
247 
248  *l = ll;
249  *r = rr;
250 
251  return;
252 }
253 //*******************************************************************************
254 
255 void d2vec_permute ( int n, double a[], int p[] )
256 
257 //*******************************************************************************
258 //
259 // Purpose:
260 //
261 // D2VEC_PERMUTE permutes an R2 vector in place.
262 //
263 // Discussion:
264 //
265 // This routine permutes an array of real "objects", but the same
266 // logic can be used to permute an array of objects of any arithmetic
267 // type, or an array of objects of any complexity. The only temporary
268 // storage required is enough to store a single object. The number
269 // of data movements made is N + the number of cycles of order 2 or more,
270 // which is never more than N + N/2.
271 //
272 // Example:
273 //
274 // Input:
275 //
276 // N = 5
277 // P = ( 2, 4, 5, 1, 3 )
278 // A = ( 1.0, 2.0, 3.0, 4.0, 5.0 )
279 // (11.0, 22.0, 33.0, 44.0, 55.0 )
280 //
281 // Output:
282 //
283 // A = ( 2.0, 4.0, 5.0, 1.0, 3.0 )
284 // ( 22.0, 44.0, 55.0, 11.0, 33.0 ).
285 //
286 // Modified:
287 //
288 // 19 February 2004
289 //
290 // Author:
291 //
292 // John Burkardt
293 //
294 // Parameters:
295 //
296 // Input, int N, the number of objects.
297 //
298 // Input/output, double A[2*N], the array to be permuted.
299 //
300 // Input, int P[N], the permutation. P(I) = J means
301 // that the I-th element of the output array should be the J-th
302 // element of the input array. P must be a legal permutation
303 // of the integers from 1 to N, otherwise the algorithm will
304 // fail catastrophically.
305 //
306 {
307  double a_temp[2];
308  int i;
309  int iget;
310  int iput;
311  int istart;
312 
313  if ( !perm_check ( n, p ) )
314  {
315  cout << "\n";
316  cout << "D2VEC_PERMUTE - Fatal error!\n";
317  cout << " The input array does not represent\n";
318  cout << " a proper permutation.\n";
319  exit ( 1 );
320  }
321 //
322 // Search for the next element of the permutation that has not been used.
323 //
324  for ( istart = 1; istart <= n; istart++ )
325  {
326  if ( p[istart-1] < 0 )
327  {
328  continue;
329  }
330  else if ( p[istart-1] == istart )
331  {
332  p[istart-1] = -p[istart-1];
333  continue;
334  }
335  else
336  {
337  a_temp[0] = a[0+(istart-1)*2];
338  a_temp[1] = a[1+(istart-1)*2];
339  iget = istart;
340 //
341 // Copy the new value into the vacated entry.
342 //
343  for ( ; ; )
344  {
345  iput = iget;
346  iget = p[iget-1];
347 
348  p[iput-1] = -p[iput-1];
349 
350  if ( iget < 1 || n < iget )
351  {
352  cout << "\n";
353  cout << "D2VEC_PERMUTE - Fatal error!\n";
354  exit ( 1 );
355  }
356 
357  if ( iget == istart )
358  {
359  a[0+(iput-1)*2] = a_temp[0];
360  a[1+(iput-1)*2] = a_temp[1];
361  break;
362  }
363  a[0+(iput-1)*2] = a[0+(iget-1)*2];
364  a[1+(iput-1)*2] = a[1+(iget-1)*2];
365  }
366  }
367  }
368 //
369 // Restore the signs of the entries.
370 //
371  for ( i = 0; i < n; i++ )
372  {
373  p[i] = -p[i];
374  }
375 
376  return;
377 }
378 //******************************************************************************
379 
380 int *d2vec_sort_heap_index_a ( int n, double a[] )
381 
382 //******************************************************************************
383 //
384 // Purpose:
385 //
386 // D2VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R2 vector.
387 //
388 // Discussion:
389 //
390 // The sorting is not actually carried out. Rather an index array is
391 // created which defines the sorting. This array may be used to sort
392 // or index the array, or to sort or index related arrays keyed on the
393 // original array.
394 //
395 // Once the index array is computed, the sorting can be carried out
396 // "implicitly:
397 //
398 // A(1:2,INDX(I)), I = 1 to N is sorted,
399 //
400 // or explicitly, by the call
401 //
402 // call D2VEC_PERMUTE ( N, A, INDX )
403 //
404 // after which A(1:2,I), I = 1 to N is sorted.
405 //
406 // Modified:
407 //
408 // 13 January 2004
409 //
410 // Author:
411 //
412 // John Burkardt
413 //
414 // Parameters:
415 //
416 // Input, int N, the number of entries in the array.
417 //
418 // Input, double A[2*N], an array to be index-sorted.
419 //
420 // Output, int D2VEC_SORT_HEAP_INDEX_A[N], the sort index. The
421 // I-th element of the sorted array is A(0:1,D2VEC_SORT_HEAP_INDEX_A(I-1)).
422 //
423 {
424  double aval[2];
425  int i;
426  int *indx;
427  int indxt;
428  int ir;
429  int j;
430  int l;
431 
432  if ( n < 1 )
433  {
434  return NULL;
435  }
436 
437  if ( n == 1 )
438  {
439  indx = new int[1];
440  indx[0] = 1;
441  return indx;
442  }
443 
444  indx = ivec_indicator ( n );
445 
446  l = n / 2 + 1;
447  ir = n;
448 
449  for ( ; ; )
450  {
451  if ( 1 < l )
452  {
453  l = l - 1;
454  indxt = indx[l-1];
455  aval[0] = a[0+(indxt-1)*2];
456  aval[1] = a[1+(indxt-1)*2];
457  }
458  else
459  {
460  indxt = indx[ir-1];
461  aval[0] = a[0+(indxt-1)*2];
462  aval[1] = a[1+(indxt-1)*2];
463  indx[ir-1] = indx[0];
464  ir = ir - 1;
465 
466  if ( ir == 1 )
467  {
468  indx[0] = indxt;
469  break;
470  }
471 
472  }
473 
474  i = l;
475  j = l + l;
476 
477  while ( j <= ir )
478  {
479  if ( j < ir )
480  {
481  if ( a[0+(indx[j-1]-1)*2] < a[0+(indx[j]-1)*2] ||
482  ( a[0+(indx[j-1]-1)*2] == a[0+(indx[j]-1)*2] &&
483  a[1+(indx[j-1]-1)*2] < a[1+(indx[j]-1)*2] ) )
484  {
485  j = j + 1;
486  }
487  }
488 
489  if ( aval[0] < a[0+(indx[j-1]-1)*2] ||
490  ( aval[0] == a[0+(indx[j-1]-1)*2] &&
491  aval[1] < a[1+(indx[j-1]-1)*2] ) )
492  {
493  indx[i-1] = indx[j-1];
494  i = j;
495  j = j + j;
496  }
497  else
498  {
499  j = ir + 1;
500  }
501  }
502  indx[i-1] = indxt;
503  }
504 
505  return indx;
506 }
507 //*****************************************************************************
508 
509 void d2vec_sort_quick_a ( int n, double a[] )
510 
511 //*****************************************************************************
512 //
513 // Purpose:
514 //
515 // D2VEC_SORT_QUICK_A ascending sorts an R2 vector using quick sort.
516 //
517 // Discussion:
518 //
519 // The data structure is a set of N pairs of real numbers.
520 // These values are stored in a one dimensional array, by pairs.
521 //
522 // Modified:
523 //
524 // 01 September 2003
525 //
526 // Author:
527 //
528 // John Burkardt
529 //
530 // Parameters:
531 //
532 // Input, int N, the number of entries in the array.
533 //
534 // Input/output, double A[N*2].
535 // On input, the array to be sorted.
536 // On output, the array has been sorted.
537 //
538 {
539 # define LEVEL_MAX 25
540 
541  int base;
542  int l_segment;
543  int level;
544  int n_segment;
545  int rsave[LEVEL_MAX];
546  int r_segment;
547 
548  if ( n < 1 )
549  {
550  cout << "\n";
551  cout << "D2VEC_SORT_QUICK_A - Fatal error!\n";
552  cout << " N < 1.\n";
553  exit ( 1 );
554  }
555 
556  if ( n == 1 )
557  {
558  return;
559  }
560 
561  level = 1;
562  rsave[level-1] = n + 1;
563  base = 1;
564  n_segment = n;
565 
566  while ( 0 < n_segment )
567  {
568 //
569 // Partition the segment.
570 //
571  d2vec_part_quick_a ( n_segment, a+2*(base-1)+0, &l_segment, &r_segment );
572 //
573 // If the left segment has more than one element, we need to partition it.
574 //
575  if ( 1 < l_segment )
576  {
577  if ( LEVEL_MAX < level )
578  {
579  cout << "\n";
580  cout << "D2VEC_SORT_QUICK_A - Fatal error!\n";
581  cout << " Exceeding recursion maximum of " << LEVEL_MAX << "\n";
582  exit ( 1 );
583  }
584 
585  level = level + 1;
586  n_segment = l_segment;
587  rsave[level-1] = r_segment + base - 1;
588  }
589 //
590 // The left segment and the middle segment are sorted.
591 // Must the right segment be partitioned?
592 //
593  else if ( r_segment < n_segment )
594  {
595  n_segment = n_segment + 1 - r_segment;
596  base = base + r_segment - 1;
597  }
598 //
599 // Otherwise, we back up a level if there is an earlier one.
600 //
601  else
602  {
603  for ( ; ; )
604  {
605  if ( level <= 1 )
606  {
607  n_segment = 0;
608  break;
609  }
610 
611  base = rsave[level-1];
612  n_segment = rsave[level-2] - rsave[level-1];
613  level = level - 1;
614 
615  if ( 0 < n_segment )
616  {
617  break;
618  }
619  }
620  }
621  }
622  return;
623 # undef LEVEL_MAX
624 }
625 //******************************************************************************
626 
627 int diaedg ( double x0, double y0, double x1, double y1, double x2, double y2,
628  double x3, double y3 )
629 
630 //******************************************************************************
631 //
632 // Purpose:
633 //
634 // DIAEDG chooses a diagonal edge.
635 //
636 // Discussion:
637 //
638 // The routine determines whether 0--2 or 1--3 is the diagonal edge
639 // that should be chosen, based on the circumcircle criterion, where
640 // (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple
641 // quadrilateral in counterclockwise order.
642 //
643 // Modified:
644 //
645 // 28 August 2003
646 //
647 // Author:
648 //
649 // Barry Joe,
650 // Department of Computing Science,
651 // University of Alberta,
652 // Edmonton, Alberta, Canada T6G 2H1
653 //
654 // Reference:
655 //
656 // Barry Joe,
657 // GEOMPACK - a software package for the generation of meshes
658 // using geometric algorithms,
659 // Advances in Engineering Software,
660 // Volume 13, pages 325-331, 1991.
661 //
662 // Parameters:
663 //
664 // Input, double X0, Y0, X1, Y1, X2, Y2, X3, Y3, the coordinates of the
665 // vertices of a quadrilateral, given in counter clockwise order.
666 //
667 // Output, int DIAEDG, chooses a diagonal:
668 // +1, if diagonal edge 02 is chosen;
669 // -1, if diagonal edge 13 is chosen;
670 // 0, if the four vertices are cocircular.
671 //
672 {
673  double ca;
674  double cb;
675  double dx10;
676  double dx12;
677  double dx30;
678  double dx32;
679  double dy10;
680  double dy12;
681  double dy30;
682  double dy32;
683  double s;
684  double tol;
685  double tola;
686  double tolb;
687  int value;
688 
689  tol = 100.0 * d_epsilon ( );
690 
691  dx10 = x1 - x0;
692  dy10 = y1 - y0;
693  dx12 = x1 - x2;
694  dy12 = y1 - y2;
695  dx30 = x3 - x0;
696  dy30 = y3 - y0;
697  dx32 = x3 - x2;
698  dy32 = y3 - y2;
699 
700  tola = tol * d_max ( fabs ( dx10 ),
701  d_max ( fabs ( dy10 ),
702  d_max ( fabs ( dx30 ), fabs ( dy30 ) ) ) );
703 
704  tolb = tol * d_max ( fabs ( dx12 ),
705  d_max ( fabs ( dy12 ),
706  d_max ( fabs ( dx32 ), fabs ( dy32 ) ) ) );
707 
708  ca = dx10 * dx30 + dy10 * dy30;
709  cb = dx12 * dx32 + dy12 * dy32;
710 
711  if ( tola < ca && tolb < cb )
712  {
713  value = -1;
714  }
715  else if ( ca < -tola && cb < -tolb )
716  {
717  value = 1;
718  }
719  else
720  {
721  tola = d_max ( tola, tolb );
722  s = ( dx10 * dy30 - dx30 * dy10 ) * cb
723  + ( dx32 * dy12 - dx12 * dy32 ) * ca;
724 
725  if ( tola < s )
726  {
727  value = -1;
728  }
729  else if ( s < -tola )
730  {
731  value = 1;
732  }
733  else
734  {
735  value = 0;
736  }
737 
738  }
739 
740  return value;
741 }
742 //******************************************************************************
743 
744 void dmat_transpose_print ( int m, int n, double a[], const char *title )
745 
746 //******************************************************************************
747 //
748 // Purpose:
749 //
750 // DMAT_TRANSPOSE_PRINT prints a real matrix, transposed.
751 //
752 // Modified:
753 //
754 // 11 August 2004
755 //
756 // Author:
757 //
758 // John Burkardt
759 //
760 // Parameters:
761 //
762 // Input, int M, N, the number of rows and columns.
763 //
764 // Input, double A[M*N], an M by N matrix to be printed.
765 //
766 // Input, const char *TITLE, an optional title.
767 //
768 {
769  dmat_transpose_print_some ( m, n, a, 1, 1, m, n, title );
770 
771  return;
772 }
773 //******************************************************************************
774 
775 void dmat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
776  int ihi, int jhi, const char *title )
777 
778 //******************************************************************************
779 //
780 // Purpose:
781 //
782 // DMAT_TRANSPOSE_PRINT_SOME prints some of a real matrix, transposed.
783 //
784 // Modified:
785 //
786 // 11 August 2004
787 //
788 // Author:
789 //
790 // John Burkardt
791 //
792 // Parameters:
793 //
794 // Input, int M, N, the number of rows and columns.
795 //
796 // Input, double A[M*N], an M by N matrix to be printed.
797 //
798 // Input, int ILO, JLO, the first row and column to print.
799 //
800 // Input, int IHI, JHI, the last row and column to print.
801 //
802 // Input, const char *TITLE, an optional title.
803 //
804 {
805 # define INCX 5
806 
807  int i;
808  int i2;
809  int i2hi;
810  int i2lo;
811  int inc;
812  int j;
813  int j2hi;
814  int j2lo;
815 
816  if ( 0 < s_len_trim ( title ) )
817  {
818  cout << "\n";
819  cout << title << "\n";
820  }
821 
822  for ( i2lo = i_max ( ilo, 1 ); i2lo <= i_min ( ihi, m ); i2lo = i2lo + INCX )
823  {
824  i2hi = i2lo + INCX - 1;
825  i2hi = i_min ( i2hi, m );
826  i2hi = i_min ( i2hi, ihi );
827 
828  inc = i2hi + 1 - i2lo;
829 
830  cout << "\n";
831  cout << " Row: ";
832  for ( i = i2lo; i <= i2hi; i++ )
833  {
834  cout << setw(7) << i << " ";
835  }
836  cout << "\n";
837  cout << " Col\n";
838  cout << "\n";
839 
840  j2lo = i_max ( jlo, 1 );
841  j2hi = i_min ( jhi, n );
842 
843  for ( j = j2lo; j <= j2hi; j++ )
844  {
845  cout << setw(5) << j << " ";
846  for ( i2 = 1; i2 <= inc; i2++ )
847  {
848  i = i2lo - 1 + i2;
849  cout << setw(14) << a[(i-1)+(j-1)*m];
850  }
851  cout << "\n";
852  }
853  }
854  cout << "\n";
855 
856  return;
857 # undef INCX
858 }
859 //******************************************************************************
860 
861 void dmat_uniform ( int m, int n, double b, double c, int *seed, double r[] )
862 
863 //******************************************************************************
864 //
865 // Purpose:
866 //
867 // DMAT_UNIFORM fills a double precision array with scaled pseudorandom values.
868 //
869 // Discussion:
870 //
871 // This routine implements the recursion
872 //
873 // seed = 16807 * seed mod ( 2**31 - 1 )
874 // unif = seed / ( 2**31 - 1 )
875 //
876 // The integer arithmetic never requires more than 32 bits,
877 // including a sign bit.
878 //
879 // Modified:
880 //
881 // 30 January 2005
882 //
883 // Author:
884 //
885 // John Burkardt
886 //
887 // Reference:
888 //
889 // Paul Bratley, Bennett Fox, Linus Schrage,
890 // A Guide to Simulation,
891 // Springer Verlag, pages 201-202, 1983.
892 //
893 // Bennett Fox,
894 // Algorithm 647:
895 // Implementation and Relative Efficiency of Quasirandom
896 // Sequence Generators,
897 // ACM Transactions on Mathematical Software,
898 // Volume 12, Number 4, pages 362-376, 1986.
899 //
900 // Peter Lewis, Allen Goodman, James Miller,
901 // A Pseudo-Random Number Generator for the System/360,
902 // IBM Systems Journal,
903 // Volume 8, pages 136-143, 1969.
904 //
905 // Parameters:
906 //
907 // Input, int M, N, the number of rows and columns.
908 //
909 // Input, double B, C, the limits of the pseudorandom values.
910 //
911 // Input/output, int *SEED, the "seed" value. Normally, this
912 // value should not be 0, otherwise the output value of SEED
913 // will still be 0, and D_UNIFORM will be 0. On output, SEED has
914 // been updated.
915 //
916 // Output, double DMAT_UNIFORM[M*N], a matrix of pseudorandom values.
917 //
918 {
919  int i;
920  int j;
921  int k;
922 
923  for ( j = 0; j < n; j++ )
924  {
925  for ( i = 0; i < m; i++ )
926  {
927  k = *seed / 127773;
928 
929  *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
930 
931  if ( *seed < 0 )
932  {
933  *seed = *seed + 2147483647;
934  }
935 //
936 // Although SEED can be represented exactly as a 32 bit integer,
937 // it generally cannot be represented exactly as a 32 bit real number!
938 //
939  r[i+j*m] = b + ( c - b ) * double( *seed ) * 4.656612875E-10;
940  }
941  }
942 
943  return;
944 }
945 //******************************************************************************
946 
947 int dtris2 ( int point_num, double point_xy[], int *tri_num,
948  int tri_vert[], int tri_nabe[] )
949 
950 //******************************************************************************
951 //
952 // Purpose:
953 //
954 // DTRIS2 constructs a Delaunay triangulation of 2D vertices.
955 //
956 // Discussion:
957 //
958 // The routine constructs the Delaunay triangulation of a set of 2D vertices
959 // using an incremental approach and diagonal edge swaps. Vertices are
960 // first sorted in lexicographically increasing (X,Y) order, and
961 // then are inserted one at a time from outside the convex hull.
962 //
963 // Modified:
964 //
965 // 15 January 2004
966 //
967 // Author:
968 //
969 // Barry Joe,
970 // Department of Computing Science,
971 // University of Alberta,
972 // Edmonton, Alberta, Canada T6G 2H1
973 //
974 // Reference:
975 //
976 // Barry Joe,
977 // GEOMPACK - a software package for the generation of meshes
978 // using geometric algorithms,
979 // Advances in Engineering Software,
980 // Volume 13, pages 325-331, 1991.
981 //
982 // Parameters:
983 //
984 // Input, int POINT_NUM, the number of vertices.
985 //
986 // Input/output, double POINT_XY[POINT_NUM*2], the coordinates of the vertices.
987 // On output, the vertices have been sorted into dictionary order.
988 //
989 // Output, int *TRI_NUM, the number of triangles in the triangulation;
990 // TRI_NUM is equal to 2*POINT_NUM - NB - 2, where NB is the number
991 // of boundary vertices.
992 //
993 // Output, int TRI_VERT[TRI_NUM*3], the nodes that make up each triangle.
994 // The elements are indices of POINT_XY. The vertices of the triangles are
995 // in counter clockwise order.
996 //
997 // Output, int TRI_NABE[TRI_NUM*3], the triangle neighbor list.
998 // Positive elements are indices of TIL; negative elements are used for links
999 // of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1)
1000 // where I, J = triangle, edge index; TRI_NABE[I,J] refers to
1001 // the neighbor along edge from vertex J to J+1 (mod 3).
1002 //
1003 // Output, int DTRIS2, is 0 for no error.
1004 {
1005  double cmax;
1006  int e;
1007  int error;
1008  int i;
1009  int *indx;
1010  int j;
1011  int k;
1012  int l;
1013  int ledg;
1014  int lr;
1015  int ltri;
1016  int m;
1017  int m1;
1018  int m2;
1019  int n;
1020  int redg;
1021  int rtri;
1022  int *stack;
1023  int t;
1024  double tol;
1025  int top;
1026 
1027  stack = new int[point_num];
1028 
1029  tol = 100.0 * d_epsilon ( );
1030 //
1031 // Sort the vertices by increasing (x,y).
1032 //
1033  indx = d2vec_sort_heap_index_a ( point_num, point_xy );
1034 
1035  d2vec_permute ( point_num, point_xy, indx );
1036 //
1037 // Make sure that the data points are "reasonably" distinct.
1038 //
1039  m1 = 1;
1040 
1041  for ( i = 2; i <= point_num; i++ )
1042  {
1043  m = m1;
1044  m1 = i;
1045 
1046  k = -1;
1047 
1048  for ( j = 0; j <= 1; j++ )
1049  {
1050  cmax = d_max ( fabs ( point_xy[2*(m-1)+j] ),
1051  fabs ( point_xy[2*(m1-1)+j] ) );
1052 
1053  if ( tol * ( cmax + 1.0 )
1054  < fabs ( point_xy[2*(m-1)+j] - point_xy[2*(m1-1)+j] ) )
1055  {
1056  k = j;
1057  break;
1058  }
1059 
1060  }
1061 
1062  if ( k == -1 )
1063  {
1064  cout << "\n";
1065  cout << "DTRIS2 - Fatal error!\n";
1066  cout << " Fails for point number I = " << i << "\n";
1067  cout << " M = " << m << "\n";
1068  cout << " M1 = " << m1 << "\n";
1069  cout << " X,Y(M) = " << point_xy[2*(m-1)+0] << " "
1070  << point_xy[2*(m-1)+1] << "\n";
1071  cout << " X,Y(M1) = " << point_xy[2*(m1-1)+0] << " "
1072  << point_xy[2*(m1-1)+1] << "\n";
1073  delete [] stack;
1074  return 224;
1075  }
1076 
1077  }
1078 //
1079 // Starting from points M1 and M2, search for a third point M that
1080 // makes a "healthy" triangle (M1,M2,M)
1081 //
1082  m1 = 1;
1083  m2 = 2;
1084  j = 3;
1085 
1086  for ( ; ; )
1087  {
1088  if ( point_num < j )
1089  {
1090  cout << "\n";
1091  cout << "DTRIS2 - Fatal error!\n";
1092  delete [] stack;
1093  return 225;
1094  }
1095 
1096  m = j;
1097 
1098  lr = lrline ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1],
1099  point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1],
1100  point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0 );
1101 
1102  if ( lr != 0 )
1103  {
1104  break;
1105  }
1106 
1107  j = j + 1;
1108 
1109  }
1110 //
1111 // Set up the triangle information for (M1,M2,M), and for any other
1112 // triangles you created because points were collinear with M1, M2.
1113 //
1114  *tri_num = j - 2;
1115 
1116  if ( lr == -1 )
1117  {
1118  tri_vert[3*0+0] = m1;
1119  tri_vert[3*0+1] = m2;
1120  tri_vert[3*0+2] = m;
1121  tri_nabe[3*0+2] = -3;
1122 
1123  for ( i = 2; i <= *tri_num; i++ )
1124  {
1125  m1 = m2;
1126  m2 = i+1;
1127  tri_vert[3*(i-1)+0] = m1;
1128  tri_vert[3*(i-1)+1] = m2;
1129  tri_vert[3*(i-1)+2] = m;
1130  tri_nabe[3*(i-1)+0] = -3 * i;
1131  tri_nabe[3*(i-1)+1] = i;
1132  tri_nabe[3*(i-1)+2] = i - 1;
1133 
1134  }
1135 
1136  tri_nabe[3*(*tri_num-1)+0] = -3 * (*tri_num) - 1;
1137  tri_nabe[3*(*tri_num-1)+1] = -5;
1138  ledg = 2;
1139  ltri = *tri_num;
1140  }
1141  else
1142  {
1143  tri_vert[3*0+0] = m2;
1144  tri_vert[3*0+1] = m1;
1145  tri_vert[3*0+2] = m;
1146  tri_nabe[3*0+0] = -4;
1147 
1148  for ( i = 2; i <= *tri_num; i++ )
1149  {
1150  m1 = m2;
1151  m2 = i+1;
1152  tri_vert[3*(i-1)+0] = m2;
1153  tri_vert[3*(i-1)+1] = m1;
1154  tri_vert[3*(i-1)+2] = m;
1155  tri_nabe[3*(i-2)+2] = i;
1156  tri_nabe[3*(i-1)+0] = -3 * i - 3;
1157  tri_nabe[3*(i-1)+1] = i - 1;
1158  }
1159 
1160  tri_nabe[3*(*tri_num-1)+2] = -3 * (*tri_num);
1161  tri_nabe[3*0+1] = -3 * (*tri_num) - 2;
1162  ledg = 2;
1163  ltri = 1;
1164  }
1165 //
1166 // Insert the vertices one at a time from outside the convex hull,
1167 // determine visible boundary edges, and apply diagonal edge swaps until
1168 // Delaunay triangulation of vertices (so far) is obtained.
1169 //
1170  top = 0;
1171 
1172  for ( i = j+1; i <= point_num; i++ )
1173  {
1174  m = i;
1175  m1 = tri_vert[3*(ltri-1)+ledg-1];
1176 
1177  if ( ledg <= 2 )
1178  {
1179  m2 = tri_vert[3*(ltri-1)+ledg];
1180  }
1181  else
1182  {
1183  m2 = tri_vert[3*(ltri-1)+0];
1184  }
1185 
1186  lr = lrline ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1],
1187  point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1],
1188  point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0 );
1189 
1190  if ( 0 < lr )
1191  {
1192  rtri = ltri;
1193  redg = ledg;
1194  ltri = 0;
1195  }
1196  else
1197  {
1198  l = -tri_nabe[3*(ltri-1)+ledg-1];
1199  rtri = l / 3;
1200  redg = (l % 3) + 1;
1201  }
1202 
1203  vbedg ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1], point_num,
1204  point_xy, *tri_num, tri_vert, tri_nabe, &ltri, &ledg, &rtri, &redg );
1205 
1206  n = *tri_num + 1;
1207  l = -tri_nabe[3*(ltri-1)+ledg-1];
1208 
1209  for ( ; ; )
1210  {
1211  t = l / 3;
1212  e = ( l % 3 ) + 1;
1213  l = -tri_nabe[3*(t-1)+e-1];
1214  m2 = tri_vert[3*(t-1)+e-1];
1215 
1216  if ( e <= 2 )
1217  {
1218  m1 = tri_vert[3*(t-1)+e];
1219  }
1220  else
1221  {
1222  m1 = tri_vert[3*(t-1)+0];
1223  }
1224 
1225  *tri_num = *tri_num + 1;
1226  tri_nabe[3*(t-1)+e-1] = *tri_num;
1227  tri_vert[3*(*tri_num-1)+0] = m1;
1228  tri_vert[3*(*tri_num-1)+1] = m2;
1229  tri_vert[3*(*tri_num-1)+2] = m;
1230  tri_nabe[3*(*tri_num-1)+0] = t;
1231  tri_nabe[3*(*tri_num-1)+1] = *tri_num - 1;
1232  tri_nabe[3*(*tri_num-1)+2] = *tri_num + 1;
1233  top = top + 1;
1234 
1235  if ( point_num < top )
1236  {
1237  cout << "\n";
1238  cout << "DTRIS2 - Fatal error!\n";
1239  cout << " Stack overflow.\n";
1240  delete [] stack;
1241  return 8;
1242  }
1243 
1244  stack[top-1] = *tri_num;
1245 
1246  if ( t == rtri && e == redg )
1247  {
1248  break;
1249  }
1250 
1251  }
1252 
1253  tri_nabe[3*(ltri-1)+ledg-1] = -3 * n - 1;
1254  tri_nabe[3*(n-1)+1] = -3 * (*tri_num) - 2;
1255  tri_nabe[3*(*tri_num-1)+2] = -l;
1256  ltri = n;
1257  ledg = 2;
1258 
1259  error = swapec ( m, &top, &ltri, &ledg, point_num, point_xy, *tri_num,
1260  tri_vert, tri_nabe, stack );
1261 
1262  if ( error != 0 )
1263  {
1264  cout << "\n";
1265  cout << "DTRIS2 - Fatal error!\n";
1266  cout << " Error return from SWAPEC.\n";
1267  delete [] stack;
1268  return error;
1269  }
1270 
1271  }
1272 //
1273 // Now account for the sorting that we did.
1274 //
1275  for ( i = 0; i < 3; i++ )
1276  {
1277  for ( j = 0; j < *tri_num; j++ )
1278  {
1279  tri_vert[i+j*3] = indx [ tri_vert[i+j*3] - 1 ];
1280  }
1281  }
1282 
1283  perm_inv ( point_num, indx );
1284 
1285  d2vec_permute ( point_num, point_xy, indx );
1286 
1287  delete [] indx;
1288  delete [] stack;
1289 
1290  return 0;
1291 }
1292 //******************************************************************************
1293 
1294 bool dvec_eq ( int n, double a1[], double a2[] )
1295 
1296 //******************************************************************************
1297 //
1298 // Purpose:
1299 //
1300 // DVEC_EQ is true if every pair of entries in two vectors is equal.
1301 //
1302 // Modified:
1303 //
1304 // 28 August 2003
1305 //
1306 // Author:
1307 //
1308 // John Burkardt
1309 //
1310 // Parameters:
1311 //
1312 // Input, int N, the number of entries in the vectors.
1313 //
1314 // Input, double A1[N], A2[N], two vectors to compare.
1315 //
1316 // Output, bool DVEC_EQ.
1317 // DVEC_EQ is TRUE if every pair of elements A1(I) and A2(I) are equal,
1318 // and FALSE otherwise.
1319 //
1320 {
1321  int i;
1322 
1323  for ( i = 0; i < n; i++ )
1324  {
1325  if ( a1[i] != a2[i] )
1326  {
1327  return false;
1328  }
1329  }
1330  return true;
1331 
1332 }
1333 //******************************************************************************
1334 
1335 bool dvec_gt ( int n, double a1[], double a2[] )
1336 
1337 //******************************************************************************
1338 //
1339 // Purpose:
1340 //
1341 // DVEC_GT == ( A1 > A2 ) for real vectors.
1342 //
1343 // Discussion:
1344 //
1345 // The comparison is lexicographic.
1346 //
1347 // A1 > A2 <=> A1(1) > A2(1) or
1348 // ( A1(1) == A2(1) and A1(2) > A2(2) ) or
1349 // ...
1350 // ( A1(1:N-1) == A2(1:N-1) and A1(N) > A2(N)
1351 //
1352 // Modified:
1353 //
1354 // 28 August 2003
1355 //
1356 // Author:
1357 //
1358 // John Burkardt
1359 //
1360 // Parameters:
1361 //
1362 // Input, int N, the dimension of the vectors.
1363 //
1364 // Input, double A1[N], A2[N], the vectors to be compared.
1365 //
1366 // Output, bool DVEC_GT, is TRUE if and only if A1 > A2.
1367 //
1368 {
1369  int i;
1370 
1371  for ( i = 0; i < n; i++ )
1372  {
1373 
1374  if ( a2[i] < a1[i] )
1375  {
1376  return true;
1377  }
1378  else if ( a1[i] < a2[i] )
1379  {
1380  return false;
1381  }
1382 
1383  }
1384 
1385  return false;
1386 }
1387 //******************************************************************************
1388 
1389 bool dvec_lt ( int n, double a1[], double a2[] )
1390 
1391 //******************************************************************************
1392 //
1393 // Purpose:
1394 //
1395 // DVEC_LT == ( A1 < A2 ) for real vectors.
1396 //
1397 // Discussion:
1398 //
1399 // The comparison is lexicographic.
1400 //
1401 // A1 < A2 <=> A1(1) < A2(1) or
1402 // ( A1(1) == A2(1) and A1(2) < A2(2) ) or
1403 // ...
1404 // ( A1(1:N-1) == A2(1:N-1) and A1(N) < A2(N)
1405 //
1406 // Modified:
1407 //
1408 // 28 August 2003
1409 //
1410 // Author:
1411 //
1412 // John Burkardt
1413 //
1414 // Parameters:
1415 //
1416 // Input, int N, the dimension of the vectors.
1417 //
1418 // Input, double A1[N], A2[N], the vectors to be compared.
1419 //
1420 // Output, bool DVEC_LT, is TRUE if and only if A1 < A2.
1421 //
1422 {
1423  int i;
1424 
1425  for ( i = 0; i < n; i++ )
1426  {
1427  if ( a1[i] < a2[i] )
1428  {
1429  return true;
1430  }
1431  else if ( a2[i] < a1[i] )
1432  {
1433  return false;
1434  }
1435 
1436  }
1437 
1438  return false;
1439 }
1440 //********************************************************************
1441 
1442 void dvec_print ( int n, double a[], const char *title )
1443 
1444 //********************************************************************
1445 //
1446 // Purpose:
1447 //
1448 // DVEC_PRINT prints a double precision vector.
1449 //
1450 // Modified:
1451 //
1452 // 16 August 2004
1453 //
1454 // Author:
1455 //
1456 // John Burkardt
1457 //
1458 // Parameters:
1459 //
1460 // Input, int N, the number of components of the vector.
1461 //
1462 // Input, double A[N], the vector to be printed.
1463 //
1464 // Input, const char *TITLE, a title to be printed first.
1465 // TITLE may be blank.
1466 //
1467 {
1468  int i;
1469 
1470  if ( 0 < s_len_trim ( title ) )
1471  {
1472  cout << "\n";
1473  cout << title << "\n";
1474  }
1475 
1476  cout << "\n";
1477  for ( i = 0; i <= n-1; i++ )
1478  {
1479  cout << setw(6) << i + 1 << " "
1480  << setw(14) << a[i] << "\n";
1481  }
1482 
1483  return;
1484 }
1485 //******************************************************************************
1486 
1487 void dvec_swap ( int n, double a1[], double a2[] )
1488 
1489 //******************************************************************************
1490 //
1491 // Purpose:
1492 //
1493 // DVEC_SWAP swaps the entries of two real vectors.
1494 //
1495 // Modified:
1496 //
1497 // 28 August 2003
1498 //
1499 // Author:
1500 //
1501 // John Burkardt
1502 //
1503 // Parameters:
1504 //
1505 // Input, int N, the number of entries in the arrays.
1506 //
1507 // Input/output, double A1[N], A2[N], the vectors to swap.
1508 //
1509 {
1510  int i;
1511  double temp;
1512 
1513  for ( i = 0; i < n; i++ )
1514  {
1515  temp = a1[i];
1516  a1[i] = a2[i];
1517  a2[i] = temp;
1518  }
1519 
1520  return;
1521 }
1522 //****************************************************************************
1523 
1524 int i_max ( int i1, int i2 )
1525 
1526 //****************************************************************************
1527 //
1528 // Purpose:
1529 //
1530 // I_MAX returns the maximum of two integers.
1531 //
1532 // Modified:
1533 //
1534 // 13 October 1998
1535 //
1536 // Author:
1537 //
1538 // John Burkardt
1539 //
1540 // Parameters:
1541 //
1542 // Input, int I1, I2, are two integers to be compared.
1543 //
1544 // Output, int I_MAX, the larger of I1 and I2.
1545 //
1546 {
1547  if ( i2 < i1 )
1548  {
1549  return i1;
1550  }
1551  else
1552  {
1553  return i2;
1554  }
1555 
1556 }
1557 //****************************************************************************
1558 
1559 int i_min ( int i1, int i2 )
1560 
1561 //****************************************************************************
1562 //
1563 // Purpose:
1564 //
1565 // I_MIN returns the smaller of two integers.
1566 //
1567 // Modified:
1568 //
1569 // 13 October 1998
1570 //
1571 // Author:
1572 //
1573 // John Burkardt
1574 //
1575 // Parameters:
1576 //
1577 // Input, int I1, I2, two integers to be compared.
1578 //
1579 // Output, int I_MIN, the smaller of I1 and I2.
1580 //
1581 {
1582  if ( i1 < i2 )
1583  {
1584  return i1;
1585  }
1586  else
1587  {
1588  return i2;
1589  }
1590 
1591 }
1592 //*********************************************************************
1593 
1594 int i_modp ( int i, int j )
1595 
1596 //*********************************************************************
1597 //
1598 // Purpose:
1599 //
1600 // I_MODP returns the nonnegative remainder of integer division.
1601 //
1602 // Formula:
1603 //
1604 // If
1605 // NREM = I_MODP ( I, J )
1606 // NMULT = ( I - NREM ) / J
1607 // then
1608 // I = J * NMULT + NREM
1609 // where NREM is always nonnegative.
1610 //
1611 // Comments:
1612 //
1613 // The MOD function computes a result with the same sign as the
1614 // quantity being divided. Thus, suppose you had an angle A,
1615 // and you wanted to ensure that it was between 0 and 360.
1616 // Then mod(A,360) would do, if A was positive, but if A
1617 // was negative, your result would be between -360 and 0.
1618 //
1619 // On the other hand, I_MODP(A,360) is between 0 and 360, always.
1620 //
1621 // Examples:
1622 //
1623 // I J MOD I_MODP I_MODP Factorization
1624 //
1625 // 107 50 7 7 107 = 2 * 50 + 7
1626 // 107 -50 7 7 107 = -2 * -50 + 7
1627 // -107 50 -7 43 -107 = -3 * 50 + 43
1628 // -107 -50 -7 43 -107 = 3 * -50 + 43
1629 //
1630 // Modified:
1631 //
1632 // 26 May 1999
1633 //
1634 // Author:
1635 //
1636 // John Burkardt
1637 //
1638 // Parameters:
1639 //
1640 // Input, int I, the number to be divided.
1641 //
1642 // Input, int J, the number that divides I.
1643 //
1644 // Output, int I_MODP, the nonnegative remainder when I is
1645 // divided by J.
1646 //
1647 {
1648  int value;
1649 
1650  if ( j == 0 )
1651  {
1652  cout << "\n";
1653  cout << "I_MODP - Fatal error!\n";
1654  cout << " I_MODP ( I, J ) called with J = " << j << "\n";
1655  exit ( 1 );
1656  }
1657 
1658  value = i % j;
1659 
1660  if ( value < 0 )
1661  {
1662  value = value + abs ( j );
1663  }
1664 
1665  return value;
1666 }
1667 //********************************************************************
1668 
1669 int i_sign ( int i )
1670 
1671 //********************************************************************
1672 //
1673 // Purpose:
1674 //
1675 // I_SIGN returns the sign of an integer.
1676 //
1677 // Discussion:
1678 //
1679 // The sign of 0 and all positive integers is taken to be +1.
1680 // The sign of all negative integers is -1.
1681 //
1682 // Modified:
1683 //
1684 // 06 May 2003
1685 //
1686 // Author:
1687 //
1688 // John Burkardt
1689 //
1690 // Parameters:
1691 //
1692 // Input, int I, the integer whose sign is desired.
1693 //
1694 // Output, int I_SIGN, the sign of I.
1695 {
1696  if ( i < 0 )
1697  {
1698  return (-1);
1699  }
1700  else
1701  {
1702  return 1;
1703  }
1704 
1705 }
1706 //*******************************************************************************
1707 
1708 int i_wrap ( int ival, int ilo, int ihi )
1709 
1710 //*******************************************************************************
1711 //
1712 // Purpose:
1713 //
1714 // I_WRAP forces an integer to lie between given limits by wrapping.
1715 //
1716 // Example:
1717 //
1718 // ILO = 4, IHI = 8
1719 //
1720 // I I_WRAP
1721 //
1722 // -2 8
1723 // -1 4
1724 // 0 5
1725 // 1 6
1726 // 2 7
1727 // 3 8
1728 // 4 4
1729 // 5 5
1730 // 6 6
1731 // 7 7
1732 // 8 8
1733 // 9 4
1734 // 10 5
1735 // 11 6
1736 // 12 7
1737 // 13 8
1738 // 14 4
1739 //
1740 // Modified:
1741 //
1742 // 19 August 2003
1743 //
1744 // Author:
1745 //
1746 // John Burkardt
1747 //
1748 // Parameters:
1749 //
1750 // Input, int IVAL, an integer value.
1751 //
1752 // Input, int ILO, IHI, the desired bounds for the integer value.
1753 //
1754 // Output, int I_WRAP, a "wrapped" version of IVAL.
1755 //
1756 {
1757  int jhi;
1758  int jlo;
1759  int value;
1760  int wide;
1761 
1762  jlo = i_min ( ilo, ihi );
1763  jhi = i_max ( ilo, ihi );
1764 
1765  wide = jhi + 1 - jlo;
1766 
1767  if ( wide == 1 )
1768  {
1769  value = jlo;
1770  }
1771  else
1772  {
1773  value = jlo + i_modp ( ival - jlo, wide );
1774  }
1775 
1776  return value;
1777 }
1778 //******************************************************************************
1779 
1780 void imat_transpose_print ( int m, int n, int a[], const char *title )
1781 
1782 //******************************************************************************
1783 //
1784 // Purpose:
1785 //
1786 // IMAT_TRANSPOSE_PRINT prints an integer matrix, transposed.
1787 //
1788 // Modified:
1789 //
1790 // 31 January 2005
1791 //
1792 // Author:
1793 //
1794 // John Burkardt
1795 //
1796 // Parameters:
1797 //
1798 // Input, int M, the number of rows in A.
1799 //
1800 // Input, int N, the number of columns in A.
1801 //
1802 // Input, int A[M*N], the M by N matrix.
1803 //
1804 // Input, const char *TITLE, a title to be printed.
1805 //
1806 {
1807  imat_transpose_print_some ( m, n, a, 1, 1, m, n, title );
1808 
1809  return;
1810 }
1811 //******************************************************************************
1812 
1813 void imat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo,
1814  int ihi, int jhi, const char *title )
1815 
1816 //******************************************************************************
1817 //
1818 // Purpose:
1819 //
1820 // IMAT_TRANSPOSE_PRINT_SOME prints some of an integer matrix, transposed.
1821 //
1822 // Modified:
1823 //
1824 // 09 February 2005
1825 //
1826 // Author:
1827 //
1828 // John Burkardt
1829 //
1830 // Parameters:
1831 //
1832 // Input, int M, the number of rows of the matrix.
1833 // M must be positive.
1834 //
1835 // Input, int N, the number of columns of the matrix.
1836 // N must be positive.
1837 //
1838 // Input, int A[M*N], the matrix.
1839 //
1840 // Input, int ILO, JLO, IHI, JHI, designate the first row and
1841 // column, and the last row and column to be printed.
1842 //
1843 // Input, const char *TITLE, a title for the matrix.
1844 {
1845 # define INCX 10
1846 
1847  int i;
1848  int i2hi;
1849  int i2lo;
1850  int j;
1851  int j2hi;
1852  int j2lo;
1853 
1854  if ( 0 < s_len_trim ( title ) )
1855  {
1856  cout << "\n";
1857  cout << title << "\n";
1858  }
1859 //
1860 // Print the columns of the matrix, in strips of INCX.
1861 //
1862  for ( i2lo = ilo; i2lo <= ihi; i2lo = i2lo + INCX )
1863  {
1864  i2hi = i2lo + INCX - 1;
1865  i2hi = i_min ( i2hi, m );
1866  i2hi = i_min ( i2hi, ihi );
1867 
1868  cout << "\n";
1869 //
1870 // For each row I in the current range...
1871 //
1872 // Write the header.
1873 //
1874  cout << " Row: ";
1875  for ( i = i2lo; i <= i2hi; i++ )
1876  {
1877  cout << setw(7) << i << " ";
1878  }
1879  cout << "\n";
1880  cout << " Col\n";
1881  cout << "\n";
1882 //
1883 // Determine the range of the rows in this strip.
1884 //
1885  j2lo = i_max ( jlo, 1 );
1886  j2hi = i_min ( jhi, n );
1887 
1888  for ( j = j2lo; j <= j2hi; j++ )
1889  {
1890 //
1891 // Print out (up to INCX) entries in column J, that lie in the current strip.
1892 //
1893  cout << setw(5) << j << " ";
1894  for ( i = i2lo; i <= i2hi; i++ )
1895  {
1896  cout << setw(6) << a[i-1+(j-1)*m] << " ";
1897  }
1898  cout << "\n";
1899  }
1900 
1901  }
1902 
1903  cout << "\n";
1904 
1905  return;
1906 # undef INCX
1907 }
1908 //********************************************************************
1909 
1910 void ivec_heap_d ( int n, int a[] )
1911 
1912 /*********************************************************************
1913 //
1914 // Purpose:
1915 //
1916 // IVEC_HEAP_D reorders an array of integers into a descending heap.
1917 //
1918 // Definition:
1919 //
1920 // A heap is an array A with the property that, for every index J,
1921 // A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices
1922 // 2*J+1 and 2*J+2 are legal).
1923 //
1924 // Diagram:
1925 //
1926 // A(0)
1927 // / \
1928 // A(1) A(2)
1929 // / \ / \
1930 // A(3) A(4) A(5) A(6)
1931 // / \ / \
1932 // A(7) A(8) A(9) A(10)
1933 //
1934 // Reference:
1935 //
1936 // Albert Nijenhuis, Herbert Wilf,
1937 // Combinatorial Algorithms,
1938 // Academic Press, 1978, second edition,
1939 // ISBN 0-12-519260-6.
1940 //
1941 // Modified:
1942 //
1943 // 30 April 1999
1944 //
1945 // Author:
1946 //
1947 // John Burkardt
1948 //
1949 // Parameters:
1950 //
1951 // Input, int N, the size of the input array.
1952 //
1953 // Input/output, int A[N].
1954 // On input, an unsorted array.
1955 // On output, the array has been reordered into a heap.
1956 */
1957 {
1958  int i;
1959  int ifree;
1960  int key;
1961  int m;
1962 //
1963 // Only nodes (N/2)-1 down to 0 can be "parent" nodes.
1964 //
1965  for ( i = (n/2)-1; 0 <= i; i-- )
1966  {
1967 //
1968 // Copy the value out of the parent node.
1969 // Position IFREE is now "open".
1970 //
1971  key = a[i];
1972  ifree = i;
1973 
1974  for ( ;; )
1975  {
1976 //
1977 // Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position
1978 // IFREE. (One or both may not exist because they equal or exceed N.)
1979 //
1980  m = 2 * ifree + 1;
1981 //
1982 // Does the first position exist?
1983 //
1984  if ( n <= m )
1985  {
1986  break;
1987  }
1988  else
1989  {
1990 //
1991 // Does the second position exist?
1992 //
1993  if ( m + 1 < n )
1994  {
1995 //
1996 // If both positions exist, take the larger of the two values,
1997 // and update M if necessary.
1998 //
1999  if ( a[m] < a[m+1] )
2000  {
2001  m = m + 1;
2002  }
2003  }
2004 //
2005 // If the large descendant is larger than KEY, move it up,
2006 // and update IFREE, the location of the free position, and
2007 // consider the descendants of THIS position.
2008 //
2009  if ( key < a[m] )
2010  {
2011  a[ifree] = a[m];
2012  ifree = m;
2013  }
2014  else
2015  {
2016  break;
2017  }
2018 
2019  }
2020 
2021  }
2022 //
2023 // When you have stopped shifting items up, return the item you
2024 // pulled out back to the heap.
2025 //
2026  a[ifree] = key;
2027 
2028  }
2029 
2030  return;
2031 }
2032 //******************************************************************************
2033 
2034 int *ivec_indicator ( int n )
2035 
2036 //******************************************************************************
2037 //
2038 // Purpose:
2039 //
2040 // IVEC_INDICATOR sets an integer vector to the indicator vector.
2041 //
2042 // Modified:
2043 //
2044 // 13 January 2004
2045 //
2046 // Author:
2047 //
2048 // John Burkardt
2049 //
2050 // Parameters:
2051 //
2052 // Input, int N, the number of elements of A.
2053 //
2054 // Output, int IVEC_INDICATOR(N), the initialized array.
2055 //
2056 {
2057  int *a;
2058  int i;
2059 
2060  a = new int[n];
2061 
2062  for ( i = 0; i < n; i++ )
2063  {
2064  a[i] = i + 1;
2065  }
2066 
2067  return a;
2068 }
2069 //********************************************************************
2070 
2071 void ivec_sort_heap_a ( int n, int a[] )
2072 
2073 //********************************************************************
2074 //
2075 // Purpose:
2076 //
2077 // IVEC_SORT_HEAP_A ascending sorts an array of integers using heap sort.
2078 //
2079 // Reference:
2080 //
2081 // Albert Nijenhuis, Herbert Wilf,
2082 // Combinatorial Algorithms,
2083 // Academic Press, 1978, second edition,
2084 // ISBN 0-12-519260-6.
2085 //
2086 // Modified:
2087 //
2088 // 30 April 1999
2089 //
2090 // Author:
2091 //
2092 // John Burkardt
2093 //
2094 // Parameters:
2095 //
2096 // Input, int N, the number of entries in the array.
2097 //
2098 // Input/output, int A[N].
2099 // On input, the array to be sorted;
2100 // On output, the array has been sorted.
2101 //
2102 {
2103  int n1;
2104  int temp;
2105 
2106  if ( n <= 1 )
2107  {
2108  return;
2109  }
2110 //
2111 // 1: Put A into descending heap form.
2112 //
2113  ivec_heap_d ( n, a );
2114 //
2115 // 2: Sort A.
2116 //
2117 // The largest object in the heap is in A[0].
2118 // Move it to position A[N-1].
2119 //
2120  temp = a[0];
2121  a[0] = a[n-1];
2122  a[n-1] = temp;
2123 //
2124 // Consider the diminished heap of size N1.
2125 //
2126  for ( n1 = n-1; 2 <= n1; n1-- )
2127  {
2128 //
2129 // Restore the heap structure of the initial N1 entries of A.
2130 //
2131  ivec_heap_d ( n1, a );
2132 //
2133 // Take the largest object from A[0] and move it to A[N1-1].
2134 //
2135  temp = a[0];
2136  a[0] = a[n1-1];
2137  a[n1-1] = temp;
2138 
2139  }
2140 
2141  return;
2142 }
2143 //******************************************************************************
2144 
2145 void ivec_sorted_unique ( int n, int a[], int *nuniq )
2146 
2147 //******************************************************************************
2148 //
2149 // Purpose:
2150 //
2151 // IVEC_SORTED_UNIQUE finds unique elements in a sorted integer array.
2152 //
2153 // Modified:
2154 //
2155 // 02 September 2003
2156 //
2157 // Author:
2158 //
2159 // John Burkardt
2160 //
2161 // Parameters:
2162 //
2163 // Input, int N, the number of elements in A.
2164 //
2165 // Input/output, int A[N]. On input, the sorted
2166 // integer array. On output, the unique elements in A.
2167 //
2168 // Output, int *NUNIQ, the number of unique elements in A.
2169 //
2170 {
2171  int i;
2172 
2173  *nuniq = 0;
2174 
2175  if ( n <= 0 )
2176  {
2177  return;
2178  }
2179 
2180  *nuniq = 1;
2181 
2182  for ( i = 1; i < n; i++ )
2183  {
2184  if ( a[i] != a[*nuniq] )
2185  {
2186  *nuniq = *nuniq + 1;
2187  a[*nuniq] = a[i];
2188  }
2189 
2190  }
2191 
2192  return;
2193 }
2194 //******************************************************************************
2195 
2196 int lrline ( double xu, double yu, double xv1, double yv1, double xv2,
2197  double yv2, double dv )
2198 
2199 //******************************************************************************
2200 //
2201 // Purpose:
2202 //
2203 // LRLINE determines where a point lies in relation to a directed line.
2204 //
2205 // Discussion:
2206 //
2207 // LRLINE determines whether a point is to the left of, right of,
2208 // or on a directed line parallel to a line through given points.
2209 //
2210 // Modified:
2211 //
2212 // 28 August 2003
2213 //
2214 // Author:
2215 //
2216 // Barry Joe,
2217 // Department of Computing Science,
2218 // University of Alberta,
2219 // Edmonton, Alberta, Canada T6G 2H1
2220 //
2221 // Reference:
2222 //
2223 // Barry Joe,
2224 // GEOMPACK - a software package for the generation of meshes
2225 // using geometric algorithms,
2226 // Advances in Engineering Software,
2227 // Volume 13, pages 325-331, 1991.
2228 //
2229 // Parameters:
2230 //
2231 // Input, double XU, YU, XV1, YV1, XV2, YV2, are vertex coordinates; the
2232 // directed line is parallel to and at signed distance DV to the left of
2233 // the directed line from (XV1,YV1) to (XV2,YV2); (XU,YU) is the vertex for
2234 // which the position relative to the directed line is to be determined.
2235 //
2236 // Input, double DV, the signed distance, positive for left.
2237 //
2238 // Output, int LRLINE, is +1, 0, or -1 depending on whether (XU,YU) is
2239 // to the right of, on, or left of the directed line. LRLINE is 0 if
2240 // the line degenerates to a point.
2241 //
2242 {
2243  double dx;
2244  double dxu;
2245  double dy;
2246  double dyu;
2247  double t;
2248  double tol = 0.0000001;
2249  double tolabs;
2250  int value = 0;
2251 
2252  dx = xv2 - xv1;
2253  dy = yv2 - yv1;
2254  dxu = xu - xv1;
2255  dyu = yu - yv1;
2256 
2257  tolabs = tol * d_max ( fabs ( dx ),
2258  d_max ( fabs ( dy ),
2259  d_max ( fabs ( dxu ),
2260  d_max ( fabs ( dyu ), fabs ( dv ) ) ) ) );
2261 
2262  t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy );
2263 
2264  if ( tolabs < t )
2265  {
2266  value = 1;
2267  }
2268  else if ( -tolabs <= t )
2269  {
2270  value = 0;
2271  }
2272  else if ( t < -tolabs )
2273  {
2274  value = -1;
2275  }
2276 
2277  return value;
2278 }
2279 //******************************************************************************
2280 
2281 bool perm_check ( int n, int p[] )
2282 
2283 //******************************************************************************
2284 //
2285 // Purpose:
2286 //
2287 // PERM_CHECK checks that a vector represents a permutation.
2288 //
2289 // Discussion:
2290 //
2291 // The routine verifies that each of the integers from 1
2292 // to N occurs among the N entries of the permutation.
2293 //
2294 // Modified:
2295 //
2296 // 13 January 2004
2297 //
2298 // Author:
2299 //
2300 // John Burkardt
2301 //
2302 // Parameters:
2303 //
2304 // Input, int N, the number of entries.
2305 //
2306 // Input, int P[N], the array to check.
2307 //
2308 // Output, bool PERM_CHECK, is TRUE if the permutation is OK.
2309 //
2310 {
2311  bool found;
2312  int i;
2313  int seek;
2314 
2315  for ( seek = 1; seek <= n; seek++ )
2316  {
2317  found = false;
2318 
2319  for ( i = 0; i < n; i++ )
2320  {
2321  if ( p[i] == seek )
2322  {
2323  found = true;
2324  break;
2325  }
2326  }
2327 
2328  if ( !found )
2329  {
2330  return false;
2331  }
2332 
2333  }
2334 
2335  return true;
2336 }
2337 //******************************************************************************
2338 
2339 void perm_inv ( int n, int p[] )
2340 
2341 //******************************************************************************
2342 //
2343 // Purpose:
2344 //
2345 // PERM_INV inverts a permutation "in place".
2346 //
2347 // Modified:
2348 //
2349 // 13 January 2004
2350 //
2351 // Parameters:
2352 //
2353 // Input, int N, the number of objects being permuted.
2354 //
2355 // Input/output, int P[N], the permutation, in standard index form.
2356 // On output, P describes the inverse permutation
2357 //
2358 {
2359  int i;
2360  int i0;
2361  int i1;
2362  int i2;
2363  int is;
2364 
2365  if ( n <= 0 )
2366  {
2367  cout << "\n";
2368  cout << "PERM_INV - Fatal error!\n";
2369  cout << " Input value of N = " << n << "\n";
2370  exit ( 1 );
2371  }
2372 
2373  if ( !perm_check ( n, p ) )
2374  {
2375  cout << "\n";
2376  cout << "PERM_INV - Fatal error!\n";
2377  cout << " The input array does not represent\n";
2378  cout << " a proper permutation.\n";
2379  exit ( 1 );
2380  }
2381 
2382  is = 1;
2383 
2384  for ( i = 1; i <= n; i++ )
2385  {
2386  i1 = p[i-1];
2387 
2388  while ( i < i1 )
2389  {
2390  i2 = p[i1-1];
2391  p[i1-1] = -i2;
2392  i1 = i2;
2393  }
2394 
2395  is = - i_sign ( p[i-1] );
2396  p[i-1] = i_sign ( is ) * abs ( p[i-1] );
2397  }
2398 
2399  for ( i = 1; i <= n; i++ )
2400  {
2401  i1 = -p[i-1];
2402 
2403  if ( 0 <= i1 )
2404  {
2405  i0 = i;
2406 
2407  for ( ; ; )
2408  {
2409  i2 = p[i1-1];
2410  p[i1-1] = i0;
2411 
2412  if ( i2 < 0 )
2413  {
2414  break;
2415  }
2416  i0 = i1;
2417  i1 = i2;
2418  }
2419  }
2420  }
2421 
2422  return;
2423 }
2424 //********************************************************************
2425 
2426 int *points_delaunay_naive_2d ( int n, double p[], int *ntri )
2427 
2428 //********************************************************************
2429 //
2430 // Purpose:
2431 //
2432 // POINTS_DELAUNAY_NAIVE_2D computes the Delaunay triangulation in 2D.
2433 //
2434 // Discussion:
2435 //
2436 // A naive and inefficient (but extremely simple) method is used.
2437 //
2438 // This routine is only suitable as a demonstration code for small
2439 // problems. Its running time is of order N^4. Much faster algorithms
2440 // are available.
2441 //
2442 // Given a set of nodes in the plane, a triangulation is a set of
2443 // triples of distinct nodes, forming triangles, so that every
2444 // point with the convex hull of the set of nodes is either one
2445 // of the nodes, or lies on an edge of one or more triangles,
2446 // or lies within exactly one triangle.
2447 //
2448 // Modified:
2449 //
2450 // 05 February 2005
2451 //
2452 // Author:
2453 //
2454 // John Burkardt
2455 //
2456 // Reference:
2457 //
2458 // Joseph O'Rourke,
2459 // Computational Geometry,
2460 // Cambridge University Press,
2461 // Second Edition, 1998, page 187.
2462 //
2463 // Parameters:
2464 //
2465 // Input, int N, the number of nodes. N must be at least 3.
2466 //
2467 // Input, double P[2*N], the coordinates of the nodes.
2468 //
2469 // Output, int *NTRI, the number of triangles.
2470 //
2471 // Output, int POINTS_DELAUNAY_NAIVE_2D[3*NTRI], the indices of the
2472 // nodes making each triangle.
2473 //
2474 {
2475  int count;
2476  int flag;
2477  int i;
2478  int j;
2479  int k;
2480  int m;
2481  int pass;
2482  int *tri = NULL;
2483  double xn;
2484  double yn;
2485  double zn;
2486  double *z;
2487 
2488  count = 0;
2489 
2490  z = new double [ n ];
2491 
2492  for ( i = 0; i < n; i++ )
2493  {
2494  z[i] = p[0+i*2] * p[0+i*2] + p[1+i*2] * p[1+i*2];
2495  }
2496 //
2497 // First pass counts triangles,
2498 // Second pass allocates triangles and sets them.
2499 //
2500  for ( pass = 1; pass <= 2; pass++ )
2501  {
2502  if ( pass == 2 )
2503  {
2504  tri = new int[3*count];
2505  }
2506  count = 0;
2507 //
2508 // For each triple (I,J,K):
2509 //
2510  for ( i = 0; i < n - 2; i++ )
2511  {
2512  for ( j = i+1; j < n; j++ )
2513  {
2514  for ( k = i+1; k < n; k++ )
2515  {
2516  if ( j != k )
2517  {
2518  xn = ( p[1+j*2] - p[1+i*2] ) * ( z[k] - z[i] )
2519  - ( p[1+k*2] - p[1+i*2] ) * ( z[j] - z[i] );
2520  yn = ( p[0+k*2] - p[0+i*2] ) * ( z[j] - z[i] )
2521  - ( p[0+j*2] - p[0+i*2] ) * ( z[k] - z[i] );
2522  zn = ( p[0+j*2] - p[0+i*2] ) * ( p[1+k*2] - p[1+i*2] )
2523  - ( p[0+k*2] - p[0+i*2] ) * ( p[1+j*2] - p[1+i*2] );
2524 
2525  flag = ( zn < 0 );
2526 
2527  if ( flag )
2528  {
2529  for ( m = 0; m < n; m++ )
2530  {
2531  flag = flag && ( ( p[0+m*2] - p[0+i*2] ) * xn
2532  + ( p[1+m*2] - p[1+i*2] ) * yn
2533  + ( z[m] - z[i] ) * zn <= 0 );
2534  }
2535  }
2536 
2537  if ( flag )
2538  {
2539  if ( pass == 2 )
2540  {
2541  tri[0+count*3] = i;
2542  tri[1+count*3] = j;
2543  tri[2+count*3] = k;
2544  }
2545  count = count + 1;
2546  }
2547 
2548  }
2549  }
2550  }
2551  }
2552  }
2553 
2554  *ntri = count;
2555  delete [] z;
2556 
2557  return tri;
2558 }
2559 //******************************************************************************
2560 
2561 int s_len_trim ( const char *s )
2562 
2563 //******************************************************************************
2564 //
2565 // Purpose:
2566 //
2567 // S_LEN_TRIM returns the length of a string to the last nonblank.
2568 //
2569 // Modified:
2570 //
2571 // 26 April 2003
2572 //
2573 // Author:
2574 //
2575 // John Burkardt
2576 //
2577 // Parameters:
2578 //
2579 // Input, const char *S, a pointer to a string.
2580 //
2581 // Output, int S_LEN_TRIM, the length of the string to the last nonblank.
2582 // If S_LEN_TRIM is 0, then the string is entirely blank.
2583 //
2584 {
2585  int n;
2586  char* t;
2587 
2588  n = strlen ( s );
2589  t = const_cast<char*>(s) + n - 1;
2590 
2591  while ( 0 < n )
2592  {
2593  if ( *t != ' ' )
2594  {
2595  return n;
2596  }
2597  t--;
2598  n--;
2599  }
2600 
2601  return n;
2602 }
2603 //******************************************************************************
2604 
2605 int swapec ( int i, int *top, int *btri, int *bedg, int point_num,
2606  double point_xy[], int tri_num, int tri_vert[], int tri_nabe[],
2607  int stack[] )
2608 
2609 //******************************************************************************
2610 //
2611 // Purpose:
2612 //
2613 // SWAPEC swaps diagonal edges until all triangles are Delaunay.
2614 //
2615 // Discussion:
2616 //
2617 // The routine swaps diagonal edges in a 2D triangulation, based on
2618 // the empty circumcircle criterion, until all triangles are Delaunay,
2619 // given that I is the index of the new vertex added to the triangulation.
2620 //
2621 // Modified:
2622 //
2623 // 03 September 2003
2624 //
2625 // Author:
2626 //
2627 // Barry Joe,
2628 // Department of Computing Science,
2629 // University of Alberta,
2630 // Edmonton, Alberta, Canada T6G 2H1
2631 //
2632 // Reference:
2633 //
2634 // Barry Joe,
2635 // GEOMPACK - a software package for the generation of meshes
2636 // using geometric algorithms,
2637 // Advances in Engineering Software,
2638 // Volume 13, pages 325-331, 1991.
2639 //
2640 // Parameters:
2641 //
2642 // Input, int I, the index of the new vertex.
2643 //
2644 // Input/output, int *TOP, the index of the top of the stack.
2645 // On output, TOP is zero.
2646 //
2647 // Input/output, int *BTRI, *BEDG; on input, if positive, are the
2648 // triangle and edge indices of a boundary edge whose updated indices
2649 // must be recorded. On output, these may be updated because of swaps.
2650 //
2651 // Input, int POINT_NUM, the number of points.
2652 //
2653 // Input, double POINT_XY[POINT_NUM*2], the coordinates of the points.
2654 //
2655 // Input, int TRI_NUM, the number of triangles.
2656 //
2657 // Input/output, int TRI_VERT[TRI_NUM*3], the triangle incidence list.
2658 // May be updated on output because of swaps.
2659 //
2660 // Input/output, int TRI_NABE[TRI_NUM*3], the triangle neighbor list;
2661 // negative values are used for links of the counter-clockwise linked
2662 // list of boundary edges; May be updated on output because of swaps.
2663 //
2664 // LINK = -(3*I + J-1) where I, J = triangle, edge index.
2665 //
2666 // Workspace, int STACK[MAXST]; on input, entries 1 through TOP
2667 // contain the indices of initial triangles (involving vertex I)
2668 // put in stack; the edges opposite I should be in interior; entries
2669 // TOP+1 through MAXST are used as a stack.
2670 //
2671 // Output, int SWAPEC, is set to 8 for abnormal return.
2672 //
2673 {
2674  int a;
2675  int b;
2676  int c;
2677  int e;
2678  int ee;
2679  int em1;
2680  int ep1;
2681  int f;
2682  int fm1;
2683  int fp1;
2684  int l;
2685  int r;
2686  int s;
2687  int swap;
2688  int t;
2689  int tt;
2690  int u;
2691  double x;
2692  double y;
2693 //
2694 // Determine whether triangles in stack are Delaunay, and swap
2695 // diagonal edge of convex quadrilateral if not.
2696 //
2697  x = point_xy[2*(i-1)+0];
2698  y = point_xy[2*(i-1)+1];
2699 
2700  for ( ; ; )
2701  {
2702  if ( *top <= 0 )
2703  {
2704  break;
2705  }
2706 
2707  t = stack[(*top)-1];
2708  *top = *top - 1;
2709 
2710  if ( tri_vert[3*(t-1)+0] == i )
2711  {
2712  e = 2;
2713  b = tri_vert[3*(t-1)+2];
2714  }
2715  else if ( tri_vert[3*(t-1)+1] == i )
2716  {
2717  e = 3;
2718  b = tri_vert[3*(t-1)+0];
2719  }
2720  else
2721  {
2722  e = 1;
2723  b = tri_vert[3*(t-1)+1];
2724  }
2725 
2726  a = tri_vert[3*(t-1)+e-1];
2727  u = tri_nabe[3*(t-1)+e-1];
2728 
2729  if ( tri_nabe[3*(u-1)+0] == t )
2730  {
2731  f = 1;
2732  c = tri_vert[3*(u-1)+2];
2733  }
2734  else if ( tri_nabe[3*(u-1)+1] == t )
2735  {
2736  f = 2;
2737  c = tri_vert[3*(u-1)+0];
2738  }
2739  else
2740  {
2741  f = 3;
2742  c = tri_vert[3*(u-1)+1];
2743  }
2744 
2745  swap = diaedg ( x, y,
2746  point_xy[2*(a-1)+0], point_xy[2*(a-1)+1],
2747  point_xy[2*(c-1)+0], point_xy[2*(c-1)+1],
2748  point_xy[2*(b-1)+0], point_xy[2*(b-1)+1] );
2749 
2750  if ( swap == 1 )
2751  {
2752  em1 = i_wrap ( e - 1, 1, 3 );
2753  ep1 = i_wrap ( e + 1, 1, 3 );
2754  fm1 = i_wrap ( f - 1, 1, 3 );
2755  fp1 = i_wrap ( f + 1, 1, 3 );
2756 
2757  tri_vert[3*(t-1)+ep1-1] = c;
2758  tri_vert[3*(u-1)+fp1-1] = i;
2759  r = tri_nabe[3*(t-1)+ep1-1];
2760  s = tri_nabe[3*(u-1)+fp1-1];
2761  tri_nabe[3*(t-1)+ep1-1] = u;
2762  tri_nabe[3*(u-1)+fp1-1] = t;
2763  tri_nabe[3*(t-1)+e-1] = s;
2764  tri_nabe[3*(u-1)+f-1] = r;
2765 
2766  if ( 0 < tri_nabe[3*(u-1)+fm1-1] )
2767  {
2768  *top = *top + 1;
2769  stack[(*top)-1] = u;
2770  }
2771 
2772  if ( 0 < s )
2773  {
2774  if ( tri_nabe[3*(s-1)+0] == u )
2775  {
2776  tri_nabe[3*(s-1)+0] = t;
2777  }
2778  else if ( tri_nabe[3*(s-1)+1] == u )
2779  {
2780  tri_nabe[3*(s-1)+1] = t;
2781  }
2782  else
2783  {
2784  tri_nabe[3*(s-1)+2] = t;
2785  }
2786 
2787  *top = *top + 1;
2788 
2789  if ( point_num < *top )
2790  {
2791  return 8;
2792  }
2793 
2794  stack[(*top)-1] = t;
2795  }
2796  else
2797  {
2798  if ( u == *btri && fp1 == *bedg )
2799  {
2800  *btri = t;
2801  *bedg = e;
2802  }
2803 
2804  l = - ( 3 * t + e - 1 );
2805  tt = t;
2806  ee = em1;
2807 
2808  while ( 0 < tri_nabe[3*(tt-1)+ee-1] )
2809  {
2810  tt = tri_nabe[3*(tt-1)+ee-1];
2811 
2812  if ( tri_vert[3*(tt-1)+0] == a )
2813  {
2814  ee = 3;
2815  }
2816  else if ( tri_vert[3*(tt-1)+1] == a )
2817  {
2818  ee = 1;
2819  }
2820  else
2821  {
2822  ee = 2;
2823  }
2824 
2825  }
2826 
2827  tri_nabe[3*(tt-1)+ee-1] = l;
2828 
2829  }
2830 
2831  if ( 0 < r )
2832  {
2833  if ( tri_nabe[3*(r-1)+0] == t )
2834  {
2835  tri_nabe[3*(r-1)+0] = u;
2836  }
2837  else if ( tri_nabe[3*(r-1)+1] == t )
2838  {
2839  tri_nabe[3*(r-1)+1] = u;
2840  }
2841  else
2842  {
2843  tri_nabe[3*(r-1)+2] = u;
2844  }
2845  }
2846  else
2847  {
2848  if ( t == *btri && ep1 == *bedg )
2849  {
2850  *btri = u;
2851  *bedg = f;
2852  }
2853 
2854  l = - ( 3 * u + f - 1 );
2855  tt = u;
2856  ee = fm1;
2857 
2858  while ( 0 < tri_nabe[3*(tt-1)+ee-1] )
2859  {
2860  tt = tri_nabe[3*(tt-1)+ee-1];
2861 
2862  if ( tri_vert[3*(tt-1)+0] == b )
2863  {
2864  ee = 3;
2865  }
2866  else if ( tri_vert[3*(tt-1)+1] == b )
2867  {
2868  ee = 1;
2869  }
2870  else
2871  {
2872  ee = 2;
2873  }
2874  }
2875  tri_nabe[3*(tt-1)+ee-1] = l;
2876  }
2877  }
2878  }
2879  return 0;
2880 }
2881 //**********************************************************************
2882 
2883 void timestamp ( void )
2884 
2885 //**********************************************************************
2886 //
2887 // Purpose:
2888 //
2889 // TIMESTAMP prints the current YMDHMS date as a time stamp.
2890 //
2891 // Example:
2892 //
2893 // May 31 2001 09:45:54 AM
2894 //
2895 // Modified:
2896 //
2897 // 21 August 2002
2898 //
2899 // Author:
2900 //
2901 // John Burkardt
2902 //
2903 // Parameters:
2904 //
2905 // None
2906 //
2907 {
2908 # define TIME_SIZE 29
2909 
2910  static char time_buffer[TIME_SIZE];
2911  const struct tm *tm;
2912  size_t len;
2913  time_t now;
2914 
2915  now = time ( NULL );
2916  tm = localtime ( &now );
2917 
2918  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
2919 
2920  if ( len != 0 )
2921  {
2922  cout << time_buffer << "\n";
2923  }
2924 
2925  return;
2926 # undef TIME_SIZE
2927 }
2928 //**********************************************************************
2929 
2930 char *timestring ( void )
2931 
2932 //**********************************************************************
2933 //
2934 // Purpose:
2935 //
2936 // TIMESTRING returns the current YMDHMS date as a string.
2937 //
2938 // Example:
2939 //
2940 // May 31 2001 09:45:54 AM
2941 //
2942 // Modified:
2943 //
2944 // 13 June 2003
2945 //
2946 // Author:
2947 //
2948 // John Burkardt
2949 //
2950 // Parameters:
2951 //
2952 // Output, char *TIMESTRING, a string containing the current YMDHMS date.
2953 //
2954 {
2955 # define TIME_SIZE 29
2956 
2957  const struct tm *tm;
2958  size_t len;
2959  time_t now;
2960  char *s;
2961 
2962  now = time ( NULL );
2963  tm = localtime ( &now );
2964 
2965  s = new char[TIME_SIZE];
2966 
2967  len = strftime ( s, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
2968 
2969  return s;
2970 # undef TIME_SIZE
2971 }
2972 //******************************************************************************
2973 
2974 double *triangle_circumcenter_2d ( double t[] )
2975 
2976 //******************************************************************************
2977 //
2978 // Purpose:
2979 //
2980 // TRIANGLE_CIRCUMCENTER_2D computes the circumcenter of a triangle in 2D.
2981 //
2982 // Discussion:
2983 //
2984 // The circumcenter of a triangle is the center of the circumcircle, the
2985 // circle that passes through the three vertices of the triangle.
2986 //
2987 // The circumcircle contains the triangle, but it is not necessarily the
2988 // smallest triangle to do so.
2989 //
2990 // If all angles of the triangle are no greater than 90 degrees, then
2991 // the center of the circumscribed circle will lie inside the triangle.
2992 // Otherwise, the center will lie outside the circle.
2993 //
2994 // The circumcenter is the intersection of the perpendicular bisectors
2995 // of the sides of the triangle.
2996 //
2997 // In geometry, the circumcenter of a triangle is often symbolized by "O".
2998 //
2999 // Modified:
3000 //
3001 // 09 February 2005
3002 //
3003 // Author:
3004 //
3005 // John Burkardt
3006 //
3007 // Parameters:
3008 //
3009 // Input, double T[2*3], the triangle vertices.
3010 //
3011 // Output, double *X, *Y, the coordinates of the circumcenter of the triangle.
3012 //
3013 {
3014 # define DIM_NUM 2
3015 
3016  double asq;
3017  double bot;
3018  double *center;
3019  double csq;
3020  double top1;
3021  double top2;
3022 
3023  center = new double[DIM_NUM];
3024 
3025  asq = ( t[0+1*2] - t[0+0*2] ) * ( t[0+1*2] - t[0+0*2] )
3026  + ( t[1+1*2] - t[1+0*2] ) * ( t[1+1*2] - t[1+0*2] );
3027 
3028  csq = ( t[0+2*2] - t[0+0*2] ) * ( t[0+2*2] - t[0+0*2] )
3029  + ( t[1+2*2] - t[1+0*2] ) * ( t[1+2*2] - t[1+0*2] );
3030 
3031  top1 = ( t[1+1*2] - t[1+0*2] ) * csq - ( t[1+2*2] - t[1+0*2] ) * asq;
3032  top2 = ( t[0+1*2] - t[0+0*2] ) * csq - ( t[0+2*2] - t[0+0*2] ) * asq;
3033 
3034  bot = ( t[1+1*2] - t[1+0*2] ) * ( t[0+2*2] - t[0+0*2] )
3035  - ( t[1+2*2] - t[1+0*2] ) * ( t[0+1*2] - t[0+0*2] );
3036 
3037  center[0] = t[0+0*2] + 0.5 * top1 / bot;
3038  center[1] = t[1+0*2] + 0.5 * top2 / bot;
3039 
3040  return center;
3041 
3042 # undef DIM_NUM
3043 }
3044 //******************************************************************************
3045 
3046 bool triangulation_plot_eps ( const char *file_out_name, int g_num, double g_xy[],
3047  int tri_num, int nod_tri[] )
3048 
3049 //******************************************************************************
3050 //
3051 // Purpose:
3052 //
3053 // TRIANGULATION_PLOT_EPS plots a triangulation of a pointset.
3054 //
3055 // Discussion:
3056 //
3057 // The triangulation is most usually a Delaunay triangulation,
3058 // but this is not necessary.
3059 //
3060 // The data can be generated by calling DTRIS2, but this is not
3061 // necessary.
3062 //
3063 // Modified:
3064 //
3065 // 08 September 2003
3066 //
3067 // Author:
3068 //
3069 // John Burkardt
3070 //
3071 // Parameters:
3072 //
3073 // Input, const char *FILE_OUT_NAME, the name of the output file.
3074 //
3075 // Input, int G_NUM, the number of points.
3076 //
3077 // Input, double G_XY[G_NUM,2], the coordinates of the points.
3078 //
3079 // Input, int TRI_NUM, the number of triangles.
3080 //
3081 // Input, int NOD_TRI[3,TRI_NUM], lists, for each triangle,
3082 // the indices of the points that form the vertices of the triangle.
3083 //
3084 // Output, bool TRIANGULATION_PLOT_EPS, is TRUE for success.
3085 //
3086 {
3087  char *date_time;
3088  int e;
3089  ofstream file_out;
3090  int g;
3091  int j;
3092  int k;
3093  int t;
3094  double x_max;
3095  double x_min;
3096  int x_ps;
3097  int x_ps_max = 576;
3098  int x_ps_max_clip = 594;
3099  int x_ps_min = 36;
3100  int x_ps_min_clip = 18;
3101  double y_max;
3102  double y_min;
3103  int y_ps;
3104  int y_ps_max = 666;
3105  int y_ps_max_clip = 684;
3106  int y_ps_min = 126;
3107  int y_ps_min_clip = 108;
3108 
3109  date_time = timestring ( );
3110 
3111  x_max = g_xy[0+0*2];
3112  x_min = g_xy[0+0*2];
3113  y_max = g_xy[1+0*2];
3114  y_min = g_xy[1+0*2];
3115 
3116  for ( g = 0; g < g_num; g++ )
3117  {
3118  x_max = d_max ( x_max, g_xy[0+g*2] );
3119  x_min = d_min ( x_min, g_xy[0+g*2] );
3120  y_max = d_max ( y_max, g_xy[1+g*2] );
3121  y_min = d_min ( y_min, g_xy[1+g*2] );
3122  }
3123 //
3124 // Plot the Delaunay triangulation.
3125 //
3126 //
3127 // Open the output file.
3128 //
3129  file_out.open ( file_out_name );
3130 
3131  if ( !file_out )
3132  {
3133  cout << "\n";
3134  cout << "TRIANGULATION_PLOT_EPS - Fatal error!\n";
3135  cout << " Cannot open the output file \"" << file_out_name << "\".\n";
3136  return false;
3137  }
3138 
3139  file_out << "%!PS-Adobe-3.0 EPSF-3.0\n";
3140  file_out << "%%Creator: triangulation_plot_eps.cc\n";
3141  file_out << "%%Title: " << file_out_name << "\n";
3142  file_out << "%%CreationDate: " << date_time << "\n";
3143  file_out << "%%Pages: 1\n";
3144  file_out << "%%Bounding Box: " << x_ps_min << " " << y_ps_min << " "
3145  << x_ps_max << " " << y_ps_max << "\n";
3146  file_out << "%%Document-Fonts: Times-Roman\n";
3147  file_out << "%%LanguageLevel: 1\n";
3148  file_out << "%%EndComments\n";
3149  file_out << "%%BeginProlog\n";
3150  file_out << "/inch {72 mul} def\n";
3151  file_out << "%%EndProlog\n";
3152  file_out << "%%Page: 1 1\n";
3153  file_out << "save\n";
3154  file_out << "%\n";
3155  file_out << "% Set the RGB line color to very light gray.\n";
3156  file_out << "%\n";
3157  file_out << "0.900 0.900 0.900 setrgbcolor\n";
3158  file_out << "%\n";
3159  file_out << "% Draw a gray border around the page.\n";
3160  file_out << "%\n";
3161  file_out << "newpath\n";
3162  file_out << " " << x_ps_min << " " << y_ps_min << " moveto\n";
3163  file_out << " " << x_ps_max << " " << y_ps_min << " lineto\n";
3164  file_out << " " << x_ps_max << " " << y_ps_max << " lineto\n";
3165  file_out << " " << x_ps_min << " " << y_ps_max << " lineto\n";
3166  file_out << " " << x_ps_min << " " << y_ps_min << " lineto\n";
3167  file_out << "stroke\n";
3168  file_out << "%\n";
3169  file_out << "% Set the RGB line color to black.\n";
3170  file_out << "%\n";
3171  file_out << "0.000 0.000 0.000 setrgbcolor\n";
3172  file_out << "%\n";
3173  file_out << "% Set the font and its size.\n";
3174  file_out << "%\n";
3175  file_out << "/Times-Roman findfont\n";
3176  file_out << "0.50 inch scalefont\n";
3177  file_out << "setfont\n";
3178  file_out << "%\n";
3179  file_out << "% Print a title.\n";
3180  file_out << "%\n";
3181  file_out << "210 702 moveto\n";
3182  file_out << "(Triangulation) show\n";
3183  file_out << "%\n";
3184  file_out << "% Define a clipping polygon.\n";
3185  file_out << "%\n";
3186  file_out << "newpath\n";
3187  file_out << " " << x_ps_min_clip << " " << y_ps_min_clip << " moveto\n";
3188  file_out << " " << x_ps_max_clip << " " << y_ps_min_clip << " lineto\n";
3189  file_out << " " << x_ps_max_clip << " " << y_ps_max_clip << " lineto\n";
3190  file_out << " " << x_ps_min_clip << " " << y_ps_max_clip << " lineto\n";
3191  file_out << " " << x_ps_min_clip << " " << y_ps_min_clip << " lineto\n";
3192  file_out << "clip newpath\n";
3193  file_out << "%\n";
3194  file_out << "% Set the RGB line color to green.\n";
3195  file_out << "%\n";
3196  file_out << "0.000 0.750 0.150 setrgbcolor\n";
3197  file_out << "%\n";
3198  file_out << "% Draw the nodes.\n";
3199  file_out << "%\n";
3200 
3201  for ( g = 0; g < g_num; g++ )
3202  {
3203  x_ps = int(
3204  ( ( x_max - g_xy[0+g*2] ) * double( x_ps_min )
3205  + ( g_xy[0+g*2] - x_min ) * double( x_ps_max ) )
3206  / ( x_max - x_min ) );
3207 
3208  y_ps = int(
3209  ( ( y_max - g_xy[1+g*2] ) * double( y_ps_min )
3210  + ( g_xy[1+g*2] - y_min ) * double( y_ps_max ) )
3211  / ( y_max - y_min ) );
3212 
3213  file_out << "newpath " << x_ps << " "
3214  << y_ps << " 5 0 360 arc closepath fill\n";
3215  }
3216 
3217  file_out << "%\n";
3218  file_out << "% Set the RGB line color to red.\n";
3219  file_out << "%\n";
3220  file_out << "0.900 0.200 0.100 setrgbcolor\n";
3221  file_out << "%\n";
3222  file_out << "% Draw the triangles.\n";
3223  file_out << "%\n";
3224 
3225  for ( t = 1; t <= tri_num; t++ )
3226  {
3227  file_out << "newpath\n";
3228 
3229  for ( j = 1; j <= 4; j++ )
3230  {
3231  e = i_wrap ( j, 1, 3 );
3232 
3233  k = nod_tri[3*(t-1)+e-1];
3234 
3235  x_ps = int(
3236  ( ( x_max - g_xy[0+(k-1)*2] ) * double( x_ps_min )
3237  + ( g_xy[0+(k-1)*2] - x_min ) * double( x_ps_max ) )
3238  / ( x_max - x_min ) );
3239 
3240  y_ps = int(
3241  ( ( y_max - g_xy[1+(k-1)*2] ) * double( y_ps_min )
3242  + ( g_xy[1+(k-1)*2] - y_min ) * double( y_ps_max ) )
3243  / ( y_max - y_min ) );
3244 
3245  if ( j == 1 )
3246  {
3247  file_out << x_ps << " " << y_ps << " moveto\n";
3248  }
3249  else
3250  {
3251  file_out << x_ps << " " << y_ps << " lineto\n";
3252  }
3253 
3254  }
3255 
3256  file_out << "stroke\n";
3257 
3258  }
3259 
3260  file_out << "restore showpage\n";
3261  file_out << "%\n";
3262  file_out << "% End of page.\n";
3263  file_out << "%\n";
3264  file_out << "%%Trailer\n";
3265  file_out << "%%EOF\n";
3266 
3267  file_out.close ( );
3268 
3269  return true;
3270 }
3271 //******************************************************************************
3272 
3273 void triangulation_print ( int point_num, double xc[], int tri_num,
3274  int tri_vert[], int tri_nabe[] )
3275 
3276 //******************************************************************************
3277 //
3278 // Purpose:
3279 //
3280 // TRIANGULATION_PRINT prints information defining a Delaunay triangulation.
3281 //
3282 // Discussion:
3283 //
3284 // Triangulations created by RTRIS include extra information encoded
3285 // in the negative values of TRI_NABE.
3286 //
3287 // Because some of the nodes counted in POINT_NUM may not actually be
3288 // used in the triangulation, I needed to compute the true number
3289 // of vertices. I added this calculation on 13 October 2001.
3290 //
3291 // Ernest Fasse pointed out an error in the indexing of VERTEX_LIST,
3292 // which was corrected on 19 February 2004.
3293 //
3294 // Modified:
3295 //
3296 // 19 February 2004
3297 //
3298 // Author:
3299 //
3300 // John Burkardt
3301 //
3302 // Parameters:
3303 //
3304 // Input, int POINT_NUM, the number of points.
3305 //
3306 // Input, double XC[2*POINT_NUM], the point coordinates.
3307 //
3308 // Input, int TRI_NUM, the number of triangles.
3309 //
3310 // Input, int TRI_VERT[3*TRI_NUM], the nodes that make up the triangles.
3311 //
3312 // Input, int TRI_NABE[3*TRI_NUM], the triangle neighbors on each side.
3313 // If there is no triangle neighbor on a particular side, the value of
3314 // TRI_NABE should be negative. If the triangulation data was created by
3315 // DTRIS2, then there is more information encoded in the negative values.
3316 //
3317 {
3318 # define DIM_NUM 2
3319 
3320  int boundary_num;
3321  int i;
3322  int j;
3323  int k;
3324  int n1;
3325  int n2;
3326  int s;
3327  int s1;
3328  int s2;
3329  bool skip;
3330  int t;
3331  int *vertex_list;
3332  int vertex_num;
3333 
3334  cout << "\n";
3335  cout << "TRIANGULATION_PRINT\n";
3336  cout << " Information defining a triangulation.\n";
3337  cout << "\n";
3338  cout << " The number of points is " << point_num << "\n";
3339 
3340  dmat_transpose_print ( DIM_NUM, point_num, xc, " Point coordinates" );
3341 
3342  cout << "\n";
3343  cout << " The number of triangles is " << tri_num << "\n";
3344  cout << "\n";
3345  cout << " Sets of three points are used as vertices of\n";
3346  cout << " the triangles. For each triangle, the points\n";
3347  cout << " are listed in counterclockwise order.\n";
3348 
3349  imat_transpose_print ( 3, tri_num, tri_vert, " Triangle nodes" );
3350 
3351  cout << "\n";
3352  cout << " On each side of a given triangle, there is either\n";
3353  cout << " another triangle, or a piece of the convex hull.\n";
3354  cout << " For each triangle, we list the indices of the three\n";
3355  cout << " neighbors, or (if negative) the codes of the\n";
3356  cout << " segments of the convex hull.\n";
3357 
3358  imat_transpose_print ( 3, tri_num, tri_nabe, " Triangle neighbors" );
3359 //
3360 // Determine VERTEX_NUM, the number of vertices. This is not
3361 // the same as the number of points!
3362 //
3363  vertex_list = new int[3*tri_num];
3364 
3365  k = 0;
3366  for ( t = 0; t < tri_num; t++ )
3367  {
3368  for ( s = 0; s < 3; s++ )
3369  {
3370  vertex_list[k] = tri_vert[s+t*3];
3371  k = k + 1;
3372  }
3373  }
3374 
3375  ivec_sort_heap_a ( 3*tri_num, vertex_list );
3376 
3377  ivec_sorted_unique ( 3*tri_num, vertex_list, &vertex_num );
3378 
3379  delete [] vertex_list;
3380 //
3381 // Determine the number of boundary points.
3382 //
3383  boundary_num = 2 * vertex_num - tri_num - 2;
3384 
3385  cout << "\n";
3386  cout << " The number of boundary points is " << boundary_num << "\n";
3387  cout << "\n";
3388  cout << " The segments that make up the convex hull can be\n";
3389  cout << " determined from the negative entries of the triangle\n";
3390  cout << " neighbor list.\n";
3391  cout << "\n";
3392  cout << " # Tri Side N1 N2\n";
3393  cout << "\n";
3394 
3395  skip = false;
3396 
3397  k = 0;
3398 
3399  for ( i = 0; i < tri_num; i++ )
3400  {
3401  for ( j = 0; j < 3; j++ )
3402  {
3403  if ( tri_nabe[j+i*3] < 0 )
3404  {
3405  s = -tri_nabe[j+i*3];
3406  t = s / 3;
3407 
3408  if ( t < 1 || tri_num < t )
3409  {
3410  cout << "\n";
3411  cout << " Sorry, this data does not use the DTRIS2\n";
3412  cout << " convention for convex hull segments.\n";
3413  skip = true;
3414  break;
3415  }
3416 
3417  s1 = ( s % 3 ) + 1;
3418  s2 = i_wrap ( s1+1, 1, 3 );
3419  k = k + 1;
3420  n1 = tri_vert[s1-1+(t-1)*3];
3421  n2 = tri_vert[s2-1+(t-1)*3];
3422  cout << setw(4) << k << " "
3423  << setw(4) << t << " "
3424  << setw(4) << s1 << " "
3425  << setw(4) << n1 << " "
3426  << setw(4) << n2 << "\n";
3427  }
3428 
3429  }
3430 
3431  if ( skip )
3432  {
3433  break;
3434  }
3435 
3436  }
3437 
3438  return;
3439 # undef DIM_NUM
3440 }
3441 //******************************************************************************
3442 
3443 void vbedg ( double x, double y, int point_num, double point_xy[], int tri_num,
3444  int tri_vert[], int tri_nabe[], int *ltri, int *ledg, int *rtri, int *redg )
3445 
3446 //******************************************************************************
3447 //
3448 // Purpose:
3449 //
3450 // VBEDG determines which boundary edges are visible to a point.
3451 //
3452 // Discussion:
3453 //
3454 // The point (X,Y) is assumed to be outside the convex hull of the
3455 // region covered by the 2D triangulation.
3456 //
3457 // Author:
3458 //
3459 // Barry Joe,
3460 // Department of Computing Science,
3461 // University of Alberta,
3462 // Edmonton, Alberta, Canada T6G 2H1
3463 //
3464 // Reference:
3465 //
3466 // Barry Joe,
3467 // GEOMPACK - a software package for the generation of meshes
3468 // using geometric algorithms,
3469 // Advances in Engineering Software,
3470 // Volume 13, pages 325-331, 1991.
3471 //
3472 // Modified:
3473 //
3474 // 02 September 2003
3475 //
3476 // Parameters:
3477 //
3478 // Input, double X, Y, the coordinates of a point outside the convex hull
3479 // of the current triangulation.
3480 //
3481 // Input, int POINT_NUM, the number of points.
3482 //
3483 // Input, double POINT_XY[POINT_NUM*2], the coordinates of the vertices.
3484 //
3485 // Input, int TRI_NUM, the number of triangles.
3486 //
3487 // Input, int TRI_VERT[TRI_NUM*3], the triangle incidence list.
3488 //
3489 // Input, int TRI_NABE[TRI_NUM*3], the triangle neighbor list; negative
3490 // values are used for links of a counter clockwise linked list of boundary
3491 // edges;
3492 // LINK = -(3*I + J-1) where I, J = triangle, edge index.
3493 //
3494 // Input/output, int *LTRI, *LEDG. If LTRI != 0 then these values are
3495 // assumed to be already computed and are not changed, else they are updated.
3496 // On output, LTRI is the index of boundary triangle to the left of the
3497 // leftmost boundary triangle visible from (X,Y), and LEDG is the boundary
3498 // edge of triangle LTRI to the left of the leftmost boundary edge visible
3499 // from (X,Y). 1 <= LEDG <= 3.
3500 //
3501 // Input/output, int *RTRI. On input, the index of the boundary triangle
3502 // to begin the search at. On output, the index of the rightmost boundary
3503 // triangle visible from (X,Y).
3504 //
3505 // Input/output, int *REDG, the edge of triangle RTRI that is visible
3506 // from (X,Y). 1 <= REDG <= 3.
3507 //
3508 {
3509  int a;
3510  double ax;
3511  double ay;
3512  int b;
3513  double bx;
3514  double by;
3515  bool done;
3516  int e;
3517  int l;
3518  int lr;
3519  int t;
3520 //
3521 // Find the rightmost visible boundary edge using links, then possibly
3522 // leftmost visible boundary edge using triangle neighbor information.
3523 //
3524  if ( *ltri == 0 )
3525  {
3526  done = false;
3527  *ltri = *rtri;
3528  *ledg = *redg;
3529  }
3530  else
3531  {
3532  done = true;
3533  }
3534 
3535  for ( ; ; )
3536  {
3537  l = -tri_nabe[3*((*rtri)-1)+(*redg)-1];
3538  t = l / 3;
3539  e = 1 + l % 3;
3540  a = tri_vert[3*(t-1)+e-1];
3541 
3542  if ( e <= 2 )
3543  {
3544  b = tri_vert[3*(t-1)+e];
3545  }
3546  else
3547  {
3548  b = tri_vert[3*(t-1)+0];
3549  }
3550 
3551  ax = point_xy[2*(a-1)+0];
3552  ay = point_xy[2*(a-1)+1];
3553 
3554  bx = point_xy[2*(b-1)+0];
3555  by = point_xy[2*(b-1)+1];
3556 
3557  lr = lrline ( x, y, ax, ay, bx, by, 0.0 );
3558 
3559  if ( lr <= 0 )
3560  {
3561  break;
3562  }
3563 
3564  *rtri = t;
3565  *redg = e;
3566 
3567  }
3568 
3569  if ( done )
3570  {
3571  return;
3572  }
3573 
3574  t = *ltri;
3575  e = *ledg;
3576 
3577  for ( ; ; )
3578  {
3579  b = tri_vert[3*(t-1)+e-1];
3580  e = i_wrap ( e-1, 1, 3 );
3581 
3582  while ( 0 < tri_nabe[3*(t-1)+e-1] )
3583  {
3584  t = tri_nabe[3*(t-1)+e-1];
3585 
3586  if ( tri_vert[3*(t-1)+0] == b )
3587  {
3588  e = 3;
3589  }
3590  else if ( tri_vert[3*(t-1)+1] == b )
3591  {
3592  e = 1;
3593  }
3594  else
3595  {
3596  e = 2;
3597  }
3598 
3599  }
3600 
3601  a = tri_vert[3*(t-1)+e-1];
3602  ax = point_xy[2*(a-1)+0];
3603  ay = point_xy[2*(a-1)+1];
3604 
3605  bx = point_xy[2*(b-1)+0];
3606  by = point_xy[2*(b-1)+1];
3607 
3608  lr = lrline ( x, y, ax, ay, bx, by, 0.0 );
3609 
3610  if ( lr <= 0 )
3611  {
3612  break;
3613  }
3614 
3615  }
3616 
3617  *ltri = t;
3618  *ledg = e;
3619 
3620  return;
3621 }