43 static const int DIM = 4;
50 void check_normal_direction(
double *normal,
double *
out);
54 double *get_centroid();
65 normal_calculated = 0;
66 centroid_calculated = 0;
67 centroid =
new double[DIM];
68 normal =
new double[DIM];
96 for (
int i=0;
i < DIM;
i++) {
97 dot_product += out[
i]*normal[
i];
101 if ( dot_product < 0 ) {
102 for (
int i=0;
i < DIM;
i++) {
103 normal[
i] = -normal[
i];
117 if ( !centroid_calculated )
118 calculate_centroid();
131 if ( !normal_calculated )
146 static const int LINE_CORNERS = 2;
147 static const int LINE_DIM = 2;
154 void calculate_centroid();
155 void calculate_normal();
159 void init(
double**,
double*,
int*);
160 void flip_start_end();
175 corners =
new double*[LINE_CORNERS];
176 for (
int i=0;
i < LINE_CORNERS;
i++) {
177 corners[
i] =
new double[DIM];
179 out =
new double[DIM];
180 const_i =
new int[DIM-LINE_DIM];
190 for (
int i=0;
i < LINE_CORNERS;
i++) {
211 for (
int i=0;
i < LINE_CORNERS;
i++) {
212 for (
int j=0;
j < DIM;
j++) {
213 corners[
i][
j] = p[
i][
j];
224 for (
int i=0;
i < DIM-LINE_DIM;
i++) {
232 }
else if ( c[1] == 2 ) {
237 }
else if ( c[0] == 1 ) {
249 normal_calculated = 0;
250 centroid_calculated = 0;
262 for (
int i=0;
i < DIM;
i++) {
263 centroid[
i] = (corners[0][
i] + corners[1][
i])/2.0;
265 centroid_calculated = 1;
276 if ( !centroid_calculated )
277 calculate_centroid();
279 normal[x1] = -(corners[1][x2] - corners[0][x2]);
280 normal[x2] = corners[1][x1] - corners[0][x1];
281 normal[const_i[0]] = 0;
282 normal[const_i[1]] = 0;
284 double *Vout =
new double[DIM];
285 for (
int j=0;
j < DIM;
j++) {
286 Vout[
j] =
out[
j] - centroid[
j];
288 check_normal_direction(normal,Vout);
290 normal_calculated = 1;
300 int temp = start_point;
301 start_point = end_point;
314 return corners[start_point];
326 return corners[end_point];
350 static const int MAX_LINES = 24;
351 static const int POLYGON_DIM = 3;
356 void calculate_centroid();
357 void calculate_normal();
362 bool add_line(
Line*,
int);
405 for (
int i=0;
i < DIM;
i++) {
409 }
else if ( x2 < 0 ) {
420 normal_calculated = 0;
421 centroid_calculated = 0;
441 if ( donotcheck || Nlines == 0 ) {
450 double *p3 =
lines[Nlines-1]->get_end();
453 for (
int j=0;
j < DIM;
j++) {
454 diff1 += fabs(p1[
j]-p3[
j]);
455 diff2 += fabs(p2[j]-p3[j]);
459 if ( diff1 < eps || diff2 < eps ) {
481 double *
mean =
new double[DIM];
482 for (
int i=0;
i < DIM;
i++ ) {
487 for (
int i=0;
i < Nlines;
i++ ) {
488 double *p1 =
lines[
i]->get_start();
489 double *p2 =
lines[
i]->get_end();
490 for (
int j=0;
j < DIM;
j++ ) {
491 mean[
j] += p1[
j] + p2[
j];
494 for (
int j=0;
j < DIM;
j++ ) {
495 mean[
j] = mean[
j]/
double(2.0*Nlines);
500 for (
int i=0;
i < DIM;
i++) {
501 centroid[
i] = mean[
i];
503 centroid_calculated = 1;
510 double *sum_up =
new double[DIM];
512 for (
int i=0;
i < DIM;
i++) {
516 double *
a =
new double[DIM];
517 double *
b =
new double[DIM];
519 double *cm_i =
new double[DIM];
520 for (
int i=0;
i < Nlines;
i++) {
521 double *p1 =
lines[
i]->get_start();
522 double *p2 =
lines[
i]->get_end();
523 for (
int j=0;
j < DIM;
j++) {
524 cm_i[
j] = (p1[
j] + p2[
j] + mean[
j])/3.0;
527 for (
int j=0;
j < DIM;
j++) {
528 a[
j] = p1[
j] - mean[
j];
529 b[
j] = p2[
j] - mean[
j];
532 double A_i = 0.5*sqrt( pow(a[x2]*b[x3]-a[x3]*b[x2],2.0) + pow(a[x1]*b[x3]-a[x3]*b[x1],2.0) + pow(a[x2]*b[x1]-a[x1]*b[x2],2.0) );
534 for (
int j=0;
j < DIM;
j++) {
535 sum_up[
j] += cm_i[
j]*A_i;
541 for (
int i=0;
i < DIM;
i++) {
542 centroid[
i] = sum_up[
i]/sum_down;
545 centroid_calculated = 1;
562 if ( !centroid_calculated )
563 calculate_centroid();
566 double **normals =
new double*[Nlines];
567 double *Vout =
new double[DIM];
568 for (
int i=0;
i < Nlines;
i++) {
569 normals[
i] =
new double[DIM];
570 for (
int j=0;
j < DIM;
j++) {
575 double *
a =
new double[DIM];
576 double *
b =
new double[DIM];
578 for (
int i=0;
i < Nlines;
i++) {
580 double *p1 =
lines[
i]->get_start();
581 double *p2 =
lines[
i]->get_end();
582 for (
int j=0;
j < DIM;
j++) {
583 a[
j] = p1[
j] - centroid[
j];
584 b[
j] = p2[
j] - centroid[
j];
587 normals[
i][x1] = 0.5*(a[x2]*b[x3]-a[x3]*b[x2]);
588 normals[
i][x2] = -0.5*(a[x1]*b[x3]-a[x3]*b[x1]);
589 normals[
i][x3] = 0.5*(a[x1]*b[x2]-a[x2]*b[x1]);
590 normals[
i][const_i] = 0;
592 double *o =
lines[
i]->get_out();
593 for (
int j=0;
j < DIM;
j++) {
594 Vout[
j] = o[
j] - centroid[
j];
597 check_normal_direction(normals[
i],Vout);
600 for (
int i=0;
i < DIM;
i++) {
603 for (
int i=0;
i < Nlines;
i++) {
604 for (
int j=0;
j < DIM;
j++) {
605 normal[
j] += normals[
i][
j];
609 normal_calculated = 1;
612 for (
int i=0;
i < Nlines;
i++) {
654 for (
int i=0;
i < Nlines;
i++) {
655 double *p1 =
lines[
i]->get_start();
656 double *p2 =
lines[
i]->get_end();
657 file << pos[x1] + p1[x1] <<
" " << pos[x2] + p1[x2] <<
" " << pos[x3] + p1[x3];
658 file <<
" " << pos[x1] + p2[x1] <<
" " << pos[x2] + p2[x2] <<
" " << pos[x3] + p2[x3];
659 file <<
" " << pos[x1] + centroid[x1] <<
" " << pos[x2] + centroid[x2] <<
" " << pos[x3] + centroid[x3] << endl;
672 static const int MAX_POLYGONS = 24;
678 void tetravolume(
double*,
double*,
double*,
double*);
679 void calculate_centroid();
680 void calculate_normal();
685 bool add_polygon(
Polygon*,
int);
697 polygons =
new Polygon*[MAX_POLYGONS];
727 normal_calculated = 0;
728 centroid_calculated = 0;
747 if ( donocheck || Npolygons == 0 ) {
748 polygons[Npolygons] =
p;
755 for (
int i=0;
i < Npolygons;
i++) {
757 int Nlines2 = polygons[
i]->get_Nlines();
759 Line **lines2 = polygons[
i]->get_lines();
760 for (
int j=0;
j < Nlines1;
j++) {
761 for (
int k=0;
k < Nlines2;
k++) {
762 if ( lines_equal(lines1[
j],lines2[
k]) ) {
763 polygons[Npolygons] =
p;
794 for (
int i=0;
i < DIM;
i++) {
795 diff1 += fabs(p11[
i]-p21[
i]);
796 diff2 += fabs(p12[i]-p21[i]);
797 if ( diff1 > eps && diff2 > eps )
817 double bc01 = b[0]*c[1]-b[1]*c[0];
818 double bc02 = b[0]*c[2]-b[2]*c[0];
819 double bc03 = b[0]*c[3]-b[3]*c[0];
820 double bc12 = b[1]*c[2]-b[2]*c[1];
821 double bc13 = b[1]*c[3]-b[3]*c[1];
822 double bc23 = b[2]*c[3]-b[3]*c[2];
823 n[0] = 1.0/6.0*(a[1]*bc23 - a[2]*bc13 + a[3]*bc12);
824 n[1] = -1.0/6.0*(a[0]*bc23 - a[2]*bc03 + a[3]*bc02);
825 n[2] = 1.0/6.0*(a[0]*bc13 - a[1]*bc03 + a[3]*bc01);
826 n[3] = -1.0/6.0*(a[0]*bc12 - a[1]*bc02 + a[2]*bc01);
844 double *
mean =
new double[DIM];
845 for (
int i=0;
i < DIM;
i++ ) {
849 for (
int i=0;
i < Npolygons;
i++ ) {
850 int Nlines = polygons[
i]->get_Nlines();
852 for (
int j=0;
j < Nlines;
j++ ) {
855 for (
int k=0;
k < DIM;
k++ ) {
856 mean[
k] += p1[
k] + p2[
k];
860 for (
int j=0;
j < DIM;
j++ ) {
861 mean[
j] = mean[
j]/
double(2.0*Ntetrahedra);
864 double *
a =
new double[DIM];
865 double *
b =
new double[DIM];
866 double *
c =
new double[DIM];
867 double *
n =
new double[DIM];
868 double *cm_i =
new double[DIM];
869 double *sum_up =
new double[DIM];
871 for (
int i=0;
i < DIM;
i++) {
874 for (
int i=0;
i < Npolygons;
i++ ) {
875 int Nlines = polygons[
i]->get_Nlines();
877 double *cent = polygons[
i]->get_centroid();
878 for (
int j=0;
j < Nlines;
j++) {
882 for (
int k=0;
k < DIM;
k++) {
883 cm_i[
k] = (p1[
k] + p2[
k] + cent[
k] + mean[
k])/4.0;
887 for (
int k=0;
k < DIM;
k++) {
888 a[
k] = p1[
k] - mean[
k];
889 b[
k] = p2[
k] - mean[
k];
890 c[
k] = cent[
k] - mean[
k];
893 tetravolume(a,b,c,n);
895 for (
int k=0;
k < DIM;
k++) {
900 for (
int k=0;
k < DIM;
k++) {
901 sum_up[
k] += cm_i[
k]*V_i;
908 for (
int i=0;
i < DIM;
i++) {
909 centroid[
i] = sum_up[
i]/sum_down;
912 centroid_calculated = 1;
932 if ( !centroid_calculated )
933 calculate_centroid();
935 double *Vout =
new double[DIM];
936 double *
a =
new double[DIM];
937 double *
b =
new double[DIM];
938 double *
c =
new double[DIM];
939 double **normals =
new double*[Ntetrahedra];
940 for (
int i=0;
i < Ntetrahedra;
i++) {
941 normals[
i] =
new double[DIM];
944 for (
int i=0;
i < Npolygons;
i++ ) {
945 int Nlines = polygons[
i]->get_Nlines();
947 double *cent = polygons[
i]->get_centroid();
948 for (
int j=0;
j < Nlines;
j++) {
952 for (
int k=0;
k < DIM;
k++) {
953 a[
k] = p1[
k] - centroid[
k];
954 b[
k] = p2[
k] - centroid[
k];
955 c[
k] = cent[
k] - centroid[
k];
958 tetravolume(a,b,c,normals[Ntetra]);
961 for (
int k=0;
k < DIM;
k++) {
962 Vout[
k] = o[
k] - centroid[
k];
964 check_normal_direction(normals[Ntetra],Vout);
969 for (
int i=0;
i < DIM;
i++) {
972 for (
int i=0;
i < Ntetrahedra;
i++) {
973 for (
int j=0;
j < DIM;
j++) {
974 normal[
j] += normals[
i][
j];
978 normal_calculated = 1;
979 for (
int i=0;
i < Ntetrahedra;
i++) {
999 static const int DIM = 4;
1000 static const int SQUARE_DIM = 2;
1001 static const int MAX_POINTS = 4;
1002 static const int MAX_LINES = 2;
1016 void ends_of_edge(
double);
1017 void find_outside(
double);
1021 void init(
double**,
int*,
double*,
double*);
1022 void construct_lines(
double);
1036 cuts =
new double*[MAX_POINTS];
1037 out =
new double*[MAX_POINTS];
1038 points =
new double*[SQUARE_DIM];
1039 const_i =
new int[DIM-SQUARE_DIM];
1040 const_value =
new double[DIM-SQUARE_DIM];
1042 for (
int i=0;
i < SQUARE_DIM;
i++) {
1043 points[
i] =
new double[SQUARE_DIM];
1045 for (
int i=0;
i < MAX_POINTS;
i++) {
1046 cuts[
i] =
new double[SQUARE_DIM];
1047 out[
i] =
new double[SQUARE_DIM];
1050 points_temp =
new double*[SQUARE_DIM];
1051 out_temp =
new double[DIM];
1052 for (
int j=0;
j < Npoints;
j++) {
1053 points_temp[
j] =
new double[DIM];
1065 for (
int i=0;
i < DIM-SQUARE_DIM;
i++) {
1069 for (
int i=0;
i < MAX_POINTS;
i++) {
1076 delete[] const_value;
1078 for (
int i=0;
i < SQUARE_DIM;
i++) {
1079 delete[] points_temp[
i];
1081 delete[] points_temp;
1099 for (
int i=0;
i < SQUARE_DIM;
i++) {
1100 for (
int j=0;
j < SQUARE_DIM;
j++) {
1101 points[
i][
j] = sq[
i][
j];
1104 for (
int i=0;
i < DIM-SQUARE_DIM;
i++) {
1105 const_i[
i] = c_i[
i];
1106 const_value[
i] = c_v[
i];
1111 for (
int i=0;
i < DIM;
i++) {
1112 if (
i != const_i[0] &&
i != const_i[1] ) {
1138 for (
int i=0;
i < DIM-SQUARE_DIM;
i++) {
1139 for (
int j=0;
j < DIM-SQUARE_DIM;
j++) {
1140 if ( points[
i][
j] >= E0 )
1146 if ( above == 0 || above == 4 ) {
1157 for (
int i=0;
i < Ncuts;
i++) {
1158 points_temp[
i%2][x1] =
cuts[
i][0];
1159 points_temp[
i%2][x2] =
cuts[
i][1];
1160 points_temp[
i%2][const_i[0]] = const_value[0];
1161 points_temp[
i%2][const_i[1]] = const_value[1];
1165 out_temp[x1] =
out[
i/2][0];
1166 out_temp[x2] =
out[
i/2][1];
1167 out_temp[const_i[0]] = const_value[0];
1168 out_temp[const_i[1]] = const_value[1];
1169 lines[
i/2].init(points_temp,out_temp,const_i);
1195 if ( (points[0][0]-E0)*(points[1][0]-E0) < 0 ) {
1196 cuts[Ncuts][0] = (points[0][0]-E0)/(points[0][0]-points[1][0])*dx[x1];
1199 }
else if ( points[0][0] == E0 && points[1][0] < E0 ) {
1200 cuts[Ncuts][0] = 1
e-9*dx[x1];
1203 }
else if ( points[1][0] == E0 && points[0][0] < E0 ) {
1204 cuts[Ncuts][0] = (1.0-1
e-9)*dx[x1];
1209 if ( (points[0][0]-E0)*(points[0][1]-E0) < 0 ) {
1211 cuts[Ncuts][1] = (points[0][0]-E0)/(points[0][0]-points[0][1])*dx[x2];
1213 }
else if ( points[0][0] == E0 && points[0][1] < E0 ) {
1215 cuts[Ncuts][1] = 1
e-9*dx[x2];
1217 }
else if ( points[0][1] == E0 && points[0][0] < E0 ) {
1219 cuts[Ncuts][1] = (1.0-1
e-9)*dx[x2];
1223 if ( (points[1][0]-E0)*(points[1][1]-E0) < 0 ) {
1224 cuts[Ncuts][0] = dx[x1];
1225 cuts[Ncuts][1] = (points[1][0]-E0)/(points[1][0]-points[1][1])*dx[x2];
1227 }
else if ( points[1][0] == E0 && points[1][1] < E0 ) {
1228 cuts[Ncuts][0] = dx[x1];
1229 cuts[Ncuts][1] = 1
e-9*dx[x2];
1231 }
else if ( points[1][1] == E0 && points[1][0] < E0 ) {
1232 cuts[Ncuts][0] = dx[x1];
1233 cuts[Ncuts][1] = (1.0-1
e-9)*dx[x2];
1237 if ( (points[0][1]-E0)*(points[1][1]-E0) < 0 ) {
1238 cuts[Ncuts][0] = (points[0][1]-E0)/(points[0][1]-points[1][1])*dx[x1];
1239 cuts[Ncuts][1] = dx[x2];
1241 }
else if ( points[0][1] == E0 && points[1][1] < E0 ) {
1242 cuts[Ncuts][0] = 1
e-9*dx[x1];
1243 cuts[Ncuts][1] = dx[x2];
1245 }
else if ( points[1][1] == E0 && points[0][1] < E0 ) {
1246 cuts[Ncuts][0] = (1.0-1
e-9)*dx[x1];
1247 cuts[Ncuts][1] = dx[x2];
1250 if ( Ncuts != 0 && Ncuts != 2 && Ncuts != 4 ) {
1251 cout <<
"Error in EndsOfEdge: Ncuts " << Ncuts << endl;
1269 for (
int i=0;
i < 2;
i++) {
1270 for (
int j=0;
j < 2;
j++) {
1271 Emid += 0.25*points[
i][
j];
1277 if ( (points[0][0] < E0 && Emid < E0) || (points[0][0] > E0 && Emid > E0) ) {
1278 for (
int i=0;
i < 2;
i++) {
1279 double temp =
cuts[1][
i];
1286 if ( (Emid-E0) < 0 ) {
1287 for (
int i=0;
i < 2;
i++) {
1288 for (
int j=0;
j < 2;
j++) {
1290 out[
i][
j] = 0.5*dx[x1];
1292 out[
i][
j] = 0.5*dx[x2];
1297 if ( (points[0][0]-E0) < 0 ) {
1298 for (
int i=0;
i < 2;
i++) {
1314 for (
int i=0;
i < 2;
i++) {
1315 for (
int j=0;
j < 2;
j++) {
1320 for (
int i=0;
i < 2;
i++) {
1321 for (
int j=0;
j < 2;
j++) {
1322 if ( points[
i][
j] < E0 ) {
1323 out[0][0] +=
i*dx[x1];
1324 out[0][1] +=
j*dx[x2];
1330 for (
int i=0;
i < 2;
i++) {
1331 for (
int j=0;
j < 2;
j++) {
1385 static const int DIM = 4;
1386 static const int CUBE_DIM = 4;
1387 static const int MAX_POLY = 8;
1388 static const int NSQUARES = 6;
1389 static const int STEPS = 2;
1401 void split_to_squares();
1402 void check_ambiguous(
int);
1406 void init(
double***&,
int,
double,
double*&);
1407 void construct_polygons(
double);
1409 int get_Npolygons();
1421 cube =
new double**[STEPS];
1422 for (
int i=0;
i < STEPS;
i++) {
1423 cube[
i] =
new double*[STEPS];
1424 for (
int j=0;
j < STEPS;
j++) {
1425 cube[
i][
j] =
new double[STEPS];
1429 polygons =
new Polygon[MAX_POLY];
1430 squares =
new Square[NSQUARES];
1440 for (
int i=0;
i < STEPS;
i++) {
1441 for (
int j=0;
j < STEPS;
j++) {
1442 delete[] cube[
i][
j];
1472 for (
int i=0;
i < DIM;
i++) {
1476 }
else if ( x2 < 0 ) {
1483 for (
int i=0;
i < STEPS;
i++) {
1484 for (
int j=0;
j < STEPS;
j++) {
1485 for (
int k=0;
k < STEPS;
k++) {
1486 cube[
i][
j][
k] = c[
i][
j][
k];
1504 double **sq =
new double*[STEPS];
1505 int *c_i =
new int[STEPS];
1506 double *c_v =
new double[STEPS];
1507 for (
int i=0;
i < STEPS;
i++) {
1508 sq[
i] =
new double[STEPS];
1511 for (
int i=0;
i < DIM;
i++) {
1514 if (
i != const_i ) {
1517 for (
int j=0;
j < STEPS;
j++) {
1518 c_v[0] = const_value;
1520 for (
int ci1=0; ci1 < STEPS; ci1++) {
1521 for (
int ci2=0; ci2 < STEPS; ci2++) {
1523 sq[ci1][ci2] = cube[
j][ci1][ci2];
1524 }
else if (
i == x2 ) {
1525 sq[ci1][ci2] = cube[ci1][
j][ci2];
1527 sq[ci1][ci2] = cube[ci1][ci2][
j];
1531 squares[Nsquares].init(sq,c_i,c_v,dx);
1536 for (
int i=0;
i < STEPS;
i++) {
1557 for (
int i=0;
i < NSQUARES;
i++) {
1558 squares[
i].construct_lines(value0);
1562 for (
int i=0;
i < NSQUARES;
i++) {
1563 int Nline = squares[
i].get_Nlines();
1564 Line *l = squares[
i].get_lines();
1565 for (
int j=0;
j < Nline;
j++) {
1567 double *
p =
lines[Nlines]->get_start();
1568 p =
lines[Nlines]->get_end();
1574 if ( Nlines == 0 ) {
1578 check_ambiguous(Nlines);
1579 if ( ambiguous > 0 ) {
1582 int *not_used =
new int[Nlines];
1583 for (
int i=0;
i < Nlines;
i++) {
1589 if ( Nlines - used < 3 ) {
1590 cout <<
"Error: cannot construct polygon from " << Nlines -used <<
" lines" << endl;
1594 polygons[Npolygons].init(const_i);
1596 for (
int i=0;
i < Nlines;
i++) {
1597 if ( not_used[
i] ) {
1599 if ( polygons[Npolygons].add_line(
lines[i],0) ) {
1610 }
while ( used < Nlines );
1615 polygons[Npolygons].init(const_i);
1616 for (
int i=0;
i < Nlines;
i++) {
1617 polygons[Npolygons].add_line(
lines[
i],1);
1631 for (
int i=0;
i < NSQUARES;
i++) {
1632 if ( squares[
i].is_ambiguous() ) {
1639 if ( ambiguous == 0 && Nlines == 6 ) {
1702 static const int DIM = 4;
1703 static const int MAX_POLY = 10;
1704 static const int NCUBES = 8;
1705 static const int STEPS = 2;
1714 void split_to_cubes();
1715 void check_ambiguous(
double);
1719 void init(
double****,
double*);
1720 void construct_polyhedrons(
double);
1721 int get_Npolyhedrons();
1732 hcube =
new double***[STEPS];
1733 for (
int i=0;
i < STEPS;
i++) {
1734 hcube[
i] =
new double**[STEPS];
1735 for (
int j=0;
j < STEPS;
j++) {
1736 hcube[
i][
j] =
new double*[STEPS];
1737 for (
int k=0;
k < STEPS;
k++) {
1738 hcube[
i][
j][
k] =
new double[STEPS];
1743 polygons =
new Polygon*[NCUBES*10];
1744 cubes =
new Cube[NCUBES];
1754 for (
int i=0;
i < STEPS;
i++) {
1755 for (
int j=0;
j < STEPS;
j++) {
1756 for (
int k=0;
k < STEPS;
k++) {
1757 delete[] hcube[
i][
j][
k];
1759 delete[] hcube[
i][
j];
1764 delete[] polyhedrons;
1786 for (
int i=0;
i < STEPS;
i++) {
1787 for (
int j=0;
j < STEPS;
j++) {
1788 for (
int k=0;
k < STEPS;
k++) {
1789 for (
int l=0; l < STEPS; l++) {
1790 hcube[
i][
j][
k][l] = c[
i][
j][
k][l];
1808 double ***cu =
new double**[STEPS];
1809 for (
int i=0;
i < STEPS;
i++) {
1810 cu[
i] =
new double*[STEPS];
1811 for (
int j=0;
j < STEPS;
j++) {
1812 cu[
i][
j] =
new double[STEPS];
1816 for (
int i=0;
i < DIM;
i++) {
1820 for (
int j=0;
j < STEPS;
j++) {
1821 double c_v =
j*dx[
i];
1822 for (
int ci1=0; ci1 < STEPS; ci1++) {
1823 for (
int ci2=0; ci2 < STEPS; ci2++) {
1824 for (
int ci3=0; ci3 < STEPS; ci3++) {
1826 cu[ci1][ci2][ci3] = hcube[
j][ci1][ci2][ci3];
1827 }
else if (
i == x2 ) {
1828 cu[ci1][ci2][ci3] = hcube[ci1][
j][ci2][ci3];
1829 }
else if (
i == x3 ) {
1830 cu[ci1][ci2][ci3] = hcube[ci1][ci2][
j][ci3];
1832 cu[ci1][ci2][ci3] = hcube[ci1][ci2][ci3][
j];
1837 cubes[Ncubes].init(cu,c_i,c_v,dx);
1841 for (
int i=0;
i < STEPS;
i++) {
1842 for (
int j=0;
j < STEPS;
j++) {
1861 for (
int i=0;
i < NCUBES;
i++) {
1862 cubes[
i].construct_polygons(value0);
1866 for (
int i=0;
i < NCUBES;
i++) {
1867 int Npoly = cubes[
i].get_Npolygons();
1869 for (
int j=0;
j < Npoly;
j++) {
1870 polygons[Npolygons] = &p[
j];
1874 check_ambiguous(value0);
1875 if ( ambiguous > 0 ) {
1878 int *not_used =
new int[Npolygons];
1879 for (
int i=0;
i < Npolygons;
i++) {
1886 polyhedrons[Npolyhedrons].init();
1888 for (
int i=0;
i < Npolygons;
i++) {
1889 if ( not_used[
i] ) {
1891 if ( polyhedrons[Npolyhedrons].add_polygon(polygons[i],0) ) {
1902 }
while ( used < Npolygons );
1910 polyhedrons[Npolyhedrons].init();
1911 for (
int i=0;
i < Npolygons;
i++) {
1912 polyhedrons[Npolyhedrons].add_polygon(polygons[
i],1);
1926 for (
int i=0;
i < NCUBES;
i++) {
1927 if ( cubes[
i].is_ambiguous() ) {
1931 if ( ambiguous == 0 ) {
1933 for (
int i=0;
i < NCUBES;
i++) {
1934 Nlines += cubes[
i].get_Nlines();
1937 for (
int i1=0; i1 < STEPS; i1++) {
1938 for (
int i2=0; i2 < STEPS; i2++) {
1939 for (
int i3=0; i3 < STEPS; i3++) {
1940 for (
int i4=0; i4 < STEPS; i4++) {
1941 if ( hcube[i1][i2][i3][i4] < value0 ) {
1959 if ( Nlines == 24 && points == 2 ) {
1974 return Npolyhedrons;
2002 static const int STEPS = 2;
2003 static const int DIM = 4;
2004 static const int MAX_ELEMENTS = 10;
2014 void surface_3d(
double***,
double*,
int);
2021 void init(
int,
double,
double*);
2022 void init_print(
string);
2023 void find_surface_2d(
double**);
2024 void find_surface_3d(
double***);
2025 void find_surface_3d_print(
double***,
double*);
2026 void find_surface_4d(
double****);
2027 int get_Nelements();
2028 double **get_normals();
2029 double **get_centroids();
2030 double **get_normals_4d();
2031 double **get_centroids_4d();
2032 double get_centroid_elem(
int,
int);
2033 double get_normal_elem(
int,
int);
2045 print_initialized = 0;
2046 normals =
new double*[MAX_ELEMENTS];
2047 centroids =
new double*[MAX_ELEMENTS];
2048 for (
int i=0;
i < MAX_ELEMENTS;
i++) {
2049 normals[
i] =
new double[DIM];
2050 centroids[
i] =
new double[DIM];
2061 for (
int i=0;
i < MAX_ELEMENTS;
i++) {
2062 delete[] normals[
i];
2063 delete[] centroids[
i];
2067 if ( initialized ) {
2071 if ( print_initialized ) {
2072 output_print.close();
2091 if ( initialized != 1 ) {
2092 dx =
new double[DIM];
2094 for (
int i=0;
i < DIM;
i++) {
2095 if (
i < DIM-cube_dim ) {
2098 dx[
i] = dex[
i-(DIM-cube_dim)];
2115 output_print.open(filename.c_str());
2116 print_initialized = 1;
2129 if ( !initialized || cube_dim != 2 ) {
2130 cout <<
"Cornelius not initialized for 2D case" << endl;
2133 int *c_i =
new int[cube_dim];
2134 double *c_v =
new double[cube_dim];
2139 cu2d.init(cube,c_i,c_v,dx);
2140 cu2d.construct_lines(value0);
2141 Nelements = cu2d.get_Nlines();
2142 Line *l = cu2d.get_lines();
2143 for (
int i=0;
i < Nelements;
i++) {
2144 for (
int j=0;
j < DIM;
j++) {
2164 surface_3d(cube,pos,0);
2179 surface_3d(cube,pos,1);
2194 if ( !initialized || cube_dim != 3 ) {
2195 cout <<
"Cornelius not initialized for 3D case" << endl;
2202 for (
int i=0;
i < STEPS;
i++) {
2203 for (
int j=0;
j < STEPS;
j++) {
2204 for (
int k=0;
k < STEPS;
k++) {
2205 if ( cube[
i][
j][
k] >= value0 )
2210 if ( above == 0 || above == 8 ) {
2219 cu3d.init(cube,c_i,c_v,dx);
2221 cu3d.construct_polygons(value0);
2223 Nelements = cu3d.get_Npolygons();
2225 for (
int i=0;
i < Nelements;
i++) {
2227 for (
int j=0;
j < DIM;
j++) {
2233 if ( print_initialized && do_print ) {
2234 p[
i].
print(output_print,pos);
2250 if ( !initialized || cube_dim != 4 ) {
2251 cout <<
"Cornelius not initialized for 4D case" << endl;
2258 for (
int i=0;
i < STEPS;
i++) {
2259 for (
int j=0;
j < STEPS;
j++) {
2260 for (
int k=0;
k < STEPS;
k++) {
2261 for (
int l=0; l < STEPS; l++) {
2262 if ( cube[
i][
j][
k][l] >= value0 )
2268 if ( above == 0 || above == 16 ) {
2277 cu4d.construct_polyhedrons(value0);
2279 Nelements = cu4d.get_Npolyhedrons();
2281 for (
int i=0;
i < Nelements;
i++) {
2282 for (
int j=0;
j < DIM;
j++) {
2314 double **vect =
new double*[Nelements];
2315 for (
int i=0;
i < Nelements;
i++) {
2316 vect[
i] =
new double[DIM];
2317 for (
int j=0;
j < DIM;
j++) {
2318 vect[
i][
j] = centroids[
i][
j];
2338 double **vect =
new double*[Nelements];
2339 for (
int i=0;
i < Nelements;
i++) {
2340 vect[
i] =
new double[DIM];
2341 for (
int j=0;
j < DIM;
j++) {
2342 vect[
i][
j] = normals[
i][
j];
2360 double **vect =
new double*[Nelements];
2361 for (
int i=0;
i < Nelements;
i++) {
2362 vect[
i] =
new double[cube_dim];
2363 for (
int j=0;
j < cube_dim;
j++) {
2364 vect[
i][
j] = centroids[
i][
j+(DIM-cube_dim)];
2383 double **vect =
new double*[Nelements];
2384 for (
int i=0;
i < Nelements;
i++) {
2385 vect[
i] =
new double[cube_dim];
2386 for (
int j=0;
j < cube_dim;
j++) {
2387 vect[
i][
j] = normals[
i][
j+(DIM-cube_dim)];
2406 if ( i >= Nelements || j >= cube_dim ) {
2407 cout <<
"Cornelius error: asking for an element which does not exist." << endl;
2410 return centroids[
i][j+(DIM-cube_dim)];
2427 if ( i >= Nelements || j >= cube_dim ) {
2428 cout <<
"Cornelius error: asking for an element which does not exist." << endl;
2431 return normals[
i][j+(DIM-cube_dim)];