From: Anthony Geay Date: Mon, 14 Aug 2023 08:43:38 +0000 (+0200) Subject: [EDF27860] : MEDCouplingUMesh.getCellsContainingPoints eps parameter specifies a... X-Git-Tag: V9_12_0a1~2 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=1823c2a86a50582f8db28250144879660a76b371;p=tools%2Fmedcoupling.git [EDF27860] : MEDCouplingUMesh.getCellsContainingPoints eps parameter specifies a balls radius relative to the conidered cell --- diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index e286f8e96..3638f9c64 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -1253,27 +1253,61 @@ namespace INTERP_KERNEL } return result; } + + /*! + * This method normalize input "tetrahedrized polyhedron" to put it around 0,0,0 point and with the normalization factor. + * + * \param [in,out] ptsOfTetrahedrizedPolyhedron nbfaces * 3 * 3 vector that store in full interlace all the 3 first points of each face of the input polyhedron + * \param [in] nbfaces number of faces in the tetrahedrized polyhedron to be normalized + * \param [out] centerPt the center of input tetrahedrized polyhedron + * \return the normalization factor corresponding to the max amplitude among all nbfaces*3 input points and among X, Y and Z axis. + */ + inline double NormalizeTetrahedrizedPolyhedron(double *ptsOfTetrahedrizedPolyhedron, int nbfaces, double centerPt[3]) + { + centerPt[0] = 0.0; centerPt[1] = 0.0; centerPt[2] = 0.0; + double minMax[6]={ std::numeric_limits::max(), -std::numeric_limits::max(), + std::numeric_limits::max(), -std::numeric_limits::max(), + std::numeric_limits::max(), -std::numeric_limits::max() }; + for(int iPt = 0 ; iPt < 3 * nbfaces ; ++iPt) + { + for(int k = 0 ; k < 3 ; ++k) + { + minMax[2*k] = std::min(minMax[2*k],ptsOfTetrahedrizedPolyhedron[3*iPt+k]); + minMax[2*k+1] = std::max(minMax[2*k+1],ptsOfTetrahedrizedPolyhedron[3*iPt+k]); + } + } + double normFact = 0.0; + for(int k = 0 ; k < 3 ; ++k) + { + centerPt[k] = (minMax[2*k] + minMax[2*k+1]) / 2.0 ; + normFact = std::max(minMax[2*k+1] - minMax[2*k], normFact); + } + // renormalize into ptsOfTetrahedrizedPolyhedron + for(int iPt = 0 ; iPt < 3 * nbfaces ; ++iPt) + { + for(int k = 0 ; k < 3 ; ++k) + { + ptsOfTetrahedrizedPolyhedron[3*iPt+k] = ( ptsOfTetrahedrizedPolyhedron[3*iPt+k] - centerPt[k] ) / normFact; + } + } + return normFact; + } - /*! Computes the triple product (XA^XB).XC (in 3D)*/ - inline double triple_product(const double* A, const double*B, const double*C, const double*X) + /*! + * Computes the triple product (XA^XB).XC/(norm(XA^XB)) (in 3D) + * Returned value is equal to the distance (positive or negative depending of the position of C relative to XAB plane) between XAB plane and C point. + */ + inline double TripleProduct(const double *A, const double *B, const double *C, const double *X) { - double XA[3]; - XA[0]=A[0]-X[0]; - XA[1]=A[1]-X[1]; - XA[2]=A[2]-X[2]; - double XB[3]; - XB[0]=B[0]-X[0]; - XB[1]=B[1]-X[1]; - XB[2]=B[2]-X[2]; - double XC[3]; - XC[0]=C[0]-X[0]; - XC[1]=C[1]-X[1]; - XC[2]=C[2]-X[2]; + double XA[3]={ A[0]-X[0], A[1]-X[1], A[2]-X[2] }; + double XB[3]={ B[0]-X[0], B[1]-X[1], B[2]-X[2] }; + double XC[3]={ C[0]-X[0], C[1]-X[1], C[2]-X[2] }; - return - (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+ - (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+ - (XA[0]*XB[1]-XA[1]*XB[0])*XC[2]; + double XA_cross_XB[3] = {XA[1]*XB[2]-XA[2]*XB[1], XA[2]*XB[0]-XA[0]*XB[2], XA[0]*XB[1]-XA[1]*XB[0]}; + // norm is equal to double the area of the triangle + double norm = std::sqrt(XA_cross_XB[0]*XA_cross_XB[0]+XA_cross_XB[1]*XA_cross_XB[1]+XA_cross_XB[2]*XA_cross_XB[2]); + + return ( XA_cross_XB[0]*XC[0]+ XA_cross_XB[1]*XC[1] + XA_cross_XB[2]*XC[2] ) / norm; } /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a common point.*/ diff --git a/src/INTERP_KERNEL/PointLocatorAlgos.txx b/src/INTERP_KERNEL/PointLocatorAlgos.txx index 1365100e3..3ccb2c391 100644 --- a/src/INTERP_KERNEL/PointLocatorAlgos.txx +++ b/src/INTERP_KERNEL/PointLocatorAlgos.txx @@ -31,6 +31,7 @@ #include #include #include +#include namespace INTERP_KERNEL { @@ -119,7 +120,7 @@ namespace INTERP_KERNEL here XA^XC and XC^XB have different signs*/ const int SPACEDIM=MyMeshType::MY_SPACEDIM; - int* sign = new int[nbEdges]; + std::unique_ptr sign( new char[nbEdges] ); for (mcIdType iedge=0; iedge sign( new char[nbfaces] ); + std::unique_ptr connOfSon( new ConnType[conn_elem_sz] ); + std::unique_ptr ptsOfTetrahedrizedPolyhedron( new double[3*3*nbfaces] ); for (int iface=0; iface::coo2C(connOfSon[0])); + std::copy(AA,AA+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface); const double* BB=coords+SPACEDIM*(OTT::coo2C(connOfSon[1])); + std::copy(BB,BB+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface+3); const double* CC=coords+SPACEDIM*(OTT::coo2C(connOfSon[2])); - double Vol=triple_product(AA,BB,CC,ptToTest); - if (Vol<-eps) + std::copy(CC,CC+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface+6); + } + double centerPt[3]; + double normalizeFact = NormalizeTetrahedrizedPolyhedron(ptsOfTetrahedrizedPolyhedron.get(),nbfaces,centerPt); + double ptToTestNorm[3] = {(ptToTest[0]-centerPt[0])/normalizeFact, + (ptToTest[1]-centerPt[1])/normalizeFact, + (ptToTest[2]-centerPt[2])/normalizeFact}; + for (int iface=0; ifaceeps) + else if( lengthNorm>eps ) sign[iface]=1; - else - sign[iface]=0; + else + sign[iface]=0; } - bool ret=decideFromSign(sign, nbfaces); - delete [] sign; - delete [] connOfSon; + bool ret=decideFromSign(sign.get(), nbfaces); return ret; } @@ -273,10 +288,10 @@ namespace INTERP_KERNEL return isElementContainsPoint(x,type,coords,conn_elem,conn_elem_sz,eps); } - static bool decideFromSign(const int* sign, mcIdType nbelem) + static bool decideFromSign(const char *sign, mcIdType nbelem) { - int min_sign = 1; - int max_sign = -1; + char min_sign = 1; + char max_sign = -1; for (int i=0; igetSpaceDimension() * \a nbOfPoints * \param [in] nbOfPoints - number of points to locate within \a this mesh. - * \param [in] eps - radius of balls (i.e. the precision). + * \param [in] eps - ball radius (i.e. the precision) relative to max dimension along X, Y or Z of the the considered cell bounding box. * \param [out] elts - vector returning ids of found cells. * \param [out] eltsIndex - an array, of length \a nbOfPoints + 1, * dividing cell ids in \a elts into groups each referring to one diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index 2a75dce56..c45f74141 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1305,6 +1305,31 @@ class MEDCouplingBasicsTest7(unittest.TestCase): c,d = DataArrayInt.ExtractFromIndexedArrays( cell, pE3[1], pE3[2] ) self.assertTrue( ref_c.isEqual(c) ) self.assertTrue( ref_d.isEqual(d) ) + + def testGetCellContainingPointRelativeEps(self): + """ + See EDF27860 : This test checks that detection of point inside a cell works by normalizing cell around origin with factor equal to the max delta of bbox along axis X, Y or Z. + """ + # in this test cell is vuluntary far from origin {15260.775604514516, 11197.646906189217, 14187.820484060947} + # and caracteritic size is ~ 1500 + coo = DataArrayDouble( [(14724.199858870656, 11928.888084722483, 14442.32726944039), (14788.407409534622, 11992.60694822231, 14453.86181555231), (15572.175148726046, 10798.586790270576, 14471.54225356788), (15643.898717334796, 10853.094666047728, 14477.233802854305), (15005.31495255754, 11573.261110174888, 13933.313698681504), (15070.29423166349, 11636.377758513776, 13946.650959030132), (15797.351350158377, 10466.40572765595, 13965.524190108257), (15869.808770928525, 10519.99285973948, 13972.419352086607), (15273.866774426142, 11216.458197414971, 13433.169979717744), (15340.421031616577, 11277.882145177837, 13446.53598386297), (16013.382514001762, 10132.795887638129, 13465.184281842226), (16086.979064572806, 10184.802292369684, 13472.147425473782)] ) + m = MEDCouplingUMesh("",3) + m.setCoords(coo) + m.allocateCells() + m.insertNextCell(NORM_TETRA4,[0,5,4,6]) + m.insertNextCell(NORM_TETRA4,[4,5,9,7]) + + ##### See EDF2760 pt is outside cell 0 (6e-4) and 1 (8e-4) + pt = DataArrayDouble([(15263.41200205526, 11314.957094727113, 13950.0)]) + a,b = m.getCellsContainingPoints(pt,1e-3) + self.assertTrue(a.isEqual(DataArrayInt([0,1]))) + self.assertTrue(b.isEqual(DataArrayInt([0,2]))) + + # by shifting pt by 10 along Z pt in only inside cell # 0 + pt += [0,0,10] + a1,b1 = m.getCellsContainingPoints(pt,1e-3) + self.assertTrue(a1.isEqual(DataArrayInt([0]))) + self.assertTrue(b1.isEqual(DataArrayInt([0,1]))) pass