1 // Copyright (C) 2007-2024 CEA, EDF
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.
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"
26 #include "VectorUtils.hxx"
28 #include "InterpKernelGeo2DNode.hxx"
29 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
37 namespace INTERP_KERNEL
39 class GenericPointLocatorAlgos
42 virtual ~GenericPointLocatorAlgos() { }
43 virtual std::list<mcIdType> locates(const double* x, double eps) = 0;
46 template<class MyMeshType>
47 class PointLocatorAlgos: public GenericPointLocatorAlgos
51 BBTree<MyMeshType::MY_SPACEDIM,typename MyMeshType::MyConnType>* _tree;
52 const MyMeshType& _mesh;
54 PointLocatorAlgos(const MyMeshType& mesh):_mesh(mesh)
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++)
66 for (int idim=0; idim<SPACEDIM; idim++)
68 _bb[2*(i*SPACEDIM+idim)]=std::numeric_limits<double>::max();
69 _bb[2*(i*SPACEDIM+idim)+1]=-std::numeric_limits<double>::max();
71 for (ConnType index= conn_index[i]; index < conn_index[i+1];index++)
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;
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++)
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];
84 _tree=new BBTree<SPACEDIM,typename MyMeshType::MyConnType>(_bb,0,0,nelem);
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)
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++)
104 ConnType ielem=candidates[i];
105 if (elementContainsPoint(ielem,x,eps))
106 retlist.push_back(OTT<ConnType,numPol>::indFC(ielem));
111 static bool isElementContainsPointAlg2DSimple(const double *ptToTest, const double *cellPts, mcIdType nbEdges, double eps)
113 /* with dimension 2, it suffices to check all the edges
114 and see if the sign of double products from the point
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++)
127 const double* A=cellPts+SPACEDIM*iedge;
128 const double* B=cellPts+SPACEDIM*((iedge+1)%nbEdges);
129 double a=mon_determinant(ptToTest, A, B);
137 bool ret=decideFromSign(sign.get(), nbEdges);
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)
146 const int SPACEDIM=MyMeshType::MY_SPACEDIM;
147 typedef typename MyMeshType::MyConnType ConnType;
148 const NumberingPolicy numPol=MyMeshType::My_numPol;
150 const CellModel& cmType=CellModel::GetCellModel(type);
153 int nbEdges=cmType.getNumberOfSons();
154 double *pts = new double[nbEdges*SPACEDIM];
155 for (int iedge=0; iedge<nbEdges; iedge++)
157 const double* a=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
158 std::copy(a,a+SPACEDIM,pts+iedge*SPACEDIM);
160 ret=isElementContainsPointAlg2DSimple(ptToTest,pts,nbEdges,eps);
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)
169 // Override precision for this method only:
170 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
172 const int SPACEDIM=MyMeshType::MY_SPACEDIM;
173 typedef typename MyMeshType::MyConnType ConnType;
174 const NumberingPolicy numPol=MyMeshType::My_numPol;
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++)
180 mcIdType nodeId(OTT<ConnType,numPol>::ind2C(conn_elem[j]));
181 nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*SPACEDIM],coords[nodeId*SPACEDIM+1]);
183 if(!INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic())
184 pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
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();
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)
197 const int SPACEDIM=MyMeshType::MY_SPACEDIM;
198 typedef typename MyMeshType::MyConnType ConnType;
199 const NumberingPolicy numPol=MyMeshType::My_numPol;
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++)
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)
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
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;
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++)
239 double lengthNorm = TripleProduct(ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 3,
240 ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 6,
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 )
246 else if( lengthNorm>eps )
251 bool ret=decideFromSign(sign.get(), nbfaces);
256 * Precondition : spacedim==meshdim. To be checked upstream to this call.
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)
262 const int SPACEDIM=MyMeshType::MY_SPACEDIM;
263 typedef typename MyMeshType::MyConnType ConnType;
264 const NumberingPolicy numPol=MyMeshType::My_numPol;
266 const CellModel& cmType=CellModel::GetCellModel(type);
270 if(type != INTERP_KERNEL::NORM_POLYGON && !cmType.isQuadratic())
271 return isElementContainsPointAlgo2DSimple2(ptToTest, type, coords, conn_elem, conn_elem_sz, eps);
273 return isElementContainsPointAlgo2DPolygon(ptToTest, type, coords, conn_elem, conn_elem_sz, eps);
277 return isElementContainsPointAlg3D(ptToTest,conn_elem,conn_elem_sz,coords,cmType,eps);
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;
288 throw INTERP_KERNEL::Exception("Invalid spacedim detected ! Managed spaceDim are 2 and 3 !");
291 bool elementContainsPoint(typename MyMeshType::MyConnType i, const double* x, double eps)
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;
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);
307 static bool decideFromSign(const char *sign, mcIdType nbelem)
311 for (int i=0; i<nbelem;i++)
313 min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
314 max_sign=(sign[i]>max_sign)?sign[i]:max_sign;
316 return (min_sign!=-1 || max_sign!=1);
320 template<class MyMeshType>
321 class PointLocatorInSimplex : public PointLocatorAlgos<MyMeshType>
323 const MyMeshType& _mesh;
325 PointLocatorInSimplex(const MyMeshType& mesh)
326 :PointLocatorAlgos<MyMeshType>(mesh),_mesh(mesh)
330 //================================================================================
332 * \brief Returns nodes composing the simplex the point x is in
334 //================================================================================
336 virtual std::list<typename MyMeshType::MyConnType> locates(const double* x, double eps)
338 typedef typename MyMeshType::MyConnType ConnType;
339 const NumberingPolicy numPol=MyMeshType::My_numPol;
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 )
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);
355 if ( cell.isQuadratic() )
356 throw Exception("P2 not implemented yet");
358 if ( cell.isSimplex())
360 for ( int n = 0; n < conn_elem_sz; ++n )
361 simplexNodes.push_back( conn_elem[ n ]);
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 )
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() );
375 std::set< std::set< ConnType > > checkedSonSimplex;
376 for ( unsigned sn = 0; sn < sonNodes.size(); ++sn )
378 std::vector< ConnType > simplexConn( cell.getDimension() + 1 );
380 for ( n = 0; n < cell.getDimension()-1; ++n )
381 simplexConn[n] = sonNodes[ (sn+n) % sonNodes.size() ];
383 for ( unsigned n2 = 0; n2 < sonNodes.size()-cell.getDimension()+1; ++n2 )
385 simplexConn[n] = sonNodes[ (sn+n+n2) % sonNodes.size() ];
386 std::set< ConnType > sonSimplex( simplexConn.begin(), --simplexConn.end());
387 if ( checkedSonSimplex.insert( sonSimplex ).second )
389 for ( unsigned cn = 0; cn < conn_elem_sz; ++cn )
390 if ( !sonNodesSet.count( conn_elem[cn] ))
392 simplexConn.back() = conn_elem[cn];
393 if ( this->isElementContainsPoint( x, simlexType, coords,
394 &simplexConn[0], simplexConn.size(), eps ))
396 simplexNodes.insert( simplexNodes.end(),
397 simplexConn.begin(), simplexConn.end());