src/INTERP_KERNEL/PointLocatorAlgos.txx

   1 // Copyright (C) 2007-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20 #ifndef __POINTLOCATORALGOS_TXX__
  21 #define __POINTLOCATORALGOS_TXX__
  22
  23 #include "InterpolationUtils.hxx"
  24 #include "CellModel.hxx"
  25 #include "BBTree.txx"
  26 #include "VectorUtils.hxx"
  27
  28 #include "InterpKernelGeo2DNode.hxx"
  29 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
  30
  31 #include <list>
  32 #include <vector>
  33 #include <set>
  34 #include <limits>
  35 #include <memory>
  36
  37 namespace INTERP_KERNEL
  38 {
  39   class GenericPointLocatorAlgos
  40   {
  41   public:
  42     virtual ~GenericPointLocatorAlgos() { }
  43     virtual std::list<mcIdType> locates(const double* x, double eps) = 0;
  44   };
  45
  46   template<class MyMeshType>
  47   class PointLocatorAlgos: public GenericPointLocatorAlgos
  48   {
  49   private :
  50     double* _bb;
  51     BBTree<MyMeshType::MY_SPACEDIM,typename MyMeshType::MyConnType>* _tree;
  52     const MyMeshType& _mesh;
  53   public:
  54     PointLocatorAlgos(const MyMeshType& mesh):_mesh(mesh)
  55     {
  56       typedef typename MyMeshType::MyConnType ConnType;
  57       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
  58       const NumberingPolicy numPol=MyMeshType::My_numPol;
  59       ConnType nelem = _mesh.getNumberOfElements();
  60       _bb = new double[SPACEDIM*2*nelem];
  61       const ConnType* conn = _mesh.getConnectivityPtr();
  62       const ConnType* conn_index = _mesh.getConnectivityIndexPtr();
  63       const double* coords=_mesh.getCoordinatesPtr();
  64       for (ConnType i=0; i<nelem; i++)
  65         {
  66           for (int idim=0; idim<SPACEDIM; idim++)
  67             {
  68               _bb[2*(i*SPACEDIM+idim)]=std::numeric_limits<double>::max();
  69               _bb[2*(i*SPACEDIM+idim)+1]=-std::numeric_limits<double>::max();
  70             }
  71           for (ConnType index= conn_index[i]; index < conn_index[i+1];index++)
  72             {
  73               //coordelem points to the coordinates of the current node of the i-th element
  74               const double* coordelem = coords+OTT<ConnType,numPol>::ind2C(conn[OTT<ConnType,numPol>::ind2C(index)])*SPACEDIM;
  75
  76               //the bounding box is updated by checking wheher the node is at the min/max in exach dimension
  77               for (int idim=0; idim<SPACEDIM;idim++)
  78                 {
  79                   _bb[2*(i*SPACEDIM+idim)]=(coordelem[idim]<_bb[2*(i*SPACEDIM+idim)])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)];
  80                   _bb[2*(i*SPACEDIM+idim)+1]=(coordelem[idim]>_bb[2*(i*SPACEDIM+idim)+1])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)+1];
  81                 }
  82             }
  83         }
  84       _tree=new BBTree<SPACEDIM,typename MyMeshType::MyConnType>(_bb,0,0,nelem);
  85     }
  86
  87     ~PointLocatorAlgos()
  88     {
  89       delete[] _bb;
  90       delete _tree;
  91     }
  92
  93     //returns the list of elements that contains
  94     //the point pointed to by x
  95     std::list<typename MyMeshType::MyConnType> locates(const double* x, double eps)
  96     {
  97       typedef typename MyMeshType::MyConnType ConnType;
  98       const NumberingPolicy numPol=MyMeshType::My_numPol;
  99       std::vector<ConnType> candidates;
 100       _tree->getElementsAroundPoint(x,candidates);
 101       std::list<ConnType> retlist;
 102       for(unsigned int i=0; i< candidates.size(); i++)
 103         {
 104           ConnType ielem=candidates[i];
 105           if (elementContainsPoint(ielem,x,eps))
 106             retlist.push_back(OTT<ConnType,numPol>::indFC(ielem));
 107         }
 108       return retlist;
 109     }
 110
 111     static bool isElementContainsPointAlg2DSimple(const double *ptToTest, const double *cellPts, mcIdType nbEdges, double eps)
 112     {
 113       /* with dimension 2, it suffices to check all the edges
 114          and see if the sign of double products from the point
 115          is always the same.
 116                          C
 117                         / \
 118                        /   \
 119              Xo       /     \
 120                      A-------B
 121
 122          here XA^XC and XC^XB have different signs*/
 123       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 124       std::unique_ptr<char[]> sign( new char[nbEdges] );
 125       for (mcIdType iedge=0; iedge<nbEdges; iedge++)
 126         {
 127           const double* A=cellPts+SPACEDIM*iedge;
 128           const double* B=cellPts+SPACEDIM*((iedge+1)%nbEdges);
 129           double a=mon_determinant(ptToTest, A, B);
 130           if(a<-eps)
 131             sign[iedge]=-1;
 132           else if(a>eps)
 133             sign[iedge]=1;
 134           else
 135             sign[iedge]=0;
 136         }
 137       bool ret=decideFromSign(sign.get(), nbEdges);
 138       return ret;
 139     }
 140
 141     // Same as isElementContainsPointAlg2DSimple() with a different input format ...
 142     static bool isElementContainsPointAlgo2DSimple2(const double *ptToTest, NormalizedCellType type,
 143                                                     const double *coords, const typename MyMeshType::MyConnType *conn_elem,
 144                                                     typename MyMeshType::MyConnType conn_elem_sz, double eps)
 145     {
 146       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 147       typedef typename MyMeshType::MyConnType ConnType;
 148       const NumberingPolicy numPol=MyMeshType::My_numPol;
 149
 150       const CellModel& cmType=CellModel::GetCellModel(type);
 151       bool ret=false;
 152
 153       int nbEdges=cmType.getNumberOfSons();
 154       double *pts = new double[nbEdges*SPACEDIM];
 155       for (int iedge=0; iedge<nbEdges; iedge++)
 156         {
 157           const double* a=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
 158           std::copy(a,a+SPACEDIM,pts+iedge*SPACEDIM);
 159         }
 160       ret=isElementContainsPointAlg2DSimple(ptToTest,pts,nbEdges,eps);
 161       delete [] pts;
 162       return ret;
 163     }
 164
 165     static bool isElementContainsPointAlgo2DPolygon(const double *ptToTest, NormalizedCellType type,
 166                                                     const double *coords, const typename MyMeshType::MyConnType *conn_elem,
 167                                                     typename MyMeshType::MyConnType conn_elem_sz, double eps)
 168     {
 169       // Override precision for this method only:
 170       INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
 171
 172       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 173       typedef typename MyMeshType::MyConnType ConnType;
 174       const NumberingPolicy numPol=MyMeshType::My_numPol;
 175
 176       std::vector<INTERP_KERNEL::Node *> nodes(conn_elem_sz);
 177       INTERP_KERNEL::QuadraticPolygon *pol(0);
 178       for(mcIdType j=0;j<conn_elem_sz;j++)
 179         {
 180           mcIdType nodeId(OTT<ConnType,numPol>::ind2C(conn_elem[j]));
 181           nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*SPACEDIM],coords[nodeId*SPACEDIM+1]);
 182         }
 183       if(!INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic())
 184         pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
 185       else
 186         pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
 187       INTERP_KERNEL::Node *n(new INTERP_KERNEL::Node(ptToTest[0],ptToTest[1]));
 188       double a(0.),b(0.),c(0.);
 189       a=pol->normalizeMe(b,c); n->applySimilarity(b,c,a);
 190       bool ret=pol->isInOrOut2(n);
 191       delete pol; n->decrRef();
 192       return ret;
 193     }
 194
 195     static bool isElementContainsPointAlg3D(const double *ptToTest, const typename MyMeshType::MyConnType *conn_elem, typename MyMeshType::MyConnType conn_elem_sz, const double *coords, const CellModel& cmType, double eps)
 196     {
 197       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 198       typedef typename MyMeshType::MyConnType ConnType;
 199       const NumberingPolicy numPol=MyMeshType::My_numPol;
 200
 201       int nbfaces = cmType.getNumberOfSons2(conn_elem,conn_elem_sz);
 202       std::unique_ptr<char[]> sign( new char[nbfaces] );
 203       std::unique_ptr<ConnType[]> connOfSon( new ConnType[conn_elem_sz] );
 204       std::unique_ptr<double[]> ptsOfTetrahedrizedPolyhedron( new double[3*3*nbfaces] );
 205       for (int iface=0; iface<nbfaces; iface++)
 206         {
 207           NormalizedCellType typeOfSon;
 208           unsigned nb_face_nodes = cmType.fillSonCellNodalConnectivity2(iface,conn_elem,conn_elem_sz,connOfSon.get(),typeOfSon);
 209           // in case of polyhedron, points can be colinear on some faces
 210           // => loop until non colinear points are found
 211           // or element type is not polyhedron
 212           // or max number of nodes is reached to avoid infinite loop on this face
 213           bool non_colinear_points_found = false;
 214           unsigned i_permutation = 0;
 215           while (!non_colinear_points_found)
 216           {
 217             const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[(0+i_permutation) % nb_face_nodes]));
 218             const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[(1+i_permutation) % nb_face_nodes]));
 219             const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[(2+i_permutation) % nb_face_nodes]));
 220             if (cmType.getEnum() == NORM_POLYHED && isColinear3DPts(AA, BB, CC, eps) && i_permutation < nb_face_nodes)
 221               // points are colinear => use next 3 other points
 222               i_permutation+=1;
 223             else
 224               {
 225               std::copy(AA,AA+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface);
 226               std::copy(BB,BB+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface+3);
 227               std::copy(CC,CC+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface+6);
 228               non_colinear_points_found = true;
 229             }
 230           }
 231         }
 232       double centerPt[3];
 233       double normalizeFact = NormalizeTetrahedrizedPolyhedron(ptsOfTetrahedrizedPolyhedron.get(),nbfaces,centerPt);
 234       double ptToTestNorm[3] = {(ptToTest[0]-centerPt[0])/normalizeFact,
 235       (ptToTest[1]-centerPt[1])/normalizeFact,
 236       (ptToTest[2]-centerPt[2])/normalizeFact};
 237       for (int iface=0; iface<nbfaces; iface++)
 238         {
 239           double lengthNorm = TripleProduct(ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 3,
 240                                             ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 6,
 241                                             ptToTestNorm,
 242                                             ptsOfTetrahedrizedPolyhedron.get() + 9*iface);
 243           // assigne sign[iface] : 0 means on the face within eps, 1 or -1 means far from tetra representing face
 244           if( lengthNorm<-eps )
 245             sign[iface]=-1;
 246           else if( lengthNorm>eps )
 247             sign[iface]=1;
 248           else
 249             sign[iface]=0;
 250         }
 251       bool ret=decideFromSign(sign.get(), nbfaces);
 252       return ret;
 253     }
 254
 255     /*!
 256      * Precondition : spacedim==meshdim. To be checked upstream to this call.
 257      */
 258     static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords,
 259                                        const typename MyMeshType::MyConnType *conn_elem,
 260                                        typename MyMeshType::MyConnType conn_elem_sz, double eps)
 261     {
 262       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 263       typedef typename MyMeshType::MyConnType ConnType;
 264       const NumberingPolicy numPol=MyMeshType::My_numPol;
 265
 266       const CellModel& cmType=CellModel::GetCellModel(type);
 267       //
 268       if (SPACEDIM==2)
 269         {
 270           if(type != INTERP_KERNEL::NORM_POLYGON && !cmType.isQuadratic())
 271             return isElementContainsPointAlgo2DSimple2(ptToTest, type, coords, conn_elem, conn_elem_sz, eps);
 272           else
 273             return isElementContainsPointAlgo2DPolygon(ptToTest, type, coords, conn_elem, conn_elem_sz, eps);
 274         }
 275       if (SPACEDIM==3)
 276         {
 277           return isElementContainsPointAlg3D(ptToTest,conn_elem,conn_elem_sz,coords,cmType,eps);
 278         }
 279
 280       if(SPACEDIM==1)
 281         {
 282           double p1=coords[(OTT<ConnType,numPol>::ind2C(conn_elem[0]))];
 283           double p2=coords[(OTT<ConnType,numPol>::ind2C(conn_elem[1]))];
 284           double delta=fabs(p1-p2)+eps;
 285           double val=*ptToTest-std::min(p1,p2);
 286           return val>-eps && val<delta;
 287         }
 288       throw INTERP_KERNEL::Exception("Invalid spacedim detected ! Managed spaceDim are 2 and 3 !");
 289     }
 290
 291     bool elementContainsPoint(typename MyMeshType::MyConnType i, const double* x, double eps)
 292     {
 293       //as i is extracted from the BBTRee, it is already in C numbering
 294       //it is not necessary to convert it from F to C
 295       typedef typename MyMeshType::MyConnType ConnType;
 296       const NumberingPolicy numPol=MyMeshType::My_numPol;
 297
 298       const double* coords= _mesh.getCoordinatesPtr();
 299       const ConnType* conn=_mesh.getConnectivityPtr();
 300       const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
 301       const ConnType* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
 302       int conn_elem_sz=conn_index[i+1]-conn_index[i];
 303       NormalizedCellType type=_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
 304       return isElementContainsPoint(x,type,coords,conn_elem,conn_elem_sz,eps);
 305     }
 306
 307     static bool decideFromSign(const char *sign, mcIdType nbelem)
 308     {
 309       char min_sign = 1;
 310       char max_sign = -1;
 311       for (int i=0; i<nbelem;i++)
 312         {
 313           min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
 314           max_sign=(sign[i]>max_sign)?sign[i]:max_sign;
 315         }
 316       return (min_sign!=-1 || max_sign!=1);
 317     }
 318   };
 319
 320   template<class MyMeshType>
 321   class PointLocatorInSimplex : public PointLocatorAlgos<MyMeshType>
 322   {
 323     const MyMeshType& _mesh;
 324   public:
 325     PointLocatorInSimplex(const MyMeshType& mesh)
 326       :PointLocatorAlgos<MyMeshType>(mesh),_mesh(mesh)
 327     {
 328     }
 329
 330     //================================================================================
 331     /*!
 332      * \brief Returns nodes composing the simplex the point x is in
 333      */
 334     //================================================================================
 335
 336     virtual std::list<typename MyMeshType::MyConnType> locates(const double* x, double eps)
 337     {
 338       typedef typename MyMeshType::MyConnType ConnType;
 339       const NumberingPolicy numPol=MyMeshType::My_numPol;
 340
 341       std::list<ConnType> simplexNodes;
 342       std::list<ConnType> candidates = PointLocatorAlgos<MyMeshType>::locates(x,eps);
 343       typename std::list<ConnType>::iterator eIt = candidates.begin();
 344       for ( ; eIt != candidates.end(); ++eIt )
 345         {
 346           const ConnType i = OTT<ConnType,numPol>::ind2C( *eIt );
 347           const double* coords= _mesh.getCoordinatesPtr();
 348           const ConnType* conn=_mesh.getConnectivityPtr();
 349           const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
 350           const ConnType* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
 351           int conn_elem_sz=conn_index[i+1]-conn_index[i];
 352           NormalizedCellType type=_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
 353           CellModel cell = CellModel::GetCellModel(type);
 354
 355           if ( cell.isQuadratic() )
 356             throw Exception("P2 not implemented yet");
 357
 358           if ( cell.isSimplex())
 359             {
 360               for ( int n = 0; n < conn_elem_sz; ++n )
 361                 simplexNodes.push_back( conn_elem[ n ]);
 362             }
 363           else
 364             {
 365               NormalizedCellType simlexType = cell.getDimension()==3 ? NORM_TETRA4 : NORM_TRI3;
 366               std::vector<mcIdType> sonNodes;
 367               NormalizedCellType sonType;
 368               const unsigned nbSons = cell.getNumberOfSons2( conn_elem, conn_elem_sz );
 369               for ( unsigned s = 0; s < nbSons; ++s )
 370                 {
 371                   sonNodes.resize( cell.getNumberOfNodesConstituentTheSon2( s, conn_elem, conn_elem_sz ));
 372                   cell.fillSonCellNodalConnectivity2( s, conn_elem, conn_elem_sz, &sonNodes[0], sonType );
 373                   std::set<mcIdType> sonNodesSet( sonNodes.begin(), sonNodes.end() );
 374
 375                   std::set< std::set< ConnType > > checkedSonSimplex;
 376                   for ( unsigned sn = 0; sn < sonNodes.size(); ++sn )
 377                     {
 378                       std::vector< ConnType > simplexConn( cell.getDimension() + 1 );
 379                       unsigned n;
 380                       for ( n = 0; n < cell.getDimension()-1; ++n )
 381                         simplexConn[n] = sonNodes[ (sn+n) % sonNodes.size() ];
 382
 383                       for ( unsigned n2 = 0; n2 < sonNodes.size()-cell.getDimension()+1; ++n2 )
 384                         {
 385                           simplexConn[n] = sonNodes[ (sn+n+n2) % sonNodes.size() ];
 386                           std::set< ConnType > sonSimplex( simplexConn.begin(), --simplexConn.end());
 387                           if ( checkedSonSimplex.insert( sonSimplex ).second )
 388                             {
 389                               for ( unsigned cn = 0; cn < conn_elem_sz; ++cn )
 390                                 if ( !sonNodesSet.count( conn_elem[cn] ))
 391                                   {
 392                                     simplexConn.back() = conn_elem[cn];
 393                                     if ( this->isElementContainsPoint( x, simlexType, coords,
 394                                                                        &simplexConn[0], simplexConn.size(), eps ))
 395                                       {
 396                                         simplexNodes.insert( simplexNodes.end(),
 397                                                              simplexConn.begin(), simplexConn.end());
 398                                         return simplexNodes;
 399                                       }
 400                                   }
 401                             }
 402                         }
 403                     }
 404                 }
 405             }
 406         }
 407       return simplexNodes;
 408     }
 409
 410   };
 411 }
 412
 413 #endif