This repository serve as a backup for my Maxwell-TD code
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 
 

545 lines
18 KiB

import mesh_utils as mu
#######################################################
## ##
## 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 = 0.01
########################################
## MAIN FUNCTION ##
########################################
# 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]
'''
# Find offset to snap vertices to a grid
regVertices, offset_trunc = snap_vertices_relative_to_center(regVertices, referenceNode, regularElementSize, unit2trunc)
reg_min = np.min(regVertices, axis=0)
reg_max = np.max(regVertices, axis=0)
print("Min Regular Corner after snapping = ", reg_min)
print("Max Regular Corner after snapping = ", reg_max)
print("Snapped regVertices to grid")
'''
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_pml( elementArray=combined_element_array, bcArray=combined_bc_array, totalNodes=totalNodes,
materialArray=temp_material_flags, region_corner=reg_corner , region_dims=reg_dims ,
excluded_bc_values=[2,3], 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")
sigma_base = 0.001
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")