diff --git a/modules/tools/CMakeLists.txt b/modules/tools/CMakeLists.txt index 9d754e467..6d5bfa280 100644 --- a/modules/tools/CMakeLists.txt +++ b/modules/tools/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(BRAILS) add_subdirectory(ShakerMaker) +add_subdirectory(DRM) diff --git a/modules/tools/DRM/CMakeLists.txt b/modules/tools/DRM/CMakeLists.txt new file mode 100644 index 000000000..502bdd41f --- /dev/null +++ b/modules/tools/DRM/CMakeLists.txt @@ -0,0 +1,3 @@ +add_subdirectory(MeshGenerator) +add_subdirectory(StructresFiles) +simcenter_add_python_script(SCRIPT modelCreator.py) diff --git a/modules/tools/DRM/MeshGenerator/CMakeLists.txt b/modules/tools/DRM/MeshGenerator/CMakeLists.txt new file mode 100644 index 000000000..a19033786 --- /dev/null +++ b/modules/tools/DRM/MeshGenerator/CMakeLists.txt @@ -0,0 +1,5 @@ +add_subdirectory(models) +simcenter_add_python_script(SCRIPT Infowriter.py) +simcenter_add_python_script(SCRIPT MeshGenrator.py) +simcenter_add_python_script(SCRIPT Structure.py) + diff --git a/modules/tools/DRM/MeshGenerator/Infowriter.py b/modules/tools/DRM/MeshGenerator/Infowriter.py new file mode 100644 index 000000000..43c226d25 --- /dev/null +++ b/modules/tools/DRM/MeshGenerator/Infowriter.py @@ -0,0 +1,153 @@ + + +def infowriter(info,meshdir): + """ + Write the information dictionary to a tcl file. + + Parameters + ---------- + info : dict + Information dictionary. + meshdir : str + Directory where the tcl file will be written + + Returns + ------- + None + + """ + # ============================================================================ + # Create a tcl file + # ============================================================================ + f = open(f"{meshdir}/Modelinfo.tcl", "w") + f.write(f"wipe\n") + f.write(f"# ============================================================================\n") + f.write(f"# Cores Information\n") + f.write(f"# ============================================================================\n") + f.write(f"set regcores {info['regcores']}\n") + f.write(f"set pmlcores {info['pmlcores']}\n") + f.write(f"set drmcores {info['drmcores']}\n") + f.write(f"set structurecores {info['structurecores']}\n") + f.write(f"set AnalysisType {info['AnalysisType']}\n") + + + f.write(f"# ============================================================================\n") + f.write(f"# Model Information\n") + f.write(f"# ============================================================================\n") + f.write(f"set StructureType \"{info['StructureType']}\"\n") + f.write(f"set NStory {info['NStory']}\n") + f.write(f"set NBay {info['NBay']}\n") + f.write(f"set NBayZ {info['NBayZ']}\n") + f.write(f"set StartNodeX {info['StartNodeX']}\n") + f.write(f"set StartNodeY {info['StartNodeY']}\n") + f.write(f"set StartNodeZ {info['StartNodeZ']}\n") + f.write(f"set LCol {info['LCol']}\n") + f.write(f"set LBeam {info['LBeam']}\n") + f.write(f"set LGird {info['LGird']}\n") + f.write(f"set SectionType {info['SectionType']}\n") + f.write(f"set HaveStructure \"{info['HaveStructure']}\"\n") + f.write(f"# ============================================================================\n") + f.write(f"# Soil Information\n") + f.write(f"# ============================================================================\n") + f.write(f"set dx {info['dx']}\n") + f.write(f"set dy {info['dy']}\n") + f.write(f"set dz {info['dz']}\n") + f.write(f"set llx {info['llx']}\n") + f.write(f"set lly {info['lly']}\n") + f.write(f"set llz {info['llz']}\n") + f.write(f"set drmthicknessx {info['drmthicknessx']}\n") + f.write(f"set drmthicknessy {info['drmthicknessy']}\n") + f.write(f"set drmthicknessz {info['drmthicknessz']}\n") + f.write(f"set numdrmlayers {info['numdrmlayers']}\n") + f.write(f"set lx {info['lx']}\n") + f.write(f"set ly {info['ly']}\n") + f.write(f"set lz {info['lz']}\n") + f.write(f"set nx {info['nx']}\n") + f.write(f"set ny {info['ny']}\n") + f.write(f"set nz {info['nz']}\n") + f.write(f"# ============================================================================\n") + f.write(f"# PML information\n") + f.write(f"# ============================================================================\n") + f.write(f"set AbsorbingElements \"{info['AbsorbingElements']}\"\n") + f.write(f"set numpmllayers {info['numpmllayers']}\n") + f.write(f"set pmlthicknessx {info['pmlthicknessx']}\n") + f.write(f"set pmlthicknessy {info['pmlthicknessy']}\n") + f.write(f"set pmlthicknessz {info['pmlthicknessz']}\n") + f.write(f"set pmltotalthickness {info['pmltotalthickness']}\n") + f.write(f"set HaveAbsorbingElements \"{info['HaveAbsorbingElements']}\"\n") + f.write(f"set Absorbing_rayleigh_alpha {info['Absorbing_rayleigh_alpha']}\n") + f.write(f"set Absorbing_rayleigh_beta {info['Absorbing_rayleigh_beta']}\n") + f.write(f"# ============================================================================\n") + f.write(f"# General information\n") + f.write(f"# ============================================================================\n") + f.write(f"set meshdir \"{info['meshdir']}\"\n") + f.write(f"set outputdir \"{info['outputdir']}\"\n") + f.write(f"# ============================================================================\n") + f.write(f"# Embedding foundation\n") + f.write(f"# ============================================================================\n") + f.write(f"set EmbeddingFoundation \"{info['EmbeddingFoundation']}\"\n") + info2 = info['EmbeddedFoundation'] + f.write(f"set EmbeddedFoundation [dict create xmax {info2['xmax']} xmin {info2['xmin']} ymax {info2['ymax']} ymin {info2['ymin']} zmax {info2['zmax']} zmin {info2['zmin']}]\n") + f.write(f"# ============================================================================\n") + f.write(f"# Fondation information\n") + f.write(f"# ============================================================================\n") + f.write(f"set HaveFoundation \"{info['HaveFoundation']}\"\n") + f.write("set foundationBlocks {}\n") + for i,block in enumerate(info['foundationBlocks']): + f.write(f"lappend foundationBlocks [dict create matTag {block['matTag']} xmax {block['xmax']} xmin {block['xmin']} ymax {block['ymax']} ymin {block['ymin']} zmax {block['zmax']} zmin {block['zmin']} Xmeshsize {block['Xmeshsize']} Ymeshsize {block['Ymeshsize']} Zmeshsize {block['Zmeshsize']}]\n") + f.write(f"# ============================================================================\n") + f.write(f"# Piles information\n") + f.write(f"# ============================================================================\n") + f.write(f"set HavePiles \"{info['HavePiles']}\"\n") + f.write("set pilelist {}\n") + for i,pile in enumerate(info['pilelist']): + f.write(f"lappend pilelist [dict create xtop {pile['xtop']} ytop {pile['ytop']} ztop {pile['ztop']} xbottom {pile['xbottom']} ybottom {pile['ybottom']} zbottom {pile['zbottom']} numberofElements {pile['numberofElements']}]\n") + f.write(f"# ============================================================================\n") + f.write(f"# cells and nodes information\n") + f.write(f"# ============================================================================\n") + f.write(f"set soilfoundation_num_cells {info['soilfoundation_num_cells']}\n") + f.write(f"set soilfoundation_num_nodes {info['soilfoundation_num_points']}\n") + f.write(f"# ============================================================================\n") + f.write(f"# DRM information\n") + f.write(f"# ============================================================================\n") + f.write(f"set DRMFile \"{info['DRMFile']}\"\n") + f.write(f"set DRM_Provider_Software \"{info['DRMinformation']['DRM_Provider_Software']}\"\n") + f.write(f"set DRM_factor \"{info['DRMinformation']['factor']}\"\n") + f.write(f"set crd_scale \"{info['DRMinformation']['crd_scale']}\"\n") + f.write(f"set distance_tolerance \"{info['DRMinformation']['distance_tolerance']}\"\n") + f.write(f"set do_coordinate_transformation \"{info['DRMinformation']['do_coordinate_transformation']}\"\n") + f.write(f"set T00 \"{info['DRMinformation']['T00']}\"\n") + f.write(f"set T01 \"{info['DRMinformation']['T01']}\"\n") + f.write(f"set T02 \"{info['DRMinformation']['T02']}\"\n") + f.write(f"set T10 \"{info['DRMinformation']['T10']}\"\n") + f.write(f"set T11 \"{info['DRMinformation']['T11']}\"\n") + f.write(f"set T12 \"{info['DRMinformation']['T12']}\"\n") + f.write(f"set T20 \"{info['DRMinformation']['T20']}\"\n") + f.write(f"set T21 \"{info['DRMinformation']['T21']}\"\n") + f.write(f"set T22 \"{info['DRMinformation']['T22']}\"\n") + f.write(f"set originX \"{info['DRMinformation']['originX']}\"\n") + f.write(f"set originY \"{info['DRMinformation']['originY']}\"\n") + f.write(f"set originZ \"{info['DRMinformation']['originZ']}\"\n") + f.write(f"set DRM_Location \"{info['DRMinformation']['DRM_Location'].lower()}\"\n") + f.write(f"# ============================================================================\n") + f.write(f"# Analysis information\n") + f.write(f"# ============================================================================\n") + f.write(f"set Analysis_dt {info['AnalysisInfo']['dt']}\n") + f.write(f"set Analysis_duration {info['AnalysisInfo']['t_final']}\n") + f.write(f"set Analysis_record_dt {info['AnalysisInfo']['recording_dt']}\n") + f.write(f"# ============================================================================\n") + f.write(f"# material information\n") + f.write(f"# ============================================================================\n") + f.write(f"set Vs {info['MaterialInfo']['Vs']}\n") + f.write(f"set rho {info['MaterialInfo']['rho']}\n") + f.write(f"set nu {info['MaterialInfo']['nu']}\n") + f.close() + + + return None + + + + + + diff --git a/modules/tools/DRM/MeshGenerator/MeshGenrator.py b/modules/tools/DRM/MeshGenerator/MeshGenrator.py new file mode 100644 index 000000000..5bf7e75ca --- /dev/null +++ b/modules/tools/DRM/MeshGenerator/MeshGenrator.py @@ -0,0 +1,736 @@ +# %% +import pyvista as pv +import numpy as np +import ctypes +import os +import sys +import h5py +from scipy.spatial import KDTree + +def partition_with_pilebox(messh,numcores, pileboxinfo,tol=1e-6): + if pileboxinfo["Havepilebox"] == False : + print("No pilebox is defined") + partition(messh,numcores) + return + + # seperate the core of the pilebox + mesh = messh.copy() + eps = 1e-6 + # check if the xmin,xmax, ymin ymax ,zmin zmax is exist in the infokeys + if "xmin" in pileboxinfo.keys(): + pilebox_xmin = pileboxinfo["xmin"] + if "xmax" in pileboxinfo.keys(): + pilebox_xmax = pileboxinfo["xmax"] + if "ymin" in pileboxinfo.keys(): + pilebox_ymin = pileboxinfo["ymin"] + if "ymax" in pileboxinfo.keys(): + pilebox_ymax = pileboxinfo["ymax"] + if "zmin" in pileboxinfo.keys(): + pilebox_zmin = pileboxinfo["zmin"] + if "zmax" in pileboxinfo.keys(): + pilebox_zmax = pileboxinfo["zmax"] + if "center" in pileboxinfo.keys(): + pilebox_xmin = pileboxinfo["center"][0] - pileboxinfo["lx"]/2 + eps + pilebox_xmax = pileboxinfo["center"][0] + pileboxinfo["lx"]/2 - eps + pilebox_ymin = pileboxinfo["center"][1] - pileboxinfo["ly"]/2 + eps + pilebox_ymax = pileboxinfo["center"][1] + pileboxinfo["ly"]/2 - eps + pilebox_zmin = pileboxinfo["center"][2] - pileboxinfo["depth"] + eps + pilebox_zmax = pileboxinfo["center"][2] + tol + eps + + + # find the cells that are inside the pilebox + cube = pv.Cube(bounds=[pilebox_xmin,pilebox_xmax,pilebox_ymin,pilebox_ymax,pilebox_zmin,pilebox_zmax]) + # extract the cells that are outside the pilebox + indices = mesh.find_cells_within_bounds(cube.bounds) + + # crete newindices for outside the pilebox + newindices = np.ones(mesh.n_cells,dtype=bool) + newindices[indices] = False + newindices = np.where(newindices)[0] + + + # extract the cells that are inside the pilebox + mesh2 = mesh.extract_cells(newindices) + # partition the mesh + + if numcores > 2: + partition(mesh2,numcores-1) + if numcores == 2: + mesh2.cell_data['partitioned'] = np.zeros(mesh2.n_cells,dtype=np.int32) + + mesh.cell_data['partitioned'] = np.zeros(mesh.n_cells,dtype=np.int32) + mesh.cell_data['partitioned'][newindices] = mesh2.cell_data['partitioned'] + 1 + messh.cell_data['partitioned'] = mesh.cell_data['partitioned'] + + + + +def partition(mesh, numcores): + # ============================================================================= + # define partioner + # ============================================================================= + + libpath = os.getcwd().split("OpenSeesProjects")[0] + "OpenSeesProjects/" + "MeshGenerator/lib" + if os.name == 'nt': + metis_partition_lib = ctypes.CDLL(f'{libpath}/Partitioner.dll') + if os.name == 'posix': + metis_partition_lib = ctypes.CDLL(f'{libpath}/libPartitioner.so') + + # Define function argument and return types + metis_partition_lib.Partition.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_int32), ctypes.c_int, ctypes.POINTER(ctypes.c_int32)] + metis_partition_lib.Partition.restype = ctypes.c_int + + numcells = mesh.n_cells + numpoints = mesh.n_points + numweights = 1 + cells = np.array(mesh.cells.reshape(-1, 9), dtype=np.int32) + cells = cells[:,1:] + cellspointer = cells.flatten().ctypes.data_as(ctypes.POINTER(ctypes.c_int32)) + partitioned = np.empty(numcells, dtype=np.int32) + partitionedpointer = partitioned.ctypes.data_as(ctypes.POINTER(ctypes.c_int32)) + metis_partition_lib.Partition(numcells, numpoints, cellspointer, numcores, partitionedpointer) + mesh.cell_data['partitioned'] = partitioned + + + +def kdtree_partition(grid, n_partitions=4): + blocks = grid.partition(n_partitions, generate_global_id=True, as_composite=True) + i = 0 + grid["partitioned"] = np.zeros(grid.n_cells,dtype=np.int32) + for block in blocks: + grid["partitioned"][block["vtkGlobalCellIds"]] = i + i += 1 + + + + + + +def DRM_PML_Foundation_Meshgenrator(meshinfo) : + xwidth, ywidth, zwidth = (meshinfo["xwidth"], meshinfo["ywidth"], meshinfo["zwidth"]) + eps = 1e-6 + Xmeshsize, Ymeshsize, Zmeshsize = (meshinfo["Xmeshsize"], meshinfo["Ymeshsize"], meshinfo["Zmeshsize"]) + PMLThickness = np.array([meshinfo["PMLThickness"][0], meshinfo["PMLThickness"][1], meshinfo["PMLThickness"][2]]) ; # thickness of the each PML layer + numPMLLayers = meshinfo["numPMLLayers"] ; # number of PML layers + PMLTotalThickness = PMLThickness * numPMLLayers ; # total thickness of the PML layers + DRMThickness = np.array([meshinfo["DRMThickness"][0], meshinfo["DRMThickness"][1], meshinfo["DRMThickness"][2]]) ; # thickness of the DRM layers + numDrmLayers = meshinfo["numDrmLayers"] ; # number of DRM layers + DRMTotalThickness = DRMThickness * numDrmLayers ; # total thickness of the DRM layers + padLayers = numPMLLayers + numDrmLayers ; # number of layers to pad the meshgrid + padThickness = PMLTotalThickness + DRMThickness ; # thickness of the padding layers + reg_num_cores = meshinfo["reg_num_cores"] + DRM_num_cores = meshinfo["DRM_num_cores"] + PML_num_cores = meshinfo["PML_num_cores"] + structurecores = meshinfo["Structure_num_cores"] + Dir = meshinfo["Dir"] + OutputDir = meshinfo["OutputDir"] + AbsorbingElements = meshinfo["AbsorbingElements"] + DRMfile = meshinfo["DRMfile"] + foundationBlocks = meshinfo["foundationBlocks"] + pilelist = meshinfo["pilelist"] + EmbeddeFoundation = meshinfo["EmbeddedFoundationDict"] + HaveEmbeddingFoundation = meshinfo["EmbeddingFoundation"] + HaveFoundation = meshinfo["HaveFoundation"] + HavePiles = meshinfo["HavePiles"] + HaveStructure = meshinfo["HaveStructure"] + AttachFoundation = meshinfo["AttachFoundation"] + DRMinformation = meshinfo["DRMinformation"] + PartitionAlgorithm = meshinfo["PartitionAlgorithm"] + + + # ============================================================================= + if HaveStructure != "YES": + structurecores = 0 + + # create a dictionary for meshing information + info = { + "Foundation": 1, + "RegularDomain": 2, + "DRMDomain": 3, + "PMLDomain": 4, + } + DomainNames = ["Foundation","Soil", "Soil", "PML"] + + + pileboxinfo = {} + + if not os.path.exists(Dir): + os.makedirs(Dir) + else : + # remove the files in the directory + for file in os.listdir(Dir): + try: + os.remove(os.path.join(Dir,file)) + except: + continue + + # adding different plots + meshplotdir = Dir + "/../meshplots" + if not os.path.exists(meshplotdir): + os.makedirs(meshplotdir) + + + # ============================================================================= + # meshing + # ============================================================================= + x = np.arange(-xwidth/2., xwidth/2.+eps, Xmeshsize) + y = np.arange(-ywidth/2., ywidth/2.+eps, Ymeshsize) + z = np.arange(-zwidth, 0+eps, Zmeshsize) + + # padding x and y for PML and DRM layers + x = np.pad(x, (numDrmLayers,numDrmLayers), "linear_ramp", end_values=(x.min()-DRMTotalThickness[0], x.max()+DRMTotalThickness[0])) + y = np.pad(y, (numDrmLayers,numDrmLayers), "linear_ramp", end_values=(y.min()-DRMTotalThickness[1], y.max()+DRMTotalThickness[1])) + z = np.pad(z, (numDrmLayers,0), "linear_ramp", end_values=(z.min()-DRMTotalThickness[2])) + + # padding the x and y for PML and PML layers + x = np.pad(x, (numPMLLayers,numPMLLayers), "linear_ramp", end_values=(x.min()-PMLTotalThickness[0], x.max()+PMLTotalThickness[0])) + y = np.pad(y, (numPMLLayers,numPMLLayers), "linear_ramp", end_values=(y.min()-PMLTotalThickness[1], y.max()+PMLTotalThickness[1])) + z = np.pad(z, (numPMLLayers,0), "linear_ramp", end_values=(z.min()-PMLTotalThickness[2])) + + + x, y, z = np.meshgrid(x, y, z, indexing='ij') + + mesh = pv.StructuredGrid(x, y, z) + mesh.cell_data['Domain'] = np.ones(mesh.n_cells,dtype=np.int8)*info["RegularDomain"] + # ============================================================================= + # sperate embedding layer + # ============================================================================= + if HaveEmbeddingFoundation=="YES": + cube = pv.Cube(bounds=[EmbeddeFoundation["xmin"],EmbeddeFoundation["xmax"],EmbeddeFoundation["ymin"],EmbeddeFoundation["ymax"],EmbeddeFoundation["zmin"],EmbeddeFoundation["zmax"]]) + mesh = mesh.clip_box(cube,invert=True,crinkle=True,progress_bar = False) + mesh.clear_data() + mesh.cell_data['Domain'] = np.ones(mesh.n_cells,dtype=np.int8)*info["RegularDomain"] + # ============================================================================= + # Add foundation blocks + # ============================================================================= + if HaveFoundation == "YES": + for i,block in enumerate(foundationBlocks): + xBLOCK = np.arange(block["xmin"], block["xmax"]+eps, block["Xmeshsize"]) + yBLOCK = np.arange(block["ymin"], block["ymax"]+eps, block["Ymeshsize"]) + zBLOCK = np.arange(block["zmin"], block["zmax"]+eps, block["Zmeshsize"]) + xBLOCK, yBLOCK, zBLOCK = np.meshgrid(xBLOCK, yBLOCK, zBLOCK, indexing='ij') + if i == 0: + foundation = pv.StructuredGrid(xBLOCK, yBLOCK, zBLOCK) + else: + foundation = foundation.merge(pv.StructuredGrid(xBLOCK, yBLOCK, zBLOCK),merge_points=False,tolerance=1e-6,progress_bar = False) + + + + # ============================================================================= + # adding piles + # ============================================================================= + pilenodes = np.zeros((len(pilelist)*2,3)) + pileelement = np.zeros((len(pilelist),3),dtype=int) + for i in range(len(pilelist)): + pilenodes[i*2] = [pilelist[i]["xtop"],pilelist[i]["ytop"],pilelist[i]["ztop"]] + pilenodes[i*2+1] = [pilelist[i]["xbottom"],pilelist[i]["ybottom"],pilelist[i]["zbottom"]] + pileelement[i] = [2,2*i,2*i+1] + celltypes = np.ones(pileelement.shape[0],dtype= int) * pv.CellType.LINE + piles = pv.UnstructuredGrid(pileelement.tolist(),celltypes.tolist(),pilenodes.tolist()) + + + pl = pv.Plotter() + if HavePiles == "YES": + pl.add_mesh(piles, show_edges=True, color = "r" , line_width=4.0,) + if HaveFoundation == "YES": + pl.add_mesh(foundation,show_edges=True, color="gray", opacity=0.5 ) + pl.add_mesh(mesh, opacity=0.5 ) + pl.camera_position = 'xz' + pl.export_html(meshplotdir+"/pile.html") + # pl.show() + + + # merge the piles and foundation + if HaveFoundation == "YES" and HavePiles == "YES": + foundation2 = foundation.merge(piles,merge_points=False,tolerance=1e-6,progress_bar = False) + foundationbounds = foundation2.bounds + pileboxinfo["xmin"] = foundationbounds[0] + pileboxinfo["xmax"] = foundationbounds[1] + pileboxinfo["ymin"] = foundationbounds[2] + pileboxinfo["ymax"] = foundationbounds[3] + pileboxinfo["zmin"] = foundationbounds[4] + pileboxinfo["zmax"] = foundationbounds[5] + pileboxinfo["Havepilebox"] = True + del foundation2 + else : + pileboxinfo["Havepilebox"] = False + if HaveFoundation == "YES" : + foundation.cell_data['Domain'] = np.ones(foundation.n_cells,dtype=np.int8)*info["Foundation"] + if AttachFoundation == "YES": + mesh = mesh.merge(foundation,merge_points=True,tolerance=1e-6,progress_bar = False) + else: + mesh = mesh.merge(foundation,merge_points=False,tolerance=1e-6,progress_bar = False) + + # delete foundation2 + # print("Havepilebox:",pileboxinfo["Havepilebox"]) + + # ============================================================================= + # sperate PML layer + # ============================================================================= + xmin = x.min() + PMLTotalThickness[0] + xmax = x.max() - PMLTotalThickness[0] + ymin = y.min() + PMLTotalThickness[1] + ymax = y.max() - PMLTotalThickness[1] + zmin = z.min() + PMLTotalThickness[2] + zmax = z.max() + PMLTotalThickness[2] + cube = pv.Cube(bounds=[xmin,xmax,ymin,ymax,zmin,1000.0]) + PML = mesh.clip_box(cube,invert=True, crinkle=True, progress_bar = False) + reg = mesh.clip_box(cube,invert=False,crinkle=True, progress_bar = False) + + + + + + + # now find DRM layer + indices = reg.find_cells_within_bounds([xmin + DRMTotalThickness[0] + eps, + xmax - DRMTotalThickness[0] - eps, + ymin + DRMTotalThickness[1] + eps, + ymax - DRMTotalThickness[1] - eps, + zmin + DRMTotalThickness[2] + eps, + zmax + DRMTotalThickness[2] + eps]) + + + + + # now create complemntary indices for DRM + DRMindices = np.ones(reg.n_cells,dtype=bool) + DRMindices[indices] = False + DRMindices = np.where(DRMindices)[0] + + + + # reg.cell_data['Domain'] = np.ones(reg.n_cells,dtype=np.int8)*info["DRMDomain"] + # reg.cell_data['Domain'][indices] = info["RegularDomain"] + + reg.cell_data['Domain'][DRMindices] = info["DRMDomain"] + PML.cell_data['Domain'] = np.ones(PML.n_cells,dtype=np.int8)*info["PMLDomain"] + reg.cell_data['partitioned'] = np.zeros(reg.n_cells,dtype=np.int32) + PML.cell_data['partitioned'] = np.zeros(PML.n_cells,dtype=np.int32) + + # partitioning regular mesh + regular = reg.extract_cells(indices,progress_bar=False) + DRM = reg.extract_cells(DRMindices,progress_bar=False) + regularvtkOriginalCellIds = regular["vtkOriginalCellIds"] + DRMvtkOriginalCellIds = DRM["vtkOriginalCellIds"] + + if reg_num_cores > 1: + if PartitionAlgorithm.lower() == "metis": + print("Using Metis partitioning") + partition_with_pilebox(regular,reg_num_cores,pileboxinfo,tol=10) + if PartitionAlgorithm.lower() == "kd-tree": + print("Using KD-Tree partitioning") + kdtree_partition(regular,reg_num_cores) + + if DRM_num_cores > 1: + if PartitionAlgorithm.lower() == "metis": + partition(DRM,DRM_num_cores) + if PartitionAlgorithm.lower() == "kd-tree": + kdtree_partition(DRM,DRM_num_cores) + + if PML_num_cores > 1: + if PartitionAlgorithm.lower() == "metis": + partition(PML,PML_num_cores) + if PartitionAlgorithm.lower() == "kd-tree": + PMLcopy = PML.copy() + kdtree_partition(PMLcopy,PML_num_cores) + PML.cell_data['partitioned'] = PMLcopy["partitioned"] + del PMLcopy + + reg.cell_data['partitioned'][regularvtkOriginalCellIds] = regular.cell_data['partitioned'] + reg.cell_data['partitioned'][DRMvtkOriginalCellIds] = DRM.cell_data['partitioned'] + reg_num_cores + PML.cell_data['partitioned'] = PML.cell_data['partitioned'] + reg_num_cores + DRM_num_cores + + + + + # regular.plot(scalars="partitioned",show_edges=True) + # DRM.plot(scalars="partitioned",show_edges=True) + # PML.plot(scalars="partitioned",show_edges=True) + + # merging PML and regular mesh to create a single mesh + mesh = reg.merge(PML,merge_points=False,tolerance=1e-6,progress_bar = False) + # stack the regular and DRM partitioned data + # mesh.plot(scalars="partitioned",show_edges=True) + + # mapping between PML and regular mesh on the boundary + mapping = mesh.clean(produce_merge_map=True)["PointMergeMap"] + regindicies = np.where(mapping[PML.n_points:]0)[0] + + + mesh["matTag"] = np.ones(mesh.n_cells,dtype=np.uint8) + + + + + + + # define the ASDA absorbing elements + if AbsorbingElements.lower() == "asda": + + mesh = mesh.clean(tolerance=1e-6,remove_unused_points=False) + mesh["ASDA_type"] = np.zeros(mesh.n_cells,dtype=np.uint8) + + ASDAelem_type = { + "B" :1 , + "L" :2 , + "R" :3 , + "F" :4 , + "K" :5 , + "BL" :6 , + "BR" :7 , + "BF" :8 , + "BK" :9 , + "LF" :10, + "LK" :11, + "RF" :12, + "RK" :13, + "BLF" :14, + "BLK" :15, + "BRF" :16, + "BRK" :17, + } + + ASDAelem_typereverse = { + 1 :"B" , + 2 :"L" , + 3 :"R" , + 4 :"F" , + 5 :"K" , + 6 :"BL" , + 7 :"BR" , + 8 :"BF" , + 9 :"BK" , + 10:"LF" , + 11:"LK" , + 12:"RF" , + 13:"RK" , + 14:"BLF" , + 15:"BLK" , + 16:"BRF" , + 17:"BRK" , + } + xmin, xmax, ymin, ymax, zmin, zmax = reg.bounds + ASDA_xwidth = xmax - xmin + ASDA_ywidth = ymax - ymin + ASDA_zwidth = zmax - zmin + # print("ASDA_xwidth", ASDA_xwidth) + # print("ASDA_ywidth", ASDA_ywidth) + # print("ASDA_zwidth", ASDA_zwidth) + + i = 0 + for ele_center in mesh.cell_centers().points: + # check if the element is in the left or rightside + if ele_center[0] < (-ASDA_xwidth/2.): + # it is in the left side + # check if it is in the front or back + if ele_center[1] < (-ASDA_ywidth/2.): + # it is in the back + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BLK"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["LK"] + elif ele_center[1] > (ASDA_ywidth/2.): + # it is in the front + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BLF"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["LF"] + else: + # it is in the middle + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BL"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["L"] + + elif ele_center[0] > (ASDA_xwidth/2.): + # it is in the right side + # check if it is in the front or back + if ele_center[1] < (-ASDA_ywidth/2.): + # it is in the back + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BRK"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["RK"] + elif ele_center[1] > (ASDA_ywidth/2.): + # it is in the front + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BRF"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["RF"] + else: + # it is in the middle + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BR"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["R"] + else: + # it is in the middle + # check if it is in the front or back + if ele_center[1] < (-ASDA_ywidth/2.): + # it is in the back + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BK"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["K"] + elif ele_center[1] > (ASDA_ywidth/2.): + # it is in the front + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BF"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["F"] + else: + # it is in the middle + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["B"] + + i += 1 + + + if AbsorbingElements.lower() == "normal": + mesh = mesh.clean(tolerance=1e-6,remove_unused_points=False) + # ============================================================================= + # write the mesh + # ============================================================================= + min_core = mesh.cell_data['partitioned'].min() + max_core = mesh.cell_data['partitioned'].max() + + + + # write the mesh nodes + for core in range(min_core,max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Nodes" + str(core + structurecores) + ".tcl", "w") + + for i in range(tmp.n_points): + f.write(f"node [expr int($StructureMaxNodeTag + 1 + {tmp['vtkOriginalPointIds'][i]})] {tmp.points[i][0]} {tmp.points[i][1]} {tmp.points[i][2]}\n") + f.close() + + + # writing the mesh elements + if AbsorbingElements == "ASDA": + # writing the mesh elements + for core in range(min_core,max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Elements" + str(core + structurecores) + ".tcl", "w") + if core >= reg_num_cores + DRM_num_cores: + for eletag in range(tmp.n_cells): + tmpeletag = "[expr int($StructureMaxEleTag + 1 + " + str(tmp['vtkOriginalCellIds'][eletag]) + ")]" + tmpNodeTags = [f"[expr int($StructureMaxNodeTag + 1 + {x})]" for x in tmp['vtkOriginalPointIds'][tmp.get_cell(eletag).point_ids]] + Domainname = DomainNames[tmp['Domain'][eletag]-1] + if Domainname == "PML": + Domainname = "Absorbing" + tmpmatTag = f"${Domainname}matTag{tmp['matTag'][eletag]}" + tmpASDA_type = ASDAelem_typereverse[tmp['ASDA_type'][eletag]] + f.write(f"eval \"element $elementType {tmpeletag} {' '.join(tmpNodeTags)} {tmpmatTag} {tmpASDA_type}\" \n") + else: + for eletag in range(tmp.n_cells): + tmpeletag = "[expr int($StructureMaxEleTag + 1 + " + str(tmp['vtkOriginalCellIds'][eletag]) + ")]" + tmpNodeTags = [f"[expr int($StructureMaxNodeTag + 1 + {x})]" for x in tmp['vtkOriginalPointIds'][tmp.get_cell(eletag).point_ids]] + Domainname = DomainNames[tmp['Domain'][eletag]-1] + tmpmatTag = f"${Domainname}matTag{tmp['matTag'][eletag]}" + f.write(f"eval \"element $elementType {tmpeletag} {' '.join(tmpNodeTags)} {tmpmatTag}\" \n") + f.close() + else : + for core in range(min_core,max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Elements" + str(core + structurecores) + ".tcl", "w") + for eletag in range(tmp.n_cells): + tmpeletag = "[expr int($StructureMaxEleTag + 1 + " + str(tmp['vtkOriginalCellIds'][eletag]) + ")]" + tmpNodeTags = [f"[expr int($StructureMaxNodeTag + 1 + {x})]" for x in tmp['vtkOriginalPointIds'][tmp.get_cell(eletag).point_ids]] + Domainname = DomainNames[tmp['Domain'][eletag]-1] + if Domainname == "PML": + Domainname = "Absorbing" + tmpmatTag = f"${Domainname}matTag{tmp['matTag'][eletag]}" + f.write(f"eval \"element $elementType {tmpeletag} {' '.join(tmpNodeTags)} {tmpmatTag}\" \n") + f.close() + + if AbsorbingElements == "PML": + # writing the boundary files + for core in range(reg_num_cores + DRM_num_cores , max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Boundary" + str(core + structurecores) + ".tcl", "w") + for i in range(tmp.n_points): + if tmp["boundary"][i] != -1: + x,y,z = tmp.points[i] + nodeTag1 = tmp['vtkOriginalPointIds'][i] + nodeTag2 = tmp['boundary'][i] + f.write(f"node [expr int($StructureMaxNodeTag + 1 +{nodeTag2})] {str(x)} {str(y)} {str(z)}\n") + f.write(f"equalDOF [expr int($StructureMaxNodeTag + 1 + {nodeTag2})] [expr int($StructureMaxNodeTag + 1 +{nodeTag1})] 1 2 3\n") + + + + # ============================================================================= + # printing information + # ============================================================================= + print(f"Number of regular cores: {reg_num_cores}") + print(f"Number of DRM cores: {DRM_num_cores}") + print(f"Number of PML cores: {PML_num_cores}") + print(f"Number of regular elements: {regular.n_cells} roughly {int(regular.n_cells/reg_num_cores)} each core") + print(f"Number of DRM elements: {DRM.n_cells} roughly {int(DRM.n_cells/DRM_num_cores)} each core") + print(f"Number of PML elements: {PML.n_cells} roughly {int(PML.n_cells/PML_num_cores)} each core") + print(f"Number of total elements: {mesh.n_cells}") + print(f"Number of total points: {mesh.n_points}") + print(f"Number of cores: {max_core-min_core+1}") + print(f"Number of PML nodes: {PML.n_points}") + print(f"Number of regular nodes: {regular.n_points}") + print(f"Number of DRM nodes: {DRM.n_points}") + if AbsorbingElements == "PML": + print(f"Number of MP constraints: {regindicies.size}") + + # calculating number of surface points on the boundaries + eps = 1e-2 + bounds = np.array(mesh.bounds)-np.array([-eps,eps,-eps,eps,-eps,-eps]) + cube = pv.Cube(bounds=bounds) + # points outside the cube + selected = mesh.select_enclosed_points(cube,inside_out = True) + pts = mesh.extract_points( + selected['SelectedPoints'].view(bool), + include_cells=False, + ) + print(f"number of sp constriants: {pts.n_points*9}") + + + + + + + # f = h5py.File('./DRMloadSmall.h5drm', 'r') + if DRMinformation["DRM_Location"].lower() == "local": + f = h5py.File(DRMfile, 'r') + pts = f["DRM_Data"]["xyz"][:] + internal = f["DRM_Data"]["internal"][:] + xyz0 = f["DRM_QA_Data"]["xyz"][:] + + + + + if DRMinformation["DRM_Provider_Software"].lower() == "shakermaker" and DRMinformation["DRM_Location"].lower() == "local": + pts = pts - xyz0 + pts[:,[0,1]] = pts[:,[1,0]] + pts[:,2] = -pts[:,2] + pts = pts * DRMinformation["crd_scale"] + + + if DRMinformation["DRM_Location"].lower() == "local": + set1 = DRM.points + set2 = pts + + # check that set2 is a subset of set1 + tree = KDTree(set1) + + tolernace = 1e-6 + issubset = True + for point in set2: + dist, ind = tree.query(point, distance_upper_bound=tolernace) + if dist > tolernace: + issubset = False + break + + if issubset: + print("The DRM nodes in the loading file are subset of the DRM mesh") + else: + print("The DRM nodes in the loading file are not subset of the DRM mesh") + print("Please check the DRM nodes in the loading file or change the DRM mesh") + else : + print("The DRM nodes are not checked for subset") + + + + + + pl = pv.Plotter() + # plot the DRM layer with the DRM nodes loaded from the DRM file + pl.add_mesh(DRM,scalars="partitioned",show_edges=True) + if DRMinformation["DRM_Location"].lower() == "local": + pl.add_points(pts, color='r') + pl.export_html(meshplotdir+"/DRM.html") + pl.clear() + + + # plot the regular mesh with the internal nodes loaded from the DRM file + pl.add_mesh(regular,scalars="partitioned",show_edges=True) + pl.export_html(meshplotdir+"/Regular.html") + pl.clear() + + # plot the PML mesh + pl.add_mesh(PML,scalars="partitioned",show_edges=True) + pl.export_html(meshplotdir+"/PML.html") + pl.clear() + + + # plot the total mesh + pl.add_mesh(mesh,scalars="partitioned",show_edges=True) + pl.export_html(meshplotdir+"/Total_partitioned.html") + pl.clear() + + # plot the total mesh with domain scalars + pl.add_mesh(mesh,scalars="Domain",show_edges=True) + pl.export_html(meshplotdir+"/Total_domain.html") + pl.clear() + + + if AbsorbingElements == "ASDA": + pl.add_mesh(mesh,scalars="ASDA_type",show_edges=True,cmap="tab20") + pl.export_html(meshplotdir+"/ASDA_total.html") + pl.clear() + + # filter the mesh with domain pml + indices = mesh['Domain'] == info["PMLDomain"] + grid = mesh.extract_cells(indices) + pl.add_mesh(grid,scalars="ASDA_type",show_edges=True,cmap="tab20b") + pl.export_html(meshplotdir+"/ASDA_PML.html") + pl.clear() + + pl.close() + + # save the mesh + mesh.save(os.path.join(OutputDir,"mesh.vtk"),binary=True) + + + + # return thenumber of elements + return (mesh.n_cells, mesh.n_points) + + + + + + + + +# %% diff --git a/modules/tools/DRM/MeshGenerator/Structure.py b/modules/tools/DRM/MeshGenerator/Structure.py new file mode 100644 index 000000000..954781558 --- /dev/null +++ b/modules/tools/DRM/MeshGenerator/Structure.py @@ -0,0 +1,774 @@ +# %% +import pyvista as pv +import numpy as np +import ctypes +import os +import sys +# ============================================================================= +# information +# ============================================================================= +# getting the meshing information from the command line + +# xwidth, ywidth, zwidth = (float(sys.argv[1]), float(sys.argv[2]), float(sys.argv[3])) +# eps = 1e-6 +# Xmeshsize, Ymeshsize, Zmeshsize = (float(sys.argv[4]), float(sys.argv[5]), float(sys.argv[6])) +# PMLThickness = np.array([float(sys.argv[7]), float(sys.argv[8]), float(sys.argv[9])]) ; # thickness of the each PML layer +# numPMLLayers = int(sys.argv[10]) ; # number of PML layers +# PMLTotalThickness = PMLThickness * numPMLLayers ; # total thickness of the PML layers +# DRMThickness = np.array([float(sys.argv[11]), float(sys.argv[12]), float(sys.argv[13])]) ; # thickness of the DRM layers +# numDrmLayers = int(sys.argv[14]) ; # number of DRM layers +# DRMTotalThickness = DRMThickness * numDrmLayers ; # total thickness of the DRM layers +# padLayers = numPMLLayers + numDrmLayers ; # number of layers to pad the meshgrid +# padThickness = PMLTotalThickness + DRMThickness ; # thickness of the padding layers +# reg_num_cores = int(sys.argv[15]) +# DRM_num_cores = int(sys.argv[16]) +# PML_num_cores = int(sys.argv[17]) +# Dir = sys.argv[18] +# OutputDir = sys.argv[19] +# pileboxcenterx = float(sys.argv[20]) +# pileboxcentery = float(sys.argv[21]) +# pileboxcenterz = float(sys.argv[22]) +# pileboxlx = float(sys.argv[23]) +# pileboxly = float(sys.argv[24]) +# pileboxdepth = float(sys.argv[25]) +# AbsorbingElements = sys.argv[26] +# DRMfile = sys.argv[27] + + + +xwidth = 100 +ywidth = 100 +zwidth = 40 +eps = 1e-6 +Xmeshsize, Ymeshsize, Zmeshsize = (2.5, 2.5, 2.5) +PMLThickness = np.array([Xmeshsize, Ymeshsize, Zmeshsize]) ; # thickness of the each PML layer +numPMLLayers = 2 ; # number of PML layers +PMLTotalThickness = PMLThickness * numPMLLayers ; # total thickness of the PML layers +DRMThickness = np.array([Xmeshsize, Ymeshsize, Zmeshsize]) ; # thickness of the DRM layers +numDrmLayers = 1 ; # number of DRM layers +DRMTotalThickness = DRMThickness * numDrmLayers ; # total thickness of the DRM layers +padLayers = numPMLLayers + numDrmLayers ; # number of layers to pad the meshgrid +padThickness = PMLTotalThickness + DRMThickness ; # thickness of the padding layers +reg_num_cores = 3 +DRM_num_cores = 3 +PML_num_cores = 3 +Dir = "OpenSeesMesh" +OutputDir = "results" +pileboxcenterx = 0 +pileboxcentery = 0 +pileboxcenterz = 0 +pileboxlx = 20 +pileboxly = 20 +pileboxdepth = 10.0 +AbsorbingElements = "PML" +DRMfile = "SurfaceWave.h5drm" + + +# # print the input information +# print(f"xwidth: {xwidth}") +# print(f"ywidth: {ywidth}") +# print(f"zwidth: {zwidth}") +# print(f"Xmeshsize: {Xmeshsize}") +# print(f"Ymeshsize: {Ymeshsize}") +# print(f"Zmeshsize: {Zmeshsize}") +# print(f"PMLThickness: {PMLThickness}") +# print(f"numPMLLayers: {numPMLLayers}") +# print(f"PMLTotalThickness: {PMLTotalThickness}") +# print(f"DRMThickness: {DRMThickness}") +# print(f"numDrmLayers: {numDrmLayers}") +# print(f"DRMTotalThickness: {DRMTotalThickness}") +# print(f"padLayers: {padLayers}") +# print(f"padThickness: {padThickness}") +# print(f"reg_num_cores: {reg_num_cores}") +# print(f"DRM_num_cores: {DRM_num_cores}") +# print(f"PML_num_cores: {PML_num_cores}") +# print(f"Dir: {Dir}") +# print(f"OutputDir: {OutputDir}") +# print(f"pileboxcenterx: {pileboxcenterx}") +# print(f"pileboxcentery: {pileboxcentery}") +# print(f"pileboxcenterz: {pileboxcenterz}") +# print(f"pileboxlx: {pileboxlx}") +# print(f"pileboxly: {pileboxly}") +# print(f"pileboxdepth: {pileboxdepth}") +# print(f"AbsorbingElements: {AbsorbingElements}") +# print(f"DRMfile: {DRMfile}") + +# create a dictionary for meshing information +info = { + "RegularDomain": 1, + "DRMDomain": 2, + "PMLDomain": 3, +} + +pileboxinfo = {} +# adding different plots +meshplotdir = OutputDir + "/meshplots" +if not os.path.exists(meshplotdir): + os.makedirs(meshplotdir) +# ============================================================================= +# buiding an structure mesh +# ============================================================================= +embededDepth = { + "xmax": 20, + "xmin": -20, + "ymax": 15, + "ymin": -15, + "zmax": 0, + "zmin": -5, +} + + + +foundationBlocks = [] + +cnetersx = np.arange(-15,15+eps,7.5) +cnetersy = np.arange(-7.5,7.5+eps,7.5) +cnetersx,cnetersy = np.meshgrid(cnetersx,cnetersy) +# make them tuples +centers = np.vstack((cnetersx.flatten(),cnetersy.flatten())).T + +for i in range(len(centers)): + blockinfo = { + "xmax" : centers[i][0] + 2, + "xmin" : centers[i][0] - 2, + "ymax" : centers[i][1] + 2, + "ymin" : centers[i][1] - 2, + "zmax" : -3, + "zmin" : -5, + "Xmeshsize": 1.0, + "Ymeshsize": 1.0, + "Zmeshsize": 1.0, + } + foundationBlocks.append(blockinfo) + +# ============================================================================= +# define the piles +# ============================================================================= +pilelist = [] +for ii in range(len(centers)): + + for tu in [(1,1),(1,-1),(-1,1),(-1,-1)]: + i,j = tu + pilelist.append({ + "xtop": centers[ii][0] + i, + "ytop": centers[ii][1] + j, + "ztop": -4, + "xbottom": centers[ii][0] + i, + "ybottom": centers[ii][1] + j, + "zbottom": -10.0, + }) + + + +# ============================================================================= +# define partioner +# ============================================================================= +# change the directory to the current directory +os.chdir(os.path.dirname(os.path.abspath(__file__))) +libpath = os.getcwd().split("OpenSeesProjects")[0] + "OpenSeesProjects/" + "MeshGenerator/lib" +print(libpath) +if os.name == 'nt': + metis_partition_lib = ctypes.CDLL(f'{libpath}/Partitioner.dll') +if os.name == 'posix': + metis_partition_lib = ctypes.CDLL(f'{libpath}/libPartitioner.so') + +# Define function argument and return types +metis_partition_lib.Partition.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_int32), ctypes.c_int, ctypes.POINTER(ctypes.c_int32)] +metis_partition_lib.Partition.restype = ctypes.c_int + +def partition(mesh, numcores): + numcells = mesh.n_cells + numpoints = mesh.n_points + numweights = 1 + cells = np.array(mesh.cells.reshape(-1, 9), dtype=np.int32) + cells = cells[:,1:] + cellspointer = cells.flatten().ctypes.data_as(ctypes.POINTER(ctypes.c_int32)) + partitioned = np.empty(numcells, dtype=np.int32) + partitionedpointer = partitioned.ctypes.data_as(ctypes.POINTER(ctypes.c_int32)) + metis_partition_lib.Partition(numcells, numpoints, cellspointer, numcores, partitionedpointer) + mesh.cell_data['partitioned'] = partitioned + + + + +def partition_with_pilebox(messh,numcores, pileboxinfo,tol=1e-6): + # seperate the core of the pilebox + mesh = messh.copy() + eps = 1e-6 + # check if the xmin,xmax, ymin ymax ,zmin zmax is exist in the infokeys + if "xmin" in pileboxinfo.keys(): + pilebox_xmin = pileboxinfo["xmin"] + if "xmax" in pileboxinfo.keys(): + pilebox_xmax = pileboxinfo["xmax"] + if "ymin" in pileboxinfo.keys(): + pilebox_ymin = pileboxinfo["ymin"] + if "ymax" in pileboxinfo.keys(): + pilebox_ymax = pileboxinfo["ymax"] + if "zmin" in pileboxinfo.keys(): + pilebox_zmin = pileboxinfo["zmin"] + if "zmax" in pileboxinfo.keys(): + pilebox_zmax = pileboxinfo["zmax"] + if "center" in pileboxinfo.keys(): + pilebox_xmin = pileboxinfo["center"][0] - pileboxinfo["lx"]/2 + eps + pilebox_xmax = pileboxinfo["center"][0] + pileboxinfo["lx"]/2 - eps + pilebox_ymin = pileboxinfo["center"][1] - pileboxinfo["ly"]/2 + eps + pilebox_ymax = pileboxinfo["center"][1] + pileboxinfo["ly"]/2 - eps + pilebox_zmin = pileboxinfo["center"][2] - pileboxinfo["depth"] + eps + pilebox_zmax = pileboxinfo["center"][2] + tol + eps + + + # find the cells that are inside the pilebox + cube = pv.Cube(bounds=[pilebox_xmin,pilebox_xmax,pilebox_ymin,pilebox_ymax,pilebox_zmin,pilebox_zmax]) + # extract the cells that are outside the pilebox + indices = mesh.find_cells_within_bounds(cube.bounds) + + # crete newindices for outside the pilebox + newindices = np.ones(mesh.n_cells,dtype=bool) + newindices[indices] = False + newindices = np.where(newindices)[0] + + + # extract the cells that are inside the pilebox + mesh2 = mesh.extract_cells(newindices) + # partition the mesh + + if numcores > 2: + partition(mesh2,numcores-1) + if numcores == 2: + mesh2.cell_data['partitioned'] = np.zeros(mesh2.n_cells,dtype=np.int32) + + mesh.cell_data['partitioned'] = np.zeros(mesh.n_cells,dtype=np.int32) + mesh.cell_data['partitioned'][newindices] = mesh2.cell_data['partitioned'] + 1 + messh.cell_data['partitioned'] = mesh.cell_data['partitioned'] + + +# ============================================================================= +# meshing +# ============================================================================= +x = np.arange(-xwidth/2., xwidth/2.+eps, Xmeshsize) +y = np.arange(-ywidth/2., ywidth/2.+eps, Ymeshsize) +z = np.arange(-zwidth, 0+eps, Zmeshsize) + +# padding x and y for PML and DRM layers +x = np.pad(x, (numDrmLayers,numDrmLayers), "linear_ramp", end_values=(x.min()-DRMTotalThickness[0], x.max()+DRMTotalThickness[0])) +y = np.pad(y, (numDrmLayers,numDrmLayers), "linear_ramp", end_values=(y.min()-DRMTotalThickness[1], y.max()+DRMTotalThickness[1])) +z = np.pad(z, (numDrmLayers,0), "linear_ramp", end_values=(z.min()-DRMTotalThickness[2])) + +# padding the x and y for PML and PML layers +x = np.pad(x, (numPMLLayers,numPMLLayers), "linear_ramp", end_values=(x.min()-PMLTotalThickness[0], x.max()+PMLTotalThickness[0])) +y = np.pad(y, (numPMLLayers,numPMLLayers), "linear_ramp", end_values=(y.min()-PMLTotalThickness[1], y.max()+PMLTotalThickness[1])) +z = np.pad(z, (numPMLLayers,0), "linear_ramp", end_values=(z.min()-PMLTotalThickness[2])) + +# %% +x, y, z = np.meshgrid(x, y, z, indexing='ij') + +mesh = pv.StructuredGrid(x, y, z) +# ============================================================================= +# sperate embedding layer +# ============================================================================= +cube = pv.Cube(bounds=[embededDepth["xmin"],embededDepth["xmax"],embededDepth["ymin"],embededDepth["ymax"],embededDepth["zmin"],embededDepth["zmax"]]) +mesh = mesh.clip_box(cube,invert=True,crinkle=True,progress_bar = True) +mesh.clear_data() +# ============================================================================= +# Add foundation blocks +# ============================================================================= +for i,block in enumerate(foundationBlocks): + xBLOCK = np.arange(block["xmin"], block["xmax"]+eps, block["Xmeshsize"]) + yBLOCK = np.arange(block["ymin"], block["ymax"]+eps, block["Ymeshsize"]) + zBLOCK = np.arange(block["zmin"], block["zmax"]+eps, block["Zmeshsize"]) + xBLOCK, yBLOCK, zBLOCK = np.meshgrid(xBLOCK, yBLOCK, zBLOCK, indexing='ij') + if i == 0: + foundation = pv.StructuredGrid(xBLOCK, yBLOCK, zBLOCK) + else: + foundation = foundation.merge(pv.StructuredGrid(xBLOCK, yBLOCK, zBLOCK),merge_points=False,tolerance=1e-6,progress_bar = True) + + + +# ============================================================================= +# adding piles +# ============================================================================= +pilenodes = np.zeros((len(pilelist)*2,3)) +pileelement = np.zeros((len(pilelist),3),dtype=int) +for i in range(len(pilelist)): + pilenodes[i*2] = [pilelist[i]["xtop"],pilelist[i]["ytop"],pilelist[i]["ztop"]] + pilenodes[i*2+1] = [pilelist[i]["xbottom"],pilelist[i]["ybottom"],pilelist[i]["zbottom"]] + pileelement[i] = [2,2*i,2*i+1] +celltypes = np.ones(pileelement.shape[0],dtype= int) * pv.CellType.LINE +piles = pv.UnstructuredGrid(pileelement.tolist(),celltypes.tolist(),pilenodes.tolist()) + +# %% +pl = pv.Plotter() +pl.add_mesh(piles, show_edges=True, color = "r" , line_width=4.0,) +pl.add_mesh(foundation, color="gray", opacity=0.5 ) +pl.add_mesh(mesh, opacity=0.5 ) +pl.camera_position = 'xz' +pl.export_html(meshplotdir+"/pile.html") +# pl.show() + +# %% +# merge the piles and foundation +foundation = foundation.merge(piles,merge_points=False,tolerance=1e-6,progress_bar = True) +foundationbounds = foundation.bounds +pileboxinfo["xmin"] = foundationbounds[0] +pileboxinfo["xmax"] = foundationbounds[1] +pileboxinfo["ymin"] = foundationbounds[2] +pileboxinfo["ymax"] = foundationbounds[3] +pileboxinfo["zmin"] = foundationbounds[4] +pileboxinfo["zmax"] = foundationbounds[5] +mesh = mesh.merge(foundation,merge_points=False,tolerance=1e-6,progress_bar = True) +# %% +# ============================================================================= +# sperate PML layer +# ============================================================================= +xmin = x.min() + PMLTotalThickness[0] +xmax = x.max() - PMLTotalThickness[0] +ymin = y.min() + PMLTotalThickness[1] +ymax = y.max() - PMLTotalThickness[1] +zmin = z.min() + PMLTotalThickness[2] +zmax = z.max() + PMLTotalThickness[2] +cube = pv.Cube(bounds=[xmin,xmax,ymin,ymax,zmin,zmax]) +PML = mesh.clip_box(cube,invert=True, crinkle=True, progress_bar = True) +reg = mesh.clip_box(cube,invert=False,crinkle=True, progress_bar = True) + + + + + + + + +# reg.plot(show_edges=True,show_grid=True,show_axes=True,show_bounds=True) +pl = pv.Plotter() +pl.add_mesh(reg, show_edges=True, style="wireframe" ) +pl.show() + + + +# %% +# now find DRM layer +indices = reg.find_cells_within_bounds([xmin + DRMTotalThickness[0] + eps, + xmax - DRMTotalThickness[0] - eps, + ymin + DRMTotalThickness[1] + eps, + ymax - DRMTotalThickness[1] - eps, + zmin + DRMTotalThickness[2] + eps, + zmax + DRMTotalThickness[2] + eps]) +# print(DRMTotalThickness) +# print(xmax, xmin, ymax, ymin, zmax, zmin) +# pl = pv.Plotter() +# pl.add_mesh(reg,show_edges=True) +# cube = pv.Cube(bounds=[xmin + DRMTotalThickness[0] + eps, +# xmax - DRMTotalThickness[0] - eps, +# ymin + DRMTotalThickness[1] + eps, +# ymax - DRMTotalThickness[1] - eps, +# zmin + DRMTotalThickness[2] + eps, +# zmax + DRMTotalThickness[2] + eps]) +# pl.add_mesh(cube,show_edges=True) +# pl.show() +# exit +# now create complemntary indices for DRM +DRMindices = np.ones(reg.n_cells,dtype=bool) +DRMindices[indices] = False +DRMindices = np.where(DRMindices)[0] + + + +reg.cell_data['Domain'] = np.ones(reg.n_cells,dtype=np.int8)*info["DRMDomain"] +reg.cell_data['Domain'][indices] = info["RegularDomain"] +PML.cell_data['Domain'] = np.ones(PML.n_cells,dtype=np.int8)*info["PMLDomain"] +reg.cell_data['partitioned'] = np.zeros(reg.n_cells,dtype=np.int32) +PML.cell_data['partitioned'] = np.zeros(PML.n_cells,dtype=np.int32) + + +# partitioning regular mesh +regular = reg.extract_cells(indices,progress_bar=True) +DRM = reg.extract_cells(DRMindices,progress_bar=True) + +if reg_num_cores > 1: + partition_with_pilebox(regular,reg_num_cores,pileboxinfo,tol=10) +if DRM_num_cores > 1: + partition(DRM,DRM_num_cores) +if PML_num_cores > 1: + partition(PML,PML_num_cores) + +reg.cell_data['partitioned'][regular["vtkOriginalCellIds"]] = regular.cell_data['partitioned'] +reg.cell_data['partitioned'][DRM["vtkOriginalCellIds"]] = DRM.cell_data['partitioned'] + reg_num_cores +PML.cell_data['partitioned'] = PML.cell_data['partitioned'] + reg_num_cores + DRM_num_cores + +# %% +# merging PML and regular mesh to create a single mesh +mesh = reg.merge(PML,merge_points=False,tolerance=1e-6,progress_bar = True) + + +# mapping between PML and regular mesh on the boundary +mapping = mesh.clean(produce_merge_map=True)["PointMergeMap"] +regindicies = np.where(mapping[PML.n_points:]0)[0] + +# %% +mesh["matTag"] = np.ones(mesh.n_cells,dtype=np.uint8) + + + + + +# %% +# define the ASDA absorbing elements +if AbsorbingElements == "ASDA": + + mesh = mesh.clean(tolerance=1e-6,remove_unused_points=False) + mesh["ASDA_type"] = np.zeros(mesh.n_cells,dtype=np.uint8) + + ASDAelem_type = { + "B" :1 , + "L" :2 , + "R" :3 , + "F" :4 , + "K" :5 , + "BL" :6 , + "BR" :7 , + "BF" :8 , + "BK" :9 , + "LF" :10, + "LK" :11, + "RF" :12, + "RK" :13, + "BLF" :14, + "BLK" :15, + "BRF" :16, + "BRK" :17, + } + + ASDAelem_typereverse = { + 1 :"B" , + 2 :"L" , + 3 :"R" , + 4 :"F" , + 5 :"K" , + 6 :"BL" , + 7 :"BR" , + 8 :"BF" , + 9 :"BK" , + 10:"LF" , + 11:"LK" , + 12:"RF" , + 13:"RK" , + 14:"BLF" , + 15:"BLK" , + 16:"BRF" , + 17:"BRK" , + } + xmin, xmax, ymin, ymax, zmin, zmax = reg.bounds + ASDA_xwidth = xmax - xmin + ASDA_ywidth = ymax - ymin + ASDA_zwidth = zmax - zmin + print("ASDA_xwidth", ASDA_xwidth) + print("ASDA_ywidth", ASDA_ywidth) + print("ASDA_zwidth", ASDA_zwidth) + + i = 0 + for ele_center in mesh.cell_centers().points: + # check if the element is in the left or rightside + if ele_center[0] < (-ASDA_xwidth/2.): + # it is in the left side + # check if it is in the front or back + if ele_center[1] < (-ASDA_ywidth/2.): + # it is in the back + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BLK"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["LK"] + elif ele_center[1] > (ASDA_ywidth/2.): + # it is in the front + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BLF"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["LF"] + else: + # it is in the middle + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BL"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["L"] + + elif ele_center[0] > (ASDA_xwidth/2.): + # it is in the right side + # check if it is in the front or back + if ele_center[1] < (-ASDA_ywidth/2.): + # it is in the back + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BRK"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["RK"] + elif ele_center[1] > (ASDA_ywidth/2.): + # it is in the front + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BRF"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["RF"] + else: + # it is in the middle + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BR"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["R"] + else: + # it is in the middle + # check if it is in the front or back + if ele_center[1] < (-ASDA_ywidth/2.): + # it is in the back + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BK"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["K"] + elif ele_center[1] > (ASDA_ywidth/2.): + # it is in the front + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["BF"] + else: + # it is in the top + mesh["ASDA_type"][i] = ASDAelem_type["F"] + else: + # it is in the middle + # check if it is in the bottom or top + if ele_center[2] < -ASDA_zwidth: + # it is in the bottom + mesh["ASDA_type"][i] = ASDAelem_type["B"] + + i += 1 +# %% +# ============================================================================= +# write the mesh +# ============================================================================= +if not os.path.exists(Dir): + os.makedirs(Dir) +else : + # remove the files in the directory + for file in os.listdir(Dir): + os.remove(os.path.join(Dir,file)) + + +min_core = mesh.cell_data['partitioned'].min() +max_core = mesh.cell_data['partitioned'].max() + + + +# write the mesh nodes +for core in range(min_core,max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Nodes" + str(core) + ".tcl", "w") + + for i in range(tmp.n_points): + f.write(f"node {tmp['vtkOriginalPointIds'][i]} {tmp.points[i][0]} {tmp.points[i][1]} {tmp.points[i][2]}\n") + f.close() +# %% + +# writing the mesh elements +if AbsorbingElements == "ASDA": + # writing the mesh elements + for core in range(min_core,max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Elements" + str(core) + ".tcl", "w") + if core >= reg_num_cores + DRM_num_cores: + for eletag in range(tmp.n_cells): + f.write(f"eval \"element $elementType {tmp['vtkOriginalCellIds'][eletag]} {' '.join(str(x) for x in tmp['vtkOriginalPointIds'][tmp.get_cell(eletag).point_ids])} $matTag{tmp['matTag'][eletag]} {ASDAelem_typereverse[tmp['ASDA_type'][eletag]]}\" \n") + else: + for eletag in range(tmp.n_cells): + f.write(f"eval \"element $elementType {tmp['vtkOriginalCellIds'][eletag]} {' '.join(str(x) for x in tmp['vtkOriginalPointIds'][tmp.get_cell(eletag).point_ids])} $matTag{tmp['matTag'][eletag]}\" \n") + f.close() +else : + for core in range(min_core,max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Elements" + str(core) + ".tcl", "w") + for eletag in range(tmp.n_cells): + f.write(f"eval \"element $elementType {tmp['vtkOriginalCellIds'][eletag]} {' '.join(str(x) for x in tmp['vtkOriginalPointIds'][tmp.get_cell(eletag).point_ids])} $matTag{tmp['matTag'][eletag]}\" \n") + f.close() + +if AbsorbingElements == "PML": + # writing the boundary files + for core in range(reg_num_cores + DRM_num_cores , max_core+1): + tmp = mesh.extract_cells(np.where(mesh.cell_data['partitioned']==core)[0]) + f = open(Dir + "/Boundary" + str(core) + ".tcl", "w") + for i in range(tmp.n_points): + if tmp["boundary"][i] != -1: + x,y,z = tmp.points[i] + nodeTag1 = tmp['vtkOriginalPointIds'][i] + nodeTag2 = tmp['boundary'][i] + f.write(f"node {nodeTag2} {str(x)} {str(y)} {str(z)}\n") + f.write(f"equalDOF {nodeTag2} {nodeTag1} 1 2 3\n") + + + + + + + + + + + +# ============================================================================= +# printing information +# ============================================================================= +print(f"Number of regular cores: {reg_num_cores}") +print(f"Number of DRM cores: {DRM_num_cores}") +print(f"Number of PML cores: {PML_num_cores}") +print(f"Number of regular elements: {regular.n_cells} roughly {int(regular.n_cells/reg_num_cores)} each core") +print(f"Number of DRM elements: {DRM.n_cells} roughly {int(DRM.n_cells/DRM_num_cores)} each core") +print(f"Number of PML elements: {PML.n_cells} roughly {int(PML.n_cells/PML_num_cores)} each core") +print(f"Number of total elements: {mesh.n_cells}") +print(f"Number of total points: {mesh.n_points}") +print(f"Number of cores: {max_core-min_core+1}") +print(f"Number of PML nodes: {PML.n_points}") +print(f"Number of regular nodes: {regular.n_points}") +print(f"Number of DRM nodes: {DRM.n_points}") +if AbsorbingElements == "PML": + print(f"Number of MP constraints: {regindicies.size}") + +# calculating number of surface points on the boundaries +eps = 1e-2 +bounds = np.array(mesh.bounds)-np.array([-eps,eps,-eps,eps,-eps,-eps]) +cube = pv.Cube(bounds=bounds) +# points outside the cube +selected = mesh.select_enclosed_points(cube,inside_out = True) +pts = mesh.extract_points( + selected['SelectedPoints'].view(bool), + include_cells=False, +) +print(f"number of sp constriants: {pts.n_points*9}") + + + + + +import h5py +# f = h5py.File('./DRMloadSmall.h5drm', 'r') +f = h5py.File(DRMfile, 'r') +pts = f["DRM_Data"]["xyz"][:] +internal = f["DRM_Data"]["internal"][:] +# %% + + + +pl = pv.Plotter() +# plot the DRM layer with the DRM nodes loaded from the DRM file +pl.add_mesh(DRM,scalars="partitioned",show_edges=True) +pl.add_points(pts, color='r') +pl.export_html(meshplotdir+"/DRM.html") +pl.clear() + + +# plot the regular mesh with the internal nodes loaded from the DRM file +pl.add_mesh(regular,scalars="partitioned",show_edges=True) +pl.export_html(meshplotdir+"/Regular.html") +pl.clear() + +# plot the PML mesh +pl.add_mesh(PML,scalars="partitioned",show_edges=True) +pl.export_html(meshplotdir+"/PML.html") +pl.clear() + + +# plot the total mesh +pl.add_mesh(mesh,scalars="partitioned",show_edges=True) +pl.export_html(meshplotdir+"/Total_partitioned.html") +pl.clear() + +# plot the total mesh with domain scalars +pl.add_mesh(mesh,scalars="Domain",show_edges=True) +pl.export_html(meshplotdir+"/Total_domain.html") +pl.clear() + + +if AbsorbingElements == "ASDA": + pl.add_mesh(mesh,scalars="ASDA_type",show_edges=True,cmap="tab20") + pl.export_html(meshplotdir+"/ASDA_total.html") + pl.clear() + + # filter the mesh with domain pml + indices = mesh['Domain'] == info["PMLDomain"] + grid = mesh.extract_cells(indices) + pl.add_mesh(grid,scalars="ASDA_type",show_edges=True,cmap="tab20b") + pl.export_html(meshplotdir+"/ASDA_PML.html") + pl.clear() + +pl.close() + +# save the mesh +mesh.save(os.path.join(OutputDir,"mesh.vtk"),binary=True) +# mesh.plot(scalars="partitioned",show_edges=True,show_grid=True,show_axes=True,show_bounds=True) +# %% +# print number of PML nodes + + +# %% +# chek for each element that the nodes go counter clockwise and first the below surface nodes and then the above surface nodes +# def iscounterclockwise(points): +# # points = np.array(points) +# points = points - points.mean(axis=0) +# signed_area = 0 +# for j in range(points.shape[0]): +# p1 = points[j,:] +# p2 = points[(j+1)%points.shape[0],:] +# x1,y1,_ = p1 +# x2,y2,_ = p2 +# signed_area += (x1*y2 - x2*y1) + +# return signed_area/2.0 + +# for i in range(mesh.n_cells): +# cell = mesh.get_cell(i) +# points = mesh.points[cell.point_ids] + +# # first surface nodes +# S1 = points[:4,:] +# S2 = points[4:,:] + +# # check that the S2 is above S1 +# if S1[:,2].mean() > S2[:,2].mean(): +# print(f"Element {i} S1 is above S2") +# # else : +# # print(f"Element {i} S2 is above S1") + +# # calcuting signed area of the S1 +# S1 = iscounterclockwise(S1) +# S2 = iscounterclockwise(S2) +# if S1 < 0: +# print(f"Element {i} S1 is not counter clockwise") +# if S2 < 0: +# print(f"Element {i} S2 is not counter clockwise") + + + + + + diff --git a/modules/tools/DRM/MeshGenerator/models/CMakeLists.txt b/modules/tools/DRM/MeshGenerator/models/CMakeLists.txt new file mode 100644 index 000000000..a2d8da276 --- /dev/null +++ b/modules/tools/DRM/MeshGenerator/models/CMakeLists.txt @@ -0,0 +1,2 @@ +simcenter_add_file(NAME Soil_with_structure.tcl) + diff --git a/modules/tools/DRM/MeshGenerator/models/Soil_with_structure.tcl b/modules/tools/DRM/MeshGenerator/models/Soil_with_structure.tcl new file mode 100644 index 000000000..ca6adaef5 --- /dev/null +++ b/modules/tools/DRM/MeshGenerator/models/Soil_with_structure.tcl @@ -0,0 +1,768 @@ +# ===================================================================== # +# 3D test model for the pml element modeling the plane strain field # +# University of Washington, Department of Civil and Environmental Eng # +# Geotechnical Eng Group, A. Pakzad, P. Arduino - Jun 2023 # +# Basic units are m, kg, s (SI units) # +# ===================================================================== # + +# ============================================================================ +# define helper functions +# ============================================================================ +proc writeNodesinFile {filename nodes} { + # check if the file exists + if {[file exists $filename] == 1} { file delete $filename } + set f [open $filename "w"] + foreach node $nodes { + puts $f "$node [nodeCoord $node]" + } + close $f +} + +proc writeElesinFile {filename eles} { + # check if the file exists + if {[file exists $filename] == 1} { file delete $filename } + set f [open $filename "w"] + foreach ele $eles { + puts $f "$ele [eleNodes $ele]" + } + close $f +} + + +proc getMaxNodeTag {np pid} { + set maxNodeTag -1 + foreach node [getNodeTags] {if {$node > $maxNodeTag} {set maxNodeTag $node}} + barrier + # puts "maxNodeTag: $maxNodeTag pid: $pid np: $np" + if {$pid == 0} {set maxNodeTaglist {}} + if {$pid == 0} { + for {set i 0} {$i < $np} {incr i} {lappend maxNodeTaglist -1} + # set the first element of the list to the maxNodeTag + set maxNodeTaglist [lreplace $maxNodeTaglist 0 0 $maxNodeTag] + } + barrier + if {$pid > 0} { + send -pid 0 $maxNodeTag + } else { + for {set i 1} {$i < $np} {incr i} { + recv -pid $i maxNodeTag + set maxNodeTaglist [lreplace $maxNodeTaglist $i $i $maxNodeTag] + } + } + barrier + if {$pid == 0} { + set maxNodeTag -1 + foreach node $maxNodeTaglist {if {$node > $maxNodeTag} {set maxNodeTag $node}} + # puts "maximum: $maxNodeTag" + } + + barrier + # return the maxNodeTag + if {$pid == 0} { + for {set i 1} {$i < $np} {incr i} { + send -pid $i $maxNodeTag + } + } else { + recv -pid 0 maxNodeTag + } + barrier + if {$maxNodeTag == -1} {set maxNodeTag 0} + return $maxNodeTag +} + + +proc getMaxEleTag {np pid} { + set maxEleTag -1 + foreach ele [getEleTags] {if {$ele > $maxEleTag} {set maxEleTag $ele}} + barrier + # puts "maxEleTag: $maxEleTag pid: $pid np: $np" + if {$pid == 0} {set maxEleTaglist {}} + if {$pid == 0} { + for {set i 0} {$i < $np} {incr i} {lappend maxEleTaglist -1} + # set the first element of the list to the maxEleTag + set maxEleTaglist [lreplace $maxEleTaglist 0 0 $maxEleTag] + } + barrier + if {$pid > 0} { + send -pid 0 $maxEleTag + } else { + for {set i 1} {$i < $np} {incr i} { + recv -pid $i maxEleTag + set maxEleTaglist [lreplace $maxEleTaglist $i $i $maxEleTag] + } + } + barrier + if {$pid == 0} { + set maxEleTag -1 + foreach ele $maxEleTaglist {if {$ele > $maxEleTag} {set maxEleTag $ele}} + # puts "maximum: $maxEleTag" + } + + barrier + # return the maxEleTag + if {$pid == 0} { + for {set i 1} {$i < $np} {incr i} { + send -pid $i $maxEleTag + } + } else { + recv -pid 0 maxEleTag + } + barrier + if {$maxEleTag == -1} {set maxEleTag 0} + return $maxEleTag +} + +proc addVartoModelInfoFile {fileName varName varValue pid writerpid} { + if {$pid != $writerpid} {return} + set f [open "$fileName" "r"] + set lines [split [read $f] "\n"] + close $f + set f [open "$fileName" "w+"] + set j 1 + foreach line $lines { + set nextline [lindex $lines $j] + if {$line == "\}"} { + puts $f "\t\"$varName\": \"$varValue\"" + puts $f "\}" + break + } + if {$nextline == "\}" && $j > 1} { + puts $f "$line," + } else { + puts $f "$line" + } + incr j + } + close $f +} + + +proc initializModelinfoFile {fileName pid writerpid } { + if {$pid != $writerpid} {return} + if {[file exists $fileName] == 1} { file delete $fileName } + set f [open "$fileName" "w+"] + puts $f "\{" + puts $f "\}" + close $f +} + + + +# ============================================================================ +# get the number of processors +# ============================================================================ + +set pid [getPID] +set np [getNP] +set ModelInfoFile "$outputdir/ModelInfo.json" + +if {$DRM_Location == "designsafe"} { + set ModelInfoFile "ModelInfo.json" + set meshdir "." + set outputdir "../Results" + set DRMFile "../$DRMFile" + +} + + +initializModelinfoFile $ModelInfoFile $pid 0 +addVartoModelInfoFile $ModelInfoFile "numberOfProcessors" $np $pid 0 + + + +# ============================================================================ +# create the structure model +# ============================================================================ +# structure model will be built on the pid=0 processor all +# other processors will be used for the for fundation, soil, DRM, and PML elements +if {$pid==0} { + + if {$HaveStructure == "YES"} { + puts "Creating the structure model" + } else { + puts "Structure model is skipped" + } +} + +if {$HaveStructure == "YES" } { + if {$pid < $structurecores} { + if {$StructureType == "STEEL"} { + source StructresFiles/SteelStructure.tcl + } + if {$StructureType == "CONCRETE"} { + source ConcreteStructure.tcl + } + if {$StructureType == "Custom"} { + source CustomStructure.tcl + } + puts "Structure model is created" + } +} +barrier + +# ============================================================================ +# update the maxNodeTag and maxEleTag +# ============================================================================ +if {$HaveStructure== "YES" } { + set maxNodeTag [getMaxNodeTag $np $pid] + set maxEleTag [getMaxEleTag $np $pid] +} else { + set maxNodeTag 0 + set maxEleTag 0 +} +barrier + +# ============================================================================ +# Setting the maxNodeTag and maxEleTag for the structure model +# ============================================================================ +# this is really important to set the maxNodeTag and maxEleTag for the structure model +set StructureMaxNodeTag $maxNodeTag +set StructureMaxEleTag $maxEleTag +barrier +addVartoModelInfoFile $ModelInfoFile "StructureMaxNodeTag" $StructureMaxNodeTag $pid 0 +addVartoModelInfoFile $ModelInfoFile "StructureMaxEleTag" $StructureMaxEleTag $pid 0 + + +# ============================================================================ +# bulding regular elements +# ============================================================================ +model BasicBuilder -ndm 3 -ndf 3 +# set Vs 200.0 +# set nu 0.25 ;# --- Poisson's Ratio +# set rho 2300.0 ;# --- Density KG/M^3 +set G [expr $Vs*$Vs*$rho] +set E [expr 2*$G*(1+$nu)] +nDMaterial ElasticIsotropic 1 $E $nu $rho; +nDMaterial ElasticIsotropic 2 [expr $E *100] $nu [expr $rho*4.0]; +set SoilmatTag1 "1 0.0 0.0 0.0"; +set FoundationmatTag1 "2 0.0 0.0 -9.81"; + + +# ============================================================================ +# create regular nodes and elements +# ============================================================================ +if {$pid>= $structurecores && $pid < [expr $regcores+$structurecores]} { + puts "StructureMaxNodeTag : $StructureMaxNodeTag\n" + model BasicBuilder -ndm 3 -ndf 3 + eval "source $meshdir/Nodes$pid.tcl" + + set matTag1 "1 0.0 0.0 0.0"; + set elementType "stdBrick" + eval "source $meshdir/Elements$pid.tcl" + + puts "Regular elements are created" + set recordList [getNodeTags] + set elerecordList [getEleTags] + set rayleighalpha 0.0 + set rayleighbeta 0.0 +} + +barrier + +# ============================================================================ +# creating DRM elements +# ============================================================================ +if {$pid >= [expr $regcores +$structurecores] && $pid < [expr $regcores + $drmcores + $structurecores]} { + model BasicBuilder -ndm 3 -ndf 3 + eval "source $meshdir/Nodes$pid.tcl" + set matTag1 "1 0.0 0.0 0.0"; + set elementType "stdBrick" + eval "source $meshdir/Elements$pid.tcl" + + set elelist [getEleTags] + # region 1 -ele $elelist -rayleigh $Damp_alpha $Damp_beta 0.0 0.0 + set rayleighalpha 0.0 + set rayleighbeta 0.0 + puts "DRM elements are created" +} +barrier + +# ============================================================================ +# Adding pile elements +# ============================================================================ +if {$pid == 0} { + if {$HavePiles == "YES"} { + puts "Creating pile elements" + } else { + puts "Pile elements are skipped" + } +} + +if {$HavePiles == "YES"} { + if {$pid == $structurecores} { + model BasicBuilder -ndm 3 -ndf 6 + set pileElements {} + set pileNodes {} + # creating pile elements + + set secTag 1 + set transfTag 1 + set diameter 1. ; # pile diameter (m) + set radius [expr $diameter/2.] + set pi 3.141593 + set Epile 1e10 + set nu 0.3 + set Gpile [expr $Epile/(2*(1+$nu))] + set Area [expr ($diameter**2)*$pi/2.] + set Iy [expr ($diameter**4)*$pi/64.] + set Iz [expr ($diameter**4)*$pi/64.] + set J [expr ($diameter**4)*$pi/32.] + set transfType "PDelta"; # PDelta, Linear, Corotational + + # section Elastic $secTag $E $A $Iz $Iy $G $J + section Elastic $secTag $Epile $Area $Iz $Iy $Gpile $J + + set numpiles [llength $pilelist] + # puts "Number of piles: $numpiles" + set j 0 + foreach pile $pilelist { + set xbottom [dict get $pile xbottom] + set ybottom [dict get $pile ybottom] + set zbottom [dict get $pile zbottom] + set xtop [dict get $pile xtop] + set ytop [dict get $pile ytop] + set ztop [dict get $pile ztop] + set pilenumelements [dict get $pile numberofElements] + set nPeri 8 + set nLong 8 + set numIntgrPts 5 + + + # creating pile nodes + set pilenumnodes [expr $pilenumelements + 1] + for {set i 1} {$i <= $pilenumnodes} {incr i} { + set xcoord [expr $xbottom + ($xtop-$xbottom)*($i-1)/$pilenumelements] + set ycoord [expr $ybottom + ($ytop-$ybottom)*($i-1)/$pilenumelements] + set zcoord [expr $zbottom + ($ztop-$zbottom)*($i-1)/$pilenumelements] + set Nodeid [expr $soilfoundation_num_nodes + $maxNodeTag + $j*$pilenumnodes + $i] + node $Nodeid $xcoord $ycoord $zcoord + lappend pileNodes $Nodeid + } + + set P1x $xtop + set P1y $ytop + set P1z $ztop + + + + # normal vector to the pile + set normalX [expr $xtop - $xbottom] + set normalY [expr $ytop - $ybottom] + set normalZ [expr $ztop - $zbottom] + set norm [expr sqrt($normalX**2 + $normalY**2 + $normalZ**2)] + set normalX [expr $normalX/$norm] + set normalY [expr $normalY/$norm] + set normalZ [expr $normalZ/$norm] + + # ax + by + cz = d + set d [expr $normalX*$xtop + $normalY*$ytop + $normalZ*$ztop] + + # find another point on the plane + set P2x [expr $xtop + 1] + set P2y [expr $ytop + 1] + set P2z [expr ($d - $normalX*$P2x - $normalY*$P2y)/$normalZ] + + set VecX_x [expr $P2x - $P1x] + set VecX_y [expr $P2y - $P1y] + set VecX_z [expr $P2z - $P1z] + set norm [expr sqrt($VecX_x**2 + $VecX_y**2 + $VecX_z**2)] + set VecX_x [expr $VecX_x/$norm] + set VecX_y [expr $VecX_y/$norm] + set VecX_z [expr $VecX_z/$norm] + + + + + + set transfTag [expr $soilfoundation_num_cells + $maxEleTag + $j*$pilenumelements + 1] + eval "geomTransf $transfType $transfTag $VecX_x $VecX_y $VecX_z" + # creating pile elements + for {set i 1} {$i < $pilenumnodes} {incr i} { + set node1 [expr $soilfoundation_num_nodes + $maxNodeTag + $j*$pilenumnodes + $i] + set node2 [expr $soilfoundation_num_nodes + $maxNodeTag + $j*$pilenumnodes + $i + 1] + set eleTag [expr $soilfoundation_num_cells + $maxEleTag + $j*$pilenumelements + $i] + element dispBeamColumn $eleTag $node1 $node2 $numIntgrPts $secTag $transfTag + lappend pileElements $eleTag + } + + # creating soil-pile interface + set interfaceElemsFile "$outputdir/PileInterfaceElements$j.dat" + set connecFile "$outputdir/PileInterfaceConnections$j.dat" + if {[file exists $interfaceElemsFile] == 1} { file delete $interfaceElemsFile } + if {[file exists $connecFile] == 1} { file delete $connecFile } + + set interfaceElems {} + set first_beam_ele [expr $soilfoundation_num_cells + $maxEleTag + $j*$pilenumelements + 1] + set last_beam_ele [expr $soilfoundation_num_cells + $maxEleTag + $j*$pilenumelements + $pilenumelements] + set interfaceElems [generateInterfacePoints -beamEleRange $first_beam_ele $last_beam_ele -gPenalty -shape circle -nP $nPeri -nL $nLong -crdTransf $transfTag -radius $radius -penaltyParam 1.0e12 -file $interfaceElemsFile -connectivity $connecFile ] + puts "interfaceElems: $interfaceElems" + + + + incr j + } + + + + + + + + # puts "writing pile nodes and elements to a file" + # write pile nodes to a file + set pileNodesFile "$outputdir/PileNodes.dat" + writeNodesinFile $pileNodesFile $pileNodes + + + # write beam elements to a file + set beamElemsFile "$outputdir/PileElements.dat" + writeElesinFile $beamElemsFile $pileElements + + puts "Pile elements are created" + } +} +barrier +set maxNodeTag [getMaxNodeTag $np $pid] +set maxEleTag [getMaxEleTag $np $pid] + +# ============================================================================ +# Creating column and foundation elements +# ============================================================================ +if {$pid == 0} { + if {$HaveStructure == "YES" } { + if {$StructureType == "STEEL" || $StructureType == "CONCRETE"} { + puts "Creating base columns and foundation elements" + } else { + puts "Base columns and foundation elements are skipped" + } + } +} +if {$HaveStructure == "YES" } { + if {$pid == $structurecores && $structurecores > 0} { + + model BasicBuilder -ndm 3 -ndf 6 + set StrucFoundConnecElems {} + set StrucFoundConnecNodes {} + + if {$StructureType == "STEEL" || $StructureType == "CONCRETE"} {set BaseColumnsFile "$outputdir/BaseColumnsNodes.dat"} + + # open BaseColumnsNodes.dat file + set f [open $BaseColumnsFile "r"] + set lines [split [read $f] "\n"] + close $f + + set EmbeddingDepth 0.75 + set j 0 + foreach line $lines { + # check if the line is empty + if {[string length $line] == 0} {continue} + incr j + set nodeTag1 [lindex $line 0] + set x [lindex $line 1] + set y [lindex $line 2] + set z [lindex $line 3] + node $nodeTag1 $x $y $z + set nodeTag2 [expr $maxNodeTag + $j] + set x [expr $x] + set y [expr $y] + set z [expr $z - $EmbeddingDepth] + node $nodeTag2 $x $y $z + set numIntgrPts 3 + lappend StrucFoundConnecNodes $nodeTag1 + lappend StrucFoundConnecNodes $nodeTag2 + + set nPeri 5 + set nLong 3 + set secTag [expr $maxNodeTag + $j] + set transfTag [expr $maxNodeTag + $j] + set diameter 0.25 ; + set radius [expr $diameter/2.] + set pi 3.141593 + set Epile 1e10 + set nu 0.3 + set Gpile [expr $Epile/(2*(1+$nu))] + set Area [expr ($diameter**2)*$pi/2.] + set Iy [expr ($diameter**4)*$pi/64.] + set Iz [expr ($diameter**4)*$pi/64.] + set J [expr ($diameter**4)*$pi/32.] + set transfType "Linear"; # PDelta, Linear, Corotational + + + section Elastic $secTag $Epile $Area $Iz $Iy $Gpile $J + geomTransf $transfType $transfTag 1 0 0 + + + set eleTag [expr $maxEleTag + $j] + element dispBeamColumn $eleTag $nodeTag2 $nodeTag1 $numIntgrPts $secTag $transfTag + lappend StrucFoundConnecElems $eleTag + + # creating soil-pile interface + set num [expr $j-1] + set interfaceElemsFile "$outputdir/StructureFoundationInterfaceElements$num.dat" + set connecFile "$outputdir/StructureFoundationInterfaceConnections$num.dat" + if {[file exists $interfaceElemsFile] == 1} { file delete $interfaceElemsFile } + if {[file exists $connecFile] == 1} { file delete $connecFile } + + set interfaceElems {} + set interfaceElems [generateInterfacePoints -beamEleRange $eleTag $eleTag -gPenalty -shape circle -nP $nPeri -nL $nLong -crdTransf $transfTag -radius $radius -penaltyParam 1.0e12 -file $interfaceElemsFile -connectivity $connecFile ] + set maxEleTag $interfaceElems + } + puts "Base columns and foundation elements are attached" + + # write pile nodes to a file + set StrucFoundConnecNodesFile "$outputdir/StructureFoundationBeamNodes.dat" + writeNodesinFile $StrucFoundConnecNodesFile $StrucFoundConnecNodes + + # write beam elements to a file + set StrucFoundConnecElemsFile "$outputdir/StructureFoundationBeamElements.dat" + writeElesinFile $StrucFoundConnecElemsFile $StrucFoundConnecElems + } +} +barrier +# ============================================================================ +# bulding PML layer +# ============================================================================ +if {$pid == 0} { + if {$HaveAbsorbingElements == "YES"} { + puts "Creating Absorbing elements" + } else { + puts "Absorbing elements are skipped" + } +} +if {$HaveAbsorbingElements == "YES" && $pid >= [expr $structurecores + $regcores + $drmcores] } { + # create PML material + set gamma 0.5 ;# --- Coefficient gamma, newmark gamma = 0.5 + set beta 0.25 ;# --- Coefficient beta, newmark beta = 0.25 + set eta [expr 1.0/12.] ;# --- Coefficient eta, newmark eta = 1/12 + set E $E ;# --- Young's modulus + set nu $nu ;# --- Poisson's Ratio + set rho $rho ;# --- Density + set EleType 6 ;# --- Element type, See line + set PML_L $pmltotalthickness ;# --- Thickness of the PML + set afp 2.0 ;# --- Coefficient m, typically m = 2 + set PML_Rcoef 1.0e-8 ;# --- Coefficient R, typically R = 1e-8 + set RD_half_width_x [expr $llx/2.] ;# --- Halfwidth of the regular domain in + set RD_half_width_y [expr $lly/2.] ;# --- Halfwidth of the regular domain in + set RD_depth [expr $llz/1.] ;# --- Depth of the regular domain + set pi 3.141593 ;# --- pi + + set Damp_alpha 0.0 ;# --- Rayleigh damping coefficient alpha + set Damp_beta 0.0 ;# --- Rayleigh damping coefficient beta + + + if {$AbsorbingElements == "PML"} { + puts "Absorbing elements are PML elements" + model BasicBuilder -ndm 3 -ndf 9; + set AbsorbingmatTag1 "$eta $beta $gamma $E $nu $rho $EleType $PML_L $afp $PML_Rcoef $RD_half_width_x $RD_half_width_y $RD_depth $Damp_alpha $Damp_beta" + # set elementType "PMLVISCOUS" + set elementType "PML" + eval "source $meshdir/Nodes$pid.tcl" + eval "source $meshdir/Elements$pid.tcl" + + # tie pml nodes to the regular nodes + model BasicBuilder -ndm 3 -ndf 3; + eval "source $meshdir/Boundary$pid.tcl" + } + if {$AbsorbingElements == "ASDA"} { + puts "Absorbing elements are ASDA elements" + model BasicBuilder -ndm 3 -ndf 3; + set AbsorbingmatTag1 "$G $nu $rho"; + set elementType "ASDAbsorbingBoundary3D" + # set elementType "PML" + eval "source $meshdir/Nodes$pid.tcl" + eval "source $meshdir/Elements$pid.tcl" + } + if {$AbsorbingElements == "Rayleigh"} { + puts "Absorbing elements are normal elements" + model BasicBuilder -ndm 3 -ndf 3; + set AbsorbingmatTag1 "1 0.0 0.0 0.0"; + set elementType "stdBrick" + eval "source $meshdir/Nodes$pid.tcl" + eval "source $meshdir/Elements$pid.tcl" + set rayleighalpha $Absorbing_rayleigh_alpha + set rayleighbeta $Absorbing_rayleigh_beta + + } +} +barrier + + +# ============================================================================ +# creating fixities +# ============================================================================ +if {$HaveAbsorbingElements == "YES"} { + set totalThickness [expr $numdrmlayers*$drmthicknessx + $numpmllayers*$pmlthicknessx] + if {$AbsorbingElements == "PML"} { + fixX [expr -$lx/2. - $totalThickness] 1 1 1 1 1 1 1 1 1; + fixX [expr $lx/2. + $totalThickness] 1 1 1 1 1 1 1 1 1; + fixY [expr -$ly/2. - $totalThickness] 1 1 1 1 1 1 1 1 1; + fixY [expr $ly/2. + $totalThickness] 1 1 1 1 1 1 1 1 1; + fixZ [expr -$lz/1. - $totalThickness] 1 1 1 1 1 1 1 1 1; + } + if {$AbsorbingElements == "ASDA" || $AbsorbingElements == "Normal"} { + fixX [expr -$lx/2. - $totalThickness] 1 1 1; + fixX [expr $lx/2. + $totalThickness] 1 1 1; + fixY [expr -$ly/2. - $totalThickness] 1 1 1; + fixY [expr $ly/2. + $totalThickness] 1 1 1; + fixZ [expr -$lz/1. - $totalThickness] 1 1 1; + } + puts "Boundary conditions are applied to the absorbing Layer nodes" +} else { + set totalThickness [expr $numdrmlayers*$drmthicknessx] + fixX [expr -$lx/2. - $totalThickness] 1 1 1; + fixX [expr $lx/2. + $totalThickness] 1 1 1; + fixY [expr -$ly/2. - $totalThickness] 1 1 1; + fixY [expr $ly/2. + $totalThickness] 1 1 1; + fixZ [expr -$lz/1. - $totalThickness] 1 1 1; + puts "Boundary conditions are applied to the DRM Layer nodes" +} +barrier + +# ============================================================================ +# printing model information for debugging +# ============================================================================ +# print "modelInfo$pid" + + +# ============================================================================ +# Gravity analysis +# ============================================================================ +domainChange +constraints Plain +numberer ParallelRCM +system Mumps -ICNTL14 400 +test NormDispIncr 1e-4 10 2 +algorithm Linear -factorOnce +# algorithm ModifiedNewton -factoronce +# algorithm Newton +integrator Newmark 0.5 0.25 +analysis Transient + + +barrier +loadConst -time 0.0 +wipeAnalysis +puts "Gravity analysis is done" +# ============================================================================ +# Gravirty recorders +# ============================================================================ +# if {$pid >= $structurecores && $pid < [expr $regcores + $structurecores]} { +# eval "recorder Node -file $outputdir/GravityNodeDisp$pid.out -node $recordList -dof 1 2 3 disp" +# eval "recorder Node -file $outputdir/GravityNodeAccl$pid.out -node $recordList -dof 1 2 3 accel" +# eval "recorder Node -file $outputdir/GravityNodeVelo$pid.out -node $recordList -dof 1 2 3 vel" +# eval "recorder Element -file $outputdir/GravityElementStress$pid.out -ele $elerecordList stresses" +# eval "recorder Element -file $outputdir/GravityElementStrain$pid.out -ele $elerecordList strains" + +# # print recordlist and elerecordlist to a file +# set f [open "$outputdir/nodeOuputTags$pid.out" "w+"] +# puts $f "$recordList" +# close $f +# set f [open "$outputdir/eleOuputTags$pid.out" "w+"] +# puts $f "$elerecordList" +# close $f +# } +# record +# record +# barrier +# remove recorders + +# ============================================================================ +# Post gravity settings +# ============================================================================ +# ASDA Absorbing layer settings +if {$pid >= [expr $regcores +$structurecores] && $pid < [expr $regcores + $drmcores + $structurecores]} { + if {$pid >= [expr $structurecores + $regcores + $drmcores] && $AbsorbingElements == "ASDA" } { + set abs_elements [getEleTags] + eval "setParameter -val 1 -ele $abs_elements stage" + } +} + +# # set the initial displacements to zero +# foreach node [getNodeTags] { +# foreach dof {1 2 3} { +# setNodeAccel $node $dof 0.0 -commit +# setNodeVel $node $dof 0.0 -commit +# setNodeDisp $node $dof 0.0 -commit +# } +# } +# ============================================================================ +# loading +# ============================================================================ + +# if {$pid < $regcores} { +if {$pid>=$regcores && $pid < [expr $regcores + $drmcores] } { + set prop1 "$DRM_factor $crd_scale $distance_tolerance $do_coordinate_transformation"; # factor crd_scale distance_tolerance do_coordinate_transformation + set prop2 "$T00 $T01 $T02"; # T00, T01, T02 + set prop3 "$T10 $T11 $T12"; # T10, T11, T12 + set prop4 "$T20 $T21 $T22"; # T20, T21, T22 + set prop5 "$originX $originY $originZ"; # x00, x01, x02 + + set DRMtag 1 + eval "pattern H5DRM $DRMtag $DRMFile $prop1 $prop2 $prop3 $prop4 $prop5" +} + +# # ============================================================================ +# # recorders +# # ============================================================================ +# find the nodetag with 0.0 0.0 0.0 +if {$pid < $regcores} { + set foundflag "No" + foreach node [getNodeTags] { + set coords [nodeCoord $node] + set diff1 [expr abs([lindex $coords 0])] + set diff2 [expr abs([lindex $coords 1])] + set diff3 [expr abs([lindex $coords 2])] + if { $diff1 < 1e-6 && $diff2 < 1e-6 && $diff3 < 1e-6 } { + set recordNode $node + set foundflag "Yes" + break + } + } + if {$foundflag == "Yes"} { + recorder Node -file $outputdir/Disp$pid.out -time -dT $Analysis_record_dt -node $recordNode -dof 1 2 3 disp + recorder Node -file $outputdir/Accel$pid.out -time -dT $Analysis_record_dt -node $recordNode -dof 1 2 3 accel + recorder Node -file $outputdir/Vel$pid.out -time -dT $Analysis_record_dt -node $recordNode -dof 1 2 3 vel + } +} +record + +# ============================================================================ +# Analysis +# ============================================================================ +domainChange + +# constraints Plain +constraints Plain +numberer ParallelRCM +system Mumps -ICNTL14 400 +test NormDispIncr 1e-4 10 2 +algorithm Linear -factorOnce +integrator Newmark 0.5 0.25 +analysis Transient + +if {$AbsorbingElements == "Rayleigh"} { + puts "Rayleigh damping is applied" + rayleigh $rayleighalpha $rayleighbeta 0.0 0.0 +} +set startTime [clock milliseconds] + +set dT $Analysis_dt +while { [getTime] < $Analysis_duration } { + if {$pid ==0 } {puts "Time: [getTime]";} + analyze 1 $dT +} + +set endTime [clock milliseconds] +set elapsedTime [expr {$endTime - $startTime}] +puts "Elapsed time: [expr $elapsedTime/1000.] seconds in $pid" + + + +# if {$pid == 0} { +# puts "natural frequency of the pile: $f Hz (assuming cantilever beam)" +# puts "wavelegth: [expr $Vs/$f] m" +# puts "Vs: $Vs" +# } +wipeAnalysis +remove recorders +remove loadPattern 2 \ No newline at end of file diff --git a/modules/tools/DRM/StructresFiles/CMakeLists.txt b/modules/tools/DRM/StructresFiles/CMakeLists.txt new file mode 100644 index 000000000..4d428ca1f --- /dev/null +++ b/modules/tools/DRM/StructresFiles/CMakeLists.txt @@ -0,0 +1 @@ +simcenter_add_file(NAME SteelStructure.tcl) diff --git a/modules/tools/DRM/StructresFiles/SteelStructure.tcl b/modules/tools/DRM/StructresFiles/SteelStructure.tcl new file mode 100644 index 000000000..70dd10e9c --- /dev/null +++ b/modules/tools/DRM/StructresFiles/SteelStructure.tcl @@ -0,0 +1,485 @@ +# -------------------------------------------------------------------------------------------------- +# Example 8. 3D Steel W-section frame +# Silvia Mazzoni & Frank McKenna, 2006 +# nonlinearBeamColumn element, inelastic fiber section +# + + + +# SET UP ---------------------------------------------------------------------------- +model BasicBuilder -ndm 3 -ndf 6; # Define the model builder, ndm=#dimension, ndf=#dofs + + +# ================================================================================== +# units +# ================================================================================== +set meter 1.; # define basic units -- input units +set Newton 1.; # define basic units -- input units +set kN 1000.; # define basic units -- input units +set MN 1.e6; # define basic units -- input units +set mm 0.001; # define basic units -- input units +set cm 0.01; # define basic units -- input units +set Pa 1.; # define basic units -- input units +set kPa 1.e3; # define basic units -- input units +set MPa 1.e6; # define basic units -- input units +set GPa 1.e9; # define basic units -- input units +set in 0.0254; # define basic units -- output units +set kip 4448.2216153; # define basic units -- output units +set sec 1.; # define basic units -- output units +set LunitTXT "meter"; # define basic-unit text for output +set FunitTXT "Newton"; # define basic-unit text for output +set TunitTXT "sec"; # define basic-unit text for output +set ft [expr 12.*$in]; # define engineering units +set ksi [expr $kip/pow($in,2)]; +set psi [expr $ksi/1000.]; +set lbf [expr $psi*$in*$in]; # pounds force +set pcf [expr $lbf/pow($ft,3)]; # pounds per cubic foot +set psf [expr $lbf/pow($ft,3)]; # pounds per square foot +set in2 [expr $in*$in]; # inch^2 +set in4 [expr $in*$in*$in*$in]; # inch^4 +set cm [expr $in/2.54]; # centimeter, needed for displacement input in MultipleSupport excitation +set PI [expr 2*asin(1.0)]; # define constants +set g [expr 32.2*$ft/pow($sec,2)]; # gravitational acceleration +set Ubig 1.e10; # a really large number +set Usmall [expr 1/$Ubig]; # a really small number + +# ================================================================================== +# fiber section W-section +# ================================================================================== +proc Wsection { secID matID d bf tf tw nfdw nftw nfbf nftf} { + # ################################################################### + # Wsection $secID $matID $d $bf $tf $tw $nfdw $nftw $nfbf $nftf + # ################################################################### + # create a standard W section given the nominal section properties + # written: Remo M. de Souza + # date: 06/99 + # modified: 08/99 (according to the new general modelbuilder) + # input parameters + # secID - section ID number + # matID - material ID number + # d = nominal depth + # tw = web thickness + # bf = flange width + # tf = flange thickness + # nfdw = number of fibers along web depth + # nftw = number of fibers along web thickness + # nfbf = number of fibers along flange width + # nftf = number of fibers along flange thickness + + set dw [expr $d - 2 * $tf] + set y1 [expr -$d/2] + set y2 [expr -$dw/2] + set y3 [expr $dw/2] + set y4 [expr $d/2] + + set z1 [expr -$bf/2] + set z2 [expr -$tw/2] + set z3 [expr $tw/2] + set z4 [expr $bf/2] + + section fiberSec $secID { + # nfIJ nfJK yI zI yJ zJ yK zK yL zL + patch quadr $matID $nfbf $nftf $y1 $z4 $y1 $z1 $y2 $z1 $y2 $z4 + patch quadr $matID $nftw $nfdw $y2 $z3 $y2 $z2 $y3 $z2 $y3 $z3 + patch quadr $matID $nfbf $nftf $y3 $z4 $y3 $z1 $y4 $z1 $y4 $z4 + } +} + + + +# ================================================================================== +# define Geometry +# ================================================================================== +set BaseColumnsIDs {}; # list of nodes at base of columns +set StructureNodes {}; # list of all nodes in the structure +# define NODAL COORDINATES +set Dlevel 10000; # numbering increment for new-level nodes +set Dframe 100; # numbering increment for new-frame nodes +set NFrame [expr $NBayZ + 1]; # number of frames + +for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} { + set Y [expr ($frame-1)*$LGird]; + set Y [expr $Y + $StartNodeY] + for {set level 1} {$level <=[expr $NStory+1]} {incr level 1} { + set Z [expr ($level-1)*$LCol]; + set Z [expr $Z + $StartNodeZ] + for {set pier 1} {$pier <= [expr $NBay+1]} {incr pier 1} { + set X [expr ($pier-1)*$LBeam]; + set X [expr $X + $StartNodeX] + + set nodeID [expr $frame*$Dframe+$level*$Dlevel+$pier] + node $nodeID $X $Y $Z; # actually define node + lappend StructureNodes $nodeID + if {$level == 1} {lappend BaseColumnsIDs $nodeID} + } + } +} + + +# ================================================================================== +# rigid diaphragm nodes +# ================================================================================== +set RigidDiaphragm OFF ; # options: ON, OFF. specify this before the analysis parameters are set the constraints are handled differently. +if {$RigidDiaphragm == "ON"} { + set Xa [expr ($NBay*$LBeam)/2]; # mid-span coordinate for rigid diaphragm + set Ya [expr ($NFrame-1)*$LGird/2]; + # move the structure to the start node + set Xa [expr $Xa + $StartNodeX] + set Ya [expr $Ya + $StartNodeY] + + set iMasterNode "" + for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} { + set Z [expr ($level-1)*$LCol]; + set Z [expr $Z + $StartNodeZ] + + + # rigid-diaphragm nodes in center of each diaphram + set MasterNodeID [expr 9900+$level] + node $MasterNodeID $Xa $Ya $Z; # master nodes for rigid diaphragm + fix $MasterNodeID 0 1 0 1 0 1; # constrain other dofs that don't belong to rigid diaphragm control + lappend iMasterNode $MasterNodeID + set perpDirn 3; # perpendicular to plane of rigid diaphragm + for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} { + for {set pier 1} {$pier <= [expr $NBay+1]} {incr pier 1} { + set nodeID [expr $level*$Dlevel+$frame*$Dframe+$pier] + rigidDiaphragm $perpDirn $MasterNodeID $nodeID; # define Rigid Diaphram, + # puts [nodeCoord $nodeID] + # puts [nodeCoord $MasterNodeID] + } + } + lappend StructureNodes $MasterNodeID + } +} + + +# ================================================================================== +# calculated MODEL PARAMETERS, particular to this model +# ================================================================================== +# Set up parameters that are particular to the model for displacement control +set IDctrlNode [expr int((($NStory+1)*$Dlevel+$NFrame*$Dframe)+1)]; # node where displacement is read for displacement control +set IDctrlDOF 1; # degree of freedom of displacement read for displacement control +set LBuilding [expr $NStory*$LCol]; # total building height + + + +# ================================================================================== +# Define SECTIONS ------------------------------------------------------------- +# ================================================================================== + +# define section tags: +set ColSecTag 1 +set BeamSecTag 2 +set GirdSecTag 3 +set ColSecTagFiber 4 +set BeamSecTagFiber 5 +set GirdSecTagFiber 6 +set SecTagTorsion 70 + + +if {$SectionType == "Elastic"} { + # material properties: + set Es [expr 29000*$ksi]; # Steel Young's Modulus + set nu 0.3; # Poisson's ratio + set Gs [expr $Es/2./[expr 1+$nu]]; # Torsional stiffness Modulus + set J $Ubig; # set large torsional stiffness + # column sections: W27x114 + set AgCol [expr 33.5*pow($in,2)]; # cross-sectional area + set IzCol [expr 4090.*pow($in,4)]; # moment of Inertia + set IyCol [expr 159.*pow($in,4)]; # moment of Inertia + # beam sections: W24x94 + set AgBeam [expr 27.7*pow($in,2)]; # cross-sectional area + set IzBeam [expr 2700.*pow($in,4)]; # moment of Inertia + set IyBeam [expr 109.*pow($in,4)]; # moment of Inertia + # girder sections: W24x94 + set AgGird [expr 27.7*pow($in,2)]; # cross-sectional area + set IzGird [expr 2700.*pow($in,4)]; # moment of Inertia + set IyGird [expr 109.*pow($in,4)]; # moment of Inertia + + section Elastic $ColSecTag $Es $AgCol $IzCol $IyCol $Gs $J + section Elastic $BeamSecTag $Es $AgBeam $IzBeam $IyBeam $Gs $J + section Elastic $GirdSecTag $Es $AgGird $IzGird $IyGird $Gs $J + + set matIDhard 1; # material numbers for recorder (this stressstrain recorder will be blank, as this is an elastic section) + +} elseif {$SectionType == "FiberSection"} { + # define MATERIAL properties + set Fy [expr 60.0*$ksi] + set Es [expr 29000*$ksi]; # Steel Young's Modulus + set nu 0.3; + set Gs [expr $Es/2./[expr 1+$nu]]; # Torsional stiffness Modulus + set Hiso 0 + set Hkin 1000 + set matIDhard 1 + uniaxialMaterial Hardening $matIDhard $Es $Fy $Hiso $Hkin + + # ELEMENT properties + # Structural-Steel W-section properties + # column sections: W27x114 + set d [expr 27.29*$in]; # depth + set bf [expr 10.07*$in]; # flange width + set tf [expr 0.93*$in]; # flange thickness + set tw [expr 0.57*$in]; # web thickness + set nfdw 16; # number of fibers along dw + set nftw 2; # number of fibers along tw + set nfbf 16; # number of fibers along bf + set nftf 4; # number of fibers along tf + Wsection $ColSecTagFiber $matIDhard $d $bf $tf $tw $nfdw $nftw $nfbf $nftf + # beam sections: W24x94 + set d [expr 24.31*$in]; # depth + set bf [expr 9.065*$in]; # flange width + set tf [expr 0.875*$in]; # flange thickness + set tw [expr 0.515*$in]; # web thickness + set nfdw 16; # number of fibers along dw + set nftw 2; # number of fibers along tw + set nfbf 16; # number of fibers along bf + set nftf 4; # number of fibers along tf + Wsection $BeamSecTagFiber $matIDhard $d $bf $tf $tw $nfdw $nftw $nfbf $nftf + # girder sections: W24x94 + set d [expr 24.31*$in]; # depth + set bf [expr 9.065*$in]; # flange width + set tf [expr 0.875*$in]; # flange thickness + set tw [expr 0.515*$in]; # web thickness + set nfdw 16; # number of fibers along dw + set nftw 2; # number of fibers along tw + set nfbf 16; # number of fibers along bf + set nftf 4; # number of fibers along tf + Wsection $GirdSecTagFiber $matIDhard $d $bf $tf $tw $nfdw $nftw $nfbf $nftf + + # assign torsional Stiffness for 3D Model + uniaxialMaterial Elastic $SecTagTorsion $Ubig + section Aggregator $ColSecTag $SecTagTorsion T -section $ColSecTagFiber + section Aggregator $BeamSecTag $SecTagTorsion T -section $BeamSecTagFiber + section Aggregator $GirdSecTag $SecTagTorsion T -section $GirdSecTagFiber +} else { + puts "No section has been defined" + return -1 +} + +set QdlCol [expr 114*$lbf/$ft]; # W-section weight per length +set QBeam [expr 94*$lbf/$ft]; # W-section weight per length +set QGird [expr 94*$lbf/$ft]; # W-section weight per length + +# ================================================================================== +# define ELEMENTS +# ================================================================================== +# set up geometric transformations of element +# separate columns and beams, in case of P-Delta analysis for columns +set IDColTransf 1; # all columns +set IDBeamTransf 2; # all beams +set IDGirdTransf 3; # all girds +set ColTransfType Linear ; # options for columns: Linear PDelta Corotational +geomTransf $ColTransfType $IDColTransf 0 1 0; # orientation of column stiffness affects bidirectional response. +geomTransf Linear $IDBeamTransf 0 -1 0; +geomTransf Linear $IDGirdTransf 1 0 0; + + +# Define Beam-Column Elements +set StructureElements {} +set numIntgrPts 5; # number of Gauss integration points for nonlinear curvature distribution +# columns +set N0col [expr 10000-1]; # column element numbers +set level 0 +for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} { + for {set level 1} {$level <=$NStory} {incr level 1} { + for {set pier 1} {$pier <= [expr $NBay+1]} {incr pier 1} { + set elemID [expr $N0col +$level*$Dlevel + $frame*$Dframe+$pier] + set nodeI [expr $level*$Dlevel + $frame*$Dframe+$pier] + set nodeJ [expr ($level+1)*$Dlevel + $frame*$Dframe+$pier] + element nonlinearBeamColumn $elemID $nodeI $nodeJ $numIntgrPts $ColSecTag $IDColTransf; # columns + lappend StructureElements $elemID + } + } +} + +# beams -- parallel to X-axis +set N0beam 1000000; # beam element numbers +for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} { + for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} { + for {set bay 1} {$bay <= $NBay} {incr bay 1} { + set elemID [expr $N0beam +$level*$Dlevel + $frame*$Dframe+ $bay] + set nodeI [expr $level*$Dlevel + $frame*$Dframe+ $bay] + set nodeJ [expr $level*$Dlevel + $frame*$Dframe+ $bay+1] + element nonlinearBeamColumn $elemID $nodeI $nodeJ $numIntgrPts $BeamSecTag $IDBeamTransf; # beams + lappend StructureElements $elemID + } + } +} + +# girders -- parallel to Z-axis +set N0gird 2000000; # gird element numbers +for {set frame 1} {$frame <=[expr $NFrame-1]} {incr frame 1} { + for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} { + for {set bay 1} {$bay <= $NBay+1} {incr bay 1} { + set elemID [expr $N0gird + $level*$Dlevel +$frame*$Dframe+ $bay] + set nodeI [expr $level*$Dlevel + $frame*$Dframe+ $bay] + set nodeJ [expr $level*$Dlevel + ($frame+1)*$Dframe+ $bay] + element nonlinearBeamColumn $elemID $nodeI $nodeJ $numIntgrPts $GirdSecTag $IDGirdTransf; # Girds + lappend StructureElements $elemID + } + } +} + +# -------------------------------------------------------------------------------------------------------------------------------- +# Define GRAVITY LOADS, weight and masses +# calculate dead load of frame, assume this to be an internal frame (do LL in a similar manner) +# calculate distributed weight along the beam length +set GammaConcrete [expr 150*$pcf]; # Reinforced-Concrete floor slabs +set Tslab [expr 6*$in]; # 6-inch slab +set Lslab [expr 2*$LGird/2]; # slab extends a distance of $LGird/2 in/out of plane +set DLfactor 2.0; # scale dead load up a little +set Qslab [expr $GammaConcrete*$Tslab*$Lslab*$DLfactor]; +set QdlBeam [expr $Qslab + $QBeam]; # dead load distributed along beam (one-way slab) +set QdlGird $QGird; # dead load distributed along girder +set WeightCol [expr $QdlCol*$LCol]; # total Column weight +set WeightBeam [expr $QdlBeam*$LBeam]; # total Beam weight +set WeightGird [expr $QdlGird*$LGird]; # total Beam weight + +# assign masses to the nodes that the columns are connected to +# each connection takes the mass of 1/2 of each element framing into it (mass=weight/$g) +set iFloorWeight "" +set WeightTotal 0.0 +set sumWiHi 0.0; # sum of storey weight times height, for lateral-load distribution + +for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} { + if {$frame == 1 || $frame == $NFrame} { + set GirdWeightFact 1; # 1x1/2girder on exterior frames + } else { + set GirdWeightFact 2; # 2x1/2girder on interior frames + } + for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} { ; + set FloorWeight 0.0 + if {$level == [expr $NStory+1]} { + set ColWeightFact 1; # one column in top story + } else { + set ColWeightFact 2; # two columns elsewhere + } + for {set pier 1} {$pier <= [expr $NBay+1]} {incr pier 1} {; + if {$pier == 1 || $pier == [expr $NBay+1]} { + set BeamWeightFact 1; # one beam at exterior nodes + } else {; + set BeamWeightFact 2; # two beams elewhere + } + set WeightNode [expr $ColWeightFact*$WeightCol/2 + $BeamWeightFact*$WeightBeam/2 + $GirdWeightFact*$WeightGird/2] + set MassNode [expr $WeightNode/$g]; + set nodeID [expr $level*$Dlevel+$frame*$Dframe+$pier] + mass $nodeID $MassNode $MassNode 0. 0. 0. 0.; # define mass + set FloorWeight [expr $FloorWeight+$WeightNode]; + } + lappend iFloorWeight $FloorWeight + set WeightTotal [expr $WeightTotal+ $FloorWeight] + set sumWiHi [expr $sumWiHi+$FloorWeight*($level-1)*$LCol]; # sum of storey weight times height, for lateral-load distribution + } +} +set MassTotal [expr $WeightTotal/$g]; # total mass + + +# # Define RECORDERS ------------------------------------------------------------- +# set FreeNodeID [expr $NFrame*$Dframe+($NStory+1)*$Dlevel+($NBay+1)]; # ID: free node +# set SupportNodeFirst [lindex $iSupportNode 0]; # ID: first support node +# set SupportNodeLast [lindex $iSupportNode [expr [llength $iSupportNode]-1]]; # ID: last support node +# set FirstColumn [expr $N0col + 1*$Dframe+1*$Dlevel +1]; # ID: first column +# recorder Node -file $dataDir/DFree.out -time -node $FreeNodeID -dof 1 2 3 disp; # displacements of free node +# recorder Node -file $dataDir/DBase.out -time -nodeRange $SupportNodeFirst $SupportNodeLast -dof 1 2 3 disp; # displacements of support nodes +# recorder Node -file $dataDir/RBase.out -time -nodeRange $SupportNodeFirst $SupportNodeLast -dof 1 2 3 reaction; # support reaction +# recorder Drift -file $dataDir/DrNode.out -time -iNode $SupportNodeFirst -jNode $FreeNodeID -dof 1 -perpDirn 2; # lateral drift +# recorder Element -file $dataDir/Fel1.out -time -ele $FirstColumn localForce; # element forces in local coordinates +# recorder Element -file $dataDir/ForceEle1sec1.out -time -ele $FirstColumn section 1 force; # section forces, axial and moment, node i +# recorder Element -file $dataDir/DefoEle1sec1.out -time -ele $FirstColumn section 1 deformation; # section deformations, axial and curvature, node i +# recorder Element -file $dataDir/ForceEle1sec$numIntgrPts.out -time -ele $FirstColumn section $numIntgrPts force; # section forces, axial and moment, node j +# recorder Element -file $dataDir/DefoEle1sec$numIntgrPts.out -time -ele $FirstColumn section $numIntgrPts deformation; # section deformations, axial and curvature, node j +# set yFiber [expr 0.]; # fiber location for stress-strain recorder, local coords +# set zFiber [expr 0.]; # fiber location for stress-strain recorder, local coords +# recorder Element -file $dataDir/SSEle1sec1.out -time -ele $FirstColumn section $numIntgrPts fiber $yFiber $zFiber stressStrain; # steel fiber stress-strain, node i + + + +# GRAVITY ------------------------------------------------------------- +# define GRAVITY load applied to beams and columns -- eleLoad applies loads in local coordinate axis +pattern Plain 101 Linear { + for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} { + for {set level 1} {$level <=$NStory} {incr level 1} { + for {set pier 1} {$pier <= [expr $NBay+1]} {incr pier 1} { + set elemID [expr $N0col + $level*$Dlevel +$frame*$Dframe+$pier] + eleLoad -ele $elemID -type -beamUniform 0. 0. -$QdlCol; # COLUMNS } + } + } + for {set frame 1} {$frame <=[expr $NFrame]} {incr frame 1} { + for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} { + for {set bay 1} {$bay <= $NBay} {incr bay 1} { + set elemID [expr $N0beam + $level*$Dlevel +$frame*$Dframe+ $bay] + eleLoad -ele $elemID -type -beamUniform -$QdlBeam 0.; # BEAMS + } + } + } + for {set frame 1} {$frame <=[expr $NFrame-1]} {incr frame 1} { + for {set level 2} {$level <=[expr $NStory+1]} {incr level 1} { + for {set bay 1} {$bay <= $NBay+1} {incr bay 1} { + set elemID [expr $N0gird + $level*$Dlevel +$frame*$Dframe+ $bay] + eleLoad -ele $elemID -type -beamUniform -$QdlGird 0.; # GIRDS + } + } + } + +} + + + +# write beam nodes to a file +set StructureNodesFile "$outputdir/StructureNodes.dat" +if {[file exists $StructureNodesFile] == 1} { file delete $StructureNodesFile } +set StructureNodesID [open $StructureNodesFile "w"] +foreach node $StructureNodes { + puts $StructureNodesID "$node [nodeCoord $node]" +} +close $StructureNodesID + +# write beam elements to a file +set StructureElemsFile "$outputdir/StructureElements.dat" +if {[file exists $StructureElemsFile] == 1} { file delete $StructureElemsFile } +set StructureElemsID [open $StructureElemsFile "w"] +foreach ele $StructureElements { + puts $StructureElemsID "$ele [eleNodes $ele]" +} +close $StructureElemsID + + + + +# write baseColumns information to a file +set BaseColumnsFile "$outputdir/BaseColumnsNodes.dat" +if {[file exists $BaseColumnsFile] == 1} { file delete $BaseColumnsFile } +set BaseColumnsID [open $BaseColumnsFile "w"] +foreach node $BaseColumnsIDs { + puts $BaseColumnsID "$node [nodeCoord $node]" +} +close $BaseColumnsID + + + + +# # Gravity-analysis parameters -- load-controlled static analysis +# set Tol 1.0e-3; # convergence tolerance for test +# variable constraintsTypeGravity Plain; # default; +# if { [info exists RigidDiaphragm] == 1} { +# if {$RigidDiaphragm=="ON"} { +# variable constraintsTypeGravity Lagrange; # large model: try Transformation +# }; # if rigid diaphragm is on +# }; # if rigid diaphragm exists + +# constraints $constraintsTypeGravity ; # how it handles boundary conditions +# numberer RCM; # renumber dof's to minimize band-width (optimization), if you want to +# system BandGeneral ; # how to store and solve the system of equations in the analysis (large model: try UmfPack) +# test EnergyIncr $Tol 6 ; # determine if convergence has been achieved at the end of an iteration step +# algorithm Newton; # use Newton's solution algorithm: updates tangent stiffness at every iteration +# set NstepGravity 10; # apply gravity in 10 steps +# set DGravity [expr 1./$NstepGravity]; # first load increment; +# integrator LoadControl $DGravity; # determine the next time step for an analysis +# analysis Static; # define type of analysis static or transient +# analyze $NstepGravity; # apply gravity +# # ------------------------------------------------- maintain constant gravity loads and reset time to zero +# loadConst -time 0.0 + +# # ------------------------------------------------------------- +# puts "Model Built" + + diff --git a/modules/tools/DRM/modelCreator.py b/modules/tools/DRM/modelCreator.py new file mode 100644 index 000000000..281c62ed1 --- /dev/null +++ b/modules/tools/DRM/modelCreator.py @@ -0,0 +1,383 @@ +# %% +import numpy as np +from MeshGenerator.MeshGenrator import * +from MeshGenerator.Infowriter import * +import os +import shutil +import argparse + +# change the current directory to the directory of the script +abspath = os.path.abspath(__file__) +dname = os.path.dirname(abspath) +os.chdir(dname) + + +parser = argparse.ArgumentParser(description='Create a model for the DRM') +parser.add_argument('--soilwidth_x', type=float, help='width of the soil in x direction') +parser.add_argument('--soilwidth_y', type=float, help='width of the soil in y direction') +parser.add_argument('--soilwidth_z', type=float, help='width of the soil in z direction') +parser.add_argument('--soilNumEle_x', type=int, help='number of elements in x direction') +parser.add_argument('--soilNumEle_y', type=int, help='number of elements in y direction') +parser.add_argument('--soilNumEle_z', type=int, help='number of elements in z direction') +parser.add_argument('--Vs', type=float, help='shear wave velocity') +parser.add_argument('--nu', type=float, help='poisson ratio') +parser.add_argument('--rho', type=float, help='density') +parser.add_argument('--DRM_filePath', type=str, help='DRM file path') +parser.add_argument('--DRM_numLayers', type=int, help='number of DRM layers') +parser.add_argument('--DRM_loadFactor', type=float, help='DRM load factor') +parser.add_argument('--DRM_crdScale', type=float, help='DRM coordinate scale') +parser.add_argument('--DRM_tolerance', type=float, help='DRM tolerance') +parser.add_argument('--DRM_T00', type=float, help='DRM T00') +parser.add_argument('--DRM_T01', type=float, help='DRM T01') +parser.add_argument('--DRM_T02', type=float, help='DRM T02') +parser.add_argument('--DRM_T10', type=float, help='DRM T10') +parser.add_argument('--DRM_T11', type=float, help='DRM T11') +parser.add_argument('--DRM_T12', type=float, help='DRM T12') +parser.add_argument('--DRM_T20', type=float, help='DRM T20') +parser.add_argument('--DRM_T21', type=float, help='DRM T21') +parser.add_argument('--DRM_T22', type=float, help='DRM T22') +parser.add_argument('--DRM_originX', type=float, help='DRM origin X') +parser.add_argument('--DRM_originY', type=float, help='DRM origin Y') +parser.add_argument('--DRM_originZ', type=float, help='DRM origin Z') +parser.add_argument('--DRM_Software', type=str, help='DRM software') +parser.add_argument('--DRM_CoordinateTransformation', type=bool, help='DRM coordinate transformation') +parser.add_argument('--Absorb_HaveAbsorbingElements', type=str, help='Absorbing elements') +parser.add_argument('--Absorb_NumAbsorbingElements', type=int, help='Number of absorbing elements') +parser.add_argument('--Absorb_rayleighAlpha', type=float, help='Rayleigh alpha') +parser.add_argument('--Absorb_rayleighBeta', type=float, help='Rayleigh beta') +parser.add_argument('--Absorb_type', type=str, help='Absorbing type') +parser.add_argument('--dt', type=float, help='time step') +parser.add_argument('--t_final', type=float, help='final time') +parser.add_argument('--recording_dt', type=float, help='recording time step') +parser.add_argument('--outputDir', type=str, help='output directory') +parser.add_argument('--PartitionAlgorithm', type=str, help='partition algorithm') +parser.add_argument('--soilcores', type=int, help='number of cores for the soil') +parser.add_argument('--absorbingcores', type=int, help='number of cores for the PML') +parser.add_argument('--drmcores', type=int, help='number of cores for the DRM') +parser.add_argument('--DRM_Location', type=str, help='DRM location') + +args = parser.parse_args() + +# ============================================================================ +# Cores information +# ============================================================================ +regcores = args.soilcores +pmlcores = args.absorbingcores +drmcores = args.drmcores +structurecores = 0 +AnalysisType = "PML" # options: "DRMFIXED", PMLDRM", "ASDADRM" +Target = "Soil-with-structure" # options: "Soil", "Structure", "Soil-with-structure" +PartitionAlgorithm = args.PartitionAlgorithm # options: "kd-tree", "metis" +# ============================================================================ +# Structure information +# ============================================================================ +StructureType = "STEEL" # Options: STEEL, CONCRETE, Custom +NStory = 10 # number of stories above ground level +NBay = 4 # number of bays in X direction +NBayZ = 4 # number of bays in Y direction +StartNodeX = -9.0 # X coordinate of the first node +StartNodeY = -9.0 # Y coordinate of the first node +StartNodeZ = -1.5 # Z coordinate of the first node +meter = 1.0 # meter to specified unit conversion (SI unit) +ft = 0.3048 # feet to meter conversion (SI unit) +LCol = 3*meter # column height (parallel to Z axis) +LBeam = 4.5*meter # beam length (parallel to X axis) +LGird = 4.5*meter # girder length (parallel to Y axis) +SectionType = "Elastic" # options: Elastic, FiberSection +HaveStructure = "NO" # options: "YES", "NO" + +# ============================================================================ +# Soil information +# ============================================================================ +dx = args.soilwidth_x/args.soilNumEle_x +dy = args.soilwidth_y/args.soilNumEle_y +dz = args.soilwidth_z/args.soilNumEle_z +lx = args.soilwidth_x +ly = args.soilwidth_y +lz = args.soilwidth_z +drmthicknessx = dx +drmthicknessy = dy +drmthicknessz = dz +numdrmlayers = args.DRM_numLayers +llx = lx + 2*numdrmlayers*drmthicknessx +lly = ly + 2*numdrmlayers*drmthicknessy +llz = lz + 1*numdrmlayers*drmthicknessz + +nx = args.soilNumEle_x +ny = args.soilNumEle_y +nz = args.soilNumEle_z + +# modelCreator.py: error: unrecognized arguments: --soilwidth_x 100 --soilwidth_y 100 --soilwidth_z 30 --soilNumEle_x 40 --soilNumEle_y 40 --soilNumEle_z 12 --Vs 200 --nu 0.3 --rho 2000 --DRM_filePath /mnt/d/Projects/OpenSeesProjects/EEUQ_DRMBackend/DRMload/DRMLoad.h5drm --DRM_numLayers 3 --DRM_laodFactor 1 --DRM_crdScale 1 --DRM_tolernace 1 --DRM_T00 0.0 --DRM_T01 1.0 --DRM_T02 0.0 --DRM_T10 1.0 --DRM_T11 0.0 --DRM_T12 0.0 --DRM_T20 0.0 --DRM_T21 0.0 --DRM_T22 -1.0 --DRM_originX 0 --DRM_originY 0 --DRM_originZ 0 --DRM_Sofrware ShakerMaker --Absorb_HaveAbsorbingElements NO --Absorb_NumAbsorbingElements 2 --Absorb_type PML --Absorb_rayleighAlpha 0.0 --Absorb_rayleighBeta 0.0 --dt 0.005 --t_final 60 --recording_dt 0.005 --outputDir /home/amnp95/Documents/EE-UQ/LocalWorkDir/DRM_Model +# ============================================================================ +# PML information +# ============================================================================ +AbsorbingElements = args.Absorb_type # could be "ASDA" or "NORMAL" +numpmllayers = args.Absorb_NumAbsorbingElements +pmlthicknessx = dx +pmlthicknessy = dy +pmlthicknessz = dz +pmltotalthickness = numpmllayers*pmlthicknessx +HaveAbsorbingElements = args.Absorb_HaveAbsorbingElements # options: "YES", "NO" +Absorbing_rayleigh_alpha = args.Absorb_rayleighAlpha +Absorbing_rayleigh_beta = args.Absorb_rayleighBeta + + +# ============================================================================ +# General information +# ============================================================================ +tmpdir = args.outputDir +meshdir = f"{tmpdir}/Mesh" +outputdir = f"{tmpdir}/Results" +DRMFile = args.DRM_filePath + +# ============================================================================ +# Embedding foundation +# ============================================================================ +EmbeddingFoundation = "NO" # options: "YES", "NO" +EmbeddedFoundation = { + "xmax": 10.0, + "xmin": -10.0, + "ymax": 10.0, + "ymin": -10.0, + "zmax": 0.0, + "zmin": -5.0, +} + +# ============================================================================ +# Fondation information +# ============================================================================ +HaveFoundation = "NO" +AttachFoundation = "NO" +foundationBlocks = [] + +# adding a block +foundationBlocks.append({ + "matTag": 1, + "xmax": 10.0, + "xmin": -10.0, + "ymax": 10.0, + "ymin": -10.0, + "zmax": -1.5, + "zmin": -4.5, + "Xmeshsize" : 1.0, + "Ymeshsize" : 1.0, + "Zmeshsize" : 1.0, +}) + + +# ============================================================================ +# piles information +# ============================================================================ +pilelist = [] + +x = np.arange(-7.0, 7.0+1e-6, 3.5) +y = np.arange(-7.0, 7.0+1e-6, 3.5) + +x, y = np.meshgrid(x, y) +x = x.flatten() +y = y.flatten() + +HavePiles = "NO" +for i in range(len(x)): + pilelist.append({ + "xtop": x[i], + "ytop": y[i], + "ztop": -3.0, + "xbottom": x[i], + "ybottom": y[i], + "zbottom": -10.0, + "numberofElements": 6, + }) + +# ============================================================================ +# DRM information +# ============================================================================ +# DRM information +DRMinformation = { + "DRM_Provider_Software" : args.DRM_Software, + "factor": args.DRM_loadFactor, + "crd_scale": args.DRM_crdScale, + "distance_tolerance": args.DRM_tolerance, + "do_coordinate_transformation": 1, + "T00":args.DRM_T00, + "T01":args.DRM_T01, + "T02":args.DRM_T02, + "T10":args.DRM_T10, + "T11":args.DRM_T11, + "T12":args.DRM_T12, + "T20":args.DRM_T20, + "T21":args.DRM_T21, + "T22":args.DRM_T22, + "originX": args.DRM_originX, + "originY": args.DRM_originY, + "originZ": args.DRM_originZ, + "DRM_Location": args.DRM_Location, +} + + +# ============================================================================ +# analysis information +# ============================================================================ +dt = args.dt +t_final = args.t_final +recording_dt = args.recording_dt + +analysisinfo = { + "dt": dt, + "t_final": t_final, + "recording_dt": recording_dt, +} + + +# ============================================================================ +# material information +# ============================================================================ +# creating the material information +materialInfo = { + "Vs": args.Vs, + "nu": args.nu, + "rho": args.rho, +} + +# ============================================================================ +# Create directories +# ============================================================================ +# delete the directories if exists +if os.path.exists(meshdir): + # removing directory + shutil.rmtree(meshdir, ignore_errors=False) + +if os.path.exists(outputdir): + shutil.rmtree(outputdir, ignore_errors=False) + +# creating directories using shutil +os.makedirs(meshdir) +os.makedirs(outputdir) + + + +# ============================================================================ +# Creating the mesh +# ============================================================================ +info = { + "xwidth": lx, + "ywidth": ly, + "zwidth": lz, + "Xmeshsize": dx, + "Ymeshsize": dy, + "Zmeshsize": dz, + "PMLThickness": [pmlthicknessx, pmlthicknessy, pmlthicknessz], + "numPMLLayers": numpmllayers, + "DRMThickness": [drmthicknessx, drmthicknessy, drmthicknessz], + "numDrmLayers": numdrmlayers, + "reg_num_cores": regcores, + "DRM_num_cores": drmcores, + "PML_num_cores": pmlcores, + "Structure_num_cores": structurecores, + "Dir": meshdir, + "OutputDir": outputdir, + "AbsorbingElements": AbsorbingElements, + "DRMfile": DRMFile, + "EmbeddingFoundation": EmbeddingFoundation, + "EmbeddedFoundationDict": EmbeddedFoundation, + "HaveFoundation": HaveFoundation, + "foundationBlocks": foundationBlocks, + "pilelist": pilelist, + "HavePiles": HavePiles, + "HaveStructure": HaveStructure, + "AttachFoundation": AttachFoundation, + "DRMinformation": DRMinformation, + "AnalysisInfo": analysisinfo, + "PartitionAlgorithm": PartitionAlgorithm, + "Absorbing_rayleigh_alpha": Absorbing_rayleigh_alpha, + "Absorbing_rayleigh_beta": Absorbing_rayleigh_beta, +} + +numcells, numpoints = DRM_PML_Foundation_Meshgenrator(info) +print(f"Number of cells: {numcells}") +print(f"Number of points: {numpoints}") +# ============================================================================ +# Writing the information file +# ============================================================================ +info = { + "soilfoundation_num_cells": numcells, + "soilfoundation_num_points": numpoints, + "AnalysisType": AnalysisType, + "regcores": regcores, + "pmlcores": pmlcores, + "drmcores": drmcores, + "structurecores": structurecores, + "StructureType": StructureType, + "NStory": NStory, + "NBay": NBay, + "NBayZ": NBayZ, + "StartNodeX": StartNodeX, + "StartNodeY": StartNodeY, + "StartNodeZ": StartNodeZ, + "meter": meter, + "ft": ft, + "LCol": LCol, + "LBeam": LBeam, + "LGird": LGird, + "SectionType": SectionType, + "dx": dx, + "dy": dy, + "dz": dz, + "llx": llx, + "lly": lly, + "llz": llz, + "drmthicknessx": drmthicknessx, + "drmthicknessy": drmthicknessy, + "drmthicknessz": drmthicknessz, + "numdrmlayers": numdrmlayers, + "lx": lx, + "ly": ly, + "lz": lz, + "nx": nx, + "ny": ny, + "nz": nz, + "AbsorbingElements": AbsorbingElements, + "numpmllayers": numpmllayers, + "pmlthicknessx": pmlthicknessx, + "pmlthicknessy": pmlthicknessy, + "pmlthicknessz": pmlthicknessz, + "pmltotalthickness": pmltotalthickness, + "HaveAbsorbingElements": HaveAbsorbingElements, + "meshdir": meshdir, + "outputdir": outputdir, + "DRMFile": args.DRM_filePath.split("/")[-1], + "EmbeddingFoundation": EmbeddingFoundation, + "EmbeddedFoundation": EmbeddedFoundation, + "HaveFoundation": HaveFoundation, + "foundationBlocks": foundationBlocks, + "pilelist": pilelist, + "HavePiles": HavePiles, + "HaveStructure": HaveStructure, + "DRMinformation": DRMinformation, + "Absorbing_rayleigh_alpha": Absorbing_rayleigh_alpha, + "Absorbing_rayleigh_beta": Absorbing_rayleigh_beta, + "AnalysisInfo": analysisinfo, + "MaterialInfo": materialInfo, + +} +infowriter(info,meshdir) + +# ============================================================================ +# copy the related file as model.tcl to the current directory +# ============================================================================ +def copy_file(source_path, destination_path): + with open(destination_path, 'wb') as dst_file: + with open(f"{meshdir}/Modelinfo.tcl", 'rb') as src_file: + dst_file.write(src_file.read()) + with open(source_path, 'rb') as src_file: + dst_file.write(src_file.read()) + + +# delete the model file if exists +if os.path.exists("./model.tcl"): + os.remove("./model.tcl") + +if Target == "Soil-with-structure": + copy_file(f"MeshGenerator/models/Soil_with_structure.tcl", f"{meshdir}/model.tcl") + +# %%