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 if not salome.sg.hasDesktop():
27 smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed
28 ([MESH], status) = smesh.CreateMeshesFromMED(meshFileIn)
30 nbAdded, MESH, addedBnd = MESH.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [])
33 aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
34 aCriteria.append(aCriterion)
35 aFilter = smesh.GetFilterFromCriteria(aCriteria)
36 aFilter.SetMesh(MESH.GetMesh())
37 FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter )
39 a = os.path.splitext(meshFileIn)
40 smesh.SetName(MESH, os.path.basename(a[0]))
43 a = os.path.splitext(meshFileIn)
44 smesh.SetName(MESH, os.path.basename(a[0]))
45 meshFileOut = a[0] + '_brd' + a[1]
48 MESH.ExportMED(meshFileOut,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1)
51 print('ExportMED() failed. Invalid file name?')
55 def explodeGroup(grp, grpName):
57 from a group of edges loaded with MEDCoupling, create ordered lists of nodes, one list for each set of connected edges.
59 grp: MEDCoupling object for the group of edges
60 grpName: name of the group
62 List of descriptors [(ordered list of nodeIds, name of the list, closed status)]
64 print(" === explodeGroup", grpName)
65 nbCells=grp.getNumberOfCells()
67 dicReverse = {} # id node --> id edges
68 for i in range(nbCells):
69 nodcell = grp.getNodeIdsOfCell(i)
70 for j in range(len(nodcell)):
71 if nodcell[j] in dicReverse:
72 dicReverse[nodcell[j]].append(i)
74 dicReverse[nodcell[j]] = [i]
77 usedCells = [False] * nbCells
78 while False in usedCells:
79 icell = usedCells.index(False)
80 usedCells[icell] = True
81 nodcell = grp.getNodeIdsOfCell(icell)
83 chain = [nodcell[0], nodcell[1]]
86 while nextnode in dicReverse:
87 nextcells = dicReverse[nextnode]
88 if len(nextcells) != 2: # end of chain(1) or "edges connector"(>2): stop
93 for i in range(len(nextcells)):
95 if not usedCells[ncell]: # the chain of nodes grows
96 usedCells[ncell] = True
98 nodcell = grp.getNodeIdsOfCell(ncell)
99 if nodcell[0] == nextnode:
100 nextnode = nodcell[1] # forward edge
102 nextnode = nodcell[0] # reversed edge ?
103 chain.append(nextnode)
104 if not newcell: # end of chain, closed
107 while prevnode in dicReverse:
108 prevcells = dicReverse[prevnode]
109 if len(prevcells) != 2: # end of chain(1) or "edges connector"(>2): stop
114 for i in range(len(prevcells)):
116 if not usedCells[ncell]: # the chain of nodes grows
117 usedCells[ncell] = True
119 nodcell = grp.getNodeIdsOfCell(ncell)
120 if nodcell[1] == prevnode:
121 prevnode = nodcell[0] # forward edge
123 prevnode = nodcell[1] # reversed edge ?
124 chain.insert(0, prevnode)
125 if not newcell: # end of chain, closed
129 chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed)
130 nodeChains.append(chainDesc)
135 def writeShapeLines(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
137 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.
139 mcMesh: mesh loaded in memory with MEDLoader
140 grpName: name associated to the group of edges
141 nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
142 outputDirectory: directory for writing the shapefile
143 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.
144 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.
146 print(" === writeShapeLines", grpName)
147 coords = mcMesh.getCoords()
148 shapeFileName = os.path.join(outputDirectory, grpName)
149 w = shapefile.Writer(shapeFileName)
152 for (chain, chainName, closed) in nodeChains:
153 print(" --- ", chainName)
156 coord = coords[node].getValues()
157 coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
158 chaincoords.append(coordLb93)
159 w.line([chaincoords])
162 return shapeFileName + '.shp'
165 def writeShapePoints(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
167 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.
169 mcMesh: mesh loaded in memory with MEDLoader
170 grpName: name associated to the group of edges
171 nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
172 outputDirectory: directory for writing the shapefile
173 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.
174 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.
176 print(" === writeShapePoints", grpName)
177 coords = mcMesh.getCoords()
178 shapeFileName = os.path.join(outputDirectory, grpName)
179 w = shapefile.Writer(shapeFileName + '_pts')
182 for (chain, chainName, closed) in nodeChains:
183 print(" --- ", chainName)
186 coord = coords[node].getValues()
187 coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
188 chaincoords.append(coordLb93)
189 w.multipoint(chaincoords)
194 def exploreEdgeGroups(meshFile, outputDirectory="", offsetX=0., offsetY=0.):
196 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)
198 meshFile: full path of the input mesh file, format MED
199 outputDirectory: directory in which the shapefiles are written (default "", if "" use the directory containing the mesh
200 offsetX: local X origin of the mesh
201 offsetY: local Y origin of the mesh
203 print(" === exploreEdgeGroups", meshFile)
204 if outputDirectory == "":
205 outputDirectory = os.path.dirname(meshFile)
207 a = os.path.splitext(meshFile)
208 prefix = os.path.basename(a[0]) # prefix = file name without extension
210 mcMesh = ml.MEDFileMesh.New(meshFile)
211 dim = mcMesh.getMeshDimension()
212 d1=-1 # when dimension 2, edges are dim -1
213 if dim == 3: # when dimension 3, edges are dim -2
216 grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
218 groups = [mcMesh.getGroup(d1, name) for name in grp_names] # list of groups in their name order
221 for (grp, grpName) in zip(groups, grp_names):
222 fullGrpName = prefix + '_' + grpName
223 nodeChains = explodeGroup(grp, fullGrpName)
224 filename = writeShapeLines(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
225 writeShapePoints(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
226 filenames.append(filename)
227 shapesListFile = os.path.join(outputDirectory, "shapesList.json")
228 with open(shapesListFile, 'w') as f:
229 json.dump(filenames, f)
231 def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust, outputDirectory="", splitFreeBorder=True, splitShapeAdjusted=True):
233 shapeFileToAdjust must be a closed line or polygon crossing freeBorderShapefile in 2 points.
234 Find in shapeFileToAdjust and in freeBorderShapefile the two closest corresponding points and move the points in shapeFileToAdjust to correspond to the points found in freeBorderShapefile. Split shapeFileToAdjust in two parts (inside or outside freeBorder). If requested, split freeBorderShapefile in two parts. Same for shapeFileToAdjust.
236 freeBorderShapefile: a set of free border lines, as generated by the functions freeBordersGroup and exploreEdgeGroups.
237 shapefileToAdjust: a closed line or polygon, supposed to be drawn in qgis to pass as close as possible to the points to be connected, on the free border.
238 outputDirectory: if empty, write the resulting shapefiles in their respective directory, with the suffix '_adj', otherwise write in the outputDirectory.
239 splitFreeBorder: boolean default True
240 splitShapeAdjusted: boolean default True
242 print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
244 # --- find domain freeBorder: bounding box englobing all others
245 # TODO: This may not be always the case, when there is not a single domain with holes, but several non connected domains.
247 fb = shapefile.Reader(freeBorderShapefile)
248 fbShapes = fb.shapes()
249 maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
251 for i,s in enumerate(fbShapes):
253 if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
256 fbs = fbShapes[outerBboxIndex] # the domain free border shape
258 # --- find the intersections of the shapefile to adjust and the domain free border:
259 # the closests points (two pairs of points)
261 sf = shapefile.Reader(shapefileToAdjust)
267 discountLastSftaPoint = 0
268 for i,p in enumerate(sfta.points):
270 x00 = p[0] # keep first point
272 d = math.sqrt((x0-p[0])*(x0-p[0]) + (y0-p[1])*(y0-p[1])) # distance to previous point
273 x0 = p[0] # keep previous point
276 print("two consecutives points of shapefile To adjust are superposed, OK")
277 continue # do not take into account consecutive superposed points in shape to adjust
278 if i == len(sfta.points) -1:
279 d = math.sqrt((x00-p[0])*(x00-p[0]) + (y00-p[1])*(y00-p[1])) # distance between first and last point
281 discountLastSftaPoint = 1
282 print("last point of shapefile To adjust is superposed to first point, last point discarded, OK")
283 continue # do not take into account last point if superposed to first point, in shape to adjust
286 for j,pfb in enumerate(fbs.points):
287 d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
291 pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
292 pdist.sort() # the 2 closest points must be set on the mesh freeBorder
294 i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
295 i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
296 if i1 == pdist[0][2]:
297 ifb1 = pdist[0][1] # index of first adjusted point in freeBorder
298 ifb2 = pdist[1][1] # index of second adjusted point in freeBorder
302 print("i1, i2, len(sfta.points)", i1, i2, len(sfta.points))
303 print("ifb1, ifb2, len(fbs.points)", ifb1, ifb2, len(fbs.points))
305 # --- write the adusted shapefile: free border closest points replaced with corresponding points
306 # on the free border. two polylines, one inside the domain, one outside
308 if outputDirectory == "":
309 outputDirectory = os.path.dirname(shapefileToAdjust)
310 a = os.path.splitext(os.path.basename(shapefileToAdjust))
311 shapefileAdjusted = os.path.join(outputDirectory, a[0] + '_adj' + a[1])
314 w = shapefile.Writer(shapefileAdjusted)
318 if splitShapeAdjusted:
320 chaincoords.append(fbs.points[ifb1])
321 for i in range(i1+1, i2):
322 chaincoords.append(sfta.points[i])
323 chaincoords.append(fbs.points[ifb2])
324 w.line([chaincoords])
325 w.record(chainName + '_0')
328 chaincoords.append(fbs.points[ifb2])
329 if i2+1 < len(sfta.points):
330 for i in range(i2+1, len(sfta.points) -discountLastSftaPoint):
331 chaincoords.append(sfta.points[i])
333 chaincoords.append(sfta.points[i])
334 chaincoords.append(fbs.points[ifb1])
335 w.line([chaincoords])
336 w.record(chainName + '_1')
339 chaincoords.append(fbs.points[ifb1])
340 for i in range(i1+1, i2):
341 chaincoords.append(sfta.points[i])
342 chaincoords.append(fbs.points[ifb2])
343 if i2+1 < len(sfta.points):
344 for i in range(i2+1, len(sfta.points) -discountLastSftaPoint):
345 chaincoords.append(sfta.points[i])
347 chaincoords.append(sfta.points[i])
348 if discountLastSftaPoint:
349 chaincoords.append(fbs.points[ifb1]) # close shape when first point if superposed with last
350 w.line([chaincoords])
356 # write the free border split in two polylines (cut by the adjusted shapefile)
358 if outputDirectory == "":
359 outputDirectory = os.path.dirname(freeBorderShapefile)
360 a = os.path.splitext(os.path.basename(freeBorderShapefile))
361 freeBorderSplit = os.path.join(outputDirectory, a[0] + '_split' + a[1])
364 w = shapefile.Writer(freeBorderSplit)
369 i = ifb1; ifb1 = ifb2; ifb2 = i
372 for i in range(ifb1, ifb2+1):
373 chaincoords.append(fbs.points[i])
374 w.line([chaincoords])
375 w.record(chainName + '_0')
378 for i in range(ifb2, len(fbs.points)):
379 chaincoords.append(fbs.points[i])
380 for i in range(ifb1+1):
381 chaincoords.append(fbs.points[i])
382 w.line([chaincoords])
383 w.record(chainName + '_1')