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..d2562c58c --- /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 cc683f10e..9b2f8c3ae 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -2498,19 +2498,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. */ @@ -2546,21 +2546,15 @@ Coords direction_vector; struct mesh_storage{ int n_facets; -int counter; -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; +int n_verts; double *normal_x; 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; }; @@ -3632,220 +3626,116 @@ 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) { -// 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++; +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; +}; - } 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. +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 (f* Dot(q,edge2) > 0){ - counter++; + intersect->h = coords_xp(intersect->rotated_velocity, intersect->edge2); + intersect->a = Dot(intersect->edge1, intersect->h); - } 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 -1; + } 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 -1; } 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++; + return 0; +} - } 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)); - } - - - - } - } - } - } +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 + // Coordinate transformation + 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 n_intersections; + 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++){ + n_intersections = 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){ + n_intersections++; + } } - - - if (counter % 2 == neg_counter % 2){ - maxC = counter; + if (n_intersections%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 { @@ -3863,112 +3753,96 @@ 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) { + // 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; - 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); - } + // 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); - - 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; - 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; + // 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) { + 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 != -1){ + t_intersect[counter] = possible_t; + facet_index[counter] = iter; + counter++; } + } *num_solutions = counter; @@ -3980,36 +3854,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; }; @@ -5346,16 +5223,12 @@ 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: + // 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); @@ -5374,10 +5247,64 @@ int mesh_overlaps_mesh(struct geometry_struct *geometry1,struct geometry_struct return 1; } } - + // 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++){ + 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_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); + 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; }; @@ -7450,7 +7377,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; diff --git a/mcstas-comps/union/Union_mesh.comp b/mcstas-comps/union/Union_mesh.comp index 04b450283..33f6fa38d 100644 --- a/mcstas-comps/union/Union_mesh.comp +++ b/mcstas-comps/union/Union_mesh.comp @@ -130,15 +130,22 @@ SHARE uint16_t attribute_byte_count; } Triangle; #pragma pack(pop) + 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) { 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 @@ -269,6 +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++) { @@ -297,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]; } @@ -308,6 +320,8 @@ SHARE 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); + // TODO: Make tests to verify the contents of FILE int n_lines = 0; char buffer[250]; @@ -432,34 +446,23 @@ 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) { // 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 +477,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; } } @@ -561,86 +564,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 = number_of_points_in_array; + mesh_shell_array.num_elements = n_verts; - 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; } @@ -889,7 +829,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); @@ -971,6 +911,21 @@ 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; @@ -978,39 +933,21 @@ 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)); + 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]]; 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); } 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)); @@ -1061,15 +998,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