]> SALOME platform Git repositories - tools/medcoupling.git/commitdiff
Salome HOME
Intersec bug fix: correct polygons with flat corners, they crash residual computation ...
authorabn <adrien.bruneton@cea.fr>
Wed, 9 Jan 2019 14:30:10 +0000 (15:30 +0100)
committerabn <adrien.bruneton@cea.fr>
Thu, 10 Jan 2019 12:52:36 +0000 (13:52 +0100)
also added test.

src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx
src/MEDCoupling/MEDCouplingUMesh_intersection.cxx
src/MEDCoupling_Swig/MEDCouplingIntersectTest.py

index 3d0f0a468aa57c1b9c712a112ddec2f83a568897..82c8f7b9caa3cfad566dc667dbaa2e0a016f2117 100644 (file)
@@ -97,6 +97,18 @@ void ComposedEdge::pushBack(ComposedEdge *elem)
   _sub_edges.insert(_sub_edges.end(),elemsOfElem->begin(),elemsOfElem->end());
 }
 
+/*!
+ * Warning! This is highly inefficient ...
+ */
+void ComposedEdge::erase(int index)
+{
+  // Not the most efficient thing to do, but rarely called ...
+  auto it = _sub_edges.begin();
+  for (int i=0; i < index; i++, it++);
+  delete (*it);
+  _sub_edges.erase(it++, it);
+}
+
 ElementaryEdge *ComposedEdge::operator[](int i) const
 {
   std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();
index 4d04aa7b9f7ed787200700297cd6ffd45ca0e362..c484b0020c972d4d038eaf0ddfcc5b3b09f11c85 100644 (file)
@@ -98,6 +98,7 @@ namespace INTERP_KERNEL
     INTERPKERNEL_EXPORT void pushBack(Edge *edge, bool direction=true);
     INTERPKERNEL_EXPORT void pushBack(ElementaryEdge *elem);
     INTERPKERNEL_EXPORT void pushBack(ComposedEdge *elem);
+    INTERPKERNEL_EXPORT void erase(int index);
     INTERPKERNEL_EXPORT int size() const { return (int)_sub_edges.size(); }
     INTERPKERNEL_EXPORT ElementaryEdge *operator[](int i) const;
     INTERPKERNEL_EXPORT Node *getEndNode() const;
index 0f62a67c0f6a5819cc422f42322b914870ec8466..3f99fcde688dba6c1af7749557fc4e8b240c3eca 100644 (file)
@@ -72,6 +72,14 @@ void ElementaryEdge::getAllNodes(std::set<Node *>& output) const
   output.insert(_ptr->getEndNode());
 }
 
+bool ElementaryEdge::hasSameExtremities(const ElementaryEdge& other) const
+{
+  std::set<Node *> s1, s2;
+  getAllNodes(s1);
+  other.getAllNodes(s2);
+  return (s1 == s2);
+}
+
 void ElementaryEdge::getBarycenter(double *bary, double& weigh) const
 {
   _ptr->getBarycenter(bary);
index 18f6efafd6f8b242c0e6d7acc8a1e9f234d67945..f585b9aa1b5758df82e9321142082dc4b6abe085 100644 (file)
@@ -50,6 +50,7 @@ namespace INTERP_KERNEL
     INTERPKERNEL_EXPORT void applySimilarity(double xBary, double yBary, double dimChar) { _ptr->applySimilarity(xBary,yBary,dimChar); }
     INTERPKERNEL_EXPORT void unApplySimilarity(double xBary, double yBary, double dimChar) { _ptr->unApplySimilarity(xBary,yBary,dimChar); }
     INTERPKERNEL_EXPORT void getAllNodes(std::set<Node *>& output) const;
+    INTERPKERNEL_EXPORT bool hasSameExtremities(const ElementaryEdge& other) const;
     INTERPKERNEL_EXPORT void getBarycenter(double *bary, double& weigh) const;
     INTERPKERNEL_EXPORT ElementaryEdge *clone() const;
     INTERPKERNEL_EXPORT void initLocations() const;
index 26c7a72231b15a7ca231edaa76fad324fed389da..b4353458dc5227cdf8602cdff8386a886ff35027 100644 (file)
@@ -463,7 +463,7 @@ void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size
  * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the
  * orientation of edge.
  */
-void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
+void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges2,
                                                 const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
                                                 const std::vector< std::vector<int> >& colinear1,
                                                 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
@@ -518,7 +518,7 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL
       if(directos)
         {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
           std::size_t oldSz=_sub_edges.size();
-          appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
+          appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges2);
           std::size_t newSz=_sub_edges.size();
           std::size_t zeSz=newSz-oldSz;
           alreadyExistingIn2[descBg[i]].resize(zeSz);
@@ -528,7 +528,7 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL
         }
       else
         {//there is subpart of edge 'edgeId' of pol2 inside pol1
-          const std::vector<int>& subEdge=intersectEdges[edgeId];
+          const std::vector<int>& subEdge=intersectEdges2[edgeId];
           std::size_t nbOfSubEdges=subEdge.size()/2;
           for(std::size_t j=0;j<nbOfSubEdges;j++)
             {
@@ -703,6 +703,33 @@ void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTE
   unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
 }
 
+/*!
+ * Remove the two 'identical' edges from the list, and drecRef the objects.
+ */
+void QuadraticPolygon::cleanDegeneratedConsecutiveEdges()
+{
+  IteratorOnComposedEdge it2ii(this);
+  ElementaryEdge * prevEdge = 0;
+  int kk = 0;
+  if  (recursiveSize() > 2)
+    for(it2ii.first();!it2ii.finished();it2ii.next())
+      {
+        ElementaryEdge * cur = it2ii.current();
+        if (prevEdge && prevEdge->hasSameExtremities(*cur))
+          {
+            // Delete the two 'identical' edges:
+            it2ii.previousLoop(); it2ii.previousLoop();
+            erase(--kk); erase(kk);
+            prevEdge = it2ii.current();
+          }
+        else
+          {
+            kk++;
+            prevEdge = cur;
+          }
+      }
+}
+
 /*!
  * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method.
  * 'other' is a QuadraticPolygon of \b non closed edges.
index a5afe97e56cbe682cd34104d1fcbe87eb92dca5c..ee286fd397178cb31896205a360cb6f4e45e9b4b 100644 (file)
@@ -93,6 +93,7 @@ namespace INTERP_KERNEL
     INTERPKERNEL_EXPORT bool haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted, bool& direction) const;
     INTERPKERNEL_EXPORT static void ComputeResidual(const QuadraticPolygon& pol1, const std::set<Edge *>& notUsedInPol1, const std::set<Edge *>& edgesInPol2OnBoundary, const std::map<INTERP_KERNEL::Node *,int>& mapp, int offset, int idThis,
                                                     std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nb1, std::vector<int>& nb2);
+    INTERPKERNEL_EXPORT void cleanDegeneratedConsecutiveEdges();
   protected:
     std::list<QuadraticPolygon *> zipConsecutiveInSegments() const;
     void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const;
index 39e7cb4eeae5b548376c62c1ad156d7387ee17fc..f61742e0b22ce0a17d5e2bfb6ce208dbbd63248d 100644 (file)
@@ -1360,6 +1360,11 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo
           pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
               pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
         }
+      // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges:
+      for (auto &p: pol2s)
+        p.cleanDegeneratedConsecutiveEdges();
+      edgesIn2ForShare.clear();  // removing temptation to use it further since it might now contain invalid edges.
+      ///
       ii=0;
       // Now rebuild intersected cells from all this:
       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
index 7694d7870bcec75c67953d22a3cafddeb97e19e4..f19780274ad4db7197f605a39c35d84ab2bc07ae 100644 (file)
@@ -548,6 +548,37 @@ class MEDCouplingIntersectTest(unittest.TestCase):
         self.assertEqual(map2.getValues(), [0])
         pass
 
+    def testIntersect2DMeshesTmp9(self):
+        """ Tricky case: two triangular shapes intersecting, but not perfectly, at their tips. Several issues fixed:
+            - Bug fix: seg seg intersector epsilon is to be taken absolutely for colinearity test (even for very small vectors 
+            we don't want to have false positive on colinearity. So go back to a comparison with an angle.)
+            - when intersecting nodes are merged, they were not properly added on pol2.
+            - bug fix in compute residual: the stop condition is really on pol1Zip only.
+            - correcting polygons with flat corners, they were crashing residual computation
+        """
+        eps = 1.0e-6  # This is the key parameter. DO NOT CHANGE IT.
+        back = MEDCouplingUMesh('crh8_rse3', 2)
+        coo = DataArrayDouble([(-31.313754538446631,-32.512813836330515),(-31.531462871779969,-32.135731941766032),(-31.422608705113298,-32.324272889048274),(-31.690836433011114,-32.295105502997181),(-31.621640616088342,-32.204927758688783),(-31.502295485728872,-32.403959669663848)])
+        back.setCoords(coo)
+        c = DataArrayInt([32, 0, 3, 1, 5, 4, 2])
+        cI = DataArrayInt([0, 7])
+        back.setConnectivity(c, cI)
+
+        tool = MEDCouplingUMesh('TA-536193G_expl_20181022_merged', 2)
+        coo = DataArrayDouble([(-29.918137808525149,-26.883223901634544),(-32.919909136264039,-26.939612990540404),(-27.866900000000001,-28.016680435212603),(-31.313800000000001,-32.512799999999999),(-27.866900000000001,-28.933918793630923)])
+        tool.setCoords(coo)
+        c = DataArrayInt([5, 1, 0, 3, 5, 0, 2, 3, 5, 4, 3, 2])
+        cI = DataArrayInt([0, 4, 8, 12])
+        tool.setConnectivity(c, cI)
+
+        inter, res2Back, res2Tool = MEDCouplingUMesh.Intersect2DMeshes(back, tool, eps)
+
+        self.assertEqual(inter.getNodalConnectivity().getValues(), [5, 14, 13, 11, 12, 5, 13, 15, 11, 32, 12, 3, 1, 14, 16, 17, 18, 19, 5, 15, 0, 11])
+        self.assertEqual(inter.getNodalConnectivityIndex().getValues(), [0, 5, 9, 18, 22])
+        self.assertEqual(res2Back.getValues(), [0, 0, 0, 0])
+        self.assertEqual(res2Tool.getValues(), [0, 1, -1, -1])
+        pass
+
     def testSwig2Intersect2DMeshWith1DLine1(self):
         """A basic test with no colinearity between m1 and m2."""
         i=MEDCouplingIMesh("mesh",2,[5,5],[0.,0.],[1.,1.])