Salome HOME
Bounding box for BBTree in MEDCouplingPointSet takes into account quadratic 2D cells
authorageay <ageay>
Tue, 10 Dec 2013 16:18:17 +0000 (16:18 +0000)
committerageay <ageay>
Tue, 10 Dec 2013 16:18:17 +0000 (16:18 +0000)
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.hxx
src/MEDCoupling/MEDCoupling1GTUMesh.cxx
src/MEDCoupling/MEDCoupling1GTUMesh.hxx
src/MEDCoupling/MEDCouplingPointSet.hxx
src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDCoupling/MEDCouplingUMesh.hxx
src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
src/MEDCoupling_Swig/MEDCouplingCommon.i

index ed3c2f305a64b0c161e1c93ccd21bdedd5588583..be590241c3ca8ebf0cfb6e4e791f090dea0ae70a 100644 (file)
@@ -44,6 +44,10 @@ namespace INTERP_KERNEL
     Bounds():_x_min(0.),_x_max(0.),_y_min(0.),_y_max(0.) { }
     double &operator[](int i);
     const double& operator[](int i) const;
+    double getXMin() const { return _x_min; }
+    double getXMax() const { return _x_max; }
+    double getYMin() const { return _y_min; }
+    double getYMax() const { return _y_max; }
     double getDiagonal() const;
     void getBarycenter(double& xBary, double& yBary) const;
     void applySimilarity(double xBary, double yBary, double dimChar);
index 3ad0c08d4838603a29cb634d307f19917a9aa259..ca92e2fa09a7705d8683bee975ebd47159153356 100644 (file)
@@ -1960,12 +1960,14 @@ MEDCoupling1DGTUMesh *MEDCoupling1SGTUMesh::computeDualMesh2D() const
 /*!
  * This method aggregate the bbox of each cell and put it into bbox 
  *
+ * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
+ *                         For all other cases this input parameter is ignored.
  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
  * 
  * \throw If \a this is not fully set (coordinates and connectivity).
  * \throw If a cell in \a this has no valid nodeId.
  */
-DataArrayDouble *MEDCoupling1SGTUMesh::getBoundingBoxForBBTree() const
+DataArrayDouble *MEDCoupling1SGTUMesh::getBoundingBoxForBBTree(double arcDetEps) const
 {
   int spaceDim(getSpaceDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()),nbOfNodesPerCell(getNumberOfNodesPerCell());
   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
@@ -3352,12 +3354,14 @@ MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::buildSetInstanceFromThis(int spaceDi
 /*!
  * This method aggregate the bbox of each cell and put it into bbox parameter.
  * 
+ * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
+ *                         For all other cases this input parameter is ignored.
  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
  * 
  * \throw If \a this is not fully set (coordinates and connectivity).
  * \throw If a cell in \a this has no valid nodeId.
  */
-DataArrayDouble *MEDCoupling1DGTUMesh::getBoundingBoxForBBTree() const
+DataArrayDouble *MEDCoupling1DGTUMesh::getBoundingBoxForBBTree(double arcDetEps) const
 {
   checkFullyDefined();
   int spaceDim(getSpaceDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
index 93ba77249f3eadcd8139b89fb5c5dc2c32325720..805ba002d0155d9b50f8b4ffd3d13d7cc67f0eb5 100644 (file)
@@ -130,7 +130,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT void renumberNodesInConn(const int *newNodeNumbersO2N);
     MEDCOUPLING_EXPORT void fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const;
     MEDCOUPLING_EXPORT int getNumberOfNodesInCell(int cellId) const;
-    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree() const;
+    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const;
     // overload of MEDCoupling1GTUMesh
     MEDCOUPLING_EXPORT void checkCoherencyOfConnectivity() const;
     MEDCOUPLING_EXPORT void allocateCells(int nbOfCells=0);
@@ -219,7 +219,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT void renumberNodesInConn(const int *newNodeNumbersO2N);
     MEDCOUPLING_EXPORT void fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const;
     MEDCOUPLING_EXPORT int getNumberOfNodesInCell(int cellId) const;
-    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree() const;
+    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const;
     // overload of MEDCoupling1GTUMesh
     MEDCOUPLING_EXPORT void checkCoherencyOfConnectivity() const;
     MEDCOUPLING_EXPORT void allocateCells(int nbOfCells=0);
index 412d83c7bb09c99d96524f1471a8d8b4256a1004..aa5622b19b926eddfab5bfe2de5c7727141b8580 100644 (file)
@@ -131,7 +131,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const;
     MEDCOUPLING_EXPORT void unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
                                             const std::vector<std::string>& littleStrings);
-    MEDCOUPLING_EXPORT virtual DataArrayDouble *getBoundingBoxForBBTree() const = 0;
+    MEDCOUPLING_EXPORT virtual DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const = 0;
     MEDCOUPLING_EXPORT virtual DataArrayInt *getCellsInBoundingBox(const double *bbox, double eps) const = 0;
     MEDCOUPLING_EXPORT virtual DataArrayInt *getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps) = 0;
     MEDCOUPLING_EXPORT virtual DataArrayInt *zipCoordsTraducer();
index 3938f391d0d76db09e5424583601daffc5d90e8a..0685848390d3dfb70bd563d7e7737d3270991965 100644 (file)
@@ -6254,12 +6254,45 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const
 /*!
  * This method aggregate the bbox of each cell and put it into bbox parameter.
  * 
+ * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
+ *                         For all other cases this input parameter is ignored.
  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
  * 
  * \throw If \a this is not fully set (coordinates and connectivity).
  * \throw If a cell in \a this has no valid nodeId.
+ * \sa MEDCouplingUMesh::getBoundingBoxForBBTreeFast, MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic
  */
-DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree() const
+DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree(double arcDetEps) const
+{
+  int mDim(getMeshDimension()),sDim(getSpaceDimension());
+  if(mDim!=2 || sDim!=2)
+    return getBoundingBoxForBBTreeFast();
+  else
+    {
+      bool presenceOfQuadratic(false);
+      for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=_types.begin();it!=_types.end();it++)
+        {
+          const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(*it));
+          if(cm.isQuadratic())
+            presenceOfQuadratic=true;
+        }
+      if(presenceOfQuadratic)
+        return getBoundingBoxForBBTree2DQuadratic(arcDetEps);
+      else
+        return getBoundingBoxForBBTreeFast();
+    }
+}
+
+/*!
+ * This method aggregate the bbox of each cell only considering the nodes constituting each cell and put it into bbox parameter.
+ * So meshes having quadratic cells the computed bounding boxes can be invalid !
+ * 
+ * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
+ * 
+ * \throw If \a this is not fully set (coordinates and connectivity).
+ * \throw If a cell in \a this has no valid nodeId.
+ */
+DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const
 {
   checkFullyDefined();
   int spaceDim(getSpaceDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
@@ -6298,6 +6331,48 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree() const
   return ret.retn();
 }
 
+/*!
+ * This method aggregate the bbox regarding foreach 2D cell in \a this the whole shape. So this method is particulary useful for 2D meshes having quadratic cells
+ * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes.
+ * 
+ * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
+ * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
+ * \throw If \a this is not fully defined.
+ * \throw If \a this is not a mesh with meshDimension equal to 2.
+ * \throw If \a this is not a mesh with spaceDimension equal to 2.
+ */
+DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const
+{
+  checkFullyDefined();
+  int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
+  if(mDim!=2 || spaceDim!=2)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic : This method should be applied on mesh with mesh dimension equal to 2 and space dimension also equal to 2!");
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
+  double *bbox(ret->getPointer());
+  const double *coords(_coords->getConstPointer());
+  const int *conn(_nodal_connec->getConstPointer()),*connI(_nodal_connec_index->getConstPointer());
+  for(int i=0;i<nbOfCells;i++,bbox+=4,connI++)
+    {
+      const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
+      int sz(connI[1]-connI[0]-1);
+      INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=1e-12;
+      std::vector<INTERP_KERNEL::Node *> nodes(sz);
+      INTERP_KERNEL::QuadraticPolygon *pol(0);
+      for(int j=0;j<sz;j++)
+        {
+          int nodeId(conn[*connI+1+j]);
+          nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*2],coords[nodeId*2+1]);
+        }
+      if(!cm.isQuadratic())
+        pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
+      else
+        pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
+      INTERP_KERNEL::Bounds b; pol->fillBounds(b); delete pol;
+      bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); 
+    }
+  return ret.retn();
+}
+
 /// @cond INTERNAL
 
 namespace ParaMEDMEMImpl
index 75b0a7ceaf24266830701fd61eef8f3c0cc1466b..204368c9f97f94ca06714849111558d776b0ba03 100644 (file)
@@ -164,7 +164,9 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT DataArrayInt *convexEnvelop2D();
     MEDCOUPLING_EXPORT DataArrayInt *findAndCorrectBadOriented3DExtrudedCells();
     MEDCOUPLING_EXPORT DataArrayInt *findAndCorrectBadOriented3DCells();
-    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree() const;
+    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const;
+    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTreeFast() const;
+    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const;
     MEDCOUPLING_EXPORT MEDCouplingUMesh *buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy);
     MEDCOUPLING_EXPORT bool isFullyQuadratic() const;
     MEDCOUPLING_EXPORT bool isPresenceOfQuadratic() const;
index 1ba9b5092981f1ecd910bd824ceab58560c36ab3..60d22ed44f48d9fe461bbfc22fff3a62073df67b 100644 (file)
@@ -14000,8 +14000,8 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertEqual([3,12,5,6,14,16,23,100,23,1,0,6],d0i.getValues())
         pass
 
-    def testSwig2GetCellsContainingPointsForNonConvexPolygon(self):
-        coo=DataArrayDouble([-0.5,-0.5,-0.5,0.5,0.5,0.5,0.5,-0.5,0.,-0.5,0.,0.,0.5,0.],7,2)
+    def testSwig2GetCellsContainingPointsForNonConvexPolygon1(self):
+        coo=DataArrayDouble([-0.5,-0.5,-0.5,0.5,0.5,0.5,0.5,-0.5,0.,-0.5,0.,0.,0.5,0.,],7,2)
         m=MEDCouplingUMesh("Intersect2D",2) ; m.setCoords(coo) ; m.allocateCells()
         m.insertNextCell(NORM_POLYGON,[6,3,4,5])
         m.insertNextCell(NORM_POLYGON,[4,0,1,2,6,5])
@@ -14011,6 +14011,19 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertTrue(m.getCellsContainingPoint((-0.4,-0.4),1e-12).isEqual(DataArrayInt([1])))
         self.assertTrue(m.getCellsContainingPoint((0.,-0.4),1e-12).isEqual(DataArrayInt([0,1])))
         pass
+    
+    def testSwig2GetCellsContainingPointsForNonConvexPolygon2(self):
+        coo=DataArrayDouble([-0.5,-0.5,-0.5,0.5,0.5,0.5,0.5,-0.5,-2.0816681711721685e-17,-2.0816681711721685e-17,-0.17677669529663687,0.1767766952966369,0.,0.5,0.5,0.,0.17677669529663684,-0.17677669529663692,0.17677669529663692,0.17677669529663684,-0.17677669529663692,-0.17677669529663687,0.,-0.5,-0.5,0.,0.33838834764831843,-0.3383883476483185,-0.33838834764831843,0.33838834764831843,-0.21213203435596423,0.21213203435596426,0.2121320343559642,-0.2121320343559643,0.21213203435596426,0.2121320343559642,-0.21213203435596423,-0.21213203435596428,0.3560660171779821,-0.35606601717798214,-0.35606601717798214,0.35606601717798214,0.19445436482630052,-0.19445436482630063,-0.19445436482630055,0.19445436482630057,0.,0.27],24,2)
+        m=MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
+        m.insertNextCell(NORM_QPOLYG,[8,5,4,9])
+        m.insertNextCell(NORM_QPOLYG,[5,8,4,10])
+        m.insertNextCell(NORM_QPOLYG,[16,8,5,15,21,9,22,17])
+        m.insertNextCell(NORM_QPOLYG,[15,1,2,3,16,20,6,7,19,17])
+        m.insertNextCell(NORM_QPOLYG,[15,5,8,16,22,10,21,18])
+        m.insertNextCell(NORM_QPOLYG,[16,3,0,1,15,19,11,12,20,18])
+        m.checkCoherency2()
+        self.assertTrue(m.getCellsContainingPoint([0.,0.27],1e-12).isEqual(DataArrayInt([2])))
+        pass
 
     def setUp(self):
         pass
index 250d034f49baa418bf5b74084f26a442a0bac734..9845712543dc131491904492896d70440473d09d 100644 (file)
@@ -267,6 +267,8 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingUMesh::ComputeRangesFromTypeDistribution;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildUnionOf2DMesh;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildUnionOf3DMesh;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTreeFast;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic;
 %newobject ParaMEDMEM::MEDCouplingUMeshCellByTypeEntry::__iter__;
 %newobject ParaMEDMEM::MEDCouplingUMeshCellEntry::__iter__;
 %newobject ParaMEDMEM::MEDCoupling1GTUMesh::New;
@@ -917,7 +919,7 @@ namespace ParaMEDMEM
       virtual void checkFullyDefined() const throw(INTERP_KERNEL::Exception);
       virtual bool isEmptyMesh(const std::vector<int>& tinyInfo) const throw(INTERP_KERNEL::Exception);
       virtual MEDCouplingPointSet *deepCpyConnectivityOnly() const throw(INTERP_KERNEL::Exception);
-      virtual DataArrayDouble *getBoundingBoxForBBTree() const throw(INTERP_KERNEL::Exception);
+      virtual DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception);
       %extend 
          {
            std::string __str__() const throw(INTERP_KERNEL::Exception)
@@ -1534,6 +1536,8 @@ namespace ParaMEDMEM
     DataArrayInt *convertNodalConnectivityToStaticGeoTypeMesh() const throw(INTERP_KERNEL::Exception);
     DataArrayInt *buildUnionOf2DMesh() const throw(INTERP_KERNEL::Exception);
     DataArrayInt *buildUnionOf3DMesh() const throw(INTERP_KERNEL::Exception);
+    DataArrayDouble *getBoundingBoxForBBTreeFast() const throw(INTERP_KERNEL::Exception);
+    DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da) throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);