From: ndjinga <ndjinga>
Date: Tue, 9 Oct 2007 13:15:26 +0000 (+0000)
Subject: addition of ConvexIntersector and generic intersector for
X-Git-Tag: trio_trio_coupling~7
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=57cb95bbf66b8c34ddd15c591973346be766b4ff;p=tools%2Fmedcoupling.git

addition of ConvexIntersector and generic intersector for
computing 2D and 3D surf intersections
---

diff --git a/src/INTERP_KERNEL/Interpolation2D.cxx b/src/INTERP_KERNEL/Interpolation2D.cxx
index 9e614e225..ee90df2c2 100755
--- a/src/INTERP_KERNEL/Interpolation2D.cxx
+++ b/src/INTERP_KERNEL/Interpolation2D.cxx
@@ -1,413 +1,183 @@
-#include "InterpolationUtils.hxx"
-#include "TranslationRotationMatrix.hxx"
-#include "Interpolation.hxx"
-#include "Interpolation2D.hxx"
-
 #include "MEDMEM_Mesh.hxx"
-#include "MEDMEM_Field.hxx"
-
-#include "MEDMEM_InterpolationHighLevelObjects.hxx"
-#include "MEDMEM_Interpolation.hxx"
+#include "Interpolation2D.hxx"
+#include "PlanarIntersector.hxx"
+#include "PlanarIntersector.H"
+#include "TriangulationIntersector.hxx"
+#include "TriangulationIntersector.cxx"
+#include "ConvexIntersector.hxx"
+#include "ConvexIntersector.cxx"
+#include "GenericIntersector.hxx"
+#include "GenericIntersector.cxx"
+#include "BBTree.H"
+#include<time.h>
 
-using namespace std;
 using namespace MED_EN;
-using namespace MEDMEM;
 
-
-namespace MEDMEM 
+namespace MEDMEM
 {
 
-  Interpolation2D::Interpolation2D():_counter(0)
-  {
+  Interpolation2D::Interpolation2D()
+  { 
     _Precision=1.0E-12;
     _DimCaracteristic=1;
-    _SecondMethodOfLocalisation=1;
     _PrintLevel=0;
+    _Intersection_type= MEDMEM::Triangulation;
   }
-
+  
   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
   /*   Options :                                        */
   /*     Precision : for geometric computation          */
-  /*     SecondMethodOfLocalisation : 0/1               */
   /*     PrintLevel : between 0 and 3                   */
   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-  void Interpolation2D::setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel)
+  void Interpolation2D::setOptions(double Precision, int PrintLevel, IntersectionType intersection_type)
   {
     _Precision=Precision;
-    _SecondMethodOfLocalisation=SecondMethodOfLocalisation;
     _PrintLevel=PrintLevel;
+    _Intersection_type=intersection_type;
   }
-
-  /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-  /*localise le barycentre des mailles de P dans le maillage S*/
-  /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ 
-
-  void Interpolation2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,INTERP_UTILS::monMaillageP& MaP,
-					INTERP_UTILS::monMaillageS& MaS,vector<int>& localis )
-  {
-
-    //calcul des coordonnées des barycentre des mailles de P
-
-
-    //calcule des coordonnees du barycentre d'une cellule
-    std::auto_ptr<SUPPORT::SUPPORT>  Support ( new SUPPORT(&myMesh_P,"monSupport",MED_CELL) );
-    std::auto_ptr<FIELD<double, FullInterlace> > Barycentre (myMesh_P.getBarycenter(Support.get() ));
-    const  double* Bary=Barycentre->getValue();  
-
-
-    //localisation des barycentres des mailles de P dans les mailles de S.
-    Meta_Wrapper<2> fromWrapper(MaS._NbNoeud,(double *)MaS._Coord,(MEDMEM::CONNECTIVITY*)myMesh_S.getConnectivityptr());
-    Meta_Wrapper<2> toWrapper(MaP._NbMaille,(double *)Bary);
-    Meta_Mapping<2> mapping(&fromWrapper,&toWrapper);
-    mapping.Cree_Mapping(0);
-    for(int i=0;i<MaP._NbMaille;i++)
-      localis.push_back(mapping.Get_Mapping()[i]+1);
-
-
-    if (_SecondMethodOfLocalisation)
-      {
-	if(_PrintLevel)
-	  cout << endl << "Option de localisations complémentaires" << endl;
-	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-	//on vérifie que tout les barycentres ont obtenu un n° de maille 
-	//non nul.
-	MEDMEM::INTERPOLATION<2> monInterpol_S(myMesh_S);
-	for(int i_P=0;i_P<MaP._NbMaille;i_P++)
-	  {
-	    if(localis[i_P]<=0)
-	      {
-		//debug cout << endl << "------------------------------------------------" << endl << "Localise : barycentre maille " << i_P << " en dehors" << endl; 
-		int Test=0;
-
-		int NP_1=MaP._Connect[3*i_P];
-		int NP_2=MaP._Connect[3*i_P+1];
-		int NP_3=MaP._Connect[3*i_P+2];
-
-		double Coord_bary_i_P[2]={Bary[2*i_P],Bary[2*i_P+1]};
-		int Vois_N=monInterpol_S.getNearestNode(Coord_bary_i_P)+1;
-		int Ni=MaS._ReversNConnectIndex[Vois_N-1];
-		int Nf=MaS._ReversNConnectIndex[Vois_N];
-		int diff=Nf-Ni;
-
-		for(int j=0;j<diff;j++)
-		  {  
-		    int N_Maille=MaS._ReversNConnect[Ni-1+j];
-		    int NS_1=MaS._Connect[3*(N_Maille-1)];
-		    int NS_2=MaS._Connect[3*(N_Maille-1)+1];
-		    int NS_3=MaS._Connect[3*(N_Maille-1)+2];
-
-
-		    //debug cout << "mailles sources testées : " << N_Maille << endl;
-		    vector<double> Inter;
-		    INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),
-						       MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),
-						       MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1),
-						       Inter,_DimCaracteristic, _Precision);
-		    _counter++;
-		    if(Inter.size()>0)
-		      {
-			//debug cout << "Localise : maille proche barycentre trouvée : " << N_Maille << endl; 
-			Test=1;
-			localis[i_P]=N_Maille;
-			break;
-		      }
-		  }
-		if(Test==0)
-		  {
-		    int NoeuMaillP[3]={NP_1,NP_2,NP_3};
-		    for(int Num=0;Num<3;Num++)
-		      {
-			int Num_noeud=NoeuMaillP[Num];
-			double coord_i_P_0[2];
-			//VB
-			coord_i_P_0[0]= MaP._Coord[2*(Num_noeud-1)];
-			coord_i_P_0[1]=MaP._Coord[2*(Num_noeud-1)+1];
-			//VB
-			int Vois_B=monInterpol_S.getNearestNode(coord_i_P_0)+1;
-			int Ni=MaS._ReversNConnectIndex[Vois_B-1];
-			int Nf=MaS._ReversNConnectIndex[Vois_B];
-			int diff=Nf-Ni;
-
-			for(int j=0;j<diff;j++)
-			  {
-			    int N_Maille=MaS._ReversNConnect[Ni-1+j];
-			    int N_1=MaS._Connect[3*(N_Maille-1)];
-			    int N_2=MaS._Connect[3*(N_Maille-1)+1];
-			    int N_3=MaS._Connect[3*(N_Maille-1)+2];
-
-			    //debug cout << "mailles sources testées : " << N_Maille << endl;
-			    vector<double> Inter;
-			    INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(N_1-1),MaS._Coord+2*(N_2-1),
-							       MaS._Coord+2*(N_3-1),MaP._Coord+2*(NP_1-1),
-							       MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1),
-							       Inter, _DimCaracteristic, _Precision);	
-			    _counter++;
-			    if(Inter.size()>0)
-			      {
-				//debug cout << "Localise : maille proche sommet " << Num+1 << " trouvée : " << N_Maille << endl; 
-				Test=1;
-				localis[i_P]=N_Maille;
-				break;//on sort de la boucle sur les maille commune à un noeud de i_S
-			      }
-			  }
-			if(Test==1)
-			  {break;}//on sort de la boucle sur les noeuds de i_P
-		      }   
-		  }
-	      }
-	  }
-      }
-  }
-
-
-
-
-
-
-
-
-
-  /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ */
-  /*  fonction qui permet de remplir le vecteur donnant la surface d'intersection           */
-  /*  de la maille i_P du maillage P (projeté) avec la maille i_S du maillage S (support)  */   
-  /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ _ _ _ _ _ _ _ _ _ _  */
-
-
-  inline void Interpolation2D::rempli_vect_d_intersect(int i_P,int i_S1,const INTERP_UTILS::monMaillageP& MaP,const INTERP_UTILS::monMaillageS& MaS,	
-						       vector<map<int,double> >& Surface_d_intersect,FIELD<double>& myField_P)
-  {
-
-    
-    //on récupere le n° des noeuds.
-
-    //debug cout << "\nintersect maille " << i_P << endl;
-    int NP_1=MaP._Connect[3*(i_P-1)];
-    int NP_2=MaP._Connect[3*(i_P-1)+1];
-    int NP_3=MaP._Connect[3*(i_P-1)+2];
-
-
-    //on calcule la surface de la maille i_P
-
-    double Surf_i_P =INTERP_UTILS::Surf_Tri(MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),
-					    MaP._Coord+2*(NP_3-1));
-
-    double Surface = 0 ;
-    vector<int> maill_deja_vu_S ;
-
-    maill_deja_vu_S.push_back(i_S1);
-
-    for (int M_S=0; M_S<maill_deja_vu_S.size(); M_S++)
-      {
-	if( abs(Surf_i_P-Surface)/Surf_i_P>_Precision && M_S!=maill_deja_vu_S.size() )
-	  {
-	    int i_S=maill_deja_vu_S[M_S];
-
-	    int NS_1=MaS._Connect[3*(i_S-1)];
-	    int NS_2=MaS._Connect[3*(i_S-1)+1];
-	    int NS_3=MaS._Connect[3*(i_S-1)+2];
-
-
-	    vector<double> Inter;
-	    INTERP_UTILS::intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),
-					       MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),
-					       MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1),
-					       Inter,_DimCaracteristic, _Precision);
-	    _counter++;
-
-	    //on teste l'intersection.
-	    int taille=Inter.size()/2;
-	    //debug cout << "  -> maille source : " << i_S << " , nb intersection=" <<  taille << endl;
-
-	    /* CAS1  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-	    // taille ==0, on ne fait rien
-
-	    /* CAS 2  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-	    if(taille==1 || taille==2)
-	      {
-		//on ne remplit pas le vecteur des correspondances n° maille-surfaces 
-		//d'intersections mais on récupère le numéro des mailles voisines de la maille i_S.
-		vector<int> M_Vois=INTERP_UTILS::touv_maill_vois(i_S,MaS);
-		for(int i=0;i< M_Vois.size();i++)
-		  {INTERP_UTILS::verif_maill_dans_vect(M_Vois[i],maill_deja_vu_S);}
-	      }
-
-	    /* CAS 3  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
-	    else if(taille==3)
-	      {
-		double Surf_inter = INTERP_UTILS::Surf_Poly(Inter) ;
-
-		//on remplit le vecteur des correspondances n° maille-surfaces d'intersections.
-		Surface_d_intersect[i_P-1][i_S]=Surf_inter;
-				
-		// on récupère le numéro des mailles voisines de la maille i_S.
-		vector<int> M_Vois=INTERP_UTILS::touv_maill_vois(i_S,MaS);
-		for(int i_M=0;i_M<(M_Vois).size();i_M++)
-		  {INTERP_UTILS::verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
-
-		Surface = Surface + Surf_inter ; 
-	      }
-			
-	    /* CAS 4  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-	    else if (taille > 3) //taille>3
-	      {
-		vector<double> Poly_Inter=INTERP_UTILS::reconstruct_polygon(Inter);
-
-		double Surf_inter =INTERP_UTILS::Surf_Poly(Poly_Inter) ;
-
-		//on remplit le vecteur des correspondances n° maille-surfaces d'intersection
-		(Surface_d_intersect[i_P-1])[i_S]=Surf_inter;
-
-		// on récupère le numéro des mailles voisines de la maille i_S.
-		vector<int> M_Vois=INTERP_UTILS::touv_maill_vois(i_S,MaS);
-		for(int i_M=0;i_M<(M_Vois).size();i_M++)
-		  {INTERP_UTILS::verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
-
-		Surface = Surface + Surf_inter ; 
-	      }
-	  }
-	//debug cout << "     surface = " << Surface << endl << flush;
-      }
-
-    /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //on rentre la fraction totale de la maille i_P qui a été considérer lors de l'interpolation.
-    double Fract=Surface/Surf_i_P;
-    myField_P.setValueIJ(i_P,1,Fract);
-
-  }
-
-
-  /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-  //fonction principale pour interpoler deux maillages triangulaires.
+  
+  
+  /***************************************************************/
+  /* Main function to interpolate triangular or quadratic meshes */
+  /***************************************************************/
   vector<map<int,double> > Interpolation2D::interpolateMeshes(const MEDMEM::MESH& myMesh_S,
-							       const MEDMEM::MESH& myMesh_P)
+								   const MEDMEM::MESH& myMesh_P)
   {
-
-    long global_start =clock();    
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    // on vérifie que les deux maillages sont formés de triangles.
-    int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_CELL);
-    int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_CELL);
-    if  ( NumberOfCellsTypes_S != 1  || NumberOfCellsTypes_P != 1)
-      { cout<<"l'un au moins des deux maillages n'est pas triangulaires"<<endl;}
-    string* Type_S = myMesh_S.getCellTypeNames(MED_CELL);
-    string* Type_P = myMesh_P.getCellTypeNames(MED_CELL);
-    if ( Type_S[0] != "MED_TRIA3" || Type_P[0] != "MED_TRIA3")
-      { cout<<"l'un au moins des deux maillages n'est pas triangulaires"<<endl;}
+    long global_start =clock(); 
+    int counter=0;   
+    /***********************************************************/
+    /* Check both meshes are made of triangles and quadrangles */
+    /***********************************************************/
+ 
+    int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_EN::MED_CELL);
+    int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_EN::MED_CELL);
+    if  ( NumberOfCellsTypes_S > 2  || NumberOfCellsTypes_P >2)
+      cout<<"Interpolation2D::intersectMeshes: one of the two meshes contains more than two types of cells"<<endl;
+    string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
+    string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
+    
+    if((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0] != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
+      cout<<"Interpolation2D::intersectMeshes: one of the two meshes is neither linear triangular nor linear quadratic"<<endl;
     //VB
     delete[] Type_S;
     delete[] Type_P;
-
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    INTERP_UTILS::monMaillageS MaS;
-    INTERP_UTILS::monMaillageP MaP;
-
-    MaS._NbNoeud = myMesh_S.getNumberOfNodes() ;
-    MaP._NbNoeud = myMesh_P.getNumberOfNodes() ;
-    MaS._NbMaille =myMesh_S.getNumberOfElements(MED_CELL,MED_TRIA3);
-    MaP._NbMaille =myMesh_P.getNumberOfElements(MED_CELL,MED_TRIA3);
-
-    //on charge les connectivités des deux maillages.
-    MaS._Connect =myMesh_S.getConnectivity(MED_FULL_INTERLACE, 
-					   MED_NODAL, MED_CELL, MED_TRIA3) ;
-    MaP._Connect =myMesh_P.getConnectivity(MED_FULL_INTERLACE, 
-					   MED_NODAL, MED_CELL, MED_TRIA3) ;
-
-    //on charge les coordonnées des noeuds.
-    MaS._Coord = myMesh_S.getCoordinates(MED_FULL_INTERLACE);
-    MaP._Coord = myMesh_P.getCoordinates(MED_FULL_INTERLACE);
-
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
-    //on a besoin de recupérer connectivité inverse de S.
-    MaS._ReversNConnect=
-      myMesh_S.getReverseConnectivity(MED_NODAL);
-    MaS._ReversNConnectIndex=
-      myMesh_S.getReverseConnectivityIndex(MED_NODAL);
-
-
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //on cherche la dimension caractéristique des maillages
+    int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+    int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+    
+    /**************************************************/
+    /* Search the characteristic size of the meshes   */
+    /**************************************************/
+    
     vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
     vector<vector<double> > BoxP=myMesh_P.getBoundingBox();
-    double AreaS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0])+(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1]));
-    MaS._DimCaracteristic=sqrt(4./sqrt(3.)*AreaS/MaS._NbMaille);
-    double AreaP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0])+(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1]));
-    MaP._DimCaracteristic=sqrt(4./sqrt(3.)*AreaP/MaP._NbMaille);
-    _DimCaracteristic=min(MaS._DimCaracteristic,MaP._DimCaracteristic);
-    if (_PrintLevel)
+    double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0])
+			 +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1]));
+    double DimCaracteristic_S=diagonalS/nbMaille_S;
+    double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0])
+			 +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1]));
+    double DimCaracteristic_P=diagonalP/nbMaille_P;
+    
+    _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P);
+    if (_PrintLevel>=1)
       {
-	cout<<"recherche de la dimension caractéristique des maillages 2D :"<<endl;
-	cout<<"  - dimension caractéristique du maillage source : "<<MaS._DimCaracteristic<<endl;
-	cout<<"  - dimension caractéristique du maillage projeté: "<<MaS._DimCaracteristic<<endl;
-	cout<<"  - d'où la dimension caractéristique: "<<_DimCaracteristic<<endl;
+	cout<<"  - Characteristic size of the source mesh : "<< DimCaracteristic_S<<endl;
+	cout<<"  - Characteristic size of the target mesh: "<< DimCaracteristic_P<<endl;
+	cout << "Interpolation2D::computation of the intersections"<<endl;
+      }
+    
+    PlanarIntersector* intersector;
+    
+    switch (_Intersection_type)
+      {
+      case  MEDMEM::Triangulation:
+	intersector=new TriangulationIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, 
+						    0, _PrintLevel);
+	break;
+      case MEDMEM::Convex:
+	intersector=new ConvexIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision,
+					     0, 0, _PrintLevel);
+	break;
+      case MEDMEM::Generic:
+	intersector=new GenericIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision,
+					      0, 0, _PrintLevel);
+	break;
       }
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //on cherche pour chaque maille i_P du maillage projetté
-    // une maille i_S du maillage S qui chevauche la maille i_P.
 
+    /****************************************************************/
+    /* Create a search tree based on the bounding boxes             */
+    /* Instanciate the intersector and initialise the result vector */
+    /****************************************************************/
+ 
     long start_filtering=clock();
-    vector<int> localis;
-    localis.reserve(MaP._NbMaille);
-    MEDMEM::MESH* ptr_S = const_cast<MEDMEM::MESH*>(&myMesh_S);
-    MEDMEM::MESH* ptr_P = const_cast<MEDMEM::MESH*>(&myMesh_P);
-    
-    cout << "Interpolation2D::local_iP_dansS"<<endl;
-    local_iP_dans_S(*ptr_S,*ptr_P,MaP,MaS,localis);
+ 
+    vector<double> bbox;
+    intersector->createBoundingBoxes<2>(myMesh_S,bbox); // create the bounding boxes
+    BBTree<2> my_tree(&bbox[0], 0, 0,nbMaille_S );//creating the search structure 
+
     long end_filtering=clock();
 
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //on crée un tableau où l'on rentrera la valeurs des intersections.
-    cout << "Interpolation2D::calcul intersec"<<endl;
     map<int,double> MAP_init;
-    vector<map<int,double> > Surface_d_intersect(MaP._NbMaille,MAP_init);//on initialise.
-
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //création d'un champ pour rentrer les fractions de mailles considérées.
-
-    std::auto_ptr<SUPPORT> mySupport_P (new SUPPORT(const_cast<MEDMEM::MESH*>(&myMesh_P),"Support on all CELLS",MED_CELL) );
-    MEDMEM::FIELD<double> myField_P(mySupport_P.get(),1);
-    string NameF=" M2 dans M1" ;
-    myField_P.setName(NameF);
-    string NameC="fracSurf";
-    string ComponentsNames[1]={NameC};
-    myField_P.setComponentsNames(ComponentsNames);
-    myField_P.setComponentsDescriptions(ComponentsNames);
-    string ComponentsUnits[1]={"%"};
-    myField_P.setMEDComponentsUnits(ComponentsUnits);
+    vector<map<int,double> > result_areas(nbMaille_P,MAP_init);//on initialise.
 
-    /*à remplacer par:
-      WriteFractSurfInFile(vector<<vector<pair<int,double> > >& FractVol,MEDMEM::MESH* myMesh1,
-      MEDMEM::MESH* myMesh2,string FileName,string NameMeshP,string NameMeshS)
-    */
+    /****************************************************/
+    /* Loop on the target cells - core of the algorithm */
+    /****************************************************/
+    const MED_EN::medGeometryElement* types = myMesh_P.getTypes(MED_EN::MED_CELL);
+    int type_max_index=0;//maximum cell number for a given type
+    int type_min_index=0;//minimum cell number for a given type
+    int i_P=0;//global index of cell
 
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //boucle sur les mailles de P.
-    //Coeur de l'algorithme
     long start_intersection=clock();
-    for(int i_P=0; i_P<MaP._NbMaille ; i_P++)
-      {
-	int i_S_init=localis[i_P];
-	if(i_S_init>0)
-	  rempli_vect_d_intersect(i_P+1,i_S_init,MaP,MaS,Surface_d_intersect,myField_P);
 
+    for (int itype = 0; itype<NumberOfCellsTypes_P; itype++)
+      {
+	int nb_nodesP=types[itype]%100;
+	int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
+	type_max_index +=  nbelem_type;
+	for( i_P=type_min_index; i_P<type_max_index ; i_P++)
+	  {
+	    vector<int> intersecting_elems;
+	    double bb[4];
+	    intersector->getElemBB<2>(bb,myMesh_P,i_P+1,nb_nodesP);
+	    my_tree.getIntersectingElems(bb, intersecting_elems);
+	    int nb_intersecting_elems = intersecting_elems.size();	    
+	    //browsing all the i_S (from mesh S) elems that can 
+	    //intersect elem i_P (from mesh P)
+	    for (int ielem=0; ielem<nb_intersecting_elems;ielem++)
+	      {
+		//BBTree structure returns numbers between 0 and n-1
+		int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
+		int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL,i_S+1)%100;
+		double surf=intersector->intersectCells(i_P+1,i_S+1,nb_nodesP,nb_nodesS); 
+		(result_areas[i_P]).insert(make_pair(i_S+1,surf));
+		counter++;
+		//rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas);
+	      }
+	    intersecting_elems.clear();
+	  }
+	type_min_index = type_max_index;
       }
-    if (_PrintLevel >= 1)
+    delete intersector;
+
+    /***********************************/
+    /*        DEBUG prints             */
+    /***********************************/
+
+    if (_PrintLevel >=1)
       {
 	long end_intersection=clock();
-	if (_PrintLevel >= 2)
+	if (_PrintLevel >=2)
 	  {
 	    cout<<endl<<"Printing intersection areas:"<<endl<<endl;
 	    cout<<"(source cell, target cell): intersection areas"<<endl;
 	    double total=0.0;
 	    double total_interm=0.0;
-	    int nb_result_areas = Surface_d_intersect.size();
+	    int nb_result_areas = result_areas.size();
 	    for(int i=0; i< nb_result_areas;i++)
 	      { 
 		map<int,double>::iterator surface;
 		total_interm=0.0;
-		for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++)
+		for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++)
 		  {
 		    cout<<"    ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
 		    total_interm +=(*surface).second;
@@ -417,54 +187,13 @@ namespace MEDMEM
 	      }
 	    cout << "total area "<<total<<endl;
 	  }
-	cout<< "Barycenter localisation time= " << end_filtering-start_filtering<< endl;
+	cout<< "Filtering time= " << end_filtering-start_filtering<< endl;
 	cout<< "Intersection time= " << end_intersection-start_intersection<< endl;
-	cout<< "counter= " << _counter << endl;
 	long global_end =clock();    
+      	cout<< "Number of computed intersections = " << counter << endl;
 	cout<< "Global time= " << global_end - global_start << endl;
       }
     
-    //    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //    //On sauvegarde le champ crée sur la maillage P.
-    //    if (_PrintLevel>=3)
-    //    {
-    ////	string NameFile="fractSurf.med";
-    ////	int id1=myMesh_P.addDriver(MED_DRIVER,NameFile,"M1");
-    ////	int id2=myMesh_S.addDriver(MED_DRIVER,NameFile,"M2");
-    ////	int id3=myField_P.addDriver(MED_DRIVER,NameFile,NameF);
-    ////	myMesh_P.write(id1);
-    ////	myMesh_S.write(id2);
-    ////	myField_P.write(id3);
-    ////
-    ////	cout<<endl<<"Pour chaque mailles de " <<myMesh_S.getName();
-    ////	cout<<"Impression de la fraction surfaciques totale dans le maillage : "<<myMesh_P.getName()<<endl; 
-    ////	for(int i=0;i<MaP._NbMaille;i++)
-    ////	    cout<<"maille n°"<<i+1<<" de "<<myMesh_P.getName()<<": "<<myField_P.getValueIJ(i+1,1)<<endl;
-    //    }
-
-    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    //   for (int i=0; i< Surface_d_intersect.size();i++)
-    //       {
-    // 	double size=0;
-    // 	for (map<int,double>::const_iterator iter = Surface_d_intersect[i].begin();
-    // 	     iter !=Surface_d_intersect[i].end();
-    // 	     iter++)
-    // 	  size+= iter->second;
-    // 	//debug cout << "\nintersect maille " << i_P << endl;
-    // 	int NP_1=MaP._Connect[3*(i)];
-    // 	int NP_2=MaP._Connect[3*(i)+1];
-    // 	int NP_3=MaP._Connect[3*(i)+2];
-
-
-    //     //on calcule la surface de la maille i_P
-
-    // 	double Surf_i_P =Surf_Tri(MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),
-    // 				  MaP._Coord+2*(NP_3-1));	
-    // 	cout <<size<<" "<<Surf_i_P<<endl;
-    //       }
-    return Surface_d_intersect;
-
+    return result_areas;
   }
-
-
 };
diff --git a/src/INTERP_KERNEL/Interpolation2D.hxx b/src/INTERP_KERNEL/Interpolation2D.hxx
index 2905f10ff..374b213ef 100755
--- a/src/INTERP_KERNEL/Interpolation2D.hxx
+++ b/src/INTERP_KERNEL/Interpolation2D.hxx
@@ -1,76 +1,33 @@
-#ifndef _INTERPOLATION2D_HXX_
-#define _INTERPOLATION2D_HXX_
+#ifndef _INTERPOLATION_2D_SURF_HXX_
+#define _INTERPOLATION_2D_SURF_HXX_
 
 #include "MEDMEM_Mesh.hxx"
-#include "MEDMEM_Field.hxx"
 #include "Interpolation.hxx"
 
-#include <cmath>
-#include <map>
-#include <sstream>
-#include <string>
-#include <set>
-#include <algorithm>
-#include <vector>
-#include <complex>
-
-// namespace MEDMEM 
-// {
-//   class MESH;
-//   template <typename T> class FIELD;
-// }
-
-namespace INTERP_UTILS
-{
-  struct monMaillageS;
-  struct monMaillageP;
-};
-
 namespace MEDMEM
 {
- 
-  
   class Interpolation2D : public Interpolation
   {
-    
   private: 
     double _Precision;
     double _DimCaracteristic;
-    int  _SecondMethodOfLocalisation;
     int  _PrintLevel;
-    IntersectionType _intersection_type;
-    int _counter;
-    // Méthodes publiques
+    IntersectionType _Intersection_type;
+
   public:
-    
     Interpolation2D();
     
-    // precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
-    void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
-    
-    //declare Intersection polygons as convex
-    void setConvex(bool flag)
-    {
-      if (flag) _intersection_type=Convex;
-    } 
- 
-    //localise le barycentre des mailles de P dans le maillage S
-    void local_iP_dans_S(MEDMEM::MESH& myMesh_S,MEDMEM::MESH& myMesh_P,
-			 INTERP_UTILS::monMaillageP& MaP,INTERP_UTILS::monMaillageS& MaS,vector<int>& localis );
-    
-    
-    //fonction qui permet de remplir le vecteur donnant la surface d'intersection 
-    //de la mailles i_P du maillage projetté avec la maille i_S du maillage support.
-    inline void rempli_vect_d_intersect(int i_P,int i_S,
-					const INTERP_UTILS::monMaillageP& MaP,const INTERP_UTILS::monMaillageS& MaS,
-					vector<map<int,double> >& Surface_d_intersect,FIELD<double>& myField_P);
-    
-    
-    //fonction principale pour interpoler deux maillages triangulaires.
-    std::vector<map<int, double> >interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+    // geometric precision, debug print level, coice of the median plane, intersection etc ...
+    void setOptions(double Precision, int PrintLevel, 
+		    IntersectionType intersectionType);
     
+    // Main function to interpolate triangular and qudadratic meshes
+    std::vector<map<int, double> > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+
+  private: 
+
   };
-  
+
 };
 
 #endif
diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in
index f0f41e270..fd38c3a8a 100644
--- a/src/INTERP_KERNEL/Makefile.in
+++ b/src/INTERP_KERNEL/Makefile.in
@@ -49,7 +49,6 @@ BoundingBox.hxx\
 TetraAffineTransform.hxx\
 TranslationRotationMatrix.hxx\
 Interpolation2D.hxx\
-Interpolation2Dbis.hxx\
 Interpolation3DSurf.hxx\
 InterpolationUtils.hxx\
 BBTree.H\
@@ -86,7 +85,6 @@ TetraAffineTransform.cxx\
 Interpolation3D.cxx\
 Interpolation3DSurf.cxx\
 Interpolation2D.cxx\
-Interpolation2Dbis.cxx\
 IntersectorTetra.cxx\
 IntersectorHexa.cxx\
 PlanarIntersector.cxx
@@ -98,7 +96,7 @@ BIN_SRC =
 BIN_SERVER_IDL = 
 BIN_CLIENT_IDL = 
 
-TEST_PROGS = testUnitTetra test_PolygonAlgorithms test_Interpol_2D test_Interpol_2Dbis  test_Interpol_3DSurf
+TEST_PROGS = testUnitTetra test_PolygonAlgorithms test_Interpol_2D test_Interpol_3DSurf
 
 LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome 
 LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
diff --git a/src/INTERP_KERNEL/test_Interpol_2D.cxx b/src/INTERP_KERNEL/test_Interpol_2D.cxx
index 0eb6f4573..d8e914a7a 100644
--- a/src/INTERP_KERNEL/test_Interpol_2D.cxx
+++ b/src/INTERP_KERNEL/test_Interpol_2D.cxx
@@ -1,12 +1,40 @@
 #include "Interpolation2D.hxx"
 #include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Field.hxx"
 
-int main()
+int main(int argc, char** argv) 
 {
-  MEDMEM::MESH source(MED_DRIVER,"/home/vb144235/resources/square128000.med","Mesh_1");
-  MEDMEM::MESH target(MED_DRIVER,"/home/vb144235/resources/square30000.med","Mesh_1");
+  if (argc <2) { cout <<"erreur nb args "<<endl; exit(1);}
+
+  IntersectionType intersectiontype;
+  int type = atoi(argv[1]);
+
+  switch(type)
+    {
+    case 0:
+      intersectiontype=MEDMEM::Triangulation;
+      break;
+    case 1 :
+      intersectiontype=MEDMEM::Convex; 
+      break;
+    case 2:
+      intersectiontype=MEDMEM::Generic; 
+      break;
+    default :
+      cout <<"error second argument should be 0, 1 or 2 (triangles or convex) "<<endl; exit(1);
+    }
+
+  MEDMEM::MESH target(MED_DRIVER,"/home/vb144235/resources/square128000.med","Mesh_1");
+  MEDMEM::MESH source(MED_DRIVER,"/home/vb144235/resources/square30000.med","Mesh_1");
+
+  MEDMEM::Interpolation2D interp; 
+ 
+  interp.setOptions(1e-6,2,intersectiontype);
+  vector<map<int,double> > surfaces= interp.interpolateMeshes(target,source); 
+
+//   map<int,double>::iterator mi;
+//    cerr<< "Résultat interpolation: celulle " << 0 << endl;
+//    for(mi=(surfaces[0]).begin(); mi!=(surfaces[0]).end(); mi++)
+//      cerr<< (*mi).first << " ; " << (*mi).second << endl;
 
-  MEDMEM::Interpolation2D interp;
-  interp.setOptions(1e-6,1,2);
-  interp.interpolateMeshes(source,target);
 }
diff --git a/src/INTERP_KERNEL/test_Interpol_3DSurf.cxx b/src/INTERP_KERNEL/test_Interpol_3DSurf.cxx
new file mode 100644
index 000000000..8e3e6a544
--- /dev/null
+++ b/src/INTERP_KERNEL/test_Interpol_3DSurf.cxx
@@ -0,0 +1,41 @@
+#include "Interpolation3DSurf.hxx"
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Field.hxx"
+
+int main(int argc, char** argv) 
+{
+  if (argc <4) { cout <<"erreur nb args "<<endl; exit(1);}
+
+  double medianPlane = atof(argv[1]);
+  IntersectionType intersectiontype;
+  int type = atoi(argv[2]);
+  bool rotate = atoi(argv[3])==1;
+
+  switch(type)
+    {
+    case 0:
+      intersectiontype=MEDMEM::Triangulation;
+      break;
+    case 1 :
+      intersectiontype=MEDMEM::Convex; 
+      break;
+    case 2:
+      intersectiontype=MEDMEM::Generic; 
+      break;
+    default :
+      cout <<"error second argument should be 0, 1 or 2 (triangles or convex) "<<endl; exit(1);
+    }
+
+  MEDMEM::MESH target(MED_DRIVER,"/data/tmpawa/ndjinga/MEDPARA/test/cylinderquad.med","cylinder");
+  MEDMEM::MESH source(MED_DRIVER,"/data/tmpawa/ndjinga/MEDPARA/test/cylindertri.med","cylinder");
+
+  MEDMEM::Interpolation3DSurf interp;
+
+   interp.setOptions(1e-6,2,medianPlane,intersectiontype, rotate); 
+   vector<map<int,double> > surfaces= interp.interpolateMeshes(source,target);
+
+//   map<int,double>::iterator mi;
+//    cerr<< "Résultat interpolation: celulle " << 0 << endl;
+//    for(mi=(surfaces[0]).begin(); mi!=(surfaces[0]).end(); mi++)
+//      cerr<< (*mi).first << " ; " << (*mi).second << endl;
+}