From b69e93f036d47b263d1840f5edae302f03ea5504 Mon Sep 17 00:00:00 2001 From: abn Date: Tue, 5 Oct 2021 15:52:43 +0200 Subject: [PATCH] Bug fix: buildInnerBoundaryAlongM1Group + separated findNodesToDuplicate() into two parts + some singular points were improperly identified + some situation where the cells around the M1 group form non-connex patterns were not handled properly. --- resources/dev/my_findNodesToDup.py | 334 +++++++++++++++++++++ src/MEDCoupling/MEDCouplingUMesh.cxx | 350 +++++++++++++++++------ src/MEDCoupling/MEDCouplingUMesh.hxx | 3 +- src/MEDCoupling_Swig/MEDCouplingCommon.i | 16 +- src/MEDLoader/MEDFileMesh.cxx | 15 +- src/MEDLoader/Swig/MEDLoaderTest3.py | 105 ++++++- 6 files changed, 730 insertions(+), 93 deletions(-) create mode 100644 resources/dev/my_findNodesToDup.py diff --git a/resources/dev/my_findNodesToDup.py b/resources/dev/my_findNodesToDup.py new file mode 100644 index 000000000..011e0928b --- /dev/null +++ b/resources/dev/my_findNodesToDup.py @@ -0,0 +1,334 @@ +""" Python version of MEDCouplingUMesh::findNodesToDuplicate() and MEDCouplingUMesh::findCellsToRenumber() methods which are at the core of the + MEDFileUMesh::buildInnerBoundaryAlongM1Group() algorithm. + This greatly helps algorithm tuning ... +""" + +from medcoupling import * + +def findNodesToDuplicate(this, otherDimM1OnSameCoords): + # Checking star-shaped M1 group: + meshM2, _,_,_,rdit0 = otherDimM1OnSameCoords.buildDescendingConnectivity() # 2D: a mesh of points, 3D: a mesh of segs + dsi = rdit0.deltaShiftIndex() + 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. + if(idsTmp0.getNumberOfTuples()): + raise ValueError("") + + # Get extreme nodes from the group (they won't be duplicated except if they also lie on bound of M0 -- see below), + # ie nodes belonging to the boundary "cells" (might be points) of M1 + xtremIdsM2 = dsi.findIdsEqual(1) + meshM2Part = meshM2[xtremIdsM2] + xtrem = meshM2Part.computeFetchedNodeIds() + # Remove from the list points on the boundary of the M0 mesh (those need duplication!) + m0desc, dt0, dit0, rdt0, rdit0 = this.buildDescendingConnectivity() + dsi = rdit0.deltaShiftIndex() + boundSegs = dsi.findIdsEqual(1) # boundary segs/faces of the M0 mesh + m0descSkin = m0desc[boundSegs] + fNodes = m0descSkin.computeFetchedNodeIds() # fNodes needs dupl + # 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) + # 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 + # although they are technically on the skin of the cube. + if this.getMeshDimension() == 3 : + m0descSkinDesc, _, _, _, _ = m0descSkin.buildDescendingConnectivity() # all segments of the skin of the 3D (M0) mesh + _, corresp = meshM2.areCellsIncludedIn(m0descSkinDesc,2) + # validIds is the list of segments which are on both the skin of *this*, and in the segments of the M1 group + # In the cube example above, this is a U shape polyline. + validIds = corresp.findIdsInRange(0, meshM2.getNumberOfCells()) + if validIds.getNumberOfTuples(): + # 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: + # (the U-shaped polyline described above) + m1IntersecSkin = m0descSkinDesc[validIds] + # Its boundary nodes should no be duplicated (this is for example the tip of the crack inside the cube described above) + notDuplSkin = m1IntersecSkin.findBoundaryNodes() + fNodes1 = fNodes.buildSubstraction(notDuplSkin) # fNodes1 needs dupl + + # Specific logic to handle singular points : + # - a point on this U-shape line used in a cell which has no face in common with M1 is deemed singular. + # - indeed, if duplicated, such a point would lead to the duplication of a cell which has no face touching M1 ! The + # algorithm would be duplicating too much ... + # 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: + meshM2Desc, _, _, _, rdit0 = meshM2.buildDescendingConnectivity() # a mesh made of node cells + dsi = rdit0.deltaShiftIndex() + singPoints = dsi.findIdsNotInRange(-1,4) # points connected to (strictly) more than 3 segments + if singPoints.getNumberOfTuples(): + print ("Hitting singular point logic") + boundNodes = m1IntersecSkin.computeFetchedNodeIds() + # If a point on this U-shape line is connected to cells which do not share any face with M1, then it + # should not be duplicated + # 1. Extract N D cells touching U-shape line: + cellsAroundBN = this.getCellIdsLyingOnNodes(boundNodes, False) # false= take cell in, even if not all nodes are in dupl + mAroundBN = this[cellsAroundBN] + mAroundBNDesc, descBN,descIBN,revDescBN,revDescIBN=mAroundBN.buildDescendingConnectivity() + # 2. Identify cells in sub-mesh mAroundBN which have a face in common with M1 + _, idsOfM1BN = mAroundBNDesc.areCellsIncludedIn(otherDimM1OnSameCoords,2) + nCells, nCellsDesc = mAroundBN.getNumberOfCells(), mAroundBNDesc.getNumberOfCells() + idsTouch = DataArrayInt.New(); idsTouch.alloc(0,1) + for v in idsOfM1BN: + if v[0] >= nCellsDesc: # Keep valid match only + continue + idx0 = revDescIBN[v[0], 0] + c1, c2 = revDescBN[idx0, 0], revDescBN[idx0+1,0] + idsTouch.pushBackSilent(c1) + idsTouch.pushBackSilent(c2) + # 3. Build complement + idsTouchCompl = idsTouch.buildComplement(nCells) + mAroundBNStrict = mAroundBN[idsTouchCompl] + nod3 = mAroundBNStrict.computeFetchedNodeIds() + inters = boundNodes.buildIntersection(nod3) + print("sing,", inters.getValues()) + fNodes1 = fNodes1.buildSubstraction(inters) # reminder: fNodes1 represent nodes that need dupl. + notDup = xtrem.buildSubstraction(fNodes1) + else: # if validIds ... + notDup = xtrem.buildSubstraction(fNodes) + else: # if 3D ... + notDup = xtrem.buildSubstraction(fNodes) + + m1Nodes = otherDimM1OnSameCoords.computeFetchedNodeIds() + dupl = m1Nodes.buildSubstraction(notDup) + return dupl + + +def findCellsToRenumber(this, otherDimM1OnSameCoords, dupl): + """ Find cells to renumber + """ + # 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 + # + 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") + + # 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,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. + 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, 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(0) + if not t.getNumberOfTuples(): + break + seed = t[0,0] + done = False + 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], neigh,neighI, -1) + done = True + for i, s in enumerate(spreadZone.getValues()): + hitCells[s] = PING_FULL + if s in toOtherSide: + other = toOtherSide[s] + if hitCells[other] != PONG_FULL: + done = False + hitCells[other] = PONG_PART + # Compute next seed, i.e. a cell on the other side of the crack + seed = other + if done: + # we might have several disjoing PONG parts in front of a single PING connex part: + idsPong = hitCells.findIdsEqual(PONG_PART) + if idsPong.getNumberOfTuples(): + seed = idsPong[0,0] + done = False + continue # continue without switching side (or break if done remains false) + else: + # Go the other side + PING_FULL, PONG_FULL = PONG_FULL, PING_FULL + PING_PART, PONG_PART = PONG_PART, PING_PART + + nonHitCells = hitCells.findIdsEqual(0) + if nonHitCells.getNumberOfTuples(): + seed = nonHitCells[0,0] + else: + break + + if nIter >= nIterMax: + raise ValueError("Too many iterations - should not happen") + + # 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: + assert(False) + neighs = neighL[neighIL[c[0]]:neighIL[c[0]+1]] + 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 + +### +### SO FAR THE LOGIC BELOW WAS NOT NEEDED .... +### + +# # 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) # Positive spread zone number +# cellsToModifyConn1_torenum = hitCellsLarge.findIdsInRange(-MAX_CP, 0) # Negative spread zone number + + +### +### C++ VERSION OF IT +### +# // +# // // 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 nbDup; +# // // Retrieve list of faces of cell c +# // mcIdType nbFac1=descILP[c+1]-descILP[c]; +# // std::vector 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 fac2(nbFac2); +# // std::copy(descLP+descILP[*n], descLP+descILP[*n+1], fac2.begin()); +# // std::sort(fac2.begin(), fac2.end()); +# // std::vector 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 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; +# // 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]; +# // } + + + 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 + + return cellIdsNeededToBeRenum, cellIdsNotModified + diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 86e466238..0762ba266 100755 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -2360,16 +2360,15 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildUnstructured() const * * \param [in] otherDimM1OnSameCoords a mesh lying on the same coords than \b this and with a mesh dimension equal to those of \b this minus 1. WARNING this input * parameter is altered during the call. - * \param [out] nodeIdsToDuplicate node ids needed to be duplicated following the algorithm explain above. - * \param [out] cellIdsNeededToBeRenum cell ids in \b this in which the renumber of nodes should be performed. - * \param [out] cellIdsNotModified cell ids mcIdType \b this that lies on \b otherDimM1OnSameCoords mesh whose connectivity do \b not need to be modified as it is the case for \b cellIdsNeededToBeRenum. + * \return node ids which need to be duplicated following the algorithm explained above. * */ -void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayIdType *& nodeIdsToDuplicate, - DataArrayIdType *& cellIdsNeededToBeRenum, DataArrayIdType *& cellIdsNotModified) const +DataArrayIdType* MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords) const { - typedef MCAuto DAInt; - typedef MCAuto MCUMesh; + // DEBUG NOTE: in case of issue with the algorithm in this method, see Python script in resources/dev + // which mimicks the C++ + using DAInt = MCAuto; + using MCUMesh = MCAuto; checkFullyDefined(); otherDimM1OnSameCoords.checkFullyDefined(); @@ -2395,7 +2394,7 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On // Remove from the list points on the boundary of the M0 mesh (those need duplication!) dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New(); MCUMesh m0desc = buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); dt0=0; dit0=0; rdt0=0; - dsi = rdit0->deltaShiftIndex(); + dsi = rdit0->deltaShiftIndex(); rdit0=0; DAInt boundSegs = dsi->findIdsEqual(1); dsi = 0; // boundary segs/faces of the M0 mesh MCUMesh m0descSkin = static_cast(m0desc->buildPartOfMySelf(boundSegs->begin(),boundSegs->end(), true)); DAInt fNodes = m0descSkin->computeFetchedNodeIds(); @@ -2410,108 +2409,297 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On dnu1=0;dnu2=0;dnu3=0;dnu4=0; DataArrayIdType * corresp=0; meshM2->areCellsIncludedIn(m0descSkinDesc,2,corresp); + // validIds is the list of segments which are on both the skin of *this*, and in the segments of the M1 group + // In the cube example above, this is a U shape polyline. DAInt validIds = corresp->findIdsInRange(0, meshM2->getNumberOfCells()); corresp->decrRef(); if (validIds->getNumberOfTuples()) { // 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: + // (the U-shaped polyline described above) MCUMesh m1IntersecSkin = static_cast(m0descSkinDesc->buildPartOfMySelf(validIds->begin(), validIds->end(), true)); // Its boundary nodes should no be duplicated (this is for example the tip of the crack inside the cube described above) DAInt notDuplSkin = m1IntersecSkin->findBoundaryNodes(); DAInt fNodes1 = fNodes->buildSubstraction(notDuplSkin); - // Also, in this (segment) mesh, nodes connected to more than 3 segs should not be dup either (singular points - see testBuildInnerBoundary6()) - dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New(); - MCUMesh meshM2Desc = meshM2->buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); dt0=0; dit0=0; rdt0=0; // a mesh made of node cells - dsi = rdit0->deltaShiftIndex(); - DAInt singPoints = dsi->findIdsNotInRange(-1,4); // points connected to (strictly) more than 3 segments - const mcIdType *cc = meshM2Desc->getNodalConnectivity()->begin(), *ccI = meshM2Desc->getNodalConnectivityIndex()->begin(); - mcIdType * singPointsP = singPoints->rwBegin(); - for (mcIdType j=0; j < singPoints->getNumberOfTuples(); j++) // replace ids in singPoints by real coordinate index (was index of cells in notDuplSkin) + // Specific logic to handle singular points : + // - a point on this U-shape line used in a cell which has no face in common with M1 is deemed singular. + // - indeed, if duplicated, such a point would lead to the duplication of a cell which has no face touching M1 ! The + // algorithm would be duplicating too much ... + // 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: + dnu1=DataArrayIdType::New(), dnu2=DataArrayIdType::New(), dnu3=DataArrayIdType::New(), rdit0=DataArrayIdType::New(); + MCUMesh meshM2Desc = meshM2->buildDescendingConnectivity(dnu1, dnu2, dnu3, rdit0); // a mesh made of node cells + dnu1=0;dnu2=0;dnu3=0; + dsi = rdit0->deltaShiftIndex(); rdit0=0; + DAInt singPoints = dsi->findIdsNotInRange(-1,4) ; dsi=0;// points connected to (strictly) more than 3 segments + if (singPoints->getNumberOfTuples()) { - mcIdType nodeCellIdx = singPointsP[j]; - singPointsP[j] = cc[ccI[nodeCellIdx]+1]; // +1 to skip type + DAInt boundNodes = m1IntersecSkin->computeFetchedNodeIds(); + // If a point on this U-shape line is connected to cells which do not share any face with M1, then it + // should not be duplicated + // 1. Extract N D cells touching U-shape line: + DAInt cellsAroundBN = getCellIdsLyingOnNodes(boundNodes->begin(), boundNodes->end(), false); // false= take cell in, even if not all nodes are in dupl + MCUMesh mAroundBN = static_cast(this->buildPartOfMySelf(cellsAroundBN->begin(), cellsAroundBN->end(), true)); + DAInt descBN=DataArrayIdType::New(), descIBN=DataArrayIdType::New(), revDescBN=DataArrayIdType::New(), revDescIBN=DataArrayIdType::New(); + MCUMesh mAroundBNDesc = mAroundBN->buildDescendingConnectivity(descBN,descIBN,revDescBN,revDescIBN); + // 2. Identify cells in sub-mesh mAroundBN which have a face in common with M1 + DataArrayIdType *idsOfM1BNt; + mAroundBNDesc->areCellsIncludedIn(&otherDimM1OnSameCoords,2, idsOfM1BNt); + DAInt idsOfM1BN(idsOfM1BNt); + mcIdType nCells=mAroundBN->getNumberOfCells(), nCellsDesc=mAroundBNDesc->getNumberOfCells(); + DAInt idsTouch=DataArrayIdType::New(); idsTouch->alloc(0,1); + const mcIdType *revDescIBNP=revDescIBN->begin(), *revDescBNP=revDescBN->begin(); + for(const auto& v: *idsOfM1BN) + { + if (v >= nCellsDesc) // Keep valid match only + continue; + mcIdType idx0 = revDescIBNP[v]; + // Keep the two cells on either side of the face v of M1: + mcIdType c1=revDescBNP[idx0], c2=revDescBNP[idx0+1]; + idsTouch->pushBackSilent(c1); idsTouch->pushBackSilent(c2); + } + // 3. Build complement + DAInt idsTouchCompl = idsTouch->buildComplement(nCells); + MCUMesh mAroundBNStrict = static_cast(mAroundBN->buildPartOfMySelf(idsTouchCompl->begin(), idsTouchCompl->end(), true)); + DAInt nod3 = mAroundBNStrict->computeFetchedNodeIds(); + DAInt inters = boundNodes->buildIntersection(nod3); + fNodes1 = fNodes1->buildSubstraction(inters); // reminder: fNodes1 represent nodes that need dupl. } - DAInt fNodes2 = fNodes1->buildSubstraction(singPoints); - notDup = xtrem->buildSubstraction(fNodes2); + notDup = xtrem->buildSubstraction(fNodes1); } - else + else // if (validIds-> ...) notDup = xtrem->buildSubstraction(fNodes); } - else + 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); - DAInt cellsAroundGroup = getCellIdsLyingOnNodes(dupl->begin(), dupl->end(), false); // false= take cell in, even if not all nodes are in notDup + return dupl.retn(); +} + + +/*! + * This method expects that \b this and \b otherDimM1OnSameCoords share the same coordinates array. + * otherDimM1OnSameCoords->getMeshDimension() is expected to be equal to this->getMeshDimension()-1. + * This method is part of the MEDFileUMesh::buildInnerBoundaryAlongM1Group() algorithm. + * Given a set of nodes to duplicate, this method identifies which cells should have their connectivity modified + * to produce the inner boundary. It is typically called after findNodesToDuplicate(). + * + * \param [in] otherDimM1OnSameCoords a mesh lying on the same coords than \b this and with a mesh dimension equal to those of \b this minus 1. WARNING this input + * parameter is altered during the call. + * \param [in] nodeIdsToDuplicateBg node ids needed to be duplicated, as returned by findNodesToDuplicate. + * \param [in] nodeIdsToDuplicateEnd node ids needed to be duplicated, as returned by findNodesToDuplicate. + * \param [out] cellIdsNeededToBeRenum cell ids in \b this in which the renumber of nodes should be performed. + * \param [out] cellIdsNotModified cell ids in \b this that lies on \b otherDimM1OnSameCoords mesh whose connectivity do \b not need to be modified as it is the case for \b cellIdsNeededToBeRenum. + * + */ +void MEDCouplingUMesh::findCellsToRenumber(const MEDCouplingUMesh& otherDimM1OnSameCoords, const mcIdType *nodeIdsToDuplicateBg, const mcIdType *nodeIdsToDuplicateEnd, + DataArrayIdType *& cellIdsNeededToBeRenum, DataArrayIdType *& cellIdsNotModified) const +{ + // DEBUG NOTE: in case of issue with the algorithm in this method, see Python script in resources/dev + // which mimicks the C++ + using DAInt = MCAuto; + using MCUMesh = MCAuto; + + checkFullyDefined(); + otherDimM1OnSameCoords.checkFullyDefined(); + if(getCoords()!=otherDimM1OnSameCoords.getCoords()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: meshes do not share the same coords array !"); + if(otherDimM1OnSameCoords.getMeshDimension()!=getMeshDimension()-1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: the mesh given in other parameter must have this->getMeshDimension()-1 !"); + + DAInt cellsAroundGroupLarge = getCellIdsLyingOnNodes(nodeIdsToDuplicateBg, nodeIdsToDuplicateEnd, false); // false= take cell in, even if not all nodes are in dupl // - MCUMesh m0Part2=static_cast(buildPartOfMySelf(cellsAroundGroup->begin(),cellsAroundGroup->end(),true)); - mcIdType nCells2 = m0Part2->getNumberOfCells(); - DAInt desc00=DataArrayIdType::New(),descI00=DataArrayIdType::New(),revDesc00=DataArrayIdType::New(),revDescI00=DataArrayIdType::New(); - MCUMesh m01=m0Part2->buildDescendingConnectivity(desc00,descI00,revDesc00,revDescI00); - - // Neighbor information of the mesh without considering the crack (serves to count how many connex pieces it is made of) - DataArrayIdType *tmp00=0,*tmp11=0; - MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00, tmp00, tmp11); - DAInt neighInit00(tmp00); - DAInt neighIInit00(tmp11); + MCUMesh mAroundGrpLarge=static_cast(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(); + + // 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; + 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 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 = 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(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): - DataArrayIdType *idsTmp=0; - m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp); - DAInt ids(idsTmp); - // 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(ids->begin(),ids->end(),desc00,descI00); - DataArrayIdType *tmp0=0,*tmp1=0; - // Compute the neighbor of each cell in m0Part2, 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. - ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00,tmp0,tmp1); - DAInt neigh00(tmp0); - DAInt neighI00(tmp1); - - // For each initial connex part of the sub-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); - hitCells->fillWithValue(-1); - DAInt cellsToModifyConn0_torenum = DataArrayIdType::New(); - cellsToModifyConn0_torenum->alloc(0,1); + // 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(),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 *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 = 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); - if (!t->getNumberOfTuples()) + DAInt t = hitCells->findIdsEqual(0); + if(!t->getNumberOfTuples()) break; - // Connex zone without the crack (to compute the next seed really) - mcIdType dnu; - DAInt connexCheck = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neighInit00,neighIInit00, -1, dnu); - mcIdType cnt(0); - for (mcIdType * ptr = connexCheck->getPointer(); cnt < connexCheck->getNumberOfTuples(); ptr++, cnt++) - hitCells->setIJ(*ptr,0,1); - // Connex zone WITH the crack (to identify cells lying on either part of the crack) - DAInt spreadZone = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neigh00,neighI00, -1, dnu); - cellsToModifyConn0_torenum = DataArrayIdType::Aggregate(cellsToModifyConn0_torenum, spreadZone, 0); - // Compute next seed, i.e. a cell in another connex part, which was not covered by the previous iterations - DAInt comple = cellsToModifyConn0_torenum->buildComplement(nCells2); - DAInt nonHitCells = hitCells->findIdsEqual(-1); - DAInt intersec = nonHitCells->buildIntersection(comple); - if (intersec->getNumberOfTuples()) - { seed = intersec->getIJ(0,0); } + mcIdType seed = t->getIJ(0,0); + bool done = false; + 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 + { + nIter++; + // Identify connex zone around the seed - this zone corresponds to some cells on the other side + // 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, neigh,neighI, -1, dnu); + done = true; + for(const mcIdType& s: *spreadZone) + { + hitCellsP[s] = PING_FULL; + const auto& it = toOtherSide.find(s); + if (it != toOtherSide.end()) + { + mcIdType other = it->second; + if (hitCellsP[other] != PONG_FULL) + { + // On the other side of the crack we hit a cell which was not fully covered previously by the + // ComputeSpreadZone process, so we are not done yet, ComputeSreadZone will need to be applied there + done = false; + hitCellsP[other] = PONG_PART; + // Compute next seed, i.e. a cell on the other side of the crack + seed = other; + } + } + } + if (done) + { + // we might have several disjoint PONG parts in front of a single PING connex part: + DAInt idsPong = hitCells->findIdsEqual(PONG_PART); + if (idsPong->getNumberOfTuples()) + { + seed = idsPong->getIJ(0,0); + done = false; + } + continue; // continue without switching side (or break if 'done' remains false) + } + else + { + // Go to the other side + std::swap(PING_FULL, PONG_FULL); + std::swap(PING_PART, PONG_PART); + } + } // while (!done ...) + DAInt nonHitCells = hitCells->findIdsEqual(0); + if (nonHitCells->getNumberOfTuples()) + seed = nonHitCells->getIJ(0,0); else - { break; } - nIter++; - } + break; + } // while (nIter < nIterMax ... if (nIter >= nIterMax) - throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate(): internal error - too many iterations."); + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: Too many iterations - should not happen"); + + // 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) (see for example testBuildInnerBoundaryAlongM1Group4) + // 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 + } + } + } + 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 - DAInt cellsToModifyConn1_torenum=cellsToModifyConn0_torenum->buildComplement(neighI00->getNumberOfTuples()-1); - cellsToModifyConn0_torenum->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end()); - cellsToModifyConn1_torenum->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end()); + if (cellsRet1->getNumberOfTuples() + cellsRet2->getNumberOfTuples() != cellsAroundGroupLarge->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: Some cells not hit - Internal error should not happen"); + cellsRet1->transformWithIndArr(cellsAroundGroupLarge->begin(),cellsAroundGroupLarge->end()); + cellsRet2->transformWithIndArr(cellsAroundGroupLarge->begin(),cellsAroundGroupLarge->end()); // - cellIdsNeededToBeRenum=cellsToModifyConn0_torenum.retn(); - cellIdsNotModified=cellsToModifyConn1_torenum.retn(); - nodeIdsToDuplicate=dupl.retn(); + cellIdsNeededToBeRenum=cellsRet1.retn(); + cellIdsNotModified=cellsRet2.retn(); } /*! diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index de9c4da21..2d5f408d7 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -144,7 +144,8 @@ namespace MEDCoupling MEDCOUPLING_EXPORT DataArrayIdType *findCellIdsOnBoundary() const; MEDCOUPLING_EXPORT void findCellIdsLyingOn(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayIdType *&cellIdsRk0, DataArrayIdType *&cellIdsRk1) const; MEDCOUPLING_EXPORT MEDCouplingUMesh *computeSkin() const; - MEDCOUPLING_EXPORT void findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayIdType *& nodeIdsToDuplicate, + MEDCOUPLING_EXPORT DataArrayIdType *findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords) const; + MEDCOUPLING_EXPORT void findCellsToRenumber(const MEDCouplingUMesh& otherDimM1OnSameCoords, const mcIdType *nodeIdsToDuplicateBg, const mcIdType *nodeIdsToDuplicateEnd, DataArrayIdType *& cellIdsNeededToBeRenum, DataArrayIdType *& cellIdsNotModified) const; MEDCOUPLING_EXPORT void duplicateNodes(const mcIdType *nodeIdsToDuplicateBg, const mcIdType *nodeIdsToDuplicateEnd); MEDCOUPLING_EXPORT void renumberNodesWithOffsetInConn(mcIdType offset); diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index ca7e73946..bfe61ebe1 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -369,6 +369,7 @@ typedef long mcPyPtrType; %newobject MEDCoupling::MEDCouplingUMesh::getRenumArrForMEDFileFrmt; %newobject MEDCoupling::MEDCouplingUMesh::convertCellArrayPerGeoType; %newobject MEDCoupling::MEDCouplingUMesh::getRenumArrForConsecutiveCellTypesSpec; +%newobject MEDCoupling::MEDCouplingUMesh::findNodesToDuplicate; %newobject MEDCoupling::MEDCouplingUMesh::buildDirectionVectorField; %newobject MEDCoupling::MEDCouplingUMesh::convertLinearCellsToQuadratic; %newobject MEDCoupling::MEDCouplingUMesh::getEdgeRatioField; @@ -2456,14 +2457,19 @@ namespace MEDCoupling return ret; } - PyObject *findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords) const + DataArrayIdType *findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords) const { - DataArrayIdType *tmp0=0,*tmp1=0,*tmp2=0; - self->findNodesToDuplicate(otherDimM1OnSameCoords,tmp0,tmp1,tmp2); - PyObject *ret=PyTuple_New(3); + DataArrayIdType *ret=self->findNodesToDuplicate(otherDimM1OnSameCoords); + return ret; + } + + PyObject *findCellsToRenumber(const MEDCouplingUMesh& otherDimM1OnSameCoords, const DataArrayIdType *dupNodes) const + { + DataArrayIdType *tmp0=0,*tmp1=0; + self->findCellsToRenumber(otherDimM1OnSameCoords,dupNodes->begin(), dupNodes->end(), tmp0,tmp1); + PyObject *ret=PyTuple_New(2); PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(tmp0),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(tmp1),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); - PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(tmp2),SWIGTITraits::TI, SWIG_POINTER_OWN | 0 )); return ret; } diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index 9ca778472..78460bf1c 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -4085,9 +4085,9 @@ void MEDFileUMesh::optimizeFamilies() * might not be duplicated at all. * After this operation a top-level cell bordering the group will loose some neighbors (typically the cell which is on the * other side of the group is no more a neighbor) - * - finally, the connectivity of (part of) the top level-cells bordering the group is also modified so that some cells + * - the connectivity of (part of) the top level-cells bordering the group is also modified so that some cells * bordering the newly created boundary use the newly computed nodes. - * Finally note that optional cell numbers are also affected by this method and might become invalid for SMESH. + * - finally note that optional cell numbers are also affected by this method and might become invalid for SMESH. * Use clearNodeAndCellNumbers() afterwards to ensure a proper SMESH loading. * * \param[in] grpNameM1 name of the (-1)-level group defining the boundary @@ -4110,11 +4110,16 @@ void MEDFileUMesh::buildInnerBoundaryAlongM1Group(const std::string& grpNameM1, MUMesh m1=getMeshAtLevel(-1); mcIdType nbNodes=m0->getNumberOfNodes(); MUMesh m11=getGroup(-1,grpNameM1); - DataArrayIdType *tmp00=0,*tmp11=0,*tmp22=0; - m0->findNodesToDuplicate(*m11,tmp00,tmp11,tmp22); - DAInt nodeIdsToDuplicate(tmp00); + DataArrayIdType *tmp00=0, *tmp11=0,*tmp22=0; + + // !!! The core of the duplication logic is in these 2 methods: + // !!! + DAInt nodeIdsToDuplicate = m0->findNodesToDuplicate(*m11); // identify nodes to duplicate + m0->findCellsToRenumber(*m11, nodeIdsToDuplicate->begin(), nodeIdsToDuplicate->end(), tmp11,tmp22); // identify cells needing renumbering DAInt cellsToModifyConn0(tmp11); DAInt cellsToModifyConn1(tmp22); + // !!!! + MUMesh tmp0=static_cast(m0->buildPartOfMySelf(cellsToModifyConn0->begin(),cellsToModifyConn0->end(),true)); // node renumbering of cells in m1 impacted by duplication of node but not in group 'grpNameM1' on level -1 DAInt descTmp0=DataArrayIdType::New(),descITmp0=DataArrayIdType::New(),revDescTmp0=DataArrayIdType::New(),revDescITmp0=DataArrayIdType::New(); diff --git a/src/MEDLoader/Swig/MEDLoaderTest3.py b/src/MEDLoader/Swig/MEDLoaderTest3.py index 1c80c68f9..e752d3544 100644 --- a/src/MEDLoader/Swig/MEDLoaderTest3.py +++ b/src/MEDLoader/Swig/MEDLoaderTest3.py @@ -1508,7 +1508,6 @@ class MEDLoaderTest3(unittest.TestCase): self.assertEqual([9,11],mfu.getGroupArr(-1,"group").getValues()) self.assertEqual([23,24],mfu.getGroupArr(-1,"group_dup").getValues()) self.assertEqual([0,1],mfu.getGroupArr(-1,"group2").getValues()) -# mfu.getMeshAtLevel(0).writeVTK("/tmp/mfu_M0.vtu") ref0 =[3, 5, 10, 12, 3, 12, 10, 11, 3, 12, 11, 13] ref1 =[3, 2, 6, 7, 3, 2, 7, 3, 3, 1, 5, 6, 3, 1, 6, 2] self.assertEqual(ref0,mfu.getMeshAtLevel(0)[[3,10,11]].getNodalConnectivity().getValues()) @@ -1625,6 +1624,110 @@ class MEDLoaderTest3(unittest.TestCase): m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999) pass + @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, 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), + (5,-3,10),(5,-3,20),(5,-3,30),(5,-3,40),(5,2,10),(5,2,20),(5,2,30),(5,2,40),(0,-3,10),(0,-3,20),(0,-3,30),(0,-3,40),(0,2,10),(0,2,20),(0,2,30),(0,2,40),(20,8,0),(0,8,0),(20,8,10),(20,8,20),(20,8,30),(20,8,40),(15,8,30),(15,8,40),(5,8,30),(5,8,40),(0,8,10),(0,8,20),(0,8,30),(0,8,40)]) + m3.setCoords(coo) + c = DataArrayInt([31, 0, 3, 2, 1, -1, 16, 20, 24, 28, -1, 0, 16, 28, 3, -1, 3, 28, 24, 2, -1, 2, 24, 20, 1, -1, 1, 20, 16, 0, 31, 16, 28, 24, 20, -1, 17, 21, 25, 29, -1, 16, 17, 29, 28, -1, 28, 29, 25, 24, -1, 24, 25, 21, 20, -1, 20, 21, 17, 16, 31, 17, 29, 25, 21, -1, 18, 22, 26, 30, -1, 17, 18, 30, 29, -1, 29, 30, 26, 25, -1, 25, 26, 22, 21, -1, 21, 22, 18, 17, 31, 18, 30, 26, 22, -1, 19, 23, 27, 31, -1, 18, 19, 31, 30, -1, 30, 31, 27, 26, -1, 26, 27, 23, 22, -1, 22, 23, 19, 18, 31, 4, 5, 3, 0, -1, 32, 16, 28, 36, -1, 4, 32, 36, 5, -1, 5, 36, 28, 3, -1, 3, 28, 16, 0, -1, 0, 16, 32, 4, 31, 32, 36, 28, 16, -1, 33, 17, 29, 37, -1, 32, 33, 37, + 36, -1, 36, 37, 29, 28, -1, 28, 29, 17, 16, -1, 16, 17, 33, 32, 31, 33, 37, 29, 17, -1, 34, 18, 30, 38, -1, 33, 34, 38, 37, -1, 37, 38, 30, 29, -1, 29, 30, 18, 17, -1, 17, 18, 34, 33, 31, 34, 38, 30, 18, -1, 35, 19, 31, 39, -1, 34, 35, 39, 38, -1, 38, 39, 31, 30, -1, 30, 31, 19, 18, -1, 18, 19, 35, 34, 31, 6, 5, 4, 7, -1, 40, 44, 32, 36, -1, 6, 40, 36, 5, -1, 5, 36, 32, 4, -1, 4, 32, 44, 7, -1, 7, 44, 40, 6, 31, 40, 36, 32, 44, -1, 41, 45, 33, 37, -1, 40, 41, 37, 36, -1, 36, 37, 33, 32, -1, 32, 33, 45, 44, -1, 44, 45, 41, 40, 31, 41, 37, 33, 45, -1, 42, 46, 34, 38, -1, 41, 42, 38, 37, -1, 37, 38, 34, 33, -1, 33, 34, 46, 45, -1, 45, 46, 42, 41, 31, + 42, 38, 34, 46, -1, 43, 47, 35, 39, -1, 42, 43, 39, 38, -1, 38, 39, 35, 34, -1, 34, 35, 47, 46, -1, 46, 47, 43, 42, 31, 80, 9, 5, 6, -1, 82, 40, 36, 52, -1, 80, 82, 52, 9, -1, 9, 52, 36, 5, -1, 5, 36, 40, 6, -1, 6, 40, 82, 80, 31, 82, 52, 36, 40, -1, 83, 41, 37, 53, -1, 82, 83, 53, 52, -1, 52, 53, 37, 36, -1, 36, 37, 41, 40, -1, 40, 41, 83, 82, 31, 83, 53, 37, 41, -1, 84, 42, 38, 86, -1, 83, 84, 86, 53, -1, 53, 86, 38, 37, -1, 37, 38, 42, 41, -1, 41, 42, 84, 83, 31, 84, 86, 38, 42, -1, 85, 43, 39, 87, -1, 84, 85, 87, 86, -1, 86, 87, 39, 38, -1, 38, 39, 43, 42, -1, 42, 43, 85, 84, 31, 10, 9, 8, 11, -1, 56, 60, 48, 52, -1, 10, 56, 52, 9, -1, 9, 52, + 48, 8, -1, 8, 48, 60, 11, -1, 11, 60, 56, 10, 31, 56, 52, + 48, 60, -1, 57, 61, 49, 53, -1, 56, 57, 53, 52, -1, 52, 53, 49, 48, -1, 48, 49, 61, 60, -1, 60, 61, 57, 56, 31, 57, 53, 49, 61, -1, 58, 62, 50, 54, -1, 57, 58, 54, 53, -1, 53, 54, 50, 49, -1, 49, 50, 62, 61, -1, 61, 62, 58, 57, 31, 58, 54, 50, 62, -1, 59, 63, 51, 55, -1, 58, 59, 55, 54, -1, 54, 55, 51, 50, -1, 50, 51, 63, 62, -1, 62, 63, 59, 58, 31, 12, 13, 9, 10, -1, 64, 56, 52, 68, -1, 12, 64, 68, 13, -1, 13, 68, 52, 9, -1, 9, 52, 56, 10, -1, 10, 56, 64, 12, 31, 64, 68, 52, 56, -1, 65, 57, 53, 69, -1, 64, 65, 69, 68, -1, 68, 69, 53, 52, -1, 52, 53, 57, 56, -1, 56, 57, 65, 64, 31, 65, 69, 53, 57, -1, 66, 58, 54, 70, -1, 65, 66, 70, 69, -1, 69, 70, + 54, 53, -1, 53, 54, 58, 57, -1, 57, 58, 66, 65, 31, 66, 70, 54, 58, -1, 67, 59, 55, 71, -1, 66, 67, 71, 70, -1, 70, 71, 55, 54, -1, 54, 55, 59, 58, -1, 58, 59, 67, 66, 31, 14, 15, 13, 12, -1, 72, 64, 68, 76, -1, 14, 72, 76, 15, -1, 15, 76, 68, 13, -1, 13, 68, 64, 12, -1, 12, 64, 72, 14, 31, 72, 76, 68, 64, -1, 73, 65, 69, 77, -1, 72, 73, 77, 76, -1, 76, 77, 69, 68, -1, 68, 69, 65, 64, -1, 64, 65, 73, 72, 31, 73, 77, 69, 65, -1, 74, 66, 70, 78, -1, 73, 74, 78, 77, -1, 77, 78, 70, 69, -1, 69, 70, 66, 65, -1, 65, 66, 74, 73, 31, 74, 78, 70, 66, -1, 75, 67, 71, 79, -1, 74, 75, 79, 78, -1, 78, 79, 71, 70, -1, 70, 71, 67, 66, -1, + 66, 67, 75, 74, 31, 2, 3, 13, 81, -1, 24, 90, 68, 28, -1, 2, 24, 28, 3, -1, 3, 28, 68, 13, -1, 13, 68, 90, 81, -1, 81, 90, 24, 2, 31, 24, 28, 68, 90, -1, 25, 91, 69, 29, -1, 24, 25, 29, 28, -1, 28, 29, 69, 68, -1, 68, 69, 91, 90, -1, 90, 91, 25, 24, 31, 25, 29, 69, 91, -1, 26, 92, 88, 30, -1, 25, 26, 30, 29, -1, 29, 30, 88, 69, -1, 69, 88, 92, 91, -1, 91, 92, 26, 25, 31, 26, 30, 88, 92, -1, 27, 93, 89, 31, -1, 26, 27, 31, 30, -1, 30, 31, 89, 88, -1, 88, 89, 93, 92, -1, 92, 93, 27, 26, 31, 13, 3, 5, 9, -1, 68, 52, 36, 28, -1, 13, 68, 28, 3, -1, 3, 28, 36, 5, -1, 5, 36, 52, 9, -1, 9, 52, 68, 13, 31, 68, 28, 36, 52, -1, 69, 53, 37, 29, -1, 68, 69, 29, + 28, -1, 28, 29, 37, 36, -1, 36, 37, 53, 52, -1, 52, 53, 69, 68, 31, 69, 29, 37, 53, -1, 88, 86, 38, 30, -1, 69, 88, 30, 29, -1, 29, 30, 38, 37, -1, 37, 38, 86, 53, -1, 53, 86, 88, 69, 31, 88, 30, 38, 86, -1, 89, 87, 39, 31, -1, 88, 89, 31, 30, -1, 30, 31, 39, 38, -1, 38, 39, 87, 86, -1, 86, 87, 89, 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]) + m3.setConnectivity(c, cI) + m3.checkConsistency() + m2, _, _, _, _ = m3.buildDescendingConnectivity() + grpIds = DataArrayInt([2,7,12,17,95,99,103,107,129,133,137,141]); 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+22, mfu.getNumberOfNodes()) + self.assertEqual(nNod+22, m3_bis.getNumberOfNodes()) + 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([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 ... + Warning: on the mesh below some points have already been doubled by a previous cut. + """ + 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+23, mfu.getNumberOfNodes()) + self.assertEqual(nNod+23, m3_bis.getNumberOfNodes()) + self.assertEqual(nNod+23, m2_bis.getNumberOfNodes()) + self.assertEqual([5, 15, 16, 35, 36, 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([0, 1, 2, 3, 20, 21, 22, 23, 34, 35, 36, 37, 38, 39]), set(cells1.getValues())) + self.assertEqual(set([4, 5, 6, 7, 42, 43, 44, 45, 46, 47]), 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) -- 2.39.2