From: vbd <vbd>
Date: Tue, 11 Mar 2008 15:02:02 +0000 (+0000)
Subject: adding PointLocator algorithms to INTERP_KERNEL
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=5cb0f0a13bd974f9b6018b795216b81b71bb45a3;p=tools%2Fmedcoupling.git

adding PointLocator algorithms to INTERP_KERNEL
---

diff --git a/src/INTERP_KERNEL/BBTree.txx b/src/INTERP_KERNEL/BBTree.txx
index dc2f10e56..597a2f08e 100644
--- a/src/INTERP_KERNEL/BBTree.txx
+++ b/src/INTERP_KERNEL/BBTree.txx
@@ -86,8 +86,8 @@ public:
 
 
       }
-    _max_left=max_left;
-    _min_right=min_right;
+    _max_left=max_left+1e-12;
+    _min_right=min_right-1e-12;
     _left=new BBTree(bbs, &(new_elems_left[0]), level+1, new_elems_left.size());
     _right=new BBTree(bbs, &(new_elems_right[0]), level+1, new_elems_right.size());
   
@@ -149,6 +149,49 @@ public:
     _right->getIntersectingElems(bb,elems);
   }
 
+ 
+  /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
+    \param xx pointer to query point coords
+    \param elems list of elements (given in 0-indexing) intersecting the bounding box
+  */
+
+  void getElementsAroundPoint(const double* xx, std::vector<int>& elems) const
+  {
+    //  terminal node : return list of elements intersecting bb
+    if (_terminal)
+      {
+        for (int i=0; i<_nbelems; i++)
+          {
+            const double* const  bb_ptr=_bb+_elems[i]*2*dim;
+            bool intersects = true;
+            for (int idim=0; idim<dim; idim++)
+              {
+                if (bb_ptr[idim*2]-xx[idim]>1e-12|| bb_ptr[idim*2+1]-xx[idim]<-1e-12)
+                  intersects=false;
+              }
+            if (intersects)
+              {
+                elems.push_back(_elems[i]);
+              }
+          }
+        return;
+      }
+
+    //non terminal node 
+    if (xx[_level%dim] < _min_right)
+      {
+        _left->getElementsAroundPoint(xx, elems);
+        return;
+      }
+    if (xx[_level%dim]> _max_left)
+      {
+        _right->getElementsAroundPoint(xx,elems);
+        return;
+      }
+    _left->getElementsAroundPoint(xx,elems);
+    _right->getElementsAroundPoint(xx,elems);
+  }
+
 
 
   int size()
@@ -157,5 +200,4 @@ public:
     return _left->size()+_right->size();
   }
 };
-
 #endif
diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx
index 350b143fe..df76196f1 100644
--- a/src/INTERP_KERNEL/Interpolation3D.txx
+++ b/src/INTERP_KERNEL/Interpolation3D.txx
@@ -46,16 +46,17 @@ namespace INTERP_KERNEL
    * intersection matrix. 
    * 
    * The matrix is partially sparse : it is a vector of maps of integer - double pairs. 
+	 * It can also be an INTERP_KERNEL::Matrix object.
    * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
    * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
    * which have a non-zero intersection volume with the target element. The vector has indices running from 
-   * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
+   * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
    * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
    * 
    
    * @param srcMesh     3-dimensional source mesh
    * @param targetMesh  3-dimesional target mesh, containing only tetraedra
-   * @param matrix      vector of maps in which the result is stored 
+   * @param result      matrix in which the result is stored 
    *
    */
   template<int SPACEDIM, int MESHDIM, class ConnType, NumberingPolicy numPol, class MatrixType, class MyMeshType>
diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx
index 7a79b962b..52d3e0615 100644
--- a/src/INTERP_KERNEL/InterpolationUtils.hxx
+++ b/src/INTERP_KERNEL/InterpolationUtils.hxx
@@ -572,4 +572,25 @@ namespace INTERP_KERNEL
 	}
 }
 
+/*! Computes the triple product (XA^XB).XC (in 3D)*/
+inline double triple_product(const double* A, const double*B, const double*C, const double*X)
+{
+	double XA[3];
+	XA[0]=A[0]-X[0];
+	XA[1]=A[1]-X[1];
+	XA[2]=A[2]-X[2];
+	double XB[3];
+	XB[0]=B[0]-X[0];
+	XB[1]=B[1]-X[1];
+	XB[2]=B[2]-X[2];
+	double XC[3];
+	XC[0]=C[0]-X[0];
+	XC[1]=C[1]-X[1];
+	XC[2]=C[2]-X[2];
+
+	return 
+		(XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
+		(XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
+		(XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
+}
 #endif
diff --git a/src/INTERP_KERNEL/IntersectorTetra.hxx b/src/INTERP_KERNEL/IntersectorTetra.hxx
index 57f90ff4a..4b4160342 100644
--- a/src/INTERP_KERNEL/IntersectorTetra.hxx
+++ b/src/INTERP_KERNEL/IntersectorTetra.hxx
@@ -12,6 +12,9 @@
 #include <map>
 #include <ext/hash_map>
 
+#include "MEDMEM_define.hxx"
+#include "MEDMEM_CellModel.hxx"
+
 using __gnu_cxx::hash_map;
 
 namespace INTERP_KERNEL
@@ -175,6 +178,7 @@ namespace INTERP_KERNEL
     inline void checkIsOutside(const double* pt, bool* isOutside) const;
     inline void calculateNode(ConnType globalNodeNum);
     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
+	
 
     /// disallow copying
     IntersectorTetra(const IntersectorTetra& t);
@@ -194,7 +198,8 @@ namespace INTERP_KERNEL
 
     /// reference to the source mesh
     const NormalizedUnstructuredMesh<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>& _srcMesh;
-    
+		
+		
   };
 
   /**
diff --git a/src/INTERP_KERNEL/IntersectorTetra.txx b/src/INTERP_KERNEL/IntersectorTetra.txx
index 867b1d050..c2ffd2db8 100644
--- a/src/INTERP_KERNEL/IntersectorTetra.txx
+++ b/src/INTERP_KERNEL/IntersectorTetra.txx
@@ -8,9 +8,6 @@
 #include "MeshUtils.hxx"
 #include "VectorUtils.hxx"
 #include "Log.hxx"
-// tony for the moment
-#include "MEDMEM_CellModel.hxx"
-#include "MEDMEM_define.hxx"
 
 #include <cmath>
 #include <cassert>
@@ -77,6 +74,7 @@ namespace INTERP_KERNEL
     delete _t;
     for(hash_map< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
       delete[] iter->second;
+	
   }
   
   /**
@@ -119,7 +117,7 @@ namespace INTERP_KERNEL
     ////const medGeometryElement type = _srcMesh.getNumberOfNodesOfElement(element);
     
     // get cell model for the element
-    const CELLMODEL cellModel(type);
+    const CELLMODEL& cellModel= CELLMODEL_Map::retrieveCellModel(type);
 
     // halfspace filtering
     bool isOutside[8] = {true, true, true, true, true, true, true, true};
@@ -155,7 +153,7 @@ namespace INTERP_KERNEL
         for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i)
           {
             const medGeometryElement faceType = cellModel.getConstituentType(1, i);
-            const CELLMODEL faceModel(faceType);
+            const CELLMODEL&  faceModel= CELLMODEL_Map::retrieveCellModel(faceType);
        
             assert(faceModel.getDimension() == 2);
 
diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am
index aee2f9b02..636df7b76 100644
--- a/src/INTERP_KERNEL/Makefile.am
+++ b/src/INTERP_KERNEL/Makefile.am
@@ -42,8 +42,8 @@ Interpolation.hxx\
 Interpolation3D.hxx\
 Interpolation2D.hxx\
 Interpolation3DSurf.hxx\
-Remapper.hxx
-
+Remapper.hxx\
+PointLocator.hxx
 
 # Libraries targets
 
@@ -53,7 +53,8 @@ TransformedTriangle_intersect.cxx\
 TransformedTriangle_math.cxx\
 BoundingBox.cxx\
 TetraAffineTransform.cxx\
-Remapper.cxx
+Remapper.cxx \
+PointLocator.cxx
 
 libinterpkernel_la_LDFLAGS= $(MED2_LIBS) $(HDF5_LIBS) ../MEDWrapper/V2_1/Core/libmed_V2_1.la $(STDLIB) ../MEDMEM/libmedmem.la  -lutil -lm
 
diff --git a/src/INTERP_KERNEL/NormalizedUnstructuredMesh.hxx b/src/INTERP_KERNEL/NormalizedUnstructuredMesh.hxx
index b14a3ea72..76a765028 100644
--- a/src/INTERP_KERNEL/NormalizedUnstructuredMesh.hxx
+++ b/src/INTERP_KERNEL/NormalizedUnstructuredMesh.hxx
@@ -9,8 +9,11 @@ namespace INTERP_KERNEL
       ALL_FORTRAN_MODE
     } NumberingPolicy;
 
+	class GenericMesh
+	{};
+
   template<int SPACEDIM, int MESHDIM, class ConnType, NumberingPolicy numPol, class MyMeshType>
-  class NormalizedUnstructuredMesh
+  class NormalizedUnstructuredMesh : public GenericMesh
   {
   public:
     void getBoundingBox(double *boundingBox) const { asLeaf().getBoundingBox(boundingBox); }
diff --git a/src/INTERP_KERNEL/PointLocator.cxx b/src/INTERP_KERNEL/PointLocator.cxx
new file mode 100644
index 000000000..9be00d3fe
--- /dev/null
+++ b/src/INTERP_KERNEL/PointLocator.cxx
@@ -0,0 +1,36 @@
+#include <list>
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Exception.hxx"
+#include "MEDNormalizedUnstructuredMesh.txx"
+#include "PointLocatorAlgos.txx"
+#include "PointLocator.hxx"
+
+namespace INTERP_KERNEL {
+PointLocator::PointLocator(const MEDMEM::MESH& mesh)
+{
+	int meshdim=mesh.getMeshDimension();
+	int spacedim=mesh.getSpaceDimension();
+	if (meshdim != spacedim) throw MEDMEM::MEDEXCEPTION("Locator is not implemented for meshdim != spacedim");
+	switch (meshdim)
+		{
+		case 2:
+			_medmesh = new MEDNormalizedUnstructuredMesh<2,2> (&mesh);
+			_point_locator=new PointLocatorAlgos<2,2,int,ALL_FORTRAN_MODE,MEDNormalizedUnstructuredMesh<2,2> >(*(static_cast<MEDNormalizedUnstructuredMesh<2,2>* >(_medmesh)));
+			break;
+		case 3:
+			_medmesh = new MEDNormalizedUnstructuredMesh<3,3> (&mesh);
+			_point_locator=new PointLocatorAlgos<3,3,int,ALL_FORTRAN_MODE,MEDNormalizedUnstructuredMesh<3,3> >(*(static_cast<MEDNormalizedUnstructuredMesh<3,3>* >(_medmesh)));
+			break;
+		}
+}
+PointLocator::~PointLocator()
+{
+	delete _medmesh;
+	delete _point_locator;
+}
+
+	std::list<int> PointLocator::locate(const double*x)
+{
+	return _point_locator->locates(x);
+}
+}
diff --git a/src/INTERP_KERNEL/PointLocator.hxx b/src/INTERP_KERNEL/PointLocator.hxx
new file mode 100644
index 000000000..43ec36fd3
--- /dev/null
+++ b/src/INTERP_KERNEL/PointLocator.hxx
@@ -0,0 +1,24 @@
+#ifndef _POINT_LOCATOR_HXX_
+#define _POINT_LOCATOR_HXX_
+#include <list>
+#include "NormalizedUnstructuredMesh.hxx"
+#include "MEDNormalizedUnstructuredMesh.hxx"
+#include "PointLocatorAlgos.txx"
+namespace INTERP_KERNEL
+{
+class PointLocator
+{
+public:
+	PointLocator(const MEDMEM::MESH& mesh);
+	virtual ~PointLocator();
+	std::list<int> locate(const double*x);
+	
+private:
+	//MEDNormalizedUnstructuredMesh<2,2>* _medmesh;
+	//GenericPointLocatorAlgos*< 2,2,int,ALL_FORTRAN_MODE,MEDNormalizedUnstructuredMesh<2,2> >* _point_locator;
+
+	GenericMesh* _medmesh;
+	GenericPointLocatorAlgos* _point_locator;
+};
+}
+#endif
diff --git a/src/INTERP_KERNEL/PointLocatorAlgos.txx b/src/INTERP_KERNEL/PointLocatorAlgos.txx
new file mode 100644
index 000000000..412f004e9
--- /dev/null
+++ b/src/INTERP_KERNEL/PointLocatorAlgos.txx
@@ -0,0 +1,181 @@
+#ifndef _POINT_LOCATOR_ALGOS_TXX_
+#define _POINT_LOCATOR_ALGOS_TXX_
+
+#include "MEDMEM_Exception.hxx"
+#include "InterpolationUtils.hxx"
+#include "MEDMEM_CellModel.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "BBTree.txx"
+#include <list>
+namespace INTERP_KERNEL{
+
+	class GenericPointLocatorAlgos
+	{
+	public:
+		virtual ~GenericPointLocatorAlgos(){};
+		virtual std::list<int> locates(const double* x)=0;
+	 
+	};
+	
+	template<int SPACEDIM, int MESHDIM, class ConnType, NumberingPolicy numPol, class MyMeshType> class PointLocatorAlgos: public GenericPointLocatorAlgos
+	{
+	private : 
+		double* _bb;
+		BBTree<SPACEDIM>* _tree;
+		const NormalizedUnstructuredMesh<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>& _mesh;
+
+	public:
+		PointLocatorAlgos<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>(const NormalizedUnstructuredMesh<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>& mesh):_mesh(mesh)
+		{
+			int nelem = _mesh.getNumberOfElements();
+			_bb = new double[SPACEDIM*2*nelem];
+			const ConnType* conn = _mesh.getConnectivityPtr();
+			const ConnType* conn_index = _mesh.getConnectivityIndexPtr();
+			const double* coords=_mesh.getCoordinatesPtr();
+			for (int i=0; i<nelem; i++)
+				{
+					for (int idim=0; idim<SPACEDIM; idim++)
+						{
+							_bb[2*(i*SPACEDIM+idim)]=HUGE;
+							_bb[2*(i*SPACEDIM+idim)+1]=-HUGE;
+						}
+					for (int index= conn_index[i]; index < conn_index[i+1];index++)
+						{
+							const double* coordelem = coords+OTT<ConnType,numPol>::ind2C(conn[OTT<ConnType,numPol>::ind2C(index)]);
+							for (int idim=0; idim<SPACEDIM;idim++)
+								{
+									_bb[2*(i*SPACEDIM+idim)]=(coordelem[idim]<_bb[2*(i*SPACEDIM+idim)])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)];
+									_bb[2*(i*SPACEDIM+idim)+1]=(coordelem[idim]>_bb[2*(i*SPACEDIM+idim)+1])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)+1];
+								}
+						}
+				}
+			
+			_tree=new BBTree<SPACEDIM>(_bb,0,0,nelem);
+		}
+		~PointLocatorAlgos<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>()
+		{
+			delete[] _bb;
+			delete _tree;
+		}
+	
+		std::list<int> locates(const double* x)
+		{
+			vector<int> candidates;
+			_tree->getElementsAroundPoint(x,candidates);
+			list<int> retlist;
+			for (int i=0; i< candidates.size(); i++)
+				{
+					if (elementContainsPoint(i,x))
+						retlist.push_back(OTT<ConnType,numPol>::indFC(i));
+				}
+			return retlist;
+		}
+	
+		bool elementContainsPoint(int i, const double*x)
+		{
+			//as i is extracted from the BBTRee, it is already in C numbering
+			//it is not necessary to convert it from F to C
+
+			const double* coords= _mesh.getCoordinatesPtr();
+			const int* conn=_mesh.getConnectivityPtr();
+			const int* conn_index= _mesh.getConnectivityIndexPtr();
+			const int* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
+
+		
+			int nbnodes = conn_index[i+1]-conn_index[i];
+		
+			// with dimension 2, it suffices to check all the edges
+			// and see if the sign of double products from the point
+			//is always the same.
+			//                 C
+			//                / \
+				//               /   \
+				//     Xo       /     \ 
+			//             A-------B
+			//
+			//here XA^XC and XC^XB have different signs
+			//
+			if (SPACEDIM==2)
+				{
+					//in 2D, nbedges==nbnodes
+					int nbedges=nbnodes;
+					int sign[nbedges];
+					for (int iedge=0; iedge<nbedges; iedge++)
+						{
+							const double* A=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
+							const double* B;
+							if (iedge+1< nbedges)
+								B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge+1]));
+							else
+								B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[0]));
+							
+							double a=mon_determinant(x, A, B);
+							if (a<-1e-12)
+								sign[iedge]=-1;
+							else if (a>1e-12)
+								sign[iedge]=1;
+							else
+								sign[iedge]=0;
+						}
+					return decide_from_sign(sign, nbedges);
+				}
+			
+			if (SPACEDIM==3)
+				{
+					MED_EN::medGeometryElement elem_type;
+					switch (nbnodes) {
+					case 4 :
+						elem_type=MED_EN::MED_TETRA4;
+						break;
+					case 5:
+						elem_type=MED_EN::MED_PYRA5;
+						break;
+					case 6 :
+						elem_type=MED_EN::MED_PENTA6;
+						break;
+					case 8: 
+						elem_type=MED_EN::MED_HEXA8;
+						break;
+					default:
+						throw MEDMEM::MEDEXCEPTION("PointLocatorAlgos : bad number of nodes in 3D locator");
+					}
+					const MEDMEM::CELLMODEL& model=MEDMEM::CELLMODEL_Map::retrieveCellModel(elem_type);
+					int nbfaces = model.getNumberOfConstituents(1);
+					int sign[nbfaces];
+					for (int iface=0; iface<nbfaces; iface++)
+						{
+							int* connface=model.getNodesConstituent(1,iface+1);
+							const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[0]-1]));
+							const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[1]-1]));
+							const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[2]-1]));
+							
+							double Vol=triple_product(AA,BB,CC,x);
+							if (Vol<-1e-12)
+								sign[iface]=-1;
+							else if (Vol>1e-12)
+								sign[iface]=1;
+							else
+								sign[iface]=0;
+						}
+					return  decide_from_sign(sign, nbfaces);
+				}
+			
+		}
+		
+		bool decide_from_sign (const int* sign, int nbelem)
+		{
+			int min_sign =1;
+			int max_sign =-1;
+			for (int i=0; i<nbelem;i++)
+				{
+					min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
+					max_sign=(sign[i]>max_sign)?sign[i]:max_sign;
+				}
+			if (min_sign==-1 && max_sign==1)
+				return false;
+			else
+				return true;	
+		}
+	};
+}
+#endif
diff --git a/src/INTERP_KERNEL/Test/BBTreeTest.cxx b/src/INTERP_KERNEL/Test/BBTreeTest.cxx
index 3180d7253..b1b0c6639 100644
--- a/src/INTERP_KERNEL/Test/BBTreeTest.cxx
+++ b/src/INTERP_KERNEL/Test/BBTreeTest.cxx
@@ -53,7 +53,12 @@ namespace INTERP_TEST
 	double bbox3[4]={5.0,6.0,7.0,9.0};
 	tree.getIntersectingElems(bbox3,elems);
 	CPPUNIT_ASSERT_EQUAL(2,(int)elems.size());
-	
+	elems.clear();
+
+	double xx[2]={1.0,1.0};
+	tree.getElementsAroundPoint(xx,elems);
+	CPPUNIT_ASSERT_EQUAL(4,(int)elems.size());
+
 	delete[] bbox;
   }
 
diff --git a/src/INTERP_KERNEL/Test/MEDMeshMaker.cxx b/src/INTERP_KERNEL/Test/MEDMeshMaker.cxx
new file mode 100644
index 000000000..d5168afdc
--- /dev/null
+++ b/src/INTERP_KERNEL/Test/MEDMeshMaker.cxx
@@ -0,0 +1,81 @@
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Meshing.hxx"
+
+MEDMEM::MESH* MEDMeshMaker(int dim, int nbedge, MED_EN::medGeometryElement type)
+{
+	MEDMEM::MESHING* mesh=new MEDMEM::MESHING();
+	mesh->setSpaceDimension(dim);
+	int nbnodes;
+	int nbelems;
+	switch (dim)
+		{
+		case 2: 
+			nbnodes=(nbedge+1)*(nbedge+1);
+			if(type==MED_EN::MED_QUAD4)
+				nbelems=(nbedge*nbedge);
+			else
+				throw MEDMEM::MEDEXCEPTION("MEDMeshMaker: type not impletmented");
+			break;
+		case 3:
+			nbnodes=(nbedge+1)*(nbedge+1)*(nbedge+1);
+			if (type==MED_EN::MED_HEXA8)
+				nbelems= nbedge*nbedge*nbedge;
+			else
+				throw MEDMEM::MEDEXCEPTION("MEDMeshMaker: type not impletmented");
+			break;
+		}
+	double* coords = new double[dim*nbnodes];
+	int nz;
+	if (dim==2) nz =1; else nz=nbedge+1;
+	{
+		for (int ix=0; ix < nbedge+1; ix++)
+			for (int iy=0; iy<nbedge+1; iy++)
+				for (int iz=0; iz<nz;iz++)
+					{
+						int inode=(ix*(nbedge+1)*nz+iy*nz+iz);
+						coords[inode*dim]=double(ix)/double(nbedge);
+						coords[inode*dim+1]=double(iy)/double(nbedge);
+						if (dim==3)
+							coords[inode*dim+2]=double(iz)/double(nbedge);
+					}
+	}
+	mesh->setCoordinates(dim, nbnodes,coords,"CARTESIAN",MED_EN::MED_FULL_INTERLACE);
+	
+	mesh->setNumberOfTypes(1,MED_EN::MED_CELL);
+	mesh->setTypes(&type,MED_EN::MED_CELL);
+	mesh->setNumberOfElements(&nbelems,MED_EN::MED_CELL);
+	
+	int* conn = new int [nbelems*(type%100)];
+	if (dim==2)
+		{
+			for (int ix=0; ix<nbedge; ix++)
+				for (int iy=0; iy<nbedge; iy++)
+					{
+						int ielem=(ix*nbedge+iy);
+						conn [ielem*4]=ix*(nbedge+1)+iy+1;
+						conn [ielem*4+1]=ix*(nbedge+1)+iy+1+1;
+						conn [ielem*4+2]=(ix+1)*(nbedge+1)+iy+1+1;
+						conn [ielem*4+3]=(ix+1)*(nbedge+1)+iy+1;															 
+					}
+		}
+	if (dim==3)
+		{
+			for (int ix=0; ix<nbedge; ix++)
+				for (int iy=0; iy<nbedge; iy++)
+					for (int iz=0; iz<nbedge; iz++)
+						{
+							int ielem=(ix*nbedge*nbedge+iy*nbedge+iz);
+							conn [ielem*8]=ix*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1;
+							conn [ielem*8+1]=(ix+1)*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1;
+							conn [ielem*8+2]=(ix+1)*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1;
+							conn [ielem*8+3]=ix*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1;
+							conn [ielem*8+4]=ix*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1+1;
+							conn [ielem*8+5]=(ix+1)*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1+1;
+							conn [ielem*8+6]=(ix+1)*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1+1;
+							conn [ielem*8+7]=ix*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1+1;
+						}
+		}
+	mesh->setConnectivity(conn, MED_EN::MED_CELL,type);
+	mesh->setMeshDimension(dim);
+	return mesh;
+}
diff --git a/src/INTERP_KERNEL/Test/MEDMeshMaker.hxx b/src/INTERP_KERNEL/Test/MEDMeshMaker.hxx
new file mode 100644
index 000000000..2c100f64c
--- /dev/null
+++ b/src/INTERP_KERNEL/Test/MEDMeshMaker.hxx
@@ -0,0 +1,7 @@
+#include "MEDMEM_define.hxx"
+namespace MEDMEM
+{
+	class MESH;
+}
+
+MEDMEM::MESH* MEDMeshMaker(int dim, int nbedge, MED_EN::medGeometryElement type);
diff --git a/src/INTERP_KERNEL/Test/Makefile.am b/src/INTERP_KERNEL/Test/Makefile.am
index 970a4ffd4..39aac6345 100644
--- a/src/INTERP_KERNEL/Test/Makefile.am
+++ b/src/INTERP_KERNEL/Test/Makefile.am
@@ -46,7 +46,9 @@ dist_libInterpKernelTest_la_SOURCES= \
      TransformedTriangleIntersectTest.cxx \
      BBTreeTest.cxx \
 		 RemapperTest.cxx \
-		 SingleElementPlanarTests.cxx
+		 SingleElementPlanarTests.cxx \
+     PointLocatorTest.cxx \
+		 MEDMeshMaker.cxx
 
 libInterpKernelTest_la_CPPFLAGS= @CPPUNIT_INCLUDES@ $(MED2_INCLUDES) $(HDF5_INCLUDES) \
 	-I$(srcdir)/.. -I$(srcdir)/../../MEDWrapper/V2_1/Core -I$(srcdir)/../../MEDMEM  -DOPTIMIZE -DLOG_LEVEL=0
@@ -72,3 +74,5 @@ if MED_ENABLE_KERNEL
 endif
 AM_CPPFLAGS= $(libInterpKernelTest_la_CPPFLAGS)
 
+check-local:TestInterpKernel
+	./TestInterpKernel
diff --git a/src/INTERP_KERNEL/Test/PointLocatorTest.cxx b/src/INTERP_KERNEL/Test/PointLocatorTest.cxx
new file mode 100644
index 000000000..ab924f125
--- /dev/null
+++ b/src/INTERP_KERNEL/Test/PointLocatorTest.cxx
@@ -0,0 +1,82 @@
+#include "PointLocatorTest.hxx"
+#include "PointLocator.hxx"
+#include "MEDMeshMaker.hxx"
+
+#include <iostream>
+#include <list>
+
+namespace INTERP_TEST
+{
+
+
+  void PointLocatorTest::setUp() 
+  {
+  }
+
+ 
+  void PointLocatorTest::tearDown() 
+  {
+  }
+
+  /**
+   * Test that creates a tree in 2D and check that 
+   * the results are correct in three
+   * cases :
+   * a non matching search
+   * a standard case
+   * a bbox overlapping the bboxes of the tree
+   */
+  void PointLocatorTest::test_PointLocator() {
+		MEDMEM::MESH* mesh2D= MEDMeshMaker(2,2,MED_EN::MED_QUAD4);
+		INTERP_KERNEL::PointLocator pl(*mesh2D) ;
+		double x[2]={0.0,0.0};
+		std::list<int> elems = pl.locate(x);
+		CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+		CPPUNIT_ASSERT_EQUAL(1,(int)(*(elems.begin())));
+		elems.clear();
+		
+		double x2[2]={0.25,0.25};
+		elems = pl.locate(x2);
+		CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+		CPPUNIT_ASSERT_EQUAL(1,(int)(*(elems.begin())));
+		elems.clear();
+		
+		double x3[2]={0.5,0.5};
+		elems = pl.locate(x3);
+		CPPUNIT_ASSERT_EQUAL(4,(int)elems.size());
+		elems.clear();
+
+		double x4[2]={-1.0,0.0};
+		elems = pl.locate(x4);
+		CPPUNIT_ASSERT_EQUAL(0,(int)elems.size());
+		elems.clear();
+
+		MEDMEM::MESH* mesh3D= MEDMeshMaker(3,2,MED_EN::MED_HEXA8);
+		INTERP_KERNEL::PointLocator pl3(*mesh3D);
+		double xx[3]={0.0,0.0,0.0};
+	  elems = pl3.locate(xx);
+		CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+		CPPUNIT_ASSERT_EQUAL(1,(int)*(elems.begin()));
+		elems.clear();
+		
+		double xx2[3]={0.25,0.25,0.25};
+		elems = pl3.locate(xx2);
+		CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+		CPPUNIT_ASSERT_EQUAL(1,(int)*(elems.begin()));
+			elems.clear();
+
+		double xx3[3]={0.5,0.5,0.5};
+		elems = pl3.locate(xx3);
+		CPPUNIT_ASSERT_EQUAL(8,(int)elems.size());
+		elems.clear();
+		
+		double xx4[3]={-1.0,0.0,0.0};
+		elems = pl3.locate(x4);
+		CPPUNIT_ASSERT_EQUAL(0,(int)elems.size());
+		elems.clear();
+
+
+  }
+
+
+}
diff --git a/src/INTERP_KERNEL/Test/PointLocatorTest.hxx b/src/INTERP_KERNEL/Test/PointLocatorTest.hxx
new file mode 100644
index 000000000..bc0027f4b
--- /dev/null
+++ b/src/INTERP_KERNEL/Test/PointLocatorTest.hxx
@@ -0,0 +1,39 @@
+#ifndef __TU_POINTLOCATOR_HXX__
+#define __TU_POINTLOCATOR_HXX__
+
+#include <cppunit/extensions/HelperMacros.h>
+#include "../PointLocator.hxx"
+
+namespace INTERP_TEST
+{
+
+  /**
+   * \brief Test suite testing some of the low level methods of TransformedTriangle.
+   *
+   */
+  class PointLocatorTest : public CppUnit::TestFixture
+  {
+
+    CPPUNIT_TEST_SUITE( PointLocatorTest );
+    CPPUNIT_TEST( test_PointLocator );
+    CPPUNIT_TEST_SUITE_END();
+
+   
+  public:
+    void setUp();
+
+    void tearDown();
+
+    // tests
+    void test_PointLocator();
+
+  };
+
+
+
+
+}
+
+
+
+#endif
diff --git a/src/INTERP_KERNEL/Test/RemapperTest.cxx b/src/INTERP_KERNEL/Test/RemapperTest.cxx
index ee78c0ff6..02c4775f7 100644
--- a/src/INTERP_KERNEL/Test/RemapperTest.cxx
+++ b/src/INTERP_KERNEL/Test/RemapperTest.cxx
@@ -60,7 +60,7 @@ namespace INTERP_TEST
 
 
 		//target square is in reverse order as compared to initial square
-		double source_integral=source_field.normL2(1,&source_areas);
+		double source_integral=source_field.normL2(1,&inv_source_areas);
 		double target_integral=target_field.normL2(1,&inv_target_areas);
 		
 		CPPUNIT_ASSERT_DOUBLES_EQUAL(source_integral,target_integral,1e-10);
diff --git a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx
index 0e2463620..1618714bb 100644
--- a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx
+++ b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx
@@ -25,6 +25,7 @@
 #include "SingleElementTetraTests.hxx"
 #include "HexaTests.hxx"
 #include "BBTreeTest.hxx"
+#include "PointLocatorTest.hxx"
 #include "RemapperTest.hxx"
 #include "MultiElement2DTests.hxx"
 #include "SingleElementPlanarTests.hxx"
@@ -39,6 +40,7 @@ CPPUNIT_TEST_SUITE_REGISTRATION( INTERP_TEST::TransformedTriangleIntersectTest )
 CPPUNIT_TEST_SUITE_REGISTRATION( INTERP_TEST::TransformedTriangleTest );
 CPPUNIT_TEST_SUITE_REGISTRATION( BBTreeTest);
 CPPUNIT_TEST_SUITE_REGISTRATION( RemapperTest);
+CPPUNIT_TEST_SUITE_REGISTRATION( PointLocatorTest);
 CPPUNIT_TEST_SUITE_REGISTRATION( MultiElement2DTests );
 CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementPlanarTests );