]> SALOME platform Git repositories - tools/medcoupling.git/blob - resources/dev/my_findNodesToDup.py
Salome HOME
Increment version: 9.12.0
[tools/medcoupling.git] / resources / dev / my_findNodesToDup.py
1 """ Python version of MEDCouplingUMesh::findNodesToDuplicate() and MEDCouplingUMesh::findCellsToRenumber() methods which are at the core of the 
2     MEDFileUMesh::buildInnerBoundaryAlongM1Group() algorithm.
3     This greatly helps algorithm tuning ...
4 """
5
6 from medcoupling import *
7
8 def findNodesToDuplicate(this, otherDimM1OnSameCoords):
9   # Checking star-shaped M1 group:
10   meshM2, _,_,_,rdit0 = otherDimM1OnSameCoords.buildDescendingConnectivity() # 2D: a mesh of points, 3D: a mesh of segs
11   dsi = rdit0.deltaShiftIndex()
12   idsTmp0 = dsi.findIdsNotInRange(-1, 3)  # for 2D: if a point is connected to more than 2 segs. For 3D: if a seg is connected to more than two faces.
13   if(idsTmp0.getNumberOfTuples()):
14     raise ValueError("")
15
16   # Get extreme nodes from the group (they won't be duplicated except if they also lie on bound of M0 -- see below),
17   # ie nodes belonging to the boundary "cells" (might be points) of M1
18   xtremIdsM2 = dsi.findIdsEqual(1)
19   meshM2Part = meshM2[xtremIdsM2]
20   xtrem = meshM2Part.computeFetchedNodeIds()
21   # Remove from the list points on the boundary of the M0 mesh (those need duplication!)
22   m0desc, dt0, dit0, rdt0, rdit0 = this.buildDescendingConnectivity()
23   dsi = rdit0.deltaShiftIndex()
24   boundSegs = dsi.findIdsEqual(1)  # boundary segs/faces of the M0 mesh
25   m0descSkin = m0desc[boundSegs]
26   fNodes = m0descSkin.computeFetchedNodeIds()  # fNodes needs dupl
27   # In 3D, some points on the boundary of M0 will NOT be duplicated (where as in 2D, points on the boundary of M0 are always duplicated)
28   # Think of a partial (plane) crack in a cube: the points at the tip of the crack and not located inside the volume of the cube are not duplicated
29   # although they are technically on the skin of the cube.
30   if this.getMeshDimension() == 3 :
31       m0descSkinDesc, _, _, _, _ = m0descSkin.buildDescendingConnectivity() # all segments of the skin of the 3D (M0) mesh
32       _, corresp = meshM2.areCellsIncludedIn(m0descSkinDesc,2)
33       # validIds is the list of segments which are on both the skin of *this*, and in the segments of the M1 group
34       # In the cube example above, this is a U shape polyline.
35       validIds = corresp.findIdsInRange(0, meshM2.getNumberOfCells())
36       if validIds.getNumberOfTuples():
37           # Build the set of segments which are: in the desc mesh of the skin of the 3D mesh (M0) **and** in the desc mesh of the M1 group:
38           # (the U-shaped polyline described above)
39           m1IntersecSkin = m0descSkinDesc[validIds]
40           # Its boundary nodes should no be duplicated (this is for example the tip of the crack inside the cube described above)
41           notDuplSkin = m1IntersecSkin.findBoundaryNodes()
42           fNodes1 = fNodes.buildSubstraction(notDuplSkin) # fNodes1 needs dupl
43
44           # Specific logic to handle singular points :
45           #   - a point on this U-shape line used in a cell which has no face in common with M1 is deemed singular.
46           #   - indeed, if duplicated, such a point would lead to the duplication of a cell which has no face touching M1 ! The
47           #   algorithm would be duplicating too much ...
48           # This is a costly algorithm so only go into it if a simple (non sufficient) criteria is met: a node connected to more than 3 segs in meshM2:
49           meshM2Desc, _, _, _, rdit0 = meshM2.buildDescendingConnectivity()  # a mesh made of node cells
50           dsi = rdit0.deltaShiftIndex()
51           singPoints = dsi.findIdsNotInRange(-1,4)   # points connected to (strictly) more than 3 segments
52           if singPoints.getNumberOfTuples():
53               print ("Hitting singular point logic")
54               boundNodes = m1IntersecSkin.computeFetchedNodeIds()
55               # If a point on this U-shape line is connected to cells which do not share any face with M1, then it
56               # should not be duplicated
57               #    1. Extract N D cells touching U-shape line:
58               cellsAroundBN = this.getCellIdsLyingOnNodes(boundNodes, False)  # false= take cell in, even if not all nodes are in dupl
59               mAroundBN = this[cellsAroundBN]
60               mAroundBNDesc, descBN,descIBN,revDescBN,revDescIBN=mAroundBN.buildDescendingConnectivity()
61               #    2. Identify cells in sub-mesh mAroundBN which have a face in common with M1
62               _, idsOfM1BN = mAroundBNDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2)
63               nCells, nCellsDesc = mAroundBN.getNumberOfCells(), mAroundBNDesc.getNumberOfCells()
64               idsTouch = DataArrayInt.New(); idsTouch.alloc(0,1)
65               for v in idsOfM1BN:
66                   if v[0] >= nCellsDesc:    # Keep valid match only
67                       continue
68                   idx0 = revDescIBN[v[0], 0]
69                   c1, c2 = revDescBN[idx0, 0], revDescBN[idx0+1,0]
70                   idsTouch.pushBackSilent(c1)
71                   idsTouch.pushBackSilent(c2)
72               #    3. Build complement
73               idsTouchCompl = idsTouch.buildComplement(nCells)
74               mAroundBNStrict = mAroundBN[idsTouchCompl]
75               nod3 = mAroundBNStrict.computeFetchedNodeIds()
76               inters = boundNodes.buildIntersection(nod3)
77               print("sing,", inters.getValues())
78               fNodes1 = fNodes1.buildSubstraction(inters)  # reminder: fNodes1 represent nodes that need dupl.
79           notDup = xtrem.buildSubstraction(fNodes1)
80       else:  # if validIds ...
81         notDup = xtrem.buildSubstraction(fNodes)
82   else:  # if 3D ...
83     notDup = xtrem.buildSubstraction(fNodes)
84
85   m1Nodes = otherDimM1OnSameCoords.computeFetchedNodeIds()
86   dupl = m1Nodes.buildSubstraction(notDup)
87   return dupl
88
89
90 def findCellsToRenumber(this, otherDimM1OnSameCoords, dupl):
91   """  Find cells to renumber
92   """
93   # All N D cells touching our group (even when this is just one point touching)
94   cellsAroundGroupLarge = this.getCellIdsLyingOnNodes(dupl, False)  # false= take cell in, even if not all nodes are in dupl
95   #
96   mAroundGrpLarge=this[cellsAroundGroupLarge]
97   mArGrpLargeDesc,descL,descIL,revDescL,revDescIL=mAroundGrpLarge.buildDescendingConnectivity()
98   mAroundGrpLarge.writeVTK("/tmp/mAr_large.vtu")
99   mArGrpLargeDesc.writeVTK("/tmp/mAr_large_desc.vtu")
100
101   # Extract now all N D cells which have a complete face in touch with the group:
102   #    1. Identify cells of M1 group in sub-mesh mAroundGrp
103   _, idsOfM1Large = mArGrpLargeDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2)
104   nL = mArGrpLargeDesc.getNumberOfCells()
105   idsStrict = DataArrayInt.New(); idsStrict.alloc(0,1)
106   #    2. Build map giving for each cell ID in mAroundGrp (not in mAroundGrpLarge) the corresponding cell
107   #       ID on the other side of the crack:
108   toOtherSide, pos = {}, {}
109   cnt = 0
110   for v in idsOfM1Large:
111       if v[0] >= nL:    # Keep valid match only
112           continue
113       idx0 = revDescIL[v[0], 0]
114       # Keep the two cells on either side of the face v of M1:
115       c1, c2 = revDescL[idx0, 0], revDescL[idx0+1,0]
116       if not c1 in idsStrict:
117           pos[c1] = cnt
118           idsStrict.pushBackSilent(c1)
119           cnt += 1
120       if not c2 in idsStrict:
121           pos[c2] = cnt
122           idsStrict.pushBackSilent(c2)
123           cnt += 1
124       k1, k2 = pos[c1], pos[c2]
125       toOtherSide[k1] = k2
126       toOtherSide[k2] = k1
127
128   cellsAroundGroup = cellsAroundGroupLarge[idsStrict]
129   mAroundGrp = this[cellsAroundGroup]
130   nCells, nCellsLarge = cellsAroundGroup.getNumberOfTuples(), cellsAroundGroupLarge.getNumberOfTuples()
131   mArGrpDesc,desc,descI,revDesc,revDescI=mAroundGrp.buildDescendingConnectivity()
132   _, idsOfM1 = mArGrpDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2) # TODO : could we avoid recomputing this??
133   mAroundGrp.writeVTK("/tmp/mAr.vtu")
134   mArGrpDesc.writeVTK("/tmp/mAr_desc.vtu")
135
136   # Neighbor information of the mesh WITH the crack (some neighbors are removed):
137   #     In the neighbor information remove the connection between high dimension cells and its low level constituents which are part
138   #     of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack):
139   DataArrayInt.RemoveIdsFromIndexedArrays(idsOfM1,desc,descI)
140   #     Compute the neighbor of each cell in mAroundGrp, taking into account the broken link above. Two
141   #     cells on either side of the crack (defined by the mesh of low dimension) are not neighbor anymore.
142   neigh, neighI = MEDCouplingUMesh.ComputeNeighborsOfCellsAdv(desc,descI,revDesc,revDescI)
143
144   # For each initial connex part of the M1 mesh (or said differently for each independent crack):
145   seed, nIter, cnt = 0, 0, 0
146   nIterMax = nCells+1 # Safety net for the loop
147   hitCells = DataArrayInt.New(); hitCells.alloc(nCells)
148   hitCells.fillWithValue(0)  # 0 : not hit, -1: one side of the crack, +1: other side of the crack
149   MAX_CP = 10000 # the choices below assume we won't have more than 10000 different connex parts ...
150   PING_FULL_init, PING_PART =  0, MAX_CP
151   PONG_FULL_init, PONG_PART = -0,-MAX_CP
152   while nIter < nIterMax:
153 #       print("dbg ", hitCells.getValues())
154       t = hitCells.findIdsEqual(0)
155       if not t.getNumberOfTuples():
156         break
157       seed = t[0,0]
158       done = False
159       cnt += 1
160       PING_FULL = PING_FULL_init+cnt
161       PONG_FULL = PONG_FULL_init-cnt
162       while not done and nIter < nIterMax:  # Start of the ping-pong
163           nIter += 1
164           # Identify connex zone around the seed
165           spreadZone, _ = MEDCouplingUMesh.ComputeSpreadZoneGraduallyFromSeed([seed],  neigh,neighI, -1)
166           done = True
167           for i, s in enumerate(spreadZone.getValues()):
168               hitCells[s] = PING_FULL
169               if s in toOtherSide:
170                   other = toOtherSide[s]
171                   if hitCells[other] != PONG_FULL:
172                       done = False
173                       hitCells[other] = PONG_PART
174                       #  Compute next seed, i.e. a cell on the other side of the crack
175                       seed = other
176           if done:
177               # we might have several disjoing PONG parts in front of a single PING connex part:
178               idsPong = hitCells.findIdsEqual(PONG_PART)
179               if idsPong.getNumberOfTuples():
180                   seed = idsPong[0,0]
181                   done = False
182               continue  # continue without switching side (or break if done remains false)
183           else:
184               # Go the other side
185               PING_FULL, PONG_FULL = PONG_FULL, PING_FULL
186               PING_PART, PONG_PART = PONG_PART, PING_PART
187
188       nonHitCells = hitCells.findIdsEqual(0)
189       if nonHitCells.getNumberOfTuples():
190         seed = nonHitCells[0,0]
191       else:
192         break
193
194   if nIter >= nIterMax:
195     raise ValueError("Too many iterations - should not happen")
196
197   # Now we have handled all N D cells which have a face touching the M1 group. It remains the cells
198   # which are just touching the group by one (or several) node(s):
199   # All those cells are in direct contact with a cell which is either PING_FULL or PONG_FULL
200   # So first reproject the PING/PONG info onto mAroundGrpLarge:
201   hitCellsLarge = DataArrayInt.New(); hitCellsLarge.alloc(nCellsLarge)
202   hitCellsLarge.fillWithValue(0)
203   hitCellsLarge[idsStrict] = hitCells
204   nonHitCells = hitCellsLarge.findIdsEqual(0)
205   # Neighbor information in mAroundGrpLarge:
206   neighL, neighIL = MEDCouplingUMesh.ComputeNeighborsOfCellsAdv(descL,descIL,revDescL,revDescIL)
207   for c in nonHitCells:
208       assert(False)
209       neighs = neighL[neighIL[c[0]]:neighIL[c[0]+1]]
210       for n in neighs:
211           neighVal = hitCellsLarge[n[0]]
212           if neighVal != 0 and abs(neighVal) < MAX_CP:  # (@test_T0) second part of the test to skip cells being assigned and target only cells assigned in the first part of the algo above
213               currVal = hitCellsLarge[c[0]]
214               if currVal != 0:   # Several neighbors have a candidate number
215                   # Unfortunately in some weird cases (see testBuildInnerBoundary8) a cell in mAroundGrpLarge
216                   # might have as neighbor two conflicting spread zone ...
217                   if currVal*neighVal < 0:
218                       # If we arrive here, the cell was already assigned a number and we found a neighbor with
219                       # a different sign ... we must swap the whole spread zone!!
220                       print("Ouch - must switch spread zones ...")
221                       ids1 = hitCellsLarge.findIdsEqual(neighVal)
222                       ids1b = hitCellsLarge.findIdsEqual(-neighVal)
223                       ids2 = hitCellsLarge.findIdsEqual(MAX_CP*neighVal)
224                       ids2b = hitCellsLarge.findIdsEqual(-MAX_CP*neighVal)
225                       hitCellsLarge[ids1] *= -1
226                       hitCellsLarge[ids1b] *= -1
227                       hitCellsLarge[ids2] *= -1
228                       hitCellsLarge[ids2b] *= -1
229               else:  # First assignation
230                   hitCellsLarge[c[0],0] = MAX_CP*neighVal   # Same sign, but different value to preserve PING_FULL and PONG_FULL
231
232 ###
233 ###  SO FAR THE LOGIC BELOW WAS NOT NEEDED ....
234 ###
235
236 #   # Now handling remaining cells not touched by the above process, called "naked" cells (see cell #20 in mArndLarge in testBuildInnerBoundary8() ...)
237 #   naked = hitCellsLarge.findIdsEqual(0)
238 #   mLargC, mLargCI = mArGrpLargeDesc.getNodalConnectivity(),mArGrpLargeDesc.getNodalConnectivityIndex()
239 #   for c in naked:
240 #       neighs = neighL[neighIL[c[0]]:neighIL[c[0]+1]]  # ExtractFromIndexedArray?
241 #       nbDup = {}
242 #       fac1 = descL[descIL[c[0]]:descIL[c[0]+1]]
243 #       for n in neighs:
244 #           if hitCellsLarge[n[0]] == 0:
245 #               continue   # this neighbour is naked too, nothing we can do for now
246 #           # Among the values found on neighbour cells, take the one from the neighbour which is connected
247 #           # with the most "economical" face, i.e. the face made of a minimal number of duplicated points.
248 #           # TODO: this is a shaky criteria ... find sth more robust ...
249 #           #   1. find face(s) making the link
250 #           fac2 = descL[descIL[n[0]]:descIL[n[0]+1]]
251 #           com = fac1.buildIntersection(fac2)
252 #           if (com.getNumberOfTuples() == 0):
253 #               raise ValueError("Internal error : no common face ?")
254 #           #   2. count number of duplicated node for this face.
255 #           for f in com: # for all common faces
256 #               faceNodes = mLargC[mLargCI[f[0]]+1:mLargCI[f[0]+1]] # first +1 to skip type
257 #               comNod = faceNodes.buildIntersection(dupl)
258 #               # in case the two cells are in contact by multiple faces, take the most conservative value
259 #               nbDup[n[0]] = max(nbDup.get(n[0],-1), comNod.getNumberOfTuples())
260 #       # Minimal value in nbDup?
261 #       cellIdx = min(nbDup, key=nbDup.get)
262 #       hitCellsLarge[c[0]] = hitCellsLarge[cellIdx]
263 #
264 #   cellsToModifyConn0_torenum = hitCellsLarge.findIdsInRange(1,MAX_CP)    # Positive spread zone number
265 #   cellsToModifyConn1_torenum = hitCellsLarge.findIdsInRange(-MAX_CP, 0)  # Negative spread zone number
266
267
268 ###
269 ###  C++ VERSION OF IT
270 ###
271 # //
272 # //  // Now handling remaining cells not touched by the for loop above, called "naked" cells (see cell #20 in mArndGrpLarge in testBuildInnerBoundary8() ...)
273 # //  DAInt naked = hitCellsLarge->findIdsEqual(0);
274 # //  const mcIdType *mLargCP=mArGrpLargeDesc->getNodalConnectivity()->begin(), *mLargCIP=mArGrpLargeDesc->getNodalConnectivityIndex()->begin();
275 # //  for (const auto &c: *naked)
276 # //    {
277 # //      std::map<mcIdType, mcIdType> nbDup;
278 # //      // Retrieve list of faces of cell c
279 # //      mcIdType nbFac1=descILP[c+1]-descILP[c];
280 # //      std::vector<mcIdType> fac1(nbFac1);
281 # //      std::copy(descLP+descILP[c], descLP+descILP[c+1], fac1.begin());
282 # //      std::sort(fac1.begin(), fac1.end());
283 # //      mcIdType cnt00 = neighILP[c];
284 # //      for (const mcIdType *n=neighLP+cnt00; cnt00 < neighILP[c+1]; n++, cnt00++)
285 # //        {
286 # //          if (hitCellsLargeP[*n] == 0)
287 # //            continue;   // this neighbour is naked too, nothing we can do for now
288 # //          // Among the values found on neighbour cells, take the one from the neighbour which is connected
289 # //          // with the most "economical" face, i.e. the face made of a minimal number of duplicated points.
290 # //          // TODO: this is a shaky criteria ... find sth more robust ...
291 # //          //   1. find face(s) making the link
292 # //          mcIdType nbFac2=descILP[*n+1]-descILP[*n];
293 # //          std::vector<mcIdType> fac2(nbFac2);
294 # //          std::copy(descLP+descILP[*n], descLP+descILP[*n+1], fac2.begin());
295 # //          std::sort(fac2.begin(), fac2.end());
296 # //          std::vector<mcIdType> comFac;
297 # //          std::set_intersection(fac1.begin(), fac1.end(),
298 # //                                fac2.begin() ,fac2.end(),
299 # //                                std::back_inserter(comFac));
300 # //          if (comFac.size() == 0)
301 # //            throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: internal error no common face between two cells should not happen");
302 # //          //   2. count number of duplicated node for this face.
303 # //          for (const auto &f : comFac) // for all common faces
304 # //            {
305 # //              std::vector<mcIdType> comNod;
306 # //              std::set_intersection(nodeIdsToDuplicateBg, nodeIdsToDuplicateEnd,
307 # //                                    mLargCP+mLargCIP[f]+1, mLargCP+mLargCIP[f+1],    // first +1 to skip type in connectivity
308 # //                                    std::back_inserter(comNod));
309 # //              // in case the two cells are in contact by multiple faces, take the most conservative value
310 # //              mcIdType val=-1;
311 # //              if(nbDup.find(*n) != nbDup.end()) val=nbDup[*n];
312 # //              nbDup[*n] = std::max(val, (mcIdType)comNod.size());
313 # //            }
314 # //        }
315 # //      // Minimal value in nbDup?
316 # //      using PairId = std::pair<mcIdType, mcIdType>;
317 # //      auto comp_fonc = [](const PairId& p1, const PairId& p2) { return p1.second < p2.second; };
318 # //      PairId zemin = *min_element(nbDup.begin(), nbDup.end(), comp_fonc);
319 # //      hitCellsLargeP[c] = hitCellsLargeP[zemin.first];
320 # //    }
321
322
323   cellsToModifyConn0_torenum = hitCellsLarge.findIdsInRange(1,MAX_CP*MAX_CP)    # Positive spread zone number
324   cellsToModifyConn1_torenum = hitCellsLarge.findIdsInRange(-MAX_CP*MAX_CP, 0)  # Negative spread zone number
325   if cellsToModifyConn0_torenum.getNumberOfTuples() + cellsToModifyConn1_torenum.getNumberOfTuples() != cellsAroundGroupLarge.getNumberOfTuples():
326       raise ValueError("Some cells not hit - Internal error should not happen")
327   cellsToModifyConn0_torenum.transformWithIndArr(cellsAroundGroupLarge)
328   cellsToModifyConn1_torenum.transformWithIndArr(cellsAroundGroupLarge)
329   #
330   cellIdsNeededToBeRenum=cellsToModifyConn0_torenum
331   cellIdsNotModified=cellsToModifyConn1_torenum
332
333   return cellIdsNeededToBeRenum, cellIdsNotModified
334