From 24c751f03aa5dbe775d68df139c65806ff6f7cff Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 28 Apr 2026 15:14:43 +0200 Subject: [PATCH 01/12] refactor moeller_Truborer intersections --- mcstas-comps/share/union-lib.c | 541 ++++++++++++----------------- mcstas-comps/union/Union_mesh.comp | 21 +- 2 files changed, 236 insertions(+), 326 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 755352aef..7253411ad 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -2598,6 +2598,9 @@ double *normal_y; double *normal_z; Coords direction_vector; Coords Bounding_Box_Center; +Coords *vertices; +int **facets; +double Bounding_Box_Extremes[6]; double Bounding_Box_Radius; }; @@ -3698,223 +3701,125 @@ This function was created by Martin Olsen at NBI on september 20, 2018. return 0; }; -int r_within_mesh(Coords pos,struct geometry_struct *geometry) { - //TODO: Make a single intersection algorithm that all the three mesh intersections - // Can use -// Unpack parameters - Coords center = geometry->center; - int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; - double *normal_x = geometry->geometry_parameters.p_mesh_storage->normal_x; - double *normal_y = geometry->geometry_parameters.p_mesh_storage->normal_y; - double *normal_z = geometry->geometry_parameters.p_mesh_storage->normal_z; - double *v1_x = geometry->geometry_parameters.p_mesh_storage->v1_x; - double *v1_y = geometry->geometry_parameters.p_mesh_storage->v1_y; - double *v1_z = geometry->geometry_parameters.p_mesh_storage->v1_z; - double *v2_x = geometry->geometry_parameters.p_mesh_storage->v2_x; - double *v2_y = geometry->geometry_parameters.p_mesh_storage->v2_y; - double *v2_z = geometry->geometry_parameters.p_mesh_storage->v2_z; - double *v3_x = geometry->geometry_parameters.p_mesh_storage->v3_x; - double *v3_y = geometry->geometry_parameters.p_mesh_storage->v3_y; - double *v3_z = geometry->geometry_parameters.p_mesh_storage->v3_z; - - double x_new,y_new,z_new; - - // Coordinate transformation - x_new = pos.x - geometry->center.x; - y_new = pos.y - geometry->center.y; - z_new = pos.z - geometry->center.z; - - Coords coordinates = coords_set(x_new,y_new,z_new); - Coords rotated_coordinates; - - rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); +struct Moeller_Trumbore{ + Coords v1; + Coords v2; + Coords v3; + Coords edge1; + Coords edge2; + Coords h; + Coords s; + Coords q; + Coords rotated_coordinates; + Coords rotated_velocity; + double a; + double f; + double V; + double u; +}; - - int verbal = 0; - // Generate unit direction vector along center axis of meshs - - // Start with vector that points along the mesh in the simple frame, and rotate to global - Coords simple_vector = coords_set(0,1,0); - Coords test_vector = coords_set(0,1,0); - // Check intersections with every single facet: - int iter =0; - int counter=0; int neg_counter=0; - Coords edge1,edge2,h,s,q,tmp,intersect_pos; - double UNION_EPSILON = 1e-27; - double this_facet_t; - double a,f,u,V; - //////printf("\n RWITHIN TEST 1ste"); - for (iter = 0 ; iter < n_facets ; iter++){ - // Intersection with face plane (Möller–Trumbore) - edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); - edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); - - vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); - - a = Dot(edge1,h); - if (a > -UNION_EPSILON && a < UNION_EPSILON){ - //////printf("\n UNION_EPSILON fail"); - } else{ - f = 1.0/a; - s = coords_sub(rotated_coordinates, coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); - u = f * (Dot(s,h)); - if (u < 0.0 || u > 1.0){ - }else{ - //q = vec_prod(s,edge1); - vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); - V = f * Dot(test_vector,q); - if (V < 0.0 || u + V > 1.0){ - } else { - // At this stage we can compute t to find out where the intersection point is on the line. - if (f* Dot(q,edge2) > 0){ - counter++; - } else { - neg_counter++; +double Moeller_Trumbore_intersection(struct Moeller_Trumbore* intersect) { + // Function to perform a Moeller Trumbore intersection between a mesh facet + // and an arbitrary vector + intersect->edge1 = coords_sub(intersect->v2, intersect->v1); + intersect->edge2 = coords_sub(intersect->v3, intersect->v1); - } - if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ - } else { //printf("\n [%f %f %f] Failed due to being close to surface, E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); - } - - - - } - } - } - } - int C1 = counter; - - int maxC; int sameNr =0; - if (counter % 2 == neg_counter % 2){ - maxC = counter; - sameNr = 1; - } else { - //printf("\n not the same intersection numbers (%i , %i)",counter,neg_counter); - maxC = counter; - sameNr = 0; - } + intersect->h = coords_xp(intersect->rotated_velocity, intersect->edge2); + intersect->a = Dot(intersect->edge1, intersect->h); - if (sameNr == 0){ - test_vector = coords_set(0,0,1); - iter =0; - counter=0; - for (iter = 0 ; iter < n_facets ; iter++){ - // Intersection with face plane (Möller–Trumbore) - edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); - edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); - - vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); - - a = Dot(edge1,h); - if (a > -UNION_EPSILON && a < UNION_EPSILON){ - //////printf("\n UNION_EPSILON fail"); - } else{ - f = 1.0/a; - s = coords_sub(rotated_coordinates , coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); - u = f * (Dot(s,h)); - if (u < 0.0 || u > 1.0){ - }else{ - //q = vec_prod(s,edge1); - vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); - V = f * Dot(test_vector,q); - if (V < 0.0 || u + V > 1.0){ - } else { - // At this stage we can compute t to find out where the intersection point is on the line. - - if (f* Dot(q,edge2) > 0){ - counter++; - - } else { - neg_counter++; + intersect->f = 1.0/intersect->a; + intersect->s = coords_sub(intersect->rotated_coordinates, intersect->v1); + intersect->u = intersect->f * (Dot(intersect->s,intersect->h)); + if (intersect->u < 0.0 || intersect->u > 1.0){ + return 0; + } else { + intersect->q = coords_xp(intersect->s, intersect->edge1); + intersect->V = intersect->f * Dot(intersect->rotated_velocity,intersect->q); - } - if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ - } else { printf("\n [%f %f %f] Failed due to being close to surface (2. iteration), E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); - } - - - - } - } - } - } - } - - if (counter % 2 == neg_counter % 2){ - maxC = counter; - sameNr = 1; + if (intersect->V < 0.0 || intersect->u + intersect->V > 1.0){ + return 0; } else { - printf("\n not the same intersection numbers (%i , %i) second iteration",counter,neg_counter); - maxC = counter; - sameNr = 0; - } + // At this stage we can compute t to find out where the intersection point is on the line. + return intersect->f * Dot(intersect->q,intersect->edge2); - if (sameNr == 0){ - test_vector = coords_set(1,0,0); - iter =0; - counter=0; - for (iter = 0 ; iter < n_facets ; iter++){ - // Intersection with face plane (Möller–Trumbore) - edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); - edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); - - vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); - - a = Dot(edge1,h); - if (a > -UNION_EPSILON && a < UNION_EPSILON){ - //////printf("\n UNION_EPSILON fail"); - } else{ - f = 1.0/a; - s = coords_sub(rotated_coordinates , coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); - u = f * (Dot(s,h)); - if (u < 0.0 || u > 1.0){ - //////printf("\n Nope 1"); - }else{ - //q = vec_prod(s,edge1); - vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); - V = f * Dot(test_vector,q); - if (V < 0.0 || u + V > 1.0){ - } else { - // At this stage we can compute t to find out where the intersection point is on the line. - - if (f* Dot(q,edge2) > 0){ - counter++; + } + } - } else { - neg_counter++; + return 0; +} - } - if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ - } else { printf("\n [%f %f %f] Failed due to being close to surface (3. iteration), E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); - } - - - - } - } - } - } +int r_within_mesh(Coords pos,struct geometry_struct *geometry) { + // r_within_mesh uses a ray casting technique to determine whether or not + // a position is within the mesh. + // It chooses a random velocity, and if the number of intersections + // is uneven, the position is inside the mesh. + // Since this can be numerically unstable, it performs this ray casting 3 times + // and then allows the majority of results to decide whether inside or outside + + Coords center = geometry->center; + double x_new,y_new,z_new; + + // Coordinate transformation + x_new = pos.x - geometry->center.x; + y_new = pos.y - geometry->center.y; + z_new = pos.z - geometry->center.z; + + Coords coordinates = coords_set(x_new,y_new,z_new); + Coords rotated_coordinates; + + rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); + + // Check intersections with every single facet: + // First do allocations: + int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; + double possible_t; + int iter =0; + int counter; + struct Moeller_Trumbore intersect_transport; + double *t_intersect=malloc(n_facets*sizeof(double)); + int *facet_index = malloc(n_facets*sizeof(int)); + if (!t_intersect || !facet_index) { + fprintf(stderr,"Failure allocating list in Union function sample_mesh_intersect - Exit!\n"); + exit(EXIT_FAILURE); + } + Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; + int **facets = geometry->geometry_parameters.p_mesh_storage->facets; + // Then loop over every facet + intersect_transport.rotated_coordinates = rotated_coordinates; + int inside_vote = 0, outside_vote = 0; + + for (int j = 0; j <3; j++){ + counter = 0; + Coords test_vector = coords_set(rand01(),rand01(),rand01()); + intersect_transport.rotated_velocity = test_vector; + for (iter = 0 ; iter < n_facets ; iter++){ + intersect_transport.v1 = verts[facets[iter][0]]; + intersect_transport.v2 = verts[facets[iter][1]]; + intersect_transport.v3 = verts[facets[iter][2]]; + possible_t = Moeller_Trumbore_intersection(&intersect_transport); + if (possible_t != 0){ + counter++; + } } - - - if (counter % 2 == neg_counter % 2){ - maxC = counter; + if (counter%2==1){ + inside_vote++; } else { - return 0; + outside_vote++; } - if ( maxC % 2 == 0) { - return 0; - }else{ - return 1; - + if (inside_vote == 2 || outside_vote == 2){ + break; } - + } + if (inside_vote == 2){ + return 1; + } else { return 0; - }; - + } +} + // Type for holding intersection and normal typedef struct { @@ -3932,114 +3837,97 @@ int compare_intersections(const void *a, const void *b) { return 0; } + + + + int sample_mesh_intersect(double *t, double *nx, double *ny, double*nz, int *surface_index, int *num_solutions,double *r,double *v, struct geometry_struct *geometry) { - - int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; - double *normal_x = geometry->geometry_parameters.p_mesh_storage->normal_x; - double *normal_y = geometry->geometry_parameters.p_mesh_storage->normal_y; - double *normal_z = geometry->geometry_parameters.p_mesh_storage->normal_z; - double *v1_x = geometry->geometry_parameters.p_mesh_storage->v1_x; - double *v1_y = geometry->geometry_parameters.p_mesh_storage->v1_y; - double *v1_z = geometry->geometry_parameters.p_mesh_storage->v1_z; - double *v2_x = geometry->geometry_parameters.p_mesh_storage->v2_x; - double *v2_y = geometry->geometry_parameters.p_mesh_storage->v2_y; - double *v2_z= geometry->geometry_parameters.p_mesh_storage->v2_z; - double *v3_x = geometry->geometry_parameters.p_mesh_storage->v3_x; - double *v3_y = geometry->geometry_parameters.p_mesh_storage->v3_y; - double *v3_z = geometry->geometry_parameters.p_mesh_storage->v3_z; - Coords Bounding_Box_Center = geometry->geometry_parameters.p_mesh_storage->Bounding_Box_Center; - double Bounding_Box_Radius = geometry->geometry_parameters.p_mesh_storage->Bounding_Box_Radius; - int i; - - double x_new,y_new,z_new; - // Coordinate transformation - x_new = r[0] - geometry->center.x; - y_new = r[1] - geometry->center.y; - z_new = r[2] - geometry->center.z; - - double x_bb,y_bb,z_bb; - x_bb = r[0] - Bounding_Box_Center.x - geometry->center.x; - y_bb = r[1] - Bounding_Box_Center.y - geometry->center.y; - z_bb = r[2] - Bounding_Box_Center.z - geometry->center.z; - - Coords coordinates = coords_set(x_new,y_new,z_new); - Coords rotated_coordinates; - - // Rotate the position of the neutron around the center of the mesh - rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); - - Coords bounding_box_coordinates = coords_set(x_bb, y_bb, z_bb); - Coords bounding_box_rotated_coordinates; - - bounding_box_rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,bounding_box_coordinates); - - Coords velocity = coords_set(v[0],v[1],v[2]); - Coords rotated_velocity; - - // Rotate the position of the neutron around the center of the mesh - rotated_velocity = rot_apply(geometry->transpose_rotation_matrix,velocity); - - int output = 0; - double tmpres[2]; - // Test intersection with bounding sphere - if ((output = sphere_intersect(&tmpres[0],&tmpres[1], - bounding_box_rotated_coordinates.x, - bounding_box_rotated_coordinates.y, - bounding_box_rotated_coordinates.z, - rotated_velocity.x, - rotated_velocity.y, - rotated_velocity.z, - Bounding_Box_Radius)) == 0) { - t[0] = -1; - t[1] = -1; - *num_solutions = 0; - return 0; - } - // Check intersections with every single facet: - int iter =0; - int counter=0; - Coords edge1,edge2,h,s,q,tmp,intersect_pos; - double UNION_EPSILON = 0.0000001; - double this_facet_t; - double a,f,u,V; - double *t_intersect=malloc(n_facets*sizeof(double)); - int *facet_index = malloc(n_facets*sizeof(int)); - if (!t_intersect || !facet_index) { - fprintf(stderr,"Failure allocating list in Union function sample_mesh_intersect - Exit!\n"); - exit(EXIT_FAILURE); - } + // Algorithm for finding the times of intersection with a mesh component. + // First, check if the neutron intersects the bounding sphere of the mesh. + // Then, if yes, loop over every single facet, and see if the neutron + // intersects with it. + Coords Bounding_Box_Center = geometry->geometry_parameters.p_mesh_storage->Bounding_Box_Center; + double Bounding_Box_Radius = geometry->geometry_parameters.p_mesh_storage->Bounding_Box_Radius; + int i; + + double x_new,y_new,z_new; + // Coordinate transformation + x_new = r[0] - geometry->center.x; + y_new = r[1] - geometry->center.y; + z_new = r[2] - geometry->center.z; + + double x_bb,y_bb,z_bb; + x_bb = r[0] - Bounding_Box_Center.x - geometry->center.x; + y_bb = r[1] - Bounding_Box_Center.y - geometry->center.y; + z_bb = r[2] - Bounding_Box_Center.z - geometry->center.z; + + Coords coordinates = coords_set(x_new,y_new,z_new); + Coords rotated_coordinates; + + // Rotate the position of the neutron around the center of the mesh + rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); + + Coords bounding_box_coordinates = coords_set(x_bb, y_bb, z_bb); + Coords bounding_box_rotated_coordinates; + + bounding_box_rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,bounding_box_coordinates); + + Coords velocity = coords_set(v[0],v[1],v[2]); + Coords rotated_velocity; + + // Rotate the position of the neutron around the center of the mesh + rotated_velocity = rot_apply(geometry->transpose_rotation_matrix,velocity); + + int output = 0; + double tmpres[2]; + // Test intersection with bounding sphere + if ((output = sphere_intersect(&tmpres[0],&tmpres[1], + bounding_box_rotated_coordinates.x, + bounding_box_rotated_coordinates.y, + bounding_box_rotated_coordinates.z, + rotated_velocity.x, + rotated_velocity.y, + rotated_velocity.z, + Bounding_Box_Radius)) == 0) { + t[0] = -1; + t[1] = -1; *num_solutions = 0; - for (iter = 0 ; iter < n_facets ; iter++){ - // Intersection with face plane (Möller–Trumbore) - edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); - edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); - vec_prod(h.x,h.y,h.z,rotated_velocity.x,rotated_velocity.y,rotated_velocity.z,edge2.x,edge2.y,edge2.z); - a = Dot(edge1,h); - //if (a > -UNION_EPSILON && a < UNION_EPSILON){ - ////printf("\n UNION_EPSILON fail"); - //} - f = 1.0/a; - s = coords_sub(rotated_coordinates, coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); - u = f * (Dot(s,h)); - if (u < 0.0 || u > 1.0){ - } else { - //q = vec_prod(s,edge1); - vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); - V = f * Dot(rotated_velocity,q); - if (V < 0.0 || u + V > 1.0){ - } else { - // At this stage we can compute t to find out where the intersection point is on the line. - t_intersect[counter] = f* Dot(q,edge2); - facet_index[counter] = iter; - //printf("\nIntersects at time: t= %f\n",t_intersect[counter] ); - counter++; - } - } + return 0; + } + // Check intersections with every single facet: + // First do allocations: + int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; + double possible_t; + int iter =0; + int counter=0; + struct Moeller_Trumbore intersect_transport; + double *t_intersect=malloc(n_facets*sizeof(double)); + int *facet_index = malloc(n_facets*sizeof(int)); + if (!t_intersect || !facet_index) { + fprintf(stderr,"Failure allocating list in Union function sample_mesh_intersect - Exit!\n"); + exit(EXIT_FAILURE); + } + Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; + int **facets = geometry->geometry_parameters.p_mesh_storage->facets; + // Then loop over every facet + intersect_transport.rotated_velocity = rotated_velocity; + intersect_transport.rotated_coordinates = rotated_coordinates; + *num_solutions = 0; + for (iter = 0 ; iter < n_facets ; iter++){ + intersect_transport.v1 = verts[facets[iter][0]]; + intersect_transport.v2 = verts[facets[iter][1]]; + intersect_transport.v3 = verts[facets[iter][2]]; + possible_t = Moeller_Trumbore_intersection(&intersect_transport); + if (possible_t != 0){ + t_intersect[counter] = possible_t; + facet_index[counter] = iter; + counter++; } + } *num_solutions = counter; @@ -4051,36 +3939,39 @@ int sample_mesh_intersect(double *t, } // Move times and normal's into structs to be sorted - Intersection *hits = malloc(*num_solutions * sizeof(Intersection)); - if (!hits) { - fprintf(stderr,"Failure allocating Intersection list struct in Union function sample_mesh_intersect - Exit!\n"); - exit(EXIT_FAILURE); - } - - for (iter=0; iter < *num_solutions; iter++){ - hits[iter].t = t_intersect[iter];; - hits[iter].nx = normal_x[facet_index[iter]]; - hits[iter].ny = normal_y[facet_index[iter]]; - hits[iter].nz = normal_z[facet_index[iter]]; - hits[iter].surface_index = 0; - } + Intersection *hits = malloc(*num_solutions * sizeof(Intersection)); + if (!hits) { + fprintf(stderr,"Failure allocating Intersection list struct in Union function sample_mesh_intersect - Exit!\n"); + exit(EXIT_FAILURE); + } - // Sort structs according to time - qsort(hits, *num_solutions, sizeof(Intersection), compare_intersections); + double *normal_x = geometry->geometry_parameters.p_mesh_storage->normal_x; + double *normal_y = geometry->geometry_parameters.p_mesh_storage->normal_y; + double *normal_z = geometry->geometry_parameters.p_mesh_storage->normal_z; + for (iter=0; iter < *num_solutions; iter++){ + hits[iter].t = t_intersect[iter];; + hits[iter].nx = normal_x[facet_index[iter]]; + hits[iter].ny = normal_y[facet_index[iter]]; + hits[iter].nz = normal_z[facet_index[iter]]; + hits[iter].surface_index = 0; + } + + // Sort structs according to time + qsort(hits, *num_solutions, sizeof(Intersection), compare_intersections); // Place the solutions into the pointers given in the function parameters for return for (int i = 0; i < *num_solutions; i++) { - t[i] = hits[i].t; - nx[i] = hits[i].nx; - ny[i] = hits[i].ny; - nz[i] = hits[i].nz; - surface_index[i] = hits[i].surface_index; + t[i] = hits[i].t; + nx[i] = hits[i].nx; + ny[i] = hits[i].ny; + nz[i] = hits[i].nz; + surface_index[i] = hits[i].surface_index; } - free(facet_index); - free(t_intersect); + free(facet_index); + free(t_intersect); free(hits); - return 1; + return 1; }; diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index 04b450283..e552e2188 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -639,7 +639,7 @@ SHARE for (i = 0; i < number_of_points_in_array; i++) { // Transpose and rotate the points such that they are in the right coordinate system mesh_shell_array.elements[i] = coords_add (mesh_shell_array.elements[i], geometry->center); - mesh_shell_array.elements[i] = rot_apply (geometry->rotation_matrix, mesh_shell_array.elements[i]); + mesh_shell_array.elements[i] = rot_apply (geometry->transpose_rotation_matrix, mesh_shell_array.elements[i]); } return mesh_shell_array; } @@ -894,6 +894,7 @@ INITIALIZE for (int i = 0; i < n_verts; i++) { verts[i] = coords_scale (verts[i], coordinate_scale); } + // Make lists that will be used to calculate the number of edges int** vert_pairs = (int**)malloc (sizeof (int*) * 3 * n_faces); @@ -971,6 +972,22 @@ INITIALIZE if (test_rad > radius) { radius = test_rad; } + if (i == 0){ + mesh_data.Bounding_Box_Extremes[0] = verts[i].x; + mesh_data.Bounding_Box_Extremes[1] = verts[i].x; + mesh_data.Bounding_Box_Extremes[2] = verts[i].y; + mesh_data.Bounding_Box_Extremes[3] = verts[i].y; + mesh_data.Bounding_Box_Extremes[4] = verts[i].z; + mesh_data.Bounding_Box_Extremes[5] = verts[i].z; + } else { + mesh_data.Bounding_Box_Extremes[0] = verts[i].x > mesh_data.Bounding_Box_Extremes[0] ? verts[i].x : mesh_data.Bounding_Box_Extremes[0]; + mesh_data.Bounding_Box_Extremes[1] = verts[i].x < mesh_data.Bounding_Box_Extremes[1] ? verts[i].x : mesh_data.Bounding_Box_Extremes[1]; + mesh_data.Bounding_Box_Extremes[2] = verts[i].y > mesh_data.Bounding_Box_Extremes[2] ? verts[i].y : mesh_data.Bounding_Box_Extremes[2]; + mesh_data.Bounding_Box_Extremes[3] = verts[i].y < mesh_data.Bounding_Box_Extremes[3] ? verts[i].y : mesh_data.Bounding_Box_Extremes[3]; + mesh_data.Bounding_Box_Extremes[4] = verts[i].z > mesh_data.Bounding_Box_Extremes[4] ? verts[i].z : mesh_data.Bounding_Box_Extremes[4]; + mesh_data.Bounding_Box_Extremes[5] = verts[i].z < mesh_data.Bounding_Box_Extremes[5] ? verts[i].z : mesh_data.Bounding_Box_Extremes[5]; + + } } mesh_data.Bounding_Box_Radius = radius; @@ -990,6 +1007,8 @@ INITIALIZE mesh_data.normal_x = (double*)malloc (n_faces * sizeof (double)); mesh_data.normal_y = (double*)malloc (n_faces * sizeof (double)); mesh_data.normal_z = (double*)malloc (n_faces * sizeof (double)); + mesh_data.facets = faces; + mesh_data.vertices = verts; for (int i = 0; i < n_faces; i++) { vert1 = verts[faces[i][0]]; vert2 = verts[faces[i][1]]; From 1dcc38a21803bac59228a3d69b05d91946f38101 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 6 May 2026 13:32:22 +0200 Subject: [PATCH 02/12] Add in prefactored r_within mesh and implement edge intersection for mesh_overlaps_mesh --- mcstas-comps/share/union-lib.c | 372 +++++++++++++++++++++++------ mcstas-comps/union/Union_mesh.comp | 2 +- 2 files changed, 304 insertions(+), 70 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 7253411ad..faa961fad 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -3733,13 +3733,13 @@ double Moeller_Trumbore_intersection(struct Moeller_Trumbore* intersect) { intersect->s = coords_sub(intersect->rotated_coordinates, intersect->v1); intersect->u = intersect->f * (Dot(intersect->s,intersect->h)); if (intersect->u < 0.0 || intersect->u > 1.0){ - return 0; + return -1; } else { intersect->q = coords_xp(intersect->s, intersect->edge1); intersect->V = intersect->f * Dot(intersect->rotated_velocity,intersect->q); if (intersect->V < 0.0 || intersect->u + intersect->V > 1.0){ - return 0; + return -1; } else { // At this stage we can compute t to find out where the intersection point is on the line. return intersect->f * Dot(intersect->q,intersect->edge2); @@ -3752,73 +3752,288 @@ double Moeller_Trumbore_intersection(struct Moeller_Trumbore* intersect) { int r_within_mesh(Coords pos,struct geometry_struct *geometry) { - // r_within_mesh uses a ray casting technique to determine whether or not - // a position is within the mesh. - // It chooses a random velocity, and if the number of intersections - // is uneven, the position is inside the mesh. - // Since this can be numerically unstable, it performs this ray casting 3 times - // and then allows the majority of results to decide whether inside or outside + //TODO: Make a single intersection algorithm that all the three mesh intersections + // Can use +// Unpack parameters - Coords center = geometry->center; - double x_new,y_new,z_new; - - // Coordinate transformation - x_new = pos.x - geometry->center.x; - y_new = pos.y - geometry->center.y; - z_new = pos.z - geometry->center.z; - - Coords coordinates = coords_set(x_new,y_new,z_new); - Coords rotated_coordinates; - - rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); + Coords center = geometry->center; + int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; + double *normal_x = geometry->geometry_parameters.p_mesh_storage->normal_x; + double *normal_y = geometry->geometry_parameters.p_mesh_storage->normal_y; + double *normal_z = geometry->geometry_parameters.p_mesh_storage->normal_z; + double *v1_x = geometry->geometry_parameters.p_mesh_storage->v1_x; + double *v1_y = geometry->geometry_parameters.p_mesh_storage->v1_y; + double *v1_z = geometry->geometry_parameters.p_mesh_storage->v1_z; + double *v2_x = geometry->geometry_parameters.p_mesh_storage->v2_x; + double *v2_y = geometry->geometry_parameters.p_mesh_storage->v2_y; + double *v2_z = geometry->geometry_parameters.p_mesh_storage->v2_z; + double *v3_x = geometry->geometry_parameters.p_mesh_storage->v3_x; + double *v3_y = geometry->geometry_parameters.p_mesh_storage->v3_y; + double *v3_z = geometry->geometry_parameters.p_mesh_storage->v3_z; - // Check intersections with every single facet: - // First do allocations: - int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; - double possible_t; - int iter =0; - int counter; - struct Moeller_Trumbore intersect_transport; - double *t_intersect=malloc(n_facets*sizeof(double)); - int *facet_index = malloc(n_facets*sizeof(int)); - if (!t_intersect || !facet_index) { - fprintf(stderr,"Failure allocating list in Union function sample_mesh_intersect - Exit!\n"); - exit(EXIT_FAILURE); - } - Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; - int **facets = geometry->geometry_parameters.p_mesh_storage->facets; - // Then loop over every facet - intersect_transport.rotated_coordinates = rotated_coordinates; - int inside_vote = 0, outside_vote = 0; + double x_new,y_new,z_new; - for (int j = 0; j <3; j++){ - counter = 0; - Coords test_vector = coords_set(rand01(),rand01(),rand01()); - intersect_transport.rotated_velocity = test_vector; + // Coordinate transformation + x_new = pos.x - geometry->center.x; + y_new = pos.y - geometry->center.y; + z_new = pos.z - geometry->center.z; + + Coords coordinates = coords_set(x_new,y_new,z_new); + Coords rotated_coordinates; + + rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); + + + int verbal = 0; + // Generate unit direction vector along center axis of meshs + + // Start with vector that points along the mesh in the simple frame, and rotate to global + Coords simple_vector = coords_set(0,1,0); + Coords test_vector = coords_set(0,1,0); + // Check intersections with every single facet: + int iter =0; + int counter=0; int neg_counter=0; + Coords edge1,edge2,h,s,q,tmp,intersect_pos; + double UNION_EPSILON = 1e-27; + double this_facet_t; + double a,f,u,V; + //////printf("\n RWITHIN TEST 1ste"); for (iter = 0 ; iter < n_facets ; iter++){ - intersect_transport.v1 = verts[facets[iter][0]]; - intersect_transport.v2 = verts[facets[iter][1]]; - intersect_transport.v3 = verts[facets[iter][2]]; - possible_t = Moeller_Trumbore_intersection(&intersect_transport); - if (possible_t != 0){ - counter++; - } + // Intersection with face plane (Möller–Trumbore) + edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); + edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); + + vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); + + a = Dot(edge1,h); + if (a > -UNION_EPSILON && a < UNION_EPSILON){ + //////printf("\n UNION_EPSILON fail"); + } else{ + f = 1.0/a; + s = coords_sub(rotated_coordinates, coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); + u = f * (Dot(s,h)); + if (u < 0.0 || u > 1.0){ + }else{ + //q = vec_prod(s,edge1); + vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); + V = f * Dot(test_vector,q); + if (V < 0.0 || u + V > 1.0){ + } else { + // At this stage we can compute t to find out where the intersection point is on the line. + if (f* Dot(q,edge2) > 0){ + counter++; + + } else { + neg_counter++; + + } + if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ + } else { //printf("\n [%f %f %f] Failed due to being close to surface, E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); + } + + + + } + } + } } - if (counter%2==1){ - inside_vote++; + int C1 = counter; + + int maxC; int sameNr =0; + if (counter % 2 == neg_counter % 2){ + maxC = counter; + sameNr = 1; } else { - outside_vote++; + //printf("\n not the same intersection numbers (%i , %i)",counter,neg_counter); + maxC = counter; + sameNr = 0; } - if (inside_vote == 2 || outside_vote == 2){ - break; + + if (sameNr == 0){ + test_vector = coords_set(0,0,1); + iter =0; + counter=0; + for (iter = 0 ; iter < n_facets ; iter++){ + // Intersection with face plane (Möller–Trumbore) + edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); + edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); + + vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); + + a = Dot(edge1,h); + if (a > -UNION_EPSILON && a < UNION_EPSILON){ + //////printf("\n UNION_EPSILON fail"); + } else{ + f = 1.0/a; + s = coords_sub(rotated_coordinates , coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); + u = f * (Dot(s,h)); + if (u < 0.0 || u > 1.0){ + }else{ + //q = vec_prod(s,edge1); + vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); + V = f * Dot(test_vector,q); + if (V < 0.0 || u + V > 1.0){ + } else { + // At this stage we can compute t to find out where the intersection point is on the line. + + if (f* Dot(q,edge2) > 0){ + counter++; + + } else { + neg_counter++; + + } + if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ + } else { printf("\n [%f %f %f] Failed due to being close to surface (2. iteration), E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); + } + + + + } + } + } } - } - if (inside_vote == 2){ - return 1; - } else { + } + + if (counter % 2 == neg_counter % 2){ + maxC = counter; + sameNr = 1; + } else { + printf("\n not the same intersection numbers (%i , %i) second iteration",counter,neg_counter); + maxC = counter; + sameNr = 0; + } + + if (sameNr == 0){ + test_vector = coords_set(1,0,0); + iter =0; + counter=0; + for (iter = 0 ; iter < n_facets ; iter++){ + // Intersection with face plane (Möller–Trumbore) + edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); + edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); + + vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); + + a = Dot(edge1,h); + if (a > -UNION_EPSILON && a < UNION_EPSILON){ + //////printf("\n UNION_EPSILON fail"); + } else{ + f = 1.0/a; + s = coords_sub(rotated_coordinates , coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); + u = f * (Dot(s,h)); + if (u < 0.0 || u > 1.0){ + //////printf("\n Nope 1"); + }else{ + //q = vec_prod(s,edge1); + vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); + V = f * Dot(test_vector,q); + if (V < 0.0 || u + V > 1.0){ + } else { + // At this stage we can compute t to find out where the intersection point is on the line. + + if (f* Dot(q,edge2) > 0){ + counter++; + + } else { + neg_counter++; + + } + if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ + } else { printf("\n [%f %f %f] Failed due to being close to surface (3. iteration), E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); + } + + + + } + } + } + } + + } + + + if (counter % 2 == neg_counter % 2){ + maxC = counter; + } else { + return 0; + } + if ( maxC % 2 == 0) { + return 0; + }else{ + return 1; + + } + return 0; - } -} + }; + // r_within_mesh uses a ray casting technique to determine whether or not + // a position is within the mesh. + // It chooses a random velocity, and if the number of intersections + // is uneven, the position is inside the mesh. + // Since this can be numerically unstable, it performs this ray casting 3 times + // and then allows the majority of results to decide whether inside or outside +// +// Coords center = geometry->center; +// double x_new,y_new,z_new; +// +// // Coordinate transformation +// x_new = pos.x - geometry->center.x; +// y_new = pos.y - geometry->center.y; +// z_new = pos.z - geometry->center.z; +// +// Coords coordinates = coords_set(x_new,y_new,z_new); +// Coords rotated_coordinates; +// +// rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); +// +// // Check intersections with every single facet: +// // First do allocations: +// int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; +// double possible_t; +// int iter =0; +// int counter; +// struct Moeller_Trumbore intersect_transport; +// double *t_intersect=malloc(n_facets*sizeof(double)); +// int *facet_index = malloc(n_facets*sizeof(int)); +// if (!t_intersect || !facet_index) { +// fprintf(stderr,"Failure allocating list in Union function sample_mesh_intersect - Exit!\n"); +// exit(EXIT_FAILURE); +// } +// Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; +// int **facets = geometry->geometry_parameters.p_mesh_storage->facets; +// // Then loop over every facet +// intersect_transport.rotated_coordinates = rotated_coordinates; +// int inside_vote = 0, outside_vote = 0; +// +// for (int j = 0; j <3; j++){ +// counter = 0; +// Coords test_vector = coords_set(rand01(),rand01(),rand01()); +// intersect_transport.rotated_velocity = test_vector; +// for (iter = 0 ; iter < n_facets ; iter++){ +// intersect_transport.v1 = verts[facets[iter][0]]; +// intersect_transport.v2 = verts[facets[iter][1]]; +// intersect_transport.v3 = verts[facets[iter][2]]; +// possible_t = Moeller_Trumbore_intersection(&intersect_transport); +// if (possible_t > 0){ +// counter++; +// } +// } +// if (counter%2==1){ +// inside_vote++; +// } else { +// outside_vote++; +// } +// if (inside_vote == 2 || outside_vote == 2){ +// break; +// } +// } +// if (inside_vote == 2){ +// return 1; +// } else { +// return 0; +// } +// } // Type for holding intersection and normal @@ -3905,6 +4120,8 @@ int sample_mesh_intersect(double *t, int iter =0; int counter=0; struct Moeller_Trumbore intersect_transport; + // TODO: Allocating T_intersect and facet index might get large with large + // meshes double *t_intersect=malloc(n_facets*sizeof(double)); int *facet_index = malloc(n_facets*sizeof(int)); if (!t_intersect || !facet_index) { @@ -3922,7 +4139,7 @@ int sample_mesh_intersect(double *t, intersect_transport.v2 = verts[facets[iter][1]]; intersect_transport.v3 = verts[facets[iter][2]]; possible_t = Moeller_Trumbore_intersection(&intersect_transport); - if (possible_t != 0){ + if (possible_t != -1){ t_intersect[counter] = possible_t; facet_index[counter] = iter; counter++; @@ -5317,15 +5534,11 @@ int cone_within_cone(struct geometry_struct *geometry_child,struct geometry_stru }; int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct *geometry2) { - // Overlap function should return 1 if the to geometries both cover some volume - // Temporary function - + // Overlap function should return 1 if the to geometries both cover some volume + // TODO: Add fast check to see if point is within largest contained sphere + // Or outside of bounding box. + // Brute force check if there is one point of geometry 1 in 2 and 2 in 1. - - // Should also have a secondary check with edges intersecting on faces. - // Could be made faster with a bounding box (or sphere) - - // Load Variables: struct pointer_to_1d_coords_list shell_points1 = geometry1->shell_points(geometry1,144); struct pointer_to_1d_coords_list shell_points2 = geometry2->shell_points(geometry2,144); @@ -5349,6 +5562,27 @@ int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct free(shell_points1.elements); free(shell_points2.elements); + + // Secondary check with edges intersecting on faces: + if (geometry1->eShape==mesh){ + int **facets = geometry1->geometry_parameters.p_mesh_storage->facets; + Coords *verts = geometry1->geometry_parameters.p_mesh_storage->vertices; + for (i = 0; i < geometry1->geometry_parameters.p_mesh_storage->n_facets; i++){ + if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][1]], geometry2) == 1) return 1; + if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][2]], geometry2) == 1) return 1; + if (existence_of_intersection(verts[facets[i][1]], verts[facets[i][2]], geometry2) == 1) return 1; + } + } + if (geometry2->eShape==mesh){ + int **facets = geometry2->geometry_parameters.p_mesh_storage->facets; + Coords *verts = geometry2->geometry_parameters.p_mesh_storage->vertices; + for (i = 0; i < geometry2->geometry_parameters.p_mesh_storage->n_facets; i++){ + if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][1]], geometry1) == 1) return 1; + if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][2]], geometry1) == 1) return 1; + if (existence_of_intersection(verts[facets[i][1]], verts[facets[i][2]], geometry1) == 1) return 1; + } + } + return 0; }; diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index e552e2188..f45be59bd 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -889,7 +889,7 @@ INITIALIZE printf ("\nERROR IN %s: File %s is read as neither an stl or an off file.", NAME_CURRENT_COMP, filename); exit (1); } - printf ("\nCOMPONENT %s: Number of faces read: %d\t Number of vertices read: %d\n", NAME_CURRENT_COMP, n_verts, n_faces); + printf ("\nCOMPONENT %s: Number of faces read: %d\t Number of vertices read: %d\n", NAME_CURRENT_COMP, n_faces, n_verts); // Loop over all vertices and multiply their positions with coordinate_scale for (int i = 0; i < n_verts; i++) { verts[i] = coords_scale (verts[i], coordinate_scale); From 06ac97b6f99350d968bdad52948b518d61416228 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 6 May 2026 13:58:45 +0200 Subject: [PATCH 03/12] Fix erroneous transposition of rotation matrix for shell points --- mcstas-comps/union/Union_mesh.comp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index f45be59bd..755890e16 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -639,7 +639,7 @@ SHARE for (i = 0; i < number_of_points_in_array; i++) { // Transpose and rotate the points such that they are in the right coordinate system mesh_shell_array.elements[i] = coords_add (mesh_shell_array.elements[i], geometry->center); - mesh_shell_array.elements[i] = rot_apply (geometry->transpose_rotation_matrix, mesh_shell_array.elements[i]); + mesh_shell_array.elements[i] = rot_apply (geometry->rotation_matrix, mesh_shell_array.elements[i]); } return mesh_shell_array; } From d30e602885f762c5e33bc40e2391cd39be765931 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 6 May 2026 14:15:14 +0200 Subject: [PATCH 04/12] Switch to clearer r_within_mesh implementation by using refactoring of Moeller Trumbore intersections --- mcstas-comps/share/union-lib.c | 329 ++++++--------------------------- 1 file changed, 55 insertions(+), 274 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index faa961fad..1230fb219 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -3752,288 +3752,72 @@ double Moeller_Trumbore_intersection(struct Moeller_Trumbore* intersect) { int r_within_mesh(Coords pos,struct geometry_struct *geometry) { - //TODO: Make a single intersection algorithm that all the three mesh intersections - // Can use -// Unpack parameters - - Coords center = geometry->center; - int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; - double *normal_x = geometry->geometry_parameters.p_mesh_storage->normal_x; - double *normal_y = geometry->geometry_parameters.p_mesh_storage->normal_y; - double *normal_z = geometry->geometry_parameters.p_mesh_storage->normal_z; - double *v1_x = geometry->geometry_parameters.p_mesh_storage->v1_x; - double *v1_y = geometry->geometry_parameters.p_mesh_storage->v1_y; - double *v1_z = geometry->geometry_parameters.p_mesh_storage->v1_z; - double *v2_x = geometry->geometry_parameters.p_mesh_storage->v2_x; - double *v2_y = geometry->geometry_parameters.p_mesh_storage->v2_y; - double *v2_z = geometry->geometry_parameters.p_mesh_storage->v2_z; - double *v3_x = geometry->geometry_parameters.p_mesh_storage->v3_x; - double *v3_y = geometry->geometry_parameters.p_mesh_storage->v3_y; - double *v3_z = geometry->geometry_parameters.p_mesh_storage->v3_z; - - double x_new,y_new,z_new; - - // Coordinate transformation - x_new = pos.x - geometry->center.x; - y_new = pos.y - geometry->center.y; - z_new = pos.z - geometry->center.z; - - Coords coordinates = coords_set(x_new,y_new,z_new); - Coords rotated_coordinates; - - rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); - - - int verbal = 0; - // Generate unit direction vector along center axis of meshs - - // Start with vector that points along the mesh in the simple frame, and rotate to global - Coords simple_vector = coords_set(0,1,0); - Coords test_vector = coords_set(0,1,0); - // Check intersections with every single facet: - int iter =0; - int counter=0; int neg_counter=0; - Coords edge1,edge2,h,s,q,tmp,intersect_pos; - double UNION_EPSILON = 1e-27; - double this_facet_t; - double a,f,u,V; - //////printf("\n RWITHIN TEST 1ste"); - for (iter = 0 ; iter < n_facets ; iter++){ - // Intersection with face plane (Möller–Trumbore) - edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); - edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); - - vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); - - a = Dot(edge1,h); - if (a > -UNION_EPSILON && a < UNION_EPSILON){ - //////printf("\n UNION_EPSILON fail"); - } else{ - f = 1.0/a; - s = coords_sub(rotated_coordinates, coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); - u = f * (Dot(s,h)); - if (u < 0.0 || u > 1.0){ - }else{ - //q = vec_prod(s,edge1); - vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); - V = f * Dot(test_vector,q); - if (V < 0.0 || u + V > 1.0){ - } else { - // At this stage we can compute t to find out where the intersection point is on the line. - if (f* Dot(q,edge2) > 0){ - counter++; - - } else { - neg_counter++; - - } - if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ - } else { //printf("\n [%f %f %f] Failed due to being close to surface, E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); - } - - - - } - } - } - } - int C1 = counter; - - int maxC; int sameNr =0; - if (counter % 2 == neg_counter % 2){ - maxC = counter; - sameNr = 1; - } else { - //printf("\n not the same intersection numbers (%i , %i)",counter,neg_counter); - maxC = counter; - sameNr = 0; - } - - if (sameNr == 0){ - test_vector = coords_set(0,0,1); - iter =0; - counter=0; - for (iter = 0 ; iter < n_facets ; iter++){ - // Intersection with face plane (Möller–Trumbore) - edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); - edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); - - vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); - - a = Dot(edge1,h); - if (a > -UNION_EPSILON && a < UNION_EPSILON){ - //////printf("\n UNION_EPSILON fail"); - } else{ - f = 1.0/a; - s = coords_sub(rotated_coordinates , coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); - u = f * (Dot(s,h)); - if (u < 0.0 || u > 1.0){ - }else{ - //q = vec_prod(s,edge1); - vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); - V = f * Dot(test_vector,q); - if (V < 0.0 || u + V > 1.0){ - } else { - // At this stage we can compute t to find out where the intersection point is on the line. - - if (f* Dot(q,edge2) > 0){ - counter++; - - } else { - neg_counter++; - - } - if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ - } else { printf("\n [%f %f %f] Failed due to being close to surface (2. iteration), E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); - } + // r_within_mesh uses a ray casting technique to determine whether or not + // a position is within the mesh. + // It chooses a random velocity, and if the number of intersections + // is uneven, the position is inside the mesh. + // Since this can be numerically unstable, it performs this ray casting 3 times + // and then allows the majority of results to decide whether inside or outside + Coords center = geometry->center; + double x_new,y_new,z_new; + // Coordinate transformation + x_new = pos.x - geometry->center.x; + y_new = pos.y - geometry->center.y; + z_new = pos.z - geometry->center.z; + Coords coordinates = coords_set(x_new,y_new,z_new); + Coords rotated_coordinates; - } - } - } - } - } + rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); - if (counter % 2 == neg_counter % 2){ - maxC = counter; - sameNr = 1; - } else { - printf("\n not the same intersection numbers (%i , %i) second iteration",counter,neg_counter); - maxC = counter; - sameNr = 0; - } + // Check intersections with every single facet: + // First do allocations: + int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; + double possible_t; + int iter =0; + int counter; + struct Moeller_Trumbore intersect_transport; + double *t_intersect=malloc(n_facets*sizeof(double)); + int *facet_index = malloc(n_facets*sizeof(int)); + if (!t_intersect || !facet_index) { + fprintf(stderr,"Failure allocating list in Union function sample_mesh_intersect - Exit!\n"); + exit(EXIT_FAILURE); + } + Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; + int **facets = geometry->geometry_parameters.p_mesh_storage->facets; + // Then loop over every facet + intersect_transport.rotated_coordinates = rotated_coordinates; + int inside_vote = 0, outside_vote = 0; - if (sameNr == 0){ - test_vector = coords_set(1,0,0); - iter =0; - counter=0; + for (int j = 0; j <3; j++){ + counter = 0; + Coords test_vector = coords_set(rand01(),rand01(),rand01()); + intersect_transport.rotated_velocity = test_vector; for (iter = 0 ; iter < n_facets ; iter++){ - // Intersection with face plane (Möller–Trumbore) - edge1 = coords_set(*(v2_x+iter)-*(v1_x+iter),*(v2_y+iter)-*(v1_y+iter),*(v2_z+iter)-*(v1_z+iter)); - edge2 = coords_set(*(v3_x+iter)-*(v1_x+iter),*(v3_y+iter)-*(v1_y+iter),*(v3_z+iter)-*(v1_z+iter)); - - vec_prod(h.x,h.y,h.z,test_vector.x,test_vector.y,test_vector.z,edge2.x,edge2.y,edge2.z); - - a = Dot(edge1,h); - if (a > -UNION_EPSILON && a < UNION_EPSILON){ - //////printf("\n UNION_EPSILON fail"); - } else{ - f = 1.0/a; - s = coords_sub(rotated_coordinates , coords_set(*(v1_x+iter),*(v1_y+iter),*(v1_z+iter))); - u = f * (Dot(s,h)); - if (u < 0.0 || u > 1.0){ - //////printf("\n Nope 1"); - }else{ - //q = vec_prod(s,edge1); - vec_prod(q.x,q.y,q.z,s.x,s.y,s.z,edge1.x,edge1.y,edge1.z); - V = f * Dot(test_vector,q); - if (V < 0.0 || u + V > 1.0){ - } else { - // At this stage we can compute t to find out where the intersection point is on the line. - - if (f* Dot(q,edge2) > 0){ - counter++; - - } else { - neg_counter++; - - } - if (fabs(f* Dot(q,edge2)) > UNION_EPSILON){ - } else { printf("\n [%f %f %f] Failed due to being close to surface (3. iteration), E = %f",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,f* Dot(q,edge2)); - } - - - - } - } - } - } - + intersect_transport.v1 = verts[facets[iter][0]]; + intersect_transport.v2 = verts[facets[iter][1]]; + intersect_transport.v3 = verts[facets[iter][2]]; + possible_t = Moeller_Trumbore_intersection(&intersect_transport); + if (possible_t > 0){ + counter++; + } } - - - if (counter % 2 == neg_counter % 2){ - maxC = counter; + if (counter%2==1){ + inside_vote++; } else { - return 0; + outside_vote++; } - if ( maxC % 2 == 0) { - return 0; - }else{ - return 1; - + if (inside_vote == 2 || outside_vote == 2){ + break; } - + } + if (inside_vote == 2){ + return 1; + } else { return 0; - }; - // r_within_mesh uses a ray casting technique to determine whether or not - // a position is within the mesh. - // It chooses a random velocity, and if the number of intersections - // is uneven, the position is inside the mesh. - // Since this can be numerically unstable, it performs this ray casting 3 times - // and then allows the majority of results to decide whether inside or outside -// -// Coords center = geometry->center; -// double x_new,y_new,z_new; -// -// // Coordinate transformation -// x_new = pos.x - geometry->center.x; -// y_new = pos.y - geometry->center.y; -// z_new = pos.z - geometry->center.z; -// -// Coords coordinates = coords_set(x_new,y_new,z_new); -// Coords rotated_coordinates; -// -// rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); -// -// // Check intersections with every single facet: -// // First do allocations: -// int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; -// double possible_t; -// int iter =0; -// int counter; -// struct Moeller_Trumbore intersect_transport; -// double *t_intersect=malloc(n_facets*sizeof(double)); -// int *facet_index = malloc(n_facets*sizeof(int)); -// if (!t_intersect || !facet_index) { -// fprintf(stderr,"Failure allocating list in Union function sample_mesh_intersect - Exit!\n"); -// exit(EXIT_FAILURE); -// } -// Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; -// int **facets = geometry->geometry_parameters.p_mesh_storage->facets; -// // Then loop over every facet -// intersect_transport.rotated_coordinates = rotated_coordinates; -// int inside_vote = 0, outside_vote = 0; -// -// for (int j = 0; j <3; j++){ -// counter = 0; -// Coords test_vector = coords_set(rand01(),rand01(),rand01()); -// intersect_transport.rotated_velocity = test_vector; -// for (iter = 0 ; iter < n_facets ; iter++){ -// intersect_transport.v1 = verts[facets[iter][0]]; -// intersect_transport.v2 = verts[facets[iter][1]]; -// intersect_transport.v3 = verts[facets[iter][2]]; -// possible_t = Moeller_Trumbore_intersection(&intersect_transport); -// if (possible_t > 0){ -// counter++; -// } -// } -// if (counter%2==1){ -// inside_vote++; -// } else { -// outside_vote++; -// } -// if (inside_vote == 2 || outside_vote == 2){ -// break; -// } -// } -// if (inside_vote == 2){ -// return 1; -// } else { -// return 0; -// } -// } + } +} // Type for holding intersection and normal @@ -4053,9 +3837,6 @@ int compare_intersections(const void *a, const void *b) { } - - - int sample_mesh_intersect(double *t, double *nx, double *ny, double*nz, int *surface_index, From f97122895acb8968c591ef86b4f2add7a2ddcb5f Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 11:33:24 +0200 Subject: [PATCH 05/12] Update rotation to be correct for shell point in mesh. Also add edge intersection checks in mesh_overlaps_mesh --- mcstas-comps/share/union-lib.c | 87 ++++++++++++++++++------------ mcstas-comps/union/Union_mesh.comp | 85 ++++------------------------- 2 files changed, 65 insertions(+), 107 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 1230fb219..a3c013f2c 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -2535,19 +2535,19 @@ struct lines_to_draw draw_circle_with_highest_priority(Coords center,Coords vect * geometry is within the same file. * * To add a new geometry one needs to: - * Write a geometry_storage_struct that contains the paramters needed to describe the geometry + * Write a geometry_storage_struct that contains the parameters needed to describe the geometry * Add a pointer to this storage type in the geometry_parameter_union * Write a function for intersection with line, using the same input scheme as for the others * Write a function checking if a point is within the geometry * Write a function checking if one instance of the geometry overlaps with another * Write a function checking if one instance of the geometry is inside another - * For each exsisting geometry: - * Write a function checking if an instance of this geometry overlaps with an instance of the exsisting - * Write a function checking if an instance of this geometry is inside an instance of the exsisting + * For each existing geometry: + * Write a function checking if an instance of this geometry overlaps with an instance of the existing + * Write a function checking if an instance of this geometry is inside an instance of the existing * Write a function checking if an instance of an existing geometry is inside an instance of this geometry * * Add these functions to geometry to the logic at the end of this file - * Write a component file similar to the exsisting ones, taking the input from the instrument file, and sending + * Write a component file similar to the existing ones, taking the input from the instrument file, and sending * it on to the master component. */ @@ -2583,7 +2583,7 @@ Coords direction_vector; struct mesh_storage{ int n_facets; -int counter; +int n_verts; double *v1_x; double *v1_y; double *v1_z; @@ -3758,25 +3758,17 @@ int r_within_mesh(Coords pos,struct geometry_struct *geometry) { // is uneven, the position is inside the mesh. // Since this can be numerically unstable, it performs this ray casting 3 times // and then allows the majority of results to decide whether inside or outside - Coords center = geometry->center; - double x_new,y_new,z_new; - // Coordinate transformation - x_new = pos.x - geometry->center.x; - y_new = pos.y - geometry->center.y; - z_new = pos.z - geometry->center.z; - - Coords coordinates = coords_set(x_new,y_new,z_new); + Coords coordinates = coords_sub(pos, geometry->center); Coords rotated_coordinates; rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates); - // Check intersections with every single facet: // First do allocations: int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; double possible_t; int iter =0; - int counter; + int n_intersections; struct Moeller_Trumbore intersect_transport; double *t_intersect=malloc(n_facets*sizeof(double)); int *facet_index = malloc(n_facets*sizeof(int)); @@ -3791,7 +3783,7 @@ int r_within_mesh(Coords pos,struct geometry_struct *geometry) { int inside_vote = 0, outside_vote = 0; for (int j = 0; j <3; j++){ - counter = 0; + n_intersections = 0; Coords test_vector = coords_set(rand01(),rand01(),rand01()); intersect_transport.rotated_velocity = test_vector; for (iter = 0 ; iter < n_facets ; iter++){ @@ -3800,10 +3792,10 @@ int r_within_mesh(Coords pos,struct geometry_struct *geometry) { intersect_transport.v3 = verts[facets[iter][2]]; possible_t = Moeller_Trumbore_intersection(&intersect_transport); if (possible_t > 0){ - counter++; + n_intersections++; } } - if (counter%2==1){ + if (n_intersections%2==1){ inside_vote++; } else { outside_vote++; @@ -5320,7 +5312,7 @@ int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct // Or outside of bounding box. // Brute force check if there is one point of geometry 1 in 2 and 2 in 1. - // Load Variables: + // Note! shell points for meshes dont generate a varying number of points: struct pointer_to_1d_coords_list shell_points1 = geometry1->shell_points(geometry1,144); struct pointer_to_1d_coords_list shell_points2 = geometry2->shell_points(geometry2,144); @@ -5339,31 +5331,60 @@ int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct return 1; } } - - free(shell_points1.elements); - free(shell_points2.elements); - - - // Secondary check with edges intersecting on faces: + // Secondary check with edges intersecting on faces if (geometry1->eShape==mesh){ int **facets = geometry1->geometry_parameters.p_mesh_storage->facets; Coords *verts = geometry1->geometry_parameters.p_mesh_storage->vertices; + Coords vert1, vert2, vert3; for (i = 0; i < geometry1->geometry_parameters.p_mesh_storage->n_facets; i++){ - if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][1]], geometry2) == 1) return 1; - if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][2]], geometry2) == 1) return 1; - if (existence_of_intersection(verts[facets[i][1]], verts[facets[i][2]], geometry2) == 1) return 1; + vert1 = shell_points2.elements[facets[i][0]]; + vert2 = shell_points2.elements[facets[i][1]]; + vert3 = shell_points2.elements[facets[i][2]]; + if (existence_of_intersection(vert1, vert2, geometry2) == 1) { + free(shell_points1.elements); + free(shell_points2.elements); + return 1; + } + if (existence_of_intersection(vert1, vert3, geometry2) == 1) { + free(shell_points1.elements); + free(shell_points2.elements); + return 1; + } + if (existence_of_intersection(vert2, vert3, geometry2) == 1) { + free(shell_points1.elements); + free(shell_points2.elements); + + return 1; + } } } if (geometry2->eShape==mesh){ int **facets = geometry2->geometry_parameters.p_mesh_storage->facets; Coords *verts = geometry2->geometry_parameters.p_mesh_storage->vertices; + Coords vert1, vert2, vert3; for (i = 0; i < geometry2->geometry_parameters.p_mesh_storage->n_facets; i++){ - if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][1]], geometry1) == 1) return 1; - if (existence_of_intersection(verts[facets[i][0]], verts[facets[i][2]], geometry1) == 1) return 1; - if (existence_of_intersection(verts[facets[i][1]], verts[facets[i][2]], geometry1) == 1) return 1; + vert1 = shell_points1.elements[facets[i][0]]; + vert2 = shell_points1.elements[facets[i][1]]; + vert3 = shell_points1.elements[facets[i][2]]; + if (existence_of_intersection(vert1, vert2, geometry1) == 1) { + free(shell_points1.elements); + free(shell_points2.elements); + return 1; + } + if (existence_of_intersection(vert1, vert3, geometry1) == 1) { + free(shell_points1.elements); + free(shell_points2.elements); + return 1; + } + if (existence_of_intersection(vert2, vert3, geometry1) == 1) { + free(shell_points1.elements); + free(shell_points2.elements); + return 1; + } } } - + free(shell_points1.elements); + free(shell_points2.elements); return 0; }; diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index 755890e16..1fb540d7f 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -561,86 +561,23 @@ SHARE struct pointer_to_1d_coords_list allocate_shell_points (struct geometry_struct* geometry, int max_number_of_points) { - // Function that returns a number (less than max) of points on the geometry surface - // Run trhough all points in list of faces, and remove dublicates - // There are three points in a face and very often these will be dublicated a few times. This removes dublicates to boost performance down stream... + // Function that returns all vertex points transported to the center of the geometry struct pointer_to_1d_coords_list mesh_shell_array; - int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets; - double* v1_x = geometry->geometry_parameters.p_mesh_storage->v1_x; - double* v1_y = geometry->geometry_parameters.p_mesh_storage->v1_y; - double* v1_z = geometry->geometry_parameters.p_mesh_storage->v1_z; - double* v2_x = geometry->geometry_parameters.p_mesh_storage->v2_x; - double* v2_y = geometry->geometry_parameters.p_mesh_storage->v2_y; - double* v2_z = geometry->geometry_parameters.p_mesh_storage->v2_z; - double* v3_x = geometry->geometry_parameters.p_mesh_storage->v3_x; - double* v3_y = geometry->geometry_parameters.p_mesh_storage->v3_y; - double* v3_z = geometry->geometry_parameters.p_mesh_storage->v3_z; - int number_of_points_in_array = 0; - mesh_shell_array.elements = malloc (3 * n_facets * sizeof (Coords)); - int is_dublicate = 0; - Coords this_vert; - int i, j; - for (i = 0; i < n_facets; i++) { - - // v1 - is_dublicate = 0; - this_vert = coords_set (*(v1_x + i), *(v1_y + i), *(v1_z + i)); - - // test if dublicate - for (j = 0; j < number_of_points_in_array; j++) { - if (this_vert.x == mesh_shell_array.elements[j].x && this_vert.y == mesh_shell_array.elements[j].y && this_vert.z == mesh_shell_array.elements[j].z) { - is_dublicate = 1; - j = number_of_points_in_array; - } - } - if (is_dublicate == 0) { - mesh_shell_array.elements[number_of_points_in_array] = this_vert; - number_of_points_in_array += 1; - } - - // v2 - is_dublicate = 0; - this_vert = coords_set (*(v2_x + i), *(v2_y + i), *(v2_z + i)); - - for (j = 0; j < number_of_points_in_array; j++) { - if (this_vert.x == mesh_shell_array.elements[j].x && this_vert.y == mesh_shell_array.elements[j].y && this_vert.z == mesh_shell_array.elements[j].z) { - is_dublicate = 1; - j = number_of_points_in_array; - } - } - if (is_dublicate == 0) { - mesh_shell_array.elements[number_of_points_in_array] = this_vert; - number_of_points_in_array += 1; - } - - // v3 - is_dublicate = 0; - this_vert = coords_set (*(v3_x + i), *(v3_y + i), *(v3_z + i)); + int n_verts = geometry->geometry_parameters.p_mesh_storage->n_verts; + Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; + mesh_shell_array.elements = malloc (n_verts * sizeof (Coords)); + for (int i = 0; i < n_verts; i++) { + mesh_shell_array.elements[i] = coords_set(verts[i].x, verts[i].y, verts[i].z); + mesh_shell_array.elements[i] = rot_apply(geometry->rotation_matrix, mesh_shell_array.elements[i]); - // test if dublicate - for (j = 0; j < number_of_points_in_array; j++) { - if (this_vert.x == mesh_shell_array.elements[j].x && this_vert.y == mesh_shell_array.elements[j].y && this_vert.z == mesh_shell_array.elements[j].z) { - is_dublicate = 1; - j = number_of_points_in_array; - } - } - if (is_dublicate == 0) { - mesh_shell_array.elements[number_of_points_in_array] = this_vert; - number_of_points_in_array += 1; - } + mesh_shell_array.elements[i] = coords_add (mesh_shell_array.elements[i], geometry->center); + // Transpose and rotate the points such that they are in the right coordinate system } - j = number_of_points_in_array - 1; // Last legal index, currently j is out of bounds. + mesh_shell_array.num_elements = n_verts; - mesh_shell_array.num_elements = number_of_points_in_array; - - for (i = 0; i < number_of_points_in_array; i++) { - // Transpose and rotate the points such that they are in the right coordinate system - mesh_shell_array.elements[i] = coords_add (mesh_shell_array.elements[i], geometry->center); - mesh_shell_array.elements[i] = rot_apply (geometry->rotation_matrix, mesh_shell_array.elements[i]); - } return mesh_shell_array; } @@ -1029,7 +966,7 @@ INITIALIZE vec_prod (mesh_data.normal_x[i], mesh_data.normal_y[i], mesh_data.normal_z[i], edge1.x, edge1.y, edge1.z, edge2.x, edge2.y, edge2.z); } mesh_data.n_facets = n_faces; - mesh_data.counter = n_verts; + mesh_data.n_verts = n_verts; mesh_vol.geometry.number_of_faces = 1; mesh_vol.geometry.surface_stack_for_each_face = malloc (mesh_vol.geometry.number_of_faces * sizeof (struct surface_stack_struct)); From 6a8344918c0b63cf99ba24bfb79d7a2d0b8f8204 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 12:04:26 +0200 Subject: [PATCH 06/12] Fix mismatch in mesh_overlap_mesh when checking for edge intersection. Also add test instrument to check for the edge intersection --- .../Test_mesh_boxes/Test_mesh_boxes.instr | 137 ++++++++++++++++++ .../Tests_union/Test_mesh_boxes/box.stl | 87 +++++++++++ .../Test_mesh_boxes/single_blade.stl | 87 +++++++++++ mcstas-comps/share/union-lib.c | 18 ++- 4 files changed, 323 insertions(+), 6 deletions(-) create mode 100644 mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr create mode 100644 mcstas-comps/examples/Tests_union/Test_mesh_boxes/box.stl create mode 100644 mcstas-comps/examples/Tests_union/Test_mesh_boxes/single_blade.stl diff --git a/mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr new file mode 100644 index 000000000..7d747ebda --- /dev/null +++ b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr @@ -0,0 +1,137 @@ +/******************************************************************************** +* +* McStas, neutron ray-tracing package +* Copyright (C) 1997-2008, All rights reserved +* Risoe National Laboratory, Roskilde, Denmark +* Institut Laue Langevin, Grenoble, France +* +* This file was written by McStasScript, which is a +* python based McStas instrument generator written by +* Mads Bertelsen in 2019 while employed at the +* European Spallation Source Data Management and +* Software Centre +* +* Instrument: Test_mesh_boxes +* +* %Identification +* Written by: Sam Lambrick with Python McStas Instrument Generator +* Date: 10:07:56 on March 05, 2026 +* Origin: ESS DMSC / ISIS +* %INSTRUMENT_SITE: Tests_union +* +* +* A small instrument to investigate mesh boxes in a user configuration +* +* %Description +* Instrument provided by Sam Lambrick from ISIS, that uses three meshes to +* make their boxes. +* +* %Example: Unused_par=1 Detector: PSD_detector_I=10.7752 +* +* %Parameters +* +* %End +********************************************************************************/ + +DEFINE INSTRUMENT Test_mesh_boxes( + Unused_par=0 +) + +DECLARE +%{ +%} + +INITIALIZE +%{ +// Start of initialize for generated test_union_filter +%} + +TRACE +COMPONENT armOrigin = Arm() +AT (0, 0, 0) ABSOLUTE + +COMPONENT Source = Moderator( + radius = 0.1, Emin = 1, + Emax = 20, dist = 2, + focus_xw = 0.1, focus_yh = 0.1, + Ec = 1000) +AT (0, 0, 0) RELATIVE armOrigin + +COMPONENT armFilter = Arm() +AT (0, 0, 2) RELATIVE armOrigin + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT Be_NCrystal = NCrystal_process( + cfg = "Be_sg194.ncmat", interact_fraction = 1) +AT (0, 0, 0) ABSOLUTE + +COMPONENT Cd_absorber = Incoherent_process( + sigma = 3.46, unit_cell_volume = 43.2, + interact_fraction = 1) +AT (0, 0, 0) ABSOLUTE + +COMPONENT Be = Union_make_material( + process_string = "Be_NCrystal", my_absorption = 0.0) +AT (0, 0, 0) ABSOLUTE + +COMPONENT Cd = Union_make_material( + process_string = "Cd_absorber", my_absorption = 11627.906976744185) +AT (0, 0, 0) ABSOLUTE + +COMPONENT Be_envelope = Union_mesh( + filename = "./box.stl", material_string = "Be", + priority = 1, coordinate_scale = 1) +AT (0, 0, 0) RELATIVE armFilter + +COMPONENT Cd_blade1 = Union_mesh( + filename = "./single_blade.stl", material_string = "Cd", + priority = 2, coordinate_scale = 1) +AT (0, -0.02, 0) RELATIVE armFilter + +COMPONENT Cd_blade2 = Union_mesh( + filename = "./single_blade.stl", material_string = "Cd", + priority = 2.5, coordinate_scale = 1) +AT (0, 0.02, 0) RELATIVE armFilter + +COMPONENT zy_logger = Union_logger_2D_space( + target_geometry = "Be_envelope,Cd_blade1,Cd_blade2", D_direction_1 = "z", + D1_min = -0.02, D1_max = 0.02, + n1 = 100, D_direction_2 = "y", + D2_min = -0.05, D2_max = 0.05, + n2 = 100, filename = "zy_logger.dat") +AT (0, 0, 0) RELATIVE armFilter + +COMPONENT zy_logger_abs = Union_abs_logger_2D_space( + target_geometry = "Be_envelope,Cd_blade1,Cd_blade2", D_direction_1 = "z", + D1_min = -0.02, D1_max = 0.02, + n1 = 100, D_direction_2 = "y", + D2_min = -0.05, D2_max = 0.05, + n2 = 100, filename = "zy_logge_abs.dat") +AT (0, 0, 0) RELATIVE armFilter + +COMPONENT master = Union_master(verbal=1) +AT (0, 0, 0) RELATIVE armFilter + +COMPONENT stop = Union_stop() +AT (0, 0, 0) RELATIVE armFilter + +COMPONENT psd_detector = PSD_monitor( + nx = 201, ny = 201, + filename = "psd_detector.dat", xwidth = 0.1, + yheight = 0.1, restore_neutron = 1) +AT (0, 0, 0.15) RELATIVE armFilter + +COMPONENT e_detector = E_monitor( + filename = "e_detector.dat", xwidth = 0.1, + yheight = 0.1, Emin = 1, + Emax = 20) +AT (0, 0, 0.15) RELATIVE armFilter + +FINALLY +%{ +// Start of finally for generated test_union_filter +%} + +END diff --git a/mcstas-comps/examples/Tests_union/Test_mesh_boxes/box.stl b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/box.stl new file mode 100644 index 000000000..8d3d07110 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/box.stl @@ -0,0 +1,87 @@ +solid +facet normal -1.0 0.0 0.0 +outer loop +vertex -0.025 -0.025 0.005 +vertex -0.025 0.025 0.005 +vertex -0.025 -0.025 -0.005 +endloop +endfacet +facet normal 0.0 -1.0 0.0 +outer loop +vertex 0.025 -0.025 -0.005 +vertex -0.025 -0.025 0.005 +vertex -0.025 -0.025 -0.005 +endloop +endfacet +facet normal -1.0 0.0 0.0 +outer loop +vertex -0.025 -0.025 -0.005 +vertex -0.025 0.025 0.005 +vertex -0.025 0.025 -0.005 +endloop +endfacet +facet normal 0.0 0.0 -1.0 +outer loop +vertex -0.025 0.025 -0.005 +vertex 0.025 -0.025 -0.005 +vertex -0.025 -0.025 -0.005 +endloop +endfacet +facet normal 0.0 0.0 1.0 +outer loop +vertex -0.025 -0.025 0.005 +vertex 0.025 0.025 0.005 +vertex -0.025 0.025 0.005 +endloop +endfacet +facet normal 0.0 -1.0 0.0 +outer loop +vertex 0.025 -0.025 0.005 +vertex -0.025 -0.025 0.005 +vertex 0.025 -0.025 -0.005 +endloop +endfacet +facet normal 0.0 0.0 1.0 +outer loop +vertex 0.025 -0.025 0.005 +vertex 0.025 0.025 0.005 +vertex -0.025 -0.025 0.005 +endloop +endfacet +facet normal 0.0 1.0 0.0 +outer loop +vertex -0.025 0.025 0.005 +vertex 0.025 0.025 0.005 +vertex -0.025 0.025 -0.005 +endloop +endfacet +facet normal 0.0 0.0 -1.0 +outer loop +vertex 0.025 0.025 -0.005 +vertex 0.025 -0.025 -0.005 +vertex -0.025 0.025 -0.005 +endloop +endfacet +facet normal 0.0 1.0 0.0 +outer loop +vertex -0.025 0.025 -0.005 +vertex 0.025 0.025 0.005 +vertex 0.025 0.025 -0.005 +endloop +endfacet +facet normal 1.0 0.0 0.0 +outer loop +vertex 0.025 0.025 -0.005 +vertex 0.025 -0.025 0.005 +vertex 0.025 -0.025 -0.005 +endloop +endfacet +facet normal 1.0 0.0 0.0 +outer loop +vertex 0.025 0.025 0.005 +vertex 0.025 -0.025 0.005 +vertex 0.025 0.025 -0.005 +endloop +endfacet + +endsolid diff --git a/mcstas-comps/examples/Tests_union/Test_mesh_boxes/single_blade.stl b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/single_blade.stl new file mode 100644 index 000000000..4c484ac4d --- /dev/null +++ b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/single_blade.stl @@ -0,0 +1,87 @@ +solid +facet normal -1.0 0.0 0.0 +outer loop +vertex -0.03 -0.00202669681901379 0.005735198340406318 +vertex -0.03 -5.708131298937399e-05 0.006082494695740179 +vertex -0.03 5.708131298937399e-05 -0.006082494695740179 +endloop +endfacet +facet normal 0.0 -0.9848077530122082 -0.17364817766693036 +outer loop +vertex 0.03 5.708131298937399e-05 -0.006082494695740179 +vertex -0.03 -0.00202669681901379 0.005735198340406318 +vertex -0.03 5.708131298937399e-05 -0.006082494695740179 +endloop +endfacet +facet normal -1.0 0.0 0.0 +outer loop +vertex -0.03 5.708131298937399e-05 -0.006082494695740179 +vertex -0.03 -5.708131298937399e-05 0.006082494695740179 +vertex -0.03 0.00202669681901379 -0.005735198340406318 +endloop +endfacet +facet normal 0.0 0.17364817766693036 -0.9848077530122082 +outer loop +vertex -0.03 0.00202669681901379 -0.005735198340406318 +vertex 0.03 5.708131298937399e-05 -0.006082494695740179 +vertex -0.03 5.708131298937399e-05 -0.006082494695740179 +endloop +endfacet +facet normal 0.0 -0.17364817766693036 0.9848077530122082 +outer loop +vertex -0.03 -0.00202669681901379 0.005735198340406318 +vertex 0.03 -5.708131298937399e-05 0.006082494695740179 +vertex -0.03 -5.708131298937399e-05 0.006082494695740179 +endloop +endfacet +facet normal 0.0 -0.9848077530122082 -0.17364817766693036 +outer loop +vertex 0.03 -0.00202669681901379 0.005735198340406318 +vertex -0.03 -0.00202669681901379 0.005735198340406318 +vertex 0.03 5.708131298937399e-05 -0.006082494695740179 +endloop +endfacet +facet normal 0.0 -0.17364817766693036 0.9848077530122082 +outer loop +vertex 0.03 -0.00202669681901379 0.005735198340406318 +vertex 0.03 -5.708131298937399e-05 0.006082494695740179 +vertex -0.03 -0.00202669681901379 0.005735198340406318 +endloop +endfacet +facet normal 0.0 0.9848077530122082 0.17364817766693036 +outer loop +vertex -0.03 -5.708131298937399e-05 0.006082494695740179 +vertex 0.03 -5.708131298937399e-05 0.006082494695740179 +vertex -0.03 0.00202669681901379 -0.005735198340406318 +endloop +endfacet +facet normal 0.0 0.17364817766693036 -0.9848077530122082 +outer loop +vertex 0.03 0.00202669681901379 -0.005735198340406318 +vertex 0.03 5.708131298937399e-05 -0.006082494695740179 +vertex -0.03 0.00202669681901379 -0.005735198340406318 +endloop +endfacet +facet normal 0.0 0.9848077530122082 0.17364817766693036 +outer loop +vertex -0.03 0.00202669681901379 -0.005735198340406318 +vertex 0.03 -5.708131298937399e-05 0.006082494695740179 +vertex 0.03 0.00202669681901379 -0.005735198340406318 +endloop +endfacet +facet normal 1.0 0.0 0.0 +outer loop +vertex 0.03 0.00202669681901379 -0.005735198340406318 +vertex 0.03 -0.00202669681901379 0.005735198340406318 +vertex 0.03 5.708131298937399e-05 -0.006082494695740179 +endloop +endfacet +facet normal 1.0 0.0 0.0 +outer loop +vertex 0.03 -5.708131298937399e-05 0.006082494695740179 +vertex 0.03 -0.00202669681901379 0.005735198340406318 +vertex 0.03 0.00202669681901379 -0.005735198340406318 +endloop +endfacet + +endsolid diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index a3c013f2c..723ed48b7 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -5337,35 +5337,39 @@ int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct Coords *verts = geometry1->geometry_parameters.p_mesh_storage->vertices; Coords vert1, vert2, vert3; for (i = 0; i < geometry1->geometry_parameters.p_mesh_storage->n_facets; i++){ - vert1 = shell_points2.elements[facets[i][0]]; - vert2 = shell_points2.elements[facets[i][1]]; - vert3 = shell_points2.elements[facets[i][2]]; + vert1 = shell_points1.elements[facets[i][0]]; + vert2 = shell_points1.elements[facets[i][1]]; + vert3 = shell_points1.elements[facets[i][2]]; if (existence_of_intersection(vert1, vert2, geometry2) == 1) { free(shell_points1.elements); free(shell_points2.elements); + printf("t1\n"); return 1; } if (existence_of_intersection(vert1, vert3, geometry2) == 1) { free(shell_points1.elements); free(shell_points2.elements); + printf("t2\n"); return 1; } if (existence_of_intersection(vert2, vert3, geometry2) == 1) { free(shell_points1.elements); free(shell_points2.elements); + printf("t3\n"); return 1; } } } + printf("t4\n"); if (geometry2->eShape==mesh){ int **facets = geometry2->geometry_parameters.p_mesh_storage->facets; Coords *verts = geometry2->geometry_parameters.p_mesh_storage->vertices; Coords vert1, vert2, vert3; for (i = 0; i < geometry2->geometry_parameters.p_mesh_storage->n_facets; i++){ - vert1 = shell_points1.elements[facets[i][0]]; - vert2 = shell_points1.elements[facets[i][1]]; - vert3 = shell_points1.elements[facets[i][2]]; + vert1 = shell_points2.elements[facets[i][0]]; + vert2 = shell_points2.elements[facets[i][1]]; + vert3 = shell_points2.elements[facets[i][2]]; if (existence_of_intersection(vert1, vert2, geometry1) == 1) { free(shell_points1.elements); free(shell_points2.elements); @@ -7491,7 +7495,9 @@ void generate_overlap_lists(struct pointer_to_1d_int_list **true_overlap_lists, if (cone_overlaps_cone(&Volumes[parent]->geometry,&Volumes[child]->geometry)) temp_list_local.elements[used_elements++] = child; } else if (strcmp("mesh",Volumes[parent]->geometry.shape) == 0 && strcmp("mesh",Volumes[child]->geometry.shape) == 0) { + printf("Just before mesh overlaps mesh\n"); if (mesh_overlaps_mesh(&Volumes[parent]->geometry,&Volumes[child]->geometry)) temp_list_local.elements[used_elements++] = child; + printf("Just after mesh overlaps mesh\n"); } else if (strcmp("box",Volumes[parent]->geometry.shape) == 0 && strcmp("cylinder",Volumes[child]->geometry.shape) == 0) { if (box_overlaps_cylinder(&Volumes[parent]->geometry,&Volumes[child]->geometry)) temp_list_local.elements[used_elements++] = child; From 12e1b15840f22e85c876300dcdd24d24d3215128 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 12:09:23 +0200 Subject: [PATCH 07/12] Fix upper case error in Test definition of Test_mesh_boxes --- .../examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr index 7d747ebda..d2562c58c 100644 --- a/mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr +++ b/mcstas-comps/examples/Tests_union/Test_mesh_boxes/Test_mesh_boxes.instr @@ -26,7 +26,7 @@ * Instrument provided by Sam Lambrick from ISIS, that uses three meshes to * make their boxes. * -* %Example: Unused_par=1 Detector: PSD_detector_I=10.7752 +* %Example: Unused_par=1 Detector: psd_detector_I=10.7752 * * %Parameters * From c32fe7daeca6f13b1b04516883a8be22e7e07709 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 13:12:15 +0200 Subject: [PATCH 08/12] Remove the last vestiges of vertices represented as 9 double arrays --- mcstas-comps/share/union-lib.c | 9 --- mcstas-comps/union/Union_mesh.comp | 92 ++++++++++-------------------- 2 files changed, 29 insertions(+), 72 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 723ed48b7..4658c339f 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -2584,15 +2584,6 @@ Coords direction_vector; struct mesh_storage{ int n_facets; int n_verts; -double *v1_x; -double *v1_y; -double *v1_z; -double *v2_x; -double *v2_y; -double *v2_z; -double *v3_x; -double *v3_y; -double *v3_z; double *normal_x; double *normal_y; double *normal_z; diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index 1fb540d7f..98ad53dce 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -443,23 +443,18 @@ SHARE // Function to call in mcdisplay section of the sample component for this volume int n_facets = Geometries[index]->geometry_parameters.p_mesh_storage->n_facets; - double* v1_x = Geometries[index]->geometry_parameters.p_mesh_storage->v1_x; - double* v1_y = Geometries[index]->geometry_parameters.p_mesh_storage->v1_y; - double* v1_z = Geometries[index]->geometry_parameters.p_mesh_storage->v1_z; - double* v2_x = Geometries[index]->geometry_parameters.p_mesh_storage->v2_x; - double* v2_y = Geometries[index]->geometry_parameters.p_mesh_storage->v2_y; - double* v2_z = Geometries[index]->geometry_parameters.p_mesh_storage->v2_z; - double* v3_x = Geometries[index]->geometry_parameters.p_mesh_storage->v3_x; - double* v3_y = Geometries[index]->geometry_parameters.p_mesh_storage->v3_y; - double* v3_z = Geometries[index]->geometry_parameters.p_mesh_storage->v3_z; - + int **facets = Geometries[index]->geometry_parameters.p_mesh_storage->facets; Coords center = Geometries[index]->center; + struct pointer_to_1d_coords_list points_struct = Geometries[index]->shell_points(Geometries[index], 1); + Coords* points = points_struct.elements; + + struct lines_to_draw lines_to_draw_temp; lines_to_draw_temp.number_of_lines = 0; Coords point1, point2, point3; - int iterate, i, j; + int i, j, k; int print1 = 0; int print2 = 0; int print3 = 0; @@ -474,55 +469,55 @@ SHARE int counter = 0; // For every triangle it should add three lines - for (iterate = 0; iterate < n_facets; iterate++) { - point1 = coords_add (rot_apply (Geometries[index]->rotation_matrix, coords_set (*(v1_x + iterate), *(v1_y + iterate), *(v1_z + iterate))), center); - point2 = coords_add (rot_apply (Geometries[index]->rotation_matrix, coords_set (*(v2_x + iterate), *(v2_y + iterate), *(v2_z + iterate))), center); - point3 = coords_add (rot_apply (Geometries[index]->rotation_matrix, coords_set (*(v3_x + iterate), *(v3_y + iterate), *(v3_z + iterate))), center); + for (i = 0; i < n_facets; i++) { + point1 = points[facets[i][0]]; + point2 = points[facets[i][1]]; + point3 = points[facets[i][2]]; print1 = 1; print2 = 1; print3 = 1; // Make sure it does not print a line if it is already printed.... (might take a while?) - for (i = 0; i < counter; i++) { - if (print1 == 1 && coord_comp (point1, list_startpoints[i])) { - for (j = 0; j < counter; j++) { - if (coord_comp (point2, list_startpoints[i])) { + for (j = 0; j < counter; j++) { + if (print1 == 1 && coord_comp (point1, list_startpoints[j])) { + for (k = 0; k < counter; k++) { + if (coord_comp (point2, list_startpoints[j])) { print1 = 0; } } } - if (print2 == 1 && coord_comp (point2, list_startpoints[i])) { - for (j = 0; j < counter; j++) { - if (coord_comp (point1, list_startpoints[i])) { + if (print2 == 1 && coord_comp (point2, list_startpoints[j])) { + for (k = 0; k < counter; k++) { + if (coord_comp (point1, list_startpoints[j])) { print1 = 0; } } } - if (print2 == 1 && coord_comp (point2, list_startpoints[i])) { - for (j = 0; j < counter; j++) { - if (coord_comp (point3, list_startpoints[i])) { + if (print2 == 1 && coord_comp (point2, list_startpoints[j])) { + for (k = 0; k < counter; k++) { + if (coord_comp (point3, list_startpoints[j])) { print2 = 0; } } } - if (print3 == 1 && coord_comp (point3, list_startpoints[i])) { - for (j = 0; j < counter; j++) { - if (coord_comp (point2, list_startpoints[i])) { + if (print3 == 1 && coord_comp (point3, list_startpoints[j])) { + for (k = 0; k < counter; k++) { + if (coord_comp (point2, list_startpoints[j])) { print2 = 0; } } } - if (print1 == 1 && coord_comp (point1, list_startpoints[i])) { - for (j = 0; j < counter; j++) { - if (coord_comp (point1, list_startpoints[i])) { + if (print1 == 1 && coord_comp (point1, list_startpoints[j])) { + for (k = 0; k < counter; k++) { + if (coord_comp (point1, list_startpoints[j])) { print3 = 0; } } } - if (print3 == 1 && coord_comp (point3, list_startpoints[i])) { - for (j = 0; j < counter; j++) { - if (coord_comp (point1, list_startpoints[i])) { + if (print3 == 1 && coord_comp (point3, list_startpoints[j])) { + for (k = 0; k < counter; k++) { + if (coord_comp (point1, list_startpoints[j])) { print3 = 0; } } @@ -932,15 +927,6 @@ INITIALIZE Coords edge1, edge2; // Allocate space vertices in mesh data // TODO: Change data format to be in a less verbose data structure, or a more effective one. - mesh_data.v1_x = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v1_y = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v1_z = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v2_x = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v2_y = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v2_z = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v3_x = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v3_y = (double*)malloc (n_faces * sizeof (double)); - mesh_data.v3_z = (double*)malloc (n_faces * sizeof (double)); mesh_data.normal_x = (double*)malloc (n_faces * sizeof (double)); mesh_data.normal_y = (double*)malloc (n_faces * sizeof (double)); mesh_data.normal_z = (double*)malloc (n_faces * sizeof (double)); @@ -950,17 +936,6 @@ INITIALIZE vert1 = verts[faces[i][0]]; vert2 = verts[faces[i][1]]; vert3 = verts[faces[i][2]]; - mesh_data.v1_x[i] = vert1.x; - mesh_data.v1_y[i] = vert1.y; - mesh_data.v1_z[i] = vert1.z; - - mesh_data.v2_x[i] = vert2.x; - mesh_data.v2_y[i] = vert2.y; - mesh_data.v2_z[i] = vert2.z; - - mesh_data.v3_x[i] = vert3.x; - mesh_data.v3_y[i] = vert3.y; - mesh_data.v3_z[i] = vert3.z; edge1 = coords_sub (vert1, vert2); edge2 = coords_sub (vert3, vert1); vec_prod (mesh_data.normal_x[i], mesh_data.normal_y[i], mesh_data.normal_z[i], edge1.x, edge1.y, edge1.z, edge2.x, edge2.y, edge2.z); @@ -1017,15 +992,6 @@ TRACE FINALLY %{ - free (mesh_data.v1_x); - free (mesh_data.v1_y); - free (mesh_data.v1_z); - free (mesh_data.v2_x); - free (mesh_data.v2_y); - free (mesh_data.v2_z); - free (mesh_data.v3_x); - free (mesh_data.v3_y); - free (mesh_data.v3_z); %} END From f635cc07ca782e67f717aa173751ab7838d2ae68 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 13:19:34 +0200 Subject: [PATCH 09/12] Add error statement if file load fails for Union mesh --- mcstas-comps/union/Union_mesh.comp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index 98ad53dce..c96126dc9 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -130,15 +130,20 @@ SHARE uint16_t attribute_byte_count; } Triangle; #pragma pack(pop) + int check_open_file_errors(FILE* fp, char* filename){ + if (fp == NULL) { + printf("\n\nERROR: %s could not be opened by Union_mesh component.\n\n", filename); + exit(1); + } + } void read_stl (char* filename, int* n_verts, int* n_faces, int* n_edges, Coords** verts, int*** faces, char* comp_name) { unsigned char buffer[1000]; Triangle* triangles; FILE* file = Open_File (filename, "r", NULL); - if (!file) { - perror ("ERROR: Failed to open file"); - exit (1); - } + + check_open_file_errors(file, filename); + // Check if the file is binary or ASCII int file_is_ascii = 0; char word[6]; // "solid" + null terminator @@ -307,7 +312,9 @@ SHARE void read_off_file_vertices_and_faces (char* filename, int* n_verts, int* n_faces, int* n_edges, Coords** verts, int*** faces, char* comp_name) { FILE* fp; - fp = Open_File (filename, "r", NULL); + fp = Open_File (filename, "r", NULL); + check_open_file_errors(fp, filename); + // TODO: Make tests to verify the contents of FILE int n_lines = 0; char buffer[250]; From b744553b54edf638d3496c6a8ebe657f987895f3 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 13:23:16 +0200 Subject: [PATCH 10/12] Use linter on Union_mesh --- mcstas-comps/union/Union_mesh.comp | 40 ++++++++++++++---------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index c96126dc9..32da3dfb9 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -130,10 +130,11 @@ SHARE uint16_t attribute_byte_count; } Triangle; #pragma pack(pop) - int check_open_file_errors(FILE* fp, char* filename){ + int + check_open_file_errors (FILE* fp, char* filename) { if (fp == NULL) { - printf("\n\nERROR: %s could not be opened by Union_mesh component.\n\n", filename); - exit(1); + printf ("\n\nERROR: %s could not be opened by Union_mesh component.\n\n", filename); + exit (1); } } void @@ -142,7 +143,7 @@ SHARE Triangle* triangles; FILE* file = Open_File (filename, "r", NULL); - check_open_file_errors(file, filename); + check_open_file_errors (file, filename); // Check if the file is binary or ASCII int file_is_ascii = 0; @@ -312,8 +313,8 @@ SHARE void read_off_file_vertices_and_faces (char* filename, int* n_verts, int* n_faces, int* n_edges, Coords** verts, int*** faces, char* comp_name) { FILE* fp; - fp = Open_File (filename, "r", NULL); - check_open_file_errors(fp, filename); + fp = Open_File (filename, "r", NULL); + check_open_file_errors (fp, filename); // TODO: Make tests to verify the contents of FILE int n_lines = 0; @@ -450,13 +451,12 @@ SHARE // Function to call in mcdisplay section of the sample component for this volume int n_facets = Geometries[index]->geometry_parameters.p_mesh_storage->n_facets; - int **facets = Geometries[index]->geometry_parameters.p_mesh_storage->facets; + int** facets = Geometries[index]->geometry_parameters.p_mesh_storage->facets; Coords center = Geometries[index]->center; - struct pointer_to_1d_coords_list points_struct = Geometries[index]->shell_points(Geometries[index], 1); + struct pointer_to_1d_coords_list points_struct = Geometries[index]->shell_points (Geometries[index], 1); Coords* points = points_struct.elements; - struct lines_to_draw lines_to_draw_temp; lines_to_draw_temp.number_of_lines = 0; @@ -568,11 +568,11 @@ SHARE struct pointer_to_1d_coords_list mesh_shell_array; int n_verts = geometry->geometry_parameters.p_mesh_storage->n_verts; - Coords *verts = geometry->geometry_parameters.p_mesh_storage->vertices; + Coords* verts = geometry->geometry_parameters.p_mesh_storage->vertices; mesh_shell_array.elements = malloc (n_verts * sizeof (Coords)); for (int i = 0; i < n_verts; i++) { - mesh_shell_array.elements[i] = coords_set(verts[i].x, verts[i].y, verts[i].z); - mesh_shell_array.elements[i] = rot_apply(geometry->rotation_matrix, mesh_shell_array.elements[i]); + mesh_shell_array.elements[i] = coords_set (verts[i].x, verts[i].y, verts[i].z); + mesh_shell_array.elements[i] = rot_apply (geometry->rotation_matrix, mesh_shell_array.elements[i]); mesh_shell_array.elements[i] = coords_add (mesh_shell_array.elements[i], geometry->center); // Transpose and rotate the points such that they are in the right coordinate system @@ -833,7 +833,6 @@ INITIALIZE for (int i = 0; i < n_verts; i++) { verts[i] = coords_scale (verts[i], coordinate_scale); } - // Make lists that will be used to calculate the number of edges int** vert_pairs = (int**)malloc (sizeof (int*) * 3 * n_faces); @@ -911,13 +910,13 @@ INITIALIZE if (test_rad > radius) { radius = test_rad; } - if (i == 0){ - mesh_data.Bounding_Box_Extremes[0] = verts[i].x; - mesh_data.Bounding_Box_Extremes[1] = verts[i].x; - mesh_data.Bounding_Box_Extremes[2] = verts[i].y; - mesh_data.Bounding_Box_Extremes[3] = verts[i].y; - mesh_data.Bounding_Box_Extremes[4] = verts[i].z; - mesh_data.Bounding_Box_Extremes[5] = verts[i].z; + if (i == 0) { + mesh_data.Bounding_Box_Extremes[0] = verts[i].x; + mesh_data.Bounding_Box_Extremes[1] = verts[i].x; + mesh_data.Bounding_Box_Extremes[2] = verts[i].y; + mesh_data.Bounding_Box_Extremes[3] = verts[i].y; + mesh_data.Bounding_Box_Extremes[4] = verts[i].z; + mesh_data.Bounding_Box_Extremes[5] = verts[i].z; } else { mesh_data.Bounding_Box_Extremes[0] = verts[i].x > mesh_data.Bounding_Box_Extremes[0] ? verts[i].x : mesh_data.Bounding_Box_Extremes[0]; mesh_data.Bounding_Box_Extremes[1] = verts[i].x < mesh_data.Bounding_Box_Extremes[1] ? verts[i].x : mesh_data.Bounding_Box_Extremes[1]; @@ -925,7 +924,6 @@ INITIALIZE mesh_data.Bounding_Box_Extremes[3] = verts[i].y < mesh_data.Bounding_Box_Extremes[3] ? verts[i].y : mesh_data.Bounding_Box_Extremes[3]; mesh_data.Bounding_Box_Extremes[4] = verts[i].z > mesh_data.Bounding_Box_Extremes[4] ? verts[i].z : mesh_data.Bounding_Box_Extremes[4]; mesh_data.Bounding_Box_Extremes[5] = verts[i].z < mesh_data.Bounding_Box_Extremes[5] ? verts[i].z : mesh_data.Bounding_Box_Extremes[5]; - } } mesh_data.Bounding_Box_Radius = radius; From 86b1679967f7ed21a569dffc8faca45fcf9ba79e Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 13:25:30 +0200 Subject: [PATCH 11/12] Remove unused function --- mcstas-comps/union/Union_mesh.comp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index 32da3dfb9..e2f630cc8 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -275,6 +275,7 @@ SHARE int vertex_is_unique; int x_are_equal, y_are_equal, z_are_equal; int* map_old_to_unique = malloc (sizeof (int) * *n_verts); + int unique_counter = 0; int vert_index_in_faces = 0; for (int i = 0; i < *n_verts; i++) { @@ -440,11 +441,6 @@ SHARE } return 0; }; - Coords - get_coords_from_string (char* line) { - Coords vert_coords = coords_set (1, 2, 1); - return vert_coords; - } void mcdisplay_mesh_function (struct lines_to_draw* lines_to_draw_output, int index, struct geometry_struct** Geometries, int number_of_volumes) { From ffe4e7809bfd2af9e8cf7c242c8f140e0c02fe80 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 7 May 2026 14:10:09 +0200 Subject: [PATCH 12/12] Fix linter complaints --- mcstas-comps/union/Union_mesh.comp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index e2f630cc8..33f6fa38d 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -130,12 +130,13 @@ SHARE uint16_t attribute_byte_count; } Triangle; #pragma pack(pop) - int + void check_open_file_errors (FILE* fp, char* filename) { if (fp == NULL) { printf ("\n\nERROR: %s could not be opened by Union_mesh component.\n\n", filename); exit (1); } + return; } void read_stl (char* filename, int* n_verts, int* n_faces, int* n_edges, Coords** verts, int*** faces, char* comp_name) { @@ -275,7 +276,10 @@ SHARE int vertex_is_unique; int x_are_equal, y_are_equal, z_are_equal; int* map_old_to_unique = malloc (sizeof (int) * *n_verts); - + if (map_old_to_unique == NULL || unique_vertices == NULL) { + printf("\nERROR ON MALLOC IN UNION MESH STL READ\n"); + exit(1); + } int unique_counter = 0; int vert_index_in_faces = 0; for (int i = 0; i < *n_verts; i++) { @@ -304,6 +308,7 @@ SHARE } *n_verts = unique_counter; *verts = (Coords*)realloc (*verts, sizeof (Coords) * unique_counter); + for (int i = 0; i < unique_counter; i++) { (*verts)[i] = unique_vertices[i]; }