src/HYDROTools/shapesGroups.py

   1
   2 import salome
   3 salome.salome_init()
   4
   5 import  SMESH, SALOMEDS
   6 from salome.smesh import smeshBuilder
   7
   8 import numpy as np
   9 import MEDLoader as ml
  10 import medcoupling as mc
  11
  12 import shapefile
  13 import math
  14 import os
  15
  16
  17 def freeBordersGroup(meshFile):
  18     print(" === freeBordersGroup", meshFile)
  19     smesh = smeshBuilder.New()
  20     smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed
  21     ([MESH], status) = smesh.CreateMeshesFromMED(meshFile)
  22     nbAdded, MESH, addedBnd = MESH.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [])
  23     groups = MESH.GetGroups()
  24     aCriteria = []
  25     aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
  26     aCriteria.append(aCriterion)
  27     aFilter = smesh.GetFilterFromCriteria(aCriteria)
  28     aFilter.SetMesh(MESH.GetMesh())
  29     FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter )
  30     a = os.path.splitext(meshFile)
  31     smesh.SetName(MESH, os.path.basename(a[0]))
  32     newMeshName = a[0] + '_brd' + a[1]
  33     print(newMeshName)
  34     try:
  35         MESH.ExportMED(newMeshName,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1)
  36         pass
  37     except:
  38         print('ExportMED() failed. Invalid file name?')
  39     return newMeshName
  40
  41
  42 def explodeGroup(grp, grpName):
  43     print(" === explodeGroup", grpName)
  44     nbCells=grp.getNumberOfCells()
  45
  46     dicReverse = {} # id noeud --> id edges
  47     for i in range(nbCells):
  48         nodcell = grp.getNodeIdsOfCell(i)
  49         for j in range(len(nodcell)):
  50             if nodcell[j] in dicReverse:
  51                 dicReverse[nodcell[j]].append(i)
  52             else:
  53                 dicReverse[nodcell[j]] = [i]
  54
  55     nodeChains = []
  56     usedCells = [False] * nbCells
  57     while False in usedCells:
  58         icell = usedCells.index(False)
  59         usedCells[icell] = True
  60         nodcell = grp.getNodeIdsOfCell(icell)
  61         closed = False
  62         chain = [nodcell[0], nodcell[1]]
  63         nextnode = nodcell[1]
  64         prevnode = nodcell[0]
  65         while nextnode in dicReverse:
  66             nextcells = dicReverse[nextnode]
  67             if len(nextcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
  68                 closed = False
  69                 nextnode = -1                   # stop
  70             else:
  71                 newcell =False
  72                 for i in range(len(nextcells)):
  73                     ncell = nextcells[i]
  74                     if not usedCells[ncell]:       # the chain of nodes grows
  75                         usedCells[ncell] = True
  76                         newcell = True
  77                         nodcell = grp.getNodeIdsOfCell(ncell)
  78                         if nodcell[0] == nextnode:
  79                             nextnode = nodcell[1]  # forward edge
  80                         else:
  81                             nextnode = nodcell[0]  # reversed edge ?
  82                         chain.append(nextnode)
  83                 if not newcell:                    # end of chain, closed
  84                     closed =True
  85                     nextnode = -1
  86         while prevnode in dicReverse:
  87             prevcells = dicReverse[prevnode]
  88             if len(prevcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
  89                 closed = False
  90                 prevnode = -1                   # stop
  91             else:
  92                 newcell =False
  93                 for i in range(len(prevcells)):
  94                     ncell = prevcells[i]
  95                     if not usedCells[ncell]:       # the chain of nodes grows
  96                         usedCells[ncell] = True
  97                         newcell = True
  98                         nodcell = grp.getNodeIdsOfCell(ncell)
  99                         if nodcell[1] == prevnode:
 100                             prevnode = nodcell[0]  # forward edge
 101                         else:
 102                             prevnode = nodcell[1]  # reversed edge ?
 103                         chain.insert(0, prevnode)
 104                 if not newcell:                    # end of chain, closed
 105                     closed =True
 106                     prevnode = -1
 107
 108         chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed)
 109         nodeChains.append(chainDesc)
 110         print(chainDesc[1:])
 111     return nodeChains
 112
 113
 114 def writeShapeLines(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
 115     print(" === writeShapeLines", grpName)
 116     coords = mcMesh.getCoords()
 117     w = shapefile.Writer(grpName)
 118     w.shapeType = 3
 119     w.field('name', 'C')
 120     for (chain, chainName, closed) in nodeChains:
 121         print("   --- ", chainName)
 122         chaincoords = []
 123         for node in chain:
 124             coord = coords[node].getValues()
 125             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 126             #print("      ", coordLb93)
 127             chaincoords.append(coordLb93)
 128         w.line([chaincoords])
 129         w.record(chainName)
 130     w.close()
 131
 132
 133 def writeShapePoints(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
 134     print(" === writeShapePoints", grpName)
 135     coords = mcMesh.getCoords()
 136     w = shapefile.Writer(grpName + '_pts')
 137     #w.shapeType = 8
 138     w.field('name', 'C')
 139     for (chain, chainName, closed) in nodeChains:
 140         print("   --- ", chainName)
 141         chaincoords = []
 142         for node in chain:
 143             coord = coords[node].getValues()
 144             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 145             #print("      ", coordLb93)
 146             chaincoords.append(coordLb93)
 147         w.multipoint(chaincoords)
 148         w.record(chainName)
 149     w.close()
 150
 151
 152 def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.):
 153     print(" === exploreEdgeGroups", meshFile)
 154     mcMesh = ml.MEDFileMesh.New(meshFile)
 155     dim = mcMesh.getSpaceDimension()
 156     d1=-1        # when dimension 2, edges are dim -1
 157     if dim == 3: # when dimension 3, edges are dim -2
 158         d1=-2
 159     a = os.path.splitext(meshFile)
 160     prefix = os.path.basename(a[0])
 161
 162     grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
 163     groups = [mcMesh.getGroup(d1, name) for name in grp_names]
 164     for (grp, grpName) in zip(groups, grp_names):
 165         fullGrpName = prefix + '_' + grpName
 166         nodeChains = explodeGroup(grp, fullGrpName)
 167         writeShapeLines(mcMesh, fullGrpName, nodeChains, offsetX, offsetY)
 168         writeShapePoints(mcMesh, fullGrpName, nodeChains, offsetX, offsetY)
 169
 170
 171 def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
 172     print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
 173
 174     # --- find domain freeBorder:  bounding box englobing all others
 175
 176     fb = shapefile.Reader(freeBorderShapefile)
 177     fbShapes = fb.shapes()
 178     maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
 179     outerBboxIndex = -1
 180     for i,s in enumerate(fbShapes):
 181         bb = s.bbox
 182         if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
 183             maxbbox = bb
 184             outerBboxIndex = i
 185     fbs = fbShapes[outerBboxIndex]
 186
 187     # --- find the intersections of the shapefile to adjust and the domain free border:
 188     #     the closests points (two pairs of points)
 189
 190     sf = shapefile.Reader(shapefileToAdjust)
 191     shapes = sf.shapes()
 192     sfta = sf.shape(0)
 193     pdist =[]
 194     for i,p in enumerate(sfta.points):
 195         dmin = 1.e30
 196         jmin = -1
 197         for j,pfb in enumerate(fbs.points):
 198             d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
 199             if d < dmin:
 200                 dmin = d
 201                 jmin = j
 202         pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
 203     pdist.sort() # the 2 closest points must be set on the mesh freeBorder
 204     print(pdist)
 205     i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
 206     i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
 207     if i1 == pdist[0][2]:
 208         ifb1 = pdist[0][1]             # index of first adjusted point in freeBorder
 209         ifb2 = pdist[1][1]             # index of second adjusted point in freeBorder
 210     else:
 211         ifb1 = pdist[1][1]
 212         ifb2 = pdist[0][1]
 213
 214     # --- write the adusted shapefile: free border closest points replaced with corresponding points
 215     #     on the free border. two polylines, one inside the domain, one outside
 216
 217     a = os.path.splitext(shapefileToAdjust)
 218     shapefileAdjusted = a[0] + '_adj' + a[1]
 219     chainName = os.path.basename(a[0])
 220
 221     w = shapefile.Writer(shapefileAdjusted)
 222     w.shapeType = 3
 223     w.field('name', 'C')
 224
 225     chaincoords = []
 226     chaincoords.append(fbs.points[ifb1])
 227     for i in range(i1+1, i2):
 228         chaincoords.append(sfta.points[i])
 229     chaincoords.append(fbs.points[ifb2])
 230     w.line([chaincoords])
 231     w.record(chainName + '_0')
 232
 233     chaincoords = []
 234     chaincoords.append(fbs.points[ifb2])
 235     if i2+1 < len(sfta.points):
 236         for i in range(i2+1, len(sfta.points)):
 237             chaincoords.append(sfta.points[i])
 238     for i in range(i1):
 239         chaincoords.append(sfta.points[i])
 240     chaincoords.append(fbs.points[ifb1])
 241     w.line([chaincoords])
 242     w.record(chainName + '_1')
 243     w.close()
 244
 245     # write the free border splited in two polylines (cut by the adjusted shapefile)
 246
 247     a = os.path.splitext(freeBorderShapefile)
 248     freeBorderSplit = a[0] + '_split' + a[1]
 249     chainName = os.path.basename(a[0])
 250
 251     w = shapefile.Writer(freeBorderSplit)
 252     w.shapeType = 3
 253     w.field('name', 'C')
 254
 255     if (ifb1 > ifb2):
 256         i = ifb1; ifb1 = ifb2; ifb2 = i
 257
 258     chaincoords = []
 259     for i in range(ifb1, ifb2+1):
 260         chaincoords.append(fbs.points[i])
 261     w.line([chaincoords])
 262     w.record(chainName + '_0')
 263
 264     chaincoords = []
 265     for i in range(ifb2, len(fbs.points)):
 266         chaincoords.append(fbs.points[i])
 267     for i in range(ifb1+1):
 268         chaincoords.append(fbs.points[i])
 269     w.line([chaincoords])
 270     w.record(chainName + '_1')
 271     w.close()
 272