diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 3d6b8e0c33..779df785c9 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; } @@ -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 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;ifaceSize;i+=m->faceArray[i]+1) + for (j=1;jfaceArray[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;itranspose_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]]; diff --git a/mcstas-comps/share/union-lib.h b/mcstas-comps/share/union-lib.h index 7ab47e5520..b670e7671d 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) @@ -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); diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 496d69aa7e..f7eac88383 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/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); diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index cdc3e6ca31..aff80d5ad2 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,78 @@ 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 +243,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 +284,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 +439,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 +470,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..9505542bc1 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); @@ -1729,6 +1727,128 @@ 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 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); + } + 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;igeometry.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);