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