From 74b7797e44b8e15f1d24e48af3e67be29abce269 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Mon, 14 Aug 2023 10:43:38 +0200 Subject: [PATCH] [EDF27860] : Improve MEDCouplingUMesh.getCellsContainingPoints detection by changing the semantic of eps. Now in 3D Eps is relative distance to the 3D cell. --- src/INTERP_KERNEL/InterpolationUtils.hxx | 64 +++++++++++++++++------- src/INTERP_KERNEL/PointLocatorAlgos.txx | 48 +++++++++++------- 2 files changed, 78 insertions(+), 34 deletions(-) diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index e286f8e96..5b58de09b 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -1253,27 +1253,57 @@ 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) + 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]}; + 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..4d4d1b8b7 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( vol>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 +287,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; i