Salome HOME
Bug fix: Intersect2DMeshWith1DLine now correctly handling closed lines ...
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh_intersection.cxx
index 629ae62b335baf49feb01c24a9bf3846b4801d08..6f29181487dc2b3382be96e39c9d0d784ba96e45 100644 (file)
@@ -120,9 +120,9 @@ void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp,
   // only the quadratic point to deal with:
   if(linOrArc)
     {
-      if(stp-start>1)
+      if(stp-start>1)  // if we are covering more than one segment we need to create a new mid point
         {
-          int tmpSrt(connBg[start]),tmpEnd(connBg[stp]);
+          int tmpSrt(connBg[start]),tmpEnd(connBg[stp % nbOfEdges]);  // % to handle last seg.
           int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
           InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
           middles.push_back(tmp3+offset);
@@ -282,7 +282,7 @@ namespace MEDCoupling
   }
 
   /**
-   * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI).
+   * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI).
    */
   void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
                                         const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
@@ -337,7 +337,7 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg
       // 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
+          for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell with one single edge
             {
               cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn,typeOfSon);
               INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
@@ -357,7 +357,7 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg
         }
       // Now move forward:
       const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt);  // the first element to be inspected going forward
-      for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++)  // 2nd condition is to avoid ending with a cell wih one single edge
+      for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++)  // 2nd condition is to avoid ending with a cell with one single edge
         {
           cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn,typeOfSon); // get edge #j's connectivity
           INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
@@ -375,9 +375,9 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg
               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!
+      // The if clauses below are (voluntary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its beginning!
       if(nbOfTurn==0)
-        // at the begining of the connectivity (insert type)
+        // at the beginning of the connectivity (insert type)
         EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
       else if((nbOfHit+nbOfTurn) != (nbs-1))
         // in the middle
@@ -423,14 +423,14 @@ 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.
  * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
  * \param m2 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.
  * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
- * \param[out] intersectEdge the same content as subDiv, but correclty oriented.
+ * \param[out] intersectEdge the same content as subDiv, but correctly oriented.
  */
 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
                                            const std::vector<double>& addCoo,
@@ -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)
 {
@@ -1206,7 +1236,7 @@ void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const M
           for(std::size_t iii=0;iii<szz;iii++,itt++)
             { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
           // end of protection
-          // Performs egde cutting:
+          // Performs edge cutting:
           pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
           delete pol2;
           delete pol1;
@@ -1252,7 +1282,7 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c
  * (newly created) nodes corresponding to the edge intersections.
  * Output params:
  * @param[out] cr, crI connectivity of the resulting mesh
- * @param[out] cNb1, cNb2 correspondance arrays giving for the merged mesh the initial cells IDs in m1 / m2
+ * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
  * TODO: describe input parameters
  */
 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
@@ -1431,7 +1461,7 @@ void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int
 }
 
 /*!
- * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additionnal nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
+ * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additional nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
  *
  * \sa MEDCouplingUMesh::split2DCells
  */
@@ -1467,7 +1497,7 @@ void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataAr
 
 
 /*!
- * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additionnal nodes can be added at the end of coordinates array object.
+ * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
  *
  * \return  int - the number of new nodes created.
  * \sa MEDCouplingUMesh::split2DCells
@@ -1619,7 +1649,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1
 /*!
  * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
  * Thus the final result contains the aggregation of nodes of \a mesh2D, then nodes of \a mesh1D, then new nodes that are the result of the intersection
- * and finaly, in case of quadratic polygon the centers of edges new nodes.
+ * and finally, in case of quadratic polygon the centers of edges new nodes.
  * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
  *
  * \param [in] mesh2D - the 2D mesh (spacedim=meshdim=2) to be intersected using \a mesh1D tool. The mesh must be so that each point in the space covered by \a mesh2D
@@ -1813,7 +1843,7 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D,
  * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
  * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this this method
  * will suppress such edges to use sub edges in \a this. So this method does not add nodes in \a this if merged edges are both linear (INTERP_KERNEL::NORM_SEG2).
- * In the other cases new nodes can be created. If any are created, they will be appended at the end of the coordinates object before the invokation of this method.
+ * In the other cases new nodes can be created. If any are created, they will be appended at the end of the coordinates object before the invocation of this method.
  *
  * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
  * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
@@ -1952,7 +1982,7 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
  * If yes, the cell is "repaired" to minimize at most its number of edges. So this method do not change the overall shape of cells in \a this (with eps precision).
  * This method do not take care of shared edges between cells, so this method can lead to a non conform mesh (\a this). If a conform mesh is required you're expected
  * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
- * This method works on any 2D geometric types of cell (even static one). If a cell is touched its type becomes dynamic automaticaly. For 2D "repaired" quadratic cells
+ * This method works on any 2D geometric types of cell (even static one). If a cell is touched its type becomes dynamic automatically. For 2D "repaired" quadratic cells
  * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
  *
  * If the returned array is empty it means that nothing has changed in \a this (as if it were a const method). If the array is not empty the connectivity of \a this is modified
@@ -2088,7 +2118,7 @@ void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxC
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
   int d = distance(startPos, endPos);
   if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
-    modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end());  // insidePoints also contains start and end node. Those dont need to be inserted.
+    modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end());  // insidePoints also contains start and end node. Those don't need to be inserted.
   else
     modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
 }
@@ -2264,7 +2294,7 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
           {
             if (packsIds[jj] == -1)
               // The below should never happen - if a face is used several times, with a different layout of the nodes
-              // it means that is is already conform, so it is *not* hit by the algorithm. The algorithm only hits
+              // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
               // faces which are actually used only once, by a single cell. This is different for edges below.
               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
             else
@@ -2401,7 +2431,7 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
             mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
             idsGoodLine->begin(), idsGoodLine->end(),
             /*out*/insidePoints, hitSegs);
-        // Optim: smaller segments completly included in eIdx and not split won't need any further treatment:
+        // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
         for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
           hit[cands2[*its]] = true;