src/INTERP_KERNEL/PointLocatorAlgos.txx

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