--- /dev/null
+#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);
+}
+};
--- /dev/null
+#ifndef _INTERPOLATION_HXX_
+#define _INTERPOLATION_HXX_
+
+
+#include <vector>
+#include <map>
+
+#include "MEDMEM_Mesh.hxx"
+
+
+namespace MEDMEM
+{
+
+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
--- /dev/null
+#include "InterpolationUtils.hxx"
+#include "Interpolation.hxx"
+#include "Interpolation2D.hxx"
+
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Field.hxx"
+
+#include "MEDMEM_InterpolationHighLevelObjects.hxx"
+#include "MEDMEM_Interpolation.hxx"
+
+using namespace std;
+using namespace MED_EN;
+using namespace MEDMEM;
+
+
+namespace MEDMEM
+{
+
+ Interpolation2D::Interpolation2D()
+ {
+ _Precision=1.0E-12;
+ _DimCaracteristic=1;
+ _SecondMethodOfLocalisation=1;
+ _PrintLevel=0;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* Options : */
+ /* Precision : for geometric computation */
+ /* SecondMethodOfLocalisation : 0/1 */
+ /* PrintLevel : between 0 and 3 */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ void Interpolation2D::setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel)
+ {
+ _Precision=Precision;
+ _SecondMethodOfLocalisation=SecondMethodOfLocalisation;
+ _PrintLevel=PrintLevel;
+ }
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /*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), _DimCaracteristic, _Precision);
+
+ 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) , _DimCaracteristic, _Precision);
+
+ 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), _DimCaracteristic, _Precision);
+
+ //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.
+ vector<map<int,double> > Interpolation2D::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;
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ 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
+ 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)
+ {
+ 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;
+ }
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //on cherche pour chaque maille i_P du maillage projetté
+ // une maille i_S du maillage S qui chevauche la maille i_P.
+
+ 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);
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //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);
+
+ /*à remplacer par:
+ WriteFractSurfInFile(vector<<vector<pair<int,double> > >& FractVol,MEDMEM::MESH* myMesh1,
+ MEDMEM::MESH* myMesh2,string FileName,string NameMeshP,string NameMeshS)
+ */
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ //boucle sur les mailles de P.
+ //Coeur de l'algorithme
+
+ 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);
+
+ }
+
+ if (_PrintLevel >= 2)
+ {
+ cout<<endl<<"Impression des surfaces d'intersection:"<<endl<<endl;
+ cout<<"(maille source, maille cible):surface d'intersection"<<endl;
+ 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;
+ }
+ }
+ }
+
+ // /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ // //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;
+
+ }
+
+
+};
--- /dev/null
+#ifndef _INTERPOLATION2D_HXX_
+#define _INTERPOLATION2D_HXX_
+
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Field.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 Interpolation;
+
+class Interpolation2D : public Interpolation
+{
+
+ private:
+ double _Precision;
+ double _DimCaracteristic;
+ int _SecondMethodOfLocalisation;
+ int _PrintLevel;
+
+ // Méthodes publiques
+ public:
+
+ Interpolation2D();
+
+ // precision geometrique, choix de la methode comlementaire de localisation, niveau d'impressions
+ void setOptions(double Precision, int SecondMethodOfLocalisation, int PrintLevel);
+
+
+ //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> > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+
+ private:
+
+
+};
+
+};
+
+#endif
--- /dev/null
+#include "MEDMEM_Field.hxx"
+#include "MEDMEM_Mesh.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "InterpolationUtils.hxx"
+#include "Interpolation.hxx"
+#include "Interpolation3DSurf.hxx"
+#include "BBTree.H"
+
+using namespace std;
+using namespace MED_EN;
+using namespace MEDMEM;
+using namespace INTERP_UTILS;
+
+namespace MEDMEM
+{
+
+ Interpolation3DSurf::Interpolation3DSurf()
+ {
+ _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 Interpolation3DSurf::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 Interpolation3DSurf::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 Interpolation3DSurf::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 =INTERP_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=INTERP_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> > Interpolation3DSurf::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;
+
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ 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 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 << "Interpolation3DSurf::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];
+ INTERP_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 Interpolation3DSurf::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 Interpolation3DSurf::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;
+
+ }
+ }
+
+};
--- /dev/null
+#ifndef _INTERPOLATION_3D_SURF_HXX_
+#define _INTERPOLATION_3D_SURF_HXX_
+
+
+
+#include "MEDMEM_Mesh.hxx"
+#include "InterpolationUtils.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "Interpolation.hxx"
+
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <string>
+#include <set>
+#include <algorithm>
+#include <vector>
+
+namespace MEDMEM
+{
+
+ class Interpolation3DSurf : public Interpolation
+ {
+
+
+ private:
+ double _Precision;
+ double _DimCaracteristic;
+ int _SecondMethodOfLocalisation;
+ int _PrintLevel;
+ double _Surf3DAdjustmentEps;
+ // Méthodes publiques
+ public:
+
+ Interpolation3DSurf();
+
+
+ // 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 INTERP_UTILS::monMaillageP& MaP,
+ const INTERP_UTILS::monMaillageS& MaS,
+ vector<map<int,double> >& Surface_d_intersect);
+
+
+ void createBoundingBoxes(const struct INTERP_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
--- /dev/null
+#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 INTERP_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=INTERP_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=INTERP_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=INTERP_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
--- /dev/null
+# MED MEDMEM : MED files in memory
+#
+# Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+# CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+#
+#
+# File : Makefile.in
+# Author : Vincent BERGEAUD (CEA/DEN/DANS/DM2S/SFME/LGLS)
+# Module : MED
+
+top_srcdir=@top_srcdir@
+top_builddir=../..
+srcdir=@srcdir@
+VPATH=.:$(srcdir):$(srcdir)/tests
+
+MACHINE=PCLINUX
+
+@COMMENCE@
+
+
+EXPORT_PYSCRIPTS = \
+
+
+EXPORT_HEADERS = \
+Interpolation2D.hxx\
+Interpolation3DSurf.hxx\
+Interpolation.hxx\
+InterpolationUtils.hxx\
+TranslationRotationMatrix.hxx\
+BBTree.H
+
+# Libraries targets
+
+LIB=libinterpkernel.la
+
+LIB_SRC = \
+Interpolation2D.cxx\
+Interpolation3DSurf.cxx
+
+
+
+# Executables targets
+BIN =
+BIN_SRC =
+BIN_SERVER_IDL =
+BIN_CLIENT_IDL =
+
+TEST_PROGS = make_plane test_INTERPOL_2D
+
+LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
+LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome
+
+CXXFLAGS+=@CXXTMPDPTHFLAGS@ $(MED2_INCLUDES)
+CPPFLAGS+=$(BOOST_CPPFLAGS)
+#LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS)
+# change motivated by the bug KERNEL4778.
+LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil
+
+#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
+# change motivated by the bug KERNEL4778.
+LDFLAGSFORBIN+= -lm $(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 -lmedmem $(BOOST_LIBS) -lutil
+
+ifeq ($(MED_WITH_KERNEL),yes)
+ CPPFLAGS+= ${KERNEL_CXXFLAGS}
+ CXXFLAGS+= ${KERNEL_CXXFLAGS}
+ LDFLAGS+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace
+ LDFLAGSFORBIN+= ${KERNEL_LDFLAGS} -lSALOMELocalTrace -lSALOMEBasics
+endif
+
+LIBSFORBIN=$(BOOSTLIBS) $(MPI_LIBS)
+
+LIBS=
+
+# build create_mesh :
+bin:
+
+@CONCLUDE@
--- /dev/null
+#ifndef _TRANSLATIONROTATIONMATRIX_HXX_
+#define _TRANSLATIONROTATIONMATRIX_HXX_
+
+namespace MEDMEM
+{
+
+ 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
--- /dev/null
+#include "Interpolation2D.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+int main()
+{
+ MEDMEM::MESH source(MED_DRIVER,"/home/vb144235/resources/square128000.med","Mesh_1");
+ MEDMEM::MESH target(MED_DRIVER,"/home/vb144235/resources/square30000.med","Mesh_1");
+
+ MEDMEM::Interpolation2D interp;
+ // interp.setOptions(1e-6,1,2);
+ interp.interpol_maillages(source,target);
+}
