From cee9ab36de294201b26f81c42bf9df54e6b2eeea Mon Sep 17 00:00:00 2001 From: abn Date: Mon, 15 May 2023 13:47:12 +0200 Subject: [PATCH] New feature: colinearizeEdges(): merge colinear edges of polyhedrons + typically to be used after simplifyPolyhedra() --- src/MEDCoupling/MEDCouplingUMesh.cxx | 80 +++++++++++++++++++ src/MEDCoupling/MEDCouplingUMesh.hxx | 1 + .../MEDCouplingBasicsTest5.py | 27 ++++++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 1 + 4 files changed, 107 insertions(+), 2 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 76d748bd3..b9a9d0304 100755 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -1340,6 +1340,8 @@ bool MEDCouplingUMesh::unPolyze() * 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. */ @@ -1376,6 +1378,84 @@ void MEDCouplingUMesh::simplifyPolyhedra(double eps) 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; + 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 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 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 diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index b47c167f7..3eca47ac7 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -108,6 +108,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT void convertExtrudedPolyhedra(); MEDCOUPLING_EXPORT bool unPolyze(); MEDCOUPLING_EXPORT void simplifyPolyhedra(double eps); + MEDCOUPLING_EXPORT void colinearizeEdges(double eps); MEDCOUPLING_EXPORT MEDCouplingUMesh *buildSpreadZonesWithPoly() const; MEDCOUPLING_EXPORT std::vector partitionBySpreadZone() const; MEDCOUPLING_EXPORT DataArrayIdType *computeFetchedNodeIds() const; diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py index c8864d1ba..86477fee5 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py @@ -2009,14 +2009,14 @@ class MEDCouplingBasicsTest5(unittest.TestCase): # Just to get a nice coords array ... mm = MEDCouplingCMesh(); arr = DataArrayDouble([0.0, 1.0,2.0]) mm.setCoords(arr, arr); mm = mm.buildUnstructured(); coo = mm.getCoords() - + mesh = MEDCouplingUMesh("M", 2) mesh.setCoords(coo) c = [NORM_POLYGON, 0,1,4,7,6,3, NORM_QUAD4, 1,2,5,4, NORM_QUAD4,4,5,8,7] cI = [0, 7,12,17] mm.setConnectivity(DataArrayInt(c),DataArrayInt(cI)) mm.checkConsistencyLight() - + mm.colinearizeKeepingConform2D(eps) c = mm.getNodalConnectivity().getValues() cI = mm.getNodalConnectivityIndex().getValues() @@ -4830,6 +4830,29 @@ class MEDCouplingBasicsTest5(unittest.TestCase): self.assertEqual(cI.getValues(), tgt_cI.getValues()) pass + def testColinearizeEdges(self): + mesh = MEDCouplingUMesh('mesh', 3) + coo = DataArrayDouble([(0,0,0), (1,0,0), (2,0,0), (3,0,0), + (0,0,3), (1,0,3), (2,0,3), (3,0,3), + (0,1,0), (3,1,0), + (0,1,3), (3,1,3)]) + mesh.setCoords(coo) + c = DataArrayInt([NORM_POLYHED, 0,1,2,3,7,6,5,4,-1, # front + 9,8,10,11, -1, # back + 0,4,10,8, -1, # left + 3,7,11,9, -1, # right + 0,1,2,3,9,8,-1, # bottom + 4,5,6,7,11,10 # top + ]) + cI = DataArrayInt([0, c.getNumberOfTuples()]) + mesh.setConnectivity(c, cI) + mesh.colinearizeEdges(1.0e-8) + c, cI = mesh.getNodalConnectivity(), mesh.getNodalConnectivityIndex() + tgt_c = DataArrayInt([NORM_POLYHED, 0, 3, 7, 4, -1, 9, 8, 10, 11, -1, 0, 4, 10, 8, -1, 3, 7, 11, 9, -1, 0, 3, 9, 8, -1, 4, 7, 11, 10]) + tgt_cI = DataArrayInt([0, 30]) + self.assertEqual(c.getValues(), tgt_c.getValues()) + self.assertEqual(cI.getValues(), tgt_cI.getValues()) + pass if __name__ == '__main__': diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index e971e5345..0812ce07e 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -3047,6 +3047,7 @@ namespace MEDCoupling void convertExtrudedPolyhedra(); bool unPolyze(); void simplifyPolyhedra(double eps); + void colinearizeEdges(double eps); MEDCouplingUMesh *buildSpreadZonesWithPoly() const; MEDCouplingUMesh *buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy); }; -- 2.39.2