src/MEDMEM/PointLocatorAlgos.txx

   1 //  Copyright (C) 2007-2008  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.
   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 #ifndef _POINT_LOCATOR_ALGOS_TXX_
  20 #define _POINT_LOCATOR_ALGOS_TXX_
  21
  22 #include "InterpolationUtils.hxx"
  23 #include "CellModel.hxx"
  24 #include "BBTree.txx"
  25 #include <list>
  26 #include <limits>
  27
  28 namespace INTERP_KERNEL{
  29
  30   class GenericPointLocatorAlgos
  31   {
  32   public:
  33     virtual ~GenericPointLocatorAlgos(){};
  34     virtual std::list<int> locates(const double* x)=0;
  35
  36   };
  37
  38   template<class MyMeshType>
  39   class PointLocatorAlgos: public GenericPointLocatorAlgos
  40   {
  41   private :
  42     double* _bb;
  43     BBTree<MyMeshType::MY_SPACEDIM,typename MyMeshType::MyConnType>* _tree;
  44     const MyMeshType& _mesh;
  45
  46   public:
  47     PointLocatorAlgos(const MyMeshType& mesh):_mesh(mesh)
  48     {
  49       typedef typename MyMeshType::MyConnType ConnType;
  50       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
  51       const NumberingPolicy numPol=MyMeshType::My_numPol;
  52       int nelem = _mesh.getNumberOfElements();
  53       _bb = new double[SPACEDIM*2*nelem];
  54       const ConnType* conn = _mesh.getConnectivityPtr();
  55       const ConnType* conn_index = _mesh.getConnectivityIndexPtr();
  56       const double* coords=_mesh.getCoordinatesPtr();
  57       for (int i=0; i<nelem; i++)
  58         {
  59           for (int idim=0; idim<SPACEDIM; idim++)
  60             {
  61               _bb[2*(i*SPACEDIM+idim)]=std::numeric_limits<double>::max();
  62               _bb[2*(i*SPACEDIM+idim)+1]=-std::numeric_limits<double>::max();
  63             }
  64           for (int index= conn_index[i]; index < conn_index[i+1];index++)
  65             {
  66               //coordelem points to the coordinates of the current node of the i-th element
  67               const double* coordelem = coords+OTT<ConnType,numPol>::ind2C(conn[OTT<ConnType,numPol>::ind2C(index)])*SPACEDIM;
  68
  69               //the bounding box is updated by checking wheher the node is at the min/max in exach dimension
  70               for (int idim=0; idim<SPACEDIM;idim++)
  71                 {
  72                   _bb[2*(i*SPACEDIM+idim)]=(coordelem[idim]<_bb[2*(i*SPACEDIM+idim)])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)];
  73                   _bb[2*(i*SPACEDIM+idim)+1]=(coordelem[idim]>_bb[2*(i*SPACEDIM+idim)+1])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)+1];
  74                 }
  75             }
  76         }
  77
  78       _tree=new BBTree<SPACEDIM,typename MyMeshType::MyConnType>(_bb,0,0,nelem);
  79     }
  80     ~PointLocatorAlgos()
  81     {
  82       delete[] _bb;
  83       delete _tree;
  84     }
  85
  86     //returns the list of elements that contains
  87     //the point pointed to by x
  88     std::list<typename MyMeshType::MyConnType> locates(const double* x)
  89     {
  90       typedef typename MyMeshType::MyConnType ConnType;
  91       const NumberingPolicy numPol=MyMeshType::My_numPol;
  92       std::vector<ConnType> candidates;
  93       _tree->getElementsAroundPoint(x,candidates);
  94       std::list<ConnType> retlist;
  95       for(int i=0; i< candidates.size(); i++)
  96         {
  97           int ielem=candidates[i];
  98           if (elementContainsPoint(ielem,x))
  99             retlist.push_back(OTT<ConnType,numPol>::indFC(ielem));
 100         }
 101       return retlist;
 102     }
 103
 104     bool elementContainsPoint(typename MyMeshType::MyConnType i, const double*x)
 105     {
 106       //as i is extracted from the BBTRee, it is already in C numbering
 107       //it is not necessary to convert it from F to C
 108       typedef typename MyMeshType::MyConnType ConnType;
 109       const int SPACEDIM=MyMeshType::MY_SPACEDIM;
 110       const NumberingPolicy numPol=MyMeshType::My_numPol;
 111
 112       const double* coords= _mesh.getCoordinatesPtr();
 113       const ConnType* conn=_mesh.getConnectivityPtr();
 114       const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
 115       const ConnType* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
 116       NormalizedCellType type=_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
 117       const CellModel& cmType=CellModel::getCellModel(type);
 118
 119       int nbnodes = cmType.getNumberOfNodes();//conn_index[i+1]-conn_index[i];
 120
 121       // with dimension 2, it suffices to check all the edges
 122       // and see if the sign of double products from the point
 123       //is always the same.
 124       //                 C
 125       //                / \
 126       //               /   \
 127       //     Xo       /     \
 128       //             A-------B
 129       //
 130       //here XA^XC and XC^XB have different signs
 131       //
 132       if (SPACEDIM==2)
 133         {
 134           //in 2D, nbedges==nbnodes
 135           int nbedges=nbnodes;
 136           int* sign = new int[nbedges];
 137           for (int iedge=0; iedge<nbedges; iedge++)
 138             {
 139               const double* A=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
 140               const double* B;
 141               if (iedge+1< nbedges)
 142                 B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge+1]));
 143               else
 144                 B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[0]));
 145
 146               double a=mon_determinant(x, A, B);
 147               if (a<-1e-12)
 148                 sign[iedge]=-1;
 149               else if (a>1e-12)
 150                 sign[iedge]=1;
 151               else
 152                 sign[iedge]=0;
 153             }
 154           bool ret=decide_from_sign(sign, nbedges);
 155           delete [] sign;
 156           return ret;
 157         }
 158
 159       if (SPACEDIM==3)
 160         {
 161           int nbfaces = cmType.getNumberOfSons();
 162           int* sign = new int[nbfaces];
 163           for (int iface=0; iface<nbfaces; iface++)
 164             {
 165               const unsigned* connface=cmType.getNodesConstituentTheSon(iface);
 166               const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[0]]));
 167               const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[1]]));
 168               const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[2]]));
 169
 170               double Vol=triple_product(AA,BB,CC,x);
 171               if (Vol<-1e-12)
 172                 sign[iface]=-1;
 173               else if (Vol>1e-12)
 174                 sign[iface]=1;
 175               else
 176                 sign[iface]=0;
 177             }
 178           bool ret=decide_from_sign(sign, nbfaces);
 179           delete [] sign;
 180           return ret;
 181         }
 182
 183     }
 184
 185     bool decide_from_sign (const int* sign, int nbelem)
 186     {
 187       int min_sign =1;
 188       int max_sign =-1;
 189       for (int i=0; i<nbelem;i++)
 190         {
 191           min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
 192           max_sign=(sign[i]>max_sign)?sign[i]:max_sign;
 193         }
 194       return (min_sign!=-1 || max_sign!=1);
 195     }
 196   };
 197 }
 198 #endif