Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cornelius.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file cornelius.cpp
1 
27 #include <iostream>
28 #include <fstream>
29 #include <math.h>
30 #include <stdlib.h>
31 
32 using namespace std;
33 
41 {
42  protected:
43  static const int DIM = 4;
44  double *centroid;
45  double *normal;
48  virtual void calculate_centroid() {};
49  virtual void calculate_normal() {};
50  void check_normal_direction(double *normal, double *out);
51  public:
53  ~GeneralElement();
54  double *get_centroid();
55  double *get_normal();
56 };
57 
64 {
65  normal_calculated = 0;
66  centroid_calculated = 0;
67  centroid = new double[DIM];
68  normal = new double[DIM];
69 }
70 
77 {
78  delete[] centroid;
79  delete[] normal;
80 }
81 
91 void GeneralElement::check_normal_direction(double *normal, double *out)
92 {
93  //We calculate the dot product, if less than 4 dimensions the elements,
94  //which are not used, are just zero and do not effect this
95  double dot_product = 0;
96  for (int i=0; i < DIM; i++) {
97  dot_product += out[i]*normal[i];
98  }
99  //If the dot product is negative, the normal is pointing in the
100  //wrong direction and we have to flip it's direction
101  if ( dot_product < 0 ) {
102  for (int i=0; i < DIM; i++) {
103  normal[i] = -normal[i];
104  }
105  }
106 }
107 
116 {
117  if ( !centroid_calculated )
118  calculate_centroid();
119  return centroid;
120 }
121 
130 {
131  if ( !normal_calculated )
132  calculate_normal();
133  return normal;
134 }
135 
136 
143 class Line : public GeneralElement
144 {
145  private:
146  static const int LINE_CORNERS = 2;
147  static const int LINE_DIM = 2;
148  int x1,x2;
151  double **corners;
152  double *out;
153  int *const_i;
154  void calculate_centroid();
155  void calculate_normal();
156  public:
157  Line();
158  ~Line();
159  void init(double**,double*,int*);
160  void flip_start_end();
161  double *get_start();
162  double *get_end();
163  double *get_out();
164 };
165 
173 {
174  GeneralElement();
175  corners = new double*[LINE_CORNERS];
176  for (int i=0; i < LINE_CORNERS; i++) {
177  corners[i] = new double[DIM];
178  }
179  out = new double[DIM];
180  const_i = new int[DIM-LINE_DIM];
181 }
182 
189 {
190  for (int i=0; i < LINE_CORNERS; i++) {
191  delete[] corners[i];
192  }
193  delete[] corners;
194  delete[] out;
195  delete[] const_i;
196 }
197 
208 void Line::init(double **p, double *o, int *c)
209 {
210  //We copy the values of the end points and the outside point
211  for (int i=0; i < LINE_CORNERS; i++) {
212  for (int j=0; j < DIM; j++) {
213  corners[i][j] = p[i][j];
214  if ( i==0 ) {
215  out[j] = o[j];
216  }
217  }
218  }
219  //Let's set indices for start and end points, so we can conveniently
220  //flip the start and end point without changing the actual values
221  start_point = 0;
222  end_point = 1;
223  //Here we copy the indices which are constant at this line
224  for (int i=0; i < DIM-LINE_DIM; i++) {
225  const_i[i] = c[i];
226  }
227  //Here we fix the non-zero indices in such a way that x1 is
228  //always smaller
229  if ( c[0] == 0 ) {
230  if ( c[1] == 1 ) {
231  x1 = 2; x2 = 3;
232  } else if ( c[1] == 2 ) {
233  x1 = 1; x2 = 3;
234  } else {
235  x1 = 1; x2 = 2;
236  }
237  } else if ( c[0] == 1 ) {
238  if ( c[1] == 2 ) {
239  x1 = 0; x2 = 3;
240  } else {
241  x1 = 0; x2 = 2;
242  }
243  } else {
244  x1 = 0; x2 = 1;
245  }
246  //We need to set these to zero so when this function is called
247  //again we know that normal and centroid are not yet calculated
248  //with new values
249  normal_calculated = 0;
250  centroid_calculated = 0;
251 }
252 
260 {
261  //Centroid is just the average of the points
262  for (int i=0; i < DIM; i++) {
263  centroid[i] = (corners[0][i] + corners[1][i])/2.0;
264  }
265  centroid_calculated = 1;
266 }
267 
274 {
275  //Centroid must be calculated before we can calculate normal
276  if ( !centroid_calculated )
277  calculate_centroid();
278  //Normal is just (-dy,dx)
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;
283  //Now we check if the normal is in the correct direction
284  double *Vout = new double[DIM];
285  for (int j=0; j < DIM; j++) {
286  Vout[j] = out[j] - centroid[j];
287  }
288  check_normal_direction(normal,Vout);
289  delete[] Vout;
290  normal_calculated = 1;
291 }
292 
299 {
300  int temp = start_point;
301  start_point = end_point;
302  end_point = temp;
303 }
304 
313 {
314  return corners[start_point];
315 }
316 
324 double *Line::get_end()
325 {
326  return corners[end_point];
327 }
328 
336 double *Line::get_out()
337 {
338  return out;
339 }
340 
347 class Polygon : public GeneralElement
348 {
349  private:
350  static const int MAX_LINES = 24;
351  static const int POLYGON_DIM = 3;
353  int Nlines;
354  int x1,x2,x3;
355  int const_i;
356  void calculate_centroid();
357  void calculate_normal();
358  public:
359  Polygon();
360  ~Polygon();
361  void init(int);
362  bool add_line(Line*,int);
363  int get_Nlines();
364  Line** get_lines();
365  void print(ofstream &file,double*);
366 };
367 
375 {
376  GeneralElement();
377  lines = new Line*[MAX_LINES];
378 }
379 
386 {
387  delete[] lines;
388 }
389 
397 void Polygon::init(int c)
398 {
399  //Let's copy the constant index
400  const_i = c;
401  //Here we fix the indices which are not constant
402  x1 = -1;
403  x2 = -1;
404  x3 = -1;
405  for (int i=0; i < DIM; i++) {
406  if ( i != c ) {
407  if ( x1 < 0 ) {
408  x1 = i;
409  } else if ( x2 < 0 ) {
410  x2 = i;
411  } else {
412  x3 = i;
413  }
414  }
415  }
416  //We need to set these to zero so when this function is called
417  //again we know that polygon has no lines added and normal and
418  //centroid are not yet calculated with new values
419  Nlines = 0;
420  normal_calculated = 0;
421  centroid_calculated = 0;
422 }
423 
436 bool Polygon::add_line(Line *l, int donotcheck)
437 {
438  double eps = 1e-10;
439  //If this is the first line or we do not want to check, line is added
440  //automatically
441  if ( donotcheck || Nlines == 0 ) {
442  lines[Nlines] = l;
443  Nlines++;
444  return true;
445  } else {
446  //Here we check if the new line is connected with the last line in
447  //polygon. This since lines are ordered.
448  double *p1 = l->get_start();
449  double *p2 = l->get_end();
450  double *p3 = lines[Nlines-1]->get_end();
451  double diff1 = 0;
452  double diff2 = 0;
453  for (int j=0; j < DIM; j++) {
454  diff1 += fabs(p1[j]-p3[j]);
455  diff2 += fabs(p2[j]-p3[j]);
456  }
457  //If start or end point is the same as the end point of previous line, line
458  //is connected to polygon and it is added
459  if ( diff1 < eps || diff2 < eps ) {
460  //We always want that start point is connected to end so we flip the
461  //start and end in the new line if needed
462  if ( diff2 < eps )
463  l->flip_start_end();
464  lines[Nlines] = l;
465  Nlines++;
466  return true;
467  }
468  //If we are here, the line was not connected to the current polygon
469  return false;
470  }
471 }
472 
479 {
480  //We need a vector for the mean of the corners.
481  double *mean = new double[DIM];
482  for (int i=0; i < DIM; i++ ) {
483  mean[i] = 0;
484  }
485  //First we determine the mean of the corner points. Every point
486  //is here twice since lines are not ordered
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];
492  }
493  }
494  for (int j=0; j < DIM; j++ ) {
495  mean[j] = mean[j]/double(2.0*Nlines);
496  }
497  //If only 3 lines, i.e. 3 corners, the polygon is always on a plane
498  //and centroid is the mean
499  if ( Nlines == 3 ) {
500  for (int i=0; i < DIM; i++) {
501  centroid[i] = mean[i];
502  }
503  centroid_calculated = 1;
504  delete[] mean;
505  return;
506  }
507  //If more than 3 corners, calculation of the centroid is more
508  //complicated
509  //Here we from triangles from the lines and the mean point
510  double *sum_up = new double[DIM]; //areas of the single triangles
511  double sum_down = 0; //area of all the triangles
512  for (int i=0; i < DIM; i++) {
513  sum_up[i] = 0;
514  }
515  //a and b are vectors which from the triangle
516  double *a = new double[DIM];
517  double *b = new double[DIM];
518  //centroid of the triangle (this is always on a plane)
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;
525  }
526  //Triangle is defined by these two vectors
527  for (int j=0; j < DIM; j++) {
528  a[j] = p1[j] - mean[j];
529  b[j] = p2[j] - mean[j];
530  }
531  //Area is calculated from cross product of these vectors
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) );
533  //Then we store the area and update the total area
534  for (int j=0; j < DIM; j++) {
535  sum_up[j] += cm_i[j]*A_i;
536  }
537  sum_down += A_i;
538  }
539  //Now we can determine centroid as a weighted average of the centroids
540  //of the triangles
541  for (int i=0; i < DIM; i++) {
542  centroid[i] = sum_up[i]/sum_down;
543  }
544  //centroid is now calculated and memory is freed
545  centroid_calculated = 1;
546  delete[] sum_up;
547  delete[] mean;
548  delete[] a;
549  delete[] b;
550  delete[] cm_i;
551 }
552 
560 {
561  //centroid must be calculated before normal can be determined
562  if ( !centroid_calculated )
563  calculate_centroid();
564  //First we find the normal for each triangle formed from
565  //one edge and centroid
566  double **normals = new double*[Nlines]; //these are the normals
567  double *Vout = new double[DIM]; //the point which is always outside
568  for (int i=0; i < Nlines; i++) {
569  normals[i] = new double[DIM];
570  for (int j=0; j < DIM; j++) {
571  normals[i][j] = 0;
572  }
573  }
574  //Normal is defined by these two vectors
575  double *a = new double[DIM];
576  double *b = new double[DIM];
577  //Loop over all triangles
578  for (int i=0; i < Nlines; i++) {
579  //First we calculate the vectors which form the triangle
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];
585  }
586  //Normal is calculated as a cross product of these vectors
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;
591  //Then we construct a vector which points out
592  double *o = lines[i]->get_out();
593  for (int j=0; j < DIM; j++) {
594  Vout[j] = o[j] - centroid[j];
595  }
596  //then we check that normal is point in the correct direction
597  check_normal_direction(normals[i],Vout);
598  }
599  //Finally the normal is a sum of the normals of the triangles
600  for (int i=0; i < DIM; i++) {
601  normal[i] = 0;
602  }
603  for (int i=0; i < Nlines; i++) {
604  for (int j=0; j < DIM; j++) {
605  normal[j] += normals[i][j];
606  }
607  }
608  //Normal is now calculated and memory can be freed
609  normal_calculated = 1;
610  delete[] a;
611  delete[] b;
612  for (int i=0; i < Nlines; i++) {
613  delete[] normals[i];
614  }
615  delete[] normals;
616  delete[] Vout;
617 }
618 
627 {
628  return Nlines;
629 }
630 
639 {
640  return lines;
641 }
642 
652 void Polygon::print(ofstream &file, double *pos)
653 {
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;
660  }
661 }
662 
670 {
671  private:
672  static const int MAX_POLYGONS = 24;
676  int x1,x2,x3,x4;
677  bool lines_equal(Line*,Line*);
678  void tetravolume(double*,double*,double*,double*);
679  void calculate_centroid();
680  void calculate_normal();
681  public:
682  Polyhedron();
683  ~Polyhedron();
684  void init();
685  bool add_polygon(Polygon*,int);
686 };
687 
695 {
696  GeneralElement();
697  polygons = new Polygon*[MAX_POLYGONS];
698 }
699 
706 {
707  delete[] polygons;
708 }
709 
716 {
717  //Here we fix the non-constant indices
718  x1 = 0;
719  x2 = 1;
720  x3 = 2;
721  x4 = 3;
722  //We need to set these to zero so when this function is called
723  //again we know that polyhedron has no polyfons added and normal
724  //and centroid are not yet calculated with new values
725  Npolygons = 0;
726  Ntetrahedra = 0;
727  normal_calculated = 0;
728  centroid_calculated = 0;
729 }
730 
743 bool Polyhedron::add_polygon(Polygon *p, int donocheck)
744 {
745  //If this is the first polygon or we want to add it in any case, we
746  //just add it automatically
747  if ( donocheck || Npolygons == 0 ) {
748  polygons[Npolygons] = p;
749  Npolygons++;
750  Ntetrahedra += p->get_Nlines();
751  return true;
752  } else {
753  //Here we check if the new polygon is connected with the last polygon in
754  //polyhedron
755  for (int i=0; i < Npolygons; i++) {
756  int Nlines1 = p->get_Nlines();
757  int Nlines2 = polygons[i]->get_Nlines();
758  Line **lines1 = p->get_lines();
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;
764  Npolygons++;
765  Ntetrahedra += p->get_Nlines();
766  return true;
767  }
768  }
769  }
770  }
771  //If we are here, the polygon was not connected to polyhedron
772  return false;
773  }
774 }
775 
787 {
788  double eps = 1e-10;
789  double *p11 = l1->get_start();
790  double *p12 = l1->get_end();
791  double *p21 = l2->get_start();
792  double diff1 = 0;
793  double diff2 = 0;
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 )
798  return false;
799  }
800  //If we are here, the lines are connected
801  return true;
802 }
803 
815 void Polyhedron::tetravolume(double *a, double *b, double *c, double *n)
816 {
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);
827  /*n[0] = 1.0/6.0*(a[1]*(b[2]*c[3]-b[3]*c[2]) - a[2]*(b[1]*c[3]-b[3]*c[1])
828  +a[3]*(b[1]*c[2]-b[2]*c[1]));
829  n[1] = -1.0/6.0*(a[0]*(b[2]*c[3]-b[3]*c[2]) - a[2]*(b[0]*c[3]-b[3]*c[0])
830  +a[3]*(b[0]*c[2]-b[2]*c[0]));
831  n[2] = 1.0/6.0*(a[0]*(b[1]*c[3]-b[3]*c[1]) - a[1]*(b[0]*c[3]-b[3]*c[0])
832  +a[3]*(b[0]*c[1]-b[1]*c[0]));
833  n[3] = -1.0/6.0*(a[0]*(b[1]*c[2]-b[2]*c[1]) - a[1]*(b[0]*c[2]-b[2]*c[0])
834  +a[2]*(b[0]*c[1]-b[1]*c[0]));*/
835 }
836 
843 {
844  double *mean = new double[DIM];
845  for (int i=0; i < DIM; i++ ) {
846  mean[i] = 0;
847  }
848  //First we determine the mean of the corners
849  for (int i=0; i < Npolygons; i++ ) {
850  int Nlines = polygons[i]->get_Nlines();
851  Line **lines = polygons[i]->get_lines();
852  for (int j=0; j < Nlines; j++ ) {
853  double *p1 = lines[j]->get_start();
854  double *p2 = lines[j]->get_end();
855  for (int k=0; k < DIM; k++ ) {
856  mean[k] += p1[k] + p2[k];
857  }
858  }
859  }
860  for (int j=0; j < DIM; j++ ) {
861  mean[j] = mean[j]/double(2.0*Ntetrahedra);
862  }
863  //Some memory allocated for temporary variables
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];
870  double sum_down = 0;
871  for (int i=0; i < DIM; i++) {
872  sum_up[i] = 0;
873  }
874  for (int i=0; i < Npolygons; i++ ) {
875  int Nlines = polygons[i]->get_Nlines();
876  Line **lines = polygons[i]->get_lines();
877  double *cent = polygons[i]->get_centroid();
878  for (int j=0; j < Nlines; j++) {
879  double *p1 = lines[j]->get_start();
880  double *p2 = lines[j]->get_end();
881  //CM of the tetrahedra
882  for (int k=0; k < DIM; k++) {
883  cm_i[k] = (p1[k] + p2[k] + cent[k] + mean[k])/4.0;
884  }
885  //Volume of a tetrahedron is determined next
886  //Tetrahedron is defined by these three vectors
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];
891  }
892  //Then we calculate the volume from the normal n
893  tetravolume(a,b,c,n);
894  double V_i = 0;
895  for (int k=0; k < DIM; k++) {
896  V_i += n[k]*n[k];
897  }
898  V_i = sqrt(V_i);
899  //Then we add contributions to sums
900  for (int k=0; k < DIM; k++) {
901  sum_up[k] += cm_i[k]*V_i;
902  }
903  sum_down += V_i;
904  }
905  }
906  //Now the centroid is the volume weighted average of individual
907  //tetrahedra
908  for (int i=0; i < DIM; i++) {
909  centroid[i] = sum_up[i]/sum_down;
910  }
911  //Centroid is now calculated and we can free memory
912  centroid_calculated = 1;
913  delete[] a;
914  delete[] b;
915  delete[] c;
916  delete[] n;
917  delete[] cm_i;
918  delete[] sum_up;
919  delete[] mean;
920 }
921 
929 {
930  //We need centroid in the following calculation and thus we
931  //need to check that it is calculated
932  if ( !centroid_calculated )
933  calculate_centroid();
934  //First we allocate memory for temporary variables
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];
942  }
943  int Ntetra = 0; //Index of tetrahedron we are dealing with
944  for (int i=0; i < Npolygons; i++ ) {
945  int Nlines = polygons[i]->get_Nlines();
946  Line **lines = polygons[i]->get_lines();
947  double *cent = polygons[i]->get_centroid();
948  for (int j=0; j < Nlines; j++) {
949  double *p1 = lines[j]->get_start();
950  double *p2 = lines[j]->get_end();
951  //Tetrahedron is defined by these three vectors
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];
956  }
957  //Normal is calculated with the same function as volume
958  tetravolume(a,b,c,normals[Ntetra]);
959  //Then we determine the direction towards lower energy
960  double *o = lines[j]->get_out();
961  for (int k=0; k < DIM; k++) {
962  Vout[k] = o[k] - centroid[k];
963  }
964  check_normal_direction(normals[Ntetra],Vout);
965  Ntetra++;
966  }
967  }
968  //Finally the element normal is a sum of the calculated normals
969  for (int i=0; i < DIM; i++) {
970  normal[i] = 0;
971  }
972  for (int i=0; i < Ntetrahedra; i++) {
973  for (int j=0; j < DIM; j++) {
974  normal[j] += normals[i][j];
975  }
976  }
977  //Normal is now determined and we can free memory
978  normal_calculated = 1;
979  for (int i=0; i < Ntetrahedra; i++) {
980  delete[] normals[i];
981  }
982  delete[] a;
983  delete[] b;
984  delete[] c;
985  delete[] normals;
986  delete[] Vout;
987 }
988 
996 class Square
997 {
998  private:
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;
1003  double **points;
1004  double **cuts;
1005  double **out;
1006  double **points_temp;
1007  double *out_temp;
1008  int *const_i;
1009  double *const_value;
1010  int x1, x2;
1011  double *dx;
1012  int Ncuts;
1013  int Nlines;
1016  void ends_of_edge(double);
1017  void find_outside(double);
1018  public:
1019  Square();
1020  ~Square();
1021  void init(double**,int*,double*,double*);
1022  void construct_lines(double);
1023  int is_ambiguous();
1024  int get_Nlines();
1025  Line* get_lines();
1026 };
1027 
1035 {
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];
1041  lines = new Line[MAX_LINES];
1042  for (int i=0; i < SQUARE_DIM; i++) {
1043  points[i] = new double[SQUARE_DIM];
1044  }
1045  for (int i=0; i < MAX_POINTS; i++) {
1046  cuts[i] = new double[SQUARE_DIM];
1047  out[i] = new double[SQUARE_DIM];
1048  }
1049  int Npoints = 2;
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];
1054  }
1055  ambiguous = 0;
1056 }
1057 
1064 {
1065  for (int i=0; i < DIM-SQUARE_DIM; i++) {
1066  delete[] points[i];
1067  }
1068  delete[] points;
1069  for (int i=0; i < MAX_POINTS; i++) {
1070  delete[] cuts[i];
1071  delete[] out[i];
1072  }
1073  delete[] cuts;
1074  delete[] out;
1075  delete[] const_i;
1076  delete[] const_value;
1077  delete[] lines;
1078  for (int i=0; i < SQUARE_DIM; i++) {
1079  delete[] points_temp[i];
1080  }
1081  delete[] points_temp;
1082  delete[] out_temp;
1083 }
1084 
1097 void Square::init(double **sq, int *c_i, double *c_v, double *dex)
1098 {
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];
1102  }
1103  }
1104  for (int i=0; i < DIM-SQUARE_DIM; i++) {
1105  const_i[i] = c_i[i];
1106  const_value[i] = c_v[i];
1107  }
1108  dx = dex;
1109  x1 = -1;
1110  x2 = -1;
1111  for (int i=0; i < DIM; i++) {
1112  if ( i != const_i[0] && i != const_i[1] ) {
1113  if ( x1 < 0 ) {
1114  x1 = i;
1115  } else {
1116  x2 = i;
1117  }
1118  }
1119  }
1120  //We need to set these to zero so when this function is called
1121  //again we know that nothing has been done for this square
1122  Ncuts = 0;
1123  Nlines = 0;
1124  ambiguous = 0;
1125 }
1126 
1135 {
1136  //First we check the corner points and see if there are any lines
1137  int above=0;
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 )
1141  above++;
1142  }
1143  }
1144  //If all corners are above or below the surface value, no lines in this
1145  //square
1146  if ( above == 0 || above == 4 ) {
1147  Nlines = 0;
1148  return;
1149  }
1150  //First we find the cut points and the points which are always outside of the
1151  //surface. Also find_outside() arranges cuts so that first two cuts form a line
1152  //as defined in the algorithm (in case there are 4 cuts)
1153  ends_of_edge(E0);
1154  if ( Ncuts > 0 )
1155  find_outside(E0);
1156  //Then we go through the cut points and form the line elements
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];
1162  //When we have inserted both endpoints we insert outside point
1163  //and we are ready to create line element.
1164  if ( i%2 == 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);
1170  Nlines++;
1171  }
1172  }
1173 }
1174 
1192 void Square::ends_of_edge(double E0)
1193 {
1194  //Edge 1
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];
1197  cuts[Ncuts][1] = 0;
1198  Ncuts++;
1199  } else if ( points[0][0] == E0 && points[1][0] < E0 ) {
1200  cuts[Ncuts][0] = 1e-9*dx[x1];
1201  cuts[Ncuts][1] = 0;
1202  Ncuts++;
1203  } else if ( points[1][0] == E0 && points[0][0] < E0 ) {
1204  cuts[Ncuts][0] = (1.0-1e-9)*dx[x1];
1205  cuts[Ncuts][1] = 0;
1206  Ncuts++;
1207  }
1208  //Edge 2
1209  if ( (points[0][0]-E0)*(points[0][1]-E0) < 0 ) {
1210  cuts[Ncuts][0] = 0;
1211  cuts[Ncuts][1] = (points[0][0]-E0)/(points[0][0]-points[0][1])*dx[x2];
1212  Ncuts++;
1213  } else if ( points[0][0] == E0 && points[0][1] < E0 ) {
1214  cuts[Ncuts][0] = 0;
1215  cuts[Ncuts][1] = 1e-9*dx[x2];
1216  Ncuts++;
1217  } else if ( points[0][1] == E0 && points[0][0] < E0 ) {
1218  cuts[Ncuts][0] = 0;
1219  cuts[Ncuts][1] = (1.0-1e-9)*dx[x2];
1220  Ncuts++;
1221  }
1222  //Edge 3
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];
1226  Ncuts++;
1227  } else if ( points[1][0] == E0 && points[1][1] < E0 ) {
1228  cuts[Ncuts][0] = dx[x1];
1229  cuts[Ncuts][1] = 1e-9*dx[x2];
1230  Ncuts++;
1231  } else if ( points[1][1] == E0 && points[1][0] < E0 ) {
1232  cuts[Ncuts][0] = dx[x1];
1233  cuts[Ncuts][1] = (1.0-1e-9)*dx[x2];
1234  Ncuts++;
1235  }
1236  //Edge 4
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];
1240  Ncuts++;
1241  } else if ( points[0][1] == E0 && points[1][1] < E0 ) {
1242  cuts[Ncuts][0] = 1e-9*dx[x1];
1243  cuts[Ncuts][1] = dx[x2];
1244  Ncuts++;
1245  } else if ( points[1][1] == E0 && points[0][1] < E0 ) {
1246  cuts[Ncuts][0] = (1.0-1e-9)*dx[x1];
1247  cuts[Ncuts][1] = dx[x2];
1248  Ncuts++;
1249  }
1250  if ( Ncuts != 0 && Ncuts != 2 && Ncuts != 4 ) {
1251  cout << "Error in EndsOfEdge: Ncuts " << Ncuts << endl;
1252  exit(1);
1253  }
1254 }
1255 
1263 void Square::find_outside(double E0)
1264 {
1265  if ( Ncuts == 4 ) { //Here the surface is ambiguous
1266  ambiguous = 1;
1267  //Let's calculate the value in the middle of the square
1268  double Emid = 0;
1269  for (int i=0; i < 2; i++) {
1270  for (int j=0; j < 2; j++) {
1271  Emid += 0.25*points[i][j];
1272  }
1273  }
1274  //The default is that cuts are connected as \\ here
1275  //If both Emid and (0,0) are above or below the criterion
1276  //the cuts should be like // and we have to switch order in cuts
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];
1280  cuts[1][i] = cuts[2][i];
1281  cuts[2][i] = temp;
1282  }
1283  }
1284  //The center is below, so the middle point is always outside
1285  //the surface
1286  if ( (Emid-E0) < 0 ) {
1287  for (int i=0; i < 2; i++) {
1288  for (int j=0; j < 2; j++) {
1289  if ( j == 0 )
1290  out[i][j] = 0.5*dx[x1];
1291  else
1292  out[i][j] = 0.5*dx[x2];
1293  }
1294  }
1295  } else { // The center is above
1296  // Cuts are \\ here so bl and tr corners are outside
1297  if ( (points[0][0]-E0) < 0 ) {
1298  for (int i=0; i < 2; i++) {
1299  out[0][i] = 0;
1300  if ( i == 0 )
1301  out[1][i] = dx[x1];
1302  else
1303  out[1][i] = dx[x2];
1304  }
1305  // Cuts are // here so br and tl corners are outside
1306  } else {
1307  out[0][0] = dx[x1];
1308  out[0][1] = 0;
1309  out[1][0] = 0;
1310  out[1][1] = dx[x2];
1311  }
1312  }
1313  } else { //This is the normal case (not ambiguous)
1314  for (int i=0; i < 2; i++) {
1315  for (int j=0; j < 2; j++) {
1316  out[i][j] = 0;
1317  }
1318  }
1319  int Nout = 0;
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];
1325  Nout++;
1326  }
1327  }
1328  }
1329  if ( Nout > 0 ) {
1330  for (int i=0; i < 2; i++) {
1331  for (int j=0; j < 2; j++) {
1332  out[i][j] = out[i][j]/double(Nout);
1333  }
1334  }
1335  }
1336  }
1337 }
1338 
1347 {
1348  return ambiguous;
1349 }
1350 
1359 {
1360  return Nlines;
1361 }
1362 
1371 {
1372  return lines;
1373 }
1374 
1375 
1382 class Cube
1383 {
1384  private:
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;
1390  double ***cube;
1394  int Nlines;
1397  int const_i;
1398  double const_value;
1399  int x1,x2,x3;
1400  double *dx;
1401  void split_to_squares();
1402  void check_ambiguous(int);
1403  public:
1404  Cube();
1405  ~Cube();
1406  void init(double***&,int,double,double*&);
1407  void construct_polygons(double);
1408  int get_Nlines();
1409  int get_Npolygons();
1410  int is_ambiguous();
1411  Polygon* get_polygons();
1412 };
1413 
1420 {
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];
1426  }
1427  }
1428  lines = new Line*[NSQUARES*2]; //Each square may have max. 2 lines
1429  polygons = new Polygon[MAX_POLY];
1430  squares = new Square[NSQUARES];
1431 }
1432 
1439 {
1440  for (int i=0; i < STEPS; i++) {
1441  for (int j=0; j < STEPS; j++) {
1442  delete[] cube[i][j];
1443  }
1444  delete[] cube[i];
1445  }
1446  delete[] cube;
1447  delete[] lines;
1448  delete[] polygons;
1449  delete[] squares;
1450 }
1451 
1463 void Cube::init(double ***&c, int c_i, double c_v, double *&dex)
1464 {
1465  const_i = c_i;
1466  const_value = c_v;
1467  dx = dex;
1468  //Here we fix the non-zero indices
1469  x1 = -1;
1470  x2 = -1;
1471  x3 = -1;
1472  for (int i=0; i < DIM; i++) {
1473  if ( i != c_i ) {
1474  if ( x1 < 0 ) {
1475  x1 = i;
1476  } else if ( x2 < 0 ) {
1477  x2 = i;
1478  } else {
1479  x3 = i;
1480  }
1481  }
1482  }
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];
1487  }
1488  }
1489  }
1490  //We need to set these to zero so when this function is called
1491  //again we know that nothing has been done for this cube
1492  Nlines = 0;
1493  Npolygons = 0;
1494  ambiguous = 0;
1495 }
1496 
1503 {
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];
1509  }
1510  int Nsquares = 0;
1511  for (int i=0; i < DIM; i++) {
1512  //i is the index which is kept constant, thus we ignore the index which
1513  //is constant in this cube
1514  if ( i != const_i ) {
1515  c_i[0] = const_i;
1516  c_i[1] = i;
1517  for (int j=0; j < STEPS; j++) {
1518  c_v[0] = const_value;
1519  c_v[1] = j*dx[i];
1520  for (int ci1=0; ci1 < STEPS; ci1++) {
1521  for (int ci2=0; ci2 < STEPS; ci2++) {
1522  if ( i == x1 ) {
1523  sq[ci1][ci2] = cube[j][ci1][ci2];
1524  } else if ( i == x2 ) {
1525  sq[ci1][ci2] = cube[ci1][j][ci2];
1526  } else {
1527  sq[ci1][ci2] = cube[ci1][ci2][j];
1528  }
1529  }
1530  }
1531  squares[Nsquares].init(sq,c_i,c_v,dx);
1532  Nsquares++;
1533  }
1534  }
1535  }
1536  for (int i=0; i < STEPS; i++) {
1537  delete[] sq[i];
1538  }
1539  delete[] sq;
1540  delete[] c_i;
1541  delete[] c_v;
1542 }
1543 
1552 void Cube::construct_polygons(double value0)
1553 {
1554  //Let's start by splitting the cube to squares and finding
1555  //the lines from them
1556  split_to_squares();
1557  for (int i=0; i < NSQUARES; i++) {
1558  squares[i].construct_lines(value0);
1559  }
1560  //Then we make a table which contains pointers to the lines
1561  Nlines = 0;
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++) {
1566  lines[Nlines] = &l[j];
1567  double *p = lines[Nlines]->get_start();
1568  p = lines[Nlines]->get_end();
1569  Nlines++;
1570  }
1571  l = NULL;
1572  }
1573  //If no lines were found we may exit. This can happen only in 4D case
1574  if ( Nlines == 0 ) {
1575  return;
1576  }
1577  //Then we check if the surface is ambiguous and continue accordingly
1578  check_ambiguous(Nlines);
1579  if ( ambiguous > 0 ) {
1580  //Surface is ambiguous, so let's connect the lines to polygons and see how
1581  //many polygons we have
1582  int *not_used = new int[Nlines];
1583  for (int i=0; i < Nlines; i++) {
1584  not_used[i] = 1;
1585  }
1586  int used = 0; //this keeps track how many lines we have used
1587  //We may found several polygons with this
1588  do {
1589  if ( Nlines - used < 3 ) {
1590  cout << "Error: cannot construct polygon from " << Nlines -used << " lines" << endl;
1591  exit(1);
1592  }
1593  //First we initialize polygon
1594  polygons[Npolygons].init(const_i);
1595  //Then we go through all lines and try to add them to polygon
1596  for (int i=0; i < Nlines; i++) {
1597  if ( not_used[i] ) {
1598  //add_line returns true if line is succesfully added
1599  if ( polygons[Npolygons].add_line(lines[i],0) ) {
1600  not_used[i]=0;
1601  used++;
1602  //if line is succesfully added we start the loop from the
1603  //beginning
1604  i=0;
1605  }
1606  }
1607  }
1608  //When we have reached this point one complete polygon is formed
1609  Npolygons++;
1610  } while ( used < Nlines );
1611  delete[] not_used;
1612  } else {
1613  //Surface is not ambiguous, so we have only one polygons and all lines
1614  //can be added to it without ordering them
1615  polygons[Npolygons].init(const_i);
1616  for (int i=0; i < Nlines; i++) {
1617  polygons[Npolygons].add_line(lines[i],1);
1618  }
1619  Npolygons++;
1620  }
1621 }
1622 
1628 void Cube::check_ambiguous(int Nlines)
1629 {
1630  //First we check if any squares may have ambiguous elements
1631  for (int i=0; i < NSQUARES; i++) {
1632  if ( squares[i].is_ambiguous() ) {
1633  ambiguous++;
1634  }
1635  }
1636  //If the surface is not ambigous already, it is still possible to
1637  //have a ambiguous case if we have exatcly 6 lines, i.e. the surface
1638  //elements are at the opposite corners
1639  if ( ambiguous == 0 && Nlines == 6 ) {
1640  ambiguous++;
1641  }
1642 }
1643 
1653 {
1654  return ambiguous;
1655 }
1656 
1665 {
1666  return Nlines;
1667 }
1668 
1677 {
1678  return Npolygons;
1679 }
1680 
1689 {
1690  return polygons;
1691 }
1692 
1700 {
1701  private:
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;
1706  double ****hcube;
1712  int x1,x2,x3,x4;
1713  double *dx;
1714  void split_to_cubes();
1715  void check_ambiguous(double);
1716  public:
1717  Hypercube();
1718  ~Hypercube();
1719  void init(double****,double*);
1720  void construct_polyhedrons(double);
1721  int get_Npolyhedrons();
1722  Polyhedron* get_polyhedrons();
1723 };
1724 
1731 {
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];
1739  }
1740  }
1741  }
1742  polyhedrons = new Polyhedron[MAX_POLY];
1743  polygons = new Polygon*[NCUBES*10];
1744  cubes = new Cube[NCUBES];
1745 }
1746 
1753 {
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];
1758  }
1759  delete[] hcube[i][j];
1760  }
1761  delete[] hcube[i];
1762  }
1763  delete[] hcube;
1764  delete[] polyhedrons;
1765  delete[] polygons;
1766  delete[] cubes;
1767 }
1768 
1778 void Hypercube::init(double ****c, double *dex)
1779 {
1780  dx = dex;
1781  //Here we fix the non-zero indices
1782  x1 = 0;
1783  x2 = 1;
1784  x3 = 2;
1785  x4 = 3;
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];
1791  }
1792  }
1793  }
1794  }
1795  //We need to set these to zero so we can use this function to initialize
1796  //many hypercubes
1797  Npolyhedrons = 0;
1798  ambiguous = 0;
1799 }
1800 
1807 {
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];
1813  }
1814  }
1815  int Ncubes = 0;
1816  for (int i=0; i < DIM; i++) {
1817  //i is the index which is kept constant, thus we ignore the index which
1818  //is constant in this cube
1819  int c_i = 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++) {
1825  if ( i == x1 ) {
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];
1831  } else {
1832  cu[ci1][ci2][ci3] = hcube[ci1][ci2][ci3][j];
1833  }
1834  }
1835  }
1836  }
1837  cubes[Ncubes].init(cu,c_i,c_v,dx);
1838  Ncubes++;
1839  }
1840  }
1841  for (int i=0; i < STEPS; i++) {
1842  for (int j=0; j < STEPS; j++) {
1843  delete[] cu[i][j];
1844  }
1845  delete[] cu[i];
1846  }
1847  delete[] cu;
1848 }
1849 
1859 {
1860  split_to_cubes();
1861  for (int i=0; i < NCUBES; i++) {
1862  cubes[i].construct_polygons(value0);
1863  }
1864  int Npolygons = 0;
1865  //then polygons are loaded to table
1866  for (int i=0; i < NCUBES; i++) {
1867  int Npoly = cubes[i].get_Npolygons();
1868  Polygon *p = cubes[i].get_polygons();
1869  for (int j=0; j < Npoly; j++) {
1870  polygons[Npolygons] = &p[j];
1871  Npolygons++;
1872  }
1873  }
1874  check_ambiguous(value0);
1875  if ( ambiguous > 0 ) {
1876  //Here surface might be ambiguous and we need to connect the polygons and
1877  //see how many polyhedrons we have
1878  int *not_used = new int[Npolygons];
1879  for (int i=0; i < Npolygons; i++) {
1880  not_used[i] = 1;
1881  }
1882  int used = 0; //this keeps track how many lines we have used
1883  //We may found several polyhedrons with this
1884  do {
1885  //First we initialize polyhedron
1886  polyhedrons[Npolyhedrons].init();
1887  //Then we go through all polygons and try to add them to polyhedron
1888  for (int i=0; i < Npolygons; i++) {
1889  if ( not_used[i] ) {
1890  //add_polygon returns true if the polygon is succesfully added
1891  if ( polyhedrons[Npolyhedrons].add_polygon(polygons[i],0) ) {
1892  not_used[i]=0;
1893  used++;
1894  //if polygon is succesfully added we start the loop from the
1895  //beginning
1896  i=0;
1897  }
1898  }
1899  }
1900  //When we have reached this point one complete polyhedron is formed
1901  Npolyhedrons++;
1902  } while ( used < Npolygons );
1903  delete[] not_used;
1904  /*if ( ambiguous == 0 && Npolyhedrons != 1 ) {
1905  cout << "error" << endl;
1906  }*/
1907  } else {
1908  //Here surface cannot be ambiguous and all polygons can be added to
1909  //the polygehdron without ordering them
1910  polyhedrons[Npolyhedrons].init();
1911  for (int i=0; i < Npolygons; i++) {
1912  polyhedrons[Npolyhedrons].add_polygon(polygons[i],1);
1913  }
1914  //When we have reached this point one complete polyhedron is formed
1915  Npolyhedrons++;
1916  }
1917 }
1918 
1924 void Hypercube::check_ambiguous(double value0)
1925 {
1926  for (int i=0; i < NCUBES; i++) {
1927  if ( cubes[i].is_ambiguous() ) {
1928  ambiguous++;
1929  }
1930  }
1931  if ( ambiguous == 0 ) {
1932  int Nlines = 0;
1933  for (int i=0; i < NCUBES; i++) {
1934  Nlines += cubes[i].get_Nlines();
1935  }
1936  int points = 0;
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 ) {
1942  points++;
1943  }
1944  }
1945  }
1946  }
1947  }
1948  if ( points > 8 ) {
1949  points = 16-points;
1950  }
1951  /* these are not needed
1952  if ( Nlines == 42 ) {
1953  //ambiguous++;
1954  } else if ( Nlines == 36 && ( points == 5 || points == 4 ) ) {
1955  //ambiguous++;
1956  } else if ( Nlines == 30 && points == 3 ) {
1957  //ambiguous++;
1958  }*/
1959  if ( Nlines == 24 && points == 2 ) {
1960  ambiguous++;
1961  }
1962  }
1963 }
1964 
1973 {
1974  return Npolyhedrons;
1975 }
1976 
1985 {
1986  return polyhedrons;
1987 }
1988 
2000 {
2001  private:
2002  static const int STEPS = 2;
2003  static const int DIM = 4;
2004  static const int MAX_ELEMENTS = 10;
2006  double **normals;
2007  double **centroids;
2011  double value0;
2012  double *dx;
2013  ofstream output_print;
2014  void surface_3d(double***,double*,int);
2018  public:
2019  Cornelius();
2020  ~Cornelius();
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);
2034 };
2035 
2042 {
2043  Nelements = 0;
2044  initialized = 0;
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];
2051  }
2052 }
2053 
2060 {
2061  for (int i=0; i < MAX_ELEMENTS; i++) {
2062  delete[] normals[i];
2063  delete[] centroids[i];
2064  }
2065  delete[] normals;
2066  delete[] centroids;
2067  if ( initialized ) {
2068  delete[] dx;
2069  }
2070  //If file for printing was opened we close it here.
2071  if ( print_initialized ) {
2072  output_print.close();
2073  }
2074 }
2075 
2087 void Cornelius::init(int d, double v0, double *dex)
2088 {
2089  cube_dim = d;
2090  value0 = v0;
2091  if ( initialized != 1 ) {
2092  dx = new double[DIM];
2093  }
2094  for (int i=0; i < DIM; i++) {
2095  if ( i < DIM-cube_dim ) {
2096  dx[i] = 1;
2097  } else {
2098  dx[i] = dex[i-(DIM-cube_dim)];
2099  }
2100  }
2101  initialized = 1;
2102 }
2103 
2114 {
2115  output_print.open(filename.c_str());
2116  print_initialized = 1;
2117 }
2118 
2127 void Cornelius::find_surface_2d(double **cube)
2128 {
2129  if ( !initialized || cube_dim != 2 ) {
2130  cout << "Cornelius not initialized for 2D case" << endl;
2131  exit(1);
2132  }
2133  int *c_i = new int[cube_dim];
2134  double *c_v = new double[cube_dim];
2135  c_i[0] = 0;
2136  c_i[1] = 1;
2137  c_v[0] = 0;
2138  c_v[1] = 0;
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++) {
2145  normals[i][j] = l[i].get_normal()[j];
2146  centroids[i][j] = l[i].get_centroid()[j];
2147  }
2148  }
2149  delete[] c_i;
2150  delete[] c_v;
2151 }
2152 
2161 void Cornelius::find_surface_3d(double ***cube)
2162 {
2163  double *pos = NULL;
2164  surface_3d(cube,pos,0);
2165 }
2166 
2177 void Cornelius::find_surface_3d_print(double ***cube, double *pos)
2178 {
2179  surface_3d(cube,pos,1);
2180 }
2181 
2192 void Cornelius::surface_3d(double ***cube, double *pos, int do_print)
2193 {
2194  if ( !initialized || cube_dim != 3 ) {
2195  cout << "Cornelius not initialized for 3D case" << endl;
2196  exit(1);
2197  }
2198  //First we check if the cube actually contains surface elements.
2199  //If all or none of the elements are below the criterion, no surface
2200  //elements exist.
2201  int above = 0;
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 )
2206  above++;
2207  }
2208  }
2209  }
2210  if ( above == 0 || above == 8 ) {
2211  //No elements in this cube
2212  Nelements = 0;
2213  return;
2214  }
2215  //This cube has surface elements, so let's start by constructing
2216  //the cube.
2217  int c_i = 0;
2218  double c_v = 0;
2219  cu3d.init(cube,c_i,c_v,dx);
2220  //Then we find the elements
2221  cu3d.construct_polygons(value0);
2222  //Now let's get the information about the elements
2223  Nelements = cu3d.get_Npolygons();
2224  Polygon *p = cu3d.get_polygons();
2225  for (int i=0; i < Nelements; i++) {
2226  //Here we always work with 4-dimensions
2227  for (int j=0; j < DIM; j++) {
2228  normals[i][j] = p[i].get_normal()[j];
2229  centroids[i][j] = p[i].get_centroid()[j];
2230  }
2231  //If we want to print the actual triangles, which are found, we must do it here
2232  //before polygons are removed from the memory.
2233  if ( print_initialized && do_print ) {
2234  p[i].print(output_print,pos);
2235  }
2236  }
2237 }
2238 
2248 void Cornelius::find_surface_4d(double ****cube)
2249 {
2250  if ( !initialized || cube_dim != 4 ) {
2251  cout << "Cornelius not initialized for 4D case" << endl;
2252  exit(1);
2253  }
2254  //First we check if the cube actually contains surface elements.
2255  //If all or none of the elements are below the criterion, no surface
2256  //elements exist.
2257  int above = 0;
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 )
2263  above++;
2264  }
2265  }
2266  }
2267  }
2268  if ( above == 0 || above == 16 ) {
2269  //No elements in this cube
2270  Nelements = 0;
2271  return;
2272  }
2273  //This cube has surface elements, so let's start by constructing
2274  //the hypercube.
2275  cu4d.init(cube,dx);
2276  //Then we find the elements
2277  cu4d.construct_polyhedrons(value0);
2278  //Now let's get the information about the elements
2279  Nelements = cu4d.get_Npolyhedrons();
2280  Polyhedron *p = cu4d.get_polyhedrons();
2281  for (int i=0; i < Nelements; i++) {
2282  for (int j=0; j < DIM; j++) {
2283  centroids[i][j] = p[i].get_centroid()[j];
2284  normals[i][j] = p[i].get_normal()[j];
2285  }
2286  }
2287 }
2288 
2297 {
2298  return Nelements;
2299 }
2300 
2313 {
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];
2319  }
2320  }
2321  return vect;
2322 }
2323 
2337 {
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];
2343  }
2344  }
2345  return vect;
2346 }
2347 
2359 {
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)];
2365  }
2366  }
2367  return vect;
2368 }
2369 
2382 {
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)];
2388  }
2389  }
2390  return vect;
2391 }
2392 
2405 {
2406  if ( i >= Nelements || j >= cube_dim ) {
2407  cout << "Cornelius error: asking for an element which does not exist." << endl;
2408  exit(1);
2409  }
2410  return centroids[i][j+(DIM-cube_dim)];
2411 }
2412 
2426 {
2427  if ( i >= Nelements || j >= cube_dim ) {
2428  cout << "Cornelius error: asking for an element which does not exist." << endl;
2429  exit(1);
2430  }
2431  return normals[i][j+(DIM-cube_dim)];
2432 }