-// Copyright (C) 2007-2022 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024 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.
* \return node ids which need to be duplicated following the algorithm explained above.
*
*/
-DataArrayIdType* MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords) const
+DataArrayIdType* MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& crackingMesh) const
{
// DEBUG NOTE: in case of issue with the algorithm in this method, see Python script in resources/dev
// which mimicks the C++
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(); rdit0=0;
- DAInt boundSegs = dsi->findIdsEqual(1); dsi = 0; // boundary segs/faces of the M0 mesh
- MCUMesh m0descSkin = static_cast<MEDCouplingUMesh *>(m0desc->buildPartOfMySelf(boundSegs->begin(),boundSegs->end(), true));
- DAInt fNodes = m0descSkin->computeFetchedNodeIds();
- // 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 shape polyline.
+ // 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:
// (the U-shaped polyline described above)
- MCUMesh m1IntersecSkin = static_cast<MEDCouplingUMesh *>(m0descSkinDesc->buildPartOfMySelf(validIds->begin(), validIds->end(), true));
+ 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);
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);
+ mAroundBNDesc->areCellsIncludedIn(otherDimM1OnSameCoords,2, idsOfM1BNt);
DAInt idsOfM1BN(idsOfM1BNt);
mcIdType nCells=mAroundBN->getNumberOfCells(), nCellsDesc=mAroundBNDesc->getNumberOfCells();
DAInt idsTouch=DataArrayIdType::New(); idsTouch->alloc(0,1);
else // if (3D ...)
notDup = xtrem->buildSubstraction(fNodes);
- DAInt m1Nodes = otherDimM1OnSameCoords.computeFetchedNodeIds();
+ DAInt m1Nodes = otherDimM1OnSameCoords->computeFetchedNodeIds();
DAInt dupl = m1Nodes->buildSubstraction(notDup);
return dupl.retn();
}
* Given a set of nodes to duplicate, this method identifies which cells should have their connectivity modified
* to produce the inner boundary. It is typically called after findNodesToDuplicate().
*
- * \param [in] otherDimM1OnSameCoords a mesh lying on the same coords than \b this and with a mesh dimension equal to those of \b this minus 1. WARNING this input
- * parameter is altered during the call.
+ * \param [in] 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.
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));
if (hitCellsLargeP[other])
{
if(val && hitCellsLargeP[other] != -val)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::findCellsToRenumber: internal error #1 - conflictint values - should not happen!");;
+ 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];
}
* 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
