From 1ddf936656b21225b59af331cde66dcd5651d4f2 Mon Sep 17 00:00:00 2001 From: abn Date: Fri, 16 May 2014 14:11:03 +0200 Subject: [PATCH] Implemented MEDCouplingUMesh::orderConsecutiveCells1D. Provides a renumbering of the cells of this (which has to be a piecewise connected 1D line), so that the segments of the line are indexed in consecutive order. --- src/MEDCoupling/MEDCouplingUMesh.cxx | 86 ++++++++++++++++++- src/MEDCoupling/MEDCouplingUMesh.hxx | 1 + .../Test/MEDCouplingBasicsTest5.cxx | 53 ++++++++++++ .../Test/MEDCouplingBasicsTest5.hxx | 2 + src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 38 +++++++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 2 + 6 files changed, 180 insertions(+), 2 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 7ff84f831..88f7073ef 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -882,7 +882,7 @@ template MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivityGen(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx, DimM1DescNbrer nbrer) const { if(!desc || !descIndx || !revDesc || !revDescIndx) - throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildDescendingConnectivityGen : present of a null pointer in input !"); + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildDescendingConnectivityGen: presence of a null pointer in input!"); checkConnectivityFullyDefined(); int nbOfCells=getNumberOfCells(); int nbOfNodes=getNumberOfNodes(); @@ -9000,6 +9000,90 @@ void MEDCouplingUMesh::Build2DCellsFrom1DCut(double eps, const MEDCouplingUMesh } } +/** + * Provides a renumbering of the cells of this (which has to be a piecewise connected 1D line), so that + * the segments of the line are indexed in consecutive order (i.e. cells \a i and \a i+1 are neighbors). + * This doesn't modify the mesh. + * The caller is to deal with the resulting DataArrayInt. + * \throw If the coordinate array is not set. + * \throw If the nodal connectivity of the cells is not defined. + * \throw If m1 is not a mesh of dimension 2, or m1 is not a mesh of dimension 1 + * \throw If m2 is not a (piecewise) line (i.e. if a point has more than 2 adjacent segments) + */ +DataArrayInt * MEDCouplingUMesh::orderConsecutiveCells1D() const +{ + checkFullyDefined(); + if(getMeshDimension()!=1 || getSpaceDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::orderConsecutiveCells1D works on unstructured mesh with (meshdim, spacedim) = (1,2)!"); + + // Check that this is a line (and not a more complex 1D mesh) - each point is used at most by 2 segments: + DataArrayInt *_d = DataArrayInt::New(), *_dI = DataArrayInt::New(); + DataArrayInt *_rD = DataArrayInt::New(), *_rDI = DataArrayInt::New(); + MEDCouplingUMesh * m_points; + m_points = buildDescendingConnectivity(_d, _dI, _rD, _rDI); + const int *d=_d->getConstPointer(), *dI=_dI->getConstPointer(); + const int *rD=_rD->getConstPointer(), *rDI=_rDI->getConstPointer(); + MEDCouplingAutoRefCountObjectPtr _dsi = _rDI->deltaShiftIndex(); + const int * dsi=_dsi->getConstPointer(); + MEDCouplingAutoRefCountObjectPtr dsii = _dsi->getIdsNotInRange(0,3); + m_points->decrRef(); + if (dsii->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::orderConsecutiveCells1D only work with a mesh being a (piecewise) connected line!"); + + int nc = getNumberOfCells(); + MEDCouplingAutoRefCountObjectPtr result(DataArrayInt::New()); + result->alloc(nc,1); + + // set of edges not used so far + std::set edgeSet; + for (int i=0; i linePiece; + // fills a list of consecutive segment linked to startSeg. This can go forward or backward. + for (int direction=0;direction<2;direction++) // direction=0 --> forward, direction=1 --> backward + { + // Fill the list forward (resp. backward) from the start segment: + int activeSeg = startSeg; + int prevPointId = -20; + int ptId; + while (!edgeSet.empty()) + { + if (!(direction == 1 && prevPointId==-20)) // prevent adding twice startSeg + { + if (direction==0) + linePiece.push_back(activeSeg); + else + linePiece.push_front(activeSeg); + edgeSet.erase(activeSeg); + } + + int ptId1 = d[dI[activeSeg]], ptId2 = d[dI[activeSeg]+1]; + ptId = direction ? (ptId1 == prevPointId ? ptId2 : ptId1) : (ptId2 == prevPointId ? ptId1 : ptId2); + if (dsi[ptId] == 1) // hitting the end of the line + break; + prevPointId = ptId; + int seg1 = rD[rDI[ptId]], seg2 = rD[rDI[ptId]+1]; + activeSeg = (seg1 == activeSeg) ? seg2 : seg1; + } + } + // Done, save final piece into DA: + std::copy(linePiece.begin(), linePiece.end(), result->getPointer()+newIdx); + newIdx += linePiece.size(); + + // identify next valid start segment (one which is not consumed) + if(!edgeSet.empty()) + startSeg = *(edgeSet.begin()); + } + while (!edgeSet.empty()); + _d->decrRef(); _dI->decrRef(); _rD->decrRef(); _rDI->decrRef(); + return result.retn(); +} + void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map& m, int forbVal0, int forbVal1, std::vector& isect) { diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 8a3c98a27..6cc91f357 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -264,6 +264,7 @@ namespace ParaMEDMEM DataArrayInt *& commonCellsArr, DataArrayInt *& commonCellsIArr); MEDCOUPLING_EXPORT DataArrayInt *buildUnionOf2DMesh() const; MEDCOUPLING_EXPORT DataArrayInt *buildUnionOf3DMesh() const; + MEDCOUPLING_EXPORT DataArrayInt *orderConsecutiveCells1D() const; private: MEDCouplingUMesh(); MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCopy); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx index 26b1ef7e7..0d0c9c89c 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx @@ -2133,3 +2133,56 @@ void MEDCouplingBasicsTest5::testIntersect2DMeshWith1DLine1() connI2->decrRef(); } +void MEDCouplingBasicsTest5::testOrderConsecutiveCells1D1() +{ + // A line in several unconnected pieces: + MEDCouplingUMesh *m2 = MEDCouplingUMesh::New("bla", 1); + const int conn2A[30] = {INTERP_KERNEL::NORM_SEG2,0,1,INTERP_KERNEL::NORM_SEG3,1,3,2, INTERP_KERNEL::NORM_SEG2,3,4, + INTERP_KERNEL::NORM_SEG3,5,7,6, INTERP_KERNEL::NORM_SEG3,7,9,8, INTERP_KERNEL::NORM_SEG2,9,10, + INTERP_KERNEL::NORM_SEG2,11,12,INTERP_KERNEL::NORM_SEG2,12,13, + INTERP_KERNEL::NORM_SEG2,14,15}; + const int connI2A[10] = {0,3,7,10,14,18,21,24,27,30}; + double coo2A[32] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,00,0,0,0,0,0,0,0,0,0,0,0}; + for(int i=0; i < 16; i++, coo2A[2*i]=0.0, coo2A[2*i+1]=double(i)); + DataArrayDouble *coord2 = DataArrayDouble::New(); + DataArrayInt *conn2 = DataArrayInt::New(); + DataArrayInt *connI2 = DataArrayInt::New(); + coord2->alloc(16,2); conn2->alloc(30,1); connI2->alloc(10, 1); + std::copy(coo2A ,coo2A+32, coord2->getPointer()); + std::copy(conn2A ,conn2A+30, conn2->getPointer()); + std::copy(connI2A,connI2A+10, connI2->getPointer()); + m2->setCoords(coord2); + m2->setConnectivity(conn2, connI2); + m2->checkCoherency2(1.0e-8); + + // Shuffle a bit :-) + const int shuf[9] = {0,3,6,8,1,4,7,5,2}; + m2->renumberCells(shuf, true); + DataArrayInt * res = m2->orderConsecutiveCells1D(); + const int expRes[9] = {0,3,6,8,1,4,2,7,5}; + CPPUNIT_ASSERT_EQUAL(m2->getNumberOfCells(),res->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expRes, expRes+9,res->getConstPointer())); + + // A closed line (should also work) + MEDCouplingUMesh *m3 = MEDCouplingUMesh::New("bla3", 1); + const int conn3A[10] = {INTERP_KERNEL::NORM_SEG2,0,1,INTERP_KERNEL::NORM_SEG3,1,3,2, INTERP_KERNEL::NORM_SEG2,3,0}; + DataArrayDouble *coord3 = coord2->selectByTupleId2(0, 5,1); + conn2->reAlloc(10); connI2->reAlloc(4); + std::copy(conn3A ,conn3A+10, conn2->getPointer()); + m3->setCoords(coord3); + m3->setConnectivity(conn2, connI2); + m3->checkCoherency2(1.0e-8); + DataArrayInt * res2 = m3->orderConsecutiveCells1D(); + const int expRes2[3] = {0,1,2}; + CPPUNIT_ASSERT_EQUAL(m3->getNumberOfCells(),res2->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expRes2, expRes2+3,res2->getConstPointer())); + + coord2->decrRef(); + coord3->decrRef(); + conn2->decrRef(); + connI2->decrRef(); + m2->decrRef(); + m3->decrRef(); + res->decrRef(); + res2->decrRef(); +} diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx index 87dcf1051..d9a3a8de6 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx @@ -77,6 +77,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( testDAIIsStrictlyMonotonic1 ); CPPUNIT_TEST( testSimplexize3 ); CPPUNIT_TEST( testIntersect2DMeshWith1DLine1 ); + CPPUNIT_TEST( testOrderConsecutiveCells1D1 ); CPPUNIT_TEST_SUITE_END(); public: void testUMeshTessellate2D1(); @@ -120,6 +121,7 @@ namespace ParaMEDMEM void testDAIIsStrictlyMonotonic1(); void testSimplexize3(); void testIntersect2DMeshWith1DLine1(); + void testOrderConsecutiveCells1D1(); }; } diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index b211bce04..b512dd297 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -14630,7 +14630,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(mu.getNodalConnectivityIndex().isEqual(DataArrayInt([0,5,10,15,20,25,30]))) coo0=DataArrayDouble([(0,0,0),(1,0,0),(2,0,0),(0,1,0),(1,1,0),(2,1,0),(0,2,0),(1,2,0),(2,2,0),(0,3,0),(1,3,0),(2,3,0)]) self.assertTrue(mu.getCoords().isEqual(coo0,1e-12)) - mu.writeVTK("tutu.vtu") + #mu.writeVTK("tutu.vtu") # m=MEDCouplingCMesh() arrX=DataArrayDouble(3) ; arrX.iota() @@ -14816,6 +14816,42 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertEqual(expConnI, m3.getNodalConnectivityIndex().getValues()) pass + def testOrderConsecutiveCells1D1(self): + # A line in several unconnected pieces: + m2 = MEDCouplingUMesh.New("bla", 1) + c = DataArrayInt([NORM_SEG2,0,1,NORM_SEG3,1,3,2, NORM_SEG2,3,4, + NORM_SEG3,5,7,6, NORM_SEG3,7,9,8, NORM_SEG2,9,10, + NORM_SEG2,11,12,NORM_SEG2,12,13, + NORM_SEG2,14,15]) + cI = DataArrayInt([0,3,7,10,14,18,21,24,27,30]) + coords2 = DataArrayDouble([float(i) for i in range(32)], 16,2) + m2.setCoords(coords2); + m2.setConnectivity(c, cI); + m2.checkCoherency2(1.0e-8); + + # Shuffle a bit :-) + m2.renumberCells(DataArrayInt([0,3,6,8,1,4,7,5,2]), True); + res = m2.orderConsecutiveCells1D() + expRes = [0,3,6,8,1,4,2,7,5] + self.assertEqual(m2.getNumberOfCells(),res.getNumberOfTuples()) + self.assertEqual(expRes, res.getValues()) + + # A closed line (should also work) + m3 = MEDCouplingUMesh.New("bla3", 1) + conn3A = DataArrayInt([NORM_SEG2,0,1,NORM_SEG3,1,3,2, NORM_SEG2,3,0]) + coord3 = coords2[0:5] + c.reAlloc(10) + cI.reAlloc(4) + + m3.setCoords(coord3) + m3.setConnectivity(conn3A, cI) + m3.checkCoherency2(1.0e-8) + res2 = m3.orderConsecutiveCells1D() + expRes2 = [0,1,2] + self.assertEqual(m3.getNumberOfCells(),res2.getNumberOfTuples()) + self.assertEqual(expRes2, res2.getValues()) + pass + def setUp(self): pass pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 467474d25..b9c814669 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -273,6 +273,7 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTreeFast; %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic; %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic; +%newobject ParaMEDMEM::MEDCouplingUMesh::orderConsecutiveCells1D; %newobject ParaMEDMEM::MEDCouplingUMeshCellByTypeEntry::__iter__; %newobject ParaMEDMEM::MEDCouplingUMeshCellEntry::__iter__; %newobject ParaMEDMEM::MEDCoupling1GTUMesh::New; @@ -1548,6 +1549,7 @@ namespace ParaMEDMEM DataArrayDouble *getBoundingBoxForBBTreeFast() const throw(INTERP_KERNEL::Exception); DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception); DataArrayDouble *getBoundingBoxForBBTree1DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception); + DataArrayInt *orderConsecutiveCells1D() const throw(INTERP_KERNEL::Exception); int split2DCells(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *midOpt=0, const DataArrayInt *midOptI=0) throw(INTERP_KERNEL::Exception); static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da) throw(INTERP_KERNEL::Exception); static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception); -- 2.39.2