1 // Copyright (C) 2007-2008 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 #ifndef _POINT_LOCATOR_ALGOS_TXX_
20 #define _POINT_LOCATOR_ALGOS_TXX_
22 #include "InterpolationUtils.hxx"
23 #include "CellModel.hxx"
28 namespace INTERP_KERNEL{
30 class GenericPointLocatorAlgos
33 virtual ~GenericPointLocatorAlgos(){};
34 virtual std::list<int> locates(const double* x)=0;
38 template<class MyMeshType>
39 class PointLocatorAlgos: public GenericPointLocatorAlgos
43 BBTree<MyMeshType::MY_SPACEDIM,typename MyMeshType::MyConnType>* _tree;
44 const MyMeshType& _mesh;
47 PointLocatorAlgos(const MyMeshType& mesh):_mesh(mesh)
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++)
59 for (int idim=0; idim<SPACEDIM; idim++)
61 _bb[2*(i*SPACEDIM+idim)]=std::numeric_limits<double>::max();
62 _bb[2*(i*SPACEDIM+idim)+1]=-std::numeric_limits<double>::max();
64 for (int index= conn_index[i]; index < conn_index[i+1];index++)
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;
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++)
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];
78 _tree=new BBTree<SPACEDIM,typename MyMeshType::MyConnType>(_bb,0,0,nelem);
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)
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++)
97 int ielem=candidates[i];
98 if (elementContainsPoint(ielem,x))
99 retlist.push_back(OTT<ConnType,numPol>::indFC(ielem));
104 bool elementContainsPoint(typename MyMeshType::MyConnType i, const double*x)
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;
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);
119 int nbnodes = cmType.getNumberOfNodes();//conn_index[i+1]-conn_index[i];
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.
130 //here XA^XC and XC^XB have different signs
134 //in 2D, nbedges==nbnodes
136 int* sign = new int[nbedges];
137 for (int iedge=0; iedge<nbedges; iedge++)
139 const double* A=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
141 if (iedge+1< nbedges)
142 B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge+1]));
144 B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[0]));
146 double a=mon_determinant(x, A, B);
154 bool ret=decide_from_sign(sign, nbedges);
161 int nbfaces = cmType.getNumberOfSons();
162 int* sign = new int[nbfaces];
163 for (int iface=0; iface<nbfaces; iface++)
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]]));
170 double Vol=triple_product(AA,BB,CC,x);
178 bool ret=decide_from_sign(sign, nbfaces);
185 bool decide_from_sign (const int* sign, int nbelem)
189 for (int i=0; i<nbelem;i++)
191 min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
192 max_sign=(sign[i]>max_sign)?sign[i]:max_sign;
194 return (min_sign!=-1 || max_sign!=1);