From bd869a08d0764d72e0f895ea7157cd85c97c048f Mon Sep 17 00:00:00 2001 From: vbd Date: Thu, 28 Jun 2007 13:00:22 +0000 Subject: [PATCH] creation of INTERP_KERNEL directory to contain the interpolation routines previously in ParaMEDMEM --- src/INTERP_KERNEL/BBTree.H | 134 +++++ src/INTERP_KERNEL/Interpolation.hxx | 26 + src/INTERP_KERNEL/Interpolation2D.cxx | 442 ++++++++++++++ src/INTERP_KERNEL/Interpolation2D.hxx | 76 +++ src/INTERP_KERNEL/Interpolation3DSurf.cxx | 396 +++++++++++++ src/INTERP_KERNEL/Interpolation3DSurf.hxx | 71 +++ src/INTERP_KERNEL/InterpolationUtils.hxx | 553 ++++++++++++++++++ src/INTERP_KERNEL/Makefile.in | 94 +++ .../TranslationRotationMatrix.hxx | 117 ++++ src/INTERP_KERNEL/test_INTERPOL_2D.cxx | 12 + 10 files changed, 1921 insertions(+) create mode 100644 src/INTERP_KERNEL/BBTree.H create mode 100644 src/INTERP_KERNEL/Interpolation.hxx create mode 100755 src/INTERP_KERNEL/Interpolation2D.cxx create mode 100755 src/INTERP_KERNEL/Interpolation2D.hxx create mode 100644 src/INTERP_KERNEL/Interpolation3DSurf.cxx create mode 100644 src/INTERP_KERNEL/Interpolation3DSurf.hxx create mode 100644 src/INTERP_KERNEL/InterpolationUtils.hxx create mode 100644 src/INTERP_KERNEL/Makefile.in create mode 100644 src/INTERP_KERNEL/TranslationRotationMatrix.hxx create mode 100644 src/INTERP_KERNEL/test_INTERPOL_2D.cxx diff --git a/src/INTERP_KERNEL/BBTree.H b/src/INTERP_KERNEL/BBTree.H new file mode 100644 index 000000000..3bc3e36bc --- /dev/null +++ b/src/INTERP_KERNEL/BBTree.H @@ -0,0 +1,134 @@ +#include +#include +#include + +const int MIN_NB_ELEMS =10; +const int MAX_LEVEL=10; +using namespace std; + +template class BBTree +{ + +private: + BBTree* _left; + BBTree* _right; + BBTree* _center; + int _level; + double _median; + double* _bb; + vector _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 nodes (nbelems*2); + for (int i=0; i::iterator median = nodes.begin()+nbelems; + nth_element< vector::iterator>(nodes.begin(), nodes.begin()+nbelems, nodes.end()); + _median=*median; + // std:: cout << *median < new_elems_left; + vector new_elems_right; + vector new_elems_center; + for (int i=0; i*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& 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; idimbb[idim*2+1] || bb_ptr[idim*2+1]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/INTERP_KERNEL/Interpolation.hxx b/src/INTERP_KERNEL/Interpolation.hxx new file mode 100644 index 000000000..787a99982 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation.hxx @@ -0,0 +1,26 @@ +#ifndef _INTERPOLATION_HXX_ +#define _INTERPOLATION_HXX_ + + +#include +#include + +#include "MEDMEM_Mesh.hxx" + + +namespace MEDMEM +{ + +class Interpolation +{ +public: + Interpolation(){} + virtual ~Interpolation(){} + + //interpolation of two triangular meshes. + virtual std::vector > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)=0; + +}; + +}; +#endif diff --git a/src/INTERP_KERNEL/Interpolation2D.cxx b/src/INTERP_KERNEL/Interpolation2D.cxx new file mode 100755 index 000000000..9bf175bd8 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation2D.cxx @@ -0,0 +1,442 @@ +#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& localis ) + { + + //calcul des coordonnées des barycentre des mailles de P + + + //calcule des coordonnees du barycentre d'une cellule + std::auto_ptr Support ( new SUPPORT(&myMesh_P,"monSupport",MED_CELL) ); + std::auto_ptr > 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 monInterpol_S(myMesh_S); + for(int i_P=0;i_P 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 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 >& Surface_d_intersect,FIELD& 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 maill_deja_vu_S ; + + maill_deja_vu_S.push_back(i_S1); + + for (int M_S=0; M_S_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 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 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 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 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 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 > 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"< > BoxS=myMesh_S.getBoundingBox(); + vector > 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 :"< localis; + localis.reserve(MaP._NbMaille); + MEDMEM::MESH* ptr_S = const_cast(&myMesh_S); + MEDMEM::MESH* ptr_P = const_cast(&myMesh_P); + + cout << "Interpolation2D::local_iP_dansS"< MAP_init; + vector > 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 mySupport_P (new SUPPORT(const_cast(&myMesh_P),"Support on all CELLS",MED_CELL) ); + MEDMEM::FIELD 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< > >& 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_P0) + rempli_vect_d_intersect(i_P+1,i_S_init,MaP,MaS,Surface_d_intersect,myField_P); + + } + + if (_PrintLevel >= 2) + { + cout<::iterator surface; + for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++) + cout<<" ("<=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<::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 < +#include +#include +#include +#include +#include +#include +#include + +// namespace MEDMEM +// { +// class MESH; +// template 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& 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 >& Surface_d_intersect, + FIELD& myField_P); + + + //fonction principale pour interpoler deux maillages triangulaires. + std::vector > interpol_maillages(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2); + + private: + + +}; + +}; + +#endif diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.cxx b/src/INTERP_KERNEL/Interpolation3DSurf.cxx new file mode 100644 index 000000000..cafa15b7b --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3DSurf.cxx @@ -0,0 +1,396 @@ +#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 >& 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 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 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 > 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"< > BoxS=myMesh_S.getBoundingBox(); + vector > 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 :"< 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"< MAP_init; + vector > 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 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= 2) + { + cout<::iterator surface; + for( surface=Surface_d_intersect[i].begin();surface!=Surface_d_intersect[i].end();surface++) + { + cout<<" ("<& 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+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+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+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& bbox) + { + /* We build the segment tree for locating possible matching intersections*/ + + int size = bbox.size()/6; + for (int i=0; i +#include +#include +#include +#include +#include +#include + +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 > 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 >& Surface_d_intersect); + + + void createBoundingBoxes(const struct INTERP_UTILS::monMaillageS mesh, vector& bbox); + + void rotate3DTriangle(const double* PP1, + const double* PP2, + const double* PP3, + TranslationRotationMatrix& matrix); + + void adjustBoundingBoxes(vector& bbox); + + + }; + +}; + +#endif diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx new file mode 100644 index 000000000..65fedf9fd --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -0,0 +1,553 @@ +#ifndef _INTERPOLATION_UTILS_HXX_ +#define _INTERPOLATION_UTILS_HXX_ + + +#include "MEDMEM_Mesh.hxx" +#include "MEDMEM_Field.hxx" + + +#include +#include +#include +#include + + +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& 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 calcul_cos_et_sin(const double* P_1, + // const double* P_2, + // const double* P_3); + // inline vector bary_poly(const vector& 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& Poly); + // vector reconstruct_polygon(const vector& 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& 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& Vect, + // double dim_caracteristic, + // double precision); + // vector 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 touv_maill_vois(int i_S,const monMaillageS& MaS); + // void verif_point_dans_vect(const double* P, vector& V,double dim_caracteristic, double precision ); + // void WriteFractSurfInFile(MEDMEM::FIELD& 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 calcul_cos_et_sin(const double* P_1, + const double* P_2, + const double* P_3) + { + + vector 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 bary_poly(const vector& V) + { + vector Bary; + int taille=V.size(); + double x=0; + double y=0; + + for (int i=0;i& 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& V,double dim_caracteristic, double precision ) + { + // int taille=V.size(); + // for(int i=0;i0.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& 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& 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 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 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 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 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 V_12; + V_12.reserve(2); + insert_iterator > ib_1(V_12,V_12.begin()); + vector V_23; + V_23.reserve(2); + insert_iterator > ib_2(V_23,V_23.begin()); + vector V_13; + V_13.reserve(2); + insert_iterator > 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) + { + int taille=V.size(); + int A=0; + for(int i=0;i reconstruct_polygon(const vector& 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 Bary=bary_poly(V); + COS[0]=1.0; + SIN[0]=0.0; + angle[0]=0.0; + for(int i=0; i 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 Pt_ordonne; + Pt_ordonne.reserve(taille); + multimap Ordre; + for(int i=0;i::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& 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]=(xbb[1])?x:bb[1]; + bb[2]=(ybb[3])?y:bb[3]; + bb[4]=(zbb[5])?z:bb[5]; + } + + } + +}; +#endif diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in new file mode 100644 index 000000000..42f8b876b --- /dev/null +++ b/src/INTERP_KERNEL/Makefile.in @@ -0,0 +1,94 @@ +# 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@ diff --git a/src/INTERP_KERNEL/TranslationRotationMatrix.hxx b/src/INTERP_KERNEL/TranslationRotationMatrix.hxx new file mode 100644 index 000000000..a5712e406 --- /dev/null +++ b/src/INTERP_KERNEL/TranslationRotationMatrix.hxx @@ -0,0 +1,117 @@ +#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 "< _rotation_coeffs; + vector _translation_coeffs; + + }; + +}; +#endif diff --git a/src/INTERP_KERNEL/test_INTERPOL_2D.cxx b/src/INTERP_KERNEL/test_INTERPOL_2D.cxx new file mode 100644 index 000000000..a2f9e7b3e --- /dev/null +++ b/src/INTERP_KERNEL/test_INTERPOL_2D.cxx @@ -0,0 +1,12 @@ +#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); +} -- 2.39.2