From 0c599cf7b1a423f9fd933cbaea6d9eb264b3dd39 Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 10 Dec 2013 16:18:17 +0000 Subject: [PATCH] Bounding box for BBTree in MEDCouplingPointSet takes into account quadratic 2D cells --- .../Geometric2D/InterpKernelGeo2DBounds.hxx | 4 + src/MEDCoupling/MEDCoupling1GTUMesh.cxx | 8 +- src/MEDCoupling/MEDCoupling1GTUMesh.hxx | 4 +- src/MEDCoupling/MEDCouplingPointSet.hxx | 2 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 77 ++++++++++++++++++- src/MEDCoupling/MEDCouplingUMesh.hxx | 4 +- src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 17 +++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 6 +- 8 files changed, 112 insertions(+), 10 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.hxx index ed3c2f305..be590241c 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.hxx @@ -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); diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx index 3ad0c08d4..ca92e2fa0 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx @@ -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 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()); diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx index 93ba77249..805ba002d 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx @@ -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); diff --git a/src/MEDCoupling/MEDCouplingPointSet.hxx b/src/MEDCoupling/MEDCouplingPointSet.hxx index 412d83c7b..aa5622b19 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.hxx +++ b/src/MEDCoupling/MEDCouplingPointSet.hxx @@ -131,7 +131,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const; MEDCOUPLING_EXPORT void unserialization(const std::vector& tinyInfoD, const std::vector& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector& 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(); diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 3938f391d..068584839 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -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::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 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 nodes(sz); + INTERP_KERNEL::QuadraticPolygon *pol(0); + for(int j=0;jfillBounds(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 diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 75b0a7cea..204368c9f 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -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; diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 1ba9b5092..60d22ed44 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -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 diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 250d034f4..984571254 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -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& 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); -- 2.39.2