9 #if defined(__APPLE__) && defined(__clang__) 10 #pragma clang fp exceptions(ignore) 191 cout <<
"D2VEC_PART_QUICK_A - Fatal error!\n";
212 for ( i = 2; i <=
n; i++ )
234 for ( i = 0; i < ll - m; i++ )
236 for ( j = 0; j < 2; j++ )
238 a[2*i+j] = a[2*(i+m)+j];
244 for ( i = ll; i < ll+m; i++ )
246 for ( j = 0; j < 2; j++ )
320 cout <<
"D2VEC_PERMUTE - Fatal error!\n";
321 cout <<
" The input array does not represent\n";
322 cout <<
" a proper permutation.\n";
328 for ( istart = 1; istart <=
n; istart++ )
330 if (
p[istart-1] < 0 )
334 else if (
p[istart-1] == istart )
336 p[istart-1] = -
p[istart-1];
341 a_temp[0] = a[0+(istart-1)*2];
342 a_temp[1] = a[1+(istart-1)*2];
352 p[iput-1] = -
p[iput-1];
354 if ( iget < 1 ||
n < iget )
357 cout <<
"D2VEC_PERMUTE - Fatal error!\n";
361 if ( iget == istart )
363 a[0+(iput-1)*2] = a_temp[0];
364 a[1+(iput-1)*2] = a_temp[1];
367 a[0+(iput-1)*2] = a[0+(iget-1)*2];
368 a[1+(iput-1)*2] = a[1+(iget-1)*2];
375 for ( i = 0; i <
n; i++ )
460 aval[0] = a[0+(indxt-1)*2];
461 aval[1] = a[1+(indxt-1)*2];
466 aval[0] = a[0+(indxt-1)*2];
467 aval[1] = a[1+(indxt-1)*2];
468 indx[ir-1] = indx[0];
486 if ( a[0+(indx[j-1]-1)*2] < a[0+(indx[j]-1)*2] ||
487 ( a[0+(indx[j-1]-1)*2] == a[0+(indx[j]-1)*2] &&
488 a[1+(indx[j-1]-1)*2] < a[1+(indx[j]-1)*2] ) )
494 if ( aval[0] < a[0+(indx[j-1]-1)*2] ||
495 ( aval[0] == a[0+(indx[j-1]-1)*2] &&
496 aval[1] < a[1+(indx[j-1]-1)*2] ) )
498 indx[i-1] = indx[j-1];
544 # define LEVEL_MAX 25 556 cout <<
"D2VEC_SORT_QUICK_A - Fatal error!\n";
567 rsave[level-1] =
n + 1;
571 while ( 0 < n_segment )
585 cout <<
"D2VEC_SORT_QUICK_A - Fatal error!\n";
586 cout <<
" Exceeding recursion maximum of " <<
LEVEL_MAX <<
"\n";
591 n_segment = l_segment;
592 rsave[level-1] = r_segment + base - 1;
598 else if ( r_segment < n_segment )
600 n_segment = n_segment + 1 - r_segment;
601 base = base + r_segment - 1;
616 base = rsave[level-1];
617 n_segment = rsave[level-2] - rsave[level-1];
632 int diaedg (
double x0,
double y0,
double x1,
double y1,
double x2,
double y2,
633 double x3,
double y3 )
705 tola = tol *
d_max ( fabs ( dx10 ),
706 d_max ( fabs ( dy10 ),
707 d_max ( fabs ( dx30 ), fabs ( dy30 ) ) ) );
709 tolb = tol *
d_max ( fabs ( dx12 ),
710 d_max ( fabs ( dy12 ),
711 d_max ( fabs ( dx32 ), fabs ( dy32 ) ) ) );
713 ca = dx10 * dx30 + dy10 * dy30;
714 cb = dx12 * dx32 + dy12 * dy32;
716 if ( tola < ca && tolb < cb )
720 else if ( ca < -tola && cb < -tolb )
726 tola =
d_max ( tola, tolb );
727 s = ( dx10 * dy30 - dx30 * dy10 ) * cb
728 + ( dx32 * dy12 - dx12 * dy32 ) * ca;
734 else if (
s < -tola )
781 int ihi,
int jhi,
const char *title )
824 cout << title <<
"\n";
827 for ( i2lo =
i_max ( ilo, 1 ); i2lo <=
i_min ( ihi, m ); i2lo = i2lo +
INCX )
829 i2hi = i2lo +
INCX - 1;
830 i2hi =
i_min ( i2hi, m );
831 i2hi =
i_min ( i2hi, ihi );
833 inc = i2hi + 1 - i2lo;
837 for ( i = i2lo; i <= i2hi; i++ )
839 cout <<
setw(7) << i <<
" ";
845 j2lo =
i_max ( jlo, 1 );
848 for ( j = j2lo; j <= j2hi; j++ )
850 cout <<
setw(5) << j <<
" ";
851 for ( i2 = 1; i2 <= inc; i2++ )
854 cout <<
setw(14) << a[(i-1)+(j-1)*m];
929 for ( j = 0; j <
n; j++ )
931 for ( i = 0; i < m; i++ )
935 *seed = 16807 * ( *seed -
k * 127773 ) -
k * 2836;
939 *seed = *seed + 2147483647;
945 r[i+j*m] =
b + (
c -
b ) *
double( *seed ) * 4.656612875E-10;
953 int dtris2 (
int point_num,
double point_xy[],
int *tri_num,
954 int tri_vert[],
int tri_nabe[] )
1034 stack =
new int[point_num];
1048 for ( i = 2; i <= point_num; i++ )
1055 for ( j = 0; j <= 1; j++ )
1057 cmax =
d_max ( fabs ( point_xy[2*(m-1)+j] ),
1058 fabs ( point_xy[2*(m1-1)+j] ) );
1060 if ( tol * ( cmax + 1.0 )
1061 < fabs ( point_xy[2*(m-1)+j] - point_xy[2*(m1-1)+j] ) )
1072 cout <<
"DTRIS2 - Fatal error!\n";
1073 cout <<
" Fails for point number I = " << i <<
"\n";
1074 cout <<
" M = " << m <<
"\n";
1075 cout <<
" M1 = " << m1 <<
"\n";
1076 cout <<
" X,Y(M) = " << point_xy[2*(m-1)+0] <<
" " 1077 << point_xy[2*(m-1)+1] <<
"\n";
1078 cout <<
" X,Y(M1) = " << point_xy[2*(m1-1)+0] <<
" " 1079 << point_xy[2*(m1-1)+1] <<
"\n";
1095 if ( point_num < j )
1098 cout <<
"DTRIS2 - Fatal error!\n";
1105 lr =
lrline ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1],
1106 point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1],
1107 point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0 );
1125 tri_vert[3*0+0] = m1;
1126 tri_vert[3*0+1] = m2;
1127 tri_vert[3*0+2] = m;
1128 tri_nabe[3*0+2] = -3;
1130 for ( i = 2; i <= *tri_num; i++ )
1134 tri_vert[3*(i-1)+0] = m1;
1135 tri_vert[3*(i-1)+1] = m2;
1136 tri_vert[3*(i-1)+2] = m;
1137 tri_nabe[3*(i-2)+0] = -3 * i;
1138 tri_nabe[3*(i-2)+1] = i;
1139 tri_nabe[3*(i-1)+2] = i - 1;
1143 tri_nabe[3*(*tri_num-1)+0] = -3 * (*tri_num) - 1;
1144 tri_nabe[3*(*tri_num-1)+1] = -5;
1150 tri_vert[3*0+0] = m2;
1151 tri_vert[3*0+1] = m1;
1152 tri_vert[3*0+2] = m;
1153 tri_nabe[3*0+0] = -4;
1155 for ( i = 2; i <= *tri_num; i++ )
1159 tri_vert[3*(i-1)+0] = m2;
1160 tri_vert[3*(i-1)+1] = m1;
1161 tri_vert[3*(i-1)+2] = m;
1162 tri_nabe[3*(i-2)+2] = i;
1163 tri_nabe[3*(i-1)+0] = -3 * i - 3;
1164 tri_nabe[3*(i-1)+1] = i - 1;
1167 tri_nabe[3*(*tri_num-1)+2] = -3 * (*tri_num);
1168 tri_nabe[3*0+1] = -3 * (*tri_num) - 2;
1179 for ( i = j+1; i <= point_num; i++ )
1182 m1 = tri_vert[3*(ltri-1)+ledg-1];
1186 m2 = tri_vert[3*(ltri-1)+ledg];
1190 m2 = tri_vert[3*(ltri-1)+0];
1193 lr =
lrline ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1],
1194 point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1],
1195 point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0 );
1205 l = -tri_nabe[3*(ltri-1)+ledg-1];
1210 vbedg ( point_xy[2*(m-1)+0], point_xy[2*(m-1)+1], point_num,
1211 point_xy, *tri_num, tri_vert, tri_nabe, <ri, &ledg, &rtri, &redg );
1214 l = -tri_nabe[3*(ltri-1)+ledg-1];
1220 l = -tri_nabe[3*(t-1)+
e-1];
1221 m2 = tri_vert[3*(t-1)+
e-1];
1225 m1 = tri_vert[3*(t-1)+
e];
1229 m1 = tri_vert[3*(t-1)+0];
1232 *tri_num = *tri_num + 1;
1233 tri_nabe[3*(t-1)+
e-1] = *tri_num;
1234 tri_vert[3*(*tri_num-1)+0] = m1;
1235 tri_vert[3*(*tri_num-1)+1] = m2;
1236 tri_vert[3*(*tri_num-1)+2] = m;
1237 tri_nabe[3*(*tri_num-1)+0] = t;
1238 tri_nabe[3*(*tri_num-1)+1] = *tri_num - 1;
1239 tri_nabe[3*(*tri_num-1)+2] = *tri_num + 1;
1242 if ( point_num < top )
1245 cout <<
"DTRIS2 - Fatal error!\n";
1246 cout <<
" Stack overflow.\n";
1251 stack[top-1] = *tri_num;
1253 if ( t == rtri &&
e == redg )
1260 tri_nabe[3*(ltri-1)+ledg-1] = -3 *
n - 1;
1261 tri_nabe[3*(
n-1)+1] = -3 * (*tri_num) - 2;
1262 tri_nabe[3*(*tri_num-1)+2] = -l;
1266 error =
swapec ( m, &top, <ri, &ledg, point_num, point_xy, *tri_num,
1267 tri_vert, tri_nabe, stack );
1272 cout <<
"DTRIS2 - Fatal error!\n";
1273 cout <<
" Error return from SWAPEC.\n";
1282 for ( i = 0; i < 3; i++ )
1284 for ( j = 0; j < *tri_num; j++ )
1286 tri_vert[i+j*3] = indx [ tri_vert[i+j*3] - 1 ];
1330 for ( i = 0; i <
n; i++ )
1332 if ( a1[i] != a2[i] )
1378 for ( i = 0; i <
n; i++ )
1381 if ( a2[i] < a1[i] )
1385 else if ( a1[i] < a2[i] )
1432 for ( i = 0; i <
n; i++ )
1434 if ( a1[i] < a2[i] )
1438 else if ( a2[i] < a1[i] )
1480 cout << title <<
"\n";
1484 for ( i = 0; i <=
n-1; i++ )
1486 cout <<
setw(6) << i + 1 <<
" " 1487 <<
setw(14) << a[i] <<
"\n";
1520 for ( i = 0; i <
n; i++ )
1660 cout <<
"I_MODP - Fatal error!\n";
1661 cout <<
" I_MODP ( I, J ) called with J = " << j <<
"\n";
1669 value = value + abs ( j );
1769 jlo =
i_min ( ilo, ihi );
1770 jhi =
i_max ( ilo, ihi );
1772 wide = jhi + 1 - jlo;
1780 value = jlo +
i_modp ( ival - jlo, wide );
1821 int ihi,
int jhi,
const char *title )
1864 cout << title <<
"\n";
1869 for ( i2lo = ilo; i2lo <= ihi; i2lo = i2lo +
INCX )
1871 i2hi = i2lo +
INCX - 1;
1872 i2hi =
i_min ( i2hi, m );
1873 i2hi =
i_min ( i2hi, ihi );
1882 for ( i = i2lo; i <= i2hi; i++ )
1884 cout <<
setw(7) << i <<
" ";
1892 j2lo =
i_max ( jlo, 1 );
1895 for ( j = j2lo; j <= j2hi; j++ )
1900 cout <<
setw(5) << j <<
" ";
1901 for ( i = i2lo; i <= i2hi; i++ )
1903 cout <<
setw(6) << a[i-1+(j-1)*m] <<
" ";
1972 for ( i = (
n/2)-1; 0 <= i; i-- )
2006 if ( a[m] < a[m+1] )
2069 for ( i = 0; i <
n; i++ )
2133 for ( n1 =
n-1; 2 <= n1; n1-- )
2189 for ( i = 1; i <
n; i++ )
2191 if ( a[i] != a[*nuniq] )
2193 *nuniq = *nuniq + 1;
2203 int lrline (
double xu,
double yu,
double xv1,
double yv1,
double xv2,
2204 double yv2,
double dv )
2255 double tol = 0.0000001;
2264 tolabs = tol *
d_max ( fabs ( dx ),
2265 d_max ( fabs ( dy ),
2266 d_max ( fabs ( dxu ),
2267 d_max ( fabs ( dyu ), fabs ( dv ) ) ) ) );
2269 t = dy * dxu - dx * dyu + dv *
sqrt ( dx * dx + dy * dy );
2275 else if ( -tolabs <= t )
2279 else if ( t < -tolabs )
2322 for ( seek = 1; seek <=
n; seek++ )
2326 for ( i = 0; i <
n; i++ )
2375 cout <<
"PERM_INV - Fatal error!\n";
2376 cout <<
" Input value of N = " <<
n <<
"\n";
2383 cout <<
"PERM_INV - Fatal error!\n";
2384 cout <<
" The input array does not represent\n";
2385 cout <<
" a proper permutation.\n";
2391 for ( i = 1; i <=
n; i++ )
2403 p[i-1] =
i_sign ( is ) * abs (
p[i-1] );
2406 for ( i = 1; i <=
n; i++ )
2497 z =
new double [
n ];
2499 for ( i = 0; i <
n; i++ )
2501 z[i] =
p[0+i*2] *
p[0+i*2] +
p[1+i*2] *
p[1+i*2];
2507 for ( pass = 1; pass <= 2; pass++ )
2511 tri =
new int[3*
count];
2517 for ( i = 0; i <
n - 2; i++ )
2519 for ( j = i+1; j <
n; j++ )
2521 for (
k = i+1;
k <
n;
k++ )
2525 xn = (
p[1+j*2] -
p[1+i*2] ) * ( z[
k] - z[i] )
2526 - (
p[1+
k*2] -
p[1+i*2] ) * ( z[j] - z[i] );
2527 yn = (
p[0+
k*2] -
p[0+i*2] ) * ( z[j] - z[i] )
2528 - (
p[0+j*2] -
p[0+i*2] ) * ( z[
k] - z[i] );
2529 zn = (
p[0+j*2] -
p[0+i*2] ) * (
p[1+
k*2] -
p[1+i*2] )
2530 - (
p[0+
k*2] -
p[0+i*2] ) * (
p[1+j*2] -
p[1+i*2] );
2536 for ( m = 0; m <
n; m++ )
2538 flag = flag && ( (
p[0+m*2] -
p[0+i*2] ) * xn
2539 + (
p[1+m*2] -
p[1+i*2] ) *
yn 2540 + ( z[m] - z[i] ) * zn <= 0 );
2596 t =
const_cast<char*
>(
s) +
n - 1;
2612 int swapec (
int i,
int *top,
int *btri,
int *bedg,
int point_num,
2613 double point_xy[],
int tri_num,
int tri_vert[],
int tri_nabe[],
2704 x = point_xy[2*(i-1)+0];
2705 y = point_xy[2*(i-1)+1];
2714 t = stack[(*top)-1];
2717 if ( tri_vert[3*(t-1)+0] == i )
2720 b = tri_vert[3*(t-1)+2];
2722 else if ( tri_vert[3*(t-1)+1] == i )
2725 b = tri_vert[3*(t-1)+0];
2730 b = tri_vert[3*(t-1)+1];
2733 a = tri_vert[3*(t-1)+
e-1];
2734 u = tri_nabe[3*(t-1)+
e-1];
2736 if ( tri_nabe[3*(u-1)+0] == t )
2739 c = tri_vert[3*(u-1)+2];
2741 else if ( tri_nabe[3*(u-1)+1] == t )
2744 c = tri_vert[3*(u-1)+0];
2749 c = tri_vert[3*(u-1)+1];
2753 point_xy[2*(a-1)+0], point_xy[2*(a-1)+1],
2754 point_xy[2*(
c-1)+0], point_xy[2*(
c-1)+1],
2755 point_xy[2*(
b-1)+0], point_xy[2*(
b-1)+1] );
2764 tri_vert[3*(t-1)+ep1-1] =
c;
2765 tri_vert[3*(u-1)+fp1-1] = i;
2766 r = tri_nabe[3*(t-1)+ep1-1];
2767 s = tri_nabe[3*(u-1)+fp1-1];
2768 tri_nabe[3*(t-1)+ep1-1] = u;
2769 tri_nabe[3*(u-1)+fp1-1] = t;
2770 tri_nabe[3*(t-1)+
e-1] =
s;
2771 tri_nabe[3*(u-1)+
f-1] = r;
2773 if ( 0 < tri_nabe[3*(u-1)+fm1-1] )
2776 stack[(*top)-1] = u;
2781 if ( tri_nabe[3*(
s-1)+0] == u )
2783 tri_nabe[3*(
s-1)+0] = t;
2785 else if ( tri_nabe[3*(
s-1)+1] == u )
2787 tri_nabe[3*(
s-1)+1] = t;
2791 tri_nabe[3*(
s-1)+2] = t;
2796 if ( point_num < *top )
2801 stack[(*top)-1] = t;
2805 if ( u == *btri && fp1 == *bedg )
2811 l = - ( 3 * t +
e - 1 );
2815 while ( 0 < tri_nabe[3*(tt-1)+ee-1] )
2817 tt = tri_nabe[3*(tt-1)+ee-1];
2819 if ( tri_vert[3*(tt-1)+0] == a )
2823 else if ( tri_vert[3*(tt-1)+1] == a )
2834 tri_nabe[3*(tt-1)+ee-1] = l;
2840 if ( tri_nabe[3*(r-1)+0] == t )
2842 tri_nabe[3*(r-1)+0] = u;
2844 else if ( tri_nabe[3*(r-1)+1] == t )
2846 tri_nabe[3*(r-1)+1] = u;
2850 tri_nabe[3*(r-1)+2] = u;
2855 if ( t == *btri && ep1 == *bedg )
2861 l = - ( 3 * u +
f - 1 );
2865 while ( 0 < tri_nabe[3*(tt-1)+ee-1] )
2867 tt = tri_nabe[3*(tt-1)+ee-1];
2869 if ( tri_vert[3*(tt-1)+0] ==
b )
2873 else if ( tri_vert[3*(tt-1)+1] ==
b )
2882 tri_nabe[3*(tt-1)+ee-1] = l;
2915 # define TIME_SIZE 29 2918 const struct tm *tm;
2922 now = time ( NULL );
2923 tm = localtime ( &now );
2925 len = strftime ( time_buffer,
TIME_SIZE,
"%d %B %Y %I:%M:%S %p", tm );
2929 cout << time_buffer <<
"\n";
2962 # define TIME_SIZE 29 2964 const struct tm *tm;
2968 now = time ( NULL );
2969 tm = localtime ( &now );
2973 strftime (
s,
TIME_SIZE,
"%d %B %Y %I:%M:%S %p", tm );
3032 asq = ( t[0+1*2] - t[0+0*2] ) * ( t[0+1*2] - t[0+0*2] )
3033 + ( t[1+1*2] - t[1+0*2] ) * ( t[1+1*2] - t[1+0*2] );
3035 csq = ( t[0+2*2] - t[0+0*2] ) * ( t[0+2*2] - t[0+0*2] )
3036 + ( t[1+2*2] - t[1+0*2] ) * ( t[1+2*2] - t[1+0*2] );
3038 top1 = ( t[1+1*2] - t[1+0*2] ) * csq - ( t[1+2*2] - t[1+0*2] ) * asq;
3039 top2 = ( t[0+1*2] - t[0+0*2] ) * csq - ( t[0+2*2] - t[0+0*2] ) * asq;
3041 bot = ( t[1+1*2] - t[1+0*2] ) * ( t[0+2*2] - t[0+0*2] )
3042 - ( t[1+2*2] - t[1+0*2] ) * ( t[0+1*2] - t[0+0*2] );
3044 center[0] = t[0+0*2] + 0.5 * top1 / bot;
3045 center[1] = t[1+0*2] + 0.5 * top2 / bot;
3054 double g_xy[],
int tri_num,
int nod_tri[] )
3105 int x_ps_max_clip = 594;
3107 int x_ps_min_clip = 18;
3112 int y_ps_max_clip = 684;
3114 int y_ps_min_clip = 108;
3118 x_max = g_xy[0+0*2];
3119 x_min = g_xy[0+0*2];
3120 y_max = g_xy[1+0*2];
3121 y_min = g_xy[1+0*2];
3123 for (
g = 0;
g < g_num;
g++ )
3125 x_max =
d_max ( x_max, g_xy[0+
g*2] );
3126 x_min =
d_min ( x_min, g_xy[0+
g*2] );
3127 y_max =
d_max ( y_max, g_xy[1+
g*2] );
3128 y_min =
d_min ( y_min, g_xy[1+
g*2] );
3136 file_out.open ( file_out_name );
3141 cout <<
"TRIANGULATION_PLOT_EPS - Fatal error!\n";
3142 cout <<
" Cannot open the output file \"" << file_out_name <<
"\".\n";
3146 file_out <<
"%!PS-Adobe-3.0 EPSF-3.0\n";
3147 file_out <<
"%%Creator: triangulation_plot_eps.cc\n";
3148 file_out <<
"%%Title: " << file_out_name <<
"\n";
3149 file_out <<
"%%CreationDate: " << date_time <<
"\n";
3150 file_out <<
"%%Pages: 1\n";
3151 file_out <<
"%%Bounding Box: " << x_ps_min <<
" " << y_ps_min <<
" " 3152 << x_ps_max <<
" " << y_ps_max <<
"\n";
3153 file_out <<
"%%Document-Fonts: Times-Roman\n";
3154 file_out <<
"%%LanguageLevel: 1\n";
3155 file_out <<
"%%EndComments\n";
3156 file_out <<
"%%BeginProlog\n";
3157 file_out <<
"/inch {72 mul} def\n";
3158 file_out <<
"%%EndProlog\n";
3159 file_out <<
"%%Page: 1 1\n";
3160 file_out <<
"save\n";
3162 file_out <<
"% Set the RGB line color to very light gray.\n";
3164 file_out <<
"0.900 0.900 0.900 setrgbcolor\n";
3166 file_out <<
"% Draw a gray border around the page.\n";
3168 file_out <<
"newpath\n";
3169 file_out <<
" " << x_ps_min <<
" " << y_ps_min <<
" moveto\n";
3170 file_out <<
" " << x_ps_max <<
" " << y_ps_min <<
" lineto\n";
3171 file_out <<
" " << x_ps_max <<
" " << y_ps_max <<
" lineto\n";
3172 file_out <<
" " << x_ps_min <<
" " << y_ps_max <<
" lineto\n";
3173 file_out <<
" " << x_ps_min <<
" " << y_ps_min <<
" lineto\n";
3174 file_out <<
"stroke\n";
3176 file_out <<
"% Set the RGB line color to black.\n";
3178 file_out <<
"0.000 0.000 0.000 setrgbcolor\n";
3180 file_out <<
"% Set the font and its size.\n";
3182 file_out <<
"/Times-Roman findfont\n";
3183 file_out <<
"0.50 inch scalefont\n";
3184 file_out <<
"setfont\n";
3186 file_out <<
"% Print a title.\n";
3188 file_out <<
"210 702 moveto\n";
3189 file_out <<
"(Triangulation) show\n";
3191 file_out <<
"% Define a clipping polygon.\n";
3193 file_out <<
"newpath\n";
3194 file_out <<
" " << x_ps_min_clip <<
" " << y_ps_min_clip <<
" moveto\n";
3195 file_out <<
" " << x_ps_max_clip <<
" " << y_ps_min_clip <<
" lineto\n";
3196 file_out <<
" " << x_ps_max_clip <<
" " << y_ps_max_clip <<
" lineto\n";
3197 file_out <<
" " << x_ps_min_clip <<
" " << y_ps_max_clip <<
" lineto\n";
3198 file_out <<
" " << x_ps_min_clip <<
" " << y_ps_min_clip <<
" lineto\n";
3199 file_out <<
"clip newpath\n";
3201 file_out <<
"% Set the RGB line color to green.\n";
3203 file_out <<
"0.000 0.750 0.150 setrgbcolor\n";
3205 file_out <<
"% Draw the nodes.\n";
3208 for (
g = 0;
g < g_num;
g++ )
3211 ( ( x_max - g_xy[0+
g*2] ) *
double( x_ps_min )
3212 + ( g_xy[0+
g*2] - x_min ) *
double( x_ps_max ) )
3213 / ( x_max - x_min ) );
3216 ( ( y_max - g_xy[1+
g*2] ) *
double( y_ps_min )
3217 + ( g_xy[1+
g*2] - y_min ) *
double( y_ps_max ) )
3218 / ( y_max - y_min ) );
3220 file_out <<
"newpath " << x_ps <<
" " 3221 << y_ps <<
" 5 0 360 arc closepath fill\n";
3225 file_out <<
"% Set the RGB line color to red.\n";
3227 file_out <<
"0.900 0.200 0.100 setrgbcolor\n";
3229 file_out <<
"% Draw the triangles.\n";
3232 for ( t = 1; t <= tri_num; t++ )
3234 file_out <<
"newpath\n";
3236 for ( j = 1; j <= 4; j++ )
3240 k = nod_tri[3*(t-1)+
e-1];
3243 ( ( x_max - g_xy[0+(
k-1)*2] ) *
double( x_ps_min )
3244 + ( g_xy[0+(
k-1)*2] - x_min ) *
double( x_ps_max ) )
3245 / ( x_max - x_min ) );
3248 ( ( y_max - g_xy[1+(
k-1)*2] ) *
double( y_ps_min )
3249 + ( g_xy[1+(
k-1)*2] - y_min ) *
double( y_ps_max ) )
3250 / ( y_max - y_min ) );
3254 file_out << x_ps <<
" " << y_ps <<
" moveto\n";
3258 file_out << x_ps <<
" " << y_ps <<
" lineto\n";
3263 file_out <<
"stroke\n";
3267 file_out <<
"restore showpage\n";
3269 file_out <<
"% End of page.\n";
3271 file_out <<
"%%Trailer\n";
3272 file_out <<
"%%EOF\n";
3281 int tri_vert[],
int tri_nabe[] )
3342 cout <<
"TRIANGULATION_PRINT\n";
3343 cout <<
" Information defining a triangulation.\n";
3345 cout <<
" The number of points is " << point_num <<
"\n";
3350 cout <<
" The number of triangles is " << tri_num <<
"\n";
3352 cout <<
" Sets of three points are used as vertices of\n";
3353 cout <<
" the triangles. For each triangle, the points\n";
3354 cout <<
" are listed in counterclockwise order.\n";
3359 cout <<
" On each side of a given triangle, there is either\n";
3360 cout <<
" another triangle, or a piece of the convex hull.\n";
3361 cout <<
" For each triangle, we list the indices of the three\n";
3362 cout <<
" neighbors, or (if negative) the codes of the\n";
3363 cout <<
" segments of the convex hull.\n";
3370 vertex_list =
new int[3*tri_num];
3373 for ( t = 0; t < tri_num; t++ )
3375 for (
s = 0;
s < 3;
s++ )
3377 vertex_list[
k] = tri_vert[
s+t*3];
3386 delete [] vertex_list;
3390 boundary_num = 2 * vertex_num - tri_num - 2;
3393 cout <<
" The number of boundary points is " << boundary_num <<
"\n";
3395 cout <<
" The segments that make up the convex hull can be\n";
3396 cout <<
" determined from the negative entries of the triangle\n";
3397 cout <<
" neighbor list.\n";
3399 cout <<
" # Tri Side N1 N2\n";
3406 for ( i = 0; i < tri_num; i++ )
3408 for ( j = 0; j < 3; j++ )
3410 if ( tri_nabe[j+i*3] < 0 )
3412 s = -tri_nabe[j+i*3];
3415 if ( t < 1 || tri_num < t )
3418 cout <<
" Sorry, this data does not use the DTRIS2\n";
3419 cout <<
" convention for convex hull segments.\n";
3425 s2 =
i_wrap ( s1+1, 1, 3 );
3427 n1 = tri_vert[s1-1+(t-1)*3];
3428 n2 = tri_vert[s2-1+(t-1)*3];
3429 cout <<
setw(4) <<
k <<
" " 3430 <<
setw(4) << t <<
" " 3431 <<
setw(4) << s1 <<
" " 3432 <<
setw(4) << n1 <<
" " 3433 <<
setw(4) << n2 <<
"\n";
3450 void vbedg (
double x,
double y,
int point_num,
double point_xy[],
int tri_num,
3451 int tri_vert[],
int tri_nabe[],
int *ltri,
int *ledg,
int *rtri,
int *redg )
3544 l = -tri_nabe[3*((*rtri)-1)+(*redg)-1];
3547 a = tri_vert[3*(t-1)+
e-1];
3551 b = tri_vert[3*(t-1)+
e];
3555 b = tri_vert[3*(t-1)+0];
3558 ax = point_xy[2*(a-1)+0];
3559 ay = point_xy[2*(a-1)+1];
3561 bx = point_xy[2*(
b-1)+0];
3562 by = point_xy[2*(
b-1)+1];
3564 lr =
lrline (
x,
y, ax, ay, bx, by, 0.0 );
3586 b = tri_vert[3*(t-1)+
e-1];
3589 while ( 0 < tri_nabe[3*(t-1)+
e-1] )
3591 t = tri_nabe[3*(t-1)+
e-1];
3593 if ( tri_vert[3*(t-1)+0] ==
b )
3597 else if ( tri_vert[3*(t-1)+1] ==
b )
3608 a = tri_vert[3*(t-1)+
e-1];
3609 ax = point_xy[2*(a-1)+0];
3610 ay = point_xy[2*(a-1)+1];
3612 bx = point_xy[2*(
b-1)+0];
3613 by = point_xy[2*(
b-1)+1];
3615 lr =
lrline (
x,
y, ax, ay, bx, by, 0.0 );
int s_len_trim(const char *s)
void imat_transpose_print(int m, int n, int a[], const char *title)
int i_wrap(int ival, int ilo, int ihi)
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)
bool dvec_lt(int n, double a1[], double a2[])
int * d2vec_sort_heap_index_a(int n, double a[])
errorManipArg< error, int > exit(error &err, const int errNo=1)
void ivec_sorted_unique(int n, int a[], int *nuniq)
dimensionedScalar yn(const int n, const dimensionedScalar &ds)
void d2vec_sort_quick_a(int n, double a[])
dimensionedScalar sqrt(const dimensionedScalar &ds)
int * ivec_indicator(int n)
double * triangle_circumcenter_2d(double t[])
dimensionedScalar y0(const dimensionedScalar &ds)
label k
Boltzmann constant.
void dmat_transpose_print(int m, int n, double a[], const char *title)
void triangulation_print(int point_num, double xc[], int tri_num, int tri_vert[], int tri_nabe[])
bool dvec_gt(int n, double a1[], double a2[])
bool triangulation_plot_eps(const char *file_out_name, int g_num, double g_xy[], int tri_num, int nod_tri[])
void perm_inv(int n, int p[])
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
void dmat_transpose_print_some(int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, const char *title)
const dimensionedScalar e
Elementary charge.
double d_max(double x, double y)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
int lrline(double xu, double yu, double xv1, double yv1, double xv2, double yv2, double dv)
void dmat_uniform(int m, int n, double b, double c, int *seed, double r[])
void dvec_print(int n, double a[], const char *title)
int i_min(int i1, int i2)
dimensionedScalar y1(const dimensionedScalar &ds)
int * points_delaunay_naive_2d(int n, double p[], int *ntri)
const uniformDimensionedVectorField & g
void imat_transpose_print_some(int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, const char *title)
void ivec_sort_heap_a(int n, int a[])
bool perm_check(int n, int p[])
int dtris2(int point_num, double point_xy[], int *tri_num, int tri_vert[], int tri_nabe[])
void d2vec_part_quick_a(int n, double a[], int *l, int *r)
void ivec_heap_d(int n, int a[])
void d2vec_permute(int n, double a[], int p[])
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
int diaedg(double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3)
const dimensionedScalar c
Speed of light in a vacuum.
Omanip< int > setw(const int i)
double d_min(double x, double y)
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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 dvec_swap(int n, double a1[], double a2[])
bool dvec_eq(int n, double a1[], double a2[])
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[])
int i_max(int i1, int i2)