src/HYDROTools/shapesGroups.py

   1
   2
   3 # -----------------------------------------------------------------------------
   4
   5 import salome
   6 salome.salome_init()
   7
   8 import  SMESH, SALOMEDS
   9 from salome.smesh import smeshBuilder
  10
  11 import numpy as np
  12 import MEDLoader as ml
  13 import medcoupling as mc
  14
  15 import shapefile
  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 def explodeGroup(grp, grpName):
  39     print(" === explodeGroup", grpName)
  40     nbCells=grp.getNumberOfCells()
  41
  42     dicReverse = {} # id noeud --> id edges
  43     for i in range(nbCells):
  44         nodcell = grp.getNodeIdsOfCell(i)
  45         for j in range(len(nodcell)):
  46             if nodcell[j] in dicReverse:
  47                 dicReverse[nodcell[j]].append(i)
  48             else:
  49                 dicReverse[nodcell[j]] = [i]
  50
  51     nodeChains = []
  52     usedCells = [False] * nbCells
  53     while False in usedCells:
  54         icell = usedCells.index(False)
  55         usedCells[icell] = True
  56         nodcell = grp.getNodeIdsOfCell(icell)
  57         closed = False
  58         chain = [nodcell[0], nodcell[1]]
  59         nextnode = nodcell[1]
  60         prevnode = nodcell[0]
  61         while nextnode in dicReverse:
  62             nextcells = dicReverse[nextnode]
  63             if len(nextcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
  64                 closed = False
  65                 nextnode = -1                   # stop
  66             else:
  67                 newcell =False
  68                 for i in range(len(nextcells)):
  69                     ncell = nextcells[i]
  70                     if not usedCells[ncell]:       # the chain of nodes grows
  71                         usedCells[ncell] = True
  72                         newcell = True
  73                         nodcell = grp.getNodeIdsOfCell(ncell)
  74                         if nodcell[0] == nextnode:
  75                             nextnode = nodcell[1]  # forward edge
  76                         else:
  77                             nextnode = nodcell[0]  # reversed edge ?
  78                         chain.append(nextnode)
  79                 if not newcell:                    # end of chain, closed
  80                     closed =True
  81                     nextnode = -1
  82         while prevnode in dicReverse:
  83             prevcells = dicReverse[prevnode]
  84             if len(prevcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
  85                 closed = False
  86                 prevnode = -1                   # stop
  87             else:
  88                 newcell =False
  89                 for i in range(len(prevcells)):
  90                     ncell = prevcells[i]
  91                     if not usedCells[ncell]:       # the chain of nodes grows
  92                         usedCells[ncell] = True
  93                         newcell = True
  94                         nodcell = grp.getNodeIdsOfCell(ncell)
  95                         if nodcell[1] == prevnode:
  96                             prevnode = nodcell[0]  # forward edge
  97                         else:
  98                             prevnode = nodcell[1]  # reversed edge ?
  99                         chain.insert(0, prevnode)
 100                 if not newcell:                    # end of chain, closed
 101                     closed =True
 102                     prevnode = -1
 103
 104         chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed)
 105         nodeChains.append(chainDesc)
 106         print(chainDesc[1:])
 107     return nodeChains
 108
 109 def writeShapeLines(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
 110     print(" === writeShapeLines", grpName)
 111     coords = mcMesh.getCoords()
 112     w = shapefile.Writer(grpName)
 113     w.shapeType = 3
 114     w.field('name', 'C')
 115     for (chain, chainName, closed) in nodeChains:
 116         print("   --- ", chainName)
 117         chaincoords = []
 118         for node in chain:
 119             coord = coords[node].getValues()
 120             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 121             #print("      ", coordLb93)
 122             chaincoords.append(coordLb93)
 123         w.line([chaincoords])
 124         w.record(chainName)
 125     w.close()
 126
 127 def writeShapePoints(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
 128     print(" === writeShapePoints", grpName)
 129     coords = mcMesh.getCoords()
 130     w = shapefile.Writer(grpName + '_pts')
 131     #w.shapeType = 8
 132     w.field('name', 'C')
 133     for (chain, chainName, closed) in nodeChains:
 134         print("   --- ", chainName)
 135         chaincoords = []
 136         for node in chain:
 137             coord = coords[node].getValues()
 138             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 139             #print("      ", coordLb93)
 140             chaincoords.append(coordLb93)
 141         w.multipoint(chaincoords)
 142         w.record(chainName)
 143     w.close()
 144
 145 def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.):
 146     print(" === exploreEdgeGroups", meshFile)
 147     mcMesh = ml.MEDFileMesh.New(meshFile)
 148     dim = mcMesh.getSpaceDimension()
 149     d1=-1        # when dimension 2, edges are dim -1
 150     if dim == 3: # when dimension 3, edges are dim -2
 151         d1=-2
 152
 153     grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
 154     groups = [mcMesh.getGroup(d1, name) for name in grp_names]
 155     for (grp, grpName) in zip(groups, grp_names):
 156         nodeChains = explodeGroup(grp, grpName)
 157         writeShapeLines(mcMesh, grpName, nodeChains, offsetX, offsetY)
 158         writeShapePoints(mcMesh, grpName, nodeChains, offsetX, offsetY)
 159
 160 # ---
 161
 162 #meshFile = freeBordersGroup('/home/paul/projets/hydro95/V9_5_BR/tests/Maill_SB1610_ALB1501_barr_eff.med')
 163 #exploreEdgeGroups(meshFile, 1000000., 6000000.)