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