6 from salome.smesh import smeshBuilder
10 import medcoupling as mc
17 def freeBordersGroup(meshFileIn, meshFileOut=""):
19 In a mesh, create a group of Edges for the borders: domain, internal holes (isles)
21 meshFileIn: full path of the input mesh file, format MED
22 meshFileOut: full path of the output mesh file, format MED (default="" : when "", output file is suffixed with "_brd.med"
23 return full path of the output mesh file
25 smesh = smeshBuilder.New()
26 smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed
27 ([MESH], status) = smesh.CreateMeshesFromMED(meshFileIn)
29 nbAdded, MESH, addedBnd = MESH.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [])
32 aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
33 aCriteria.append(aCriterion)
34 aFilter = smesh.GetFilterFromCriteria(aCriteria)
35 aFilter.SetMesh(MESH.GetMesh())
36 FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter )
38 a = os.path.splitext(meshFileIn)
39 smesh.SetName(MESH, os.path.basename(a[0]))
42 a = os.path.splitext(meshFileIn)
43 smesh.SetName(MESH, os.path.basename(a[0]))
44 meshFileOut = a[0] + '_brd' + a[1]
47 MESH.ExportMED(meshFileOut,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1)
50 print('ExportMED() failed. Invalid file name?')
54 def explodeGroup(grp, grpName):
56 from a group of edges loaded with MEDCoupling, create ordered lists of nodes, one list for each set of connected edges.
58 grp: MEDCoupling object for the group of edges
59 grpName: name of the group
61 List of descriptors [(ordered list of nodeIds, name of the list, closed status)]
63 print(" === explodeGroup", grpName)
64 nbCells=grp.getNumberOfCells()
66 dicReverse = {} # id node --> id edges
67 for i in range(nbCells):
68 nodcell = grp.getNodeIdsOfCell(i)
69 for j in range(len(nodcell)):
70 if nodcell[j] in dicReverse:
71 dicReverse[nodcell[j]].append(i)
73 dicReverse[nodcell[j]] = [i]
76 usedCells = [False] * nbCells
77 while False in usedCells:
78 icell = usedCells.index(False)
79 usedCells[icell] = True
80 nodcell = grp.getNodeIdsOfCell(icell)
82 chain = [nodcell[0], nodcell[1]]
85 while nextnode in dicReverse:
86 nextcells = dicReverse[nextnode]
87 if len(nextcells) != 2: # end of chain(1) or "edges connector"(>2): stop
92 for i in range(len(nextcells)):
94 if not usedCells[ncell]: # the chain of nodes grows
95 usedCells[ncell] = True
97 nodcell = grp.getNodeIdsOfCell(ncell)
98 if nodcell[0] == nextnode:
99 nextnode = nodcell[1] # forward edge
101 nextnode = nodcell[0] # reversed edge ?
102 chain.append(nextnode)
103 if not newcell: # end of chain, closed
106 while prevnode in dicReverse:
107 prevcells = dicReverse[prevnode]
108 if len(prevcells) != 2: # end of chain(1) or "edges connector"(>2): stop
113 for i in range(len(prevcells)):
115 if not usedCells[ncell]: # the chain of nodes grows
116 usedCells[ncell] = True
118 nodcell = grp.getNodeIdsOfCell(ncell)
119 if nodcell[1] == prevnode:
120 prevnode = nodcell[0] # forward edge
122 prevnode = nodcell[1] # reversed edge ?
123 chain.insert(0, prevnode)
124 if not newcell: # end of chain, closed
128 chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed)
129 nodeChains.append(chainDesc)
134 def writeShapeLines(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
136 from a mesh loaded in memory with MEDLoader, and a list of list of connected nodes associated to a group of edges, write a shapefile of type line, with one record per list of connected nodes.
138 mcMesh: mesh loaded in memory with MEDLoader
139 grpName: name associated to the group of edges
140 nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
141 outputDirectory: directory for writing the shapefile
142 offsetX : offset of origin X in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
143 offsetY : offset of origin Y in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
145 print(" === writeShapeLines", grpName)
146 coords = mcMesh.getCoords()
147 shapeFileName = os.path.join(outputDirectory, grpName)
148 w = shapefile.Writer(shapeFileName)
151 for (chain, chainName, closed) in nodeChains:
152 print(" --- ", chainName)
155 coord = coords[node].getValues()
156 coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
157 chaincoords.append(coordLb93)
158 w.line([chaincoords])
163 def writeShapePoints(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
165 from a mesh loaded in memory with MEDLoader, and a list of list of connected nodes associated to a group of edges, write a shapefile of type multi points, with one record per list of connected nodes.
167 mcMesh: mesh loaded in memory with MEDLoader
168 grpName: name associated to the group of edges
169 nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
170 outputDirectory: directory for writing the shapefile
171 offsetX : offset of origin X in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
172 offsetY : offset of origin Y in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
174 print(" === writeShapePoints", grpName)
175 coords = mcMesh.getCoords()
176 shapeFileName = os.path.join(outputDirectory, grpName)
177 w = shapefile.Writer(shapeFileName + '_pts')
180 for (chain, chainName, closed) in nodeChains:
181 print(" --- ", chainName)
184 coord = coords[node].getValues()
185 coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
186 chaincoords.append(coordLb93)
187 w.multipoint(chaincoords)
192 def exploreEdgeGroups(meshFile, outputDirectory="", offsetX=0., offsetY=0.):
194 Find all the groups of edges in a mesh and, for each group, create one shapefile of lines and one of points. The shapefiles are created in the same system of coordinates as the mesh (For instance Lambert 93), but without origin offset (to be ready for a direct load in Qgis)
196 meshFile: full path of the input mesh file, format MED
197 outputDirectory: directory in which the shapefiles are written (default "", if "" use the directory containing the mesh
198 offsetX: local X origin of the mesh
199 offsetY: local Y origin of the mesh
201 print(" === exploreEdgeGroups", meshFile)
202 if outputDirectory == "":
203 outputDirectory = os.path.dirname(meshFile)
205 a = os.path.splitext(meshFile)
206 prefix = os.path.basename(a[0]) # prefix = file name without extension
208 mcMesh = ml.MEDFileMesh.New(meshFile)
209 dim = mcMesh.getSpaceDimension()
210 d1=-1 # when dimension 2, edges are dim -1
211 if dim == 3: # when dimension 3, edges are dim -2
214 grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
216 groups = [mcMesh.getGroup(d1, name) for name in grp_names] # list of groups in their name order
218 for (grp, grpName) in zip(groups, grp_names):
219 fullGrpName = prefix + '_' + grpName
220 nodeChains = explodeGroup(grp, fullGrpName)
221 writeShapeLines(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
222 writeShapePoints(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
225 def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
226 print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
228 # --- find domain freeBorder: bounding box englobing all others
230 fb = shapefile.Reader(freeBorderShapefile)
231 fbShapes = fb.shapes()
232 maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
234 for i,s in enumerate(fbShapes):
236 if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
239 fbs = fbShapes[outerBboxIndex]
241 # --- find the intersections of the shapefile to adjust and the domain free border:
242 # the closests points (two pairs of points)
244 sf = shapefile.Reader(shapefileToAdjust)
248 for i,p in enumerate(sfta.points):
251 for j,pfb in enumerate(fbs.points):
252 d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
256 pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
257 pdist.sort() # the 2 closest points must be set on the mesh freeBorder
259 i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
260 i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
261 if i1 == pdist[0][2]:
262 ifb1 = pdist[0][1] # index of first adjusted point in freeBorder
263 ifb2 = pdist[1][1] # index of second adjusted point in freeBorder
268 # --- write the adusted shapefile: free border closest points replaced with corresponding points
269 # on the free border. two polylines, one inside the domain, one outside
271 a = os.path.splitext(shapefileToAdjust)
272 shapefileAdjusted = a[0] + '_adj' + a[1]
273 chainName = os.path.basename(a[0])
275 w = shapefile.Writer(shapefileAdjusted)
280 chaincoords.append(fbs.points[ifb1])
281 for i in range(i1+1, i2):
282 chaincoords.append(sfta.points[i])
283 chaincoords.append(fbs.points[ifb2])
284 w.line([chaincoords])
285 w.record(chainName + '_0')
288 chaincoords.append(fbs.points[ifb2])
289 if i2+1 < len(sfta.points):
290 for i in range(i2+1, len(sfta.points)):
291 chaincoords.append(sfta.points[i])
293 chaincoords.append(sfta.points[i])
294 chaincoords.append(fbs.points[ifb1])
295 w.line([chaincoords])
296 w.record(chainName + '_1')
299 # write the free border splited in two polylines (cut by the adjusted shapefile)
301 a = os.path.splitext(freeBorderShapefile)
302 freeBorderSplit = a[0] + '_split' + a[1]
303 chainName = os.path.basename(a[0])
305 w = shapefile.Writer(freeBorderSplit)
310 i = ifb1; ifb1 = ifb2; ifb2 = i
313 for i in range(ifb1, ifb2+1):
314 chaincoords.append(fbs.points[i])
315 w.line([chaincoords])
316 w.record(chainName + '_0')
319 for i in range(ifb2, len(fbs.points)):
320 chaincoords.append(fbs.points[i])
321 for i in range(ifb1+1):
322 chaincoords.append(fbs.points[i])
323 w.line([chaincoords])
324 w.record(chainName + '_1')