Union: polyhedron display (#2898) combined with Union_mesh BVH display (#2893)

mcstas-comps/share/union-lib.c

@@ -1390,38 +1390,36 @@ void print_intersection_table(struct intersection_time_table_struct *intersectio
     }
 
 void merge_lines_to_draw(struct lines_to_draw *lines_master,struct lines_to_draw *lines_new) {
-    if (lines_master->number_of_lines == 0) {
-    lines_master->number_of_lines = lines_new->number_of_lines;
-    if (!lines_master->number_of_lines) {
-      return;
-    }
-    lines_master->lines = malloc(lines_master->number_of_lines*sizeof(struct line_segment));
-    if (!lines_master->lines) {
-      fprintf(stderr,"Failure allocating list in Union function merge_lines_to_draw 2 - Exit!\n");
-      exit(EXIT_FAILURE);
-    }
-    lines_master->lines = lines_new->lines;
-    // One could free lines_new->lines;
-    } else {
-    int iterate;
-    struct line_segment *temp_lines;
-    temp_lines = malloc(lines_master->number_of_lines*sizeof(struct line_segment));
-    if (!temp_lines) {
-      fprintf(stderr,"Failure allocating list in Union function merge_lines_to_draw 3 - Exit!\n");
-      exit(EXIT_FAILURE);
-    }
-    for (iterate = 0;iterate < lines_master->number_of_lines;iterate++) temp_lines[iterate] = lines_master->lines[iterate];
-    free(lines_master->lines);
-    lines_master->lines = malloc((lines_master->number_of_lines+lines_new->number_of_lines)*sizeof(struct line_segment));
-    if (!lines_master->lines) {
-      fprintf(stderr,"Failure allocating list in Union function merge_lines_to_draw 4 - Exit!\n");
-      exit(EXIT_FAILURE);
-    }
-    for (iterate = 0;iterate < lines_master->number_of_lines;iterate++) lines_master->lines[iterate] = temp_lines[iterate];
-    for (iterate = 0;iterate < lines_new->number_of_lines;iterate++) lines_master->lines[iterate+lines_master->number_of_lines] = lines_new->lines[iterate];
-    lines_master->number_of_lines = lines_master->number_of_lines + lines_new->number_of_lines;
-    free(temp_lines);
+    // Appends the lines of lines_new to lines_master, then frees lines_new->lines
+    // and leaves lines_new empty.
+    // The allocation behind lines_master->lines is only ever made here, and its
+    // capacity is always the smallest power of two holding number_of_lines, so it
+    // need not be stored: appending one line at a time costs amortized O(1).
+    // A master is started by setting number_of_lines = 0.
+    int capacity, needed, iterate;
+    if (lines_new->number_of_lines <= 0) return;
+    needed = lines_master->number_of_lines + lines_new->number_of_lines;
+    capacity = 0; // an empty master owns no allocation (lines may be unset)
+    if (lines_master->number_of_lines > 0) {
+      capacity = 1;
+      while (capacity < lines_master->number_of_lines) capacity *= 2;
+    }
+    if (needed > capacity) {
+      int new_capacity = capacity > 0 ? capacity : 1;
+      struct line_segment *grown;
+      while (new_capacity < needed) new_capacity *= 2;
+      grown = realloc(capacity > 0 ? lines_master->lines : NULL, new_capacity*sizeof(struct line_segment));
+      if (!grown) {
+        fprintf(stderr,"Failure allocating list in Union function merge_lines_to_draw - Exit!\n");
+        exit(EXIT_FAILURE);
+      }
+      lines_master->lines = grown;
     }
+    for (iterate = 0;iterate < lines_new->number_of_lines;iterate++) lines_master->lines[lines_master->number_of_lines+iterate] = lines_new->lines[iterate];
+    lines_master->number_of_lines = needed;
+    free(lines_new->lines);
+    lines_new->lines = NULL;
+    lines_new->number_of_lines = 0;
     }
 
 int r_has_highest_priority(Coords point,int N,struct geometry_struct **Geometries,int number_of_volumes) {
@@ -1589,6 +1587,7 @@ struct lines_to_draw draw_line_with_highest_priority(Coords position1,Coords pos
         }
     }
     free(double_dummy);
+    free(int_dummy);
     free(temp_intersection);
     // Now we have a list of intersection distances between r1 and r2 and all volumes.
     // This list needs to be sorted before we continue!
@@ -1652,9 +1651,9 @@ struct lines_to_draw draw_line_with_highest_priority(Coords position1,Coords pos
               draw_order.lines[draw_order.number_of_lines++] = lines[iterate];
             }
         }
-        if (intersection_list.num_elements > 0) free(intersection_list.elements);
     }
-
+    if (intersection_list.num_elements > 0) free(intersection_list.elements);
+
     free(points);
     free(lines);
     free(draw_logic);
@@ -1805,14 +1804,14 @@ union geometry_parameter_union allocate_cone_storage_copy(union geometry_paramet
 
 union geometry_parameter_union allocate_mesh_copy(union geometry_parameter_union *union_input) {
   union geometry_parameter_union union_output;
-  // Allocate the space for a gmesh_struct in the new union_output (union as the c structre)
-  union_output.p_mesh = malloc(sizeof(gmesh_struct));
-  if (!union_output.p_mesh) {
+  // Allocate the space for a mesh_storage in the new union_output (union as the c structre)
+  union_output.p_mesh_storage = malloc(sizeof(struct mesh_storage));
+  if (!union_output.p_mesh_storage) {
     fprintf(stderr,"Failure allocating list in Union function allocate_mesh_copy - Exit!\n");
     exit(EXIT_FAILURE);
   }
   // Copy the input storage to the output
-  *union_output.p_mesh = *union_input->p_mesh;
+  *union_output.p_mesh_storage = *union_input->p_mesh_storage;
 
   return union_output;
 }
@@ -1848,6 +1847,136 @@ Coords point_on_circle(Coords center, Coords direction, double radius, int point
     return output;
 }
 
+// Polyhedron display: triangulate the volume surface and keep the triangles where it has highest priority
+struct union_tri_list {
+  Coords *v; // 3 vertices per triangle
+  int n, size;
+  int N, number_of_volumes;
+  struct geometry_struct **Geometries;
+};
+
+void union_tri_emit(struct union_tri_list *t, Coords a, Coords b, Coords c) {
+  if (t->n == t->size) {
+    t->size = t->size ? 2*t->size : 256;
+    t->v = realloc(t->v, 3*t->size*sizeof(Coords));
+    if (!t->v) {
+      fprintf(stderr,"Failure allocating list in Union function union_tri_emit - Exit!\n");
+      exit(EXIT_FAILURE);
+    }
+  }
+  t->v[3*t->n] = a; t->v[3*t->n+1] = b; t->v[3*t->n+2] = c;
+  t->n++;
+}
+
+#define UNION_TRI_DEPTH 6
+// Splits while the priority test disagrees between the 4 sub-triangles, down to depth 0. Depth above UNION_TRI_DEPTH forces a split.
+// ponytail: priority sampled only at sub-triangle centroids, an intersecting feature that hits none of them is missed
+void union_tri_add(struct union_tri_list *t, Coords a, Coords b, Coords c, int depth) {
+  Coords ab = coords_scalar_mult(coords_add(a,b),0.5);
+  Coords bc = coords_scalar_mult(coords_add(b,c),0.5);
+  Coords ca = coords_scalar_mult(coords_add(c,a),0.5);
+  Coords sub[4][3] = {{a,ab,ca},{ab,b,bc},{ca,bc,c},{ab,bc,ca}};
+  int keep[4], i;
+  for (i=0;i<4;i++)
+    keep[i] = r_has_highest_priority(coords_scalar_mult(coords_add(coords_add(sub[i][0],sub[i][1]),sub[i][2]),1.0/3.0),
+                                     t->N,t->Geometries,t->number_of_volumes) != 0;
+  int mixed = keep[0]!=keep[1] || keep[0]!=keep[2] || keep[0]!=keep[3];
+  if (!mixed && depth <= UNION_TRI_DEPTH) {
+    if (keep[0]) union_tri_emit(t,a,b,c);
+    return;
+  }
+  for (i=0;i<4;i++) {
+    if (depth > 0) union_tri_add(t,sub[i][0],sub[i][1],sub[i][2],depth-1);
+    else if (keep[i]) union_tri_emit(t,sub[i][0],sub[i][1],sub[i][2]);
+  }
+}
+
+// Surface of revolution through profile points (z[k], r[k]) along axis, with caps where r != 0 at the ends
+void union_tri_lathe(struct union_tri_list *t, Coords center, Coords axis, double *z, double *r, int np) {
+  int nseg = 32, i, k;
+  Coords c0 = coords_add(center,coords_scalar_mult(axis,z[0]));
+  Coords c1 = coords_add(center,coords_scalar_mult(axis,z[np-1]));
+  for (i=0;i<nseg;i++) {
+    for (k=0;k<np-1;k++) {
+      Coords zk = coords_add(center,coords_scalar_mult(axis,z[k]));
+      Coords zk1 = coords_add(center,coords_scalar_mult(axis,z[k+1]));
+      Coords p00 = point_on_circle(zk,axis,r[k],i,nseg),   p01 = point_on_circle(zk,axis,r[k],i+1,nseg);
+      Coords p10 = point_on_circle(zk1,axis,r[k+1],i,nseg), p11 = point_on_circle(zk1,axis,r[k+1],i+1,nseg);
+      if (r[k] > 0)   union_tri_add(t,p00,p01,p11,UNION_TRI_DEPTH);
+      if (r[k+1] > 0) union_tri_add(t,p00,p11,p10,UNION_TRI_DEPTH);
+    }
+    if (r[0] > 0) union_tri_add(t,c0,point_on_circle(c0,axis,r[0],i+1,nseg),point_on_circle(c0,axis,r[0],i,nseg),UNION_TRI_DEPTH);
+    if (r[np-1] > 0) union_tri_add(t,c1,point_on_circle(c1,axis,r[np-1],i,nseg),point_on_circle(c1,axis,r[np-1],i+1,nseg),UNION_TRI_DEPTH);
+  }
+}
+
+// Draws volume N as one polyhedron in the frame of the master component, returns 0 if the volume should be drawn as lines instead
+int union_mcdisplay_polyhedron(int N, struct geometry_struct **Geometries, int number_of_volumes) {
+  struct geometry_struct *g = Geometries[N];
+  struct union_tri_list t = {NULL,0,0,N,number_of_volumes,Geometries};
+  int i, j;
+
+  if (g->eShape == box) {
+    struct box_storage *b = g->geometry_parameters.p_box_storage;
+    Coords s[2][4];
+    double w[2] = {b->x_width1,b->x_width2}, h[2] = {b->y_height1,b->y_height2}, zz[2] = {-0.5,0.5};
+    for (j=0;j<2;j++) {
+      s[j][0] = coords_add(coords_add(coords_add(g->center,coords_scalar_mult(b->z_vector,zz[j]*b->z_depth)),
+                                      coords_scalar_mult(b->x_vector,-0.5*w[j])),coords_scalar_mult(b->y_vector,-0.5*h[j]));
+      s[j][1] = coords_add(s[j][0],coords_scalar_mult(b->x_vector,w[j]));
+      s[j][2] = coords_add(s[j][1],coords_scalar_mult(b->y_vector,h[j]));
+      s[j][3] = coords_add(s[j][0],coords_scalar_mult(b->y_vector,h[j]));
+    }
+    for (j=0;j<2;j++) {
+      union_tri_add(&t,s[j][0],s[j][1],s[j][2],UNION_TRI_DEPTH+2);
+      union_tri_add(&t,s[j][0],s[j][2],s[j][3],UNION_TRI_DEPTH+2);
+    }
+    for (i=0;i<4;i++) {
+      union_tri_add(&t,s[0][i],s[0][(i+1)%4],s[1][(i+1)%4],UNION_TRI_DEPTH+2);
+      union_tri_add(&t,s[0][i],s[1][(i+1)%4],s[1][i],UNION_TRI_DEPTH+2);
+    }
+  } else if (g->eShape == cylinder) {
+    struct cylinder_storage *c = g->geometry_parameters.p_cylinder_storage;
+    double z[2] = {-0.5*c->height,0.5*c->height}, r[2] = {c->cyl_radius,c->cyl_radius};
+    union_tri_lathe(&t,g->center,c->direction_vector,z,r,2);
+  } else if (g->eShape == cone) {
+    struct cone_storage *c = g->geometry_parameters.p_cone_storage;
+    double z[2] = {-0.5*c->height,0.5*c->height}, r[2] = {c->cone_radius_bottom,c->cone_radius_top};
+    union_tri_lathe(&t,g->center,c->direction_vector,z,r,2);
+  } else if (g->eShape == sphere) {
+    double R = g->geometry_parameters.p_sphere_storage->sph_radius, z[17], r[17];
+    for (i=0;i<17;i++) { z[i] = -R*cos(PI*i/16); r[i] = (i==0 || i==16) ? 0 : R*sin(PI*i/16); }
+    union_tri_lathe(&t,g->center,coords_set(0,1,0),z,r,17);
+  } else if (g->eShape == mesh) {
+    if (g->geometry_parameters.p_mesh_storage->display_layer > 0) return 0; // coarse BVH boxes are drawn as lines
+    gmesh_struct *m = &g->geometry_parameters.p_mesh_storage->mesh;
+    struct pointer_to_1d_coords_list pts = g->shell_points(g,1); // all vertices, master frame
+    for (i=0;i<m->faceSize;i+=m->faceArray[i]+1)
+      for (j=1;j<m->faceArray[i]-1;j++) // fan over each polygon
+        union_tri_add(&t,pts.elements[m->faceArray[i+1]],pts.elements[m->faceArray[i+1+j]],pts.elements[m->faceArray[i+2+j]],UNION_TRI_DEPTH);
+    free(pts.elements);
+  }
+  if (t.n == 0) return 1;
+
+  // Same json layout as mcdis_polygon, bounds per vertex / face as there
+  char *json = malloc(64 + (size_t)t.n*(3*48 + 64));
+  if (!json) {
+    fprintf(stderr,"Failure allocating list in Union function union_mcdisplay_polyhedron - Exit!\n");
+    exit(EXIT_FAILURE);
+  }
+  char *p = json + sprintf(json,"{ \"vertices\": [");
+  for (i=0;i<3*t.n;i++)
+    p += sprintf(p,"%s[%g, %g, %g]",i ? ", " : "",t.v[i].x,t.v[i].y,t.v[i].z);
+  p += sprintf(p,"], \"faces\": [");
+  for (i=0;i<t.n;i++)
+    p += sprintf(p,"%s{ \"face\": [%d, %d, %d]}",i ? ", " : "",3*i,3*i+1,3*i+2);
+  sprintf(p,"]}");
+  mcdis_polyhedron(json);
+  free(json);
+  free(t.v);
+  return 1;
+}
+
 void points_on_circle(Coords *output, Coords center, Coords direction, double radius, int number_of_points) {
 
     Coords cross_input = coords_set(0,1,0);
@@ -2760,7 +2889,7 @@ This function was created by Martin Olsen at NBI on september 20, 2018.
 int r_within_mesh(Coords pos,struct geometry_struct *geometry) {
   // Point in mesh test (parity of ray crossings) in the frame of the mesh component
   Coords local_position = rot_apply(geometry->transpose_rotation_matrix, coords_sub(pos, geometry->center));
-  return gmesh_inside(geometry->geometry_parameters.p_mesh, local_position);
+  return gmesh_inside(&geometry->geometry_parameters.p_mesh_storage->mesh, local_position);
 }
 
 int sample_mesh_intersect(double *t,
@@ -2771,7 +2900,7 @@ int sample_mesh_intersect(double *t,
   // Times of intersection with a mesh component, nearest first, with the unit
   // outward normals rotated back to the frame of the master component.
   // Hits slightly before r are kept, so that the surface the ray leaves is found.
-  gmesh_struct *triangle_mesh = geometry->geometry_parameters.p_mesh;
+  gmesh_struct *triangle_mesh = &geometry->geometry_parameters.p_mesh_storage->mesh;
   gmesh_hit hits[GMESH_MAX_HITS];
   Coords position = rot_apply(geometry->transpose_rotation_matrix, coords_sub(coords_set(r[0], r[1], r[2]), geometry->center));
   Coords velocity = rot_apply(geometry->transpose_rotation_matrix, coords_set(v[0], v[1], v[2]));
@@ -4163,7 +4292,7 @@ int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct
     }
     // Secondary check with edges intersecting on faces
     if (geometry1->eShape==mesh){
-      gmesh_struct *mesh1 = geometry1->geometry_parameters.p_mesh;
+      gmesh_struct *mesh1 = &geometry1->geometry_parameters.p_mesh_storage->mesh;
       Coords vert1, vert2, vert3;
       for (i = 0; i < mesh1->ntris; i++){
         vert1 = shell_points1.elements[mesh1->tri[3*i]];
@@ -4187,7 +4316,7 @@ int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct
       }
     }
     if (geometry2->eShape==mesh){
-      gmesh_struct *mesh2 = geometry2->geometry_parameters.p_mesh;
+      gmesh_struct *mesh2 = &geometry2->geometry_parameters.p_mesh_storage->mesh;
       Coords vert1, vert2, vert3;
       for (i = 0; i < mesh2->ntris; i++){
         vert1 = shell_points2.elements[mesh2->tri[3*i]];

mcstas-comps/share/union-lib.h

@@ -225,7 +225,7 @@ union geometry_parameter_union{
     struct cylinder_storage *p_cylinder_storage;
     struct box_storage      *p_box_storage;
     struct cone_storage     *p_cone_storage;
-    gmesh_struct            *p_mesh;
+    struct mesh_storage     *p_mesh_storage;
     // add as many pointers to structs as wanted, without increasing memory footprint.
 };
 
@@ -1074,6 +1074,11 @@ double height;
 Coords direction_vector;
 };
 
+struct mesh_storage{
+gmesh_struct mesh; // triangle mesh and BVH, in the frame of the mesh component
+int display_layer; // BVH layer whose boxes mcdisplay draws, 0 to draw the triangle edges
+};
+
 
 // A number of functions below use Dot() as scalar product, replace by coords_sp define
 #define Dot(a, b) coords_sp(a, b)
@@ -1099,6 +1104,7 @@ int r_within_surroundings(Coords pos,struct geometry_struct *geometry);
 // -------------    General geometry ------------------------------------------------------------
 
 Coords point_on_circle(Coords center, Coords direction, double radius, int point_nr, int number_of_points);
+int union_mcdisplay_polyhedron(int N, struct geometry_struct **Geometries, int number_of_volumes);
 
 void points_on_circle(Coords *output, Coords center, Coords direction, double radius, int number_of_points);
 

mcstas-comps/union/Union_master.comp

@@ -42,6 +42,7 @@
 *   enable_conditionals:                 [0/1] Use conditionals with this master
 *   inherit_number_of_scattering_events: [0/1] Inherit the number of scattering events from last master
 *   weight_ratio_limit:                  [1]   Limit on the ratio between initial weight and current, ray absorbed if ratio becomes lower than limit
+*   display_mode:                        [0/1/2] Volume display in mcdisplay: 0 follows --trace, 1 forces lines, 2 polyhedra (only with --trace=2, --trace=1 always gives lines)
 * init:              [string]    Deprecated and unused. Accepted so that instruments written for McStas/McXtrace 3.8.7 and earlier, which name the Union_init component here, still compile.
 *
 * CALCULATED PARAMETERS:
@@ -63,7 +64,7 @@ SETTING PARAMETERS(int enable_refraction = 1,
                    history_limit=300000,
                    enable_conditionals=1,
                    inherit_number_of_scattering_events=0,
-                    weight_ratio_limit=1e-90, string init=""
+                    weight_ratio_limit=1e-90, int display_mode=0, string init=""
                    )
 
 NOACC
@@ -3132,6 +3133,9 @@ MCDISPLAY
   for (volume_index = 1; volume_index < number_of_volumes; volume_index++) {
 
     if (Volumes[volume_index]->geometry.visualization_on == 1) {
+      if (mcdotrace == 2 && display_mode != 1 && Geometries[volume_index]->is_mask_volume == 0 // polyhedra, masks stay wireframe
+          && union_mcdisplay_polyhedron (volume_index, Geometries, number_of_volumes))
+        continue;
       lines_to_draw_master.number_of_lines = 0;
 
       Volumes[volume_index]->geometry.mcdisplay_function (&lines_to_draw_master, volume_index, Geometries, number_of_volumes);