From: abn Date: Wed, 29 Apr 2015 08:52:31 +0000 (+0200) Subject: 1) Fixed colinearize2D to handle degenerated cases (colinear SEG3 forming a circle). X-Git-Tag: V7_6_0rc1~4 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=afe1d4b0731be1e93acd208b1c83ada45fe50753;p=tools%2Fmedcoupling.git 1) Fixed colinearize2D to handle degenerated cases (colinear SEG3 forming a circle). Introduced getMiddleOfPointsOriented() in the Edge classes. 2) GetAbsoluteAngleOfNormalizedVect() actually does exactly the same as atan2() Micro correction of ref test. Valgrinification of MEDCouplingBasicsTest.py remaining. --- 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..491ff1578 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -529,25 +529,7 @@ double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect) */ double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy) { - //When arc is lower than 0.707 Using Asin - if(fabs(ux)<0.707) - { - double ret=SafeAcos(ux); - if(uy>0.) - return ret; - ret=-ret; - return ret; - } - else - { - double ret=SafeAsin(uy); - if(ux>0.) - return ret; - if(ret>0.) - return M_PI-ret; - else - return -M_PI-ret; - } + return atan2(uy, ux); } void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end, @@ -649,6 +631,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); @@ -656,14 +644,42 @@ void EdgeArcCircle::getMiddleOfPoints(const double *p1, const double *p2, double // double myDelta1(angle1-_angle0),myDelta2(angle2-_angle0); if(_angle>0.) - { myDelta1=myDelta1>=0.?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>=0.?myDelta2:myDelta2+2.*M_PI; } + { myDelta1=myDelta1>-QUADRATIC_PLANAR::_precision?myDelta1:myDelta1+2.*M_PI; myDelta2=myDelta2>-QUADRATIC_PLANAR::_precision?myDelta2:myDelta2+2.*M_PI; } else - { myDelta1=myDelta1<=0.?myDelta1:myDelta1-2.*M_PI; myDelta2=myDelta2<=0.?myDelta2:myDelta2-2.*M_PI; } + { myDelta1=myDelta1_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..615969ffb 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(0.37969180470645592, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.4372640310451197, mid[1], 1.e-7); + + e.getMiddleOfPoints(p2, p1, mid); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.37969180470645592, mid[0], 1.e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.4372640310451197, 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 b8421f251..ab7425894 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -242,6 +242,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; @@ -249,20 +257,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()); } } @@ -11358,6 +11366,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) @@ -11371,7 +11394,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 @@ -11388,7 +11411,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 @@ -11398,13 +11421,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 @@ -11426,22 +11450,47 @@ 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 + for(unsigned j=fwdStart+1;jareColinears()); @@ -11454,37 +11503,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 e8681a913..1675b5a8e 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -14747,10 +14747,59 @@ class MEDCouplingBasicsTest(unittest.TestCase): m.colinearize2D(1e-12) m.checkCoherency2() self.assertEqual(refPtr,m.getCoords().getHiddenCppPointer()) - self.assertTrue(m.getNodalConnectivity().isEqual(DataArrayInt([5,0,2,3,4]))) + self.assertTrue(m.getNodalConnectivity().isEqual(DataArrayInt([NORM_POLYGON,0,2,3,4]))) 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 6907c90e1..b3d9a9ee2 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -3294,7 +3294,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]; @@ -3308,7 +3308,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]; @@ -4331,7 +4331,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); @@ -4367,7 +4367,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); @@ -4401,7 +4401,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);