From 11f4044bee88716316b121e3f5ef6bb4010bfd67 Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 10 Dec 2013 10:37:33 +0000 Subject: [PATCH] Some missing tiny services for implementation of a general implementation of MEDCouplingUMesh::getCellsContainingPoints --- .../InterpKernelGeo2DComposedEdge.cxx | 93 +++++++++++++++++++ .../InterpKernelGeo2DComposedEdge.hxx | 2 + src/MEDCoupling/MEDCouplingUMesh.cxx | 40 +++++++- src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 12 +++ 4 files changed, 143 insertions(+), 4 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx index 8242278d7..fc64730db 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx @@ -246,6 +246,17 @@ void ComposedEdge::getBarycenterGeneral(double *bary) const _sub_edges.back()->getBarycenter(bary,w); } +double ComposedEdge::normalizeMe(double& xBary, double& yBary) +{ + Bounds b; + b.prepareForAggregation(); + fillBounds(b); + double dimChar=b.getCaracteristicDim(); + b.getBarycenter(xBary,yBary); + applyGlobalSimilarity(xBary,yBary,dimChar); + return dimChar; +} + double ComposedEdge::normalize(ComposedEdge *other, double& xBary, double& yBary) { Bounds b; @@ -415,6 +426,9 @@ void ComposedEdge::getBarycenter(double *bary, double& weigh) const bary[1]/=weigh; } +/*! + * \sa ComposedEdge::isInOrOut2 + */ bool ComposedEdge::isInOrOut(Node *nodeToTest) const { Bounds b; b.prepareForAggregation(); @@ -481,6 +495,85 @@ bool ComposedEdge::isInOrOut(Node *nodeToTest) const return ret; } +/*! + * This method is close to ComposedEdge::isInOrOut behaviour except that here EPSILON is taken into account to detect if it is IN or OUT. + * If \a nodeToTest is close to an edge in \a this, true will be returned even if it is outside informatically from \a this. + * This method makes the hypothesis that + * + * \sa ComposedEdge::isInOrOut + */ +bool ComposedEdge::isInOrOut2(Node *nodeToTest) const +{ + Bounds b; b.prepareForAggregation(); + fillBounds(b); + if(b.nearlyWhere((*nodeToTest)[0],(*nodeToTest)[1])==OUT) + return false; + // searching for e1 + std::set nodes; + getAllNodes(nodes); + std::set radialDistributionOfNodes; + std::set::const_iterator iter; + for(iter=nodes.begin();iter!=nodes.end();iter++) + radialDistributionOfNodes.insert(nodeToTest->getSlope(*(*iter))); + std::vector radialDistrib(radialDistributionOfNodes.begin(),radialDistributionOfNodes.end()); + radialDistributionOfNodes.clear(); + std::vector radialDistrib2(radialDistrib.size()); + copy(radialDistrib.begin()+1,radialDistrib.end(),radialDistrib2.begin()); + radialDistrib2.back()=M_PI+radialDistrib.front(); + std::vector radialDistrib3(radialDistrib.size()); + std::transform(radialDistrib2.begin(),radialDistrib2.end(),radialDistrib.begin(),radialDistrib3.begin(),std::minus()); + std::vector::iterator iter3=max_element(radialDistrib3.begin(),radialDistrib3.end()); + int i=iter3-radialDistrib3.begin(); + // ok for e1 - Let's go. + EdgeInfLin *e1=new EdgeInfLin(nodeToTest,radialDistrib[i]+radialDistrib3[i]/2.); + double ref=e1->getCharactValue(*nodeToTest); + std::set< IntersectElement > inOutSwitch; + for(std::list::const_iterator iter4=_sub_edges.begin();iter4!=_sub_edges.end();iter4++) + { + ElementaryEdge *val=(*iter4); + if(val) + { + Edge *e=val->getPtr(); + std::auto_ptr intersc(Edge::BuildIntersectorWith(e1,e)); + bool obviousNoIntersection,areOverlapped; + intersc->areOverlappedOrOnlyColinears(0,obviousNoIntersection,areOverlapped); + if(obviousNoIntersection) + { + continue; + } + if(!areOverlapped) + { + std::list< IntersectElement > listOfIntesc=intersc->getIntersectionsCharacteristicVal(); + for(std::list< IntersectElement >::iterator iter2=listOfIntesc.begin();iter2!=listOfIntesc.end();iter2++) + if((*iter2).isIncludedByBoth()) + inOutSwitch.insert(*iter2); + } + //if overlapped we can forget + } + else + throw Exception("Invalid use of ComposedEdge::isInOrOut : only one level supported !"); + } + e1->decrRef(); + bool ret=false; + for(std::set< IntersectElement >::iterator iter4=inOutSwitch.begin();iter4!=inOutSwitch.end();iter4++) + { + double val((*iter4).getVal1()); + if(fabs(val-ref)>=QUADRATIC_PLANAR::_precision) + { + if(valgetLoc()==ON_1) + ret=!ret; + } + else + break; + } + else + return true; + } + return ret; +} + /*bool ComposedEdge::isInOrOut(Node *aNodeOn, Node *nodeToTest) const { diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx index 206a7017e..545cc4544 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx @@ -58,6 +58,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT double getHydraulicDiameter() const; INTERPKERNEL_EXPORT void getBarycenter(double *bary) const; INTERPKERNEL_EXPORT void getBarycenterGeneral(double *bary) const; + INTERPKERNEL_EXPORT double normalizeMe(double& xBary, double& yBary); INTERPKERNEL_EXPORT double normalize(ComposedEdge *other, double& xBary, double& yBary); INTERPKERNEL_EXPORT double normalizeExt(ComposedEdge *other, double& xBary, double& yBary); INTERPKERNEL_EXPORT void unApplyGlobalSimilarityExt(ComposedEdge& other, double xBary, double yBary, double fact); @@ -89,6 +90,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT bool changeStartNodeWith(Node *node) const; INTERPKERNEL_EXPORT void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; INTERPKERNEL_EXPORT bool isInOrOut(Node *nodeToTest) const; + INTERPKERNEL_EXPORT bool isInOrOut2(Node *nodeToTest) const; INTERPKERNEL_EXPORT bool getDirection() const; INTERPKERNEL_EXPORT bool intresincEqCoarse(const Edge *other) const; private: diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 28490ba3b..3938f391d 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -4021,6 +4021,9 @@ void MEDCouplingUMesh::DistanceToPoint2DCurveAlg(const double *pt, const int *ce /*! * Finds cells in contact with a ball (i.e. a point with precision). + * For speed reasons, the INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6 and INTERP_KERNEL::NORM_QUAD8 cells are considered as convex cells to detect if a point is IN or OUT. + * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method. + * * \warning This method is suitable if the caller intends to evaluate only one * point, for more points getCellsContainingPoints() is recommended as it is * faster. @@ -4042,6 +4045,8 @@ int MEDCouplingUMesh::getCellContainingPoint(const double *pos, double eps) cons /*! * Finds cells in contact with a ball (i.e. a point with precision). + * For speed reasons, the INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6 and INTERP_KERNEL::NORM_QUAD8 cells are considered as convex cells to detect if a point is IN or OUT. + * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method. * \warning This method is suitable if the caller intends to evaluate only one * point, for more points getCellsContainingPoints() is recommended as it is * faster. @@ -4211,10 +4216,35 @@ void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const d myTree.getIntersectingElems(bb,candidates); for(std::vector::const_iterator iter=candidates.begin();iter!=candidates.end();iter++) { - int sz=connI[(*iter)+1]-connI[*iter]-1; - if(INTERP_KERNEL::PointLocatorAlgos >::isElementContainsPoint(pos+i*SPACEDIM, - (INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]], - coords,conn+connI[*iter]+1,sz,eps)) + int sz(connI[(*iter)+1]-connI[*iter]-1); + INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]]); + bool status(false); + if(ct!=INTERP_KERNEL::NORM_POLYGON && ct!=INTERP_KERNEL::NORM_QPOLYG) + status=INTERP_KERNEL::PointLocatorAlgos >::isElementContainsPoint(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps); + else + { + if(SPACEDIM!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPointsAlg : not implemented yet for POLYGON and QPOLYGON in spaceDim 3 !"); + INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; + INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; + std::vector nodes(sz); + INTERP_KERNEL::QuadraticPolygon *pol(0); + for(int j=0;jnormalizeMe(b,c); n->applySimilarity(b,c,a); + status=pol->isInOrOut2(n); + delete pol; n->decrRef(); + } + if(status) { eltsIndexPtr[i+1]++; elts->pushBackSilent(*iter); @@ -4226,6 +4256,8 @@ void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const d * Finds cells in contact with several balls (i.e. points with precision). * This method is an extension of getCellContainingPoint() and * getCellsContainingPoint() for the case of multiple points. + * For speed reasons, the INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6 and INTERP_KERNEL::NORM_QUAD8 cells are considered as convex cells to detect if a point is IN or OUT. + * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method. * \param [in] pos - an array of coordinates of points in full interlace mode : * X0,Y0,Z0,X1,Y1,Z1,... Size of the array must be \a * this->getSpaceDimension() * \a nbOfPoints diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index cb3168eed..1ba9b5092 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -14000,6 +14000,18 @@ 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) + 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]) + m.checkCoherency2() + # + self.assertTrue(m.getCellsContainingPoint((0.4,-0.4),1e-12).isEqual(DataArrayInt([0]))) + 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 setUp(self): pass pass -- 2.39.2