initial optimization of revised intersect_line_segment_with_prism (#21)

* update point_in_prism to avoid double-counting of polygon edge intersections when plumb lines pass through vertices

* revised normal_to_object implementation for prisms to account for in-plane distance from point to prism face when determining closest face

* re-updated normal_to_prism to behave correctly when the projection of the point into the plane does not lie within the prism face in question

* updates

* revise algorithm for intersecting line segment with prism

* overhauled intersect_line_segment_with_prism to do the full exact calculation with no approximations or assumptions about proximity

* update acquisition of array slices for DFT fields to behave properly in the presence of symmetries (2)

* update to intersect_line_segment_with_prism to fix failing unit test

* remove memory allocation  from intersect_line_segment_with_prism

Author: Homer Reid

utils/geom.c

@@ -2281,7 +2281,7 @@ int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double *slist)
 /***************************************************************/
 double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b)
 {
-  double *slist = (double *)malloc( (2+prsm->vertices.num_items)*sizeof(double) );
+  double *slist=prsm->workspace.items;
   int num_intersections=intersect_line_with_prism(prsm, p, d, slist);
 
   // na=smallest index such that slist[na] > a
@@ -2292,9 +2292,7 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, doub
     na=ns;
 
   if (na==-1)
-   { free(slist);
-     return 0.0;
-   }
+   return 0.0;
 
   int inside = ( (na%2)==0 ? 0 : 1);
   double last_s=a;
@@ -2306,7 +2304,6 @@ double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, doub
      inside = (1-inside);
      last_s = this_s;
    }
-  free(slist);
   return ds > 0.0 ? ds : 0.0;
 }
 
@@ -2577,6 +2574,11 @@ geometric_object make_prism(material_type material,
   for(nv=0; nv<num_vertices; nv++)
    prsm->vertices.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]);
 
+  // workspace is an internally-stored double-valued array of length num_vertices+2
+  // that is used by some geometry routines
+  prsm->workspace.num_items = num_vertices+2;
+  prsm->workspace.items     = (double *)malloc( (num_vertices+2)*sizeof(double) );
+
   // note the center and centroid are different!
   vector3 center = vector3_plus(centroid, vector3_scale(0.5*height,zhat) );
   geometric_object o=make_geometric_object(material, center);

utils/geom.scm

@@ -93,6 +93,7 @@
   (define-property vertices '() (make-list-type 'vector3))
   (define-property centroid (vector3 0 0 0) 'vector3)
   (define-property height 0 'number)
+  (define-property workspace '() (make-list-type 'number))
   (define-property m_c2p identity_matrix 'matrix3x3)
   (define-property m_p2c identity_matrix 'matrix3x3))