src/INTERPOLATION/MEDMEM_Mapping.hxx

Go to the documentation of this file.
00001 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00002 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00003 // 
00004 // This library is free software; you can redistribute it and/or
00005 // modify it under the terms of the GNU Lesser General Public
00006 // License as published by the Free Software Foundation; either 
00007 // version 2.1 of the License.
00008 // 
00009 // This library is distributed in the hope that it will be useful 
00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 
00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
00012 // Lesser General Public License for more details.
00013 //
00014 // You should have received a copy of the GNU Lesser General Public  
00015 // License along with this library; if not, write to the Free Software 
00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00017 //
00018 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00019 //
00020 #ifndef MEDMEM_MAPPING_HXX
00021 #define MEDMEM_MAPPING_HXX
00022 
00023 #include "MEDMEM_MappingTools.hxx"
00024 #include "MEDMEM_dTree.hxx"
00025 
00026 #define NBR_MAX_MAILLES_EXAMINEES 100
00027 
00028 #ifndef  NBR_FACES_MAX
00029 #define NBR_FACES_MAX 6
00030 #endif
00031 
00032 #define _TEMPLATE_ template <class MAILLAGE, class NUAGEMAILLE, class NUAGENOEUD, class NOEUD, int DIMENSION>
00033 #define _MAPPING_ Mapping<MAILLAGE,NUAGEMAILLE,NUAGENOEUD,NOEUD,DIMENSION>
00034 
00035 
00041 
00042 /*********************************************************/
00043 /*                                                       */
00044 /*                   Classe  Mapping                     */
00045 /*                                                       */
00046 /*********************************************************/
00047 
00048 // ATTENTION LE NUAGE DE NOEUD EST SUPPOSE NON REDONDANT ET AUCUNE VERIFICATION N'EST FAITE !
00049 
00050 template <class MAILLAGE, class NUAGEMAILLE, class NUAGENOEUD, class NOEUD, int DIMENSION> class Mapping
00051 {
00052 protected :
00053      MAILLAGE * maillage_back;
00054      NUAGEMAILLE * mailles_back;
00055      NUAGENOEUD * noeuds_back;
00056      NUAGENOEUD * noeuds_front;
00057      Coordonnees_Barycentriques<NUAGEMAILLE,NUAGENOEUD,NOEUD,DIMENSION> * CB;
00058      dTree<NOEUD,NUAGENOEUD,DIMENSION> * my_dTree;
00059      vector<int> resultat_mapping;
00060      vector<int> point_le_plus_proche;
00061      
00062 public :
00063      Mapping():maillage_back(NULL),mailles_back(NULL),noeuds_back(NULL),noeuds_front(NULL),CB(NULL),my_dTree(NULL) {}
00064      Mapping(MAILLAGE * mb,NUAGENOEUD * nb,NUAGENOEUD * nf); // le dTree est crée à l'initialisation, par contre, le mapping lui meme doit etre invoqué
00065      ~Mapping() {if (CB) delete CB;if (my_dTree) delete my_dTree;}
00066         dTree<NOEUD,NUAGENOEUD,DIMENSION> * Donne_dTree() {return my_dTree;}
00067      enum { INTERIEUR = 1, EXTERIEUR_AU_MILIEU = -1, EXTERIEUR_AU_BORD = -2 };
00068      int Donne_Directions(int num_maille,const NOEUD &n,int etat_face[NBR_FACES_MAX]);
00069      // Méthode interne de localisation
00070      int Trouve_Maille_Contenant_Point_Mth_Co(const NOEUD &n,int num_maille,int num_maille_interdit,int max_loop,int &nbr_mailles_examinees,int flag_convexe);
00071      void Cree_Mapping(int flag_convexe=0);                                             // SUPPOSE NON CONVEXE PAR DEFAUT
00072      void Cree_Mapping(NUAGENOEUD * nf, int flag_convexe=0);                            // SUPPOSE NON CONVEXE PAR DEFAUT
00073      inline int operator[](int i) const {return resultat_mapping[i];}                   // Renvoie la valeur mappé, si le mapping a été fait, sinon, n'importe quoi
00074      inline vector<int> & Get_Mapping() {return resultat_mapping;}                        // Renvoie le vector contenant le mapping
00075      inline int Get_Noeud_Le_Plus_Proche(int i) const {return point_le_plus_proche[i];} // Invoque la méthode de d-Tree qui donne le noeud le plus proche
00076      inline int Exist_dTree() const {return (my_dTree);}                                // Teste si le dTree existe
00077      void affiche()
00078           {
00079           for (int i=0;i<resultat_mapping.size();i++) cout<<"Noeud "<<i<<" de noeud plus proche "<<point_le_plus_proche[i]<<" mappé dans maille "<<resultat_mapping[i]<<endl;       
00080           }
00081 };
00082 
00088 
00089 _TEMPLATE_ void _MAPPING_::Cree_Mapping(NUAGENOEUD * nf, int flag_convexe)
00090      {
00091      noeuds_front=nf;
00092      Cree_Mapping(flag_convexe);
00093      }
00094 
00095 _TEMPLATE_ void _MAPPING_::Cree_Mapping(int flag_convexe)
00096      {
00097      
00098      if (resultat_mapping.size()==0) 
00099           {
00100           if (!noeuds_front) 
00101                 {
00102                 cerr<<"Mapping : Pas de noeuds à mapper !"<<endl;
00103                 exit(-1);
00104                 }
00105            
00106           int i;
00107           int nbr_noeuds=noeuds_front->SIZE();
00108           int num_maille_depart;
00109           int nma=0;
00110           resultat_mapping     = vector<int>(nbr_noeuds,UNDEFINED);
00111           point_le_plus_proche = vector<int>(nbr_noeuds,UNDEFINED);
00112      
00113           // noeuds_back->affiche();
00114           
00115           for (i=0;i<nbr_noeuds;i++)
00116                {
00117                point_le_plus_proche[i]=my_dTree->trouve_plus_proche_point((*noeuds_front)[i]);
00118                num_maille_depart=maillage_back->DONNE_PREMIERE_MAILLE_CONTENANT(point_le_plus_proche[i]);
00119                resultat_mapping[i]=Trouve_Maille_Contenant_Point_Mth_Co((*noeuds_front)[i],num_maille_depart,num_maille_depart,NBR_MAX_MAILLES_EXAMINEES,nma,flag_convexe);
00120                }
00121           }
00122           
00123      else
00124           {
00125           cout<<"Le mapping semble déja existé, interrogation sur l'existant"<<endl;
00126           }
00127           
00128      }
00129 _TEMPLATE_ _MAPPING_::Mapping(MAILLAGE * mb,NUAGENOEUD * nb,NUAGENOEUD * nf):maillage_back(mb),noeuds_back(nb),noeuds_front(nf),my_dTree(NULL)
00130      {
00131      
00132      if (!maillage_back)
00133           {
00134           cerr<<"Mapping : constructeur appelé avec Maillage Vide"<<endl;
00135           exit(-1);
00136           }
00137           
00138      if (!noeuds_back)
00139           {
00140           cerr<<"Mapping : constructeur appelé avec Nuage Noeuds Vide"<<endl;
00141           exit(-1);
00142           }
00143           
00144      mailles_back=maillage_back->DONNE_POINTEUR_NUAGEMAILLE();
00145      
00146      CB=new Coordonnees_Barycentriques<NUAGEMAILLE,NUAGENOEUD,NOEUD,DIMENSION>(mailles_back,noeuds_back);
00147 
00148      // TEST REDONDANCE
00149      /*
00150      int nnb=noeuds_back->SIZE();
00151      if (nnb<20000) 
00152           {
00153           cout<<"MAPPING : VERIFICATION REDONDANCES DANS NUAGE NOEUD BACK"<<endl;
00154           noeuds_back->affiche();
00155           int i,j;       
00156           vector<int> redondance(nnb,0);
00157           for (i=0;i<nnb;i++) 
00158                {
00159                for (j=i+1;j<nnb;j++) if ((*noeuds_back)[i]==(*noeuds_back)[j]) 
00160                     {
00161                     cerr<<"Le Noeud "<<j<<" est identique au Noeud "<<i<<endl;
00162                     exit(-1);
00163                     }
00164                }
00165           cout<<"FIN TEST VERIFICATION REDONDANCES"<<endl;
00166           }
00167      // FIN TEST */
00168      
00169      my_dTree=new dTree<NOEUD,NUAGENOEUD,DIMENSION>(noeuds_back);
00170 
00171      }
00172 // Renvoie :
00173 //     1 si le point est intérieur
00174 //    -1 si le point est extérieur à la maille via uniquement des faces qui ne sont pas au bord
00175 //    -2 si le point est extérieur à la maille par au moins une face de bord
00176 // Et modifie etat_face de telle sorte que :
00177 // etat_face[i] = -1 s'il n'existe pas de voisin via la face i
00178 // etat_face[i] =  0 si le point est intérieur via la face i et que le voisin i existe
00179 // etat_face[i] =  1 si le point est extérieur via la face i et que le voisin i existe
00180 _TEMPLATE_ int _MAPPING_::Donne_Directions(int num_maille,const NOEUD &n,int etat_face[NBR_FACES_MAX])
00181      {
00182      vector<double> ef=CB->Donne_Pseudo_Coord_Baryc(num_maille,n);
00183      int etat_int=VRAI;
00184      int etat_ext_bord=FAUX;
00185      int tf,tv,tb;
00186      int nbr_faces=(*mailles_back)[num_maille].DONNE_NBR_FACES();
00187      for (int i=0;i<nbr_faces;i++)
00188           {
00189           tf=(ef[i]<0);
00190           tv=(maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,i)==UNDEFINED);
00191           tb=(maillage_back->EST_AU_BORD_FACE_DE_MAILLE(num_maille,i));
00192           if (tf) // extérieur
00193                {
00194                etat_int=FAUX;
00195                if (tb) etat_ext_bord=VRAI;
00196                }
00197           if (tv) etat_face[i]=-1; // ya pas de voisin
00198           else
00199                {
00200                if (tf) etat_face[i]=1;
00201                else etat_face[i]=0;
00202                }
00203           }
00204      if (etat_int) return INTERIEUR;
00205      if (etat_ext_bord) return EXTERIEUR_AU_BORD;
00206      return EXTERIEUR_AU_MILIEU;
00207      }
00208 _TEMPLATE_ int _MAPPING_::Trouve_Maille_Contenant_Point_Mth_Co(const NOEUD &n,int num_maille,int num_maille_interdit,int max_loop,int &nbr_mailles_examinees,int flag_convexe)
00209      {
00210 
00211      int etat_face[NBR_FACES_MAX];
00212      int i,tmp,nbr_rnd;
00213      int indirection[NBR_FACES_MAX];
00214      int ind_reel;
00215      int num_reel;
00216      int new_num=UNDEFINED;
00217      
00218      int test=Donne_Directions(num_maille,n,etat_face);
00219      
00220      int nbr_faces=maillage_back->DONNE_NBR_FACES_MAILLE(num_maille);
00221      
00222      if ( test != INTERIEUR ) { // EAP, for PAL11458
00223           // check neighbors
00224           int etat_face_for_check[NBR_FACES_MAX];
00225           for (i=0;i<nbr_faces;i++) {
00226             int num_neighbor=maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,i);
00227             if ( num_neighbor != UNDEFINED &&
00228                  Donne_Directions(num_neighbor,n,etat_face_for_check) == INTERIEUR )
00229               return num_neighbor;
00230             indirection[i]=i;
00231           }
00232      }
00233      
00234      nbr_mailles_examinees=0;
00235      
00236      while (nbr_mailles_examinees<max_loop)
00237           {
00238           if (test==INTERIEUR) 
00239                {
00240                return num_maille;
00241                }
00242           if ((test==EXTERIEUR_AU_BORD)&&(flag_convexe)) 
00243                {
00244                return 2*UNDEFINED;
00245                }
00246           nbr_mailles_examinees++;
00247           for (i=0;i<nbr_faces;i++)
00248                {
00249                tmp=indirection[i];
00250                nbr_rnd=rand()%nbr_faces;
00251                indirection[i]=indirection[nbr_rnd];
00252                indirection[nbr_rnd]=tmp;
00253                }
00254           for (i=0;(i<nbr_faces)&&(new_num==UNDEFINED);i++) 
00255                {
00256                ind_reel=indirection[i];
00257                num_reel=maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,ind_reel);
00258                if ((etat_face[ind_reel]==1)&&(num_reel!=num_maille_interdit)) 
00259                     {
00260                     new_num=num_reel;
00261                     }
00262                }
00263           for (i=0;(i<nbr_faces)&&(new_num==UNDEFINED);i++) 
00264                {
00265                ind_reel=indirection[i];
00266                num_reel=maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,ind_reel);
00267                if ((etat_face[ind_reel]==0)&&(num_reel!=num_maille_interdit)) 
00268                     {
00269                     new_num=num_reel;
00270                     }
00271                }
00272           if (new_num==UNDEFINED) 
00273                {
00274                new_num=num_maille_interdit;
00275                }
00276           num_maille_interdit=num_maille;
00277           num_maille=new_num;
00278           new_num=UNDEFINED;
00279           test=Donne_Directions(num_maille,n,etat_face);
00280           }
00281      return UNDEFINED;
00282      }
00283 
00284 #undef _TEMPLATE_
00285 #undef _MAPPING_
00286 
00287 #endif