6 from salome.smesh import smeshBuilder
10 import medcoupling as mc
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 nbAdded, MESH, addedBnd = MESH.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [])
23 groups = MESH.GetGroups()
25 aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
26 aCriteria.append(aCriterion)
27 aFilter = smesh.GetFilterFromCriteria(aCriteria)
28 aFilter.SetMesh(MESH.GetMesh())
29 FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter )
30 a = os.path.splitext(meshFile)
31 smesh.SetName(MESH, os.path.basename(a[0]))
32 newMeshName = a[0] + '_brd' + a[1]
35 MESH.ExportMED(newMeshName,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1)
38 print('ExportMED() failed. Invalid file name?')
42 def explodeGroup(grp, grpName):
43 print(" === explodeGroup", grpName)
44 nbCells=grp.getNumberOfCells()
46 dicReverse = {} # id noeud --> id edges
47 for i in range(nbCells):
48 nodcell = grp.getNodeIdsOfCell(i)
49 for j in range(len(nodcell)):
50 if nodcell[j] in dicReverse:
51 dicReverse[nodcell[j]].append(i)
53 dicReverse[nodcell[j]] = [i]
56 usedCells = [False] * nbCells
57 while False in usedCells:
58 icell = usedCells.index(False)
59 usedCells[icell] = True
60 nodcell = grp.getNodeIdsOfCell(icell)
62 chain = [nodcell[0], nodcell[1]]
65 while nextnode in dicReverse:
66 nextcells = dicReverse[nextnode]
67 if len(nextcells) != 2: # end of chain(1) or "edges connector"(>2): stop
72 for i in range(len(nextcells)):
74 if not usedCells[ncell]: # the chain of nodes grows
75 usedCells[ncell] = True
77 nodcell = grp.getNodeIdsOfCell(ncell)
78 if nodcell[0] == nextnode:
79 nextnode = nodcell[1] # forward edge
81 nextnode = nodcell[0] # reversed edge ?
82 chain.append(nextnode)
83 if not newcell: # end of chain, closed
86 while prevnode in dicReverse:
87 prevcells = dicReverse[prevnode]
88 if len(prevcells) != 2: # end of chain(1) or "edges connector"(>2): stop
93 for i in range(len(prevcells)):
95 if not usedCells[ncell]: # the chain of nodes grows
96 usedCells[ncell] = True
98 nodcell = grp.getNodeIdsOfCell(ncell)
99 if nodcell[1] == prevnode:
100 prevnode = nodcell[0] # forward edge
102 prevnode = nodcell[1] # reversed edge ?
103 chain.insert(0, prevnode)
104 if not newcell: # end of chain, closed
108 chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed)
109 nodeChains.append(chainDesc)
114 def writeShapeLines(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
115 print(" === writeShapeLines", grpName)
116 coords = mcMesh.getCoords()
117 w = shapefile.Writer(grpName)
120 for (chain, chainName, closed) in nodeChains:
121 print(" --- ", chainName)
124 coord = coords[node].getValues()
125 coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
126 #print(" ", coordLb93)
127 chaincoords.append(coordLb93)
128 w.line([chaincoords])
133 def writeShapePoints(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
134 print(" === writeShapePoints", grpName)
135 coords = mcMesh.getCoords()
136 w = shapefile.Writer(grpName + '_pts')
139 for (chain, chainName, closed) in nodeChains:
140 print(" --- ", chainName)
143 coord = coords[node].getValues()
144 coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
145 #print(" ", coordLb93)
146 chaincoords.append(coordLb93)
147 w.multipoint(chaincoords)
152 def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.):
153 print(" === exploreEdgeGroups", meshFile)
154 mcMesh = ml.MEDFileMesh.New(meshFile)
155 dim = mcMesh.getSpaceDimension()
156 d1=-1 # when dimension 2, edges are dim -1
157 if dim == 3: # when dimension 3, edges are dim -2
159 a = os.path.splitext(meshFile)
160 prefix = os.path.basename(a[0])
162 grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
163 groups = [mcMesh.getGroup(d1, name) for name in grp_names]
164 for (grp, grpName) in zip(groups, grp_names):
165 fullGrpName = prefix + '_' + grpName
166 nodeChains = explodeGroup(grp, fullGrpName)
167 writeShapeLines(mcMesh, fullGrpName, nodeChains, offsetX, offsetY)
168 writeShapePoints(mcMesh, fullGrpName, nodeChains, offsetX, offsetY)
171 def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust):
172 print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
174 # --- find domain freeBorder: bounding box englobing all others
176 fb = shapefile.Reader(freeBorderShapefile)
177 fbShapes = fb.shapes()
178 maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
180 for i,s in enumerate(fbShapes):
182 if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
185 fbs = fbShapes[outerBboxIndex]
187 # --- find the intersections of the shapefile to adjust and the domain free border:
188 # the closests points (two pairs of points)
190 sf = shapefile.Reader(shapefileToAdjust)
194 for i,p in enumerate(sfta.points):
197 for j,pfb in enumerate(fbs.points):
198 d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
202 pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
203 pdist.sort() # the 2 closest points must be set on the mesh freeBorder
205 i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
206 i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
207 if i1 == pdist[0][2]:
208 ifb1 = pdist[0][1] # index of first adjusted point in freeBorder
209 ifb2 = pdist[1][1] # index of second adjusted point in freeBorder
214 # --- write the adusted shapefile: free border closest points replaced with corresponding points
215 # on the free border. two polylines, one inside the domain, one outside
217 a = os.path.splitext(shapefileToAdjust)
218 shapefileAdjusted = a[0] + '_adj' + a[1]
219 chainName = os.path.basename(a[0])
221 w = shapefile.Writer(shapefileAdjusted)
226 chaincoords.append(fbs.points[ifb1])
227 for i in range(i1+1, i2):
228 chaincoords.append(sfta.points[i])
229 chaincoords.append(fbs.points[ifb2])
230 w.line([chaincoords])
231 w.record(chainName + '_0')
234 chaincoords.append(fbs.points[ifb2])
235 if i2+1 < len(sfta.points):
236 for i in range(i2+1, len(sfta.points)):
237 chaincoords.append(sfta.points[i])
239 chaincoords.append(sfta.points[i])
240 chaincoords.append(fbs.points[ifb1])
241 w.line([chaincoords])
242 w.record(chainName + '_1')
245 # write the free border splited in two polylines (cut by the adjusted shapefile)
247 a = os.path.splitext(freeBorderShapefile)
248 freeBorderSplit = a[0] + '_split' + a[1]
249 chainName = os.path.basename(a[0])
251 w = shapefile.Writer(freeBorderSplit)
256 i = ifb1; ifb1 = ifb2; ifb2 = i
259 for i in range(ifb1, ifb2+1):
260 chaincoords.append(fbs.points[i])
261 w.line([chaincoords])
262 w.record(chainName + '_0')
265 for i in range(ifb2, len(fbs.points)):
266 chaincoords.append(fbs.points[i])
267 for i in range(ifb1+1):
268 chaincoords.append(fbs.points[i])
269 w.line([chaincoords])
270 w.record(chainName + '_1')