diff --git a/meshing_scripts/create_full_mesh.py b/meshing_scripts/create_full_mesh.py index 7b872dee..4ca8166f 100644 --- a/meshing_scripts/create_full_mesh.py +++ b/meshing_scripts/create_full_mesh.py @@ -14,6 +14,7 @@ from one_dim_mesh import OneDimMesh from natural_order_mesh import NatOrdMeshRow from mesh_utils import * +from remove_cells_step import * # ------------------------------------------------ def checkDomainBoundsCoordinates(xValues, yValues): @@ -44,6 +45,8 @@ def main(workDir, debug, Nx, Ny, \ # figure out if domain is a plain rectangle or has a step in it plainDomain = True + print(len(xBd)) + print(xBd) if dim!=1 and len(xBd) > 4: plainDomain = False @@ -59,22 +62,39 @@ def main(workDir, debug, Nx, Ny, \ axFM = figFM.gca() # --------------------------------------------- + x = None + y = None + G = None + gids = None if dim == 1: meshObj = OneDimMesh(Nx, dx, xBd[0], xBd[1], stencilSize, enablePeriodicBc) + # get mesh coordinates, gids and graph + [x, y] = meshObj.getCoordinates() + gids = meshObj.getGIDs() + G = meshObj.getGraph() elif dim==2 and plainDomain: meshObj = NatOrdMeshRow(Nx,Ny, dx,dy,\ xBd[0], xBd[1], yBd[0], yBd[1], \ stencilSize, enablePeriodicBc) - else: - print("invalid domain: with step not impl yet") - sys.exit(1) + # get mesh coordinates, gids and graph + [x, y] = meshObj.getCoordinates() + gids = meshObj.getGIDs() + G = meshObj.getGraph() + + elif dim==2 and not plainDomain: + minX, maxX = min(xBd), max(xBd) + minY, maxY = min(yBd), max(yBd) + L = [maxX-minX, maxY-minY] + dx,dy = L[0]/Nx, L[1]/Ny + + meshObj = NatOrdMeshRow(Nx, Ny, dx, dy,\ + minX, maxX, minY, maxY, \ + stencilSize, False) + + x,y,G,gids = removeStepCells(meshObj, xBd, yBd) - # get mesh coordinates, gids and graph - [x, y] = meshObj.getCoordinates() - gids = meshObj.getGIDs() - G = meshObj.getGraph() if debug: print("full mesh connectivity") printDicPretty(G) @@ -136,12 +156,13 @@ def main(workDir, debug, Nx, Ny, \ f.write("stencilSize %2d\n" % stencilSize) f.write("nx %8d\n" % Nx) f.close() + else: f.write("dim %1d\n" % 2) - f.write("xMin %.14f\n" % xBd[0]) - f.write("xMax %.14f\n" % xBd[1]) - f.write("yMin %.14f\n" % yBd[0]) - f.write("yMax %.14f\n" % yBd[1]) + f.write("xMin %.14f\n" % min(xBd)) + f.write("xMax %.14f\n" % max(xBd)) + f.write("yMin %.14f\n" % min(yBd)) + f.write("yMax %.14f\n" % max(yBd)) f.write("dx %.14f\n" % dx) f.write("dy %.14f\n" % dy) f.write("sampleMeshSize %8d\n" % meshSize) @@ -174,9 +195,9 @@ def main(workDir, debug, Nx, Ny, \ if plotting != "none": plotLabels(x, y, dx, dy, gids, axFM, fontSz=plotFontSize) axFM.set_aspect(1.0) - axFM.set_xlim(xBd[0], xBd[1]) + axFM.set_xlim(min(xBd), max(xBd)) if dim!=1: - axFM.set_ylim(yBd[0], yBd[1]) + axFM.set_ylim(min(yBd), max(yBd)) if plotting == "show": plt.show() @@ -273,6 +294,17 @@ def main(workDir, debug, Nx, Ny, \ yCoords = [args.bounds[2], args.bounds[3]] dim = 2 + elif nx!=1 and ny!=1 and len(args.bounds) == 12: + xCoords = [] + yCoords = [] + print(args.bounds) + for i,v in enumerate(args.bounds): + if i % 2 == 0: + xCoords.append(v) + else: + yCoords.append(v) + dim = 2 + # other things plotFontSize = 0 if len(args.plottingInfo) == 2: diff --git a/meshing_scripts/create_sample_mesh.py b/meshing_scripts/create_sample_mesh.py index 10013763..d46b294e 100644 --- a/meshing_scripts/create_sample_mesh.py +++ b/meshing_scripts/create_sample_mesh.py @@ -36,13 +36,9 @@ def readFullMeshConnec(fullMeshDir): lineList = line.split() numNeigh = len(lineList[1:]) connections = [np.int64(i) for i in lineList[1:]] - #np.zeros(numNeigh, dtype=np.int64) - #for i in range(numNeigh): - #connections[i] = np.int64(lineList[1+i]) pt = np.int64(lineList[0]) gids.append(pt) G[pt] = np.copy(connections) - line = fp.readline() cnt += 1 @@ -51,6 +47,7 @@ def readFullMeshConnec(fullMeshDir): #==================================================================== def readFullMeshInfo(fullMeshDir): + dim=2 bounds=[None]*4 dx,dy = 0., 0. sampleMeshSize = 0 @@ -104,6 +101,7 @@ def readFullMeshCoordinates(fullMeshDir): cnt += 1 return np.array(x), np.array(y) +#==================================================================== #==================================================================== def main(workDir, debug, fullMeshDir, tilingDir, plotting, plotFontSize): dx,dy,numCells,domainBounds,stencilSize = readFullMeshInfo(fullMeshDir) @@ -119,8 +117,8 @@ def main(workDir, debug, fullMeshDir, tilingDir, plotting, plotFontSize): if plotting != "none": plotLabels(x, y, dx, dy, gids, axFM, fontSz=plotFontSize) axFM.set_aspect(1.0) - axFM.set_xlim(np.min(x)-0.5,np.max(x)+0.5) - axFM.set_ylim(np.min(y)-0.5,np.max(y)+0.5) + axFM.set_xlim(np.min(x)-dx*0.5, np.max(x)+dx*0.5) + axFM.set_ylim(np.min(y)-dy*0.5, np.max(y)+dy*0.5) if debug: print("natural order full mesh connectivity") @@ -141,7 +139,7 @@ def main(workDir, debug, fullMeshDir, tilingDir, plotting, plotFontSize): # ----------------------------------------------------- smGraph0 = collections.OrderedDict() for rPt in sampleMeshGIDs: - smGraph0[rPt] = G[rPt] + smGraph0[rPt] = G[rPt] print("\n") if debug: print("sample mesh graph0 (IDs wrt full mesh)") @@ -150,10 +148,11 @@ def main(workDir, debug, fullMeshDir, tilingDir, plotting, plotFontSize): stencilMeshGIDs = [] # loop over target cells wherw we want residual for k,v in smGraph0.items(): - # append the GID of this cell - stencilMeshGIDs.append(k) - # append GID of stencils/neighborin cells - for j in v: + # append the GID of this cell + stencilMeshGIDs.append(k) + # append GID of stencils/neighborin cells + for j in v: + if j != -1: stencilMeshGIDs.append(np.int64(j)) # remove duplicates and sort @@ -180,6 +179,7 @@ def main(workDir, debug, fullMeshDir, tilingDir, plotting, plotFontSize): if tilingDir!=None: gidsfiles = glob.glob(tilingDir+"/cell_gids_p_*.txt") + # sort based on the ID, so need to extract ID which is last of dir name def func(elem): return int(elem.split('_')[-1].split('.')[0]) gidsfiles = sorted(gidsfiles, key=func) @@ -211,14 +211,14 @@ def func(elem): return int(elem.split('_')[-1].split('.')[0]) if debug: for k,v in fm_to_sm_map.items(): print(k, v) - print("doing now the sm -> fm gids mapping ") sm_to_fm_map = collections.OrderedDict() for k,v in fm_to_sm_map.items(): sm_to_fm_map[v] = k print("Done with sm_to_fm_map") if debug: - for k,v in sm_to_fm_map.items(): print(k, v) + for k,v in sm_to_fm_map.items(): + print(k, v) # ----------------------------------------------------- # Here we have a list of unique GIDs for the sample mesh. @@ -227,12 +227,13 @@ def func(elem): return int(elem.split('_')[-1].split('.')[0]) # ----------------------------------------------------- sampleMeshGraph = collections.OrderedDict() for rGidFM, v in smGraph0.items(): - smGID = fm_to_sm_map[rGidFM] - smStencilGIDs = v - for i in range(len(smStencilGIDs)): - thisGID = smStencilGIDs[i] - smStencilGIDs[i] = fm_to_sm_map[thisGID] - sampleMeshGraph[smGID] = smStencilGIDs + smGID = fm_to_sm_map[rGidFM] + smStencilGIDs = v + for i in range(len(smStencilGIDs)): + thisGID = smStencilGIDs[i] + if thisGID != -1: + smStencilGIDs[i] = fm_to_sm_map[thisGID] + sampleMeshGraph[smGID] = smStencilGIDs print("\n") print("Done with sampleMeshGraph") @@ -245,8 +246,9 @@ def func(elem): return int(elem.split('_')[-1].split('.')[0]) x2, y2 = x[ list(fm_to_sm_map.keys()) ], y[ list(fm_to_sm_map.keys()) ] plotLabels(x2, y2, dx, dy, gids_sm, axSM, 's', 'r', 0, fontSz=plotFontSize) axSM.set_aspect(1.0) - axSM.set_xlim(np.min(x)-0.5,np.max(x)+0.5) - axSM.set_ylim(np.min(y)-0.5,np.max(y)+0.5) + axSM.set_xlim(np.min(x)-dx*0.5, np.max(x)+dx*0.5) + axSM.set_ylim(np.min(y)-dy*0.5, np.max(y)+dy*0.5) + # ----------------------------------------------------- sampleMeshSize = len(sampleMeshGraph) diff --git a/meshing_scripts/natural_order_mesh.py b/meshing_scripts/natural_order_mesh.py index 07f82f02..84fdaf38 100644 --- a/meshing_scripts/natural_order_mesh.py +++ b/meshing_scripts/natural_order_mesh.py @@ -49,6 +49,9 @@ def __init__(self, Nx, Ny, dx, dy, \ self.buildGraph(enablePeriodicBc) #self.spMat_ = convertGraphDicToSparseMatrix(self.G_) + def getTotCells(self): + return self.Nx_ * self.Ny_ + def getCoordinates(self): return [self.x_, self.y_] diff --git a/meshing_scripts/remove_cells_step.py b/meshing_scripts/remove_cells_step.py new file mode 100644 index 00000000..459c00b7 --- /dev/null +++ b/meshing_scripts/remove_cells_step.py @@ -0,0 +1,60 @@ + +import matplotlib.pyplot as plt +import sys, os, time +import numpy as np +from numpy import linspace, meshgrid +from matplotlib import cm +import collections +from argparse import ArgumentParser +import random +import scipy.sparse as sp +from scipy.sparse import csr_matrix +from scipy.sparse.csgraph import reverse_cuthill_mckee + +from mesh_utils import printDicPretty + +def removeStepCells(meshObj, xBd, yBd): + + stepXb = [xBd[2], xBd[3]] + stepYb = [yBd[1], yBd[2]] + print("stepXb ", stepXb) + print("stepYb ", stepYb) + + totNumCells = meshObj.getTotCells() + [x, y] = meshObj.getCoordinates() + gids = meshObj.getGIDs() + G = meshObj.getGraph() + + a = np.where(x>stepXb[0]) + b = np.where(xstepYb[0]) + d = np.where(y