]> SALOME platform Git repositories - tools/medcoupling.git/commitdiff
Salome HOME
Handling another funny case.
authorabn <adrien.bruneton@cea.fr>
Fri, 22 Oct 2021 11:49:22 +0000 (13:49 +0200)
committerabn <adrien.bruneton@cea.fr>
Fri, 22 Oct 2021 11:49:22 +0000 (13:49 +0200)
resources/dev/my_findNodesToDup.py
src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDLoader/Swig/MEDLoaderTest3.py

index 6ba2b41eaedbfc643285ab12f6037374f4f15bca..cbc53afcd800cdd93da9074a9dade1bc5c82bcf8 100644 (file)
@@ -72,7 +72,6 @@ def findNodesToDuplicate(this, otherDimM1OnSameCoords):
                     break
               if nIter >= nIterMax:
                 raise ValueError("Internal error should not happen")
-              print("dbg", hitCells.getValues())
               if hitCells.buildUniqueNotSorted().getNumberOfTuples() > 1:
                   # All points common to two (or more) different spread zones are singular:
                   cc, ccI = m1Bound.getReverseNodalConnectivity()
@@ -99,60 +98,79 @@ def findNodesToDuplicate(this, otherDimM1OnSameCoords):
 def findCellsToRenumber(this, otherDimM1OnSameCoords, dupl):
   """  Find cells to renumber
   """
-  cellsAroundGroup = this.getCellIdsLyingOnNodes(dupl, False)  # false= take cell in, even if not all nodes are in dupl
-
+  # All N D cells touching our group (even when this is just one point touching)
+  cellsAroundGroupLarge = this.getCellIdsLyingOnNodes(dupl, False)  # false= take cell in, even if not all nodes are in dupl
   #
-  mAroundGrp=this[cellsAroundGroup]
-  nCells2 = mAroundGrp.getNumberOfCells()
-  mArGrpDesc,desc00,descI00,revDesc00,revDescI00=mAroundGrp.buildDescendingConnectivity()
+  mAroundGrpLarge=this[cellsAroundGroupLarge]
+  mArGrpLargeDesc,descL,descIL,revDescL,revDescIL=mAroundGrpLarge.buildDescendingConnectivity()
+  mAroundGrpLarge.writeVTK("/tmp/mAr_large.vtu")
+  mArGrpLargeDesc.writeVTK("/tmp/mAr_large_desc.vtu")
+
+  # Extract now all N D cells which have a complete face in touch with the group: 
+  #    1. Identify cells of M1 group in sub-mesh mAroundGrp
+  _, idsOfM1Large = mArGrpLargeDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2)
+  nL = mArGrpLargeDesc.getNumberOfCells()
+  idsStrict = DataArrayInt.New(); idsStrict.alloc(0,1)
+  #    2. Build map giving for each cell ID in mAroundGrp (not in mAroundGrpLarge) the corresponding cell 
+  #       ID on the other side of the crack:
+  toOtherSide, pos = {}, {}
+  cnt = 0
+  for v in idsOfM1Large:
+      if v[0] >= nL:    # Keep valid match only
+          continue
+      idx0 = revDescIL[v[0], 0]
+      # Keep the two cells on either side of the face v of M1:
+      c1, c2 = revDescL[idx0, 0], revDescL[idx0+1,0]
+      if not c1 in idsStrict:
+          pos[c1] = cnt
+          idsStrict.pushBackSilent(c1)
+          cnt += 1
+      if not c2 in idsStrict:
+          pos[c2] = cnt
+          idsStrict.pushBackSilent(c2)
+          cnt += 1
+      k1, k2 = pos[c1], pos[c2]
+      toOtherSide[k1] = k2
+      toOtherSide[k2] = k1
 
+  cellsAroundGroup = cellsAroundGroupLarge[idsStrict]
+  mAroundGrp = this[cellsAroundGroup]
+  nCells, nCellsLarge = cellsAroundGroup.getNumberOfTuples(), cellsAroundGroupLarge.getNumberOfTuples()
+  mArGrpDesc,desc,descI,revDesc,revDescI=mAroundGrp.buildDescendingConnectivity()
+  _, idsOfM1 = mArGrpDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2) # TODO : could we avoid recomputing this??
   mAroundGrp.writeVTK("/tmp/mAr.vtu")
   mArGrpDesc.writeVTK("/tmp/mAr_desc.vtu")
 
-  # Identify cells of M1 group in sub-mesh mAroundGrp
-  _, idsOfM1 = mArGrpDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2)
-  nCellsArGrpDesc = mArGrpDesc.getNumberOfCells()
-#   print(idsOfM1.getValues())
-
-  # Build map giving for each cell ID in mAroundGrp the corresponding cell ID on the other side of the crack:
-  # Note that this does not cover all cells around the crack (a cell like a triangle might touche the crack only with its tip)
-  toOtherSide = {}
-  for i, v in enumerate(idsOfM1):
-    if v[0] >= nCellsArGrpDesc:    # Keep valid match only
-        continue
-    idx0 = revDescI00[v[0], 0]
-    c1, c2 = revDesc00[idx0, 0], revDesc00[idx0+1,0]
-    toOtherSide[c1] = c2
-    toOtherSide[c2] = c1
-
   # Neighbor information of the mesh WITH the crack (some neighbors are removed):
   #     In the neighbor information remove the connection between high dimension cells and its low level constituents which are part
   #     of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack):
-  DataArrayInt.RemoveIdsFromIndexedArrays(idsOfM1,desc00,descI00)
+  DataArrayInt.RemoveIdsFromIndexedArrays(idsOfM1,desc,descI)
   #     Compute the neighbor of each cell in mAroundGrp, taking into account the broken link above. Two
   #     cells on either side of the crack (defined by the mesh of low dimension) are not neighbor anymore.
-  neigh00, neighI00 = MEDCouplingUMesh.ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00)
+  neigh, neighI = MEDCouplingUMesh.ComputeNeighborsOfCellsAdv(desc,descI,revDesc,revDescI)
 
   # For each initial connex part of the M1 mesh (or said differently for each independent crack):
-  seed, nIter = 0, 0
-  nIterMax = nCells2+1 # Safety net for the loop
-  hitCells = DataArrayInt.New(); hitCells.alloc(nCells2)
-  hitCells.fillWithValue(-1)  # -1 : not hit, 0: one side of the crack, 1: other side of the crack
-  cellsToModifyConn0_torenum = DataArrayInt.New()
-  cellsToModifyConn0_torenum.alloc(0,1)
+  seed, nIter, cnt = 0, 0, 0
+  nIterMax = nCells+1 # Safety net for the loop
+  hitCells = DataArrayInt.New(); hitCells.alloc(nCells)
+  hitCells.fillWithValue(0)  # 0 : not hit, -1: one side of the crack, +1: other side of the crack
+  MAX_CP = 10000 # the choices below assume we won't have more than 10000 different connex parts ...
+  PING_FULL_init, PING_PART =  0, MAX_CP
+  PONG_FULL_init, PONG_PART = -0,-MAX_CP
   while nIter < nIterMax:
 #       print("dbg ", hitCells.getValues())
-      t = hitCells.findIdsEqual(-1)
+      t = hitCells.findIdsEqual(0)
       if not t.getNumberOfTuples():
         break
       seed = t[0,0]
       done = False
-      PING_FULL, PING_PART = 1,11
-      PONG_FULL, PONG_PART = 2,22
+      cnt += 1
+      PING_FULL = PING_FULL_init+cnt
+      PONG_FULL = PONG_FULL_init-cnt
       while not done and nIter < nIterMax:  # Start of the ping-pong
           nIter += 1
           # Identify connex zone around the seed
-          spreadZone, _ = MEDCouplingUMesh.ComputeSpreadZoneGraduallyFromSeed([seed],  neigh00,neighI00, -1)
+          spreadZone, _ = MEDCouplingUMesh.ComputeSpreadZoneGraduallyFromSeed([seed],  neigh,neighI, -1)
           done = True
           for i, s in enumerate(spreadZone.getValues()):
               hitCells[s] = PING_FULL
@@ -175,20 +193,88 @@ def findCellsToRenumber(this, otherDimM1OnSameCoords, dupl):
               PING_FULL, PONG_FULL = PONG_FULL, PING_FULL
               PING_PART, PONG_PART = PONG_PART, PING_PART
 
-      nonHitCells = hitCells.findIdsEqual(-1)
+      nonHitCells = hitCells.findIdsEqual(0)
       if nonHitCells.getNumberOfTuples():
         seed = nonHitCells[0,0]
       else:
         break
+      
   if nIter >= nIterMax:
     raise ValueError("Too many iterations - should not happen")
 
-  cellsToModifyConn0_torenum = hitCells.findIdsEqual(PONG_FULL)
-  cellsToModifyConn1_torenum = hitCells.findIdsEqual(PING_FULL)
-  if cellsToModifyConn0_torenum.getNumberOfTuples() + cellsToModifyConn1_torenum.getNumberOfTuples() != cellsAroundGroup.getNumberOfTuples():
-    raise ValueError("Some cells not hit - Internal error should not happen")
-  cellsToModifyConn0_torenum.transformWithIndArr(cellsAroundGroup)
-  cellsToModifyConn1_torenum.transformWithIndArr(cellsAroundGroup)
+  # Now we have handled all N D cells which have a face touching the M1 group. It remains the cells
+  # which are just touching the group by one (or several) node(s):
+  # All those cells are in direct contact with a cell which is either PING_FULL or PONG_FULL
+  # So first reproject the PING/PONG info onto mAroundGrpLarge:
+  hitCellsLarge = DataArrayInt.New(); hitCellsLarge.alloc(nCellsLarge)
+  hitCellsLarge.fillWithValue(0)
+  hitCellsLarge[idsStrict] = hitCells
+  nonHitCells = hitCellsLarge.findIdsEqual(0)
+  # Neighbor information in mAroundGrpLarge:
+  neighL, neighIL = MEDCouplingUMesh.ComputeNeighborsOfCellsAdv(descL,descIL,revDescL,revDescIL)
+  for c in nonHitCells:
+      neighs = neighL[neighIL[c[0]]:neighIL[c[0]+1]]
+      ok = False
+      for n in neighs:
+          neighVal = hitCellsLarge[n[0]]
+          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
+              currVal = hitCellsLarge[c[0]]
+              if currVal != 0:   # Several neighbors have a candidate number
+                  # Unfortunately in some weird cases (see testBuildInnerBoundary8) a cell in mAroundGrpLarge
+                  # might have as neighbor two conflicting spread zone ...
+                  if currVal*neighVal < 0:
+                      # If we arrive here, the cell was already assigned a number and we found a neighbor with
+                      # a different sign ... we must swap the whole spread zone!!
+                      print("Ouch - must switch spread zones ...")
+                      ids1 = hitCellsLarge.findIdsEqual(neighVal)
+                      ids1b = hitCellsLarge.findIdsEqual(-neighVal)
+                      ids2 = hitCellsLarge.findIdsEqual(MAX_CP*neighVal)
+                      ids2b = hitCellsLarge.findIdsEqual(-MAX_CP*neighVal)
+                      hitCellsLarge[ids1] *= -1
+                      hitCellsLarge[ids1b] *= -1
+                      hitCellsLarge[ids2] *= -1
+                      hitCellsLarge[ids2b] *= -1
+              else:  # First assignation
+                  hitCellsLarge[c[0],0] = MAX_CP*neighVal   # Same sign, but different value to preserve PING_FULL and PONG_FULL
+              ok = True
+      if not ok: 
+          print("Hitting 'naked' cell")
+#           raise ValueError("Some cells in mAroundGrpLarge not computable!! Should not happen.")
+  
+  # Now handling remaining cells not touched by the above process, called "naked" cells (see cell #20 in mArndLarge in testBuildInnerBoundary8() ...)
+  naked = hitCellsLarge.findIdsEqual(0)
+  mLargC, mLargCI = mArGrpLargeDesc.getNodalConnectivity(),mArGrpLargeDesc.getNodalConnectivityIndex()
+  for c in naked:
+      neighs = neighL[neighIL[c[0]]:neighIL[c[0]+1]]  # ExtractFromIndexedArray?
+      nbDup = {}
+      fac1 = descL[descIL[c[0]]:descIL[c[0]+1]]
+      for n in neighs:
+          if hitCellsLarge[n[0]] == 0:
+              continue   # this neighbour is naked too, nothing we can do for now
+          # Among the values found on neighbour cells, take the one from the neighbour which is connected
+          # with the most "economical" face, i.e. the face made of a minimal number of duplicated points.
+          # TODO: this is a shaky criteria ... find sth more robust ...
+          #   1. find face(s) making the link
+          fac2 = descL[descIL[n[0]]:descIL[n[0]+1]]
+          com = fac1.buildIntersection(fac2)
+          if (com.getNumberOfTuples() == 0):
+              raise ValueError("Internal error : no common face ?")
+          #   2. count number of duplicated node for this face.
+          for f in com: # for all common faces
+              faceNodes = mLargC[mLargCI[f[0]]+1:mLargCI[f[0]+1]] # first +1 to skip type
+              comNod = faceNodes.buildIntersection(dupl)
+              # in case the two cells are in contact by multiple faces, take the most conservative value
+              nbDup[n[0]] = max(nbDup.get(n[0],-1), comNod.getNumberOfTuples())
+      # Minimal value in nbDup?
+      cellIdx = min(nbDup, key=nbDup.get)
+      hitCellsLarge[c[0]] = hitCellsLarge[cellIdx]
+  
+  cellsToModifyConn0_torenum = hitCellsLarge.findIdsInRange(1,MAX_CP*MAX_CP)    # Positive spread zone number
+  cellsToModifyConn1_torenum = hitCellsLarge.findIdsInRange(-MAX_CP*MAX_CP, 0)  # Negative spread zone number
+  if cellsToModifyConn0_torenum.getNumberOfTuples() + cellsToModifyConn1_torenum.getNumberOfTuples() != cellsAroundGroupLarge.getNumberOfTuples():
+      raise ValueError("Some cells not hit - Internal error should not happen")
+  cellsToModifyConn0_torenum.transformWithIndArr(cellsAroundGroupLarge)
+  cellsToModifyConn1_torenum.transformWithIndArr(cellsAroundGroupLarge)
   #
   cellIdsNeededToBeRenum=cellsToModifyConn0_torenum
   cellIdsNotModified=cellsToModifyConn1_torenum
index f4e1789e3e7a42c6e1a77ac6d531d2cfa4620f7e..f7ca33b4c8fc6d1f0a6707361d961ecadf58c395 100755 (executable)
@@ -2494,7 +2494,6 @@ DataArrayIdType* MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh&
   else  // if (3D ...)
     notDup = xtrem->buildSubstraction(fNodes);
 
-  // Now compute cells around group (i.e. cells where we will do the propagation to identify the two sub-sets delimited by the group)
   DAInt m1Nodes = otherDimM1OnSameCoords.computeFetchedNodeIds();
   DAInt dupl = m1Nodes->buildSubstraction(notDup);
   return dupl.retn();
@@ -2529,59 +2528,85 @@ void MEDCouplingUMesh::findCellsToRenumber(const MEDCouplingUMesh& otherDimM1OnS
   if(otherDimM1OnSameCoords.getMeshDimension()!=getMeshDimension()-1)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: the mesh given in other parameter must have this->getMeshDimension()-1 !");
 
-  DAInt cellsAroundGroup = getCellIdsLyingOnNodes(nodeIdsToDuplicateBg, nodeIdsToDuplicateEnd, false);  // false= take cell in, even if not all nodes are in dupl
+  DAInt cellsAroundGroupLarge = getCellIdsLyingOnNodes(nodeIdsToDuplicateBg, nodeIdsToDuplicateEnd, false);  // false= take cell in, even if not all nodes are in dupl
 
   //
-  MCUMesh mAroundGrp=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(cellsAroundGroup->begin(),cellsAroundGroup->end(),true));
-  mcIdType nCells2 = mAroundGrp->getNumberOfCells();  
-  DAInt desc00=DataArrayIdType::New(),descI00=DataArrayIdType::New(),revDesc00=DataArrayIdType::New(),revDescI00=DataArrayIdType::New();
-  MCUMesh mArGrpDesc=mAroundGrp->buildDescendingConnectivity(desc00,descI00,revDesc00,revDescI00);
+  MCUMesh mAroundGrpLarge=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(cellsAroundGroupLarge->begin(),cellsAroundGroupLarge->end(),true));
+  DAInt descL=DataArrayIdType::New(),descIL=DataArrayIdType::New(),revDescL=DataArrayIdType::New(),revDescIL=DataArrayIdType::New();
+  MCUMesh mArGrpLargeDesc=mAroundGrpLarge->buildDescendingConnectivity(descL,descIL,revDescL,revDescIL);
+  const mcIdType *descILP=descIL->begin(), *descLP=descL->begin();
 
-  // Identify cells of M1 group in sub-mesh mAroundGrp
+  // Extract now all N D cells which have a complete face in touch with the group:
+  //   1. Identify cells of M1 group in sub-mesh mAroundGrp
   DataArrayIdType *idsOfM1t;
-  mArGrpDesc->areCellsIncludedIn(&otherDimM1OnSameCoords,2, idsOfM1t);
-  DAInt idsOfM1(idsOfM1t);
-
-  // Build map giving for each cell ID in mAroundGrp the corresponding cell ID on the other side of the crack:
-  // Note that this does not cover all cells around the crack (a cell like a triangle might touche the crack only with its tip)
-  std::map<mcIdType, mcIdType> toOtherSide;
-  const mcIdType *revDP=revDesc00->begin(), *revDIP=revDescI00->begin();
-  const mcIdType nCellsArGrpDesc = mArGrpDesc->getNumberOfCells();
-  for (const auto& v: *idsOfM1)
-    {
-      if (v >= nCellsArGrpDesc)   // Keep valid match only
+  mArGrpLargeDesc->areCellsIncludedIn(&otherDimM1OnSameCoords,2, idsOfM1t);
+  DAInt idsOfM1Large(idsOfM1t);
+  mcIdType nL = mArGrpLargeDesc->getNumberOfCells();
+  DAInt idsStrict = DataArrayIdType::New(); idsStrict->alloc(0,1);
+  //  2. Build map giving for each cell ID in mAroundGrp (not in mAroundGrpLarge) the corresponding cell
+  //     ID on the other side of the crack:
+  std::map<mcIdType, mcIdType> toOtherSide, pos;
+  mcIdType cnt = 0;
+  const mcIdType *revDescILP=revDescIL->begin(), *revDescLP=revDescL->begin();
+  for(const auto& v: *idsOfM1Large)
+    {
+      if (v >= nL)    // Keep valid match only
         continue;
-      mcIdType idx0 = revDIP[v];
-      mcIdType c1=revDP[idx0], c2=revDP[idx0+1];
-      toOtherSide[c1] = c2;
-      toOtherSide[c2] = c1;
-    }
+      mcIdType idx0 = revDescILP[v];
+      // Keep the two cells on either side of the face v of M1:
+      mcIdType c1=revDescLP[idx0], c2=revDescLP[idx0+1];
+      DAInt t1=idsStrict->findIdsEqual(c1), t2=idsStrict->findIdsEqual(c2);
+
+      if (!t1->getNumberOfTuples())
+        {  pos[c1] = cnt++; idsStrict->pushBackSilent(c1);   }
+      if (!t2->getNumberOfTuples())
+        {  pos[c2] = cnt++; idsStrict->pushBackSilent(c2);   }
+
+      mcIdType k1 = pos[c1], k2=pos[c2];
+      toOtherSide[k1] = k2;
+      toOtherSide[k2] = k1;
+    }
+
+  DAInt cellsAroundGroup = cellsAroundGroupLarge->selectByTupleId(idsStrict->begin(), idsStrict->end());
+  MCUMesh mAroundGrp = static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(cellsAroundGroup->begin(), cellsAroundGroup->end(), true));
+  mcIdType nCells=cellsAroundGroup->getNumberOfTuples(), nCellsLarge=cellsAroundGroupLarge->getNumberOfTuples();
+  DAInt desc=DataArrayIdType::New(),descI=DataArrayIdType::New(),revDesc=DataArrayIdType::New(),revDescI=DataArrayIdType::New();  
+  MCUMesh mArGrpDesc=mAroundGrp->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
+  DataArrayIdType *idsOfM1t2;
+  mArGrpDesc->areCellsIncludedIn(&otherDimM1OnSameCoords,2, idsOfM1t2);  // TODO can we avoid recomputation here?
+  DAInt idsOfM1(idsOfM1t2);
 
   // Neighbor information of the mesh WITH the crack (some neighbors are removed):
   //     In the neighbor information remove the connection between high dimension cells and its low level constituents which are part
   //     of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack):
-  DataArrayIdType::RemoveIdsFromIndexedArrays(idsOfM1->begin(), idsOfM1->end(),desc00,descI00);
+  DataArrayIdType::RemoveIdsFromIndexedArrays(idsOfM1->begin(), idsOfM1->end(),desc,descI);
   //     Compute the neighbor of each cell in mAroundGrp, taking into account the broken link above. Two
   //     cells on either side of the crack (defined by the mesh of low dimension) are not neighbor anymore.
-  DataArrayIdType *neigh00t=0, *neighI00t=0;
-  MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00, neigh00t, neighI00t);
-  DAInt neigh00(neigh00t), neighI00(neighI00t);
+  DataArrayIdType *neight=0, *neighIt=0;
+  MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(desc,descI,revDesc,revDescI, neight, neighIt);
+  DAInt neigh(neight), neighI(neighIt);
 
   // For each initial connex part of the M1 mesh (or said differently for each independent crack):
   mcIdType seed=0, nIter=0;
-  mcIdType nIterMax = nCells2+1; // Safety net for the loop
-  DAInt hitCells = DataArrayIdType::New(); hitCells->alloc(nCells2,1);
-  hitCells->fillWithValue(-1);  // -1 : not hit, 0: one side of the crack, 1: other side of the crack
+  mcIdType nIterMax = nCells+1; // Safety net for the loop
+  DAInt hitCells = DataArrayIdType::New(); hitCells->alloc(nCells,1);
+  mcIdType* hitCellsP = hitCells->rwBegin();
+  hitCells->fillWithValue(0);  // 0 : not hit, +x: one side of the crack, -x: other side of the crack, with 'x' the index of the connex component
   mcIdType PING_FULL, PONG_FULL;
+  mcIdType MAX_CP = 10000;  // the choices below assume we won't have more than 10000 different connex parts ...
+  mcIdType PING_FULL_init = 0, PING_PART = MAX_CP;
+  mcIdType PONG_FULL_init = 0, PONG_PART = -MAX_CP;
+  cnt=0;
   while (nIter < nIterMax)
     {
-      DAInt t = hitCells->findIdsEqual(-1);
+      DAInt t = hitCells->findIdsEqual(0);
       if(!t->getNumberOfTuples())
         break;
       mcIdType seed = t->getIJ(0,0);
       bool done = false;
-      PING_FULL         =1;  PONG_FULL=2;
-      mcIdType PING_PART=11, PONG_PART=22;
+      cnt++;
+      PING_FULL = PING_FULL_init+cnt;
+      PONG_FULL = PONG_FULL_init-cnt;
       // while the connex bits in correspondance on either side of the crack are not fully covered
       while(!done && nIter < nIterMax)  // Start of the ping-pong
         {
@@ -2590,9 +2615,8 @@ void MEDCouplingUMesh::findCellsToRenumber(const MEDCouplingUMesh& otherDimM1OnS
           // of the crack that might extend further away. So we will need to compute spread zone on the other side
           // too ... and this process can repeat, hence the "ping-pong" logic.
           mcIdType dnu;
-          DAInt spreadZone = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neigh00,neighI00, -1, dnu);
+          DAInt spreadZone = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neigh,neighI, -1, dnu);
           done = true;
-          mcIdType* hitCellsP = hitCells->rwBegin();
           for(const mcIdType& s: *spreadZone)
             {
               hitCellsP[s] = PING_FULL;
@@ -2629,7 +2653,7 @@ void MEDCouplingUMesh::findCellsToRenumber(const MEDCouplingUMesh& otherDimM1OnS
               std::swap(PING_PART, PONG_PART);
             }
         } // while (!done ...)
-      DAInt nonHitCells = hitCells->findIdsEqual(-1);
+      DAInt nonHitCells = hitCells->findIdsEqual(0);
       if (nonHitCells->getNumberOfTuples())
         seed = nonHitCells->getIJ(0,0);
       else
@@ -2638,12 +2662,110 @@ void MEDCouplingUMesh::findCellsToRenumber(const MEDCouplingUMesh& otherDimM1OnS
   if (nIter >= nIterMax)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: Too many iterations - should not happen");
 
-  DAInt cellsRet1 = hitCells->findIdsEqual(PONG_FULL);
-  DAInt cellsRet2 = hitCells->findIdsEqual(PING_FULL);
-  if (cellsRet1->getNumberOfTuples() + cellsRet2->getNumberOfTuples() != cellsAroundGroup->getNumberOfTuples())
+  // Now we have handled all N D cells which have a face touching the M1 group. It remains the cells
+  // which are just touching the group by one (or several) node(s):
+  // All those cells are in direct contact with a cell which is either PING_FULL or PONG_FULL
+  // So first reproject the PING/PONG info onto mAroundGrpLarge:
+  DAInt hitCellsLarge = DataArrayIdType::New(); hitCellsLarge->alloc(nCellsLarge,1);
+  hitCellsLarge->fillWithValue(0);
+  mcIdType *hitCellsLargeP=hitCellsLarge->rwBegin(), tt=0;
+  for(const auto &i: *idsStrict)
+    { hitCellsLargeP[i] = hitCellsP[tt++]; }
+  DAInt nonHitCells = hitCellsLarge->findIdsEqual(0);
+  // Neighbor information in mAroundGrpLarge:
+  DataArrayIdType *neighLt=0, *neighILt=0;
+  MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(descL,descIL,revDescL,revDescIL, neighLt, neighILt);
+  DAInt neighL(neighLt), neighIL(neighILt);
+  const mcIdType *neighILP=neighIL->begin(), *neighLP=neighL->begin();
+  for(const auto& c : *nonHitCells)
+    {
+      mcIdType cnt00 = neighILP[c];
+      for (const mcIdType *n=neighLP+cnt00; cnt00 < neighILP[c+1]; n++, cnt00++)
+        {
+          mcIdType neighVal = hitCellsLargeP[*n];
+          if (neighVal != 0 && std::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
+            {
+              mcIdType currVal = hitCellsLargeP[c];
+              if (currVal != 0)   // Several neighbors have a candidate number
+                {
+                  // Unfortunately in some weird cases (see testBuildInnerBoundary8) a cell in mAroundGrpLarge
+                  // might have as neighbor two conflicting spread zone ...
+                  if (currVal*neighVal < 0)
+                    {
+                      // If we arrive here, the cell was already assigned a number and we found a neighbor with
+                      // a different sign ... we must swap the whole spread zone!!
+                      DAInt ids1 = hitCellsLarge->findIdsEqual(neighVal), ids1b = hitCellsLarge->findIdsEqual(-neighVal);
+                      DAInt ids2 = hitCellsLarge->findIdsEqual(MAX_CP*neighVal), ids2b = hitCellsLarge->findIdsEqual(-MAX_CP*neighVal);
+                      // A nice little lambda to multiply part of a DAInt by -1 ...
+                      auto mul_part_min1 = [hitCellsLargeP](const DAInt& ids) { for(const auto& i: *ids) hitCellsLargeP[i] *= -1; };
+                      mul_part_min1(ids1);
+                      mul_part_min1(ids1b);
+                      mul_part_min1(ids2);
+                      mul_part_min1(ids2b);
+                    }
+                }
+              else  // First assignation
+                hitCellsLargeP[c] = MAX_CP*neighVal;  // Same sign, but different value to preserve PING_FULL and PONG_FULL
+            }
+        }
+    }
+
+  // Now handling remaining cells not touched by the for loop above, called "naked" cells (see cell #20 in mArndGrpLarge in testBuildInnerBoundary8() ...)
+  DAInt naked = hitCellsLarge->findIdsEqual(0);
+  const mcIdType *mLargCP=mArGrpLargeDesc->getNodalConnectivity()->begin(), *mLargCIP=mArGrpLargeDesc->getNodalConnectivityIndex()->begin();
+  for (const auto &c: *naked)
+    {
+      std::map<mcIdType, mcIdType> nbDup;
+      // Retrieve list of faces of cell c
+      mcIdType nbFac1=descILP[c+1]-descILP[c];
+      std::vector<mcIdType> fac1(nbFac1);
+      std::copy(descLP+descILP[c], descLP+descILP[c+1], fac1.begin());
+      std::sort(fac1.begin(), fac1.end());
+      mcIdType cnt00 = neighILP[c];
+      for (const mcIdType *n=neighLP+cnt00; cnt00 < neighILP[c+1]; n++, cnt00++)
+        {
+          if (hitCellsLargeP[*n] == 0)
+            continue;   // this neighbour is naked too, nothing we can do for now
+          // Among the values found on neighbour cells, take the one from the neighbour which is connected
+          // with the most "economical" face, i.e. the face made of a minimal number of duplicated points.
+          // TODO: this is a shaky criteria ... find sth more robust ...
+          //   1. find face(s) making the link
+          mcIdType nbFac2=descILP[*n+1]-descILP[*n];
+          std::vector<mcIdType> fac2(nbFac2);
+          std::copy(descLP+descILP[*n], descLP+descILP[*n+1], fac2.begin());
+          std::sort(fac2.begin(), fac2.end());
+          std::vector<mcIdType> comFac;
+          std::set_intersection(fac1.begin(), fac1.end(),
+                             fac2.begin() ,fac2.end(),
+                             std::back_inserter(comFac));
+          if (comFac.size() == 0)
+            throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: internal error - no common face between two cells - should not happen");
+          //   2. count number of duplicated node for this face.
+          for (const auto &f : comFac) // for all common faces
+            {
+              std::vector<mcIdType> comNod;
+              std::set_intersection(nodeIdsToDuplicateBg, nodeIdsToDuplicateEnd,
+                                    mLargCP+mLargCIP[f]+1, mLargCP+mLargCIP[f+1],    // first +1 to skip type in connectivity
+                                    std::back_inserter(comNod));
+              // in case the two cells are in contact by multiple faces, take the most conservative value
+              mcIdType val=-1;
+              if(nbDup.find(*n) != nbDup.end()) val=nbDup[*n];
+              nbDup[*n] = std::max(val, (mcIdType)comNod.size());
+            }
+        }
+      // Minimal value in nbDup?
+      using PairId = std::pair<mcIdType, mcIdType>;
+      auto comp_fonc = [](const PairId& p1, const PairId& p2) { return p1.second < p2.second; };
+      PairId zemin = *min_element(nbDup.begin(), nbDup.end(), comp_fonc);
+      hitCellsLargeP[c] = hitCellsLargeP[zemin.first];
+    }
+
+  DAInt cellsRet1 = hitCellsLarge->findIdsInRange(1,MAX_CP*MAX_CP);   // Positive spread zone number
+  DAInt cellsRet2 = hitCellsLarge->findIdsInRange(-MAX_CP*MAX_CP, 0); // Negative spread zone number
+  if (cellsRet1->getNumberOfTuples() + cellsRet2->getNumberOfTuples() != cellsAroundGroupLarge->getNumberOfTuples())
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: Some cells not hit - Internal error should not happen");
-  cellsRet1->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end());
-  cellsRet2->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end());
+  cellsRet1->transformWithIndArr(cellsAroundGroupLarge->begin(),cellsAroundGroupLarge->end());
+  cellsRet2->transformWithIndArr(cellsAroundGroupLarge->begin(),cellsAroundGroupLarge->end());
   //
   cellIdsNeededToBeRenum=cellsRet1.retn();
   cellIdsNotModified=cellsRet2.retn();
index 5960607e471c8c0fcae973e7c03185ebd9250061..cf11334d1427004469c17e8210b7c27eba6dd543 100644 (file)
@@ -1628,7 +1628,7 @@ class MEDLoaderTest3(unittest.TestCase):
     @WriteInTmpDir
     def testBuildInnerBoundary7(self):
         """ 3D test where the crack has another funny shape with another singular point (i.e. two faces of the M1 group are only connected by one point, not a full segment)
-        Once the crack is inserted 
+        Once the crack is inserted, the cells on either side of the crack do not necessarily form a connex spread zone. This was not properly handled either. 
         """
         m3 = MEDCouplingUMesh('box', 3)
         coo = DataArrayDouble([(5,17,0),(0,17,0),(0,12,0),(5,12,0),(15,17,0),(15,12,0),(20,12,0),(20,17,0),(20,2,0),(15,2,0),(15,-3,0),(20,-3,0),(5,-3,0),(5,2,0),(0,-3,0),(0,2,0),(5,17,10),(5,17,20),(5,17,30),(5,17,40),(0,17,10),(0,17,20),(0,17,30),(0,17,40),(0,12,10),(0,12,20),(0,12,30),(0,12,40),(5,12,10),(5,12,20),(5,12,30),(5,12,40),(15,17,10),(15,17,20),(15,17,30),(15,17,40),(15,12,10),(15,12,20),(15,12,30),(15,12,40),(20,12,10),(20,12,20),(20,12,30),(20,12,40),(20,17,10),(20,17,20),(20,17,30),(20,17,40),(20,2,10),(20,2,20),(20,2,30),(20,2,40),(15,2,10),(15,2,20),(15,2,30),(15,2,40),(15,-3,10),(15,-3,20),(15,-3,30),(15,-3,40),(20,-3,10),(20,-3,20),(20,-3,30),(20,-3,40),
@@ -1651,7 +1651,6 @@ class MEDLoaderTest3(unittest.TestCase):
         mfu.setMeshAtLevel(0, m3)
         mfu.setMeshAtLevel(-1, m2)
         mfu.setGroupsAtLevel(-1, [grpIds])
-        mfu.write("/tmp/tst7.med", 2)
         nNod = m3.getNumberOfNodes()
         nodesDup, cells1, cells2 = mfu.buildInnerBoundaryAlongM1Group("group")
         m3_bis = mfu.getMeshAtLevel(0)
@@ -1663,14 +1662,72 @@ class MEDLoaderTest3(unittest.TestCase):
         self.assertEqual(nNod+22, m2_bis.getNumberOfNodes())
         self.assertEqual([0, 3, 12, 13, 16, 17, 18, 19, 28, 29, 30, 31, 64, 65, 66, 67, 68, 69, 70, 71, 88, 89], nodesDup.getValues())
         self.assertEqual(m3_bis.getCoords()[nodesDup].getValues(), m3_bis.getCoords()[nNod:].getValues())
-        self.assertEqual(set([4, 5, 6, 7, 20, 21, 22, 23, 32, 33, 34, 35]), set(cells1.getValues()))
-        self.assertEqual(set([0, 1, 2, 3, 24, 25, 26, 27, 28, 29, 30, 31]), set(cells2.getValues()))
+        self.assertEqual(set([0, 1, 2, 3, 24, 25, 26, 27, 28, 29, 30, 31]), set(cells1.getValues()))
+        self.assertEqual(set([4, 5, 6, 7, 20, 21, 22, 23, 32, 33, 34, 35]), set(cells2.getValues()))
         self.assertEqual([2, 7, 12, 17, 95, 99, 103, 107, 129, 133, 137, 141],mfu.getGroupArr(-1,"group").getValues())
         self.assertEqual([151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162],mfu.getGroupArr(-1,"group_dup").getValues())  # here only one cell has been duplicated
         m_desc, _, _, _, _ = m3_bis.buildDescendingConnectivity()
         m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999)
         pass
 
+    def testBuildInnerBoundary8(self):
+        """ 3D test where the crack leaves 'naked' cells. If we call a 'close-to-crack cell' a cell which shares a face with the M1 group,
+         a 'naked cell' is a cell that has some node duplicated, but which do not share any face with a 'close-to-crack cell'. In this case
+         it is tricky to decide whether this cell should be renumbered or not ... 
+        """
+        m3 = MEDCouplingUMesh('box', 3)
+        coo = DataArrayDouble([(0,15,0),(0,5,0),(3,5,0),(5,5,0),(5,15,0),(5,20,0),(0,20,0),(15,20,0),(15,15,0),(20,15,0),(20,20,0),(20,5,0),(15,5,0),(15,0,0),(20,0,0),(5,-1.60551e-25,0),(5,3,0),(3,0,0),
+        (3,3,0),(0,0,0),(0,3,0),(0,15,10),(0,15,20),(0,15,30),(0,15,40),(0,5,10),(0,5,20),(0,5,30),(0,5,40),(3,5,10),(3,5,20),(3,5,30),(3,5,40),(5,5,10),(5,5,20),(5,5,30),(5,5,40),(5,15,10),(5,15,20),(5,15,30),
+        (5,15,40),(5,20,10),(5,20,20),(5,20,30),(5,20,40),(0,20,10),(0,20,20),(0,20,30),(0,20,40),(15,20,10),(15,20,20),(15,20,30),(15,20,40),(15,15,10),(15,15,20),(15,15,30),(15,15,40),(20,15,10),(20,15,20),
+        (20,15,30),(20,15,40),(20,20,10),(20,20,20),(20,20,30),(20,20,40),(20,5,10),(20,5,20),(20,5,30),(20,5,40),(15,5,10),(15,5,20),(15,5,30),(15,5,40),(15,0,10),(15,0,20),(15,0,30),(15,0,40),(20,0,10),
+        (20,0,20),(20,0,30),(20,0,40),(5,-1.60551e-25,10),(5,-1.60551e-25,20),(5,-1.60551e-25,30),(5,-1.60551e-25,40),(5,3,10),(5,3,20),(5,3,30),(5,3,40),(3,0,10),(3,0,20),(3,0,30),(3,0,40),(3,3,10),(3,3,20),
+        (3,3,30),(3,3,40),(0,0,10),(0,0,20),(0,0,30),(0,0,40),(0,3,10),(0,3,20),(0,3,30),(0,3,40),(0,9,0),(3,9,0),(20,9,0),(0,9,10),(0,9,20),(0,9,30),(0,9,40),(3,9,10),(3,9,20),(3,9,30),(3,9,40),(5,9,30),
+        (5,9,40),(20,9,10),(20,9,20),(20,9,30),(20,9,40),(15,9,30),(15,9,40),(0,15,0),(20,15,0),(0,15,10),(0,15,20),(0,15,30),(0,15,40),(5,15,30),(5,15,40),(15,15,30),(15,15,40),(20,15,10),(20,15,20),(20,15,30),
+        (20,15,40)])
+        m3.setCoords(coo)
+        c = DataArrayInt([31, 5, 4, 124, 6, -1, 41, 45, 126, 37, -1, 5, 41, 37, 4, -1, 4, 37, 126, 124, -1, 124, 126, 45, 6, -1, 6, 45, 41, 5, 31, 41, 37, 126, 45, -1, 42, 46, 127, 38, -1, 41, 42, 38, 37, -1, 37, 38, 127, 126, -1, 126, 127, 46, 45, -1, 45, 46, 42, 41, 31, 42, 38, 127, 46, -1, 43, 47, 128, 130, -1, 42, 43, 130, 38, -1, 38, 130, 128, 127, -1, 127, 128, 47, 46, -1, 46, 47, 43, 42, 31, 43, 130, 128, 47,
+        -1, 44, 48, 129, 131, -1, 43, 44, 131, 130, -1, 130, 131, 129, 128, -1, 128, 129, 48, 47, -1, 47, 48, 44, 43, 31, 7, 8, 4, 5, -1, 49, 41, 37, 53, -1, 7, 49, 53, 8, -1, 8, 53, 37, 4, -1, 4, 37, 41, 5, -1, 5, 41, 49, 7, 31, 49, 53, 37, 41, -1, 50, 42, 38, 54, -1, 49, 50, 54, 53, -1, 53, 54, 38, 37, -1, 37, 38, 42, 41, -1, 41, 42, 50, 49, 31, 50, 54, 38, 42, -1, 51, 43, 130, 132, -1, 50, 51, 132, 54, -1, 54, 132,
+        130, 38, -1, 38, 130, 43, 42, -1, 42, 43, 51, 50, 31, 51, 132, 130, 43, -1, 52, 44, 131, 133, -1, 51, 52, 133, 132, -1, 132, 133, 131, 130, -1, 130, 131, 44, 43, -1, 43, 44, 52, 51, 31, 125, 8, 7, 10, -1, 134, 61, 49, 53, -1, 125, 134, 53, 8, -1, 8, 53, 49, 7, -1, 7, 49, 61, 10, -1, 10, 61, 134, 125, 31, 134, 53, 49, 61, -1, 135, 62, 50, 54, -1, 134, 135, 54, 53, -1, 53, 54, 50, 49, -1, 49, 50, 62, 61, -1,
+        61, 62, 135, 134, 31, 135, 54, 50, 62, -1, 136, 63, 51, 132, -1, 135, 136, 132, 54, -1, 54, 132, 51, 50, -1, 50, 51, 63, 62, -1, 62, 63, 136, 135, 31, 136, 132, 51, 63, -1, 137, 64, 52, 133, -1, 136, 137, 133, 132, -1, 132, 133, 52, 51, -1, 51, 52, 64, 63, -1, 63, 64, 137, 136, 31, 107, 12, 8, 9, -1, 118, 57, 53, 69, -1, 107, 118, 69, 12, -1, 12, 69, 53, 8, -1, 8, 53, 57, 9, -1, 9, 57, 118, 107, 31, 118, 69,
+        53, 57, -1, 119, 58, 54, 70, -1, 118, 119, 70, 69, -1, 69, 70, 54, 53, -1, 53, 54, 58, 57, -1, 57, 58, 119, 118, 31, 119, 70, 54, 58, -1, 120, 59, 55, 122, -1, 119, 120, 122, 70, -1, 70, 122, 55, 54, -1, 54, 55, 59, 58, -1, 58, 59, 120, 119, 31, 120, 122, 55, 59, -1, 121, 60, 56, 123, -1, 120, 121, 123, 122, -1, 122, 123, 56, 55, -1, 55, 56, 60, 59, -1, 59, 60, 121, 120, 31, 13, 12, 11, 14, -1, 73, 77, 65, 69,
+        -1, 13, 73, 69, 12, -1, 12, 69, 65, 11, -1, 11, 65, 77, 14, -1, 14, 77, 73, 13, 31, 73, 69, 65, 77, -1, 74, 78, 66, 70, -1, 73, 74, 70, 69, -1, 69, 70, 66, 65, -1, 65, 66, 78, 77, -1, 77, 78, 74, 73, 31, 74, 70, 66, 78, -1, 75, 79, 67, 71, -1, 74, 75, 71, 70, -1, 70, 71, 67, 66, -1, 66, 67, 79, 78, -1, 78, 79, 75, 74, 31, 75, 71, 67, 79, -1, 76, 80, 68, 72, -1, 75, 76, 72, 71, -1, 71, 72, 68, 67, -1, 67, 68, 80,
+        79, -1, 79, 80, 76, 75, 31, 17, 18, 16, 15, -1, 89, 81, 85, 93, -1, 17, 89, 93, 18, -1, 18, 93, 85, 16, -1, 16, 85, 81, 15, -1, 15, 81, 89, 17, 31, 89, 93, 85, 81, -1, 90, 82, 86, 94, -1, 89, 90, 94, 93, -1, 93, 94, 86, 85, -1, 85, 86, 82, 81, -1, 81, 82, 90, 89, 31, 90, 94, 86, 82, -1, 91, 83, 87, 95, -1, 90, 91, 95, 94, -1, 94, 95, 87, 86, -1, 86, 87, 83, 82, -1, 82, 83, 91, 90, 31, 91, 95, 87, 83, -1, 92, 84,
+        88, 96, -1, 91, 92, 96, 95, -1, 95, 96, 88, 87, -1, 87, 88, 84, 83, -1, 83, 84, 92, 91, 31, 19, 20, 18, 17, -1, 97, 89, 93, 101, -1, 19, 97, 101, 20, -1, 20, 101, 93, 18, -1, 18, 93, 89, 17, -1, 17, 89, 97, 19, 31, 97, 101, 93, 89, -1, 98, 90, 94, 102, -1, 97, 98, 102, 101, -1, 101, 102, 94, 93, -1, 93, 94, 90, 89, -1, 89, 90, 98, 97, 31, 98, 102, 94, 90, -1, 99, 91, 95, 103, -1, 98, 99, 103, 102, -1, 102, 103,
+        95, 94, -1, 94, 95, 91, 90, -1, 90, 91, 99, 98, 31, 99, 103, 95, 91, -1, 100, 92, 96, 104, -1, 99, 100, 104, 103, -1, 103, 104, 96, 95, -1, 95, 96, 92, 91, -1, 91, 92, 100, 99, 31, 1, 2, 18, 20, -1, 25, 101, 93, 29, -1, 1, 25, 29, 2, -1, 2, 29, 93, 18, -1, 18, 93, 101, 20, -1, 20, 101, 25, 1, 31, 25, 29, 93, 101, -1, 26, 102, 94, 30, -1, 25, 26, 30, 29, -1, 29, 30, 94, 93, -1, 93, 94, 102, 101, -1, 101, 102,
+        26, 25, 31, 26, 30, 94, 102, -1, 27, 103, 95, 31, -1, 26, 27, 31, 30, -1, 30, 31, 95, 94, -1, 94, 95, 103, 102, -1, 102, 103, 27, 26, 31, 27, 31, 95, 103, -1, 28, 104, 96, 32, -1, 27, 28, 32, 31, -1, 31, 32, 96, 95, -1, 95, 96, 104, 103, -1, 103, 104, 28, 27, 31, 3, 4, 8, 12, -1, 33, 69, 53, 37, -1, 3, 33, 37, 4, -1, 4, 37, 53, 8, -1, 8, 53, 69, 12, -1, 12, 69, 33, 3, 31, 33, 37, 53, 69, -1, 34, 70, 54, 38, -1,
+        33, 34, 38, 37, -1, 37, 38, 54, 53, -1, 53, 54, 70, 69, -1, 69, 70, 34, 33, 31, 34, 38, 54, 70, -1, 116, 122, 55, 39, -1, 34, 116, 39, 38, -1, 38, 39, 55, 54, -1, 54, 55, 122, 70, -1, 70, 122, 116, 34, 31, 116, 39, 55, 122, -1, 117, 123, 56, 40, -1, 116, 117, 40, 39, -1, 39, 40, 56, 55, -1, 55, 56, 123, 122, -1, 122, 123, 117, 116, 31, 16, 18, 2, 3, -1, 85, 33, 29, 93, -1, 16, 85, 93, 18, -1, 18, 93, 29, 2,
+        -1, 2, 29, 33, 3, -1, 3, 33, 85, 16, 31, 85, 93, 29, 33, -1, 86, 34, 30, 94, -1, 85, 86, 94, 93, -1, 93, 94, 30, 29, -1, 29, 30, 34, 33, -1, 33, 34, 86, 85, 31, 86, 94, 30, 34, -1, 87, 35, 31, 95, -1, 86, 87, 95, 94, -1, 94, 95, 31, 30, -1, 30, 31, 35, 34, -1, 34, 35, 87, 86, 31, 87, 95, 31, 35, -1, 88, 36, 32, 96, -1, 87, 88, 96, 95, -1, 95, 96, 32, 31, -1, 31, 32, 36, 35, -1, 35, 36, 88, 87, 31, 4, 3, 106,
+        105, 0, -1, 37, 21, 108, 112, 33, -1, 3, 4, 37, 33, -1, 106, 3, 33, 112, -1, 105, 106, 112, 108, -1, 0, 105, 108, 21, -1, 4, 0, 21, 37, 31, 37, 33, 112, 108, 21, -1, 38, 22, 109, 113, 34, -1, 33, 37, 38, 34, -1, 112, 33, 34, 113, -1, 108, 112, 113, 109, -1, 21, 108, 109, 22, -1, 37, 21, 22, 38, 31, 38, 34, 113, 109, 22, -1, 39, 23, 110, 114, 116, -1, 34, 38, 39, 116, -1, 113, 34, 116, 114, -1, 109, 113, 114, 110,
+        -1, 22, 109, 110, 23, -1, 38, 22, 23, 39, 31, 39, 116, 114, 110, 23, -1, 40, 24, 111, 115, 117, -1, 116, 39, 40, 117, -1, 114, 116, 117, 115, -1, 110, 114, 115, 111, -1, 23, 110, 111, 24, -1, 39, 23, 24, 40, 31, 16, 3, 12, 13, 15, -1, 85, 81, 73, 69, 33, -1, 3, 16, 85, 33, -1, 12, 3, 33, 69, -1, 13, 12, 69, 73, -1, 15, 13, 73, 81, -1, 16, 15, 81, 85, 31, 85, 33, 69, 73, 81, -1, 86, 82, 74, 70, 34, -1, 33, 85,
+        86, 34, -1, 69, 33, 34, 70, -1, 73, 69, 70, 74, -1, 81, 73, 74, 82, -1, 85, 81, 82, 86, 31, 86, 34, 70, 74, 82, -1, 87, 83, 75, 71, 35, -1, 34, 86, 87, 35, -1, 70, 34, 35, 71, -1, 74, 70, 71, 75, -1, 82, 74, 75, 83, -1, 86, 82, 83, 87, 31, 87, 35, 71, 75, 83, -1, 88, 84, 76, 72, 36, -1, 35, 87, 88, 36, -1, 71, 35, 36, 72, -1, 75, 71, 72, 76, -1, 83, 75, 76, 84, -1, 87, 83, 84, 88])
+        cI = DataArrayInt([0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360, 390, 420, 450, 480, 510, 540, 570, 600, 630, 660, 690, 720, 750, 780, 810, 840, 870, 900, 930, 960, 990, 1020, 1050, 1080, 1110, 1140, 1170, 1200, 1237, 1274, 1311, 1348, 1385, 1422, 1459, 1496])
+        m3.setConnectivity(c, cI)
+        m3.checkConsistency()
+        m2, _, _, _, _ = m3.buildDescendingConnectivity()
+        grpIds = DataArrayInt([2,7,12,17,101,106,111,116,160,164,170,173,176,179]); grpIds.setName("group")
+        mfu = MEDFileUMesh()
+        mfu.setMeshAtLevel(0, m3)
+        mfu.setMeshAtLevel(-1, m2)
+        mfu.setGroupsAtLevel(-1, [grpIds])
+        nNod = m3.getNumberOfNodes()
+        nodesDup, cells1, cells2 = mfu.buildInnerBoundaryAlongM1Group("group")
+        m3_bis = mfu.getMeshAtLevel(0)
+        m3_bis.checkConsistency()
+        m2_bis = mfu.getMeshAtLevel(-1)
+        m2_bis.checkConsistency()
+        self.assertEqual(nNod+27, mfu.getNumberOfNodes())
+        self.assertEqual(nNod+27, m3_bis.getNumberOfNodes())
+        self.assertEqual(nNod+27, m2_bis.getNumberOfNodes())
+        self.assertEqual([3, 4, 5, 15, 16, 33, 35, 36, 37, 39, 40, 41, 42, 43, 44, 81, 82, 83, 84, 85, 86, 87, 88, 116, 117, 130, 131], nodesDup.getValues())
+        self.assertEqual(m3_bis.getCoords()[nodesDup].getValues(), m3_bis.getCoords()[nNod:].getValues())
+        self.assertEqual(set([4, 5, 6, 7, 32, 33, 34, 35, 44, 45, 46, 47]), set(cells1.getValues()))
+        self.assertEqual(set([0, 1, 2, 3, 20, 21, 22, 23, 36, 37, 38, 39, 40, 41, 42, 43]), set(cells2.getValues()))
+        self.assertEqual([2, 7, 12, 17, 101, 106, 111, 116, 160, 164, 170, 173, 176, 179],mfu.getGroupArr(-1,"group").getValues())
+        self.assertEqual([212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225],mfu.getGroupArr(-1,"group_dup").getValues())  # here only one cell has been duplicated
+        m_desc, _, _, _, _ = m3_bis.buildDescendingConnectivity()
+        m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999)
+        pass
+
     @WriteInTmpDir
     def testBasicConstructors(self):
         GeneratePyfile18(self)