]> SALOME platform Git repositories - tools/medcoupling.git/commitdiff
Salome HOME
[EDF27860] : MEDCouplingUMesh.getCellsContainingPoints eps parameter specifies a... agy/edf27860
authorAnthony Geay <anthony.geay@edf.fr>
Mon, 14 Aug 2023 13:54:18 +0000 (15:54 +0200)
committerAnthony Geay <anthony.geay@edf.fr>
Mon, 14 Aug 2023 13:54:18 +0000 (15:54 +0200)
src/INTERP_KERNEL/InterpolationUtils.hxx
src/INTERP_KERNEL/PointLocatorAlgos.txx
src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py

index 5b58de09bb6e962d729a862e3441841ec2707f8f..3638f9c64840479888a28839297c8c5a2cf87c4b 100644 (file)
@@ -1293,14 +1293,18 @@ namespace INTERP_KERNEL
     return normFact;
   }
   
-  /*! Computes the triple product (XA^XB).XC (in 3D)*/
-  inline double TripleProduct(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]={ 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] };
     
     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;
index 4d4d1b8b766ac28c43abf2e727259057cb2c7d43..3ccb2c391623a0d1aaa651fe47aedb9edf19888a 100644 (file)
@@ -219,14 +219,15 @@ namespace INTERP_KERNEL
       (ptToTest[2]-centerPt[2])/normalizeFact};
       for (int iface=0; iface<nbfaces; iface++)
         {
-          double vol = ( TripleProduct(ptToTestNorm,
+          double lengthNorm = ( TripleProduct(
           ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 3,
           ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 6,
-          ptsOfTetrahedrizedPolyhedron.get() + 9*iface) ) / normalizeFact;
+          ptToTestNorm,
+          ptsOfTetrahedrizedPolyhedron.get() + 9*iface) );
           // assigne sign[iface] : 0 means on the face within eps, 1 or -1 means far from tetra representing face
-          if( vol<-eps )
+          if( lengthNorm<-eps )
             sign[iface]=-1;
-          else if( vol>eps )
+          else if( lengthNorm>eps )
             sign[iface]=1;
           else
             sign[iface]=0;
index c06cb28386ad625c81412162a7c014db9f2d4164..6d2ad7a88fe6025c15336e6cb6086beb81449664 100755 (executable)
@@ -4356,7 +4356,7 @@ mcIdType MEDCouplingUMesh::getCellContainingPoint(const double *pos, double eps)
  *          point, for more points getCellsContainingPoints() is recommended as it is
  *          faster.
  *  \param [in] pos - array of coordinates of the ball central point.
- *  \param [in] eps - ball radius.
+ *  \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 the found cells. It is cleared
  *         before inserting ids.
  *  \throw If the coordinates array is not set.
@@ -4423,7 +4423,7 @@ void MEDCouplingUMesh::getCellsContainingPointsZeAlg(const double *pos, mcIdType
  *         X0,Y0,Z0,X1,Y1,Z1,... Size of the array must be \a
  *         this->getSpaceDimension() * \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
index 2a75dce56becf801ec0017baedcc8cbd11f12cfb..c45f74141b3e8a558b556b8bfa4ea15a985a826d 100644 (file)
@@ -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