From: vbd <vbd>
Date: Mon, 30 Apr 2007 14:57:44 +0000 (+0000)
Subject: improvement to the ParaMEDMEM classes with the possibility to address
X-Git-Tag: trio_trio_coupling~105
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=473c7e6d9ff83b5ac997f94ac6c2ac317da47efb;p=tools%2Fmedcoupling.git

improvement to the ParaMEDMEM classes with the possibility to address
coupling of 3D surfaces

correction of a bug in the reconstruct_polygon algorithm
---

diff --git a/src/ParaMEDMEM/BBTree.H b/src/ParaMEDMEM/BBTree.H
new file mode 100644
index 000000000..3bc3e36bc
--- /dev/null
+++ b/src/ParaMEDMEM/BBTree.H
@@ -0,0 +1,134 @@
+#include <iostream> 
+#include <vector>
+#include <algorithm>
+
+const int MIN_NB_ELEMS =10;
+const int MAX_LEVEL=10;
+using namespace std;
+
+template <int dim> class BBTree
+{
+
+private:
+  BBTree* _left;
+  BBTree* _right;
+  BBTree* _center;
+  int _level;
+  double _median;
+  double* _bb;
+  vector <int> _elems;
+  bool _terminal;
+
+public:
+
+
+BBTree(double* bbs, int* elems, int level, int nbelems):
+  _left(0), _right(0), _center(0), _level(level), _median(0.0), _bb(bbs), _terminal(false)
+{
+  if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
+    {
+      _terminal=true;
+		
+    }
+  vector<double> nodes (nbelems*2);
+  for (int i=0; i<nbelems; i++)
+    {
+      int elem;
+      if (elems!=0)
+	elem= elems[i];
+      else
+	elem=i;
+
+      _elems.push_back(elem);
+      nodes[i*2]=bbs[elem*dim*2+(level%dim)*2];
+      nodes[i*2+1]=bbs[elem*dim*2+(level%dim)*2+1];
+    }
+  if (_terminal) return;
+
+	//elem list is not useful if the node is not terminal
+	_elems.clear();
+
+  vector<double>::iterator median = nodes.begin()+nbelems;
+   nth_element< vector<double>::iterator>(nodes.begin(), nodes.begin()+nbelems, nodes.end());
+   _median=*median;
+   // std:: cout << *median <<std::endl;
+
+  vector<int> new_elems_left;
+  vector<int> new_elems_right;
+  vector<int> new_elems_center;
+  for (int i=0; i<nbelems;i++)
+    {
+      int elem;
+      if (elems!=0)
+	elem= elems[i];
+      else
+	elem=i;
+      double max=bbs[elem*dim*2+(level%dim)*2+1];
+      if (max < *median)
+	{
+	  new_elems_left.push_back(elem);
+	  continue;
+	}
+      double min = bbs[elem*dim*2+(level%dim)*2];
+      if (min>*median)
+	{
+	  new_elems_right.push_back(elem);
+	  continue;
+	}
+      new_elems_center.push_back(elem);
+    }
+
+  _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());
+  _center=new BBTree(bbs, &(new_elems_center[0]), level+1, new_elems_center.size());
+ 
+  
+}
+
+~BBTree()
+{
+  if (_left!=0)  delete _left;
+  if (_right!=0) delete _right;
+  if (_center!=0) delete _center;
+}
+
+void getIntersectingElems(double* bb, vector<int>& elems)
+{
+  //  terminal node : return list of elements intersecting bb
+  if (_terminal)
+    {
+      for (int i=0; i<_elems.size(); i++)
+	{
+	  double* bb_ptr=_bb+_elems[i]*2*dim;
+	  bool intersects = true;
+	  for (int idim=0; idim<dim; idim++)
+	    {
+	      if (bb_ptr[idim*2]>bb[idim*2+1] || bb_ptr[idim*2+1]<bb[idim*2])
+		intersects=false;
+	    }
+	  if (intersects)
+	    elems.push_back(_elems[i]);
+	}
+      return;
+    }
+
+  //non terminal node 
+  double min = bb[(_level%dim)*2];
+  double max = bb[(_level%dim)*2+1];
+  if (max < _median)
+    {
+      _left->getIntersectingElems(bb, elems);
+			_center->getIntersectingElems(bb,elems);
+      return;
+    }
+  if (min > _median)
+    {
+      _right->getIntersectingElems(bb,elems);
+			_center->getIntersectingElems(bb,elems);
+      return;
+    }
+  _left->getIntersectingElems(bb,elems);
+  _center->getIntersectingElems(bb,elems);
+  _right->getIntersectingElems(bb,elems);
+}
+};
diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx
index 74e344b79..60281fe22 100644
--- a/src/ParaMEDMEM/ElementLocator.cxx
+++ b/src/ParaMEDMEM/ElementLocator.cxx
@@ -19,7 +19,8 @@ const double HUGE = 1e300;
 ElementLocator::ElementLocator(const ParaMESH& mesh, const ProcessorGroup& distant_group) 
 :_local_mesh(mesh.getMesh()),
 _local_group(*mesh.getBlockTopology()->getProcGroup()),
-_distant_group(distant_group)
+ _distant_group(distant_group),
+ _adjustment_eps(0.1)
 { 
   _union_group = _local_group.fuse(distant_group);
   _computeBoundingBoxes();
@@ -162,9 +163,22 @@ bool ElementLocator::_intersectsBoundingBox(int irank)
 
 bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim)
 {
+	double bbtemp[2*dim];
+	double deltamax=0.0;
+	for (int i=0; i< dim; i++)
+		{
+			double delta = bb1[2*i+1]-bb1[2*i];
+			deltamax = (delta>deltamax)?delta:deltamax;
+		}
+	for (int i=0; i<dim; i++)
+		{
+			bbtemp[i*2]=bb1[i*2]-deltamax*_adjustment_eps;
+			bbtemp[i*2+1]=bb1[i*2+1]+deltamax*_adjustment_eps;
+		}
+	
   for (int idim=0; idim < dim; idim++)
   {
-    bool intersects = bb1[idim*2]<bb2[idim*2+1] && bb2[idim*2]<bb1[idim*2+1];
+    bool intersects = bbtemp[idim*2]<bb2[idim*2+1] && bb2[idim*2]<bbtemp[idim*2+1];
     if (!intersects) return false; 
   }
   return true;
@@ -255,6 +269,7 @@ void ElementLocator::_exchangeMesh(MEDMEM::MESH* local_mesh, MEDMEM::MESH*& dist
     MEDMEM::MESHING* meshing = new MEDMEM::MESHING ();
     int* recv_buffer_ptr = recv_buffer;
     meshing->setCoordinates(distant_space_dim, distant_nnodes, recv_coords, "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
+		
     meshing->setNumberOfTypes(distant_nb_types,MED_EN::MED_CELL);
 
     // converting the types from int to medGeometryElement
@@ -281,7 +296,8 @@ void ElementLocator::_exchangeMesh(MEDMEM::MESH* local_mesh, MEDMEM::MESH*& dist
     {
       distant_ids_recv[i]=*recv_buffer_ptr++;
     }
-    
+    meshing->setMeshDimension(distant_mesh_dim);
+
     distant_mesh=meshing;  
     delete[] recv_buffer;
     delete[] nbtypes;
@@ -358,13 +374,16 @@ MEDMEM::MESH* ElementLocator::_meshFromElems(set<int>& elems)
   meshing->setNumberOfTypes(nbelems_per_type.size(),MED_EN::MED_CELL);
   meshing->setTypes(new_types,MED_EN::MED_CELL);
   meshing->setNumberOfElements(nbtypes,MED_EN::MED_CELL);
- 
+
+	int dimmax=0;
  int* conn_ptr= conn;
   for (int i=0; i<nbelems_per_type.size(); i++)
   {
     meshing->setConnectivity(conn_ptr, MED_EN::MED_CELL,new_types[i]);
     conn_ptr+=nbtypes[i]*(new_types[i]%100);
+		if (new_types[i]/100>dimmax) dimmax=new_types[i]/100;
   }
+	meshing->setMeshDimension(dimmax);
   delete [] small2big;
   delete [] nbtypes;
   delete [] conn;
diff --git a/src/ParaMEDMEM/ElementLocator.hxx b/src/ParaMEDMEM/ElementLocator.hxx
index af3243c9e..2a9ea9cf3 100644
--- a/src/ParaMEDMEM/ElementLocator.hxx
+++ b/src/ParaMEDMEM/ElementLocator.hxx
@@ -32,6 +32,7 @@ private:
   const ProcessorGroup& _local_group;
   ProcessorGroup* _union_group;
   std::vector<int> _distant_proc_ids;
+	double _adjustment_eps;
   //InterpolationMatrix _matrix;
   //MxN_Mapping _mapping; 
   
diff --git a/src/ParaMEDMEM/INTERPOLATION.hxx b/src/ParaMEDMEM/INTERPOLATION.hxx
new file mode 100644
index 000000000..2ffbea81c
--- /dev/null
+++ b/src/ParaMEDMEM/INTERPOLATION.hxx
@@ -0,0 +1,26 @@
+#ifndef _INTERPOLATION_HXX_
+#define _INTERPOLATION_HXX_
+
+
+#include <vector>
+#include <map>
+
+#include "MEDMEM_Mesh.hxx"
+
+
+namespace ParaMEDMEM
+{
+
+class INTERPOLATION
+{
+public:
+	INTERPOLATION(){}
+	virtual ~INTERPOLATION(){}
+
+	//interpolation of two triangular meshes.
+	virtual std::vector<std::map<int, double> > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)=0;
+
+};
+
+};
+#endif
diff --git a/src/ParaMEDMEM/INTERPOLATION_2D.cxx b/src/ParaMEDMEM/INTERPOLATION_2D.cxx
index 27eadc667..dc23b8f16 100755
--- a/src/ParaMEDMEM/INTERPOLATION_2D.cxx
+++ b/src/ParaMEDMEM/INTERPOLATION_2D.cxx
@@ -1,31 +1,12 @@
+#include "INTERPOLATION_Utils.hxx"
 #include "INTERPOLATION_2D.hxx"
 using namespace std;
 using namespace MED_EN;
 using namespace MEDMEM;
 
 
-
-struct monMaillageS
-{
-    int _NbMaille;
-    int _NbNoeud;
-    const int* _Connect;
-    const double* _Coord;
-    const int* _ReversNConnect;
-    const int*  _ReversNConnectIndex;
-    double _DimCaracteristic;
-};
-
-struct monMaillageP
+namespace ParaMEDMEM 
 {
-    int _NbMaille;
-    int _NbNoeud;
-    const int* _Connect;
-    const double* _Coord;
-    double _DimCaracteristic;
-};
-
-
 
 INTERPOLATION_2D::INTERPOLATION_2D()
 {
@@ -48,428 +29,12 @@ void INTERPOLATION_2D::setOptions(double Precision, int SecondMethodOfLocalisati
     _PrintLevel=PrintLevel;
 }
 
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*   calcul la surface d'un triangle                  */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-double INTERPOLATION_2D::Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
-{
-    double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]);
-    double Surface = 1/2.0*abs(A);
-    return Surface;
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*     fonction qui calcul le déterminant            */
-/*      de deux vecteur(cf doc CGAL).                */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-
-//fonction qui calcul le déterminant des vecteurs: P3P1 et P3P2
-//(cf doc CGAL).
-
-double INTERPOLATION_2D::mon_determinant(const double* P_1,const double*  P_2,
-	const double* P_3)
-{
-    double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
-    return mon_det;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-//calcul la norme du vecteur P1P2
-
-double INTERPOLATION_2D::norme_vecteur(const double* P_1,const double* P_2)
-{
-    double X=P_1[0]-P_2[0];
-    double Y=P_1[1]-P_2[1];
-    double norme=sqrt(X*X+Y*Y);
-    return norme;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*         calcul le cos et le sin de l'angle P1P2,P1P3     */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-vector<double> INTERPOLATION_2D::calcul_cos_et_sin(const double* P_1,const double* P_2,
-	const double* P_3)
-{
-
-    vector<double> Vect;
-    double P1_P2=norme_vecteur(P_1,P_2);
-    double P2_P3=norme_vecteur(P_2,P_3);
-    double P3_P1=norme_vecteur(P_3,P_1);
-
-    double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
-    double D=2.0*P1_P2*P3_P1;
-    double COS=N/D;
-    Vect.push_back(COS);
-    double V=mon_determinant(P_2,P_3,P_1);
-    double D_1=P1_P2*P3_P1;
-    double SIN=V/D_1;
-    Vect.push_back(SIN);
-
-    return Vect;
-
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*fonction pour vérifier qu'un point n'a pas déja été considérer dans   */ 
-/*      le vecteur et le rajouter au vecteur sinon.                     */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-void INTERPOLATION_2D::verif_point_dans_vect(const double* P, vector<double>& V)
-{
-  
-  int taille=V.size();
-  //  bool isPresent=false;
-  for(int i=0;i<taille/2;i++) 
-    {
-      double dx=P[0]-V[2*i];
-      dx=dx>0.0?dx:-dx;
-      if (dx>_Precision)
-	continue;
-      double dy=P[1]-V[2*i+1];
-      dy=dy>0.0?dy:-dy;
-      if (dy<_Precision)
-	
-	return;
-      //	if( sqrt((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))/_DimCaracteristic<_Precision)
-      //  isPresent=true;
-      
-    }
-  //    if(!isPresent)
-  //{
-  
-  V.push_back(P[0]);
-  V.push_back(P[1]);
-  
-  //}
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*     calcul les coordonnées du barycentre d'un polygone   */ 
-/*     le vecteur en entrée est constitué des coordonnées   */
-/*     des sommets du polygone                              */                             
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-vector<double> INTERPOLATION_2D::bary_poly(const vector<double>& V)
-{
-    vector<double> Bary;
-    int taille=V.size();
-    double x=0;
-    double y=0;
-
-    for (int i=0;i<taille/2;i++)
-    {
-	x=x+V[2*i];
-	y=y+V[2*i+1];
-    }
-    double A=2*x/taille;
-    double B=2*y/taille;
-    Bary.push_back(A);//taille vecteur=2*nb de points.
-    Bary.push_back(B);
-
-
-    return Bary;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*         calcul la surface d'un polygone.                 */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-double  INTERPOLATION_2D::Surf_Poly(const vector<double>& Poly)
-{ 
-
-    double Surface=0;
-    for (int i=0; i<(Poly.size())/2-2; i++)
-    {
-	double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] ); 
-	Surface=Surface + Surf ;
-    }
-    return Surface ;
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/
-/*  calcul de l'intersection de deux segments: segments P1P2 avec P3P4      */
-/*  . Si l'intersection est non nulle et si celle-ci n'est                  */
-/*  n'est pas déjà contenue dans Vect on la rajoute à Vect                  */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/ 
-
-void INTERPOLATION_2D::inters_de_segment(const double * P_1,const double * P_2,
-	const double * P_3,const double * P_4,vector<double>&  Vect)
-{
-
-
-    // calcul du déterminant de P1_1P_2 et P_3P_4.
-    double det=(P_2[0]-P_1[0])*(P_4[1]-P_3[1])-(P_4[0]-P_3[0])*(P_2[1]-P_1[1]);
-
-
-    if(abs(det)/_DimCaracteristic>_Precision)
-    {
-
-	double k_1=1/((P_2[0]-P_1[0])*(P_3[1]-P_4[1])+(P_4[0]-P_3[0])*(P_2[1]-P_1[1]))
-	    *((P_3[1]-P_4[1])*(P_3[0]-P_1[0])+(P_4[0]-P_3[0])*(P_3[1]-P_1[1]));
-	
-	if( k_1>=0 &&  k_1<=1)
-	  {
-	    double k_2=1/((P_4[0]-P_3[0])*(P_1[1]-P_2[1])+(P_2[0]-P_1[0])*(P_4[1]-P_3[1]))
-	      *((P_1[1]-P_2[1])*(P_1[0]-P_3[0])+(P_2[0]-P_1[0])*(P_1[1]-P_3[1]));
-	    
-
-	    if( k_2>=0 &&  k_2<=1)
-	    {
-		double P_0[2];
-		P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]);
-		P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]);
-		verif_point_dans_vect(P_0,Vect);
-
-	    }
-	}
-    }
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*   fonction qui teste si un point est dans une maille     */
-/*   point: P_0                                             */
-/*   P_1, P_2, P_3 sommet des mailles                       */   
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-int INTERPOLATION_2D::point_dans_triangle(const double* P_0,const double* P_1,
-	const double* P_2,const double* P_3)
-{
-
-    int A=0;
-    double det_1=mon_determinant(P_1,P_3,P_0);
-    double det_2=mon_determinant(P_3,P_2,P_0);
-    double det_3=mon_determinant(P_2,P_1,P_0);
-    if( (det_1>=0 && det_2>=0 && det_3>=0) || (det_1<=0 && det_2<=0 && det_3<=0) )
-    {
-	A=1;
-    }
-
-    return A;
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/* fonction qui rajoute les sommet du triangle P dans le vecteur V        */ 
-/* si ceux-ci sont compris dans le triangle S et ne sont pas déjà dans    */
-/* V.                                                                     */
-/*sommets de P: P_1, P_2, P_3                                             */
-/*sommets de S: P_4, P_5, P_6                                             */  
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ 
-
-void INTERPOLATION_2D::rajou_sommet_triangl(const double* P_1,const double* P_2,const double* P_3,
-	const double* P_4,const double* P_5,const double* P_6,vector<double>& V)
-{
-
-    int A_1=point_dans_triangle(P_1,P_4,P_5,P_6);
-    if(A_1==1)
-    {verif_point_dans_vect(P_1,V);
-  //   if (V.size()==1)
-//       {
-// 	// all nodes are necessarily in the triangle
-// 	verif_point_dans_vect(P_2,V);
-// 	verif_point_dans_vect(P_3,V);
-// 	return;
-//       }
-    }
-    int A_2=point_dans_triangle(P_2,P_4,P_5,P_6);
-    if(A_2==1)
-    {verif_point_dans_vect(P_2,V);}
-
-    int A_3=point_dans_triangle(P_3,P_4,P_5,P_6);
-    if(A_3==1)
-    {verif_point_dans_vect(P_3,V);}
-    
-	
-      
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
-/*      calcul l'intersection de deux triangles            */
-/* P_1, P_2, P_3: sommets du premier triangle              */
-/* P_4, P_5, P_6: sommets du deuxième triangle             */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
-
-vector<double> INTERPOLATION_2D::intersec_de_triangle(const double* P_1,const double* P_2,
-	const double* P_3,const double* P_4,const double* P_5,const double* P_6)
-{
-
-    //debug cout << "  T1 : " << *P_1 << " , " << *(P_1+1) << " , " << *P_2 << " , " << *(P_2+1) << " , " << *P_3 << " , " << *(P_3+1)<< endl;
-    //debug cout << "  T2 : " << *P_4 << " , " << *(P_4+1) << " , " << *P_5 << " , " << *(P_5+1) << " , " << *P_6 << " , " << *(P_6+1)<< endl;
-    vector<double> Vect;
-
-    inters_de_segment(P_1,P_2,P_4,P_5,Vect);
-    inters_de_segment(P_1,P_2,P_5,P_6,Vect);
-    inters_de_segment(P_1,P_2,P_6,P_4,Vect);
-    inters_de_segment(P_2,P_3,P_4,P_5,Vect);
-    inters_de_segment(P_2,P_3,P_5,P_6,Vect);
-    inters_de_segment(P_2,P_3,P_6,P_4,Vect);
-    inters_de_segment(P_3,P_1,P_4,P_5,Vect);
-    inters_de_segment(P_3,P_1,P_5,P_6,Vect);
-    inters_de_segment(P_3,P_1,P_6,P_4,Vect);
-    rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect);
-    rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect);
-    //debug cout << "  Inter : ";
-    //debug for (int i=0; i<Vect.size(); ++i)
-	//debug cout << Vect[i] << "  ";
-    //debug cout << endl << endl;
-
-    return Vect;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */  
-/* fonction pour reconstituer un polygone convexe à partir  */
-/*              d'un nuage de point.                        */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */  
-
-vector<double> INTERPOLATION_2D::reconstruct_polygon(const vector<double>& V)
-{
-
-    int taille=V.size();
-    if(taille<=6)
-    {return V;}
-    else
-    {
-	double COS[taille/2];
-	double SIN[taille/2];
-	double angle[taille/2];
-	vector<double> Bary=bary_poly(V);
-	COS[0]=1.0;
-	SIN[0]=0.0;
-	angle[0]=0.0;
-	for(int i=0; i<taille/2-1;i++)
-	{
-	    vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
-	    COS[i+1]=Trigo[0];
-	    SIN[i+1]=Trigo[1];
-	    if(SIN[i+1]>=0)
-	    {angle[i+1]=acos(COS[i+1]);}
-	    else
-	    {angle[i+1]=-acos(COS[i+1]);}
-	}
-
-	//ensuite on ordonne les angles.
-	vector<double> Pt_ordonne;
-	Pt_ordonne.reserve(taille);
-	map<double,int> Ordre;
-	for(int i=0;i<taille/2;i++)	
-	{Ordre[angle[i]]=i;}
-	map <double,int>::iterator mi;
-	for(mi=Ordre.begin();mi!=Ordre.end();mi++)
-	{
-	    int j=(*mi).second;
-	    Pt_ordonne.push_back(V[2*j]);
-	    Pt_ordonne.push_back(V[2*j+1]);
-	}
-	return Pt_ordonne;
-    }
-}
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */ 
-/*fonction pour trouver les mailles voisines d'une maille triangle.*/
-/* (mailles voisines par les faces)                                */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */ 
-
-vector<int> INTERPOLATION_2D::touv_maill_vois(int i_S,const monMaillageS& MaS)
-{
-
-    int N_1=MaS._Connect[3*(i_S-1)];
-    int N_2=MaS._Connect[3*(i_S-1)+1];
-    int N_3=MaS._Connect[3*(i_S-1)+2];
-    vector<int> tab;
-    int I_1=MaS._ReversNConnectIndex[N_1-1]-1;// Dans MED le n°des noeuds commencent à 1. 
-    int I_2=MaS._ReversNConnectIndex[N_2-1]-1;
-    int I_3=MaS._ReversNConnectIndex[N_3-1]-1;
-    int F_1=MaS._ReversNConnectIndex[N_1]-1;
-    int F_2=MaS._ReversNConnectIndex[N_2]-1;
-    int F_3=MaS._ReversNConnectIndex[N_3]-1;
-
-
-    vector<int> V_12;
-    V_12.reserve(2);
-    insert_iterator<vector<int> > ib_1(V_12,V_12.begin());
-    vector<int> V_23;
-    V_23.reserve(2);
-    insert_iterator<vector<int> > ib_2(V_23,V_23.begin());
-    vector<int> V_13;
-    V_13.reserve(2);
-    insert_iterator<vector<int> > ib_3(V_13,V_13.begin());
-
-
-    set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_2,
-	    MaS._ReversNConnect+F_2,ib_1);
-    set_intersection(MaS._ReversNConnect+I_2,MaS._ReversNConnect+F_2,MaS._ReversNConnect+I_3,
-	    MaS._ReversNConnect+F_3,ib_2);
-    set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_3,
-	    MaS._ReversNConnect+F_3,ib_3);
-
-
-    for(int i=0;i<(V_12).size();i++)
-    {
-	if(V_12[i]!=i_S)
-	{
-	    tab.push_back(V_12[i]);
-	    break;
-	}
-    }
-
-    for(int i=0;i<V_23.size();i++)
-    {
-	if(V_23[i]!=i_S)
-	{
-	    tab.push_back(V_23[i]);
-	    break;
-	}
-    }
-
-    for(int i=0;i<V_13.size();i++)
-    {
-	if(V_13[i]!=i_S)
-	{
-	    tab.push_back(V_13[i]);
-	    break;
-	}
-    }
-    return tab;
-}
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-/* fonction pour vérifier qu'un n°de maille n'a pas déja été considérer  */
-/*  dans le vecteur et le rajouter au vecteur sinon.                     */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
-
-void INTERPOLATION_2D::verif_maill_dans_vect(int Num, vector<int>& V)
-{
-    int taille=V.size();
-    int A=0;
-    for(int i=0;i<taille;i++)
-    {
-	if(Num==V[i])
-	{
-	    A=1;
-	    break;
-	} 
-    }
-    if(A==0)
-    {V.push_back(Num); }
-}
-
-
-
 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 /*localise le barycentre des mailles de P dans le maillage S*/
 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ 
 
-void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,monMaillageP& MaP,
-	monMaillageS& MaS,vector<int>& localis )
+void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,ParaMEDMEM_utils::monMaillageP& MaP,
+																			 ParaMEDMEM_utils::monMaillageS& MaS,vector<int>& localis )
 {
 
     //calcul des coordonnées des barycentre des mailles de P
@@ -510,7 +75,7 @@ void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,monMaillage
 		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);
+		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;
@@ -524,9 +89,9 @@ void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,monMaillage
 
 
 		    //debug cout << "mailles sources testées : " << N_Maille << endl;
-		    vector<double> Inter=intersec_de_triangle(MaS._Coord+2*(NS_1-1),MaS._Coord+2*(NS_2-1),
+		    vector<double> Inter=ParaMEDMEM_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));
+																																		MaS._Coord+2*(NS_3-1),MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),MaP._Coord+2*(NP_3-1), _DimCaracteristic, _Precision);
 
 		    if(Inter.size()>0)
 		    {
@@ -560,8 +125,8 @@ void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,monMaillage
 			    int N_3=MaS._Connect[3*(N_Maille-1)+2];
 
 			    //debug cout << "mailles sources testées : " << N_Maille << endl;
-			    vector<double> Inter=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) );	
+			    vector<double> Inter=ParaMEDMEM_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) , _DimCaracteristic, _Precision);	
 
 			    if(Inter.size()>0)
 			    {
@@ -594,7 +159,7 @@ void INTERPOLATION_2D::local_iP_dans_S(MESH& myMesh_S,MESH& myMesh_P,monMaillage
 /* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ _ _ _ _ _ _ _ _ _ _  */
 
 
-inline void INTERPOLATION_2D::rempli_vect_d_intersect(int i_P,int i_S1,const monMaillageP& MaP,const monMaillageS& MaS,	
+inline void INTERPOLATION_2D::rempli_vect_d_intersect(int i_P,int i_S1,const ParaMEDMEM_utils::monMaillageP& MaP,const ParaMEDMEM_utils::monMaillageS& MaS,	
 	vector<map<int,double> >& Surface_d_intersect,FIELD<double>& myField_P)
 {
 
@@ -609,7 +174,7 @@ inline void INTERPOLATION_2D::rempli_vect_d_intersect(int i_P,int i_S1,const mon
 
     //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),
+    double Surf_i_P =ParaMEDMEM_utils::Surf_Tri(MaP._Coord+2*(NP_1-1),MaP._Coord+2*(NP_2-1),
 	    MaP._Coord+2*(NP_3-1));
 
     double Surface = 0 ;
@@ -628,61 +193,56 @@ inline void INTERPOLATION_2D::rempli_vect_d_intersect(int i_P,int i_S1,const mon
 	    int NS_3=MaS._Connect[3*(i_S-1)+2];
 
 
-	    vector<double> Inter=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));
-
-
-
+	    vector<double> Inter=ParaMEDMEM_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), _DimCaracteristic, _Precision);
 
 	    //on teste l'intersection.
 	    int taille=Inter.size()/2;
 	    //debug cout << "  -> maille source : " << i_S << " , nb intersection=" <<  taille << endl;
 
 	    /* CAS1  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-	    if(taille==0)
-	    {int Rien=0;}
+	    // taille ==0, on ne fait rien
 
 	    /* CAS 2  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-	    else if(taille==1 || taille==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=touv_maill_vois(i_S,MaS);
+				vector<int> M_Vois=ParaMEDMEM_utils::touv_maill_vois(i_S,MaS);
 		for(int i=0;i< M_Vois.size();i++)
-		{verif_maill_dans_vect(M_Vois[i],maill_deja_vu_S);}
+			{ParaMEDMEM_utils::verif_maill_dans_vect(M_Vois[i],maill_deja_vu_S);}
 	    }
 
 	    /* CAS 3  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
 	    else if(taille==3)
 	    {
-		double Surf_inter = Surf_Poly(Inter) ;
+				double Surf_inter = ParaMEDMEM_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=touv_maill_vois(i_S,MaS);
-		for(int i_M=0;i_M<(M_Vois).size();i_M++)
-		{verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
-
-		Surface = Surface + Surf_inter ; 
+				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=ParaMEDMEM_utils::touv_maill_vois(i_S,MaS);
+				for(int i_M=0;i_M<(M_Vois).size();i_M++)
+					{ParaMEDMEM_utils::verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
+
+				Surface = Surface + Surf_inter ; 
 	    }
-
+			
 	    /* CAS 4  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 
-	    else //taille>3
+	    else if (taille > 3) //taille>3
 	    {
-		vector<double> Poly_Inter=reconstruct_polygon(Inter);
+				vector<double> Poly_Inter=ParaMEDMEM_utils::reconstruct_polygon(Inter);
 
-		double Surf_inter =Surf_Poly(Poly_Inter) ;
+				double Surf_inter =ParaMEDMEM_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=touv_maill_vois(i_S,MaS);
+		vector<int> M_Vois=ParaMEDMEM_utils::touv_maill_vois(i_S,MaS);
 		for(int i_M=0;i_M<(M_Vois).size();i_M++)
-		{verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
+		{ParaMEDMEM_utils::verif_maill_dans_vect(M_Vois[i_M],maill_deja_vu_S);}
 
 		Surface = Surface + Surf_inter ; 
 	    }
@@ -720,8 +280,8 @@ vector<map<int,double> > INTERPOLATION_2D::interpol_maillages(const MEDMEM::MESH
     delete[] Type_P;
 
     /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-    monMaillageS MaS;
-    monMaillageP MaP;
+		ParaMEDMEM_utils::monMaillageS MaS;
+		ParaMEDMEM_utils::monMaillageP MaP;
 
     MaS._NbNoeud = myMesh_S.getNumberOfNodes() ;
     MaP._NbNoeud = myMesh_P.getNumberOfNodes() ;
@@ -748,10 +308,6 @@ vector<map<int,double> > INTERPOLATION_2D::interpol_maillages(const MEDMEM::MESH
 
     /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 
-
-    int SpaceDimension_S = myMesh_S.getNumberOfNodes() ;
-    int SpaceDimension_P = myMesh_P.getNumberOfNodes() ;
-
     /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
     //on cherche la dimension caractéristique des maillages
     vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
@@ -851,29 +407,28 @@ vector<map<int,double> > INTERPOLATION_2D::interpol_maillages(const MEDMEM::MESH
 //    }
 
     /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+  //   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;
 
 }
 
 
-
-
-
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-/*  fonction qui écrit les résultats d'annelise dans un fichier:      */
-/*  pour chaque maille du maillage 1 on stoque la fraction volumique  */
-/*          considéré lors de l'algorithme                            */
-/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
-
-void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
-	MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S)
-{
-    //On sauvegarde le champ crée sur la maillage P.
-    string FileName="FracSurf"+NameMesh_S+"_to_"+NameMesh_P;
-    int id1=(*myMesh_P).addDriver(MED_DRIVER,FileName,"MP");
-    int id2=(*myMesh_S).addDriver(MED_DRIVER,FileName,"MS");
-    int id3=FractSurf.addDriver(MED_DRIVER,FileName,"% MP");
-    (*myMesh_P).write(id1);
-    (*myMesh_S).write(id2);
-    FractSurf.write(id3);
-}
+};
diff --git a/src/ParaMEDMEM/INTERPOLATION_2D.hxx b/src/ParaMEDMEM/INTERPOLATION_2D.hxx
index 58026ff48..941c24979 100755
--- a/src/ParaMEDMEM/INTERPOLATION_2D.hxx
+++ b/src/ParaMEDMEM/INTERPOLATION_2D.hxx
@@ -1,7 +1,6 @@
 #ifndef _INTERPOLATION_2D_HXX_
 #define _INTERPOLATION_2D_HXX_
 
-
 #include "MEDMEM_InterpolationHighLevelObjects.hxx"
 #include "MEDMEM_Mesh.hxx"
 #include "MEDMEM_Field.hxx"
@@ -14,6 +13,7 @@
 #include "MEDMEM_DriverFactory.hxx"
 #include "MEDMEM_Meshing.hxx"
 #include "MEDMEM_define.hxx"
+#include "INTERPOLATION.hxx"
 #include "MEDMEM_Interpolation.hxx"
 
 
@@ -27,97 +27,43 @@
 #include <complex>
 
 
-struct monMaillageS;
-
+namespace ParaMEDMEM_utils
+{
+	struct monMaillageS;
+	struct monMaillageP;
+};
 
-struct monMaillageP;
+namespace ParaMEDMEM
+{
 
+	class INTERPOLATION;
 
-class INTERPOLATION_2D
+class INTERPOLATION_2D : public INTERPOLATION
 {
 
-
     private: 
 	double _Precision;
 	double _DimCaracteristic;
 	int  _SecondMethodOfLocalisation;
 	int  _PrintLevel;
+
 	// Méthodes publiques
     public:
 
 	INTERPOLATION_2D();
 
-
 	// precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
 	void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
 
-	//pour calculer la surface d'un triangle.
-	double Surf_Tri(const double* P_1,const double* P_2,
-		const double* P_3);
-
-	//calcul déterminant de trois points (cf doc CGAL)
-	double mon_determinant(const double* P_1,const double*  P_2,
-		const double* P_3);
-
-	//calcul la norme d'un vecteur.
-	double norme_vecteur(const double* P_1,const double* P_2);
-
-	//calcul le cos et le sin.
-	vector<double> calcul_cos_et_sin(const double* P_1,
-		const double* P_2,const double* P_3);
-
-	//fonction pour vérifier qu'un point n'a pas déja été considérer dans 
-	//le vecteur et le rajouter au vecteur sinon.
-	void verif_point_dans_vect(const double* P, vector<double>& V);
-
-	//calcul les coordonnées du barycentre
-	//d'un polygone
-	vector<double> bary_poly(const vector<double>& V);
-
-	//fonction qui calcul la surface d'un polygone
-	double Surf_Poly(const vector<double>& Poly);
-
-	//calcul de l'intersection de deux segments.
-	void inters_de_segment(const double* P_1,const double* P_2,
-		const double* P_3,const double* P_4,vector<double>& Vect);
-
-	//fonction qui teste si un point est dans une maille
-	int point_dans_triangle(const double* P_0,const double* P_1,
-		const double* P_2,const double* P_3);
-
-	//fonction qui rajouter les sommets du triangle P2P2P3 si ceux-ci sont compris
-	// dans le triangle P4P5P6 et n'ont pas encore été considérer.
-	//Rq:point du triangle donné dans ordre quelconque.
-	void rajou_sommet_triangl(const double* P_1,const double* P_2,
-		const double* P_3,const double* P_4,const double* P_5,
-		const double* P_6,vector<double>& V);
-
-	//calcul l'intersection de deux triangles (sommets du triangle donné dans
-	//ordre quelconque)
-	vector<double> intersec_de_triangle(const double* P_1,const double* P_2,
-		const double* P_3,const double * P_4,const double * P_5,
-		const double * P_6);  
-
-	//fonction pour reconstituer un polygon convexe à partir
-	//d'un nuage de point.
-	vector<double> reconstruct_polygon(const vector<double>& V);
-
-	//fonction pour trouver les mailles voisines d'une maille triangle.
-	vector<int> touv_maill_vois(int i_S,const monMaillageS& MaS);
-
-	//fonction pour vérifier qu'un n°de maille n'a pas déja été considérer
-	// dans le vecteur et le rajouter au vecteur sinon.
-	void INTERPOLATION_2D::verif_maill_dans_vect(int Num, vector<int>& V);
-
 
 	//localise le barycentre des mailles de P dans le maillage S
 	void local_iP_dans_S(MEDMEM::MESH& myMesh_S,MEDMEM::MESH& myMesh_P,
-		monMaillageP& MaP,monMaillageS& MaS,vector<int>& localis );
+		ParaMEDMEM_utils::monMaillageP& MaP,ParaMEDMEM_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 monMaillageP& MaP,const monMaillageS& MaS,
+	inline void rempli_vect_d_intersect(int i_P,int i_S,const ParaMEDMEM_utils::monMaillageP& MaP,const ParaMEDMEM_utils::monMaillageS& MaS,
 		vector<map<int,double> >& Surface_d_intersect,
 		FIELD<double>& myField_P);
 
@@ -125,26 +71,11 @@ class INTERPOLATION_2D
 	//fonction principale pour interpoler deux maillages triangulaires.
 	std::vector<map<int, double> > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
 
-
-	//  fonction qui écrit les résultats d'annelise dans un fichier:      
-	//  pour chaque maille du maillage 1 on stoque la fraction volumique  
-	//          considéré lors de l'algorithme                            
-
-
-
-
-
-	void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
-		MEDMEM::MESH* myMesh_S,string FileName,string NameMesh_P,string NameMesh_S);
-
-
-
-	// MEDMEM::FIELD<double>* createField();
-	// ...
-
     private: 
 
 
+};
+
 };
 
 #endif
diff --git a/src/ParaMEDMEM/INTERPOLATION_3D_surf.cxx b/src/ParaMEDMEM/INTERPOLATION_3D_surf.cxx
new file mode 100644
index 000000000..cf897a169
--- /dev/null
+++ b/src/ParaMEDMEM/INTERPOLATION_3D_surf.cxx
@@ -0,0 +1,395 @@
+#include "MEDMEM_Field.hxx"
+#include "MEDMEM_Mesh.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "INTERPOLATION_Utils.hxx"
+#include "INTERPOLATION_3D_surf.hxx"
+#include "BBTree.H"
+
+using namespace std;
+using namespace MED_EN;
+using namespace MEDMEM;
+using namespace ParaMEDMEM_utils;
+
+namespace ParaMEDMEM
+{
+
+	INTERPOLATION_3D_surf::INTERPOLATION_3D_surf()
+	{
+    _Precision=1.0E-12;
+    _DimCaracteristic=1;
+    _SecondMethodOfLocalisation=1;
+    _PrintLevel=0;
+		_Surf3DAdjustmentEps=0.1;
+	}
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*   Options :                                        */
+	/*     Precision : for geometric computation          */
+	/*     SecondMethodOfLocalisation : 0/1               */
+	/*     PrintLevel : between 0 and 3                   */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	void INTERPOLATION_3D_surf::setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel)
+	{
+    _Precision=Precision;
+    _SecondMethodOfLocalisation=SecondMethodOfLocalisation;
+    _PrintLevel=PrintLevel;
+	}
+
+
+
+	/*! method computing the translation/rotation matrix
+		necessary to transform a triangle into a triangle located inside the Oxy plane. The method fails only if two of the three points are coincident
+		\param P1 pointer to three coordinates defining the first point
+		\param P2 pointer to three coordinates defining the second point
+		\param P3 pointer to three coordinates defining the third point
+		\return TranslationRotationMatrix object containing the transform that must be applied to the triangle for putting it inside the Oxy plane.
+	*/
+	void INTERPOLATION_3D_surf::rotate3DTriangle(const double* PP1, const double*PP2, const double*PP3, TranslationRotationMatrix& rotation_matrix)
+	{
+		//initializes
+		rotation_matrix.translate(PP1);
+  
+		double P2w[3];
+		double P3w[3];
+		P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
+		P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
+
+		// translating to set P1 at the origin
+		for (int i=0; i<2; i++)
+			{
+				P2w[i]-=PP1[i];
+				P3w[i]-=PP1[i];
+			}
+   
+		// rotating to set P2 on the Oxy plane
+		TranslationRotationMatrix A;
+		A.rotate_x(P2w);
+		A.rotate_vector(P3w);
+		rotation_matrix.multiply(A);
+
+		//rotating to set P2 on the Ox axis
+		TranslationRotationMatrix B;
+		B.rotate_z(P2w);
+		B.rotate_vector(P3w);
+		rotation_matrix.multiply(B);
+  
+		//rotating to set P3 on the Oxy plane
+		TranslationRotationMatrix C;
+		C.rotate_x(P3w);
+		rotation_matrix.multiply(C);
+  
+	}
+
+
+
+	/* _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ */
+	/*  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 INTERPOLATION_3D_surf::rempli_vect_d_intersect(int i_P,int i_S,const monMaillageP& MaP,const monMaillageS& MaS,	
+																														 vector<map<int,double> >& Surface_d_intersect
+																														 )
+ 
+	{
+
+
+    //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];
+
+    
+    TranslationRotationMatrix rotation;
+		rotate3DTriangle(MaP._Coord+3*(NP_1-1),MaP._Coord+3*(NP_2-1), MaP._Coord+3*(NP_3-1), rotation);
+    
+    //on calcule la surface de la maille i_P
+
+		//    double Surf_i_P =ParaMEDMEM_utils::Surf_Tri(MaP._Coord+3*(NP_1-1),MaP._Coord+3*(NP_2-1),
+		//																						MaP._Coord+3*(NP_3-1));
+
+		// double Surface = 0 ;
+ 
+    
+    
+    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];
+    
+    double PS1[3];
+    double PS2[3];
+    double PS3[3];
+		double PP1[3];
+		double PP2[3];
+		double PP3[3];
+		
+    for (int i=0; i<3; i++)
+      {
+				PS1[i] = *(MaS._Coord+3*(NS_1-1)+i);
+				PS2[i] = *(MaS._Coord+3*(NS_2-1)+i);
+				PS3[i] = *(MaS._Coord+3*(NS_3-1)+i);
+				PP1[i] = *(MaP._Coord+3*(NP_1-1)+i);
+				PP2[i] = *(MaP._Coord+3*(NP_2-1)+i);
+				PP3[i] = *(MaP._Coord+3*(NP_3-1)+i);	
+      }
+
+    rotation.transform_vector(PS1);
+    rotation.transform_vector(PS2);
+    rotation.transform_vector(PS3);
+    rotation.transform_vector(PP1);
+    rotation.transform_vector(PP2);
+    rotation.transform_vector(PP3);
+
+
+
+    //intersects 3D triangle only considering the two first coordinates
+    // 
+    //> MaP is located in the Oxy plane
+    //> PS is not and is therefore projected on MaP orthogonally along the z
+    //   direction 
+    vector<double> Inter=ParaMEDMEM_utils::
+			intersec_de_triangle(PS1,PS2,PS3,
+													 PP1,PP2,PP3,
+													 _DimCaracteristic,
+													 _Precision);
+    
+    
+    //on teste l'intersection.
+    int taille=Inter.size()/2;
+    //debug cout << "  -> maille source : " << i_S << " , nb intersection=" <<  taille << endl;
+    
+    
+    /* CAS 3  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
+    if(taille==3)
+      {
+				double Surf_inter = Surf_Poly(Inter) ;
+	
+				//on remplit le vecteur des correspondances n° maille-surfaces d'intersections.
+				Surface_d_intersect[i_P-1][i_S]=Surf_inter;
+	
+	
+				//	Surface = Surface + Surf_inter ; 
+      }
+    
+    /* CAS 4  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+    
+    else if (taille>3) //taille>3
+      {
+				vector<double> Poly_Inter=reconstruct_polygon(Inter);
+	
+				double Surf_inter =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;
+	
+				//	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.
+	vector<map<int,double> > INTERPOLATION_3D_surf::interpol_maillages(const MEDMEM::MESH& myMesh_S,
+																																		 const MEDMEM::MESH& myMesh_P)
+	{
+
+
+    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+    // 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;}
+    //VB
+    delete[] Type_S;
+    delete[] Type_P;
+
+    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+		ParaMEDMEM_utils::monMaillageS MaS;
+		ParaMEDMEM_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 cherche la dimension caractéristique des maillages
+		/*______________________________________________________________*/
+
+    vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+    vector<vector<double> > BoxP=myMesh_P.getBoundingBox();
+    double VolumeS=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]))+(BoxS[1][2]-BoxS[0][2])*(BoxS[1][2]-BoxS[0][2]);
+    MaS._DimCaracteristic=sqrt(VolumeS/MaS._NbMaille);
+    double VolumeP=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]))+(BoxP[1][2]-BoxP[0][2])*(BoxP[1][2]-BoxP[0][2]);
+    MaP._DimCaracteristic=sqrt(VolumeP/MaP._NbMaille);
+    
+    _DimCaracteristic=min(MaS._DimCaracteristic,MaP._DimCaracteristic);
+    if (_PrintLevel)
+			{
+				cout<<"recherche de la dimension caractéristique des maillages 3D_surf :"<<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;
+			}
+
+		//creating a search structure based on the bounding boxes 
+		//of the elements to know 
+
+		vector<double> bbox;
+		createBoundingBoxes(MaS,bbox); // create the bounding boxes
+		adjustBoundingBoxes(bbox); // expand them so that no intersection is missed
+ 
+		BBTree<3> my_tree(&bbox[0], 0, 0,MaS._NbMaille );//creating the search structure 
+
+    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+    //on crée un tableau où l'on rentrera la valeurs des intersections.
+    cout << "INTERPOLATION_3D_surf::calcul intersec"<<endl;
+    map<int,double> MAP_init;
+    vector<map<int,double> > Surface_d_intersect(MaP._NbMaille,MAP_init);//on initialise.
+
+
+    /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+    //boucle sur les mailles de P.
+    //Coeur de l'algorithme
+
+    for(int i_P=0; i_P<MaP._NbMaille ; i_P++)
+			{
+				vector<int> intersecting_elems;
+				double bb[6];
+				ParaMEDMEM_utils::getElemBB(bb,MaP,i_P);
+				my_tree.getIntersectingElems(bb, intersecting_elems);
+				
+				//browsing all the i_S (from mesh S) elems that can 
+				//intersect elem i_P (from mesh P)
+				for (int ielem=0; ielem<intersecting_elems.size();ielem++)
+					{
+						//BBTree structure returns numbers between 0 and n-1
+						int i_S=intersecting_elems[ielem]+1; 
+	
+						rempli_vect_d_intersect(i_P+1,i_S,MaP,MaS,Surface_d_intersect);
+					}
+				intersecting_elems.clear();
+			}
+
+
+		//DEBUG prints
+    if (_PrintLevel >= 2)
+			{
+				cout<<endl<<"Impression des surfaces d'intersection:"<<endl<<endl;
+				cout<<"(maille source, maille cible):surface d'intersection"<<endl;
+				double total=0.0;
+				for(int i=0; i< Surface_d_intersect.size();i++)
+					{ 
+						if(Surface_d_intersect[i].size()!=0)
+							{
+								map<int,double>::iterator surface;
+								for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++)
+									{
+										cout<<"    ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+										total+=(*surface).second;
+									}
+							}
+					}  
+				cout << "surface totale "<<total<<endl;
+			}
+
+    return Surface_d_intersect;
+
+	}
+
+
+
+
+
+	/*! Creates a set of bounding boxes that correspond to mesh \a mesh. The vector returned is of dimension 6*nb_elems with bounding boxes stored as xmin1, xmax1, ymin1, ymax1, zmin1, zmax1, xmin2, xmax2, ymin2,... 
+		The returned pointer must be deleted by the calling code.
+
+		\param mesh structure pointing to the mesh
+		\param bbox vector containing the bounding boxes
+	*/
+
+	void INTERPOLATION_3D_surf::createBoundingBoxes(const struct monMaillageS Ma_S, vector<double>& bbox)
+	{
+		/* We build the segment tree for locating possible matching intersections*/
+		bbox.resize(6* Ma_S._NbMaille);
+		for (int i=0; i< Ma_S._NbMaille; i++)
+			{
+				bbox[i*6]=HUGE;
+				bbox[i*6+1]=-HUGE;
+				bbox[i*6+2]=HUGE;
+				bbox[i*6+3]=-HUGE;
+				bbox[i*6+4]=HUGE;
+				bbox[i*6+5]=-HUGE;
+			
+				for (int j=0; j<3; j++)
+					{
+						double x=Ma_S._Coord[(Ma_S._Connect[i*3+j]-1)*3];
+						bbox[i*6]=(bbox[i*6]<x)?bbox[i*6]:x;
+						bbox[i*6+1]=(bbox[i*6+1]>x)?bbox[i*6+1]:x;
+						double y=Ma_S._Coord[(Ma_S._Connect[i*3+j]-1)*3+1];
+						bbox[i*6+2]=(bbox[i*6+2]<y)?bbox[i*6+2]:y;
+						bbox[i*6+3]=(bbox[i*6+3]>y)?bbox[i*6+3]:y;
+						double z=Ma_S._Coord[(Ma_S._Connect[i*3+j]-1)*3+2];
+						bbox[i*6+4]=(bbox[i*6+4]<z)?bbox[i*6+4]:z;
+						bbox[i*6+5]=(bbox[i*6+5]>z)?bbox[i*6+5]:z;
+					}
+			}
+	}
+
+
+	/*! Readjusts a set of bounding boxes so that they are extended
+		in all dimensions for avoiding missing interesting intersections
+
+		\param bbox vector containing the bounding boxes
+	*/
+
+	void INTERPOLATION_3D_surf::adjustBoundingBoxes(vector<double>& bbox)
+	{
+		/* We build the segment tree for locating possible matching intersections*/
+
+		int size = bbox.size()/6;
+		for (int i=0; i<size; i++)
+			{
+				double Dx=bbox[i*6+1]-bbox[i*6];
+				double Dy=bbox[i*6+3]-bbox[i*6+2];
+				double Dz=bbox[i*6+5]-bbox[i*6+4];
+				double max=(Dx<Dy)?Dy:Dx;
+				max=(max<Dz)?Dz:max;
+				bbox[i*6]-=_Surf3DAdjustmentEps*max;
+				bbox[i*6+1]+=_Surf3DAdjustmentEps*max;
+				bbox[i*6+2]-=_Surf3DAdjustmentEps*max;
+				bbox[i*6+3]+=_Surf3DAdjustmentEps*max;
+				bbox[i*6+4]-=_Surf3DAdjustmentEps*max;
+				bbox[i*6+5]+=_Surf3DAdjustmentEps*max;
+
+			}
+	}
+
+};
diff --git a/src/ParaMEDMEM/INTERPOLATION_3D_surf.hxx b/src/ParaMEDMEM/INTERPOLATION_3D_surf.hxx
new file mode 100644
index 000000000..b2a279065
--- /dev/null
+++ b/src/ParaMEDMEM/INTERPOLATION_3D_surf.hxx
@@ -0,0 +1,71 @@
+#ifndef _INTERPOLATION_3D_SURF_HXX_
+#define _INTERPOLATION_3D_SURF_HXX_
+
+
+
+#include "MEDMEM_Mesh.hxx"
+#include "INTERPOLATION_Utils.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "INTERPOLATION.hxx"
+
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <string>
+#include <set>
+#include <algorithm>
+#include <vector>
+
+namespace ParaMEDMEM
+{
+	
+	class INTERPOLATION_3D_surf : public INTERPOLATION
+{
+
+
+    private: 
+	double _Precision;
+	double _DimCaracteristic;
+	int  _SecondMethodOfLocalisation;
+	int  _PrintLevel;
+	double _Surf3DAdjustmentEps;
+	// Méthodes publiques
+    public:
+
+	INTERPOLATION_3D_surf();
+
+
+	// precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
+	void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
+
+
+
+	//fonction principale pour interpoler deux maillages triangulaires.
+	std::vector<map<int, double> > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+
+
+    private: 
+
+	//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 ParaMEDMEM_utils::monMaillageP& MaP,
+																			const ParaMEDMEM_utils::monMaillageS& MaS,
+																			vector<map<int,double> >& Surface_d_intersect);
+
+
+	void createBoundingBoxes(const struct ParaMEDMEM_utils::monMaillageS mesh, vector<double>& bbox);
+
+	void  rotate3DTriangle(const double* PP1, 
+											 const double* PP2,
+											 const double* PP3,
+											 TranslationRotationMatrix& matrix);
+
+	void adjustBoundingBoxes(vector<double>& bbox);
+
+
+};
+
+};
+
+#endif
diff --git a/src/ParaMEDMEM/INTERPOLATION_Utils.hxx b/src/ParaMEDMEM/INTERPOLATION_Utils.hxx
new file mode 100644
index 000000000..fd2a66996
--- /dev/null
+++ b/src/ParaMEDMEM/INTERPOLATION_Utils.hxx
@@ -0,0 +1,553 @@
+#ifndef _INTERPOLATION_UTILS_HXX_
+#define _INTERPOLATION_UTILS_HXX_
+
+
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Field.hxx"
+
+
+#include <cmath>
+#include <map>
+#include <algorithm>
+#include <vector>
+
+
+namespace ParaMEDMEM_utils
+{
+	//	const double HUGE=1e300;
+	struct monMaillageS
+	{
+    int _NbMaille;
+    int _NbNoeud;
+    const int* _Connect;
+    const double* _Coord;
+    const int* _ReversNConnect;
+    const int*  _ReversNConnectIndex;
+    double _DimCaracteristic;
+	};
+
+	struct monMaillageP
+	{
+    int _NbMaille;
+    int _NbNoeud;
+    const int* _Connect;
+    const double* _Coord;
+    double _DimCaracteristic;
+	};
+
+// 	void verif_point_dans_vect(const double* P, vector<double>& V,double dim_caracteristic, double precision );
+// 	inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3);
+// 	inline double norme_vecteur(const double* P_1,const double* P_2);
+// 	inline double mon_determinant(const double* P_1,
+// 																const double*  P_2,
+// 																const double* P_3);
+// 	vector<double> calcul_cos_et_sin(const double* P_1,
+// 																	 const double* P_2,
+// 																	 const double* P_3);
+// 	inline vector<double> bary_poly(const vector<double>& V);
+// 	bool point_dans_triangle(const double* P_0,const double* P_1,
+// 													 const double* P_2,const double* P_3);
+// 	double  Surf_Poly(const vector<double>& Poly);
+// 	vector<double> reconstruct_polygon(const vector<double>& V);
+	
+// 	void rajou_sommet_triangl(const double* P_1,
+// 														const double* P_2,
+// 														const double* P_3,
+// 														const double* P_4,
+// 														const double* P_5,
+// 														const double* P_6,
+// 														vector<double>& V, 
+// 														double dim_caracteristic, 
+// 														double precision);
+	
+	
+// 	inline void inters_de_segment(const double * P_1,
+// 																const double * P_2,
+// 																const double * P_3,
+// 																const double * P_4,
+// 																vector<double>&  Vect, 
+// 																double dim_caracteristic, 
+// 																double precision);
+// 	vector<double> intersec_de_triangle(const double* P_1,const double* P_2,
+// 																			const double* P_3,const double* P_4,const double* P_5,const double* P_6, double dim_caracteristic, double precision);
+// 	vector<int> touv_maill_vois(int i_S,const monMaillageS& MaS);
+// 	void verif_point_dans_vect(const double* P, vector<double>& V,double dim_caracteristic, double precision );
+// 	void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
+// 														MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S);
+
+
+
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*   calcul la surface d'un triangle                  */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+	inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
+	{
+    double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]);
+    double Surface = 1/2.0*abs(A);
+    return Surface;
+	}
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*     fonction qui calcul le déterminant            */
+	/*      de deux vecteur(cf doc CGAL).                */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+	//fonction qui calcul le déterminant des vecteurs: P3P1 et P3P2
+	//(cf doc CGAL).
+
+	inline double mon_determinant(const double* P_1,
+																const double*  P_2,
+																const double* P_3)
+	{
+    double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
+    return mon_det;
+	}
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+	//calcul la norme du vecteur P1P2
+
+	inline double norme_vecteur(const double* P_1,const double* P_2)
+	{
+		double X=P_1[0]-P_2[0];
+		double Y=P_1[1]-P_2[1];
+		double norme=sqrt(X*X+Y*Y);
+		return norme;
+	}
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*         calcul le cos et le sin de l'angle P1P2,P1P3     */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+	inline vector<double> calcul_cos_et_sin(const double* P_1,
+																	 const double* P_2,
+																	 const double* P_3)
+	{
+	
+		vector<double> Vect;
+		double P1_P2=norme_vecteur(P_1,P_2);
+		double P2_P3=norme_vecteur(P_2,P_3);
+		double P3_P1=norme_vecteur(P_3,P_1);
+	
+		double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
+		double D=2.0*P1_P2*P3_P1;
+		double COS=N/D;
+		Vect.push_back(COS);
+		double V=mon_determinant(P_2,P_3,P_1);
+		double D_1=P1_P2*P3_P1;
+		double SIN=V/D_1;
+		Vect.push_back(SIN);
+	
+		return Vect;
+	
+	}
+
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*     calcul les coordonnées du barycentre d'un polygone   */ 
+	/*     le vecteur en entrée est constitué des coordonnées   */
+	/*     des sommets du polygone                              */                             
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+	inline vector<double> bary_poly(const vector<double>& V)
+	{
+    vector<double> Bary;
+    int taille=V.size();
+    double x=0;
+    double y=0;
+
+    for (int i=0;i<taille/2;i++)
+			{
+				x=x+V[2*i];
+				y=y+V[2*i+1];
+			}
+    double A=2*x/taille;
+    double B=2*y/taille;
+    Bary.push_back(A);//taille vecteur=2*nb de points.
+    Bary.push_back(B);
+
+
+    return Bary;
+	}
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*         calcul la surface d'un polygone.                 */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+	inline double  Surf_Poly(const vector<double>& Poly)
+	{ 
+
+    double Surface=0;
+    for (int i=0; i<(Poly.size())/2-2; i++)
+			{
+				double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] ); 
+				Surface=Surface + Surf ;
+			}
+    return Surface ;
+	}
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*   fonction qui teste si un point est dans une maille     */
+	/*   point: P_0                                             */
+	/*   P_1, P_2, P_3 sommet des mailles                       */   
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+	inline bool point_dans_triangle(const double* P_0,const double* P_1,
+													 const double* P_2,const double* P_3)
+	{
+
+    bool A=false;
+    double det_1=mon_determinant(P_1,P_3,P_0);
+    double det_2=mon_determinant(P_3,P_2,P_0);
+    double det_3=mon_determinant(P_2,P_1,P_0);
+    if( (det_1>=0 && det_2>=0 && det_3>=0) || (det_1<=0 && det_2<=0 && det_3<=0) )
+			{
+				A=true;
+			}
+
+    return A;
+	}
+
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*fonction pour vérifier qu'un point n'a pas déja été considérer dans   */ 
+	/*      le vecteur et le rajouter au vecteur sinon.                     */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+	inline void verif_point_dans_vect(const double* P, vector<double>& V,double dim_caracteristic, double precision )
+	{
+		//    int taille=V.size();
+		//   for(int i=0;i<taille/2;i++) 
+		//     {
+		//       double dx=P[0]-V[2*i];
+		//       dx=dx>0.0?dx:-dx;
+		//       if (dx>_Precision)
+		// 	continue;
+		//       double dy=P[1]-V[2*i+1];
+		//       dy=dy>0.0?dy:-dy;
+		//       if (dy<_Precision)
+	
+		// 	return;
+      
+		//     }
+		//   V.push_back(P[0]);
+		//   V.push_back(P[1]);
+  
+		int taille=V.size();
+		bool isPresent=false;
+		for(int i=0;i<taille/2;i++) 
+			{
+				if (sqrt(((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))/dim_caracteristic<precision)
+					isPresent=true;
+      
+			}
+		if(!isPresent)
+			{
+      
+				V.push_back(P[0]);
+				V.push_back(P[1]);    
+			}
+	}
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/* fonction qui rajoute les sommet du triangle P dans le vecteur V        */ 
+	/* si ceux-ci sont compris dans le triangle S et ne sont pas déjà dans    */
+	/* V.                                                                     */
+	/*sommets de P: P_1, P_2, P_3                                             */
+	/*sommets de S: P_4, P_5, P_6                                             */  
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ 
+
+	inline void rajou_sommet_triangl(const double* P_1,const double* P_2,const double* P_3,
+														const double* P_4,const double* P_5,const double* P_6,vector<double>& V, double dim_caracteristic, double precision)
+	{
+
+		bool A_1=ParaMEDMEM_utils::point_dans_triangle(P_1,P_4,P_5,P_6);
+    if(A_1)
+			{verif_point_dans_vect(P_1,V,dim_caracteristic, precision);
+				//   if (V.size()==1)
+				//       {
+				// 	// all nodes are necessarily in the triangle
+				// 	verif_point_dans_vect(P_2,V);
+				// 	verif_point_dans_vect(P_3,V);
+				// 	return;
+				//       }
+			}
+    bool A_2=ParaMEDMEM_utils::point_dans_triangle(P_2,P_4,P_5,P_6);
+    if(A_2)
+			{verif_point_dans_vect(P_2,V,dim_caracteristic,precision);}
+
+    bool A_3=ParaMEDMEM_utils::point_dans_triangle(P_3,P_4,P_5,P_6);
+    if(A_3)
+			{verif_point_dans_vect(P_3,V,dim_caracteristic, precision);}
+          
+	}
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/
+	/*  calcul de l'intersection de deux segments: segments P1P2 avec P3P4      */
+	/*  . Si l'intersection est non nulle et si celle-ci n'est                  */
+	/*  n'est pas déjà contenue dans Vect on la rajoute à Vect                  */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/ 
+
+	inline void inters_de_segment(const double * P_1,const double * P_2,
+																const double * P_3,const double * P_4,vector<double>&  Vect, 
+																double dim_caracteristic, double precision)
+	{
+
+
+    // calcul du déterminant de P1_1P_2 et P_3P_4.
+    double det=(P_2[0]-P_1[0])*(P_4[1]-P_3[1])-(P_4[0]-P_3[0])*(P_2[1]-P_1[1]);
+
+
+    if(abs(det)/dim_caracteristic>precision)
+			{
+
+				double k_1=1/((P_2[0]-P_1[0])*(P_3[1]-P_4[1])+(P_4[0]-P_3[0])*(P_2[1]-P_1[1]))
+					*((P_3[1]-P_4[1])*(P_3[0]-P_1[0])+(P_4[0]-P_3[0])*(P_3[1]-P_1[1]));
+	
+				if( k_1>=0 &&  k_1<=1)
+					{
+						double k_2=1/((P_4[0]-P_3[0])*(P_1[1]-P_2[1])+(P_2[0]-P_1[0])*(P_4[1]-P_3[1]))
+							*((P_1[1]-P_2[1])*(P_1[0]-P_3[0])+(P_2[0]-P_1[0])*(P_1[1]-P_3[1]));
+	    
+
+						if( k_2>=0 &&  k_2<=1)
+							{
+								double P_0[2];
+								P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]);
+								P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]);
+								verif_point_dans_vect(P_0,Vect, dim_caracteristic, precision);
+
+							}
+					}
+			}
+	}
+
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
+	/*      calcul l'intersection de deux triangles            */
+	/* P_1, P_2, P_3: sommets du premier triangle              */
+	/* P_4, P_5, P_6: sommets du deuxième triangle             */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
+
+	inline vector<double> intersec_de_triangle(const double* P_1,const double* P_2,
+																			const double* P_3,const double* P_4,const double* P_5,const double* P_6, double dim_caracteristic, double precision)
+	{
+
+		//		cout << "  T1 : " << *P_1 << " , " << *(P_1+1) << " , " << *P_2 << " , " << *(P_2+1) << " , " << *P_3 << " , " << *(P_3+1)<< endl;
+		// cout << "  T2 : " << *P_4 << " , " << *(P_4+1) << " , " << *P_5 << " , " << *(P_5+1) << " , " << *P_6 << " , " << *(P_6+1)<< endl;
+    vector<double> Vect;
+
+    inters_de_segment(P_1,P_2,P_4,P_5,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_1,P_2,P_5,P_6,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_1,P_2,P_6,P_4,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_2,P_3,P_4,P_5,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_2,P_3,P_5,P_6,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_2,P_3,P_6,P_4,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_3,P_1,P_4,P_5,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_3,P_1,P_5,P_6,Vect, dim_caracteristic, precision);
+    inters_de_segment(P_3,P_1,P_6,P_4,Vect, dim_caracteristic, precision);
+    rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect, dim_caracteristic, precision);
+    rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect, dim_caracteristic, precision);
+    //debug cout << "  Inter : ";
+    //debug for (int i=0; i<Vect.size(); ++i)
+		//debug cout << Vect[i] << "  ";
+    //debug cout << endl << endl;
+
+    return Vect;
+	}
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */ 
+	/*fonction pour trouver les mailles voisines d'une maille triangle.*/
+	/* (mailles voisines par les faces)                                */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ __ _ _ _ */ 
+
+	inline vector<int> touv_maill_vois(int i_S,const monMaillageS& MaS)
+	{
+
+    int N_1=MaS._Connect[3*(i_S-1)];
+    int N_2=MaS._Connect[3*(i_S-1)+1];
+    int N_3=MaS._Connect[3*(i_S-1)+2];
+    vector<int> tab;
+    int I_1=MaS._ReversNConnectIndex[N_1-1]-1;// Dans MED le n°des noeuds commencent à 1. 
+    int I_2=MaS._ReversNConnectIndex[N_2-1]-1;
+    int I_3=MaS._ReversNConnectIndex[N_3-1]-1;
+    int F_1=MaS._ReversNConnectIndex[N_1]-1;
+    int F_2=MaS._ReversNConnectIndex[N_2]-1;
+    int F_3=MaS._ReversNConnectIndex[N_3]-1;
+
+
+    vector<int> V_12;
+    V_12.reserve(2);
+    insert_iterator<vector<int> > ib_1(V_12,V_12.begin());
+    vector<int> V_23;
+    V_23.reserve(2);
+    insert_iterator<vector<int> > ib_2(V_23,V_23.begin());
+    vector<int> V_13;
+    V_13.reserve(2);
+    insert_iterator<vector<int> > ib_3(V_13,V_13.begin());
+
+
+    set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_2,
+										 MaS._ReversNConnect+F_2,ib_1);
+    set_intersection(MaS._ReversNConnect+I_2,MaS._ReversNConnect+F_2,MaS._ReversNConnect+I_3,
+										 MaS._ReversNConnect+F_3,ib_2);
+    set_intersection(MaS._ReversNConnect+I_1,MaS._ReversNConnect+F_1,MaS._ReversNConnect+I_3,
+										 MaS._ReversNConnect+F_3,ib_3);
+
+
+    for(int i=0;i<(V_12).size();i++)
+			{
+				if(V_12[i]!=i_S)
+					{
+						tab.push_back(V_12[i]);
+						break;
+					}
+			}
+
+    for(int i=0;i<V_23.size();i++)
+			{
+				if(V_23[i]!=i_S)
+					{
+						tab.push_back(V_23[i]);
+						break;
+					}
+			}
+
+    for(int i=0;i<V_13.size();i++)
+			{
+				if(V_13[i]!=i_S)
+					{
+						tab.push_back(V_13[i]);
+						break;
+					}
+			}
+    return tab;
+	}
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+	/* fonction pour vérifier qu'un n°de maille n'a pas déja été considérer  */
+	/*  dans le vecteur et le rajouter au vecteur sinon.                     */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+
+	inline void verif_maill_dans_vect(int Num, vector<int>& V)
+	{
+    int taille=V.size();
+    int A=0;
+    for(int i=0;i<taille;i++)
+			{
+				if(Num==V[i])
+					{
+						A=1;
+						break;
+					} 
+			}
+    if(A==0)
+			{V.push_back(Num); }
+	}
+
+
+
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */  
+	/* fonction pour reconstituer un polygone convexe à partir  */
+	/*              d'un nuage de point.                        */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */  
+
+	inline vector<double> reconstruct_polygon(const vector<double>& V)
+	{
+
+		int taille=V.size();
+
+		//VB : why 6 ?
+
+		if(taille<=6)
+			{return V;}
+		else
+			{
+				double COS[taille/2];
+				double SIN[taille/2];
+				double angle[taille/2];
+				vector<double> Bary=bary_poly(V);
+				COS[0]=1.0;
+				SIN[0]=0.0;
+				angle[0]=0.0;
+				for(int i=0; i<taille/2-1;i++)
+					{
+						vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
+						COS[i+1]=Trigo[0];
+						SIN[i+1]=Trigo[1];
+						if(SIN[i+1]>=0)
+							{angle[i+1]=acos(COS[i+1]);}
+						else
+							{angle[i+1]=-acos(COS[i+1]);}
+					}
+			
+				//ensuite on ordonne les angles.
+				vector<double> Pt_ordonne;
+				Pt_ordonne.reserve(taille);
+				multimap<double,int> Ordre;
+				for(int i=0;i<taille/2;i++)	
+					{Ordre.insert(make_pair(angle[i],i));}
+				multimap <double,int>::iterator mi;
+				for(mi=Ordre.begin();mi!=Ordre.end();mi++)
+					{
+						int j=(*mi).second;
+						Pt_ordonne.push_back(V[2*j]);
+						Pt_ordonne.push_back(V[2*j+1]);
+					}
+				return Pt_ordonne;
+			}
+	}
+
+
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+	/*  fonction qui écrit les résultats d'annelise dans un fichier:      */
+	/*  pour chaque maille du maillage 1 on stoque la fraction volumique  */
+	/*          considéré lors de l'algorithme                            */
+	/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+
+	inline void WriteFractSurfInFile(MEDMEM::FIELD<double>& FractSurf,MEDMEM::MESH* myMesh_P,
+														MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S)
+	{
+    //On sauvegarde le champ crée sur la maillage P.
+    string FileName="FracSurf"+NameMesh_S+"_to_"+NameMesh_P;
+    int id1=(*myMesh_P).addDriver(MED_DRIVER,FileName,"MP");
+    int id2=(*myMesh_S).addDriver(MED_DRIVER,FileName,"MS");
+    int id3=FractSurf.addDriver(MED_DRIVER,FileName,"% MP");
+    (*myMesh_P).write(id1);
+    (*myMesh_S).write(id2);
+    FractSurf.write(id3);
+	}
+
+	inline void getElemBB(double* bb, const monMaillageP& MaP, int iP)
+	{
+		bb[0]=HUGE;bb[1]=-HUGE;
+		bb[2]=HUGE;bb[3]=-HUGE;
+		bb[4]=HUGE;bb[5]=-HUGE;
+		
+		for (int i=0; i<3; i++)
+			{
+				double x = MaP._Coord[(MaP._Connect[3*iP+i]-1)*3];
+				double y = MaP._Coord[(MaP._Connect[3*iP+i]-1)*3+1];
+				double z = MaP._Coord[(MaP._Connect[3*iP+i]-1)*3+2];
+				bb[0]=(x<bb[0])?x:bb[0];
+				bb[1]=(x>bb[1])?x:bb[1];
+				bb[2]=(y<bb[2])?y:bb[2];
+				bb[3]=(y>bb[3])?y:bb[3];
+				bb[4]=(z<bb[4])?z:bb[4];
+				bb[5]=(z>bb[5])?z:bb[5];
+			}
+		
+	}
+
+};
+#endif
diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx
index 508760181..fa99bc41a 100644
--- a/src/ParaMEDMEM/InterpolationMatrix.cxx
+++ b/src/ParaMEDMEM/InterpolationMatrix.cxx
@@ -3,6 +3,7 @@
 #include "MxN_Mapping.hxx"
 #include "InterpolationMatrix.hxx"
 #include "INTERPOLATION_2D.hxx"
+#include "INTERPOLATION_3D_surf.hxx"
 
 /*! \class InterpolationMatrix
 This class enables the storage of an interpolation matrix Wij mapping 
@@ -51,10 +52,19 @@ adds the contribution of a distant subdomain to the interpolation matrix
  */
 void InterpolationMatrix::addContribution(MEDMEM::MESH& distant_support, int iproc_distant, int* distant_elems)
 {
-  INTERPOLATION_2D interpolator;
-  int nbelems=  _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
-  //_row_offsets.resize(nbelems);
-  vector<map<int,double> > surfaces = interpolator.interpol_maillages(distant_support,_source_support);
+
+	// computing the intersections between distant and local elements
+  INTERPOLATION* interpolator;
+	if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==3)
+		interpolator=new INTERPOLATION_3D_surf();
+	else if (distant_support.getMeshDimension()==2 && distant_support.getSpaceDimension()==2)
+		interpolator=new INTERPOLATION_2D();
+
+  int nbelems=  _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);  
+  vector<map<int,double> > surfaces = interpolator->interpol_maillages(distant_support,_source_support);
+	delete interpolator;
+
+
   if (surfaces.size() != nbelems)
     throw MEDEXCEPTION("uncoherent number of rows in interpolation matrix");
 
@@ -69,23 +79,31 @@ void InterpolationMatrix::addContribution(MEDMEM::MESH& distant_support, int ipr
          iter!=surfaces[ielem].end();
          iter++)
          {
+	   //surface of the source triangle
 	   double surf = triangle_surf->getValueIJ(iter->first,1);
+	   //oriented surfaces are not considered
 	   surf = (surf>0)?surf:-surf;
+	   //locating the (iproc, itriangle) pair in the list of columns
            vector<pair<int,int> >::iterator iter2 = find (_col_offsets.begin(), _col_offsets.end(),make_pair(iproc_distant,iter->first));
            int col_id;
+	   
+	   //(iproc, itriangle) is not registered in the list
+	   //of distant elements
            if (iter2==_col_offsets.end())
            {
              _col_offsets.push_back(make_pair(iproc_distant,iter->first));
              col_id =_col_offsets.size();
-	     _mapping.addElementFromSource(col_id,iproc_distant,distant_elems[iter->first-1]);
+	     _mapping.addElementFromSource(iproc_distant,distant_elems[iter->first-1]);
            } 
-           else
+           else 
            {
             col_id=iter2-_col_offsets.begin()+1;
            }
+	   //the column indirection number is registered
+	   //with the interpolation coefficient
+	   //actual column number is obtained with _col_offsets[col_id]
            _col_numbers.push_back(col_id);
            _coeffs.push_back(iter->second/surf);
-         
          }
   }
   delete triangle_surf;
@@ -154,7 +172,6 @@ void InterpolationMatrix::multiply(MEDMEM::FIELD<double>& field) const
     {
       target_value[i]=0.0;
     }
-  // int nbrows = _source_indices.size();
   int nbrows=  _source_support.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
   for (int irow=0; irow<nbrows; irow++)
   {
diff --git a/src/ParaMEDMEM/IntersectionDEC.cxx b/src/ParaMEDMEM/IntersectionDEC.cxx
index 346aede57..78a15ab71 100644
--- a/src/ParaMEDMEM/IntersectionDEC.cxx
+++ b/src/ParaMEDMEM/IntersectionDEC.cxx
@@ -49,8 +49,8 @@ void IntersectionDEC::synchronize()
       int* distant_ids=0;
       for (int i=0; i<_target_group->size(); i++)
       {
-	int idistant_proc = (i+_source_group->myRank())%_target_group->size();
-	//	idistant_proc=i;
+	//	int idistant_proc = (i+_source_group->myRank())%_target_group->size();
+	int	idistant_proc=i;
 
         //gathers pieces of the target meshes that can intersect the local mesh
         locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
@@ -81,8 +81,8 @@ void IntersectionDEC::synchronize()
       int* distant_ids=0;
       for (int i=0; i<_source_group->size(); i++)
       {
-	int idistant_proc = (i+_target_group->myRank())%_source_group->size();
-       
+	//	int idistant_proc = (i+_target_group->myRank())%_source_group->size();
+	int  idistant_proc=i;
         //gathers pieces of the target meshes that can intersect the local mesh
         locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
 	cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<endl;
diff --git a/src/ParaMEDMEM/Makefile.in b/src/ParaMEDMEM/Makefile.in
index 2778fa0e8..1aaddbdb4 100644
--- a/src/ParaMEDMEM/Makefile.in
+++ b/src/ParaMEDMEM/Makefile.in
@@ -54,13 +54,18 @@ ParaFIELD.hxx\
 DEC.hxx\
 MxN_Mapping.hxx\
 INTERPOLATION_2D.hxx\
+INTERPOLATION_3D_surf.hxx\
+INTERPOLATION.hxx\
+INTERPOLATION_Utils.hxx\
 StructuredCoincidentDEC.hxx\
 InterpolationMatrix.hxx\
 IntersectionDEC.hxx\
 UnstructuredParaSUPPORT.hxx\
 ExplicitCoincidentDEC.hxx\
 ElementLocator.hxx\
-ExplicitMapping.hxx
+ExplicitMapping.hxx\
+TranslationRotationMatrix.hxx\
+BBTree.H
 
 # Libraries targets
 
@@ -77,6 +82,7 @@ ComponentTopology.cxx\
 ParaFIELD.cxx\
 DEC.cxx\
 INTERPOLATION_2D.cxx\
+INTERPOLATION_3D_surf.cxx\
 MxN_Mapping.cxx\
 InterpolationMatrix.cxx\
 StructuredCoincidentDEC.cxx\
@@ -95,7 +101,7 @@ BIN_SERVER_IDL =
 BIN_CLIENT_IDL = 
 
 TEST_PROGS = test_ProcessorGroup test_BlockTopology test_ParaStructuredSupport \
-test_ParaField test_DEC test_UnstructuredDEC test_ExplicitDEC test_IntersectionDEC
+test_ParaField test_DEC test_UnstructuredDEC test_ExplicitDEC test_IntersectionDEC test_SeqIntersectionDEC test_Seq3DsurfIntersection test_MEDMEMConstructor make_cylinder make_plane test_INTERPOL_2D
 
 LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome 
 LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
diff --git a/src/ParaMEDMEM/MxN_Mapping.cxx b/src/ParaMEDMEM/MxN_Mapping.cxx
index ea93fa0c4..8dc6fa749 100644
--- a/src/ParaMEDMEM/MxN_Mapping.cxx
+++ b/src/ParaMEDMEM/MxN_Mapping.cxx
@@ -19,10 +19,15 @@ MxN_Mapping::~MxN_Mapping()
   delete _union_group;
 }
 
-void MxN_Mapping::addElementFromSource(int local_element, int distant_proc, int distant_element)
+
+/*!
+Method registering a new element for correspondence with a distant element
+\param distant_proc proc rank of the distant processor (in terms of the union group)
+\param distant_element id of the element on the distant processor
+ */
+void MxN_Mapping::addElementFromSource(int distant_proc, int distant_element)
 {
   _sending_ids.push_back(make_pair(distant_proc,distant_element));
- // _local_ids.push_back(local_element);
   for (int i=distant_proc; i<_union_group->size(); i++)
     _send_proc_offsets[i+1]++;
 }
diff --git a/src/ParaMEDMEM/MxN_Mapping.hxx b/src/ParaMEDMEM/MxN_Mapping.hxx
index af58fae4f..e38360bf3 100644
--- a/src/ParaMEDMEM/MxN_Mapping.hxx
+++ b/src/ParaMEDMEM/MxN_Mapping.hxx
@@ -16,7 +16,7 @@ public:
 	MxN_Mapping();
   MxN_Mapping(const ProcessorGroup& local_group, const ProcessorGroup& distant_group);
 	virtual ~MxN_Mapping();
-  void addElementFromSource(int local_elem, int distant_proc, int distant_elem);
+  void addElementFromSource(int distant_proc, int distant_elem);
   void prepareSendRecv();
   void sendRecv(MEDMEM::FIELD<double>& field);
   void sendRecv(double* field, MEDMEM::FIELD<double>& field) const ;
diff --git a/src/ParaMEDMEM/ParaFIELD.cxx b/src/ParaMEDMEM/ParaFIELD.cxx
index 07aada33d..50b4988bc 100644
--- a/src/ParaMEDMEM/ParaFIELD.cxx
+++ b/src/ParaMEDMEM/ParaFIELD.cxx
@@ -52,15 +52,24 @@ ParaFIELD::ParaFIELD(const ParaSUPPORT* para_support, const ComponentTopology& c
 	_field->setNumberOfValues(para_support->getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS));
 	string* compnames=new string[nb_components];
 	string* compdesc=new string[nb_components];
+		string* compunit=new string[nb_components];
 	for (int i=0; i<nb_components; i++)
 	{
 		ostringstream stream(compnames[i]);
 		ostringstream stream2(compdesc[i]);
+		ostringstream stream3(compunit[i]);
 		stream<<"component "<<i;
 		stream2<<"component description "<<i;
+				stream3<<"component unit "<<i;
+
+		compnames[i]=stream.str();
+		compdesc[i]=stream2.str();
+				compunit[i]=stream3.str();
+
 	}
 	_field->setComponentsNames(compnames);
 	_field->setComponentsDescriptions(compdesc);
+	_field->setMEDComponentsUnits(compunit);
 	_field->setIterationNumber(0);
 	_field->setOrderNumber(0);
 	_field->setTime(0.0);
diff --git a/src/ParaMEDMEM/ParaMESH.cxx b/src/ParaMEDMEM/ParaMESH.cxx
index 6b9a5f064..153cb79e1 100644
--- a/src/ParaMEDMEM/ParaMESH.cxx
+++ b/src/ParaMEDMEM/ParaMESH.cxx
@@ -200,10 +200,13 @@ throw (MEDMEM::MEDEXCEPTION){
 
 ParaMESH::ParaMESH(MEDMEM::MESH& subdomain_mesh, const ProcessorGroup& proc_group, const string& name):
 _mesh(&subdomain_mesh),
-_name (name),
 _my_domain_id(proc_group.myRank()),
 _block_topology (new BlockTopology(proc_group, subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS)))
 {
+	ostringstream stream;
+	stream<<name<<"_"<<_my_domain_id+1;
+	_name=stream.str();
+
   _cellglobal = new int[subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS)];
   int offset = _block_topology->localToGlobal(make_pair(_my_domain_id,0));
   for (int i=0; i<subdomain_mesh.getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS); i++)
diff --git a/src/ParaMEDMEM/TranslationRotationMatrix.hxx b/src/ParaMEDMEM/TranslationRotationMatrix.hxx
new file mode 100644
index 000000000..3cabb74ec
--- /dev/null
+++ b/src/ParaMEDMEM/TranslationRotationMatrix.hxx
@@ -0,0 +1,115 @@
+#ifndef _TRANSLATIONROTATIONMATRIX_HXX_
+#define _TRANSLATIONROTATIONMATRIX_HXX_
+
+namespace ParaMEDMEM
+{
+
+  const double EPS=1e-12;
+
+class TranslationRotationMatrix
+{
+
+public:
+
+  TranslationRotationMatrix():_rotation_coeffs(9,0.0), _translation_coeffs(3,0.0)
+  { 
+    _rotation_coeffs[0]=1.0;
+    _rotation_coeffs[4]=1.0;
+    _rotation_coeffs[8]=1.0;
+  }
+
+ void multiply(const TranslationRotationMatrix& A)
+  {
+    TranslationRotationMatrix result;
+    //setting output matrix to zero
+    for (int i=0; i<3; i++)
+      result._rotation_coeffs[i*4]=0.0;
+    //multiplying
+    for (int i=0; i<3;i++)
+      for (int j=0; j<3; j++)
+	for (int k=0; k<3; k++)
+	  result._rotation_coeffs[j+i*3]+=A._rotation_coeffs[3*i+k]*_rotation_coeffs[j+k*3];
+    
+    for (int i=0;i<9; i++)
+      _rotation_coeffs[i]=result._rotation_coeffs[i];
+
+	// 	cout << "Matrix "<<endl;
+// 		for (int i=0; i<3; i++)
+// 			{
+// 				cout<<"| ";
+// 				for (int j=0; j<3;j++)
+// 					cout <<_rotation_coeffs[i*3+j]<<" ";
+// 				cout<<" "<<endl;
+// 			}
+  }
+  
+  void rotate_vector(double* P)
+  {
+    double temp[3]={0.0, 0.0, 0.0};
+    
+    for (int i=0; i<3;i++)
+      for (int j=0; j<3; j++)
+	  temp[i] +=_rotation_coeffs[3*i+j]*P[j];
+	
+    P[0]=temp[0];P[1]=temp[1];P[2]=temp[2];
+  }
+ 
+  void transform_vector(double*P)
+  {
+    P[0]+=_translation_coeffs[0];
+    P[1]+=_translation_coeffs[1];
+    P[2]+=_translation_coeffs[2];
+    rotate_vector(P);
+  }
+  void translate(const double* P)
+  {
+    _translation_coeffs[0]=P[0];
+    _translation_coeffs[1]=P[1];
+    _translation_coeffs[2]=P[2];
+  }
+  
+  void  rotate_x (double* P)
+  {
+    _rotation_coeffs[0]=1.0;
+    double r_sqr = P[1]*P[1]+P[2]*P[2];
+    if (r_sqr < EPS)
+      {_rotation_coeffs[4]=1.0; _rotation_coeffs[8]=1.0; return;}
+    double r = sqrt(r_sqr);
+    double cos =P[1]/r;
+    double sin =P[2]/r;
+
+		_rotation_coeffs[4]=cos;
+		_rotation_coeffs[5]=sin;
+		_rotation_coeffs[7]=-sin;
+		_rotation_coeffs[8]=cos;
+
+
+    rotate_vector(P);
+  }
+   void  rotate_z (double* P)
+  {
+    _rotation_coeffs[8]=1.0;
+    double r_sqr = P[0]*P[0]+P[1]*P[1];
+    if (r_sqr < EPS)
+      {_rotation_coeffs[4]=1.0; _rotation_coeffs[0]=1.0; return;}
+    double r = sqrt(r_sqr);
+    double cos =P[0]/r;
+    double sin =P[1]/r;
+    
+    _rotation_coeffs[0]=cos;
+    _rotation_coeffs[1]=sin; 
+    _rotation_coeffs[3]=-sin;
+    _rotation_coeffs[4]=cos;
+    
+    rotate_vector(P);
+  }
+			
+	
+private:
+  vector<double> _rotation_coeffs;
+  vector<double> _translation_coeffs;
+  
+};
+
+};
+#endif