-// Copyright (C) 2007-2020 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2023 CEA, EDF
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
return buildDescendingConnectivityGen<MinusOneSonsGenerator>(desc,descIndx,revDesc,revDescIndx,MEDCouplingFastNbrer);
}
+/*!
+ * This method is a generalization of buildDescendingConnectivity method. This method return mesh containing subcells of this at level specified by \a targetDeltaLevel
+ * and return descending and reverse descending correspondances to this.
+ *
+ * \param [in] targetDeltaLevel target delta level compared to \a this mesh dimension. This parameter is expected to be lower than zero.
+ *
+ * \throw If targetDeltaLevel is greater or equal to zero
+ * \throw If targetDeltaLevel is lower than -meshDim
+ * \sa MEDCouplingUMesh::buildDescendingConnectivity, MEDCouplingUMesh::explode3DMeshTo1D
+ */
+MCAuto<MEDCouplingUMesh> MEDCouplingUMesh::explodeMeshTo(int targetDeltaLevel, MCAuto<DataArrayIdType>& desc, MCAuto<DataArrayIdType>& descIndx, MCAuto<DataArrayIdType>& revDesc, MCAuto<DataArrayIdType>& revDescIndx) const
+{
+ this->checkConsistencyLight();
+ if( targetDeltaLevel >= 0 )
+ THROW_IK_EXCEPTION("Input parameter targetDeltaLevel is expected to be lower than zero !");
+ if( targetDeltaLevel == -1 )
+ {
+ desc = DataArrayIdType::New(); descIndx = DataArrayIdType::New(); revDesc = DataArrayIdType::New(); revDescIndx = DataArrayIdType::New();
+ MCAuto<MEDCouplingUMesh> ret( this->buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx) );
+ return ret;
+ }
+ if( targetDeltaLevel == -2 && this->getMeshDimension() == 3 )
+ {
+ desc = DataArrayIdType::New(); descIndx = DataArrayIdType::New(); revDesc = DataArrayIdType::New(); revDescIndx = DataArrayIdType::New();
+ MCAuto<MEDCouplingUMesh> ret( this->explode3DMeshTo1D(desc,descIndx,revDesc,revDescIndx) );
+ return ret;
+ }
+ if( targetDeltaLevel == -this->getMeshDimension() )
+ {
+ MCAuto<MEDCouplingUMesh> ret = MEDCouplingUMesh::Build0DMeshFromCoords( const_cast<DataArrayDouble *>( this->getCoords() ) );
+ MEDCouplingUMesh::DeleteCellTypeInIndexedArray(getNodalConnectivity(),getNodalConnectivityIndex(),desc,descIndx);
+ revDesc = DataArrayIdType::New(); revDescIndx = DataArrayIdType::New();
+ this->getReverseNodalConnectivity(revDesc,revDescIndx);
+ return ret;
+ }
+ THROW_IK_EXCEPTION("Not valid input targetDeltaLevel regarding mesh dimension");
+}
+
/*!
* \a this has to have a mesh dimension equal to 3. If it is not the case an INTERP_KERNEL::Exception will be thrown.
* This behaves exactly as MEDCouplingUMesh::buildDescendingConnectivity does except that this method compute directly the transition from mesh dimension 3 to sub edges (dimension 1)
* This method performs operation only on polyhedrons in \b this. If no polyhedrons exists in \b this, \b this remains unchanged.
* This method allows to merge if any coplanar 3DSurf cells that may appear in some polyhedrons cells.
*
+ * \b WARNING: this method will not modify edges connectivity! Take a look at colinearizeEdges for that.
+ *
* \param [in] eps is a relative precision that allows to establish if some 3D plane are coplanar or not. This epsilon is used to recenter around origin to have maximal
* precision.
*/
setConnectivity(connNew,connINew,false);
}
+/*!
+ * This method expects that spaceDimension is equal to 3 and meshDimension equal to 3.
+ * This method performs operation only on polyhedrons in \b this. If no polyhedrons exists in \b this, \b this remains unchanged.
+ * This method allows to simplify edges of polyhedron cells so that consecutive colinear segments (with intermediate points
+ * not used by any other cell) are merged together.
+ *
+ * \param [in] eps is a relative precision that allows to establish if two consecutive 3D segments are colinear or not.
+ *
+ * \sa simplifyPolyhedra
+ */
+void MEDCouplingUMesh::colinearizeEdges(double eps)
+{
+ //
+ // Thanks to Antoine Gerschenfeld (CEA) for contributing this method!
+ //
+ using DAI = MCAuto<DataArrayIdType>;
+ checkFullyDefined();
+ if(getMeshDimension()!=3 || getSpaceDimension()!=3)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearizeEdges() : works with meshdim=3 and spaceDim=3!");
+ double seps = sqrt(1-eps);
+ // Computing connectivities and correspondances : elements -> segments -> points
+ DAI E_Fi(DataArrayIdType::New()), E_F(DataArrayIdType::New()), F_Ei(DataArrayIdType::New()), F_E(DataArrayIdType::New()),
+ F_Si(DataArrayIdType::New()), F_S(DataArrayIdType::New()), S_Fi(DataArrayIdType::New()), S_F(DataArrayIdType::New()),
+ S_Pi(DataArrayIdType::New()), S_P(DataArrayIdType::New()), P_Si(DataArrayIdType::New()), P_S(DataArrayIdType::New());
+ MCAuto<MEDCouplingUMesh> m_f(buildDescendingConnectivity(E_F, E_Fi, F_E, F_Ei)),
+ m_s(m_f->buildDescendingConnectivity(F_S, F_Si, S_F, S_Fi)),
+ m_p(m_s->buildDescendingConnectivity(S_P, S_Pi, P_S, P_Si)); // E: elem, F: faces, S: segments (edges), P: points (vertices)
+ const mcIdType *S_Pp(S_P->begin()), *S_Pip(S_Pi->begin()), *P_Sp(P_S->begin()), *P_Sip(P_Si->begin());
+ std::set<mcIdType> pt_rem;
+ const mcIdType *m_pi = m_p->getNodalConnectivityIndex()->begin(),
+ *m_pc = m_p->getNodalConnectivity()->begin();
+ double (*coord)[3] = (double (*)[3]) getCoords()->begin();
+ // Find all points only connected to exaclty 2 segments - they are the candidates for elimination
+ // Note that in 3D this can only happen for polyhedrons (when this happens at all)
+ DAI dsi = P_Si->deltaShiftIndex();
+ DAI cand = dsi->findIdsEqual(2);
+ for (const mcIdType& i: *cand) // i is a point to be potentially eliminated, shared by 2 segs only
+ {
+ double n2[2] = {0., 0.}, scal = 0.; // n2 is a squared norm, scal is a scalar product
+ mcIdType p[2][2]; // p[j][k] is the ID (in the coord array) of the k-th point of the j-th segment
+ for (mcIdType j = 0; j < 2; j++)
+ for (mcIdType k = 0; k < 2; k++)
+ {
+ mcIdType off1 = P_Sip[i] + j; // offset to get ID of the j-th seg (around the i-th point) in the point->seg correspondance
+ mcIdType pt_id = P_Sp[off1] + k; // ID of the k-th point of the j-th seg in the point->seg correspondance
+ mcIdType pt_id2 = S_Pp[S_Pip[pt_id]]; // ID of the point in the point mesh
+ p[j][k] = m_pc[m_pi[pt_id2] + 1]; // Absolute ID, as read from the connectvity (+1 to skip type: NORM_POINT1)
+ // Just for fun, as initially written by Antoine :-)
+ // p[j][k] = m_pc[m_pi[S_P->getIJ(S_Pi->getIJ(P_S->getIJ(P_Si->getIJ(i, 0) + j, 0), 0) + k, 0)] + 1];
+ }
+ // Geometric test on scalar product
+ for (int d = 0; d < 3; d++) // dimension
+ {
+ for (int j = 0; j < 2; j++)
+ n2[j] += std::pow(coord[p[j][1]][d] - coord[p[j][0]][d], 2);
+ scal += (coord[p[1][1]][d] - coord[p[1][0]][d]) * (coord[p[0][1]][d] - coord[p[0][0]][d]);
+ }
+ if (scal * scal > seps * n2[0] * n2[1]) // seps is a sqrt for homogeneity
+ pt_rem.insert(m_pc[m_pi[i] + 1]); // point should be removed
+ }
+ // Clean connectivity by filtering points to be removed:
+ DataArrayIdType *old_index = getNodalConnectivityIndex(), *old_conn = getNodalConnectivity();
+ DAI new_index(DataArrayIdType::New()), new_conn(DataArrayIdType::New());
+ const mcIdType *old_index_p(old_index->begin()), *old_conn_p(old_conn->begin());
+ for (mcIdType i = 0; i < getNumberOfCells(); i++)
+ {
+ new_index->pushBackSilent(new_conn->getNbOfElems());
+ for (mcIdType j = old_index_p[i]; j < old_index_p[i + 1]; j++)
+ {
+ // Keep point if it is not to be removed, or if is first in connectivity (TODO this last check could be removed?)
+ if (std::find(pt_rem.begin(), pt_rem.end(), old_conn_p[j]) == pt_rem.end() || j == old_index_p[i])
+ new_conn->pushBackSilent(old_conn_p[j]);
+ }
+ }
+ new_index->pushBackSilent(new_conn->getNbOfElems());
+ setConnectivity(new_conn, new_index);
+}
+
/*!
* This method returns all node ids used in the connectivity of \b this. The data array returned has to be dealt by the caller.
* The returned node ids are sorted ascendingly. This method is close to MEDCouplingUMesh::getNodeIdsInUse except
DataArrayIdType *MEDCouplingUMesh::findCellIdsOnBoundary() const
{
checkFullyDefined();
- MCAuto<DataArrayIdType> desc=DataArrayIdType::New();
- MCAuto<DataArrayIdType> descIndx=DataArrayIdType::New();
- MCAuto<DataArrayIdType> revDesc=DataArrayIdType::New();
- MCAuto<DataArrayIdType> revDescIndx=DataArrayIdType::New();
+ MCAuto<DataArrayIdType> ret2(DataArrayIdType::New());
+
+ if (getNumberOfCells() == 0)
+ {
+ ret2->alloc(0,1);
+ return ret2.retn();
+ }
+
+ MCAuto<DataArrayIdType> desc(DataArrayIdType::New()), descIndx(DataArrayIdType::New()), revDesc(DataArrayIdType::New()), revDescIndx(DataArrayIdType::New());
//
buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx)->decrRef();
desc=(DataArrayIdType*)0; descIndx=(DataArrayIdType*)0;
if(!ret1[revDescPtr[revDescIndxPtr[*pt]]])
{ ret1[revDescPtr[revDescIndxPtr[*pt]]]=true; sz++; }
//
- DataArrayIdType *ret2=DataArrayIdType::New();
ret2->alloc(sz,1);
mcIdType *ret2Ptr=ret2->getPointer();
sz=0;
if(*it)
*ret2Ptr++=sz;
ret2->setName("BoundaryCells");
- return ret2;
+ return ret2.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 searches for nodes needed to be duplicated. These nodes are nodes fetched by \b otherDimM1OnSameCoords which are not part of the boundary of \b otherDimM1OnSameCoords.
- * If a node is in the boundary of \b this \b and in the boundary of \b otherDimM1OnSameCoords this node is considered as needed to be duplicated.
+ * This method searches for nodes needed to be duplicated. Those are the nodes fetched by \b otherDimM1OnSameCoords which are not part of the boundary of \b otherDimM1OnSameCoords.
+ * If a node is in the boundary of \b this \b and in the boundary of \b otherDimM1OnSameCoords, it will be duplicated.
* When the set of node ids \b nodeIdsToDuplicate is computed, cell ids in \b this is searched so that their connectivity includes at least 1 node in \b nodeIdsToDuplicate.
*
- * \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
+ * \param [in] crackingMesh 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& crackingMesh) const
{
- typedef MCAuto<DataArrayIdType> DAInt;
- typedef MCAuto<MEDCouplingUMesh> 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<DataArrayIdType>;
+ using MCUMesh = MCAuto<MEDCouplingUMesh>;
checkFullyDefined();
- otherDimM1OnSameCoords.checkFullyDefined();
- if(getCoords()!=otherDimM1OnSameCoords.getCoords())
+ crackingMesh.checkFullyDefined();
+ if(getCoords()!=crackingMesh.getCoords())
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate : meshes do not share the same coords array !");
- if(otherDimM1OnSameCoords.getMeshDimension()!=getMeshDimension()-1)
+ if(crackingMesh.getMeshDimension()!=getMeshDimension()-1)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findNodesToDuplicate : the mesh given in other parameter must have this->getMeshDimension()-1 !");
+ // Clean the M1 group (cracking mesh): the M1 cells which are part of M0 boundary are irrelevant (we can't create a crack on the boundary of M0!)
+ MCUMesh m0skin = computeSkin();
+ DataArrayIdType *idsToKeepP;
+ m0skin->areCellsIncludedIn(&crackingMesh,2, idsToKeepP);
+ DAInt idsToKeep(idsToKeepP);
+ DAInt ids2 = idsToKeep->findIdsNotInRange(0, m0skin->getNumberOfCells()); // discard cells on the skin of M0
+ MCUMesh otherDimM1OnSameCoords =static_cast<MEDCouplingUMesh *>(crackingMesh.buildPartOfMySelf(ids2->begin(), ids2->end(), true));
+
+ if (!otherDimM1OnSameCoords->getNumberOfCells())
+ return MCAuto<DataArrayIdType>(DataArrayIdType::New()).retn();
+
// Checking star-shaped M1 group:
DAInt dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New();
- MCUMesh meshM2 = otherDimM1OnSameCoords.buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); // 2D: a mesh of points, 3D: a mesh of segs
+ MCUMesh meshM2 = otherDimM1OnSameCoords->buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); // 2D: a mesh of points, 3D: a mesh of segs
DAInt dsi = rdit0->deltaShiftIndex();
DAInt 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())
DAInt xtremIdsM2 = dsi->findIdsEqual(1); dsi = 0;
MCUMesh meshM2Part = static_cast<MEDCouplingUMesh *>(meshM2->buildPartOfMySelf(xtremIdsM2->begin(), xtremIdsM2->end(),true));
DAInt xtrem = meshM2Part->computeFetchedNodeIds();
- // 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();
- 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();
- // 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.
+ // Remove from the list points on the boundary of the M0 mesh (those need duplication!).
+ // 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.
+ DAInt fNodes = m0skin->computeFetchedNodeIds();
DAInt notDup = 0;
if (getMeshDimension() == 3)
{
DAInt dnu1=DataArrayIdType::New(), dnu2=DataArrayIdType::New(), dnu3=DataArrayIdType::New(), dnu4=DataArrayIdType::New();
- MCUMesh m0descSkinDesc = m0descSkin->buildDescendingConnectivity(dnu1, dnu2, dnu3, dnu4); // all segments of the skin of the 3D (M0) mesh
+ MCUMesh m0skinDesc = m0skin->buildDescendingConnectivity(dnu1, dnu2, dnu3, dnu4); // all segments of the skin of the 3D (M0) mesh
dnu1=0;dnu2=0;dnu3=0;dnu4=0;
DataArrayIdType * corresp=0;
- meshM2->areCellsIncludedIn(m0descSkinDesc,2,corresp);
+ meshM2->areCellsIncludedIn(m0skinDesc,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-shaped 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:
- MCUMesh m1IntersecSkin = static_cast<MEDCouplingUMesh *>(m0descSkinDesc->buildPartOfMySelf(validIds->begin(), validIds->end(), true));
+ // (the U-shaped polyline described above)
+ MCUMesh m1IntersecSkin = static_cast<MEDCouplingUMesh *>(m0skinDesc->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<MEDCouplingUMesh *>(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<MEDCouplingUMesh *>(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 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();
+}
- //
- MCUMesh m0Part2=static_cast<MEDCouplingUMesh *>(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);
- // 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);
- while (nIter < nIterMax)
- {
- DAInt t = hitCells->findIdsEqual(-1);
- 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); }
- else
- { break; }
- nIter++;
+
+/*!
+ * 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.
+ * \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 !");
+
+ // Degenerated case - no nodes to duplicate
+ if (nodeIdsToDuplicateBg == nodeIdsToDuplicateEnd)
+ {
+ cellIdsNeededToBeRenum = DataArrayIdType::New(); cellIdsNeededToBeRenum->alloc(0,1);
+ cellIdsNotModified = DataArrayIdType::New(); cellIdsNotModified->alloc(0,1);
+ return;
+ }
+
+ // Compute cell IDs of the mesh with cells that touch the M1 group with a least one node:
+ DAInt cellsAroundGroupLarge = getCellIdsLyingOnNodes(nodeIdsToDuplicateBg, nodeIdsToDuplicateEnd, false); // false= take cell in, even if not all nodes are in dupl
+ MCUMesh mAroundGrpLarge=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(cellsAroundGroupLarge->begin(),cellsAroundGroupLarge->end(),true));
+ mcIdType nCellsLarge=cellsAroundGroupLarge->getNumberOfTuples();
+ 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();
+ DataArrayIdType *idsOfM1t;
+ mArGrpLargeDesc->areCellsIncludedIn(&otherDimM1OnSameCoords,2, idsOfM1t);
+ DAInt idsOfM1Large(idsOfM1t);
+ mcIdType nL = mArGrpLargeDesc->getNumberOfCells();
+
+ // Computation of the neighbor information of the mesh WITH the crack (some neighbor links 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):
+ DAInt descLTrunc = descL->deepCopy(), descILTrunc = descIL->deepCopy();
+ DataArrayIdType::RemoveIdsFromIndexedArrays(idsOfM1Large->begin(), idsOfM1Large->end(),descLTrunc,descILTrunc);
+ DataArrayIdType *neight=0, *neighIt=0;
+ MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(descLTrunc,descILTrunc,revDescL,revDescIL, neight, neighIt);
+ DAInt neighL(neight), neighIL(neighIt);
+
+ DAInt hitCellsLarge = DataArrayIdType::New(); hitCellsLarge->alloc(nCellsLarge,1);
+ hitCellsLarge->fillWithValue(0); // 0 : not hit, +1: one side of the crack, -1: other side of the crack,
+ mcIdType* hitCellsLargeP = hitCellsLarge->rwBegin();
+
+ // Now loop on the faces of the M1 group and fill spread zones on either side of the crack:
+ const mcIdType *revDescILP=revDescIL->begin(), *revDescLP=revDescL->begin();
+ for(const auto& v: *idsOfM1Large)
+ {
+ if (v >= nL) continue; // Keep valid match only - see doc of areCellsIncludedIn()
+ mcIdType idx0 = revDescILP[v];
+ // Retrieve the two cells on either side of the face v of M1:
+ mcIdType c1=revDescLP[idx0], c2=revDescLP[idx0+1];
+ std::map<mcIdType, mcIdType> toOther = {{c1, c2}, {c2, c1}};
+ // Handle the spread zones on the two sides of the crack:
+ for (const auto c: {c1, c2})
+ {
+ if (hitCellsLargeP[c]) continue;
+ // Identify connex zone around this cell - if we find a value already assigned there, use it.
+ mcIdType dnu;
+ DAInt spreadZone = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&c, &c+1, neighL,neighIL, -1, dnu);
+ std::set<mcIdType> sv;
+ for (const mcIdType& s: *spreadZone)
+ if (hitCellsLargeP[s]) sv.insert(hitCellsLargeP[s]);
+ if (sv.size() > 1)
+ // Strange: we find in the same spread zone a +1 and -1 !
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: internal error #0 - conflicting values - should not happen!");
+ // If a valid value was found, use it:
+ mcIdType val = sv.size()==1 ? *sv.begin() : 0;
+ // Hopefully this does not conflict with an potential value on the other side:
+ mcIdType other = toOther[c];
+ if (hitCellsLargeP[other])
+ {
+ if(val && hitCellsLargeP[other] != -val)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: internal error #1 - conflicting values - should not happen!");;
+ // We do not yet have a value, but other side has one. Use it!
+ if(!val) val = -hitCellsLargeP[other];
+ }
+ // Cover first initialisation:
+ if (!val) val = 1;
+ // And finally, fill the current spread zone:
+ for(const mcIdType& s: *spreadZone) hitCellsLargeP[s] = val;
+ }
}
- 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());
+ DAInt cellsRet1 = hitCellsLarge->findIdsEqual(1);
+ DAInt cellsRet2 = hitCellsLarge->findIdsEqual(-1);
+
+ if (cellsRet1->getNumberOfTuples() + cellsRet2->getNumberOfTuples() != cellsAroundGroupLarge->getNumberOfTuples())
+ {
+ DAInt nonHitCells = hitCellsLarge->findIdsEqual(0); // variable kept for debug ...
+ 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();
}
/*!
revDesc2=0; revDescIndx2=0;
MCAuto<MEDCouplingUMesh> mDesc1=mDesc2->buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1);//meshDim==1 spaceDim==3
revDesc1=0; revDescIndx1=0;
- mDesc1->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),true,cellIds1D);
+ //Marking all 1D cells that contained at least one node located on the plane
+ //the intersection between those cells and the plane, which consist of the nodes previously tagged, thus don't need to be computed afterwards
+ //(if said intersection is computed in MEDCouplingUMesh::split3DCurveWithPlane, then we might create additional nodes
+ //due to accuracy errors when the needed nodes already exist)
+ mDesc1->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),false,cellIds1D);
MCAuto<DataArrayIdType> cellIds1DTmp(cellIds1D);
//
std::vector<mcIdType> cut3DCurve(mDesc1->getNumberOfCells(),-2);
* point, for more points getCellsContainingPoints() is recommended as it is
* faster.
* \param [in] pos - array of coordinates of the ball central point.
- * \param [in] eps - ball radius.
+ * \param [in] eps - ball radius (i.e. the precision) relative to max dimension along X, Y or Z of the the considered cell bounding box.
* \param [out] elts - vector returning ids of the found cells. It is cleared
* before inserting ids.
* \throw If the coordinates array is not set.
* X0,Y0,Z0,X1,Y1,Z1,... Size of the array must be \a
* this->getSpaceDimension() * \a nbOfPoints
* \param [in] nbOfPoints - number of points to locate within \a this mesh.
- * \param [in] eps - radius of balls (i.e. the precision).
+ * \param [in] eps - ball radius (i.e. the precision) relative to max dimension along X, Y or Z of the the considered cell bounding box.
* \param [out] elts - vector returning ids of found cells.
* \param [out] eltsIndex - an array, of length \a nbOfPoints + 1,
* dividing cell ids in \a elts into groups each referring to one
ptId = direction ? (ptId1 == prevPointId ? ptId2 : ptId1) : (ptId2 == prevPointId ? ptId1 : ptId2);
if (dsi[ptId] == 1) // hitting the end of the line
break;
+
prevPointId = ptId;
mcIdType seg1 = rD[rDI[ptId]], seg2 = rD[rDI[ptId]+1];
activeSeg = (seg1 == activeSeg) ? seg2 : seg1;
+
+ //for piecewise meshes made up of closed parts
+ bool segmentAlreadyTreated = (std::find(linePiece.begin(), linePiece.end(), activeSeg) != linePiece.end());
+ if(segmentAlreadyTreated)
+ break;
}
}
// Done, save final piece into DA:
// identify next valid start segment (one which is not consumed)
if(!edgeSet.empty())
startSeg = *(edgeSet.begin());
+
}
while (!edgeSet.empty());
return result.retn();
else
return 0;
}
+
+/// @cond INTERNAL
+
+namespace MEDCouplingImpl
+{
+ const mcIdType theUndefID = std::numeric_limits< mcIdType >::max(); //!< undefined cell id
+
+ //================================================================================
+ /*!
+ * \brief Encode a cell id and a mesh index into a code
+ * \param [in] id - cell id
+ * \param [in] iMesh - mesh index [0,1]
+ * \return mcIdType - code
+ */
+ //================================================================================
+
+ mcIdType encodeID( mcIdType id, int iMesh )
+ {
+ return ( id + 1 ) * ( iMesh ? -1 : 1 );
+ }
+ //================================================================================
+ /*!
+ * \brief Return cell id and mesh index by a given id
+ * \param [in] id - code of a cell in a mesh
+ * \param [out] iMesh - returned mesh index
+ * \return mcIdType - cell id
+ */
+ //================================================================================
+
+ mcIdType decodeID( mcIdType id, int& iMesh )
+ {
+ iMesh = ( id < 0 );
+ return std::abs( id ) - 1;
+ }
+
+ //================================================================================
+ /*!
+ * \brief return another face sharing two given nodes of a face edge
+ * \param [in] n0 - 1st node of the edge
+ * \param [in] n1 - 2nd node of the edge
+ * \param [in] inputFaceID - face including \a n0 andf \a n2
+ * \param [in] mesh - object and reference meshes
+ * \param [in] revNodal - reverse nodal connectivity of the two meshes
+ * \param [in] revNodalIndx - index of reverse nodal connectivity of the two meshes
+ * \param [out] facesByEdge - return another face including \a n0 andf \a n2
+ * \param [out] equalFaces - return faces equal to facesByEdge
+ */
+ //================================================================================
+
+ void getFacesOfEdge( mcIdType n0,
+ mcIdType n1,
+ mcIdType inputFaceID,
+ MEDCouplingUMesh* mesh[],
+ MCAuto<DataArrayIdType> revNodal[],
+ MCAuto<DataArrayIdType> revNodalIndx[],
+ std::vector< mcIdType >& facesByEdge,
+ std::vector< mcIdType >& equalFaces)
+ {
+ // find faces sharing the both nodes of edge
+
+ facesByEdge.clear();
+ size_t prevNbF; // nb faces found in 0-th mesh
+ for ( int iM = 0; iM < 2; ++iM )
+ {
+ const mcIdType * revInd = revNodalIndx[ iM ]->begin();
+ const mcIdType * rev = revNodal [ iM ]->begin();
+
+ mcIdType nbRevFaces0 = revInd[ n0 + 1 ] - revInd[ n0 ];
+ mcIdType nbRevFaces1 = revInd[ n1 + 1 ] - revInd[ n1 ];
+
+ prevNbF = facesByEdge.size();
+ facesByEdge.resize( prevNbF + std::max( nbRevFaces0, nbRevFaces1 ));
+
+ auto it = std::set_intersection( rev + revInd[ n0 ],
+ rev + revInd[ n0 ] + nbRevFaces0,
+ rev + revInd[ n1 ],
+ rev + revInd[ n1 ] + nbRevFaces1,
+ facesByEdge.begin() + prevNbF );
+ facesByEdge.resize( it - facesByEdge.begin() );
+ }
+
+ // facesByEdge now contains at least the 'inputFaceID'
+ // check if there are other faces
+
+ size_t nbF = facesByEdge.size();
+ if ( nbF > 1 )
+ {
+ if ( prevNbF > 0 && prevNbF < nbF ) // faces found in both meshes
+ {
+ // remove from facesByEdge equal faces in different meshes
+ const mcIdType *conn [2] = { mesh[0]->getNodalConnectivity()->getConstPointer(),
+ mesh[1]->getNodalConnectivity()->getConstPointer() };
+ const mcIdType *connI[2] = { mesh[0]->getNodalConnectivityIndex()->getConstPointer(),
+ mesh[1]->getNodalConnectivityIndex()->getConstPointer() };
+ for ( size_t i0 = 0; i0 < prevNbF; ++i0 )
+ {
+ if ( facesByEdge[ i0 ] == theUndefID )
+ continue;
+ mcIdType objFaceID = MEDCouplingImpl::encodeID( facesByEdge[ i0 ], 0 );
+ bool isInputFace = ( objFaceID == inputFaceID );
+
+ for ( size_t i1 = prevNbF; i1 < facesByEdge.size(); ++i1 )
+ {
+ if ( facesByEdge[ i1 ] == theUndefID )
+ continue;
+
+ mcIdType f0 = facesByEdge[ i0 ];
+ mcIdType f1 = facesByEdge[ i1 ];
+ size_t nbNodes0 = connI[0][ f0 + 1 ] - connI[0][ f0 ] - 1;
+ size_t nbNodes1 = connI[1][ f1 + 1 ] - connI[1][ f1 ] - 1;
+ if ( nbNodes0 != nbNodes1 )
+ continue;
+
+ const mcIdType * fConn0 = conn[0] + connI[0][ f0 ] + 1;
+ const mcIdType * fConn1 = conn[1] + connI[1][ f1 ] + 1;
+ if ( std::equal( fConn0, fConn0 + nbNodes0, fConn1 ))
+ {
+ // equal faces; remove an object one
+ mcIdType refFaceID = MEDCouplingImpl::encodeID( facesByEdge[ i1 ], 1 );
+ if ( refFaceID == inputFaceID )
+ isInputFace = true;
+
+ if ( std::find( equalFaces.begin(),
+ equalFaces.end(), objFaceID ) == equalFaces.end() )
+ equalFaces.push_back( objFaceID );
+
+ facesByEdge[ i0 ] = theUndefID;
+ if ( isInputFace )
+ facesByEdge[ i1 ] = theUndefID;
+ break;
+ }
+ }
+ if ( isInputFace )
+ facesByEdge[ i0 ] = theUndefID;
+ }
+ }
+ }
+
+ nbF = facesByEdge.size();
+ for ( size_t i = 0; i < facesByEdge.size(); ++i )
+ {
+ if ( facesByEdge[ i ] != theUndefID )
+ {
+ facesByEdge[ i ] = MEDCouplingImpl::encodeID( facesByEdge[ i ], i >= prevNbF );
+ if ( facesByEdge[ i ] == inputFaceID )
+ facesByEdge[ i ] = theUndefID;
+ }
+ nbF -= ( facesByEdge[ i ] == theUndefID );
+ }
+
+ if ( nbF > 1 )
+ return; // non-manifold
+
+ if ( nbF < 1 )
+ {
+ facesByEdge.clear();
+ }
+ else // nbF == 1, set a found face first
+ {
+ if ( facesByEdge[ 0 ] == theUndefID )
+ {
+ for ( size_t i = 1; i < facesByEdge.size(); ++i )
+ if ( facesByEdge[ i ] != theUndefID )
+ {
+ facesByEdge[ 0 ] = facesByEdge[ i ];
+ break;
+ }
+ }
+ facesByEdge.resize( 1 );
+ }
+ return;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Remove a face from nodal reversed connectivity
+ * \param [in] node - a node of the face
+ * \param [in] face - the face
+ * \param [in.out] revNodal - reverse nodal connectivity
+ * \param [in,out] revNodalIndx - reverse nodal connectivity index
+ */
+ //================================================================================
+
+ void removeFromRevNodal( mcIdType node,
+ mcIdType face,
+ MCAuto<DataArrayIdType>& revNodal,
+ MCAuto<DataArrayIdType>& revNodalIndx)
+ {
+ mcIdType* fBeg = revNodal->getPointer() + revNodalIndx->getIJ( node, 0 );
+ mcIdType* fEnd = revNodal->getPointer() + revNodalIndx->getIJ( node + 1, 0);
+ auto it = std::find( fBeg, fEnd, face );
+ if ( it != fEnd )
+ {
+ for ( auto it2 = it + 1; it2 < fEnd; ++it2 ) // keep faces sorted
+ *( it2 - 1 ) = *it2;
+
+ *( fEnd - 1 ) = theUndefID;
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Check order of two nodes in a given face
+ * \param [inout] n0 - node 1
+ * \param [inout] n1 - node 2
+ * \param [inout] iFEnc - face
+ * \param [inout] mesh - mesh
+ * \return bool - true if the nodes are in [ .., n1, n0, ..] order in face
+ */
+ //================================================================================
+
+ bool isReverseOrder( mcIdType n0,
+ mcIdType n1,
+ mcIdType iFEnc,
+ MEDCouplingUMesh* mesh[] )
+ {
+ int iMesh;
+ mcIdType iF = decodeID( iFEnc, iMesh );
+
+ const mcIdType *conn = mesh[ iMesh ]->getNodalConnectivity()->getConstPointer();
+ const mcIdType *connI = mesh[ iMesh ]->getNodalConnectivityIndex()->getConstPointer();
+
+ auto it0 = std::find( conn + connI[ iF ] + 1,
+ conn + connI[ iF + 1 ],
+ n0 );
+ auto it1 = std::find( conn + connI[ iF ] + 1,
+ conn + connI[ iF + 1 ],
+ n1 );
+ long i0 = it0 - conn;
+ long i1 = it1 - conn;
+
+ bool isRev = ( std::abs( i1 - i0 ) == 1 ) ? i1 < i0 : i0 < i1;
+ return isRev;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Change orientation of a face in one of given meshes
+ * \param [in] iFEnc - face ID also encoding a mesh index
+ * \param [in,out] mesh - object and reference meshes
+ */
+ //================================================================================
+
+ void reverseFace( mcIdType iFEnc, MEDCouplingUMesh* mesh[] )
+ {
+ int iMesh;
+ mcIdType face = decodeID( iFEnc, iMesh );
+
+ mcIdType *conn = mesh[ iMesh ]->getNodalConnectivity()->getPointer();
+ mcIdType *connI = mesh[ iMesh ]->getNodalConnectivityIndex()->getPointer();
+
+ const INTERP_KERNEL::CellModel& cm =
+ INTERP_KERNEL::CellModel::GetCellModel( mesh[iMesh]->getTypeOfCell( face ));
+
+ cm.changeOrientationOf2D( conn + connI[ face ] + 1,
+ (unsigned int)( connI[ face + 1 ] - connI[ face ] - 1 ));
+ return;
+ }
+}
+
+/// @endcond
+
+//================================================================================
+/*!
+ * \brief Orient cells of \a this 2D mesh equally to \a refFaces
+ * \param [in] refFaces - 2D mesh containing correctly oriented faces. It is optional.
+ * If there are no cells in \a refFaces or it is nullptr, then any face
+ * in \a this mesh is used as a reference
+ * \throw If \a this mesh is not well defined.
+ * \throw If \a this mesh or \refFaces are not 2D.
+ * \throw If \a this mesh and \refFaces do not share nodes.
+ * \throw If \a refFaces are not equally oriented.
+ * \throw If \a this mesh plus \a refFaces together form a non-manifold mesh.
+ *
+ * \if ENABLE_EXAMPLES
+ * \ref cpp_mcumesh_are2DCellsNotCorrectlyOriented "Here is a C++ example".<br>
+ * \ref py_mcumesh_are2DCellsNotCorrectlyOriented "Here is a Python example".
+ * \endif
+ */
+//================================================================================
+
+void MEDCouplingUMesh::orientCorrectly2DCells(const MEDCouplingUMesh* refFaces)
+{
+ checkConsistencyLight();
+ if ( getMeshDimension() != 2 )
+ throw INTERP_KERNEL::Exception("The mesh dimension must be 2");
+ if ( refFaces )
+ {
+ refFaces->checkConsistencyLight();
+ if ( refFaces->getMeshDimension() != 2 )
+ throw INTERP_KERNEL::Exception("The reference mesh dimension must be 2");
+ if ( getCoords() != refFaces->getCoords() )
+ throw INTERP_KERNEL::Exception("Object and reference meshes must share nodes ");
+ if ( refFaces->getNumberOfCells() == 0 )
+ refFaces = nullptr;
+ }
+ if ( getNumberOfCells() == 0 )
+ return;
+
+ enum { _OBJ, _REF };
+ MEDCouplingUMesh* mesh[2] = { this, const_cast< MEDCouplingUMesh* >( refFaces ) };
+ MCAuto<MEDCouplingUMesh> meshPtr;
+ if ( !mesh[_REF] )
+ {
+ meshPtr = mesh[_REF] = MEDCouplingUMesh::New();
+ mesh[_REF]->setCoords( mesh[_OBJ]->getCoords() );
+ mesh[_REF]->allocateCells(0);
+ mesh[_REF]->finishInsertingCells();
+ }
+ mcIdType nbFacesToCheck[2] = { mesh[_OBJ]->getNumberOfCells(),
+ mesh[_REF]->getNumberOfCells() };
+ std::vector< bool > isFaceQueued[ 2 ]; // enqueued faces of 2 meshes
+ isFaceQueued[_OBJ].resize( nbFacesToCheck[_OBJ] );
+ isFaceQueued[_REF].resize( nbFacesToCheck[_REF] );
+
+ MCAuto<DataArrayIdType> revNodal [2] = { DataArrayIdType::New(), DataArrayIdType::New() };
+ MCAuto<DataArrayIdType> revNodalIndx[2] = { DataArrayIdType::New(), DataArrayIdType::New() };
+ mesh[_OBJ]->getReverseNodalConnectivity( revNodal[_OBJ], revNodalIndx[_OBJ] );
+ mesh[_REF]->getReverseNodalConnectivity( revNodal[_REF], revNodalIndx[_REF] );
+
+ std::vector< mcIdType > faceNodes(4);
+ std::vector< mcIdType > facesByEdge(4), equalFaces;
+ std::vector< mcIdType > faceQueue; // starting faces with IDs counted from 1; negative ID mean a face in ref mesh
+
+ while ( nbFacesToCheck[_OBJ] + nbFacesToCheck[_REF] > 0 ) // until all faces checked
+ {
+ if ( faceQueue.empty() ) // all neighbors checked, find more faces to check
+ {
+ for ( int iMesh = 1; iMesh >= 0; --iMesh ) // on [ _REF, _OBJ ]
+ if ( nbFacesToCheck[iMesh] > 0 )
+ for ( mcIdType f = 0, nbF = mesh[iMesh]->getNumberOfCells(); f < nbF; ++f )
+ if ( !isFaceQueued[iMesh][f] )
+ {
+ faceQueue.push_back( MEDCouplingImpl::encodeID( f, iMesh ));
+ isFaceQueued[ iMesh ][ f ] = true;
+ iMesh = 0;
+ break;
+ }
+ if ( faceQueue.empty() )
+ break;
+ }
+
+ mcIdType fID = faceQueue.back();
+ faceQueue.pop_back();
+
+ int iMesh, iMesh2;
+ mcIdType refFace = MEDCouplingImpl::decodeID( fID, iMesh );
+
+ nbFacesToCheck[iMesh]--;
+
+ equalFaces.clear();
+ faceNodes.clear();
+ mesh[iMesh]->getNodeIdsOfCell( refFace, faceNodes );
+ const INTERP_KERNEL::CellModel& cm = INTERP_KERNEL::CellModel::GetCellModel( mesh[iMesh]->getTypeOfCell( refFace ));
+ const int nbEdges = cm.getNumberOfSons();
+
+ // loop on edges of the refFace
+ mcIdType n0 = faceNodes[ nbEdges - 1 ]; // 1st node of edge
+ for ( int edge = 0; edge < nbEdges; ++edge )
+ {
+ mcIdType n1 = faceNodes[ edge ]; // 2nd node of edge
+
+ // get faces sharing the edge
+ MEDCouplingImpl::getFacesOfEdge( n0, n1, fID, mesh, revNodal, revNodalIndx,
+ facesByEdge, equalFaces );
+
+ if ( facesByEdge.size() > 1 )
+ THROW_IK_EXCEPTION("Non-manifold mesh at edge " << n0+1 << " - " << n1+1);
+
+ if ( facesByEdge.size() == 1 )
+ {
+ // compare orientation of two faces
+ //
+ if ( !MEDCouplingImpl::isReverseOrder( n0, n1, facesByEdge[0], mesh ))
+ {
+ if ( facesByEdge[0] < 0 ) // in the ref mesh
+ throw INTERP_KERNEL::Exception("Different orientation of reference faces");
+
+ MEDCouplingImpl::reverseFace( facesByEdge[0], mesh );
+ }
+ mcIdType face2 = MEDCouplingImpl::decodeID( facesByEdge[0], iMesh2 );
+ if ( !isFaceQueued[iMesh2][face2] )
+ {
+ isFaceQueued[iMesh2][face2] = true;
+ faceQueue.push_back( facesByEdge[0] );
+ }
+ }
+ n0 = n1;
+ }
+
+ // remove face and equalFaces from revNodal in order not to treat them again
+ equalFaces.push_back( fID );
+ for ( mcIdType face : equalFaces )
+ {
+ mcIdType f = MEDCouplingImpl::decodeID( face, iMesh2 );
+ const mcIdType *conn = mesh[iMesh2]->getNodalConnectivity()->getConstPointer();
+ const mcIdType *connI = mesh[iMesh2]->getNodalConnectivityIndex()->getConstPointer();
+ mcIdType nbNodes = connI[ f + 1 ] - connI[ f ] - 1;
+ for ( const mcIdType* n = conn + connI[ f ] + 1, *nEnd = n + nbNodes; n < nEnd; ++n )
+
+ MEDCouplingImpl::removeFromRevNodal( *n, f, // not to treat f again
+ revNodal[ iMesh2 ], revNodalIndx[ iMesh2 ] );
+ }
+
+ } // while() until all faces checked
+
+ return;
+}
