]> SALOME platform Git repositories - tools/medcoupling.git/commitdiff
Salome HOME
Bug fix: Intersect2DMeshWith1DLine now correctly handling closed lines ...
authorabn <adrien.bruneton@cea.fr>
Mon, 9 Apr 2018 14:16:16 +0000 (16:16 +0200)
committerabn <adrien.bruneton@cea.fr>
Wed, 18 Apr 2018 14:46:40 +0000 (16:46 +0200)
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx
src/MEDCoupling/MEDCouplingUMesh_intersection.cxx
src/MEDCoupling_Swig/MEDCouplingIntersectTest.py

index 43eed2a8727c519d5f58b26da8a5beebb160cef4..c847530f261315bd29d4ed2e89c7ac587b355a71 100644 (file)
@@ -124,10 +124,10 @@ namespace INTERP_KERNEL
     Node *getNodeAndReleaseIt() { Node *tmp=_node; _node=0; return tmp; }
     ~IntersectElement();
   private:
-    bool _1S;
-    bool _1E;
-    bool _2S;
-    bool _2E;
+    bool _1S;  // true if starting point of edge 1 is located exactly on edge 2 (not nearby)
+    bool _1E;  // true if ending point of edge 1 is located exactly on edge 2 (not nearby)
+    bool _2S;  // true if starting point of edge 2 is located exactly on edge 1 (not nearby)
+    bool _2E;  // true if ending point of edge 2 is located exactly on edge 1 (not nearby)
     double _chararct_val_for_e1;
     double _chararct_val_for_e2;
     Node *_node;
index 3cb86c501b6c21266fc3fe9ccd0c4b740b9956c9..ef169e3d6a616bef3f6ee495bc21e433c7db0d25 100644 (file)
@@ -293,6 +293,7 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other,
 {
   double xBaryBB, yBaryBB;
   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
+
   //
   IteratorOnComposedEdge it1(this),it3(&other);
   MergePoints merge;
@@ -300,13 +301,13 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other,
   ComposedEdge *c2=new ComposedEdge;
   int i=0;
   std::map<INTERP_KERNEL::Node *,int> mapAddCoo;
-  for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other' _sub_edges
+  for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other->_sub_edges'
     {
       QuadraticPolygon otherTmp;
       ElementaryEdge* curE3=it3.current();
       otherTmp.pushBack(new ElementaryEdge(curE3->getPtr(),curE3->getDirection())); curE3->getPtr()->incrRef();
       IteratorOnComposedEdge it2(&otherTmp);
-      for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'other->_sub_edge'
+      for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'otherTmp->_sub_edge'
         {
           ElementaryEdge* curE2=it2.current();
           if(!curE2->isThereStartPoint())
index 9aaae278b67bf0b8b2775b9025da69f56c2dea6a..6f29181487dc2b3382be96e39c9d0d784ba96e45 100644 (file)
@@ -423,7 +423,7 @@ bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1,
  * This method has 4 inputs :
  *  - a mesh 'm1' with meshDim==1 and a SpaceDim==2
  *  - a mesh 'm2' with meshDim==1 and a SpaceDim==2
- *  - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids randomly sorted.
+ *  - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm2' the splitting node ids randomly sorted.
  * The aim of this method is to sort the splitting nodes, if any, and to put them in 'intersectEdge' output parameter based on edges of mesh 'm2'
  * Nodes end up lying consecutively on a cutted edge.
  * \param m1 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method.
@@ -998,9 +998,9 @@ CellInfo& VectorOfCellInfo::get(int pos)
  *
  * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
  *
- * \param [in] allEdges a list of pairs (beginNode, endNode). Linked with \a allEdgesPtr to get the equation of edge.
+ * \param [in] allEdges a list of pairs (beginNode, endNode). Represents all edges (already cut) in the single 2D cell being handled here. Linked with \a allEdgesPtr to get the equation of edge.
  */
-MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
+MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
                                          MCAuto<DataArrayInt>& idsLeftRight)
 {
   int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
@@ -1020,6 +1020,32 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *spl
   idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
   int *idsLeftRightPtr(idsLeftRight->getPointer());
   VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
+
+  // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
+  // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
+  // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
+  // of the connectivity.
+  MCAuto <DataArrayInt> renumb(DataArrayInt::New());
+  renumb->alloc(nbCellsInSplitMesh1D,1);
+  const int * renumbP(renumb->begin());
+
+  int i, first=cSplitPtr[1];
+  // Follow 1D line backward as long as it is connected:
+  for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
+    first=cSplitPtr[ciSplitPtr[i]+1];
+  if (i < nbCellsInSplitMesh1D-1)
+    {
+      // Build circular permutation to shift consecutive edges together
+      renumb->iota(i+1);
+      renumb->applyModulus(nbCellsInSplitMesh1D);
+      splitMesh1D->renumberCells(renumbP, false);
+      cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
+      ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
+    }
+  else
+    renumb->iota();
+  //
+
   for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
     {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
       int iEnd(iStart);
@@ -1033,7 +1059,7 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *spl
         }
       if(iEnd<nbCellsInSplitMesh1D)
         iEnd++;
-      //
+
       MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
       int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
       //
@@ -1052,11 +1078,15 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *spl
       iStart=iEnd;
     }
   for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
-    pool.feedEdgeInfoAt(eps,mm,offset,idsLeftRightPtr+2*mm);
+    pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
+
   return pool.getZeMesh().retn();
 }
 
-MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, const MEDCouplingUMesh *splitMesh1D,
+/*
+ * splitMesh1D is an input parameter but might have its cells renumbered.
+ */
+MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
                                      const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
                                      MCAuto<DataArrayInt>& idsLeftRight)
 {
index 36a1e55736a30e31440948e3069483b670ab7ecc..7c2eae6f7fa796c58bb33d2d633bf69720258c9f 100644 (file)
@@ -892,6 +892,33 @@ class MEDCouplingIntersectTest(unittest.TestCase):
         self.assertTrue([0,1,2], c.getValues())
         self.assertEqual([2,1], d.getValues())
 
+    def testSwig2Intersect2DMeshWith1DLine18(self):
+        """ Rare case: a *closed* line used as a tool, with the closing point inside a 2D cell ... """
+        tool = MEDCouplingUMesh('circle', 1)
+        coo = DataArrayDouble([(39.35,0),(27.8247,27.8247),(2.40949e-15,39.35),(-27.8247,27.8247),(-39.35,4.81899e-15),(-27.8247,-27.8247),(-7.22848e-15,-39.35),(27.8247,-27.8247),(39.35,7.39805e-15)])
+        tool.setCoords(coo)
+        c = DataArrayInt([2, 3, 5, 8,    2, 5, 3, 4])
+        cI = DataArrayInt([0, 4, 8])
+        tool.setConnectivity(c, cI)
+
+        meh = MEDCouplingUMesh('meh', 2)
+        coo = DataArrayDouble([(-26.4275,36.6199),(-23.5868,31.6996),(-34.1861,41.0993),(-30.3383,25.0214),(-40.1861,30.707),(-35.2622,27.8642),(-37.1861,35.9032),(-30.3068,38.8596),(-25.0071,34.1598),(-26.9625,28.3605),(-25.7138,32.5128),(-27.354,36.4726),(-36.9138,32.5128),(-27.354,28.553),(-26.8908,36.5462),(-28.8461,26.7872)])
+        meh.setCoords(coo)
+        c = DataArrayInt([32, 0, 1, 3, 13, 11, 8, 9, 15, 10, 14,
+                             32, 3, 4, 2, 0, 11, 13, 5, 6, 7, 14, 12, 15])
+        cI = DataArrayInt([0, 11, 24])
+        meh.setConnectivity(c, cI)
+
+        res2D, res1D, m1, m2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(meh, tool, 1e-12)
+        self.assertEqual(4, res2D.getNumberOfCells())
+        self.assertEqual(res2D.getNodalConnectivity().getValues(),[32, 13, 11, 0, 1, 25, 19, 26, 33, 34, 35, 36, 37, 38, 39, 32, 3, 26, 19, 25, 40, 41, 42, 43,
+                                                                   32, 4, 2, 0, 11, 13, 26, 27, 44, 45, 46, 47, 48, 49, 50, 32, 3, 27, 26, 51, 52, 53])
+        self.assertEqual(res2D.getNodalConnectivityIndex().getValues(),[0, 15, 24, 39, 46])
+        self.assertEqual(res1D.getNodalConnectivity().getValues(),[2, 19, 25, 28, 2, 25, 21, 29, 2, 21, 27, 30, 2, 27, 26, 31, 2, 26, 19, 32])
+        self.assertEqual(res1D.getNodalConnectivityIndex().getValues(),[0, 4, 8, 12, 16, 20])
+        self.assertEqual(m1.getValues(), [0,0,1,1])
+        self.assertEqual(m2.getValues(), [0,1,   -1,-1,  -1,-1,   2,3,  0,1])
+
     def testSwig2Conformize2D1(self):
         eps = 1.0e-8
         coo = [0.,-0.5,0.,0.,0.5,0.,0.5,-0.5,0.25,