*
* \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 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.
+ * \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
{
using DAInt = MCAuto<DataArrayIdType>;
using MCUMesh = MCAuto<MEDCouplingUMesh>;
// 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<MEDCouplingUMesh *>(m0desc->buildPartOfMySelf(boundSegs->begin(),boundSegs->end(), true));
DAInt fNodes = m0descSkin->computeFetchedNodeIds();
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 :
+ // - take the sub-mesh of M1 (dim N-1) built keeping only elements touching the m1IntersecSkin
+ // - build its spread zones
+ // - a node common to two different spread zones actually connects two cells of M1 by just one point : this is a singular point
+ // which should not be duplicated.
+ // 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();
+ DAInt boundCells = otherDimM1OnSameCoords.getCellIdsLyingOnNodes(boundNodes->begin(), boundNodes->end(), false); // false= take cell in, even if not all nodes are in dupl
+ MCUMesh m1Bound = static_cast<MEDCouplingUMesh *>(otherDimM1OnSameCoords.buildPartOfMySelf(boundCells->begin(), boundCells->end(), true));
+ DAInt d=DataArrayIdType::New(), dI=DataArrayIdType::New(), revD=DataArrayIdType::New(), revDI=DataArrayIdType::New();
+ MCUMesh dnuM = m1Bound->buildDescendingConnectivity(d,dI,revD,revDI); dnuM = 0;
+ DataArrayIdType *tmp3=0, *tmp4=0;
+ MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(d,dI,revD,revDI, tmp3, tmp4 );
+ DAInt neighSing(tmp3), neighISing(tmp4);
+ mcIdType seed=0;
+ DAInt hitCells = DataArrayIdType::New(); hitCells->alloc(m1Bound->getNumberOfCells(),1); hitCells->fillWithValue(-1);
+ mcIdType nIter=0, nIterMax=m1Bound->getNumberOfCells()+1;
+ while (nIter < nIterMax)
+ {
+ mcIdType dnuId;
+ DAInt sprdZone = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed,&seed+1, neighSing,neighISing, -1, dnuId);
+ const mcIdType * ptr = sprdZone->begin();
+ for (mcIdType cnt = 0; cnt < sprdZone->getNumberOfTuples(); ptr++, cnt++)
+ hitCells->setIJSilent(*ptr,0,nIter);
+ nIter++;
+ DAInt noHit = hitCells->findIdsEqual(-1);
+ if (noHit->getNumberOfTuples())
+ seed = noHit->getIJ(0,0);
+ else
+ break;
+ }
+ if (nIter >= nIterMax)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate(): internal error - too many iterations.");
+ DAInt uniq = hitCells->buildUniqueNotSorted();
+ if (uniq->getNumberOfTuples() > 1)
+ {
+ // All points common to two (or more) different spread zones are singular:
+ DAInt cc=DataArrayIdType::New(), ccI=DataArrayIdType::New();
+ m1Bound->getReverseNodalConnectivity(cc, ccI);
+ DAInt spreadCC = cc->deepCopy();
+ spreadCC->transformWithIndArr(hitCells->begin(), hitCells->end());
+ DAInt singPoints2 = DataArrayIdType::New(); singPoints2->alloc(0,1);
+ const mcIdType* ccIP = ccI->begin();
+ const mcIdType* spreadCCP = spreadCC->begin();
+ for(mcIdType j=0; j < ccI->getNumberOfTuples()-1; j++, ccIP++)
+ {
+ mcIdType sz = *(ccIP+1)-*ccIP;
+ if (sz == 0) continue; // skip empty packs
+ mcIdType val = spreadCCP[*ccIP];
+ for (mcIdType k=1; k < sz; k++)
+ if (val != spreadCCP[*ccIP + k])
+ {
+ singPoints2->pushBackSilent(j);
+ break;
+ }
+ }
+ fNodes1 = fNodes1->buildSubstraction(singPoints2);
+ }
}
- DAInt fNodes2 = fNodes1->buildSubstraction(singPoints);
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 dupl
+ 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
+{
+ using DAInt = MCAuto<DataArrayIdType>;
+ using MCUMesh = MCAuto<MEDCouplingUMesh>;
+
+ 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 cellsAroundGroup = 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();
+ 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);
- // 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);
+ // 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();
+ for (const auto& v: *idsOfM1)
+ {
+ mcIdType idx0 = revDIP[v];
+ mcIdType c1=revDP[idx0], c2=revDP[idx0+1];
+ toOtherSide[c1] = c2;
+ toOtherSide[c2] = c1;
+ }
+
// Neighbor information of the mesh WITH the crack (some neighbors are removed):
- DataArrayIdType *idsTmp=0;
- mArGrpDesc->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;
+ DataArrayIdType::RemoveIdsFromIndexedArrays(idsOfM1->begin(), idsOfM1->end(),desc00,descI00);
// 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.
- ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00,tmp0,tmp1);
- DAInt neigh00(tmp0);
- DAInt neighI00(tmp1);
+ DataArrayIdType *neigh00t=0, *neighI00t=0;
+ MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(desc00,descI00,revDesc00,revDescI00, neigh00t, neighI00t);
+ DAInt neigh00(neigh00t), neighI00(neighI00t);
// For each initial connex part of the M1 mesh (or said differently for each independent crack):
- mcIdType seed = 0, nIter = 0;
+ 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);
+ 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 PING_FULL, PONG_FULL;
while (nIter < nIterMax)
{
DAInt t = hitCells->findIdsEqual(-1);
- if (!t->getNumberOfTuples())
+ 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);
+ mcIdType seed = t->getIJ(0,0);
+ bool done = false;
+ PING_FULL =1; PONG_FULL=2;
+ mcIdType PING_PART=11, PONG_PART=22;
+ // 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, neigh00,neighI00, -1, dnu);
+ done = true;
+ mcIdType* hitCellsP = hitCells->rwBegin();
+ 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(-1);
- DAInt intersec = nonHitCells->buildIntersection(comple);
- if (intersec->getNumberOfTuples())
- { seed = intersec->getIJ(0,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.");
-
- DAInt cellsToModifyConn1_torenum=cellsToModifyConn0_torenum->buildComplement(neighI00->getNumberOfTuples()-1);
- cellsToModifyConn0_torenum->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end());
- cellsToModifyConn1_torenum->transformWithIndArr(cellsAroundGroup->begin(),cellsAroundGroup->end());
+ 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())
+ 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());
//
- cellIdsNeededToBeRenum=cellsToModifyConn0_torenum.retn();
- cellIdsNotModified=cellsToModifyConn1_torenum.retn();
- nodeIdsToDuplicate=dupl.retn();
+ cellIdsNeededToBeRenum=cellsRet1.retn();
+ cellIdsNotModified=cellsRet2.retn();
}
/*!
