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