import sys sys.path.append('.') import Geometry as geo import numpy as np import pandas as pd from numpy import linalg as la ####################################################### ####################################################### ## ## ## PART TO MODIFY ## ## ## ####################################################### ####################################################### directory = "/home/user/Projects/problems/testCylinder/" # Absolute path to the model files fileName = "test_cylinder" # 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.08 # Mesh size in the unit used unit2trunc = 4 # Number of significant digits in the regular mesh PEC = 1 # PEC in the .bc file ABC = 2 # PW in the .bc file ####################################################### ####################################################### ## ## ## DO NOT CHANGE FROM HERE ## ## ## ####################################################### ####################################################### ######################################## ## AUXILIARY FUNCTIONS ## ######################################## def floor(number, sig_figs): return np.trunc(number * (10 ** sig_figs)) / (10 ** sig_figs) 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)) ######################################## ## MAIN FUNCTION ## ######################################## def main(): name = f"{directory}/{fileName}" regularElementSize = truncate_to_significant_figures(mesh_size, 4) print(f"Using {directory}/{fileName}.g") mesh = geo.VolumeMesh.from_exodus(f"{directory}/{fileName}.g") config = mesh.apply_config(f"{directory}/{fileName}.json") BC = geo.BoundaryCondition regular = mesh.make_grid_mesh(regularElementSize, mesh.face_sets_with_bc(BC.PEC)) irregular = mesh.make_hybrid_mesh(regularElementSize, mesh.face_sets_with_bc(BC.PEC), max_volume=pow(regularElementSize,3)*np.sqrt(2)/12) surfs = irregular.extract_surfaces(irregular.face_sets_with_bc(BC.PEC)) surfs.save_stl(f"{name}_pec_surfaces.stl") surfs = regular.extract_surfaces(regular.face_sets_with_bc(BC.NCB)) surfs.save_stl(f"{name}_ncbc_surfaces.stl") surfs = regular.extract_surfaces(regular.face_sets_with_bc(BC.ABC)) surfs.save_stl(f"{name}_acb_surfaces.stl") regularVertices = regular.vertices.shape[1] regVertices = np.array(regular.vertices) referenceNode = regular.vertices[:,0] center = referenceNode - np.floor((np.absolute(referenceNode) + (regularElementSize / 2.0)) / regularElementSize) * regularElementSize * np.sign(referenceNode) center = floor(center, unit2trunc) regVertices = np.floor((np.absolute(regVertices - center.reshape((3,1))) + (regularElementSize / 2.0)) / regularElementSize) * regularElementSize * np.sign(regVertices) + center.reshape((3,1)) regVertices = regVertices.T totalNodes = np.vstack((np.array(regular.vertices).T, np.array(irregular.vertices).T)) totalNodes, indices, verticesMap, _ = np.unique(totalNodes, axis=0, return_index=True, return_inverse=True, return_counts=True) np.savetxt(f'{name}.node', totalNodes, newline="\n3 0 0\n" , fmt='%f', header=f"{unit}\n{totalNodes.shape[0]}", comments='') totalNodes = np.vstack((regVertices, np.array(irregular.vertices).T))[indices, :] ######################################## ## REGULAR ELEMENTS ## ######################################## dimRegular = (regular.info.num_elems, 4) regularFacesBC = np.zeros(dimRegular) NCBC = 1000 for i in range(regular.info.num_face_sets): bcValue = 0 if (regular.get_face_set_bc(i) == BC.PEC): bcValue = PEC elif (regular.get_face_set_bc(i) == BC.ABC): bcValue = ABC elif (regular.get_face_set_bc(i) == BC.NCB): bcValue = NCBC else: bcValue = 0 for j in range(len(regular.face_set(i))): regularFacesBC[regular.face_set(i)[j] // 4][regular.face_set(i)[j] % 4] = int(bcValue) elementArray = np.empty((regular.info.num_elems, 4), dtype=np.int64) bcArray = np.zeros((regular.info.num_elems, 4), dtype=np.int64) nElements = 0 for i in range(regular.info.num_elem_sets): elemSet = np.array(regular.elem_set(i)) num_elements = elemSet.shape[1] row_range = range(nElements, nElements + num_elements) elementArray[row_range, :] = verticesMap[elemSet].T bcArray[row_range, :] = regularFacesBC[row_range][:, [3, 1, 2, 0]] nElements += num_elements sorted_indices = np.argsort(elementArray, axis=1) elementArray = np.take_along_axis(elementArray, sorted_indices, axis=1) bcArray = np.take_along_axis(bcArray, sorted_indices, axis=1) ################################# ## GROUP REGULARS ## ################################# nodesArray = np.empty((regular.info.num_elems, 4, 3)) nodesArray[:,0,:] = [totalNodes[i][:3] for i in elementArray[:, 0]] nodesArray[:,1,:] = [totalNodes[i][:3] for i in elementArray[:, 1]] nodesArray[:,2,:] = [totalNodes[i][:3] for i in elementArray[:, 2]] nodesArray[:,3,:] = [totalNodes[i][:3] for i in elementArray[:, 3]] edgesArray = np.empty((regular.info.num_elems, 3, 3)) edgesArray[:,0,:] = nodesArray[:,0,:] - nodesArray[:,1,:] edgesArray[:,1,:] = nodesArray[:,0,:] - nodesArray[:,2,:] edgesArray[:,2,:] = nodesArray[:,0,:] - nodesArray[:,3,:] tolerance = 1e-10 _, indices = np.unique(np.round(edgesArray / tolerance) * tolerance, axis=0, return_inverse=True) indices = indices + np.ones_like(indices) indices[np.any(bcArray != 0, axis=1)] = 0 uniqueValues, indices = np.unique(indices, axis=0, return_inverse=True) indices = np.hstack((indices, np.zeros((irregular.info.num_elems,),dtype=np.int64))) ######################################## ## IRREGULAR ELEMENTS ## ######################################## dimIrregular = (irregular.info.num_elems, 4) irregularFacesBC = np.zeros(dimIrregular) for i in range(irregular.info.num_face_sets): bcValue = 0 if (irregular.get_face_set_bc(i) == BC.PEC): bcValue = PEC elif (irregular.get_face_set_bc(i) == BC.ABC): bcValue = ABC elif (irregular.get_face_set_bc(i) == BC.NCB): bcValue = NCBC else: bcValue = 0 for j in range(len(irregular.face_set(i))): irregularFacesBC[irregular.face_set(i)[j] // 4][irregular.face_set(i)[j] % 4] = int(bcValue) irregularElementArray = np.empty((irregular.info.num_elems, 4), dtype=np.int64) irregularBcArray = np.zeros((irregular.info.num_elems, 4), dtype=np.int64) nElements = 0 for i in range(irregular.info.num_elem_sets): elemSet = np.array(irregular.elem_set(i)) num_elements = elemSet.shape[1] row_range = range(nElements, nElements + num_elements) irregularElementArray[row_range, :] = verticesMap[elemSet + regularVertices].T irregularBcArray[row_range, :] = irregularFacesBC[row_range][:, [3, 1, 2, 0]] nElements += num_elements ######################################## ## JOIN ## ######################################## elementArray = np.vstack((elementArray, irregularElementArray)) bcArray = np.vstack((bcArray, irregularBcArray)) np.savetxt(f'{name}.regular', indices, fmt='%d', header=str(uniqueValues.shape[0]), comments='') bcArray[bcArray == NCBC] = 0 print(f"The number of elements is {bcArray.shape[0]}") ######################################## ## PRINT ## ######################################## c_with_extra = np.ones((2 * (regular.info.num_elems + irregular.info.num_elems), 4 + 1), dtype=np.int64) * (-1) c_with_extra[::2] = np.hstack((np.ones(((regular.info.num_elems + irregular.info.num_elems),1), dtype=np.int64) * 0, elementArray)) c_with_extra[1::2, 1:] = bcArray df = pd.DataFrame(c_with_extra, dtype=np.int64) mask = df != -1 df = df.where(mask, other=np.nan) new_row = pd.Series([(regular.info.num_elems + irregular.info.num_elems)] + [np.nan] * (df.shape[1] - 1), index=df.columns) # Insert the new row at the top df = pd.concat([new_row.to_frame().T, df], ignore_index=True) df = df.astype('Int64') df.to_csv(f'{name}.tetra', sep=' ', index=False, header=False) if __name__ == '__main__': main()