/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright(C) 2007 Risoe National Laboratory. * * %I * Written by: Martin Olsen, and expanded upon by Daniel Lomholt Christensen * Date: 17.09.18 * Origin: University of Copenhagen / ACTNXT project * * Component for including 3D mesh in Union geometry * * %D * Part of the Union components, a set of components that work together and thus * separates geometry and physics within McStas. * The use of this component requires other components to be used. * * 1) One specifies a number of processes using process components * 2) These are gathered into material definitions using Union_make_material * 3) Geometries are placed using Union_box/cylinder/sphere/mesh, assigned a material * 4) A Union_master component placed after all the above * * Only in step 4 will any simulation happen, and per default all geometries * defined before this master, but after the previous will be simulated here. * * There is a dedicated manual available for the Union_components * * The mesh component loads a 3D off or stl files as the geometry. * The mesh geometry that is loaded must be a watertight mesh. * In order to check this, the component implements a simple Euler Pointcare value, * To check if the given geometry is watertight. This check is sometimes too rigid, * so a user can set the skip_convex_check parameter in order to not have this * check performed. * * * If you load an off file, we currently only support rank 3 polygons, i.e * triangular meshes. * * It is allowed to overlap components, but it is not allowed to have two * parallel planes that coincides. This will crash the code on run time. * * * * %P * INPUT PARAMETERS: * filename: [str] Name of file that contains the 3D geometry. Can be either .OFF, .off, .STL, or .stl * material_string: [string] material name of this volume, defined using Union_make_material * priority: [1] priotiry of the volume (can not be the same as another volume) A high priority is on top of low. * p_interact: [1] probability to interact with this geometry [0-1] * visualize: [1] set to 0 if you wish to hide this geometry in mcdisplay * number_of_activations: [1] Number of subsequent Union_master components that will simulate this geometry * mask_string: [string] Comma seperated list of geometry names which this geometry should mask * mask_setting: [string] "All" or "Any", should the masked volume be simulated when the ray is in just one mask, or all. * target_index: [1] Focuses on component a component this many steps further in the component sequence * target_x: [m] * target_y: [m] Position of target to focus at * target_z: [m] * focus_aw: [deg] horiz. angular dimension of a rectangular area * focus_ah: [deg] vert. angular dimension of a rectangular area * focus_xw: [m] horiz. dimension of a rectangular area * focus_xh: [m] vert. dimension of a rectangular area * focus_r: [m] focusing on circle with this radius * skip_convex_check: [1] when set to 1, skips the check for whether input .off file is convex. * coordinate_scale [1e-3] Multiplier that the vertex positions are multiplied by. Set to 1 for input coordinates in meters. Set to 1e-3 for input coordinates in mm. * init: [string] name of Union_init component (typically "init", default) * * CALCULATED PARAMETERS: * * %L * * %E ******************************************************************************/ DEFINE COMPONENT Union_mesh SETTING PARAMETERS(string filename = 0, string material_string=0, priority, visualize=1, string all_surfaces=0, string cut_surface=0, int target_index=0, target_x=0, target_y=0, target_z=0, focus_aw=0, focus_ah=0, focus_xw=0, focus_xh=0, focus_r=0, p_interact=0, string mask_string=0, string mask_setting=0, number_of_activations=1, int skip_convex_check = 0, coordinate_scale = 1e-3, string init="init") /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ //============================================================================= // SHARE SECTION //============================================================================= SHARE %{ %include "read_table-lib" %include "interoff-lib" #ifndef Union #error "The Union_init component must be included before this Union_mesh component" #endif #ifndef ANY_GEOMETRY_DETECTOR_DECLARE #define ANY_GEOMETRY_DETECTOR_DECLARE dummy // struct pointer_to_global_geometry_list global_geometry_list = {0,NULL}; #endif #ifndef MAX_VERT_DIST #define MAX_VERT_DIST 1e-30 // Maximum distance between vertex points, before they are the same point. #endif // define struct to import triangles from stl file. #pragma pack(push, 1) // Disable padding typedef struct { float normal[3]; float vertex1[3]; float vertex2[3]; float vertex3[3]; 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); 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 if (fscanf (file, "%5s", word) == 1) { file_is_ascii = strcmp (word, "solid") == 0; } fclose (file); if (file_is_ascii) { FILE* file = Open_File (filename, "r", NULL); char line[256]; int capacity = 100; // initial allocation int count = 0; // Allocate memory for triangles triangles = (Triangle*)malloc (capacity * sizeof (Triangle)); if (triangles == NULL) { free (triangles); printf ("\nERROR: malloc failed for triangles"); exit (1); } Triangle current; int vertex_index = 0; while (fgets (line, sizeof (line), file)) { char* trimmed = line; while (*trimmed == ' ' || *trimmed == '\t') { trimmed++; } // If line is empty after trimming, skip it if (*trimmed == '\0') { continue; } if (strncmp (trimmed, "facet normal", 12) == 0) { memset (¤t, 0, sizeof (Triangle)); sscanf (trimmed, "facet normal %f %f %f", ¤t.normal[0], ¤t.normal[1], ¤t.normal[2]); } else if (strncmp (trimmed, "vertex", 6) == 0) { float x, y, z; sscanf (trimmed, "vertex %f %f %f", &x, &y, &z); if (vertex_index == 0) { current.vertex1[0] = x; current.vertex1[1] = y; current.vertex1[2] = z; } else if (vertex_index == 1) { current.vertex2[0] = x; current.vertex2[1] = y; current.vertex2[2] = z; } else if (vertex_index == 2) { current.vertex3[0] = x; current.vertex3[1] = y; current.vertex3[2] = z; } vertex_index++; if (vertex_index == 3) { vertex_index = 0; } } else if (strncmp (trimmed, "endfacet", 8) == 0) { current.attribute_byte_count = 0; // ASCII STL always 0 if (count >= capacity) { capacity *= 2; void* tmp = realloc (triangles, capacity * sizeof (Triangle)); if (tmp) { free (triangles); free (tmp); perror ("ERROR: Memory reallocation failed"); fclose (file); exit (1); } triangles = tmp; } triangles[count++] = current; } } *n_faces = count; fclose (file); } else { FILE* file = Open_File (filename, "rb", NULL); // Read header char header[80]; fread (header, sizeof (char), 80, file); // Read number of triangles fread (n_faces, sizeof (uint32_t), 1, file); // Allocate memory for triangles triangles = (Triangle*)malloc (*n_faces * sizeof (Triangle)); if (!triangles) { perror ("ERROR: Memory allocation failed"); fclose (file); exit (1); } // Read triangles for (int i = 0; i < *n_faces; i++) { fread (&triangles[i], sizeof (Triangle), 1, file); } fclose (file); } *n_verts = 3 * *n_faces; *verts = (Coords*)malloc (sizeof (Coords) * *n_verts); *faces = (int**)malloc (sizeof (int*) * *n_faces); // Now, we make a list of the triangles, and then give them each an index. int j = 0; for (int i = 0; i < *n_faces; i++) { (*faces)[i] = (int*)malloc (sizeof (int) * 3); j = i * 3; (*verts)[j] = coords_set ((double)triangles[i].vertex1[0], (double)triangles[i].vertex1[1], (double)triangles[i].vertex1[2]); (*verts)[j + 1] = coords_set ((double)triangles[i].vertex2[0], (double)triangles[i].vertex2[1], (double)triangles[i].vertex2[2]); (*verts)[j + 2] = coords_set ((double)triangles[i].vertex3[0], (double)triangles[i].vertex3[1], (double)triangles[i].vertex3[2]); (*faces)[i][0] = j; (*faces)[i][1] = j + 1; (*faces)[i][2] = j + 2; } free (triangles); // Loop over vertices, and make sure we only use unique vertices Coords* unique_vertices = (Coords*)malloc (sizeof (Coords) * *n_verts); 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++) { vertex_is_unique = 1; // First we check if the vertex exist for (j = 0; j < i; j++) { x_are_equal = fabs ((*verts)[i].x - (*verts)[j].x) < MAX_VERT_DIST; y_are_equal = fabs ((*verts)[i].y - (*verts)[j].y) < MAX_VERT_DIST; z_are_equal = fabs ((*verts)[i].z - (*verts)[j].z) < MAX_VERT_DIST; if (x_are_equal && y_are_equal && z_are_equal) { vertex_is_unique = 0; break; } } if (vertex_is_unique) { // Add vertex to unique vertex list map_old_to_unique[j] = unique_counter; unique_vertices[unique_counter] = (*verts)[i]; unique_counter += 1; } // If the vertex isn't unique, we move the faces index, to the // first vertex vert_index_in_faces = floor (i / 3); (*faces)[vert_index_in_faces][i % 3] = map_old_to_unique[j]; } *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]; } free (map_old_to_unique); free (unique_vertices); } 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); // TODO: Make tests to verify the contents of FILE int n_lines = 0; char buffer[250]; while (fgets (buffer, sizeof (buffer), fp)) { if (!strncmp (buffer, "#", 1)) { continue; } if (n_lines == 1) { // in .off files the second line contains faces, vertices and edges sscanf (buffer, "%d %d %d", n_verts, n_faces, n_edges); break; } ++n_lines; } fclose (fp); *verts = (Coords*)malloc (sizeof (Coords) * *n_verts); *faces = (int**)malloc (sizeof (int*) * *n_faces); // Using the Euler-Poincare characteristic, we check if the mesh is closed // To do this, we count the number of edges (unique vertice pairs) fp = Open_File (filename, "r", NULL); n_lines = 0; double x, y, z; int vert_index = 0; int face_vertices = 0; int face_index = 0; while (fgets (buffer, sizeof (buffer), fp)) { // Skip any comments and the first two lines if (!strncmp (buffer, "#", 1)) { continue; } else if (n_lines < 2) { ++n_lines; continue; } // If we are in the vertex area of the file, record the vertex if (n_lines < *n_verts + 2) { // Some .off files have more than three values // these are usually normal vectors or colors. ignore these sscanf (buffer, "%lf %lf %lf", &(*verts)[vert_index].x, &(*verts)[vert_index].y, &(*verts)[vert_index].z); vert_index++; } // If we are in the facet area of the file, record the facet if (n_lines >= (*n_verts + 2)) { (*faces)[face_index] = (int*)malloc (sizeof (int) * 3); // The faces are always after the vertices in an off file. sscanf (buffer, "%d %d %d %d", &face_vertices, &(*faces)[face_index][0], &(*faces)[face_index][1], &(*faces)[face_index][2]); if (face_vertices > 3) { printf ("\nERROR: Facets use more than 3 vertices." "\n Please correct your .off file used for Component %s", comp_name); exit (1); } ++face_index; } ++n_lines; } fclose (fp); } int mesh_is_not_closed (int n_verts, int n_faces, int n_edges) { if (n_verts - n_edges + n_faces == 2) { return 0; } return 1; } void generate_vertex_vertex_pair_list (int** faces, int** vert_pairs, int n_faces) { int vert1, vert2, vert3, vert_pair_0; // Make list of vertex vertex pairs for (int i = 0; i < n_faces; i++) { vert1 = faces[i][0]; vert2 = faces[i][1]; vert3 = faces[i][2]; vert_pairs[i][0] = vert1; vert_pairs[i][1] = vert2; vert_pairs[i + n_faces][0] = vert1; vert_pairs[i + n_faces][1] = vert3; vert_pairs[i + 2 * n_faces][0] = vert2; vert_pairs[i + 2 * n_faces][1] = vert3; } } void find_unique_vertex_vertex_pairs (int* unique_index, int** unique_verts, int** vert_pairs, int* n_faces) { int vert1, vert2; int pair_is_unique; // Make a copy of only the unique pairs for (int i = 0; i < 3 * (*n_faces); i++) { // Check if the first vertex exist in the unique list vert1 = vert_pairs[i][0]; vert2 = vert_pairs[i][1]; pair_is_unique = 1; for (int j = 0; j < (*unique_index); j++) { if (unique_verts[j][0] == vert1) { if (unique_verts[j][1] == vert2) { pair_is_unique = 0; break; } } else if (unique_verts[j][0] == vert2) { if (unique_verts[j][1] == vert1) { pair_is_unique = 0; break; } } } if (pair_is_unique) { unique_verts[*unique_index][0] = vert1; unique_verts[*unique_index][1] = vert2; *unique_index += 1; } } } int coord_comp (Coords A, Coords B) { if (A.x == B.x && A.y == B.y && A.z == B.z) { return 1; } return 0; }; 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; 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 i, j, k; int print1 = 0; int print2 = 0; int print3 = 0; Coords* list_startpoints = (Coords*)calloc (n_facets * 3, sizeof (Coords)); Coords* list_endpoints = (Coords*)calloc (n_facets * 3, sizeof (Coords)); // Check for failed allocation if (list_startpoints == NULL || list_endpoints == NULL) { free (list_startpoints); free (list_endpoints); } int counter = 0; // For every triangle it should add three lines 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 (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[j])) { for (k = 0; k < counter; k++) { if (coord_comp (point1, list_startpoints[j])) { print1 = 0; } } } 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[j])) { for (k = 0; k < counter; k++) { if (coord_comp (point2, list_startpoints[j])) { print2 = 0; } } } 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[j])) { for (k = 0; k < counter; k++) { if (coord_comp (point1, list_startpoints[j])) { print3 = 0; } } } } // Create lines // Line 1 if (print1 == 1) { lines_to_draw_temp = draw_line_with_highest_priority (point1, point2, index, Geometries, number_of_volumes, 100); merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp); list_startpoints[counter] = point1; list_endpoints[counter] = point2; counter++; } // Line 2 if (print2 == 1) { lines_to_draw_temp = draw_line_with_highest_priority (point2, point3, index, Geometries, number_of_volumes, 100); merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp); list_startpoints[counter] = point2; list_endpoints[counter] = point3; counter++; } // Line 3 if (print3 == 1) { lines_to_draw_temp = draw_line_with_highest_priority (point3, point1, index, Geometries, number_of_volumes, 100); merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp); list_startpoints[counter] = point3; list_endpoints[counter] = point1; counter++; } } free (list_startpoints); free (list_endpoints); }; struct pointer_to_1d_coords_list allocate_shell_points (struct geometry_struct* geometry, int max_number_of_points) { // Function that returns all vertex points transported to the center of the geometry 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; 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_add (mesh_shell_array.elements[i], geometry->center); // Transpose and rotate the points such that they are in the right coordinate system } mesh_shell_array.num_elements = n_verts; return mesh_shell_array; } void initialize_mesh_geometry_from_main_component (struct geometry_struct* mesh) { // Function to be called in initialize of the main component // This is done as the rotation matrix needs to be relative to the main component instead of global // Everything done in initialize in this component file has the rotation matrix relative to global Coords simple_vector; Coords mesh_vector; // Start with vector that points along the mesh in the local frame simple_vector = coords_set (0, 1, 0); // Rotate the direction vector of the mesh to the master component frame of reference mesh_vector = rot_apply (mesh->rotation_matrix, simple_vector); NORM (mesh_vector.x, mesh_vector.y, mesh_vector.z); mesh->geometry_parameters.p_mesh_storage->direction_vector.x = mesh_vector.x; mesh->geometry_parameters.p_mesh_storage->direction_vector.y = mesh_vector.y; mesh->geometry_parameters.p_mesh_storage->direction_vector.z = mesh_vector.z; mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center = rot_apply (mesh->rotation_matrix, mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center); /* // Works for pure translation print_position(mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center, "BB before adjustment"); mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center = coords_add(mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center, mesh->center); print_position(mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center, "BB after adjustment"); */ } double square (double x) { return x * x; } %} //============================================================================= // DECLARE SECTION //============================================================================= DECLARE %{ // Needed for transport to the main component struct global_geometry_element_struct global_geometry_element; int loop_index; int loop_2_index; int material_index; off_struct offdata; // off struct isnt currently in use // Volume struct contains name, geometry struct, physics struct // loggers struct and absorption loggers struct. struct Volume_struct mesh_vol; // Mesh storage contains the number of facets, and for each facet // the three vertexes coordinates. Not in an optimal way. // However, the Bounding box center and radius is also stored here. struct mesh_storage mesh_data; // Structs for surface effects struct surface_stack_struct surface_stack; struct surface_stack_struct cut_surface_stack; %} //============================================================================= // INITIALIZE SECTION //============================================================================= INITIALIZE %{ geometry_struct_init (&(mesh_vol.geometry)); mesh_vol.geometry.skip_hierarchy_optimization = skip_convex_check; // Initializes the focusing system for this volume including input sanitation. focus_initialize (&mesh_vol.geometry, POS_A_COMP_INDEX (INDEX_CURRENT_COMP + target_index), POS_A_CURRENT_COMP, ROT_A_CURRENT_COMP, target_index, target_x, target_y, target_z, focus_aw, focus_ah, focus_xw, focus_xh, focus_r, NAME_CURRENT_COMP); if (_getcomp_index (init) < 0) { fprintf (stderr, "Union_mesh:%s: Error identifying Union_init component," "%s is not a known component name.\n", NAME_CURRENT_COMP, init); exit (-1); } struct pointer_to_global_material_list* global_material_list = COMP_GETPAR3 (Union_init, init, global_material_list); // Use sanitation #ifdef MATERIAL_DETECTOR if (global_material_list->num_elements == 0) { // Here if the user have defined a material, but only after this material printf ("\nERROR: Need to define a material using Union_make_material" " before using a Union geometry component. \n"); printf ("%s was defined before first use of Union_make_material.\n", NAME_CURRENT_COMP); exit (1); } #endif #ifndef MATERIAL_DETECTOR printf ("\nERROR: Need to define a material using Union_make_material" " before using a Union geometry component. \n"); exit (1); #endif mesh_vol.geometry.is_masked_volume = 0; mesh_vol.geometry.is_exit_volume = 0; mesh_vol.geometry.is_mask_volume = 0; struct pointer_to_global_geometry_list* global_geometry_list = COMP_GETPAR3 (Union_init, init, global_geometry_list); //============================================================================== // Find out what this component masks //============================================================================== // Read the material input, or if it lacks, use automatic linking. if (mask_string && strlen (mask_string) && strcmp (mask_string, "NULL") && strcmp (mask_string, "0")) { // A mask volume is used to limit the extend of other volumes, // called the masked volumes. These are specified in the mask_string. // In order for a ray to enter a masked volume, // it needs to be both in the region covered by that volume AND the mask volume. // When more than mesh_vol.geometry.mask_mode = 1; // Default mask mode is ALL if (mask_setting && strlen (mask_setting) && strcmp (mask_setting, "NULL") && strcmp (mask_setting, "0")) { if (strcmp (mask_setting, "ALL") == 0 || strcmp (mask_setting, "All") == 0) { mesh_vol.geometry.mask_mode = 1; } else if (strcmp (mask_setting, "ANY") == 0 || strcmp (mask_setting, "Any") == 0) { mesh_vol.geometry.mask_mode = 2; } else { printf ("The mask_mode of component %s is set to %s," " but must be either ALL or ANY.\n", NAME_CURRENT_COMP, mask_setting); exit (1); } } int found_geometries = 0; for (loop_index = 0; loop_index < global_geometry_list->num_elements; loop_index++) { // Add mask list if (1 == manual_linking_function (global_geometry_list->elements[loop_index].name, mask_string)) { add_element_to_int_list (&mesh_vol.geometry.mask_list, global_geometry_list->elements[loop_index].component_index); add_element_to_int_list (&global_geometry_list->elements[loop_index].Volume->geometry.masked_by_list, INDEX_CURRENT_COMP); global_geometry_list->elements[loop_index].Volume->geometry.is_masked_volume = 1; if (mesh_vol.geometry.mask_mode == 2) global_geometry_list->elements[loop_index].Volume->geometry.mask_mode = 2; if (mesh_vol.geometry.mask_mode == 1) { if (global_geometry_list->elements[loop_index].Volume->geometry.is_masked_volume == 1 && global_geometry_list->elements[loop_index].Volume->geometry.mask_mode != 2) // If more than one mask is added to one volume, the ANY mode overwrites the (default) ALL mode. global_geometry_list->elements[loop_index].Volume->geometry.mask_mode = 1; } found_geometries = 1; } } if (found_geometries == 0) { printf ("The mask_string in geometry: " "%s did not find any of the specified volumes in the mask_string " "%s \n", NAME_CURRENT_COMP, mask_string); exit (1); } mesh_vol.p_physics = malloc (sizeof (struct physics_struct)); mesh_vol.p_physics->is_vacuum = 0; // Makes this volume a vacuum mesh_vol.p_physics->number_of_processes = (int)0; // Should not be used. mesh_vol.p_physics->my_a = 0; // Should not be used. sprintf (mesh_vol.p_physics->name, "Mask"); mesh_vol.geometry.is_mask_volume = 1; // Read the material input, or if it lacks, use automatic linking. } else if (material_string && strlen (material_string) && strcmp (material_string, "NULL") && strcmp (material_string, "0")) { // A geometry string was given, use it to determine which material if (0 == strcmp (material_string, "vacuum") || 0 == strcmp (material_string, "Vacuum")) { // One could have a global physics struct for vacuum instead of creating one for each mesh_vol.p_physics = malloc (sizeof (struct physics_struct)); mesh_vol.p_physics->is_vacuum = 1; // Makes this volume a vacuum mesh_vol.p_physics->number_of_processes = (int)0; mesh_vol.p_physics->my_a = 0; // Should not be used. sprintf (mesh_vol.p_physics->name, "Vacuum"); } else if (0 == strcmp (material_string, "exit") || 0 == strcmp (material_string, "Exit")) { // One could have a global physics struct for exit instead of creating one for each mesh_vol.p_physics = malloc (sizeof (struct physics_struct)); mesh_vol.p_physics->is_vacuum = 1; // Makes this volume a vacuum mesh_vol.p_physics->number_of_processes = (int)0; mesh_vol.p_physics->my_a = 0; // Should not be used. mesh_vol.geometry.is_exit_volume = 1; sprintf (mesh_vol.p_physics->name, "Exit"); } else { for (loop_index = 0; loop_index < global_material_list->num_elements; loop_index++) { if (0 == strcmp (material_string, global_material_list->elements[loop_index].name)) { mesh_vol.p_physics = global_material_list->elements[loop_index].physics; break; } if (loop_index == global_material_list->num_elements - 1) { printf ("\n"); printf ("ERROR: The material string \"%s\" in Union geometry \"%s\"" " did not match a specified material. \n", material_string, NAME_CURRENT_COMP); printf (" The materials available at this point" " (need to be defined before the geometry): \n"); for (loop_index = 0; loop_index < global_material_list->num_elements; loop_index++) printf (" %s\n", global_material_list->elements[loop_index].name); printf ("\n"); printf (" It is also possible to use one of the defualt materials avaiable: \n"); printf (" Vacuum (for a Volume without scattering or absorption)\n"); printf (" Exit (for a Volume where the ray exits the component if it enters)\n"); printf (" Mask (for a Volume that masks existing volumes specified in the mask_string\n"); exit (1); } } } } else { // Automatic linking, simply using the last defined material. #ifndef MATERIAL_DETECTOR printf ("Need to define a material before the geometry to use automatic linking %s.\n", NAME_CURRENT_COMP); exit (1); #endif mesh_vol.p_physics = global_material_list->elements[global_material_list->num_elements - 1].physics; } //============================================================================== // READ FIle //============================================================================== // Read input file and put into storage. // We assume that the input file is using triangular faces, otherwise // the check for closed mesh using the Euler-Poincare characteristic // is wrong. int n_verts = 0; int n_faces = 0; int n_edges = 0; Coords* verts; Coords tmp_vert = coords_set (0, 0, 0); verts = &tmp_vert; // This tmp vert is just for linter play... int** faces; int unique_index = 0; const char* dot = strrchr (filename, '.'); if (!dot || dot == filename) { printf ("ERROR: given file has no extension %s", NAME_CURRENT_COMP); exit (1); } if (strcmp (dot, ".stl") == 0 || strcmp (dot, ".STL") == 0) { printf ("\n%s, Given file is in stl format. Union_mesh assumes three vertices per face.\n", NAME_CURRENT_COMP); read_stl (filename, &n_verts, &n_faces, &n_edges, &verts, &faces, NAME_CURRENT_COMP); } else if (strcmp (dot, ".off") == 0 || strcmp (dot, ".OFF") == 0) { printf ("\n%s, Given file is in off format. Union_mesh assumes three vertices per face.\n", NAME_CURRENT_COMP); read_off_file_vertices_and_faces (filename, &n_verts, &n_faces, &n_edges, &verts, &faces, NAME_CURRENT_COMP); } else { 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_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); } // Make lists that will be used to calculate the number of edges int** vert_pairs = (int**)malloc (sizeof (int*) * 3 * n_faces); int** unique_verts = (int**)malloc (sizeof (int*) * 3 * n_faces); // Check for failing mallocs if (vert_pairs == NULL || unique_verts == NULL) { free (vert_pairs); free (unique_verts); printf ("ERROR: malloc failed on vert pairs or unique verts"); exit (1); } for (int i = 0; i < 3 * n_faces; i++) { vert_pairs[i] = (int*)malloc (sizeof (int) * 2); unique_verts[i] = (int*)malloc (sizeof (int) * 2); if (vert_pairs[i] == NULL || unique_verts[i] == NULL) { // Should probably be a loop over all vert pairs free (vert_pairs[i]); free (unique_verts[i]); printf ("\nERROR: malloc failed on specific vert pairs or unique verts"); exit (1); } } generate_vertex_vertex_pair_list (faces, vert_pairs, n_faces); find_unique_vertex_vertex_pairs (&unique_index, unique_verts, vert_pairs, &n_faces); n_edges = unique_index; // The number of edges is the number of unique // vertex vertex pairs. if (mesh_is_not_closed (n_faces, n_verts, n_edges) && !skip_convex_check) { printf ("\nNumber of edges is the unique index, %d", unique_index); printf ("\nEuler Pointcare characteristic is %d", n_faces + n_verts - unique_index); printf ("\n\nERROR: Input mesh from %s, is not a closed shape according to Euler Pointcare." "\nNumber of faces = %d, Number of vertices = %d, Number of Edges = %d" "\nUnion_mesh_off requires the user input to be a convex hull.\n" "If you are sure that it is, set skip_convex_check to 1\n\n", NAME_CURRENT_COMP, n_faces, n_verts, n_edges); exit (1); } // Calculate bounding box double max_dist = 0; double dist = 0; if (n_verts <= 0) { printf ("\nERROR: Number of vertices read is 0\n"); exit (1); } Coords B_sphere_x = verts[0]; Coords B_sphere_y = B_sphere_x; for (int i = 0; i < n_verts; i++) { dist = sqrt (B_sphere_x.x - verts[i].x + B_sphere_x.y - verts[i].y + B_sphere_x.z - verts[i].z); if (dist > max_dist) { max_dist = dist; B_sphere_y = verts[i]; } } Coords B_sphere_z = B_sphere_y; max_dist = 0; for (int i = 0; i < n_verts; i++) { dist = sqrt (square (B_sphere_y.x - verts[i].x) + square (B_sphere_y.y - verts[i].y) + square (B_sphere_y.z - verts[i].z)); if (dist > max_dist) { max_dist = dist; B_sphere_z = verts[i]; } } double radius = sqrt (square (B_sphere_y.x - B_sphere_z.x) + square (B_sphere_y.y - B_sphere_z.y) + square (B_sphere_y.z - B_sphere_z.z)) / 2; Coords bbcenter = coords_set ((B_sphere_y.x + B_sphere_z.x) / 2, (B_sphere_y.y + B_sphere_z.y) / 2, (B_sphere_y.z + B_sphere_z.z) / 2); mesh_data.Bounding_Box_Center = bbcenter; double test_rad = 0; for (int i = 0; i < n_verts; i++) { test_rad = sqrt (square (bbcenter.x - verts[i].x) + square (bbcenter.y - verts[i].y) + square (bbcenter.z - verts[i].z)); 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; Coords vert1, vert2, vert3; 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.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]]; 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.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)); mesh_vol.geometry.surface_stack_for_each_face[0] = &surface_stack; mesh_vol.geometry.internal_cut_surface_stack = &cut_surface_stack; struct pointer_to_global_surface_list* global_surface_list = COMP_GETPAR3 (Union_init, init, global_surface_list); fill_surface_stack (all_surfaces, global_surface_list, NAME_CURRENT_COMP, &surface_stack); fill_surface_stack (cut_surface, global_surface_list, NAME_CURRENT_COMP, &cut_surface_stack); sprintf (mesh_vol.name, "%s", NAME_CURRENT_COMP); sprintf (mesh_vol.geometry.shape, "mesh"); mesh_vol.geometry.eShape = mesh; mesh_vol.geometry.priority_value = priority; mesh_vol.geometry.visualization_on = visualize; mesh_vol.geometry.geometry_p_interact = p_interact; mesh_vol.geometry.within_function = &r_within_mesh; mesh_vol.geometry.geometry_parameters.p_mesh_storage = &mesh_data; mesh_vol.geometry.intersect_function = &sample_mesh_intersect; mesh_vol.geometry.mcdisplay_function = &mcdisplay_mesh_function; mesh_vol.geometry.shell_points = &allocate_shell_points; mesh_vol.geometry.initialize_from_main_function = &initialize_mesh_geometry_from_main_component; mesh_vol.geometry.process_rot_allocated = 0; mesh_vol.geometry.copy_geometry_parameters = &allocate_mesh_storage_copy; mesh_vol.geometry.center = POS_A_CURRENT_COMP; rot_copy (mesh_vol.geometry.rotation_matrix, ROT_A_CURRENT_COMP); rot_transpose (ROT_A_CURRENT_COMP, mesh_vol.geometry.transpose_rotation_matrix); // Initialize loggers mesh_vol.loggers.num_elements = 0; mesh_vol.abs_loggers.num_elements = 0; // packing the information into the global_geometry_element, which is then included in the global_geometry_list. sprintf (global_geometry_element.name, "%s", NAME_CURRENT_COMP); global_geometry_element.activation_counter = number_of_activations; global_geometry_element.component_index = INDEX_CURRENT_COMP; global_geometry_element.Volume = &mesh_vol; // Would be nicer if this m was a pointer, now we have the (small) data two places add_element_to_geometry_list (global_geometry_list, global_geometry_element); %} TRACE %{ %} FINALLY %{ %} END