From b3378de3546e624de2f4454a7d5adc409bddef9d Mon Sep 17 00:00:00 2001 From: abn Date: Wed, 29 Apr 2015 10:52:31 +0200 Subject: [PATCH] Fixed colinearize2D to handle degenerated cases (colinear SEG3 forming a circle). Introduced getMiddleOfPointsOriented() in the Edge classes. --- .../Geometric2D/InterpKernelGeo2DEdge.cxx | 5 + .../Geometric2D/InterpKernelGeo2DEdge.hxx | 3 + .../InterpKernelGeo2DEdgeArcCircle.cxx | 34 ++++++ .../InterpKernelGeo2DEdgeArcCircle.hxx | 1 + .../QuadraticPlanarInterpTest.hxx | 4 + .../QuadraticPlanarInterpTest5.cxx | 92 +++++++++++++++ src/MEDCoupling/MEDCouplingUMesh.cxx | 107 +++++++++++------- src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 55 +++++++++ src/MEDLoader/MEDFileMesh.cxx | 10 +- 9 files changed, 268 insertions(+), 43 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx index 99edac74b..166252287 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx @@ -726,6 +726,11 @@ void Edge::unApplySimilarity(double xBary, double yBary, double dimChar) _bounds.unApplySimilarity(xBary,yBary,dimChar); } +void Edge::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const +{ + return getMiddleOfPoints(p1, p2, mid); +} + bool Edge::Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, const Bounds *whereToFind, MergePoints& commonNode, ComposedEdge& outValForF1, ComposedEdge& outValForF2) { diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx index 3ef275c78..006fb2565 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx @@ -248,7 +248,10 @@ namespace INTERP_KERNEL virtual double getCurveLength() const = 0; virtual void getBarycenter(double *bary) const = 0; virtual void getBarycenterOfZone(double *bary) const = 0; + //! return the middle of two points virtual void getMiddleOfPoints(const double *p1, const double *p2, double *mid) const = 0; + //! return the middle of two points respecting the orientation defined by this (relevant for arc of circle). By default same as getMiddleOfPoints() + virtual void getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const; //! Retrieves a point that is owning to this, well placed for IN/OUT detection of this. Typically midlle of this is returned. virtual Node *buildRepresentantOfMySelf() const = 0; //! Given a magnitude specified by sub-type returns if in or not. See getCharactValue method. diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index 92f43b6a9..66c6d2ea4 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -649,6 +649,12 @@ void EdgeArcCircle::getBarycenterOfZone(double *bary) const +tmp2*(tmp4-tmp3+(tmp3*tmp3*tmp3-tmp4*tmp4*tmp4)/3.)/2.; } +/** + * Compute the "middle" of two points on the arc of circle. + * The order (p1,p2) or (p2,p1) doesn't matter. p1 and p2 have to be localized on the edge defined by this. + * \param[out] mid the point located half-way between p1 and p2 on the arc defined by this. + * \sa getMiddleOfPointsOriented() a generalisation working also when p1 and p2 are not on the arc. + */ void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const { double dx1((p1[0]-_center[0])/_radius),dy1((p1[1]-_center[1])/_radius),dx2((p2[0]-_center[0])/_radius),dy2((p2[1]-_center[1])/_radius); @@ -664,6 +670,34 @@ void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double mid[1]=_center[1]+_radius*sin(_angle0+(myDelta1+myDelta2)/2.); } +/** + * Compute the "middle" of two points on the arc of circle. + * Walk on the circle from p1 to p2 using the rotation direction indicated by this->_angle (i.e. by the orientation of the arc). + * This function is sensitive to the ordering of p1 and p2. + * \param[out] mid the point located half-way between p1 and p2 + * \sa getMiddleOfPoints() to be used when the order of p1 and p2 is not relevant. + */ +void EdgeArcCircle::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const +{ + double dx1((p1[0]-_center[0])/_radius),dy1((p1[1]-_center[1])/_radius),dx2((p2[0]-_center[0])/_radius),dy2((p2[1]-_center[1])/_radius); + double angle1(GetAbsoluteAngleOfNormalizedVect(dx1,dy1)),angle2(GetAbsoluteAngleOfNormalizedVect(dx2,dy2)); + + if (angle1 <= 0.0) + angle1 += 2.*M_PI; + if (angle2 <= 0.0) + angle2 += 2.*M_PI; + + double avg; + if((_angle>0. && angle1 <= angle2) || (_angle<=0. && angle1 >= angle2)) + avg = (angle1+angle2)/2.; + else + avg = (angle1+angle2)/2. - M_PI; + + mid[0]=_center[0]+_radius*cos(avg); + mid[1]=_center[1]+_radius*sin(avg); +} + + /*! * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x. */ diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx index 3b46d7935..382077bc4 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx @@ -83,6 +83,7 @@ namespace INTERP_KERNEL void getBarycenter(double *bary) const; void getBarycenterOfZone(double *bary) const; void getMiddleOfPoints(const double *p1, const double *p2, double *mid) const; + void getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const; bool isIn(double characterVal) const; Node *buildRepresentantOfMySelf() const; bool isLower(double val1, double val2) const; diff --git a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx index 63c457c86..991e47b28 100644 --- a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx +++ b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx @@ -108,6 +108,8 @@ namespace INTERP_TEST CPPUNIT_TEST( checkMakePartitionAbs1 ); // CPPUNIT_TEST( checkIsInOrOut ); + CPPUNIT_TEST( checkGetMiddleOfPoints ); + CPPUNIT_TEST( checkGetMiddleOfPointsOriented ); CPPUNIT_TEST_SUITE_END(); public: void setUp(); @@ -198,6 +200,8 @@ namespace INTERP_TEST void checkMakePartitionAbs1(); // From Adrien: void checkIsInOrOut(); + void checkGetMiddleOfPoints(); + void checkGetMiddleOfPointsOriented(); private: INTERP_KERNEL::QuadraticPolygon *buildQuadraticPolygonCoarseInfo(const double *coords, const int *conn, int lgth); diff --git a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest5.cxx b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest5.cxx index 6d4862337..92e9a9591 100644 --- a/src/INTERP_KERNELTest/QuadraticPlanarInterpTest5.cxx +++ b/src/INTERP_KERNELTest/QuadraticPlanarInterpTest5.cxx @@ -1190,5 +1190,97 @@ void QuadraticPlanarInterpTest::checkIsInOrOut() delete pol1; } +void QuadraticPlanarInterpTest::checkGetMiddleOfPoints() +{ + { // from testIntersect2DMeshWith1DLine6() + double p1[] = {0.51641754716735844, 2.0}; + double p2[] = {0.0, 1.0}; + double e_center[] = {-0.71, 2.0}; + double mid[] = {0.0,0.0}; // out + double mide[] = {0.0,0.0}; // expected + + Node * start = new Node(0.,0.); Node * end = new Node(0.,0.); // unused + // start, end, center_x, center_y, radius, angle0, angle + EdgeArcCircle e(start, end, e_center, 1.2264175471673588, -0.9533904350433241, 0.95339043504332388); + + e.getMiddleOfPoints(p1, p2, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.7996918047064556, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.5627359689548808, mid[1], 1.e-7); + + e.getMiddleOfPoints(p2, p1, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.7996918047064556, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.5627359689548808, mid[1], 1.e-7); + + start->decrRef(); end->decrRef(); + } + { // from testSwig2Intersect2DMeshWith1DLine11() + double p1[] = {-1., 0.23453685964236054}; + double p2[] = {-0.23453685964235979, 1.0}; + double e_center[] = {-4.85, 4.85}; + double mid[] = {0.0,0.0}; // out + + Node * start = new Node(0.,0.); Node * end = new Node(0.,0.); // unused + // start, end, center_x, center_y, radius, angle0, angle + EdgeArcCircle e(start, end, e_center, 6.0104076400856474, -0.69522150912422953, -0.18035330854643861); + + e.getMiddleOfPoints(p1, p2, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.6, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.6, mid[1], 1.e-7); + + e.getMiddleOfPoints(p2, p1, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.6, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.6, mid[1], 1.e-7); + + start->decrRef(); end->decrRef(); + } + { // from testSwig2Intersect2DMeshWith1DLine11() + double p1[] = {-0.1303327636866019, -1.0}; + double p2[] = {-1.0, -0.1303327636866019}; + double e_center[] = {-1.9833333333333298, -1.9833333333333298}; + double mid[] = {0.0,0.0}; // out + + Node * start = new Node(0.,0.); Node * end = new Node(0.,0.); // unused + // start, end, center_x, center_y, radius, angle0, angle + EdgeArcCircle e(start, end, e_center, 2.0977501175200861, 1.0829141821052615, -0.59503203741562627); + + e.getMiddleOfPoints(p1, p2, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.5, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.5, mid[1], 1.e-7); + + e.getMiddleOfPoints(p2, p1, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.5, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.5, mid[1], 1.e-7); + + start->decrRef(); end->decrRef(); + } +} + +void QuadraticPlanarInterpTest::checkGetMiddleOfPointsOriented() +{ + { // from testSwig2Colinearize2D3() + double p1[] = {-0.70710678118654746, 0.70710678118654757}; + double p2[] = {-0.70710678118654768, -0.70710678118654746}; + double e_center[] = {0., 0.}; + double mid[] = {0.0,0.0}; // out + + Node * start = new Node(0.,0.); Node * end = new Node(0.,0.); // unused + // start, end, center_x, center_y, radius, angle0, angle + EdgeArcCircle e(start, end, e_center, 1.0, -0.7853981633974485, -1.5707963267948966); + + e.getMiddleOfPointsOriented(p1, p2, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1., mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0., mid[1], 1.e-7); + + e.getMiddleOfPoints(p1, p2, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0., mid[1], 1.e-7); + + e.getMiddleOfPointsOriented(p2, p1, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0., mid[1], 1.e-7); + + start->decrRef(); end->decrRef(); + } +} } diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index f213a9d42..9e5c5b950 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -241,6 +241,14 @@ void MEDCouplingUMesh::checkCoherency1(double eps) const oss << " nodes whereas in connectivity there is " << nbOfNodesInCell << " nodes ! Looks very bad !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } + if(cm.isQuadratic() && cm.isDynamic() && meshDim == 2) + if (nbOfNodesInCell % 2 || nbOfNodesInCell < 4) + { + std::ostringstream oss; + oss << "MEDCouplingUMesh::checkCoherency1 : cell #" << i << " with quadratic type '" << cm.getRepr() << "' has " << nbOfNodesInCell; + oss << " nodes. This should be even, and greater or equal than 4!! Looks very bad!"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } for(const int *w=ptr+ptrI[i]+1;w!=ptr+ptrI[i+1];w++) { int nodeId=*w; @@ -248,20 +256,20 @@ void MEDCouplingUMesh::checkCoherency1(double eps) const { if(nodeId>=nbOfNodes) { - std::ostringstream oss; oss << "Cell #" << i << " is consituted of node #" << nodeId << " whereas there are only " << nbOfNodes << " nodes !"; + std::ostringstream oss; oss << "Cell #" << i << " is built with node #" << nodeId << " whereas there are only " << nbOfNodes << " nodes in the mesh !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } } else if(nodeId<-1) { - std::ostringstream oss; oss << "Cell #" << i << " is consituted of node #" << nodeId << " in connectivity ! sounds bad !"; + std::ostringstream oss; oss << "Cell #" << i << " is built with node #" << nodeId << " in connectivity ! sounds bad !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } else { if((INTERP_KERNEL::NormalizedCellType)(ptr[ptrI[i]])!=INTERP_KERNEL::NORM_POLYHED) { - std::ostringstream oss; oss << "Cell #" << i << " is consituted of node #-1 in connectivity ! sounds bad !"; + std::ostringstream oss; oss << "Cell #" << i << " is built with node #-1 in connectivity ! sounds bad !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } } @@ -11338,6 +11346,21 @@ int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, in } } +int InternalAddPointOriented(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter) +{ + if(id!=-1) + return id; + else + { + int ret(nodesCnter++); + double newPt[2]; + e->getMiddleOfPointsOriented(coo+2*startId,coo+2*endId,newPt); + addCoo.insertAtTheEnd(newPt,newPt+2); + return ret; + } +} + + /// @cond INTERNAL void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) @@ -11351,7 +11374,7 @@ void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int st if(stp-start>1) { int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); - InternalAddPoint(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2); + InternalAddPointOriented(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2); middles.push_back(tmp3+offset); } else @@ -11368,7 +11391,7 @@ void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int s if(stp-start>1) { int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); - InternalAddPoint(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2); + InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2); middles.push_back(tmp3+offset); } else @@ -11378,13 +11401,14 @@ void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int s void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) { + // only the quadratic point to deal with: if(linOrArc) { if(stp-start>1) { int tmpSrt(connBg[start]),tmpEnd(connBg[stp]); int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); - InternalAddPoint(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2); + InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2); middles.push_back(tmp3+offset); } else @@ -11406,22 +11430,48 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg sz--; INTERP_KERNEL::AutoPtr tmpConn(new int[sz]); const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0])); - unsigned nbs(cm.getNumberOfSons2(connBg+1,sz)),nbOfHit(0); + unsigned nbs(cm.getNumberOfSons2(connBg+1,sz)); + unsigned nbOfHit(0); // number of fusions operated int posBaseElt(0),posEndElt(0),nbOfTurn(0); + const unsigned int maxNbOfHit = cm.isQuadratic() ? nbs-2 : nbs-3; // a quad cell is authorized to end up with only two edges, a linear one has to keep 3 at least INTERP_KERNEL::NormalizedCellType typeOfSon; std::vector middles; bool ret(false); - for(;nbOfHit,int> m; INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m)); - posEndElt++; - nbOfHit++; - unsigned endI(nbs-nbOfHit); - for(unsigned i=0;iareColinears(); + if(isColinear) + { + nbOfHit++; + posBaseElt--; + ret=true; + } + delete eint; + eCand->decrRef(); + if(!isColinear) + break; + } + } + // Now move forward: + const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward + const unsigned endJ(nbs-nbOfHit); // put in a constant because modified in the loop + for(unsigned j=1;jareColinears()); @@ -11434,37 +11484,18 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg delete eint; eCand->decrRef(); if(!isColinear) - { - if(nbOfTurn==0) - {//look if the first edge of cell is not colinear with last edges in this case the start of nodal connectivity is shifted back - unsigned endII(nbs-nbOfHit-1);//warning nbOfHit can be modified, so put end condition in a variable. - for(unsigned ii=0;iiareColinears(); - if(isColinear) - { - nbOfHit++; - posBaseElt--; - ret=true; - } - delete eint; - eCand->decrRef(); - if(!isColinear) - break; - } - } break; - } } //push [posBaseElt,posEndElt) in newConnOfCell using e + // The if clauses below are (volontary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its begining! if(nbOfTurn==0) + // at the begining of the connectivity (insert type) EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); - else if(nbOfHit!=nbs) + else if((nbOfHit+nbOfTurn) != (nbs-1)) + // in the middle EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); - else + if ((nbOfHit+nbOfTurn) == (nbs-1)) + // at the end (only quad points to deal with) EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); posBaseElt=posEndElt; e->decrRef(); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index c6bc0d7d0..b0cbd8274 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -14622,8 +14622,14 @@ class MEDCouplingBasicsTest(unittest.TestCase): coo2=DataArrayDouble([(-5,0),(-1,0),(7.,6.),(7,0),(1,6),(1,0),(-3,0),(8.2426406871192839,3),(5,0),(3,0), (2,0),(4,0),(6,0),(7.9196888946291288,1.3764116995614091),(7.9196888946291288,4.6235883004385911),(4,7.2426406871192848),(-2,3),(-4,0),(-2,0),(0,0)]) m=MEDCouplingUMesh("mesh",2) ; m.setCoords(coo2) ; m.allocateCells() m.insertNextCell(NORM_QPOLYG,[5,9,8,3,7,2,4,0,6,1,10,11,12,13,14,15,16,17,18,19]) + mm = m.deepCpy() + mm.tessellate2D(0.1) + mm.writeVTK("/tmp/colin_1_a.vtu") refPtr=m.getCoords().getHiddenCppPointer() self.assertTrue(m.colinearize2D(1e-12).isEqual(DataArrayInt([0]))) + mm = m.deepCpy() + mm.tessellate2D(0.1) + mm.writeVTK("/tmp/colin_1_b.vtu") self.assertNotEqual(refPtr,m.getCoords().getHiddenCppPointer())#not same coordinates here self.assertTrue(coo2.isEqual(m.getCoords()[:20],1e-12)) self.assertTrue(m.getCoords()[20:].isEqual(DataArrayDouble([(1.,0.),(7.,6.)]),1e-12)) @@ -14751,6 +14757,55 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(m.getNodalConnectivityIndex().isEqual(DataArrayInt([0,5]))) pass + def testSwig2Colinearize2D3(self): + """ colinearize was too agressive, potentially producing cells with one edge """ + # Flat polygon with 3 edges - nothing should happen (min number of edges for a linear polyg) + coo = DataArrayDouble([0.0,0.0, 2.0,0.0, 1.5,0.0, 1.0,0.0, 0.5,0.0], 5,2) + m = MEDCouplingUMesh("m", 2) + c, cI = [DataArrayInt(l) for l in [[NORM_POLYGON, 0,1,2], [0,4]] ] + m.setCoords(coo); m.setConnectivity(c, cI) + m.colinearize2D(1e-10) + m.checkCoherency2() + self.assertEqual(c.getValues(), m.getNodalConnectivity().getValues()) + self.assertEqual(cI.getValues(), m.getNodalConnectivityIndex().getValues()) + + # Flat quad polygon, 2 edges - nothing should happen (min number of edges for a quad polyg) + m = MEDCouplingUMesh("m", 2) + c, cI = [DataArrayInt(l) for l in [[NORM_QPOLYG, 0,1, 2,3], [0,5]] ] + m.setCoords(coo); m.setConnectivity(c, cI) + m.colinearize2D(1e-10) + m.checkCoherency2() + self.assertEqual(c.getValues(), m.getNodalConnectivity().getValues()) + self.assertEqual(cI.getValues(), m.getNodalConnectivityIndex().getValues()) + + # Flat polygon, 4 edges - one reduction should happen + m = MEDCouplingUMesh("m", 2) + c, cI = [DataArrayInt(l) for l in [[NORM_POLYGON, 0,1,2,3], [0,5]] ] + m.setCoords(coo); m.setConnectivity(c, cI) + m.colinearize2D(1e-10) + m.checkCoherency2() + self.assertEqual([NORM_POLYGON, 3,1,2], m.getNodalConnectivity().getValues()) + self.assertEqual([0,4], m.getNodalConnectivityIndex().getValues()) + + # Flat quad polygon, 3 edges - one reduction expected + m = MEDCouplingUMesh("m", 2) + c, cI = [DataArrayInt(l) for l in [[NORM_QPOLYG, 0,1,3, 3,2,4], [0,7]] ] + m.setCoords(coo); m.setConnectivity(c, cI) + m.colinearize2D(1e-10) + m.checkCoherency2() + self.assertEqual([NORM_QPOLYG, 3,1, 5,2], m.getNodalConnectivity().getValues()) + self.assertTrue( m.getCoords()[5].isEqual( DataArrayDouble([(1.5,0.0)]), 1.0e-12 ) ) + self.assertEqual([0,5], m.getNodalConnectivityIndex().getValues()) + + # Now an actual (neutronic) case: circle made of 4 SEG3. Should be reduced to 2 SEG3 + m = MEDCouplingDataForTest.buildCircle2(0.0, 0.0, 1.0) + c, cI = [DataArrayInt(l) for l in [[NORM_QPOLYG, 7,5,3,1, 6,4,2,0], [0,9]] ] + m.colinearize2D(1e-10) + m.checkCoherency2() + self.assertEqual([NORM_QPOLYG, 3,5, 8,4], m.getNodalConnectivity().getValues()) + self.assertTrue( m.getCoords()[8].isEqual( DataArrayDouble([(1.0,0.0)]), 1.0e-12 ) ) + self.assertEqual([0,5], m.getNodalConnectivityIndex().getValues()) + def testSwig2CheckAndPreparePermutation2(self): a=DataArrayInt([10003,9999999,5,67]) self.assertTrue(DataArrayInt.CheckAndPreparePermutation(a).isEqual(DataArrayInt([2,3,0,1]))) diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index f90a4603a..454d90174 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -3219,7 +3219,7 @@ const MEDFileUMeshSplitL1 *MEDFileUMesh::getMeshAtLevSafe(int meshDimRelToMaxExt throw INTERP_KERNEL::Exception("Dimension request is invalid (>1) !"); int tracucedRk=-meshDimRelToMaxExt; if(tracucedRk>=(int)_ms.size()) - throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! To low !"); + throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! Too low !"); if((const MEDFileUMeshSplitL1 *)_ms[tracucedRk]==0) throw INTERP_KERNEL::Exception("On specified lev (or entity) no cells exists !"); return _ms[tracucedRk]; @@ -3233,7 +3233,7 @@ MEDFileUMeshSplitL1 *MEDFileUMesh::getMeshAtLevSafe(int meshDimRelToMaxExt) throw INTERP_KERNEL::Exception("Dimension request is invalid (>1) !"); int tracucedRk=-meshDimRelToMaxExt; if(tracucedRk>=(int)_ms.size()) - throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! To low !"); + throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! Too low !"); if((const MEDFileUMeshSplitL1 *)_ms[tracucedRk]==0) throw INTERP_KERNEL::Exception("On specified lev (or entity) no cells exists !"); return _ms[tracucedRk]; @@ -4247,7 +4247,7 @@ void MEDFileUMesh::setFamilyFieldArr(int meshDimRelToMaxExt, DataArrayInt *famAr throw INTERP_KERNEL::Exception("MEDFileUMesh::setFamilyFieldArr : Dimension request is invalid (>1) !"); int traducedRk=-meshDimRelToMaxExt; if(traducedRk>=(int)_ms.size()) - throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! To low !"); + throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! Too low !"); if((MEDFileUMeshSplitL1 *)_ms[traducedRk]==0) throw INTERP_KERNEL::Exception("On specified lev (or entity) no cells exists !"); return _ms[traducedRk]->setFamilyArr(famArr); @@ -4283,7 +4283,7 @@ void MEDFileUMesh::setRenumFieldArr(int meshDimRelToMaxExt, DataArrayInt *renumA throw INTERP_KERNEL::Exception("MEDFileUMesh::setRenumArr : Dimension request is invalid (>1) !"); int traducedRk=-meshDimRelToMaxExt; if(traducedRk>=(int)_ms.size()) - throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! To low !"); + throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! Too low !"); if((MEDFileUMeshSplitL1 *)_ms[traducedRk]==0) throw INTERP_KERNEL::Exception("On specified lev (or entity) no cells exists !"); return _ms[traducedRk]->setRenumArr(renumArr); @@ -4317,7 +4317,7 @@ void MEDFileUMesh::setNameFieldAtLevel(int meshDimRelToMaxExt, DataArrayAsciiCha throw INTERP_KERNEL::Exception("MEDFileUMesh::setNameFieldAtLevel : Dimension request is invalid (>1) !"); int traducedRk=-meshDimRelToMaxExt; if(traducedRk>=(int)_ms.size()) - throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! To low !"); + throw INTERP_KERNEL::Exception("Invalid mesh dim relative to max given ! Too low !"); if((MEDFileUMeshSplitL1 *)_ms[traducedRk]==0) throw INTERP_KERNEL::Exception("On specified lev (or entity) no cells exists !"); return _ms[traducedRk]->setNameArr(nameArr); -- 2.39.2