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     fb = shapefile.Reader(freeBorderShapefile)
 168     fbShapes = fb.shapes()
 169     maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
 170     outerBboxIndex = -1
 171     for i,s in enumerate(fbShapes):
 172         bb = s.bbox
 173         if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
 174             maxbbox = bb
 175             outerBboxIndex = i
 176     fbs = fbShapes[outerBboxIndex]
 177
 178     sf = shapefile.Reader(shapefileToAdjust)
 179     shapes = sf.shapes()
 180     sfta = sf.shape(0)
 181     pdist =[]
 182     for i,p in enumerate(sfta.points):
 183         dmin = 1.e30
 184         jmin = -1
 185         for j,pfb in enumerate(fbs.points):
 186             d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
 187             if d < dmin:
 188                 dmin = d
 189                 jmin = j
 190         pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
 191     pdist.sort() # the 2 closest points must be set on the mesh freeBorder
 192     print(pdist)
 193     i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
 194     i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
 195     if i1 == pdist[0][2]:
 196         ifb1 = pdist[0][1]             # index of first adjusted point in freeBorder
 197         ifb2 = pdist[1][1]             # index of second adjusted point in freeBorder
 198     else:
 199         ifb1 = pdist[1][1]
 200         ifb2 = pdist[0][1]
 201
 202     a = os.path.splitext(shapefileToAdjust)
 203     shapefileAdjusted = a[0] + '_adj' + a[1]
 204     chainName = os.path.basename(a[0])
 205
 206     w = shapefile.Writer(shapefileAdjusted)
 207     w.shapeType = 3
 208     w.field('name', 'C')
 209
 210     chaincoords = []
 211     chaincoords.append(fbs.points[ifb1])
 212     for i in range(i1+1, i2):
 213         chaincoords.append(sfta.points[i])
 214     chaincoords.append(fbs.points[ifb2])
 215     w.line([chaincoords])
 216     w.record(chainName + '_0')
 217
 218     chaincoords = []
 219     chaincoords.append(fbs.points[ifb2])
 220     for i in range(i2+1, len(sfta.points)):
 221         chaincoords.append(sfta.points[i])
 222     for i in range(i1):
 223         chaincoords.append(sfta.points[i])
 224     chaincoords.append(fbs.points[ifb1])
 225     w.line([chaincoords])
 226     w.record(chainName + '_1')
 227     w.close()
 228