1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
20 #ifndef __POINTLOCATORALGOS_TXX__
21 #define __POINTLOCATORALGOS_TXX__
23 #include "InterpolationUtils.hxx"
24 #include "CellModel.hxx"
31 namespace INTERP_KERNEL
33 class GenericPointLocatorAlgos
36 virtual ~GenericPointLocatorAlgos() { }
37 virtual std::list<int> locates(const double* x, double eps) = 0;
40 template<class MyMeshType>
41 class PointLocatorAlgos: public GenericPointLocatorAlgos
45 BBTree<MyMeshType::MY_SPACEDIM,typename MyMeshType::MyConnType>* _tree;
46 const MyMeshType& _mesh;
48 PointLocatorAlgos(const MyMeshType& mesh):_mesh(mesh)
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++)
60 for (int idim=0; idim<SPACEDIM; idim++)
62 _bb[2*(i*SPACEDIM+idim)]=std::numeric_limits<double>::max();
63 _bb[2*(i*SPACEDIM+idim)+1]=-std::numeric_limits<double>::max();
65 for (int index= conn_index[i]; index < conn_index[i+1];index++)
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;
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++)
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];
78 _tree=new BBTree<SPACEDIM,typename MyMeshType::MyConnType>(_bb,0,0,nelem);
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)
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++)
98 int ielem=candidates[i];
99 if (elementContainsPoint(ielem,x,eps))
100 retlist.push_back(OTT<ConnType,numPol>::indFC(ielem));
105 static bool isElementContainsPointAlg2D(const double *ptToTest, const double *cellPts, int nbEdges, double eps)
107 /* with dimension 2, it suffices to check all the edges
108 and see if the sign of double products from the point
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++)
121 const double* A=cellPts+SPACEDIM*iedge;
122 const double* B=cellPts+SPACEDIM*((iedge+1)%nbEdges);
123 double a=mon_determinant(ptToTest, A, B);
131 bool ret=decideFromSign(sign, nbEdges);
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)
138 const int SPACEDIM=MyMeshType::MY_SPACEDIM;
139 typedef typename MyMeshType::MyConnType ConnType;
140 const NumberingPolicy numPol=MyMeshType::My_numPol;
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++)
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);
160 bool ret=decideFromSign(sign, nbfaces);
166 static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords, const typename MyMeshType::MyConnType *conn_elem, int conn_elem_sz, double eps)
168 const int SPACEDIM=MyMeshType::MY_SPACEDIM;
169 typedef typename MyMeshType::MyConnType ConnType;
170 const NumberingPolicy numPol=MyMeshType::My_numPol;
172 const CellModel& cmType=CellModel::GetCellModel(type);
176 int nbEdges=cmType.getNumberOfSons();
177 double *pts = new double[nbEdges*SPACEDIM];
178 for (int iedge=0; iedge<nbEdges; iedge++)
180 const double* a=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
181 std::copy(a,a+SPACEDIM,pts+iedge*SPACEDIM);
183 bool ret=isElementContainsPointAlg2D(ptToTest,pts,nbEdges,eps);
190 return isElementContainsPointAlg3D(ptToTest,conn_elem,conn_elem_sz,coords,cmType,eps);
195 double p1=coords[(OTT<ConnType,numPol>::ind2C(conn_elem[0]))];
196 double p2=coords[(OTT<ConnType,numPol>::ind2C(conn_elem[1]))];
197 double delta=fabs(p1-p2)+eps;
198 double val=*ptToTest-std::min(p1,p2);
199 return val>-eps && val<delta;
201 throw INTERP_KERNEL::Exception("Invalid spacedim detected ! Managed spaceDim are 2 and 3 !");
204 bool elementContainsPoint(typename MyMeshType::MyConnType i, const double* x, double eps)
206 //as i is extracted from the BBTRee, it is already in C numbering
207 //it is not necessary to convert it from F to C
208 typedef typename MyMeshType::MyConnType ConnType;
209 const NumberingPolicy numPol=MyMeshType::My_numPol;
211 const double* coords= _mesh.getCoordinatesPtr();
212 const ConnType* conn=_mesh.getConnectivityPtr();
213 const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
214 const ConnType* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
215 int conn_elem_sz=conn_index[i+1]-conn_index[i];
216 NormalizedCellType type=_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
217 return isElementContainsPoint(x,type,coords,conn_elem,conn_elem_sz,eps);
220 static bool decideFromSign(const int* sign, int nbelem)
224 for (int i=0; i<nbelem;i++)
226 min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
227 max_sign=(sign[i]>max_sign)?sign[i]:max_sign;
229 return (min_sign!=-1 || max_sign!=1);
233 template<class MyMeshType>
234 class PointLocatorInSimplex : public PointLocatorAlgos<MyMeshType>
236 const MyMeshType& _mesh;
238 PointLocatorInSimplex(const MyMeshType& mesh)
239 :PointLocatorAlgos<MyMeshType>(mesh),_mesh(mesh)
243 //================================================================================
245 * \brief Returns nodes composing the simplex the point x is in
247 //================================================================================
249 virtual std::list<int> locates(const double* x, double eps)
251 typedef typename MyMeshType::MyConnType ConnType;
252 const NumberingPolicy numPol=MyMeshType::My_numPol;
254 std::list<int> simplexNodes;
255 std::list<int> candidates = PointLocatorAlgos<MyMeshType>::locates(x,eps);
256 std::list<int>::iterator eIt = candidates.begin();
257 for ( ; eIt != candidates.end(); ++eIt )
259 const int i = OTT<ConnType,numPol>::ind2C( *eIt );
260 const double* coords= _mesh.getCoordinatesPtr();
261 const ConnType* conn=_mesh.getConnectivityPtr();
262 const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
263 const ConnType* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
264 int conn_elem_sz=conn_index[i+1]-conn_index[i];
265 NormalizedCellType type=_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
266 CellModel cell = CellModel::GetCellModel(type);
268 if ( cell.isQuadratic() )
269 throw Exception("P2 not implemented yet");
271 if ( cell.isSimplex())
273 for ( int n = 0; n < conn_elem_sz; ++n )
274 simplexNodes.push_back( conn_elem[ n ]);
278 NormalizedCellType simlexType = cell.getDimension()==3 ? NORM_TETRA4 : NORM_TRI3;
279 std::vector<int> sonNodes;
280 NormalizedCellType sonType;
281 const unsigned nbSons = cell.getNumberOfSons2( conn_elem, conn_elem_sz );
282 for ( unsigned s = 0; s < nbSons; ++s )
284 sonNodes.resize( cell.getNumberOfNodesConstituentTheSon2( s, conn_elem, conn_elem_sz ));
285 cell.fillSonCellNodalConnectivity2( s, conn_elem, conn_elem_sz, &sonNodes[0], sonType );
286 std::set<int> sonNodesSet( sonNodes.begin(), sonNodes.end() );
288 std::set< std::set< ConnType > > checkedSonSimplex;
289 for ( unsigned sn = 0; sn < sonNodes.size(); ++sn )
291 std::vector< ConnType > simplexConn( cell.getDimension() + 1 );
293 for ( n = 0; n < cell.getDimension()-1; ++n )
294 simplexConn[n] = sonNodes[ (sn+n) % sonNodes.size() ];
296 for ( unsigned n2 = 0; n2 < sonNodes.size()-cell.getDimension()+1; ++n2 )
298 simplexConn[n] = sonNodes[ (sn+n+n2) % sonNodes.size() ];
299 std::set< ConnType > sonSimplex( simplexConn.begin(), --simplexConn.end());
300 if ( checkedSonSimplex.insert( sonSimplex ).second )
302 for ( unsigned cn = 0; cn < conn_elem_sz; ++cn )
303 if ( !sonNodesSet.count( conn_elem[cn] ))
305 simplexConn.back() = conn_elem[cn];
306 if ( this->isElementContainsPoint( x, simlexType, coords,
307 &simplexConn[0], simplexConn.size(), eps ))
309 simplexNodes.insert( simplexNodes.end(),
310 simplexConn.begin(), simplexConn.end());