diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c
index 3d6b8e0c33..2f1c55346d 100644
--- a/mcstas-comps/share/union-lib.c
+++ b/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;
 }
@@ -2760,7 +2759,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 +2770,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 +4162,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 +4186,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]];
diff --git a/mcstas-comps/share/union-lib.h b/mcstas-comps/share/union-lib.h
index 7ab47e5520..d0d298e17b 100755
--- a/mcstas-comps/share/union-lib.h
+++ b/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)
diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp
index cdc3e6ca31..903150111a 100644
--- a/mcstas-comps/union/Union_mesh.comp
+++ b/mcstas-comps/union/Union_mesh.comp
@@ -67,6 +67,7 @@
 * focus_xw:      [m]   horiz. dimension of a rectangular area
 * focus_xh:      [m]   vert. dimension of a rectangular area
 * focus_r:       [m]   focusing on circle with this radius
+* bvh_display_layers: [1] mcdisplay draws the bounding boxes of this layer of the mesh BVH (the root is layer 1) when the BVH has more layers than this, otherwise the triangle edges. Set to 0 to always draw the triangle edges.
 * coordinate_scale [1e-3] Multiplier that the vertex positions are multiplied by. Set to 1 for input coordinates in meters. Set to 1e-3 for input coordinates in mm.
 * all_surfaces:       [string] Comma seperated string of Union surface definitions defined prior to this component
 * cut_surface:   [string] Comma seperated string of Union surface definitions defined prior to this component, applied to all internal cuts
@@ -100,7 +101,8 @@ SETTING PARAMETERS(string filename = 0,
    string mask_string=0,
    string mask_setting=0,
    number_of_activations=1,
-    coordinate_scale = 1e-3, string init=""
+    coordinate_scale = 1e-3, string init="",
+    int bvh_display_layers = 12
    )
 
 
@@ -131,10 +133,79 @@ SHARE
     return 0;
   }
 
+  int
+  mesh_bvh_layers (gmesh_struct* triangle_mesh) {
+    // Number of layers of the BVH, the root alone is one layer
+    int node_stack[GMESH_BVH_STACK + 2];
+    int layer_stack[GMESH_BVH_STACK + 2];
+    int top = 0, layers = 0;
+    node_stack[0] = 0;
+    layer_stack[0] = 1;
+    top = 1;
+    while (top > 0) {
+      top--;
+      gmesh_node* node = &triangle_mesh->node[node_stack[top]];
+      int layer = layer_stack[top];
+      if (layer > layers)
+        layers = layer;
+      if (node->count == 0) {
+        node_stack[top] = node->first;
+        layer_stack[top++] = layer + 1;
+        node_stack[top] = node->first + 1;
+        layer_stack[top++] = layer + 1;
+      }
+    }
+    return layers;
+  }
+
+  void
+  mcdisplay_mesh_bvh_boxes (struct lines_to_draw* lines_to_draw_output, int index, struct geometry_struct** Geometries, int number_of_volumes) {
+    // Draws the bounding box of every BVH node in the display layer of the mesh storage, and of the leaves above it,
+    // in the frame of the master component
+    struct geometry_struct* geometry = Geometries[index];
+    gmesh_struct* triangle_mesh = &geometry->geometry_parameters.p_mesh_storage->mesh;
+    int node_stack[GMESH_BVH_STACK + 2];
+    int layer_stack[GMESH_BVH_STACK + 2];
+    // corner bits (x, y, z) of the 12 box edges
+    const int box_edges[12][2] = { { 0, 1 }, { 2, 3 }, { 4, 5 }, { 6, 7 }, { 0, 2 }, { 1, 3 }, { 4, 6 }, { 5, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 } };
+    int top = 1;
+    int corner, edge;
+    struct lines_to_draw lines_to_draw_temp;
+    node_stack[0] = 0;
+    layer_stack[0] = 1;
+    while (top > 0) {
+      top--;
+      gmesh_node* node = &triangle_mesh->node[node_stack[top]];
+      int layer = layer_stack[top];
+      if (node->count == 0 && layer < geometry->geometry_parameters.p_mesh_storage->display_layer) {
+        node_stack[top] = node->first;
+        layer_stack[top++] = layer + 1;
+        node_stack[top] = node->first + 1;
+        layer_stack[top++] = layer + 1;
+        continue;
+      }
+      Coords corners[8];
+      for (corner = 0; corner < 8; corner++) {
+        Coords local = coords_set ((corner & 1) ? node->hi.x : node->lo.x, (corner & 2) ? node->hi.y : node->lo.y, (corner & 4) ? node->hi.z : node->lo.z);
+        corners[corner] = coords_add (rot_apply (geometry->rotation_matrix, local), geometry->center);
+      }
+      for (edge = 0; edge < 12; edge++) {
+        lines_to_draw_temp
+            = draw_line_with_highest_priority (corners[box_edges[edge][0]], corners[box_edges[edge][1]], index, Geometries, number_of_volumes, 100);
+        merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp);
+      }
+    }
+  }
+
   void
   mcdisplay_mesh_function (struct lines_to_draw* lines_to_draw_output, int index, struct geometry_struct** Geometries, int number_of_volumes) {
-    // Draws each edge of the polygons in the file once, in the frame of the master component
-    gmesh_struct* triangle_mesh = Geometries[index]->geometry_parameters.p_mesh;
+    // Draws each edge of the polygons in the file once, in the frame of the master component,
+    // or the BVH boxes of the display layer when that is set
+    if (Geometries[index]->geometry_parameters.p_mesh_storage->display_layer > 0) {
+      mcdisplay_mesh_bvh_boxes (lines_to_draw_output, index, Geometries, number_of_volumes);
+      return;
+    }
+    gmesh_struct* triangle_mesh = &Geometries[index]->geometry_parameters.p_mesh_storage->mesh;
     struct pointer_to_1d_coords_list points_struct = Geometries[index]->shell_points (Geometries[index], 1);
     Coords* points = points_struct.elements;
     long* edges = malloc (2 * triangle_mesh->faceSize * sizeof (long));
@@ -173,7 +244,7 @@ SHARE
   struct pointer_to_1d_coords_list
   allocate_shell_points (struct geometry_struct* geometry, int max_number_of_points) {
     // Function that returns all vertex points in the frame of the master component
-    gmesh_struct* triangle_mesh = geometry->geometry_parameters.p_mesh;
+    gmesh_struct* triangle_mesh = &geometry->geometry_parameters.p_mesh_storage->mesh;
     struct pointer_to_1d_coords_list mesh_shell_array;
 
     mesh_shell_array.elements = malloc (triangle_mesh->nverts * sizeof (Coords));
@@ -214,7 +285,7 @@ DECLARE
   struct Volume_struct mesh_vol;
 
   // Mesh storage contains the triangle mesh read with geometry-mesh-lib
-  gmesh_struct mesh_data;
+  struct mesh_storage this_mesh_storage;
 
   // Structs for surface effects
   struct surface_stack_struct surface_stack;
@@ -369,19 +440,19 @@ INITIALIZE
   //                       READ FILE
   //==============================================================================
   // Read the STL/OFF/PLY file in the frame of this component, scaled by coordinate_scale
-  if (!gmesh_init (filename, coordinate_scale, GMESH_PROPS_IGNORE, &mesh_data)) {
+  if (!gmesh_init (filename, coordinate_scale, GMESH_PROPS_IGNORE, &this_mesh_storage.mesh)) {
     printf ("\nERROR IN %s: could not read the geometry file %s\n", NAME_CURRENT_COMP, filename);
     exit (1);
   }
-  if (!mesh_data.closed) {
+  if (!this_mesh_storage.mesh.closed) {
     printf ("\n\nERROR: Input mesh %s of %s is not closed: %d boundary and %d non-manifold edges."
             "\nUnion_mesh requires a closed (watertight) mesh.\n\n",
-            filename, NAME_CURRENT_COMP, mesh_data.nboundary, mesh_data.nnonmanifold);
+            filename, NAME_CURRENT_COMP, this_mesh_storage.mesh.nboundary, this_mesh_storage.mesh.nnonmanifold);
     exit (1);
   }
 
   // surface_index of a mesh hit is the triangle index, all triangles share one surface stack
-  mesh_vol.geometry.number_of_faces = mesh_data.ntris;
+  mesh_vol.geometry.number_of_faces = this_mesh_storage.mesh.ntris;
   mesh_vol.geometry.surface_stack_for_each_face = malloc (mesh_vol.geometry.number_of_faces * sizeof (struct surface_stack_struct*));
   if (!mesh_vol.geometry.surface_stack_for_each_face) {
     printf ("\nERROR IN %s: could not allocate surface stacks\n", NAME_CURRENT_COMP);
@@ -400,9 +471,19 @@ INITIALIZE
   mesh_vol.geometry.eShape = mesh;
   mesh_vol.geometry.priority_value = priority;
   mesh_vol.geometry.visualization_on = visualize;
+  this_mesh_storage.display_layer = 0;
+  if (bvh_display_layers > 0) {
+    int layers = mesh_bvh_layers (&this_mesh_storage.mesh);
+    if (layers > bvh_display_layers) {
+      this_mesh_storage.display_layer = bvh_display_layers;
+      if (mcdotrace)
+        MPI_MASTER (printf ("Union_mesh %s: mcdisplay draws the BVH boxes of layer %d of %d instead of the %d triangles\n", NAME_CURRENT_COMP, bvh_display_layers,
+                            layers, this_mesh_storage.mesh.ntris);)
+    }
+  }
   mesh_vol.geometry.geometry_p_interact = p_interact;
   mesh_vol.geometry.within_function = &r_within_mesh;
-  mesh_vol.geometry.geometry_parameters.p_mesh = &mesh_data;
+  mesh_vol.geometry.geometry_parameters.p_mesh_storage = &this_mesh_storage;
   mesh_vol.geometry.intersect_function = &sample_mesh_intersect;
   mesh_vol.geometry.mcdisplay_function = &mcdisplay_mesh_function;
   mesh_vol.geometry.shell_points = &allocate_shell_points;
diff --git a/mcxtrace-comps/share/union-lib.c b/mcxtrace-comps/share/union-lib.c
index 04c4554ada..edf44fc676 100644
--- a/mcxtrace-comps/share/union-lib.c
+++ b/mcxtrace-comps/share/union-lib.c
@@ -1293,39 +1293,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) {
-      fprintf(stderr,"Failure allocating list in Union function merge_lines_to_draw 1 - Exit!\n");
-      exit(EXIT_FAILURE);
-    }
-    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) {
@@ -1473,6 +1470,7 @@ struct lines_to_draw draw_line_with_highest_priority(Coords position1,Coords pos
 
         }
     }
+    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!
 
@@ -1534,9 +1532,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);