From 89350c4ade6c812fa42c7dcfb7c349b3c64a8f1c Mon Sep 17 00:00:00 2001 From: Paul RASCLE Date: Wed, 8 Jul 2020 17:54:08 +0200 Subject: [PATCH] fonctions Python conversion groupes en shapefiles --- src/HYDROTools/CMakeLists.txt | 1 + src/HYDROTools/shapesGroups.py | 163 +++++++++++++++++++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 src/HYDROTools/shapesGroups.py diff --git a/src/HYDROTools/CMakeLists.txt b/src/HYDROTools/CMakeLists.txt index cceff83c..1074e2d9 100644 --- a/src/HYDROTools/CMakeLists.txt +++ b/src/HYDROTools/CMakeLists.txt @@ -24,6 +24,7 @@ SET(PYFILES controls.py interpolZ.py interpolS.py + shapesGroups.py ) # --- rules --- diff --git a/src/HYDROTools/shapesGroups.py b/src/HYDROTools/shapesGroups.py new file mode 100644 index 00000000..14f6e8ee --- /dev/null +++ b/src/HYDROTools/shapesGroups.py @@ -0,0 +1,163 @@ + + +# ----------------------------------------------------------------------------- + +import salome +salome.salome_init() + +import SMESH, SALOMEDS +from salome.smesh import smeshBuilder + +import numpy as np +import MEDLoader as ml +import medcoupling as mc + +import shapefile + +def freeBordersGroup(meshFile): + print(" === freeBordersGroup", meshFile) + smesh = smeshBuilder.New() + smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed + ([MESH], status) = smesh.CreateMeshesFromMED(meshFile) + groups = MESH.GetGroups() + aCriteria = [] + aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0) + aCriteria.append(aCriterion) + aFilter = smesh.GetFilterFromCriteria(aCriteria) + aFilter.SetMesh(MESH.GetMesh()) + FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter ) + smesh.SetName(MESH, 'MESH') + newMeshName = '/tmp/freeBorders.med' + try: + MESH.ExportMED(newMeshName,auto_groups=0,minor=41,overwrite=1,meshPart=None,autoDimension=1) + pass + except: + print('ExportMED() failed. Invalid file name?') + return newMeshName + +def explodeGroup(grp, grpName): + print(" === explodeGroup", grpName) + nbCells=grp.getNumberOfCells() + + dicReverse = {} # id noeud --> id edges + for i in range(nbCells): + nodcell = grp.getNodeIdsOfCell(i) + for j in range(len(nodcell)): + if nodcell[j] in dicReverse: + dicReverse[nodcell[j]].append(i) + else: + dicReverse[nodcell[j]] = [i] + + nodeChains = [] + usedCells = [False] * nbCells + while False in usedCells: + icell = usedCells.index(False) + usedCells[icell] = True + nodcell = grp.getNodeIdsOfCell(icell) + closed = False + chain = [nodcell[0], nodcell[1]] + nextnode = nodcell[1] + prevnode = nodcell[0] + while nextnode in dicReverse: + nextcells = dicReverse[nextnode] + if len(nextcells) != 2: # end of chain(1) or "edges connector"(>2): stop + closed = False + nextnode = -1 # stop + else: + newcell =False + for i in range(len(nextcells)): + ncell = nextcells[i] + if not usedCells[ncell]: # the chain of nodes grows + usedCells[ncell] = True + newcell = True + nodcell = grp.getNodeIdsOfCell(ncell) + if nodcell[0] == nextnode: + nextnode = nodcell[1] # forward edge + else: + nextnode = nodcell[0] # reversed edge ? + chain.append(nextnode) + if not newcell: # end of chain, closed + closed =True + nextnode = -1 + while prevnode in dicReverse: + prevcells = dicReverse[prevnode] + if len(prevcells) != 2: # end of chain(1) or "edges connector"(>2): stop + closed = False + prevnode = -1 # stop + else: + newcell =False + for i in range(len(prevcells)): + ncell = prevcells[i] + if not usedCells[ncell]: # the chain of nodes grows + usedCells[ncell] = True + newcell = True + nodcell = grp.getNodeIdsOfCell(ncell) + if nodcell[1] == prevnode: + prevnode = nodcell[0] # forward edge + else: + prevnode = nodcell[1] # reversed edge ? + chain.insert(0, prevnode) + if not newcell: # end of chain, closed + closed =True + prevnode = -1 + + chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed) + nodeChains.append(chainDesc) + print(chainDesc[1:]) + return nodeChains + +def writeShapeLines(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.): + print(" === writeShapeLines", grpName) + coords = mcMesh.getCoords() + w = shapefile.Writer(grpName) + w.shapeType = 3 + w.field('name', 'C') + for (chain, chainName, closed) in nodeChains: + print(" --- ", chainName) + chaincoords = [] + for node in chain: + coord = coords[node].getValues() + coordLb93=[coord[0] + offsetX, coord[1] + offsetY] + #print(" ", coordLb93) + chaincoords.append(coordLb93) + w.line([chaincoords]) + w.record(chainName) + w.close() + +def writeShapePoints(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.): + print(" === writeShapePoints", grpName) + coords = mcMesh.getCoords() + w = shapefile.Writer(grpName + '_pts') + #w.shapeType = 8 + w.field('name', 'C') + for (chain, chainName, closed) in nodeChains: + print(" --- ", chainName) + chaincoords = [] + for node in chain: + coord = coords[node].getValues() + coordLb93=[coord[0] + offsetX, coord[1] + offsetY] + #print(" ", coordLb93) + chaincoords.append(coordLb93) + w.multipoint(chaincoords) + w.record(chainName) + w.close() + +def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.): + print(" === exploreEdgeGroups", meshFile) + mcMesh = ml.MEDFileMesh.New(meshFile) + dim = mcMesh.getSpaceDimension() + d1=-1 # when dimension 2, edges are dim -1 + if dim == 3: # when dimension 3, edges are dim -2 + d1=-2 + + grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups + groups = [mcMesh.getGroup(d1, name) for name in grp_names] + for (grp, grpName) in zip(groups, grp_names): + nodeChains = explodeGroup(grp, grpName) + writeShapeLines(mcMesh, grpName, nodeChains, offsetX, offsetY) + writeShapePoints(mcMesh, grpName, nodeChains, offsetX, offsetY) + +# --- + +#meshFile = freeBordersGroup('/home/paul/projets/hydro95/V9_5_BR/tests/Maill_SB1610_ALB1501_barr_eff.med') +#exploreEdgeGroups(meshFile, 1000000., 6000000.) -- 2.39.2