Salome HOME
New functionality: colinearizeKeepingConform2D().
authorabn <adrien.bruneton@cea.fr>
Mon, 6 Aug 2018 08:39:51 +0000 (10:39 +0200)
committerabn <adrien.bruneton@cea.fr>
Mon, 6 Aug 2018 08:40:40 +0000 (10:40 +0200)
src/INTERP_KERNEL/CellModel.cxx
src/MEDCoupling/MEDCouplingUMesh.hxx
src/MEDCoupling/MEDCouplingUMesh_intersection.cxx
src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py
src/MEDCoupling_Swig/MEDCouplingCommon.i

index 9d7e6e29a96c31d07a76501f8256efd7c2c4621b..8b45eaadc306f0b279289ae2e87b51aeb9740f09 100644 (file)
@@ -573,7 +573,7 @@ namespace INTERP_KERNEL
                 sonNodalConn[1]=nodalConn[(sonId+1)%lgth];
                 return 2;
               }
-            else
+            else  // NORM_QPOLYG
               {
                 sonNodalConn[0]=nodalConn[sonId];
                 sonNodalConn[1]=nodalConn[(sonId+1)%(lgth/2)];
index 32dd8bd236a5aa302d2078fafa78033678dc7ca8..a7be645d8e3104dbed85713059b89026701b1059 100644 (file)
@@ -236,6 +236,7 @@ namespace MEDCoupling
     MEDCOUPLING_EXPORT DataArrayDouble *computePlaneEquationOf3DFaces() const;
     MEDCOUPLING_EXPORT DataArrayInt *conformize2D(double eps);
     MEDCOUPLING_EXPORT DataArrayInt *colinearize2D(double eps);
+    MEDCOUPLING_EXPORT DataArrayInt *colinearizeKeepingConform2D(double eps);
     MEDCOUPLING_EXPORT DataArrayInt *conformize3D(double eps);
     MEDCOUPLING_EXPORT int split2DCells(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *midOpt=0, const DataArrayInt *midOptI=0);
     MEDCOUPLING_EXPORT static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da);
@@ -351,7 +352,7 @@ namespace MEDCoupling
     void buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf, const int *desc, const int *descIndx, const double *coords, double eps, std::vector<std::vector<int> >& res) const;
     void split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI);
     int split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI);
-    static bool Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords);
+    static bool Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, const std::map<int, bool>& forbiddenPoints, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords);
     static void ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayInt *nodalConnec, const DataArrayInt *nodalConnecIndex);
     static bool OrderPointsAlongLine(const double * coo, int startNode, int endNode,
                                                 const int * c, const int * cI, const int *idsBg, const int *endBg,
@@ -359,6 +360,7 @@ namespace MEDCoupling
     static void ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
                                       const std::vector<int>& insidePoints, std::vector<int>& modifiedFace);
     void attractSeg3MidPtsAroundNodesUnderground(double ratio, const int *nodeIdsBg, const int *nodeIdsEnd);
+    DataArrayInt *internalColinearize2D(double eps, bool stayConform);
   public:
     MEDCOUPLING_EXPORT static DataArrayInt *ComputeRangesFromTypeDistribution(const std::vector<int>& code);
     MEDCOUPLING_EXPORT static const int N_MEDMEM_ORDER=25;
index 6c863910ad4fbb62857dc54d451132091450c662..ba461e5a47ae9fb4fd0f4c0050df4023092004ef 100644 (file)
@@ -310,14 +310,18 @@ namespace MEDCoupling
 /*!
  * Returns true if a colinearization has been found in the given cell. If false is returned the content pushed in \a newConnOfCell is equal to [ \a connBg , \a connEnd ) .
  * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
+ * \param forbiddenPoints the list of points that should not be removed in the process
  */
-bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
+bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset,
+                                         const std::map<int, bool>& forbiddenPoints,
+                                         DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
 {
   std::size_t sz(std::distance(connBg,connEnd));
   if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
   sz--;
   INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
+  INTERP_KERNEL::AutoPtr<int> tmpConn2(new int[sz]);
   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
   unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
   unsigned nbOfHit(0); // number of fusions operated
@@ -339,8 +343,13 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg
         {
           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));
+              cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn2,typeOfSon);
+              // Identify common point:
+              int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
+              auto itE(forbiddenPoints.end());
+              if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
+                break;
+              INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
               INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
               bool isColinear=eint->areColinears();
               if(isColinear)
@@ -353,14 +362,21 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg
               eCand->decrRef();
               if(!isColinear)
                 break;
+              // Update last connectivity
+              std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
             }
         }
       // 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 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));
+          cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn2,typeOfSon); // get edge #j's connectivity
+          // Identify common point:
+          int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
+          auto itE(forbiddenPoints.end());
+          if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
+            break;
+          INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
           INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
           bool isColinear(eint->areColinears());
           if(isColinear)
@@ -373,6 +389,8 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg
           eCand->decrRef();
           if(!isColinear)
               break;
+          // Update last connectivity
+          std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
         }
       //push [posBaseElt,posEndElt) in newConnOfCell using e
       // 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!
@@ -1997,6 +2015,27 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
  * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
  */
 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
+{
+  internalColinearize2D(eps, false);
+}
+
+/*!
+ * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
+ * In a given 2D cell, if two edges are colinear and the junction point is used by a third edge, the two edges will not be
+ * merged, contrary to colinearize2D().
+ *
+ * \sa MEDCouplingUMesh::colinearize2D
+ */
+DataArrayInt *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
+{
+  internalColinearize2D(eps, true);
+}
+
+
+/*!
+ * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
+ */
+DataArrayInt *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
 {
   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
   checkConsistencyLight();
@@ -2006,12 +2045,27 @@ DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
   int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
   const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
   MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
+  std::map<int, bool> forbiddenPoints;  // list of points that can not be removed (or it will break conformity)
+  if(stayConform) //
+    {
+      // A point that is used by more than 2 edges can not be removed without breaking conformity:
+      MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descI1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescI1(DataArrayInt::New());
+      MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
+      MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
+      MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
+      MCAuto<DataArrayInt> dsi(revDescI2->deltaShiftIndex());
+      MCAuto<DataArrayInt> ids(dsi->findIdsGreaterThan(2));
+      const int * cPtr(mDesc0D->getNodalConnectivity()->begin());
+      for(auto it = ids->begin(); it != ids->end(); it++)
+         forbiddenPoints[cPtr[2*(*it)+1]] = true;  // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
+    }
+
   MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
   const double *coords(_coords->begin());
   int *newciptr(newci->getPointer());
   for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
     {
-      if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,newc,appendedCoords))
+      if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
         ret->pushBackSilent(i);
       newciptr[1]=newc->getNumberOfTuples();
     }
index ebb5bfe4d1ca5cf713a780946ce11d764937dcf4..b5b326bc48c68a98f60aca749294c7def4b359d0 100644 (file)
@@ -2001,6 +2001,26 @@ class MEDCouplingBasicsTest5(unittest.TestCase):
         self.assertTrue(m.getNodalConnectivityIndex().isEqual(DataArrayInt([0,7])))
         pass
 
+    def testSwig2ColinearizeKeepingConform2D1(self):
+        eps = 1.0e-6
+        # Just to get a nice coords array ...
+        mm = MEDCouplingCMesh(); arr = DataArrayDouble([0.0, 1.0,2.0])
+        mm.setCoords(arr, arr);  mm = mm.buildUnstructured();   coo = mm.getCoords()
+         
+        mesh = MEDCouplingUMesh("M", 2)
+        mesh.setCoords(coo)
+        c = [NORM_POLYGON, 0,1,4,7,6,3,  NORM_QUAD4, 1,2,5,4,  NORM_QUAD4,4,5,8,7]
+        cI = [0, 7,12,17]
+        mm.setConnectivity(DataArrayInt(c),DataArrayInt(cI))
+        mm.checkConsistencyLight()
+        
+        mm.colinearizeKeepingConform2D(eps)
+        c = mm.getNodalConnectivity().getValues()
+        cI = mm.getNodalConnectivityIndex().getValues()
+        self.assertEqual(c, [NORM_POLYGON, 0, 1, 4, 7, 6, NORM_POLYGON, 1, 2, 5, 4, NORM_POLYGON, 4, 5, 8, 7])
+        self.assertEqual(cI, [0,6,11,16])
+        pass
+
     def testSwig2BoundingBoxForBBTree1(self):
         """ This test appears simple but it checks that bounding box are correctly computed for quadratic polygons. It can help a lot to reduce the amount of intersections !
         """
index 52461a4fc2eb88802d626bea3a02a55857dfc96e..40dcdafcdd7f7c56899067e8f3fac6ee1e910ce8 100644 (file)
@@ -325,6 +325,7 @@ using namespace INTERP_KERNEL;
 %newobject MEDCoupling::MEDCouplingUMesh::conformize2D;
 %newobject MEDCoupling::MEDCouplingUMesh::conformize3D;
 %newobject MEDCoupling::MEDCouplingUMesh::colinearize2D;
+%newobject MEDCoupling::MEDCouplingUMesh::colinearizeKeepingConform2D;
 %newobject MEDCoupling::MEDCouplingUMesh::rearrange2ConsecutiveCellTypes;
 %newobject MEDCoupling::MEDCouplingUMesh::sortCellsInMEDFileFrmt;
 %newobject MEDCoupling::MEDCouplingUMesh::getRenumArrForMEDFileFrmt;
@@ -1950,6 +1951,7 @@ namespace MEDCoupling
     DataArrayInt *conformize2D(double eps) throw(INTERP_KERNEL::Exception);
     DataArrayInt *conformize3D(double eps) throw(INTERP_KERNEL::Exception);
     DataArrayInt *colinearize2D(double eps) throw(INTERP_KERNEL::Exception);
+    DataArrayInt *colinearizeKeepingConform2D(double eps) throw(INTERP_KERNEL::Exception);
     void shiftNodeNumbersInConn(int delta) throw(INTERP_KERNEL::Exception);
     std::vector<bool> getQuadraticStatus() const throw(INTERP_KERNEL::Exception);
     DataArrayInt *findCellIdsOnBoundary() const throw(INTERP_KERNEL::Exception);