]> SALOME platform Git repositories - tools/medcoupling.git/commitdiff
Salome HOME
Fixed colinearize2D to handle degenerated cases (colinear SEG3 forming a circle).
authorabn <adrien.bruneton@cea.fr>
Wed, 29 Apr 2015 08:52:31 +0000 (10:52 +0200)
committerabn <adrien.bruneton@cea.fr>
Wed, 29 Apr 2015 08:54:52 +0000 (10:54 +0200)
Introduced getMiddleOfPointsOriented() in the Edge classes.

src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx
src/INTERP_KERNELTest/QuadraticPlanarInterpTest.hxx
src/INTERP_KERNELTest/QuadraticPlanarInterpTest5.cxx
src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
src/MEDLoader/MEDFileMesh.cxx

index 99edac74b5a635b824dc73a73e8f6a70e87301cb..16625228794d106c2a78d875e0e932208a6d1dee 100644 (file)
@@ -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)
 {
index 3ef275c782c38ae13b8534769511bd0751b77f7c..006fb25656a1a2fe95b1bf10657a5148ea90896c 100644 (file)
@@ -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.
index 92f43b6a949d27acb1415b4762b0387c20bafb2b..66c6d2ea4732532d69b45fa2cc738c314e96ddc0 100644 (file)
@@ -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.
  */
index 3b46d7935f4385eef165a3d587fb00f00a1fdd2c..382077bc4348b68f35ec84b7af53bc4579c75c21 100644 (file)
@@ -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;
index 63c457c86f61df156e98897e64a606613b8a3f06..991e47b281b8d0864b3d236d2ffe7e4d1a369351 100644 (file)
@@ -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);
index 6d4862337fa2dbf66f73d9848ee47feeb5f44fef..92e9a959128643536bd19ba1631ebd29e70fec9b 100644 (file)
@@ -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();
+  }
+}
 
 }
index f213a9d42d53de8b084d6e3fb19c1913bfbc112e..9e5c5b950c3d03e69373c427908d532a2cca8c33 100644 (file)
@@ -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<int>& 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<int>& 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<int> 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<int> middles;
   bool ret(false);
-  for(;nbOfHit<nbs;nbOfTurn++)
+  for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
     {
       cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
       std::map<MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int> m;
       INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
-      posEndElt++;
-      nbOfHit++;
-      unsigned endI(nbs-nbOfHit);
-      for(unsigned i=0;i<endI;i++)
+      posEndElt = posBaseElt+1;
+
+      // Look backward first: are the final edges of the cells colinear with the first ones?
+      // This initializes posBaseElt.
+      if(nbOfTurn==0)
+        {
+          for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell wih one single edge
+            {
+              cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn,typeOfSon);
+              INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
+              INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
+              bool isColinear=eint->areColinears();
+              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;j<endJ && nbOfHit<maxNbOfHit;j++)  // 2nd condition is to avoid ending with a cell wih one single edge
         {
-          cm.fillSonCellNodalConnectivity2(posBaseElt+(int)i+1,connBg+1,sz,tmpConn,typeOfSon);
+          cm.fillSonCellNodalConnectivity2(fwdStart+(int)j,connBg+1,sz,tmpConn,typeOfSon); // get edge #(posBaseElt+i)'s connectivity
           INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
           INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
           bool isColinear(eint->areColinears());
@@ -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;ii<endII;ii++)
-                    {
-                      cm.fillSonCellNodalConnectivity2(nbs-ii-1,connBg+1,sz,tmpConn,typeOfSon);
-                      eCand=MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m);
-                      eint=INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand);
-                      isColinear=eint->areColinears();
-                      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();
index c6bc0d7d085687dc7978453e9ae27eca387c7171..b0cbd8274012590787897a1ceae886211c867930 100644 (file)
@@ -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])))
index f90a4603a677a8f886a7dce149bd23a667e8c27c..454d90174806314c12c0ec5398695d272ea0c830 100644 (file)
@@ -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);