From 24c751f03aa5dbe775d68df139c65806ff6f7cff Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 28 Apr 2026 15:14:43 +0200 Subject: [PATCH] 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]];