from mesh_utils import * ####################################################### ####################################################### ## ## ## PART TO MODIFY ## ## ## ####################################################### ####################################################### directory = "./" # Absolute path to the model files fileName = "PW" # Name of the .g and .json files unit = 1.0 # Unit of the input (in meters) -> 1.0 for meters or 2.54e-2 for inches mesh_size = 0.02 # Mesh size in the unit used unit2trunc = 4 # Number of significant digits in the regular mesh number_of_layers = 5 # Number of PML layers to generate sigma_base = 0.005 # Setup full file path and mesh size name = f"{directory}/{fileName}" regularElementSize = truncate_to_significant_figures(mesh_size, 4) # Load and configure volume mesh mesh = geo.VolumeMesh.from_exodus(f"{name}.g") config = mesh.apply_config(f"{name}.json") BC = geo.BoundaryCondition print_boundary_conditions() # Check PEC boundaries pec_faces = mesh.face_sets_with_bc(BC.PEC) if len(pec_faces) == 0: print("Error: No PEC faces found.") print("PEC face set size:", len(pec_faces)) ########################################### ## Generate Regular Airbox Mesh ## ########################################### regular = mesh.make_Airbox_pml_grid_mesh(regularElementSize, pec_faces) regVertices = np.array(regular.vertices) #print("Regular Vertices = ", regular.vertices) print("------------------------------------------------") print("Shape of regular vertices = ", regVertices.shape) print("------------------------------------------------") reg_min = np.min(regVertices, axis=0) reg_max = np.max(regVertices, axis=0) print("########################################") print(" Regular Grid Bounding Box") print("########################################") print(f"X range: {reg_min[0]:.6f} to {reg_max[0]:.6f}") print(f"Y range: {reg_min[1]:.6f} to {reg_max[1]:.6f}") print(f"Z range: {reg_min[2]:.6f} to {reg_max[2]:.6f}") print(f"Box Size (dx, dy, dz): {reg_max - reg_min}") print("########################################") # First vertex is used as the reference node referenceNode = regular.vertices[:, 0] print("First Node from regular tets = ", referenceNode) regularVertices = regular.vertices.shape[0] export_bc_surfaces(regular, "Regular", BC) ######################################## ## Generate Irregular Mesh ## ######################################## irregular = mesh.make_hybrid_mesh_Airbox( regularElementSize, pec_faces, max_volume=(regularElementSize ** 3) * (np.sqrt(2) / 12) ) print("------------------------------------------------") #print("Irregular Vertices = ", irregular.vertices) print("Irregular Vertices Shape = ", irregular.vertices.shape) print("------------------------------------------------") export_bc_surfaces(irregular, "Irregular", BC) ######################################## ## Generate Box-Shaped PML ## ######################################## # Get bounding box of the regular mesh reg_min = np.min(regVertices, axis=1) reg_max = np.max(regVertices, axis=1) print("Min Corner = ", reg_min) print("Max Corner = ", reg_max) # ---- Inner PML (existing) ---- sc_layer, sc_vert = generate_cuboid_pml_layer( regular, regularElementSize, reg_min, # inner (core) min reg_max, # inner (core) max layer_thickness=1 ) export_bc_surfaces(sc_layer, "sc_layer", BC) # Collect min/max of inner PML (pml_vert is (N,3)) pml1_min = np.min(sc_vert, axis=1) pml1_max = np.max(sc_vert, axis=1) all_sc_vertices = np.vstack(sc_vert) if len(sc_vert) > 0 else np.empty((0, 3)) # ---- Outer PML (new) ---- # thickness of the outer shell (cells). Feel free to change: outer_number_of_layers = number_of_layers pml_layer, pml_vert = generate_cuboid_pml_layer( regular, regularElementSize, pml1_min, # inner hole of outer PML is exactly the outer bounds of the inner PML pml1_max, layer_thickness=outer_number_of_layers ) export_bc_surfaces(pml_layer, "pml_layer", BC) all_pml_vertices = np.vstack(pml_vert) if len(pml_vert) > 0 else np.empty((0, 3)) pml_min = np.min(all_pml_vertices, axis=0) pml_max = np.max(all_pml_vertices, axis=0) export_bc_surfaces(pml_layer, "pml", BC) # ------------------------------- # Deduplicate and export node list # The unique set of nodes # The indices of the first occurrences # A mapping array from original to deduplicated nodes # The count of how many times each unique row occurred (but unused here) # ------------------------------- print("Shape of all_pml_vertices:", (pml_vert).T.shape) print("Shape of all_sc_vertices:", (sc_vert).T.shape) print("Shape of regular vertices:", np.array(regular.vertices).T.shape) print("Shape of irregular vertices:", np.array(irregular.vertices).T.shape) # Stack all vertex coordinates for deduplication totalNodes = np.vstack(( np.array(irregular.vertices).T, np.array(regular.vertices).T, np.array(sc_vert).T, np.array(pml_vert).T )) print("🔢 Number of nodes before rounding/deduplication:", totalNodes.shape[0]) print("TotalNodes = ",totalNodes.shape) # 2. Round coordinates to a precision (e.g., 1e-8) precision = 8 totalNodes = np.round(totalNodes, decimals=precision) # Deduplicate and get mapping to global IDs totalNodes, indices, verticesMap, _ = np.unique(totalNodes, axis=0, return_index=True, return_inverse=True, return_counts=True) print("✅ Number of unique nodes after rounding:", totalNodes.shape[0]) # Save to .node file np.savetxt(f'{name}.node', totalNodes, newline="\n3 0 0\n" , fmt='%f', header=f"{unit}\n{totalNodes.shape[0]}", comments='') print("Total deduplicated nodes:", totalNodes.shape[0]) ######################################## ## IRREGULAR ELEMENTS ## ######################################## offset = 0 # Build face BC array for irregular mesh selected_BC = [BC.PEC] output_BC_number = [1] irregularFacesBC = build_face_bc_array_custom(irregular, BC, selected_BC, output_BC_number) # Build element and face BC arrays for irregular mesh irregularElementArray, irregularBcArray = build_element_and_bc_arrays( irregular, irregularFacesBC, verticesMap, offset=offset) print("Unique face BC values in irregularFacesBC:") print(np.unique(irregularFacesBC)) offset += irregular.vertices.shape[1] ######################################## ## REGULAR ELEMENTS ## ######################################## ### This section processes regular mesh elements by assigning boundary conditions (BCs), ### building element connectivity arrays, and canonically ordering nodes for consistency. selected_BC = [BC.PML_Excitation] output_BC_number = [2] regularFacesBC = build_face_bc_array_custom(regular, BC, selected_BC, output_BC_number) elementArray, bcArray = build_element_and_bc_arrays(regular, regularFacesBC, verticesMap,offset=offset) offset += regular.vertices.shape[1] print("Unique face BC values in regularFacesBC:") print(np.unique(regularFacesBC)) print("Shape of elementArray:", elementArray.shape) print("Max vertex index used:", np.max(elementArray)) surfs = irregular.extract_surfaces(irregular.face_sets_with_bc(BC.PEC)) surfs.save_stl(f"{name}_regular_pec_surfaces.stl") surfs = regular.extract_surfaces(regular.face_sets_with_bc(BC.NCB)) surfs.save_stl(f"{name}_regular_ncbc_surfaces.stl") surfs = regular.extract_surfaces(regular.face_sets_with_bc(BC.ABC)) surfs.save_stl(f"{name}_regular_abc_surfaces.stl") ################################### ## SC ELEMENTS ## ################################### # Accumulators for all PML elements and BCs selected_BC = [BC.PML_Excitation,BC.PML_ABC] output_BC_number = [2,0] scFacesBC = build_face_bc_array_custom(sc_layer, BC, selected_BC, output_BC_number) scElementArrays, scBcArrays = build_element_and_bc_arrays(sc_layer, scFacesBC, verticesMap, offset = offset) offset += sc_layer.vertices.shape[1] print("Unique face BC values in scFacesBC:") print(np.unique(scFacesBC)) print("SC Shape of elementArray:", scElementArrays.shape) print("SC Max vertex index used:", np.max(scElementArrays)) #################################### ## PML ELEMENTS ## #################################### # Accumulators for all PML elements and BCs selected_BC = [BC.PML_Excitation,BC.PML_ABC] output_BC_number = [0,3] pmlFacesBC = build_face_bc_array_custom(pml_layer, BC, selected_BC, output_BC_number) pmlElementArrays, pmlBcArrays = build_element_and_bc_arrays(pml_layer, pmlFacesBC, verticesMap, offset = offset) offset += pml_layer.vertices.shape[1] print("Unique face BC values in pmlFacesBC:") print(np.unique(pmlFacesBC)) print("PML Shape of elementArray:", pmlElementArrays.shape) print("PML Max vertex index used:", np.max(pmlElementArrays)) ######################################## ## BUILD FULL ARRAYS ## ######################################## # Concatenate tetrahedral elements in order: PML, regular, irregular full_element_array = np.vstack([irregularElementArray, elementArray, scElementArrays, pmlElementArrays]) # Total element counts Nirreg = irregularElementArray.shape[0] Nreg = elementArray.shape[0] Nsc = scElementArrays.shape[0] Npml = pmlElementArrays.shape[0] Ntotal = Nsc + Npml + Nreg + Nirreg assert full_element_array.shape[0] == Ntotal ######################################## ## MATERIAL ASSIGNMENT ## ######################################## material_flags = np.zeros(Ntotal, dtype=np.int64) # SC: material 1 material_flags[Nirreg+Nreg:Nirreg+Nreg+Nsc] = 1 # PML: material 2 material_flags[Nirreg+Nreg+Nsc::] = 2 ######################################## ## BUILD BC ARRAY ## ######################################## # Ensure BC arrays match the stacked order bc_full_array = np.vstack([ irregularBcArray, # shape: (Nirreg, 4) bcArray, # regular, shape: (Nreg, 4) scBcArrays, # shape: (Npml, 4) pmlBcArrays # shape: (Nsc, 4) ]) assert bc_full_array.shape == (Ntotal, 4) print("Unique BC codes in full BC array:",np.unique(bc_full_array)) ################################# ## GROUPING REGULARS ## ################################# # Concatenate tetrahedral elements in order: PML, regular, irregular combined_element_array = np.vstack([elementArray, scElementArrays, pmlElementArrays]) combined_bc_array = np.vstack([ bcArray, # regular, shape: (Nreg, 4) scFacesBC, pmlBcArrays # shape: (Npml, 4) ]) temp_material_flags = np.zeros(Nreg+Nsc+Npml, dtype=np.int64) temp_material_flags[Nreg:Nreg+Nsc] = 1 temp_material_flags[Nreg+Nsc:Nreg+Nsc+Npml] = 2 print("Unique bc IDS = ", np.unique(combined_bc_array)) print("==================================================") print("Grouping regular elements by geometry...") print("==================================================") # Bounding box of regular airbox minX = np.min(regVertices[0, :]) maxX = np.max(regVertices[0, :]) minY = np.min(regVertices[1, :]) maxY = np.max(regVertices[1, :]) minZ = np.min(regVertices[2, :]) maxZ = np.max(regVertices[2, :]) reg_corner = np.array([minX, minY, minZ]) reg_dims = np.array([maxX - minX, maxY - minY, maxZ - minZ]) print("Bounding box for regular region grouping:") print(reg_corner, reg_corner + reg_dims) Excitation_minZ = minZ Excitation_maxZ = maxZ Excitation_minX = minX Excitation_maxX = maxX Excitation_minY = minY Excitation_maxY = maxY # Bounding box of SC layer minX = np.min(all_sc_vertices[0, :]) maxX = np.max(all_sc_vertices[0, :]) minY = np.min(all_sc_vertices[1, :]) maxY = np.max(all_sc_vertices[1, :]) minZ = np.min(all_sc_vertices[2, :]) maxZ = np.max(all_sc_vertices[2, :]) reg_corner = np.array([minX, minY, minZ]) reg_dims = np.array([maxX - minX, maxY - minY, maxZ - minZ]) # Bounding box of PML layer PMLminX = np.min(all_pml_vertices[0, :]) thickness = abs(PMLminX - minX) pml_bounds=(thickness,thickness,thickness) reg_indices = group_tets_by_box_and_conductivity_zero_bc_only( elementArray=combined_element_array, bcArray=combined_bc_array, totalNodes=totalNodes, materialArray=temp_material_flags, region_corner=reg_corner , region_dims=reg_dims , pml_bounds=pml_bounds, min_group_size=100) print("=====================================") vals, counts = np.unique(reg_indices, return_counts=True) print("Counts (including 0):") for v, c in zip(vals, counts): print(f"{int(v)}: {int(c)}") print("=====================================") ################################# ## CONSTRUCT FINAL INDEX ## ################################# print("==================================================") print("Constructing overall global regular indices ...") print("==================================================") # Combine all tets and group indices full_element_array = np.vstack([irregularElementArray, elementArray, scElementArrays, pmlElementArrays]) final_indices = np.hstack([ np.zeros(irregularElementArray.shape[0], dtype=np.int64), # Irregular reg_indices # Regular ]) print("=====================================") vals, counts = np.unique(final_indices, return_counts=True) print("Final Counts (including 0):") for v, c in zip(vals, counts): print(f"{int(v)}: {int(c)}") print("=====================================") print("final_indices shape:", final_indices.shape) print("Max group ID = ", np.sort(np.unique(final_indices))) assert final_indices.shape[0] == full_element_array.shape[0], \ "Mismatch between final_indices and total element count!" # Save .regular group file (now matches .tetra) total_group_count = np.max(final_indices) + 1 np.savetxt(f"{name}.regular", final_indices, fmt='%d', header=str(total_group_count), comments='') print(f"Saved {final_indices.shape[0]} group indices to {name}.regular") ######################################## ## WRITE TO .tetra FILE ## ######################################## c_with_extra = np.ones((2 * Ntotal, 5), dtype=np.int64) * (-1) # Even rows: material + connectivity c_with_extra[::2] = np.hstack(( material_flags.reshape(-1, 1), full_element_array )) # Odd rows: face BCs c_with_extra[1::2, 1:] = bc_full_array # Format as DataFrame for easier output df = pd.DataFrame(c_with_extra, dtype=np.int64).where(lambda x: x != -1, np.nan) # Insert element count header header_row = pd.Series([Ntotal] + [np.nan] * 4, index=df.columns) df = pd.concat([header_row.to_frame().T, df], ignore_index=True) # Save as .tetra file df = df.astype("Int64") df.to_csv(f"{name}.tetra", sep=' ', index=False, header=False) print(f"✅ Saved {Ntotal} tetrahedra to {name}.tetra") ################################# ## WRITE FILES ## ################################# print("Writing freespace.m...") with open("freespace.m", "w") as f: # First 6 rows: identity tensors f.write("standard\n") f.write("1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00\n") f.write("1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00\n") # 3 rows of zero vectors (dispersion parameters) for _ in range(3): f.write("0.000000e+00 0.000000e+00 0.000000e+00\n") print("Writing scaterring.m...") with open("scattering.m", "w") as f: # First 6 rows: identity tensors f.write("standard\n") f.write("1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00\n") f.write("1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00\n") # 3 rows of zero vectors (dispersion parameters) for _ in range(3): f.write("0.000000e+00 0.000000e+00 0.000000e+00\n") print("Writing .bc file...") with open(f"{name}.bc", "w") as f: f.write("4\n") f.write("0 none\n") f.write("1 pec\n") f.write(f"2 pml 1 pmlExci 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 {Excitation_minZ} \n") f.write("3 abc 1\n") with open(f"{name}.prop", "w") as f: f.write("./\n") f.write("freespace\n") f.write("scattering\n") f.write("pml\n") print("Writing pml.m...") with open("pml.m", "w") as f: f.write("standard\n") f.write("1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00\n") f.write("1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00\n") f.write("0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00\n") f.write(f"{sigma_base:.6e} 0.000000e+00 0.000000e+00\n") f.write(f"0.000000e+00 {sigma_base:.6e} 0.000000e+00\n") f.write(f"0.000000e+00 0.000000e+00 {sigma_base:.6e}\n") f.write(f"3 {thickness:.3f} -1.0 -1.0 -1.0\n") # Use expanded inner cavity bounds (aligned to grid) f.write(f"{pml1_min[0]:.3f} {pml1_max[0]:.3f} " f"{pml1_min[1]:.3f} {pml1_max[1]:.3f} " f"{pml1_min[2]:.3f} {pml1_max[2]:.3f}\n")