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, splitFreeBorder=True):
227 shapeFileToAdjust must be a closed line or polygon crossing freeBorderShapefile in 2 points.
228 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.
230 freeBorderShapefile: a set of free border lines, as generated by the functions freeBordersGroup and exploreEdgeGroups.
231 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.
232 splitFreeBorder: boolean default True
234 print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
236 # --- find domain freeBorder: bounding box englobing all others
237 # TODO: This may not be always the case, when there is not a single domain with holes, but several non connected domains.
239 fb = shapefile.Reader(freeBorderShapefile)
240 fbShapes = fb.shapes()
241 maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
243 for i,s in enumerate(fbShapes):
245 if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
248 fbs = fbShapes[outerBboxIndex] # the domain free border shape
250 # --- find the intersections of the shapefile to adjust and the domain free border:
251 # the closests points (two pairs of points)
253 sf = shapefile.Reader(shapefileToAdjust)
257 for i,p in enumerate(sfta.points):
260 for j,pfb in enumerate(fbs.points):
261 d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
265 pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
266 pdist.sort() # the 2 closest points must be set on the mesh freeBorder
268 i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
269 i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
270 if i1 == pdist[0][2]:
271 ifb1 = pdist[0][1] # index of first adjusted point in freeBorder
272 ifb2 = pdist[1][1] # index of second adjusted point in freeBorder
277 # --- write the adusted shapefile: free border closest points replaced with corresponding points
278 # on the free border. two polylines, one inside the domain, one outside
280 a = os.path.splitext(shapefileToAdjust)
281 shapefileAdjusted = a[0] + '_adj' + a[1]
282 chainName = os.path.basename(a[0])
284 w = shapefile.Writer(shapefileAdjusted)
289 chaincoords.append(fbs.points[ifb1])
290 for i in range(i1+1, i2):
291 chaincoords.append(sfta.points[i])
292 chaincoords.append(fbs.points[ifb2])
293 w.line([chaincoords])
294 w.record(chainName + '_0')
297 chaincoords.append(fbs.points[ifb2])
298 if i2+1 < len(sfta.points):
299 for i in range(i2+1, len(sfta.points)):
300 chaincoords.append(sfta.points[i])
302 chaincoords.append(sfta.points[i])
303 chaincoords.append(fbs.points[ifb1])
304 w.line([chaincoords])
305 w.record(chainName + '_1')
309 # write the free border split in two polylines (cut by the adjusted shapefile)
311 a = os.path.splitext(freeBorderShapefile)
312 freeBorderSplit = a[0] + '_split' + a[1]
313 chainName = os.path.basename(a[0])
315 w = shapefile.Writer(freeBorderSplit)
320 i = ifb1; ifb1 = ifb2; ifb2 = i
323 for i in range(ifb1, ifb2+1):
324 chaincoords.append(fbs.points[i])
325 w.line([chaincoords])
326 w.record(chainName + '_0')
329 for i in range(ifb2, len(fbs.points)):
330 chaincoords.append(fbs.points[i])
331 for i in range(ifb1+1):
332 chaincoords.append(fbs.points[i])
333 w.line([chaincoords])
334 w.record(chainName + '_1')