src/INTERP_KERNEL/PointLocatorAlgos.txx

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