}
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<double>::max(), -std::numeric_limits<double>::max(),
+ std::numeric_limits<double>::max(), -std::numeric_limits<double>::max(),
+ std::numeric_limits<double>::max(), -std::numeric_limits<double>::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.*/
#include <vector>
#include <set>
#include <limits>
+#include <memory>
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<char[]> sign( new char[nbEdges] );
for (mcIdType iedge=0; iedge<nbEdges; iedge++)
{
const double* A=cellPts+SPACEDIM*iedge;
else
sign[iedge]=0;
}
- bool ret=decideFromSign(sign, nbEdges);
- delete [] sign;
+ bool ret=decideFromSign(sign.get(), nbEdges);
return ret;
}
const NumberingPolicy numPol=MyMeshType::My_numPol;
int nbfaces = cmType.getNumberOfSons2(conn_elem,conn_elem_sz);
- int *sign = new int[nbfaces];
- ConnType *connOfSon = new ConnType[conn_elem_sz];
+ std::unique_ptr<char[]> sign( new char[nbfaces] );
+ std::unique_ptr<ConnType[]> connOfSon( new ConnType[conn_elem_sz] );
+ std::unique_ptr<double[]> ptsOfTetrahedrizedPolyhedron( new double[3*3*nbfaces] );
for (int iface=0; iface<nbfaces; iface++)
{
NormalizedCellType typeOfSon;
- cmType.fillSonCellNodalConnectivity2(iface,conn_elem,conn_elem_sz,connOfSon,typeOfSon);
+ cmType.fillSonCellNodalConnectivity2(iface,conn_elem,conn_elem_sz,connOfSon.get(),typeOfSon);
const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[0]));
+ std::copy(AA,AA+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface);
const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[1]));
+ std::copy(BB,BB+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface+3);
const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::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; iface<nbfaces; iface++)
+ {
+ double lengthNorm = ( TripleProduct(
+ ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 3,
+ ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 6,
+ 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( lengthNorm<-eps )
sign[iface]=-1;
- else if (Vol>eps)
+ 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;
}
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; i<nbelem;i++)
{
min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
* 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.
* 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
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
