import sys sys.path.append('/media/digitalstorm/edf2449b-a566-40a2-8d0f-393bbe8ea1c0/maxwell_td_hybrid_mesher/Maxwell_PML_Regular/HybridMesher') import Geometry as geo # Custom mesh/geometry handling module import numpy as np import pandas as pd from numpy import linalg as la from Geometry import Vector # pybind11 C++ binding from scipy.spatial import cKDTree ####################################################### # CONSTANTS # ####################################################### PEC = 1 ABC_PW = 2 ABC = 3 Interface_PW = 4 NCBC = 1000 unit2trunc = 4 # Number of significant digits in the regular mesh ######################################## ## AUXILIARY FUNCTIONS ## ######################################## def floor(number, sig_figs): return np.trunc(number * (10 ** sig_figs)) / (10 ** sig_figs) # Handles scalars or arrays and truncates each to the desired significant figures def truncate_to_significant_figures(number, sig_figs): isNumber = False if not isinstance(number, np.ndarray): number = np.array([number]) isNumber = True maskLower = number < 0 maskGreater = number > 0 order_of_magnitude = np.zeros_like(number) order_of_magnitude[maskGreater] = np.floor(np.log10(number[maskGreater])) order_of_magnitude[maskLower] = np.floor(np.log10(-number[maskLower])) if isNumber: return floor(number, -(order_of_magnitude - sig_figs + 1))[0] return floor(number, -(order_of_magnitude - sig_figs + 1)) def build_face_bc_array(mesh, BC, PEC, ABC, NCBC): number_mesh_elems = mesh.info.num_elems number_mesh_faces = mesh.info.num_face_sets faceBC = np.zeros((number_mesh_elems, 4), dtype=np.int64) print(number_mesh_faces, " face sets found in the mesh.") # Loops over each face set in the regular mesh (each face set groups faces sharing the same BC). # Assign face BCs for i in range(number_mesh_faces): bc_val = 0 bc_type = mesh.get_face_set_bc(i) if bc_type == BC.PEC: bc_val = PEC elif bc_type == BC.ABC: bc_val = ABC elif bc_type == BC.NCB: bc_val = NCBC # Each face in the set is identified by a global face ID. # faceID // 4 gives the element index. # faceID % 4 gives the local face number (0–3). # This maps the BC code to the correct element face in regularFacesBC. for face_id in mesh.face_set(i): elem_idx = face_id // 4 local_face = face_id % 4 faceBC[elem_idx][local_face] = bc_val return faceBC def build_element_and_bc_arrays(mesh, faceBC, verticesMap, offset=0): # Pre-allocate space for all tetrahedral elements and face BCs # Each tetrahedron has 4 vertices and 4 faces number_mesh_elems = mesh.info.num_elems number_mesh_faces = mesh.info.num_face_sets number_mesh_sets = mesh.info.num_elem_sets # Element and BC arrays for regular mesh # elementArray[i] = [v0, v1, v2, v3] stores the vertex indices (global IDs) of the # ith tetrahedral element. # bcArray[i] = [bc0, bc1, bc2, bc3] stores boundary condition codes for each face (0–3) of element i. elementArray = np.empty((number_mesh_elems, 4), dtype=np.int64) bcArray = np.zeros((number_mesh_elems, 4), dtype=np.int64) nElements = 0 for i in range(number_mesh_sets): # mesh.elem_set(i) returns a (4, num_elems) array: vertex indices of each tetrahedron in set i elemSet = np.array(mesh.elem_set(i)) num_elems = elemSet.shape[1] # verticesMap[...] remaps vertex indices to global deduplicated indices # + offset shifts local indices (e.g., when merging meshes) row_range = range(nElements, nElements + num_elems) elementArray[row_range, :] = verticesMap[elemSet + offset].T # faceBC[...] holds boundary condition codes per face (e.g., PEC, ABC, etc.) # The face permutation [3, 1, 2, 0] is used to reorder face BCs consistently bcArray[row_range, :] = faceBC[row_range][:, [3, 1, 2, 0]] nElements += num_elems sorted_indices = np.argsort(elementArray, axis=1) # Sort vertex indices within each tetrahedron to normalize ordering (for hashing, comparison, deduplication) # Also permute boundary conditions accordingly so they stay matched to the correct face elementArray = np.take_along_axis(elementArray, sorted_indices, axis=1) bcArray = np.take_along_axis(bcArray, sorted_indices, axis=1) return elementArray, bcArray def snap_to_grid(reference_node, grid_size): reference_node = np.asarray(reference_node) abs_node = np.abs(reference_node) shifted = abs_node + (grid_size / 2.0) floored = np.floor(shifted / grid_size) snapped_unsigned = floored * grid_size snapped = snapped_unsigned * np.sign(reference_node) return snapped def snap_vertices_relative_to_center(regVertices, referenceNode, regularElementSize, unit2trunc): """ Snap regVertices to a grid aligned with a coarsely-snapped referenceNode. Parameters: - regVertices: np.ndarray, shape (N, 3) The vertices to be snapped to a regular grid. - referenceNode: np.ndarray, shape (3,) The point around which snapping is aligned. - regularElementSize: float The fine grid spacing. - unit2trunc: float The coarser unit to truncate the offset from referenceNode. Returns: - snappedVertices: np.ndarray, shape (N, 3) The grid-aligned version of regVertices. """ regVertices = np.asarray(regVertices) referenceNode = np.asarray(referenceNode) # Step 1: Snap referenceNode to regular grid snapped = snap_to_grid(referenceNode, regularElementSize) # Step 2: Offset from actual reference to snapped offset = referenceNode - snapped print("Offset from referenceNode to snapped grid:", offset) # Step 3: Truncate that offset to a coarser grid offset_trunc = np.floor(offset / unit2trunc) * unit2trunc # Step 4: Center vertices relative to this truncated offset, then snap centered_vertices = regVertices - offset_trunc[:,None] # shape (N, 3) # Apply grid snapping around the center snapped_vertices = ( np.floor((np.abs(centered_vertices) + regularElementSize / 2) / regularElementSize) * regularElementSize * np.sign(centered_vertices) + offset_trunc[:,None] ).T return snapped_vertices, offset_trunc def snap_offset(vertices, regularElementSize, offset_trunc): # Step 4: Center vertices relative to this truncated offset, then snap centered_vertices = vertices - offset_trunc[:,None] # shape (N, 3) # Apply grid snapping around the center snapped_vertices = ( np.floor((np.abs(centered_vertices) + regularElementSize / 2) / regularElementSize) * regularElementSize * np.sign(centered_vertices) + offset_trunc[:,None] ).T return snapped_vertices def generate_stacked_pml_layers(base_mesh, number_of_layers, base_size, normal, corner, length, width): """ Generate multiple stacked PML layers without using geo.Vector. Parameters: base_mesh : Initial VolumeMesh (starting layer) number_of_layers : Number of stacked layers to generate base_size : Grid cube size normal : (3,) numpy array of floats corner : (3,) numpy array of floats length : Length in in-plane axis 1 width : Width in in-plane axis 2 Returns: layers : List of VolumeMesh layers layer_verts : List of (N_i, 3) numpy vertex arrays per layer """ import numpy as np layers = [] layer_verts = [] axis = int(np.argmax(np.abs(normal))) direction = np.sign(normal[axis]) current_mesh = base_mesh for n in range(number_of_layers): adjusted_corner = corner.copy() adjusted_corner[axis] += direction * base_size * n print(f"adjusted_corner: {adjusted_corner}") # Convert to geo.Vector normal_vec = Vector(*normal) corner_vec = Vector(*adjusted_corner) # Call C++ layer generator layer = base_mesh.makeGridPMLLayer(base_size, normal_vec, corner_vec, length, width) verts = np.array(layer.vertices).T verts = verts.T print(f" Vertices: min = {verts.min(axis=1)}, max = {verts.max(axis=1)}") layers.append(layer) layer_verts.append(verts) current_mesh = layer return layers, layer_verts def generate_pml_layer(mesh, base_size, reg_min, reg_max): """ Generate 6-directional PML layers using bounding box of regular region. Parameters: mesh : Initial VolumeMesh base_size : Cube size for grid elements reg_min : Tuple (x, y, z), lower corner of regular mesh reg_max : Tuple (x, y, z), upper corner of regular mesh Returns: layers: list of VolumeMesh objects (one per layer) layer_verts: list of (N_i, 3) numpy arrays of vertices (one per layer) """ max_x, max_y, max_z = reg_max min_x, min_y, min_z = reg_min Length_X = max_x - min_x Length_Y = max_y - min_y Length_Z = max_z - min_z number_of_layers = 10 layers = [] layer_verts = [] directions = [ ("+X", [1.0, 0.0, 0.0], [max_x, min_y, min_z], Length_Y, Length_Z), ("-X", [-1.0, 0.0, 0.0], [min_x, min_y, min_z], Length_Y, Length_Z), ("+Y", [0.0, 1.0, 0.0], [min_x, max_y, min_z], Length_X, Length_Z), ("-Y", [0.0, -1.0, 0.0], [min_x, min_y, min_z], Length_X, Length_Z), ("+Z", [0.0, 0.0, 1.0], [min_x, min_y, max_z], Length_X, Length_Y), ("-Z", [0.0, 0.0, -1.0], [min_x, min_y, min_z], Length_X, Length_Y), ] for label, normal, corner, length, width in directions: print(f"\n########## GENERATING PML LAYER {label} ##########") print(f"Normal: {normal}, Corner: {corner}, Length: {length}, Width: {width}") sub_layers, sub_verts = generate_stacked_pml_layers( mesh, number_of_layers, base_size, np.array(normal, dtype=float), np.array(corner, dtype=float), length, width ) layers.extend(sub_layers) layer_verts.extend(sub_verts) return layers, layer_verts def generate_cuboid_pml_layer(mesh, base_size, reg_min, reg_max, layer_thickness=10): """ Generate a full PML box (6-directional shell) with a hollow core surrounding the regular region. Parameters: mesh : Initial VolumeMesh object (with method makeGridPMLBox) base_size : Cube side length reg_min : Tuple (x, y, z) — lower corner of the regular region reg_max : Tuple (x, y, z) — upper corner of the regular region layer_thickness: Number of grid cells for PML thickness in each direction Returns: layers : list with one VolumeMesh PML box layer_verts : list with one (N, 3) numpy array of vertices """ reg_min = np.array(reg_min) reg_max = np.array(reg_max) reg_dims = reg_max - reg_min # Compute number of grid cells along each axis num_cells = np.ceil((reg_dims + 2 * layer_thickness * base_size) / base_size).astype(int) total_min = reg_min - layer_thickness * base_size total_max = reg_max + layer_thickness * base_size # Compute indices of regular (hollow) region in the grid inner_start = np.array([layer_thickness, layer_thickness, layer_thickness]) inner_end = inner_start + np.floor(reg_dims / base_size).astype(int) # Create the PML box using the C++ function via pybind11 pml_mesh = mesh.makeGridPMLBox( size=base_size, numCells=tuple(num_cells), innerStart=tuple(inner_start), innerEnd=tuple(inner_end), corner=tuple(total_min) ) # Extract vertices as (N, 3) numpy array verts = np.array(pml_mesh.vertices) # verts() → shape (3, N), transpose to (N, 3) return pml_mesh, verts def export_bc_surfaces(mesh, name, BC): """ Extracts and saves STL files for selected boundary conditions. Parameters: - mesh: VolumeMesh object (e.g., regular or irregular) - name: output name prefix - BC: BoundaryCondition enum from geo.BoundaryCondition """ print(BC) bc_map = { BC.PEC: "pec", BC.NCB: "ncbc", BC.ABC: "abc", BC.PML_Excitation: "pml_excitation", BC.PML_ABC: "pml_abc" } for bc_type, label in bc_map.items(): face_ids = mesh.face_sets_with_bc(bc_type) if face_ids: surfs = mesh.extract_surfaces(face_ids) surfs.save_stl(f"{name}_{label}_surfaces.stl") print(f"Saved: {name}_{label}_surfaces.stl") else: print(f"[Warning] No faces found for BC.{label.upper()}") def print_boundary_conditions(): print("Available Boundary Conditions:") for name in dir(geo.BoundaryCondition): if not name.startswith("__"): value = getattr(geo.BoundaryCondition, name) print(f"BC.{name} = {value}") def build_face_bc_array_custom(mesh, BC, selected_BCs, output_BC_Number): """ Builds an array of boundary condition codes for each element face with custom output values. Parameters: - mesh: the VolumeMesh object. - BC: the BoundaryCondition enum class. - selected_BCs: list of BC enum types to assign, e.g., [BC.PEC, BC.ABC] - output_BC_Number: list of integer codes corresponding to selected_BCs, e.g., [1, 2] Returns: - faceBC: (num_elems, 4) NumPy array of integer BC codes. """ assert len(selected_BCs) == len(output_BC_Number), "Length mismatch between BC list and output values" number_mesh_elems = mesh.info.num_elems number_mesh_faces = mesh.info.num_face_sets faceBC = np.zeros((number_mesh_elems, 4), dtype=np.int64) print(number_mesh_faces, "face sets found in the mesh.") # Build dictionary for fast lookup bc_to_value = dict(zip(selected_BCs, output_BC_Number)) for i in range(number_mesh_faces): bc_type = mesh.get_face_set_bc(i) if bc_type not in bc_to_value: continue # skip unselected BCs bc_val = bc_to_value[bc_type] for face_id in mesh.face_set(i): elem_idx = face_id // 4 local_face = face_id % 4 faceBC[elem_idx][local_face] = bc_val return faceBC def auto_merge_duplicate_nodes(totalNodes, tetra_data, tol=1e-8): """ Merges nodes that are within `tol` distance and updates tetra connectivity. Parameters: - totalNodes: (N, 3) array of node coordinates - tetra_data: list of tetrahedra, each as [n0, n1, n2, n3] - tol: tolerance to treat nodes as duplicates Returns: - new_totalNodes: deduplicated node coordinates - new_tetra_data: tetra connectivity updated to deduplicated nodes """ tree = cKDTree(totalNodes) close_pairs = tree.query_pairs(r=tol) # Step 1: Build a mapping from old index → representative (canonical) index parent = np.arange(len(totalNodes)) # Initially, each node is its own parent def find(i): while parent[i] != i: parent[i] = parent[parent[i]] i = parent[i] return i def union(i, j): pi, pj = find(i), find(j) if pi != pj: parent[max(pi, pj)] = min(pi, pj) for i, j in close_pairs: union(i, j) canonical = np.array([find(i) for i in range(len(totalNodes))]) # Step 2: Get unique canonical indices and build remapping _, new_indices, inverse = np.unique(totalNodes[canonical], axis=0, return_index=True, return_inverse=True) new_totalNodes = totalNodes[canonical][new_indices] # Step 3: Update tetra connectivity using canonical + inverse mapping new_tetra_data = [] for tet in tetra_data: updated = [inverse[canonical[v]] for v in tet] new_tetra_data.append(updated) return new_totalNodes, new_tetra_data def parse_node_file(lines): """ Parses the `.node` file. Each node spans two lines: a header and a line with x y z coordinates. """ num_nodes = int(lines[1].strip()) totalNodes = [] for i in range(num_nodes): coords_line = lines[3 + i * 2].strip() # Skip unit, node count, and header lines coords = list(map(float, coords_line.split())) totalNodes.append(coords) return np.array(totalNodes) def parse_tetra_file(lines): """ Parses the `.tetra` file: - First line: number of tets (optional) - Even lines: [flag node0 node1 node2 node3] - Odd lines: face BCs (ignored) Returns: - tetra_data: list of [n0, n1, n2, n3] """ num_tets = int(lines[0].strip()) tetra_data = [] for i in range(num_tets): parts = lines[1 + 2 * i].strip().split() if len(parts) != 5: continue _, n0, n1, n2, n3 = map(int, parts) tetra_data.append([n0, n1, n2, n3]) return tetra_data import numpy as np def group_tets_by_shape_material_bc( elementArray, bcArray, totalNodes, materialArray, excluded_bc_values=[], tolerance=1e-10): """ Groups tetrahedra by shape, material, and BC pattern. Tets with any face BC in excluded_bc_values are excluded (group ID = 0). Parameters: - elementArray: (N, 4) vertex indices per tet - bcArray: (N, 4) per-face BC codes - totalNodes: (M, 3) node coordinates - materialArray: (N,) material ID per tet - excluded_bc_values: list of BCs to exclude - tolerance: float, edge length quantization threshold Returns: - group_indices: (N,) array of group IDs (0 = ungrouped) - uniqueGroups: array of unique group IDs """ num_elems = elementArray.shape[0] # Step 1: Get node positions for each tet nodesArray = totalNodes[elementArray] # shape: (N, 4, 3) # Step 2: Create exclusion mask exclude_mask = np.isin(bcArray, excluded_bc_values).any(axis=1) valid_mask = ~exclude_mask valid_indices = np.where(valid_mask)[0] # Step 3: Extract shape info via edge vectors valid_nodes = nodesArray[valid_indices] edgesArray = np.empty((valid_nodes.shape[0], 3, 3)) edgesArray[:, 0, :] = valid_nodes[:, 0, :] - valid_nodes[:, 1, :] edgesArray[:, 1, :] = valid_nodes[:, 0, :] - valid_nodes[:, 2, :] edgesArray[:, 2, :] = valid_nodes[:, 0, :] - valid_nodes[:, 3, :] rounded_edges = np.round(edgesArray / tolerance) * tolerance shape_key = rounded_edges.reshape(valid_nodes.shape[0], -1) # Step 4: Combine with material and BC valid_materials = materialArray[valid_indices].reshape(-1, 1) valid_bcs = bcArray[valid_indices] full_key = np.hstack([shape_key, valid_materials, valid_bcs]) # Step 5: Group by unique key _, group_ids = np.unique(full_key, axis=0, return_inverse=True) group_ids = group_ids + 1 # reserve group ID 0 for excluded # Step 6: Final group array group_indices = np.zeros(num_elems, dtype=np.int64) group_indices[valid_indices] = group_ids # Step 7: Recompress uniqueGroups, compressed = np.unique(group_indices, return_inverse=True) group_indices = compressed # === Diagnostics === print(f"🔢 Total tets: {num_elems}") print(f"🚫 Excluded by BC list: {exclude_mask.sum()}") print(f"✅ Grouped tets: {len(valid_indices)}") print(f"🔁 Unique groups (excluding group 0): {len(uniqueGroups) - 1}") return group_indices, uniqueGroups import numpy as np def group_tets_by_box_and_conductivity( elementArray, bcArray, totalNodes, materialArray, region_corner, region_dims, pml_bounds, excluded_bc_values=[], sigma_base=0.01, m=3, tolerance=1e-10, min_group_size=1, ): num_elems = elementArray.shape[0] nodesArray = totalNodes[elementArray] centroids = nodesArray.mean(axis=1) core_lower = np.asarray(region_corner) core_upper = core_lower + np.asarray(region_dims) inside_mask = np.all((centroids >= core_lower) & (centroids <= core_upper), axis=1) outside_mask = ~inside_mask group_indices = np.zeros(num_elems, dtype=np.int64) inside_gid_set = set() outside_gid_set = set() # ======================== # Inside grouping # ======================== inside_indices = np.where(inside_mask)[0] inside_nodes = nodesArray[inside_indices] inside_bcs = bcArray[inside_indices] inside_materials = materialArray[inside_indices] exclude_mask = np.isin(inside_bcs, excluded_bc_values).any(axis=1) valid_inside_mask = ~exclude_mask valid_inside_indices = inside_indices[valid_inside_mask] valid_nodes = nodesArray[valid_inside_indices] edges = np.empty((valid_nodes.shape[0], 3, 3)) edges[:, 0, :] = valid_nodes[:, 0, :] - valid_nodes[:, 1, :] edges[:, 1, :] = valid_nodes[:, 0, :] - valid_nodes[:, 2, :] edges[:, 2, :] = valid_nodes[:, 0, :] - valid_nodes[:, 3, :] shape_key = np.round(edges / tolerance) * tolerance shape_key = shape_key.reshape(valid_nodes.shape[0], -1) valid_materials = materialArray[valid_inside_indices].reshape(-1, 1) valid_bcs = bcArray[valid_inside_indices] key_inside = np.hstack([shape_key, valid_materials, valid_bcs]) _, group_ids_inside = np.unique(key_inside, axis=0, return_inverse=True) group_ids_inside = group_ids_inside + 1 group_indices[valid_inside_indices] = group_ids_inside next_group_id = group_ids_inside.max() + 1 inside_gid_set.update(group_ids_inside) # ======================== # Outside grouping # ======================== outside_indices = np.where(outside_mask)[0] if outside_indices.size > 0: outside_nodes = nodesArray[outside_indices] outside_materials = materialArray[outside_indices] outside_bcs = bcArray[outside_indices] edges_out = np.empty((outside_nodes.shape[0], 3, 3)) edges_out[:, 0, :] = outside_nodes[:, 0, :] - outside_nodes[:, 1, :] edges_out[:, 1, :] = outside_nodes[:, 0, :] - outside_nodes[:, 2, :] edges_out[:, 2, :] = outside_nodes[:, 0, :] - outside_nodes[:, 3, :] shape_key_out = np.round(edges_out / tolerance) * tolerance shape_key_out = shape_key_out.reshape(outside_nodes.shape[0], -1) rc = centroids[outside_indices] dx = np.maximum.reduce([region_corner[0] - rc[:, 0], np.zeros_like(rc[:, 0]), rc[:, 0] - (region_corner[0] + region_dims[0])]) dy = np.maximum.reduce([region_corner[1] - rc[:, 1], np.zeros_like(rc[:, 1]), rc[:, 1] - (region_corner[1] + region_dims[1])]) dz = np.maximum.reduce([region_corner[2] - rc[:, 2], np.zeros_like(rc[:, 2]), rc[:, 2] - (region_corner[2] + region_dims[2])]) s1 = sigma_base * (dx / pml_bounds[0])**m s2 = sigma_base * (dy / pml_bounds[1])**m s3 = sigma_base * (dz / pml_bounds[2])**m sigmas = np.stack([np.round(s1, 6), np.round(s2, 6), np.round(s3, 6)], axis=1) materials_out = outside_materials.reshape(-1, 1) key_outside = np.hstack([shape_key_out, materials_out, outside_bcs, sigmas]) _, group_ids_outside = np.unique(key_outside, axis=0, return_inverse=True) group_ids_outside = group_ids_outside + next_group_id group_indices[outside_indices] = group_ids_outside outside_gid_set.update(group_ids_outside) # ======================== # Filter small groups # ======================== final = group_indices.copy() surviving_inside_gid_set = set() surviving_outside_gid_set = set() for gid in np.unique(group_indices): if gid == 0: continue count = np.sum(group_indices == gid) if count < min_group_size: final[group_indices == gid] = 0 else: if gid in inside_gid_set: surviving_inside_gid_set.add(gid) elif gid in outside_gid_set: surviving_outside_gid_set.add(gid) # ======================== # Remap surviving group IDs to be contiguous # ======================== uniqueGroups = np.unique(final) remap = {old_gid: new_gid for new_gid, old_gid in enumerate(uniqueGroups)} group_indices_remapped = np.array([remap[gid] for gid in final], dtype=np.int64) # Compute remapped surviving GIDs remapped_inside_ids = sorted([remap[gid] for gid in surviving_inside_gid_set if gid in remap]) remapped_outside_ids = sorted([remap[gid] for gid in surviving_outside_gid_set if gid in remap]) # ======================== # Diagnostics # ======================== def print_gid_range(name, gid_list): if gid_list: print(f"🔢 {name} group IDs: min={min(gid_list)}, max={max(gid_list)}") else: print(f"🔢 {name} group IDs: No valid groups") print(f"🔢 Total tets: {num_elems}") print(f"📦 Inside box: {len(inside_indices)} (grouped by shape+material+BC)") print(f"🌐 Outside box: {len(outside_indices)} (grouped by +conductivity)") print(f"🚫 Groups filtered out (<{min_group_size} tets)") print(f"🔁 Unique groups (excluding group 0): {len(uniqueGroups) - (1 if 0 in uniqueGroups else 0)}") print_gid_range("Inside", remapped_inside_ids) print_gid_range("Outside", remapped_outside_ids) return group_indices_remapped def group_tets_by_box_and_pml( elementArray, bcArray, totalNodes, materialArray, region_corner, region_dims, excluded_bc_values=[], tolerance=1e-10, min_group_size=1, ): num_elems = elementArray.shape[0] nodesArray = totalNodes[elementArray] centroids = nodesArray.mean(axis=1) core_lower = np.asarray(region_corner) core_upper = core_lower + np.asarray(region_dims) inside_mask = np.all((centroids >= core_lower) & (centroids <= core_upper), axis=1) outside_mask = ~inside_mask group_indices = np.zeros(num_elems, dtype=np.int64) inside_gid_set = set() # ======================== # Inside grouping only # ======================== inside_indices = np.where(inside_mask)[0] inside_nodes = nodesArray[inside_indices] inside_bcs = bcArray[inside_indices] inside_materials = materialArray[inside_indices] exclude_mask = np.isin(inside_bcs, excluded_bc_values).any(axis=1) valid_inside_mask = ~exclude_mask valid_inside_indices = inside_indices[valid_inside_mask] valid_nodes = nodesArray[valid_inside_indices] edges = np.empty((valid_nodes.shape[0], 3, 3)) edges[:, 0, :] = valid_nodes[:, 0, :] - valid_nodes[:, 1, :] edges[:, 1, :] = valid_nodes[:, 0, :] - valid_nodes[:, 2, :] edges[:, 2, :] = valid_nodes[:, 0, :] - valid_nodes[:, 3, :] shape_key = np.round(edges / tolerance) * tolerance shape_key = shape_key.reshape(valid_nodes.shape[0], -1) valid_materials = materialArray[valid_inside_indices].reshape(-1, 1) valid_bcs = bcArray[valid_inside_indices] key_inside = np.hstack([shape_key, valid_materials, valid_bcs]) _, group_ids_inside = np.unique(key_inside, axis=0, return_inverse=True) group_ids_inside = group_ids_inside + 1 group_indices[valid_inside_indices] = group_ids_inside inside_gid_set.update(group_ids_inside) # ======================== # Filter small groups # ======================== final = group_indices.copy() surviving_inside_gid_set = set() for gid in np.unique(group_indices): if gid == 0: continue count = np.sum(group_indices == gid) if count < min_group_size: final[group_indices == gid] = 0 else: surviving_inside_gid_set.add(gid) # ======================== # Remap surviving group IDs to be contiguous # ======================== uniqueGroups = np.unique(final) remap = {old_gid: new_gid for new_gid, old_gid in enumerate(uniqueGroups)} group_indices_remapped = np.array([remap[gid] for gid in final], dtype=np.int64) # Compute remapped surviving GIDs remapped_inside_ids = sorted([remap[gid] for gid in surviving_inside_gid_set if gid in remap]) # ======================== # Diagnostics # ======================== def print_gid_range(name, gid_list): if gid_list: print(f"🔢 {name} group IDs: min={min(gid_list)}, max={max(gid_list)}") else: print(f"🔢 {name} group IDs: No valid groups") print(f"🔢 Total tets: {num_elems}") print(f"📦 Inside box: {len(inside_indices)} (grouped by shape+material+BC)") print(f"🌐 Outside box: {len(np.where(outside_mask)[0])} (set to group 0)") print(f"🚫 Groups filtered out (<{min_group_size} tets)") print(f"🔁 Unique groups (excluding group 0): {len(uniqueGroups) - (1 if 0 in uniqueGroups else 0)}") print_gid_range("Inside", remapped_inside_ids) return group_indices_remapped def group_tets_by_box_and_conductivity_zero_bc_only( elementArray, bcArray, totalNodes, materialArray, region_corner, region_dims, pml_bounds, sigma_base=0.01, m=3, tolerance=1e-10, min_group_size=1, ): """ Group tets by shape+material (+BC==all-zero) inside the box, and by shape+material+BC==all-zero+PML conductivity outside the box. Only tets with per-face BCs == [0,0,0,0] are eligible for any group. Others are forced to group 0. """ num_elems = elementArray.shape[0] nodesArray = totalNodes[elementArray] centroids = nodesArray.mean(axis=1) # ---- core box masks ---- core_lower = np.asarray(region_corner) core_upper = core_lower + np.asarray(region_dims) inside_mask = np.all((centroids >= core_lower) & (centroids <= core_upper), axis=1) outside_mask = ~inside_mask # ---- NEW: require all-zero BCs ---- bc_all_zero_mask = np.all(bcArray == 0, axis=1) # result container group_indices = np.zeros(num_elems, dtype=np.int64) inside_gid_set = set() outside_gid_set = set() # ======================== # Inside grouping (only BC==0 rows) # ======================== inside_indices = np.where(inside_mask & bc_all_zero_mask)[0] if inside_indices.size > 0: valid_nodes = nodesArray[inside_indices] # Geometry shape key (three edges from node 0) edges = np.empty((valid_nodes.shape[0], 3, 3)) edges[:, 0, :] = valid_nodes[:, 0, :] - valid_nodes[:, 1, :] edges[:, 1, :] = valid_nodes[:, 0, :] - valid_nodes[:, 2, :] edges[:, 2, :] = valid_nodes[:, 0, :] - valid_nodes[:, 3, :] shape_key = np.round(edges / tolerance) * tolerance shape_key = shape_key.reshape(valid_nodes.shape[0], -1) valid_materials = materialArray[inside_indices].reshape(-1, 1) # BCs are known to be zeros; you can omit them in the key or keep them for clarity. key_inside = np.hstack([shape_key, valid_materials]) # BCs omitted because all-zero _, group_ids_inside = np.unique(key_inside, axis=0, return_inverse=True) group_ids_inside = group_ids_inside + 1 # start at 1 group_indices[inside_indices] = group_ids_inside next_group_id = group_ids_inside.max() + 1 inside_gid_set.update(group_ids_inside) else: next_group_id = 1 # ======================== # Outside grouping (only BC==0 rows) # ======================== outside_indices = np.where(outside_mask & bc_all_zero_mask)[0] if outside_indices.size > 0: outside_nodes = nodesArray[outside_indices] outside_materials = materialArray[outside_indices] # Geometry key (same as inside) edges_out = np.empty((outside_nodes.shape[0], 3, 3)) edges_out[:, 0, :] = outside_nodes[:, 0, :] - outside_nodes[:, 1, :] edges_out[:, 1, :] = outside_nodes[:, 0, :] - outside_nodes[:, 2, :] edges_out[:, 2, :] = outside_nodes[:, 0, :] - outside_nodes[:, 3, :] shape_key_out = np.round(edges_out / tolerance) * tolerance shape_key_out = shape_key_out.reshape(outside_nodes.shape[0], -1) # Simple “distance outside box” PML conductivity profile (component-wise) rc = centroids[outside_indices] dx = np.maximum.reduce([ region_corner[0] - rc[:, 0], np.zeros_like(rc[:, 0]), rc[:, 0] - (region_corner[0] + region_dims[0]), ]) dy = np.maximum.reduce([ region_corner[1] - rc[:, 1], np.zeros_like(rc[:, 1]), rc[:, 1] - (region_corner[1] + region_dims[1]), ]) dz = np.maximum.reduce([ region_corner[2] - rc[:, 2], np.zeros_like(rc[:, 2]), rc[:, 2] - (region_corner[2] + region_dims[2]), ]) s1 = sigma_base * (dx / pml_bounds[0])**m s2 = sigma_base * (dy / pml_bounds[1])**m s3 = sigma_base * (dz / pml_bounds[2])**m sigmas = np.stack([np.round(s1, 6), np.round(s2, 6), np.round(s3, 6)], axis=1) materials_out = outside_materials.reshape(-1, 1) # BCs are all-zero here too; key uses geometry + material + sigmas key_outside = np.hstack([shape_key_out, materials_out, sigmas]) _, group_ids_outside = np.unique(key_outside, axis=0, return_inverse=True) group_ids_outside = group_ids_outside + next_group_id group_indices[outside_indices] = group_ids_outside outside_gid_set.update(group_ids_outside) # ======================== # Filter small groups # ======================== final = group_indices.copy() surviving_inside_gid_set = set() surviving_outside_gid_set = set() for gid in np.unique(group_indices): if gid == 0: continue count = np.sum(group_indices == gid) if count < min_group_size: final[group_indices == gid] = 0 else: if gid in inside_gid_set: surviving_inside_gid_set.add(gid) elif gid in outside_gid_set: surviving_outside_gid_set.add(gid) # ======================== # Remap surviving group IDs to be contiguous (0..K) # ======================== uniqueGroups = np.unique(final) remap = {old_gid: new_gid for new_gid, old_gid in enumerate(uniqueGroups)} group_indices_remapped = np.array([remap[gid] for gid in final], dtype=np.int64) # Compute remapped surviving GIDs remapped_inside_ids = sorted([remap[gid] for gid in surviving_inside_gid_set if gid in remap]) remapped_outside_ids = sorted([remap[gid] for gid in surviving_outside_gid_set if gid in remap]) # ======================== # Diagnostics # ======================== def print_gid_range(name, gid_list): if gid_list: print(f"🔢 {name} group IDs: min={min(gid_list)}, max={max(gid_list)}") else: print(f"🔢 {name} group IDs: No valid groups") print(f"🔢 Total tets: {num_elems}") print(f"✅ Eligible (BC==[0,0,0,0]): {int(bc_all_zero_mask.sum())}") print(f"🚫 Ineligible (any non-zero BC): {int((~bc_all_zero_mask).sum())}") print(f"📦 Inside box (eligible only): {len(inside_indices)}") print(f"🌐 Outside box (eligible only): {len(outside_indices)}") print(f"🚫 Groups filtered out (<{min_group_size} tets)") print(f"🔁 Unique groups (excluding group 0): {len(uniqueGroups) - (1 if 0 in uniqueGroups else 0)}") print_gid_range("Inside", remapped_inside_ids) print_gid_range("Outside", remapped_outside_ids) return group_indices_remapped