-// 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
/*!
* Checks if \a this mesh is well defined. If no exception is thrown by this method,
- * then \a this mesh is most probably is writable, exchangeable and available for all
+ * then \a this mesh is informatically clean, most probably is writable, exchangeable and available for all
* algorithms. <br> In addition to the checks performed by checkConsistencyLight(), this
- * method thoroughly checks the nodal connectivity.
+ * method thoroughly checks the nodal connectivity. For more geometrical checking
+ * checkGeomConsistency method is better than this.
+ *
+ * \sa MEDCouplingUMesh::checkGeomConsistency
+ *
* \param [in] eps - a not used parameter.
* \throw If the mesh dimension is not set.
* \throw If the coordinates array is not set (if mesh dimension != -1 ).
}
}
+/*!
+ * This method adds some geometrical checks in addition to the informatical check of checkConsistency method.
+ * This method in particular checks that a same node is not repeated several times in a cell.
+ *
+ * \throw If there is a presence a multiple same node ID in nodal connectivity of cell.
+ */
+void MEDCouplingUMesh::checkGeomConsistency(double eps) const
+{
+ this->checkConsistency(eps);
+ auto nbOfCells(getNumberOfCells());
+ const mcIdType *ptr(_nodal_connec->begin()),*ptrI(_nodal_connec_index->begin());
+ for(auto icell = 0 ; icell < nbOfCells ; ++icell)
+ {
+ std::set<mcIdType> s(ptr+ptrI[icell]+1,ptr+ptrI[icell+1]);
+ if(ToIdType(s.size())==ptrI[icell+1]-ptrI[icell]-1)
+ continue;
+ std::ostringstream oss; oss << "MEDCouplingUMesh::checkGeomConsistency : for cell #" << icell << " presence of multiple same nodeID !";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+}
+
+
/*!
* Sets dimension of \a this mesh. The mesh dimension in general depends on types of
* elements contained in the mesh. For more info on the mesh dimension see
const mcIdType *endNdlConnOfCurCell=conn+connIndex[eltId+1];
for(const mcIdType *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
if(*iter>=0)//for polyhedrons
- *std::find_if(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1],std::bind2nd(std::equal_to<mcIdType>(),-1))=eltId;
+ *std::find_if(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1],std::bind(std::equal_to<mcIdType>(),std::placeholders::_1,-1))=eltId;
}
}
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)
* \b WARNING this method do the assumption that connectivity lies on the coordinates set.
* For speed reasons no check of this will be done. This method calls
* MEDCouplingUMesh::buildDescendingConnectivity to compute the result.
- * This method lists cell by cell in \b this which are its neighbors. To compute the result
+ * This method lists for every cell in \b this its neighbor \b cells. To compute the result
* only connectivities are considered.
- * The neighbor cells of cell having id 'cellId' are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
+ * The neighbor cells of a given cell having id 'cellId' are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
* The format of return is hence \ref numbering-indirect.
*
* \param [out] neighbors is an array storing all the neighbors of all cells in \b this. This array is newly
ComputeNeighborsOfCellsAdv(desc,descIndx,revDesc,revDescIndx,neighbors,neighborsIndx);
}
-void MEDCouplingUMesh::computeCellNeighborhoodFromNodesOne(const DataArrayIdType *nodeNeigh, const DataArrayIdType *nodeNeighI, MCAuto<DataArrayIdType>& cellNeigh, MCAuto<DataArrayIdType>& cellNeighIndex) const
+/**
+ * Given a set of identifiers indexed by the node IDs of the mesh (and given in the (\ref numbering-indirect format) ,
+ * re-arrange the data to produce a set indexed by cell IDs. The mapping between a node ID and a cell ID is done using the connectivity
+ * of the mesh (e.g. a triangular element will receive the information from its three vertices).
+ * Doublons are eliminated. If present in the inital dataset, the ID of the cell itself is also remooved.
+ *
+ * \param [in] nodeNeigh a set of identifiers (mcIdType) stored by node index (\ref numbering-indirect format)
+ * \param [in] nodeNeighI a set of identifiers (mcIdType) stored by node index (\ref numbering-indirect format)
+ * \param [out] cellNeigh This array is newly allocated and should be dealt by the caller. It contains the initial identifiers
+ * provided in the input parameters but stored now by cell index (See 2nd output parameter and \ref numbering-indirect).
+ * \param [out] cellNeighI is an array of size this->getNumberOfCells()+1 newly allocated and should be
+ * dealt by the caller. This arrays allow to use the first output parameter \b neighbors (\ref numbering-indirect).
+ *
+ * \raise if the number of tuples in nodeNeighI is not equal to the number of nodes in the mesh.
+ */
+void MEDCouplingUMesh::computeCellNeighborhoodFromNodesOne(const DataArrayIdType *nodeNeigh, const DataArrayIdType *nodeNeighI,
+ MCAuto<DataArrayIdType>& cellNeigh, MCAuto<DataArrayIdType>& cellNeighIndex) const
{
if(!nodeNeigh || !nodeNeighI)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::computeCellNeighborhoodFromNodesOne : null pointer !");
{
std::set<mcIdType> s;
for(const mcIdType *it=c+ci[i]+1;it!=c+ci[i+1];it++)
- if(*it>=0)
+ if(*it>=0) // avoid -1 in polygons or polyedrons
s.insert(ne+nei[*it],ne+nei[*it+1]);
s.erase(i);
cellNeigh->insertAtTheEnd(s.begin(),s.end());
* \b WARNING this method do the assumption that connectivity lies on the coordinates set.
* For speed reasons no check of this will be done. This method calls
* MEDCouplingUMesh::buildDescendingConnectivity to compute the result.
- * This method lists node by node in \b this which are its neighbors. To compute the result
+ * This method lists for every node in \b this its neighbor \b nodes. To compute the result
* only connectivities are considered.
* The neighbor nodes of node having id 'nodeId' are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
*
INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[posOfCurCell];
const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
INTERP_KERNEL::NormalizedCellType newType=INTERP_KERNEL::NORM_ERROR;
- mcIdType newLgth;
+ mcIdType newLgth=0;
if(cm.isDynamic())
{
switch(cm.getDimension())
newType=INTERP_KERNEL::CellSimplify::tryToUnPoly3D(zipFullReprOfPolyh,nbOfFaces,lgthOfPolyhConn,conn+newPos+1,newLgth);
break;
}
- case 1:
+ /* case 1: // Not supported yet
{
newType=(lgthOfCurCell==3)?INTERP_KERNEL::NORM_SEG2:INTERP_KERNEL::NORM_POLYL;
break;
}
+ */
}
ret=ret || (newType!=type);
conn[newPos]=newType;
* 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
else
return 0;
}
-
- return work!=tmp+sz1?1:0;
}
else
{//case of SEG2 and SEG3
/*!
- * This method find cells that are equal (regarding \a compType) in \a this. The comparison is specified
- * by \a compType.
- * This method keeps the coordiantes of \a this. This method is time consuming.
+ * This method find cells that are equal (regarding \a compType) in \a this. The comparison is specified by \a compType (see zipConnectivityTraducer).
+ * This method keeps the coordinates of \a this. The comparison starts at rank \a startCellId cell id (included).
+ * If \a startCellId is equal to 0 algorithm starts at cell #0 and for each cell candidates being searched have cell id higher than current cellId.
+ * If \a startCellId is greater than 0 algorithm starts at cell #startCellId but for each cell all candidates are considered.
+ * This method is time consuming.
*
* \param [in] compType input specifying the technique used to compare cells each other.
* - 0 : exactly. A cell is detected to be the same if and only if the connectivity is exactly the same without permutation and types same too. This is the strongest policy.
* \param [in] startCellId specifies the cellId starting from which the equality computation will be carried out. By default it is 0, which it means that all cells in \a this will be scanned.
* \param [out] commonCellsArr common cells ids (\ref numbering-indirect)
* \param [out] commonCellsIArr common cells ids (\ref numbering-indirect)
- * \return the correspondence array old to new in a newly allocated array.
*
*/
void MEDCouplingUMesh::findCommonCells(int compType, mcIdType startCellId, DataArrayIdType *& commonCellsArr, DataArrayIdType *& commonCellsIArr) const
std::vector<bool> isFetched(nbOfCells,false);
if(startCellId==0)
{
- for(mcIdType i=0;i<nbOfCells;i++)
+ for(mcIdType i=startCellId;i<nbOfCells;i++)
{
if(!isFetched[i])
{
- const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind2nd(std::not_equal_to<mcIdType>(),-1));
+ const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,-1));
std::vector<mcIdType> v,v2;
if(connOfNode!=connPtr+connIPtr[i+1])
{
{
if(!isFetched[i])
{
- const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind2nd(std::not_equal_to<mcIdType>(),-1));
+ const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,-1));
+ // v2 contains the result of successive intersections using rev nodal on on each node of cell #i
std::vector<mcIdType> v,v2;
if(connOfNode!=connPtr+connIPtr[i+1])
{
std::vector<mcIdType>::iterator it=std::set_intersection(v.begin(),v.end(),revNodalPtr+revNodalIPtr[*connOfNode],revNodalPtr+revNodalIPtr[*connOfNode+1],v2.begin());
v2.resize(std::distance(v2.begin(),it));
}
+ // v2 contains now candidates. Problem candidates are sorted using id rank.
if(v2.size()>1)
{
+ if(v2[0]!=i)
+ {
+ auto it(std::find(v2.begin(),v2.end(),i));
+ std::swap(*v2.begin(),*it);
+ }
if(AreCellsEqualInPool(v2,compType,connPtr,connIPtr,commonCells))
{
- mcIdType pos=commonCellsI->back();
- commonCellsI->pushBackSilent(commonCells->getNumberOfTuples());
+ mcIdType newPos(commonCells->getNumberOfTuples());
+ mcIdType pos(commonCellsI->back());
+ std::sort(commonCells->getPointerSilent()+pos,commonCells->getPointerSilent()+newPos);
+ commonCellsI->pushBackSilent(newPos);
for(const mcIdType *it=commonCells->begin()+pos;it!=commonCells->end();it++)
isFetched[*it]=true;
}
oss << " !";
throw INTERP_KERNEL::Exception(oss.str());
}
- MCAuto<DataArrayIdType> o2n=mesh->zipConnectivityTraducer(compType,nbOfCells);
- arr=o2n->subArray(nbOfCells);
- arr->setName(other->getName());
- mcIdType tmp;
+ //
if(other->getNumberOfCells()==0)
+ {
+ MCAuto<DataArrayIdType> dftRet(DataArrayIdType::New()); dftRet->alloc(0,1); arr=dftRet.retn(); arr->setName(other->getName());
return true;
- return arr->getMaxValue(tmp)<nbOfCells;
+ }
+ DataArrayIdType *commonCells(nullptr),*commonCellsI(nullptr);
+ mesh->findCommonCells(compType,nbOfCells,commonCells,commonCellsI);
+ MCAuto<DataArrayIdType> commonCellsTmp(commonCells),commonCellsITmp(commonCellsI);
+ mcIdType newNbOfCells=-1;
+ MCAuto<DataArrayIdType> o2n = DataArrayIdType::ConvertIndexArrayToO2N(ToIdType(mesh->getNumberOfCells()),commonCells->begin(),commonCellsI->begin(),commonCellsI->end(),newNbOfCells);
+ MCAuto<DataArrayIdType> p0(o2n->selectByTupleIdSafeSlice(0,nbOfCells,1));
+ mcIdType maxPart(p0->getMaxValueInArray());
+ bool ret(maxPart==newNbOfCells-1);
+ MCAuto<DataArrayIdType> p1(p0->invertArrayO2N2N2O(newNbOfCells));
+ // fill p1 array in case of presence of cells in other not in this
+ mcIdType *pt(p1->getPointer());
+ for(mcIdType i = maxPart ; i < newNbOfCells-1 ; ++i )
+ pt[i+1] = i+1;
+ //
+ MCAuto<DataArrayIdType> p2(o2n->subArray(nbOfCells));
+ p2->transformWithIndArr(p1->begin(),p1->end()); p2->setName(other->getName());
+ arr = p2.retn();
+ return ret;
}
/*!
* By default coordinates are kept. This method is close to MEDCouplingUMesh::buildPartOfMySelf except that here input
* cellIds is not given explicitly but by a range python like.
*
- * \param start
- * \param end
- * \param step
+ * \param start starting ID
+ * \param end end ID (excluded)
+ * \param step step size
* \param keepCoords that specifies if you want or not to keep coords as this or zip it (see MEDCoupling::MEDCouplingUMesh::zipCoords). If true zipCoords is \b NOT called, if false, zipCoords is called.
* \return a newly allocated
*
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();
}
/*!
* \throw if \b otherDimM1OnSameCoords is not part of constituent of \b this, or if coordinate pointer of \b this and \b otherDimM1OnSameCoords
* are not same, or if this->getMeshDimension()-1!=otherDimM1OnSameCoords.getMeshDimension()
*
- * \param [in] otherDimM1OnSameCoords
+ * \param [in] otherDimM1OnSameCoords other mesh
* \param [out] cellIdsRk0 a newly allocated array containing the cell ids of s0 (which are cell ids of \b this) in the above algorithm.
* \param [out] cellIdsRk1 a newly allocated array containing the cell ids of s1 \b indexed into the \b cellIdsRk0 subset. To get the absolute ids of s1, simply invoke
* cellIdsRk1->transformWithIndArr(cellIdsRk0->begin(),cellIdsRk0->end());
/*!
* 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.
*
- * \warning This method modifies param \b otherDimM1OnSameCoords (for speed reasons).
*/
-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);
+ 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);
+ 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())
throw INTERP_KERNEL::Exception("MEDFileUMesh::buildInnerBoundaryAlongM1Group: group is too complex: some points (or edges) have more than two connected segments (or faces)!");
dt0=0; dit0=0; rdt0=0; rdit0=0; idsTmp0=0;
- // Get extreme nodes from the group (they won't be duplicated), ie nodes belonging to boundary cells of M1
+ // 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
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); // 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 still need duplication:
+ // 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);
+ 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())
{
- MCUMesh m1IntersecSkin = static_cast<MEDCouplingUMesh *>(m0descSkinDesc->buildPartOfMySelf(validIds->begin(), validIds->end(), true));
+ // 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<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);
+
+ // 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())
+ {
+ 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.
+ }
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();
}
/*!
* should be contained in[0;this->getNumberOfCells()).
*
* \param [in] old2NewBg is expected to be a dynamically allocated pointer of size at least equal to this->getNumberOfCells()
- * \param check
+ * \param check whether to check content of old2NewBg
*/
void MEDCouplingUMesh::renumberCells(const mcIdType *old2NewBg, bool check)
{
if(pt[ptI[cellId]]!=INTERP_KERNEL::NORM_POLYHED)
return ptI[cellId+1]-ptI[cellId]-1;
else
- return ToIdType(std::count_if(pt+ptI[cellId]+1,pt+ptI[cellId+1],std::bind2nd(std::not_equal_to<mcIdType>(),-1)));
+ return ToIdType(std::count_if(pt+ptI[cellId]+1,pt+ptI[cellId+1],std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,-1)));
}
/*!
/*!
* Second step of serialization process.
* \param tinyInfo must be equal to the result given by getTinySerializationInformation method.
- * \param a1
- * \param a2
- * \param littleStrings
+ * \param a1 DataArrayDouble
+ * \param a2 DataArrayDouble
+ * \param littleStrings string vector
*/
void MEDCouplingUMesh::resizeForUnserialization(const std::vector<mcIdType>& tinyInfo, DataArrayIdType *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
{
area_vol[iel]=INTERP_KERNEL::computeVolSurfOfCell2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(type,connec+ipt+1,connec_index[iel+1]-ipt-1,coords,dim_space);
}
if(isAbs)
- std::transform(area_vol,area_vol+nbelem,area_vol,std::ptr_fun<double,double>(fabs));
+ std::transform(area_vol,area_vol+nbelem,area_vol,[](double c){return fabs(c);});
}
else
{
*area_vol++=INTERP_KERNEL::computeVolSurfOfCell2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(type,connec+ipt+1,connec_index[*iel+1]-ipt-1,coords,dim_space);
}
if(isAbs)
- std::transform(array->getPointer(),area_vol,array->getPointer(),std::ptr_fun<double,double>(fabs));
+ std::transform(array->getPointer(),area_vol,array->getPointer(),[](double c){return fabs(c);});
}
else
{
mcIdType nbNodes=getNumberOfNodes();
MCAuto<DataArrayDouble> nnpc;
{
- MCAuto<DataArrayIdType> tmp(computeNbOfNodesPerCell());
- nnpc=tmp->convertToDblArr();
+ MCAuto<DataArrayIdType> tmp2(computeNbOfNodesPerCell());
+ nnpc=tmp2->convertToDblArr();
}
std::for_each(nnpc->rwBegin(),nnpc->rwEnd(),[](double& v) { v=1./v; });
const double *nnpcPtr(nnpc->begin());
mcIdType offset=connI[i];
INTERP_KERNEL::crossprod<3>(locPtr+3*i,coords+3*conn[offset+1],coords+3*conn[offset+2],vals);
double n=INTERP_KERNEL::norm<3>(vals);
- std::transform(vals,vals+3,vals,std::bind2nd(std::multiplies<double>(),1./n));
+ std::transform(vals,vals+3,vals,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
}
}
else
mcIdType offset=connI[i];
std::transform(coords+2*conn[offset+2],coords+2*conn[offset+2]+2,coords+2*conn[offset+1],tmp,std::minus<double>());
double n=INTERP_KERNEL::norm<2>(tmp);
- std::transform(tmp,tmp+2,tmp,std::bind2nd(std::multiplies<double>(),1./n));
+ std::transform(tmp,tmp+2,tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
*vals++=-tmp[1];
*vals++=tmp[0];
}
mcIdType offset=connI[*i];
INTERP_KERNEL::crossprod<3>(locPtr,coords+3*conn[offset+1],coords+3*conn[offset+2],vals);
double n=INTERP_KERNEL::norm<3>(vals);
- std::transform(vals,vals+3,vals,std::bind2nd(std::multiplies<double>(),1./n));
+ std::transform(vals,vals+3,vals,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
}
}
else
mcIdType offset=connI[*i];
std::transform(coords+2*conn[offset+2],coords+2*conn[offset+2]+2,coords+2*conn[offset+1],tmp,std::minus<double>());
double n=INTERP_KERNEL::norm<2>(tmp);
- std::transform(tmp,tmp+2,tmp,std::bind2nd(std::multiplies<double>(),1./n));
+ std::transform(tmp,tmp+2,tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
*vals++=-tmp[1];
*vals++=tmp[0];
}
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
* The result is put in the array \a idsPerType. In the returned parameter \a code, foreach i \a code[3*i+2] refers (if different from -1) to a location into the \a idsPerType.
* This method has 1 input \a profile and 3 outputs \a code \a idsInPflPerType and \a idsPerType.
*
- * \param [in] profile
+ * \param [in] profile list of IDs constituing the profile
* \param [out] code is a vector of size 3*n where n is the number of different geometric type in \a this \b reduced to the profile \a profile. \a code has exactly the same semantic than in MEDCouplingUMesh::checkTypeConsistencyAndContig method.
* \param [out] idsInPflPerType is a vector of size of different geometric type in the subpart defined by \a profile of \a this ( equal to \a code.size()/3). For each i,
* \a idsInPflPerType[i] stores the tuple ids in \a profile that correspond to the geometric type code[3*i+0]
/*!
* Convert the nodal connectivity of the mesh so that all the cells are of dynamic types (polygon or quadratic
* polygon). This returns the corresponding new nodal connectivity in \ref numbering-indirect format.
- * \param nodalConn
- * \param nodalConnI
+ * \param nodalConn nodal connectivity
+ * \param nodalConnIndex nodal connectivity indices
*/
void MEDCouplingUMesh::convertNodalConnectivityToDynamicGeoTypeMesh(DataArrayIdType *&nodalConn, DataArrayIdType *&nodalConnIndex) const
{
for(mcIdType j=0;j<code[3*i+1];j++)
*idsPtr++=offset+j;
else
- idsPtr=std::transform(idsPerGeoTypeBg,idsPerGeoTypeEnd,idsPtr,std::bind2nd(std::plus<mcIdType>(),offset));
+ idsPtr=std::transform(idsPerGeoTypeBg,idsPerGeoTypeEnd,idsPtr,std::bind(std::plus<mcIdType>(),std::placeholders::_1,offset));
offset+=code[3*i+1];
}
MCAuto<MEDCouplingUMesh> ret=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(idsTokeep->begin(),idsTokeep->end(),true));
}
mcIdType nbOfNodesInCell=nodalI[i+1]-nodalI[i]-1;
if(nbOfNodesInCell>0)
- std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)nbOfNodesInCell));
+ std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind(std::multiplies<double>(),std::placeholders::_1,1./(double)nbOfNodesInCell));
else
{
std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of cell with no nodes !";
}
}
if(!s.empty())
- std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)s.size()));
+ std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind(std::multiplies<double>(),std::placeholders::_1,1./(double)s.size()));
else
{
std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on polyhedron cell #" << i << " there are no nodes !";
for(mcIdType offset=nodalI[0]+1;offset<nodalI[1];offset++)
std::transform(coor+3*nodal[offset],coor+3*(nodal[offset]+1),dd,dd,std::plus<double>());
mcIdType nbOfNodesInCell(nodalI[1]-nodalI[0]-1);
- std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies<double>(),1./(double)nbOfNodesInCell));
+ std::transform(dd,dd+3,dd,std::bind(std::multiplies<double>(),std::placeholders::_1,1./(double)nbOfNodesInCell));
std::copy(dd,dd+3,matrix+4*2);
INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
mcIdType meshLgth2=(*iter)->getNodalConnectivityArrayLen();
nodalPtr=std::copy(nod,nod+meshLgth2,nodalPtr);
if(iter!=meshes.begin())
- nodalIndexPtr=std::transform(index+1,index+nbOfCells+1,nodalIndexPtr,std::bind2nd(std::plus<mcIdType>(),offset));
+ nodalIndexPtr=std::transform(index+1,index+nbOfCells+1,nodalIndexPtr,std::bind(std::plus<mcIdType>(),std::placeholders::_1,offset));
else
nodalIndexPtr=std::copy(index,index+nbOfCells+1,nodalIndexPtr);
offset+=meshLgth2;
*
* \param [in] eps is a relative precision that allows to establish if some 3D plane are coplanar or not.
* \param [in] coords the coordinates with nb of components exactly equal to 3
- * \param [in] begin begin of the nodal connectivity (geometric type included) of a single polyhedron cell
- * \param [in] end end of nodal connectivity of a single polyhedron cell (excluded)
+ * \param [in] index begin of the nodal connectivity (geometric type included) of a single polyhedron cell
* \param [out] res the result is put at the end of the vector without any alteration of the data.
*/
void MEDCouplingUMesh::SimplifyPolyhedronCell(double eps, const DataArrayDouble *coords, mcIdType index, DataArrayIdType *res, MEDCouplingUMesh *faces,
nbOfEdgesInFace=std::distance(bgFace,endFace);
if(!isPerm[i])
{
- bool b;
+ bool b=false;
for(std::size_t j=0;j<nbOfEdgesInFace;j++)
{
std::pair<mcIdType,mcIdType> p1(bgFace[j],bgFace[(j+1)%nbOfEdgesInFace]);
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;
+}
