src/INTERPOLATION/MEDMEM_Mapping.hxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 #ifndef MEDMEM_MAPPING_HXX
  23 #define MEDMEM_MAPPING_HXX
  24
  25 #include "MEDMEM_MappingTools.hxx"
  26 #include "MEDMEM_dTree.hxx"
  27
  28 #define NBR_MAX_MAILLES_EXAMINEES 100
  29
  30 #ifndef  NBR_FACES_MAX
  31 #define NBR_FACES_MAX 6
  32 #endif
  33
  34 #define _TEMPLATE_ template <class MAILLAGE, class NUAGEMAILLE, class NUAGENOEUD, class NOEUD, int DIMENSION>
  35 #define _MAPPING_ Mapping<MAILLAGE,NUAGEMAILLE,NUAGENOEUD,NOEUD,DIMENSION>
  36
  37
  38 //////////////////////////////////////////////////////////////////
  39 ///                                                            ///
  40 ///                        DECLARATIONS                        ///
  41 ///                                                            ///
  42 //////////////////////////////////////////////////////////////////
  43
  44 /*********************************************************/
  45 /*                                                       */
  46 /*                   Classe  Mapping                     */
  47 /*                                                       */
  48 /*********************************************************/
  49
  50 // ATTENTION LE NUAGE DE NOEUD EST SUPPOSE NON REDONDANT ET AUCUNE VERIFICATION N'EST FAITE !
  51
  52 template <class MAILLAGE, class NUAGEMAILLE, class NUAGENOEUD, class NOEUD, int DIMENSION> class Mapping
  53 {
  54 protected :
  55         MAILLAGE * maillage_back;
  56         NUAGEMAILLE * mailles_back;
  57         NUAGENOEUD * noeuds_back;
  58         NUAGENOEUD * noeuds_front;
  59         Coordonnees_Barycentriques<NUAGEMAILLE,NUAGENOEUD,NOEUD,DIMENSION> * CB;
  60         dTree<NOEUD,NUAGENOEUD,DIMENSION> * my_dTree;
  61         vector<int> resultat_mapping;
  62         vector<int> point_le_plus_proche;
  63
  64 public :
  65         Mapping():maillage_back(NULL),mailles_back(NULL),noeuds_back(NULL),noeuds_front(NULL),CB(NULL),my_dTree(NULL) {}
  66         Mapping(MAILLAGE * mb,NUAGENOEUD * nb,NUAGENOEUD * nf); // le dTree est crée à l'initialisation, par contre, le mapping lui meme doit etre invoqué
  67         ~Mapping() {if (CB) delete CB;if (my_dTree) delete my_dTree;}
  68         dTree<NOEUD,NUAGENOEUD,DIMENSION> * Donne_dTree() {return my_dTree;}
  69         enum { INTERIEUR = 1, EXTERIEUR_AU_MILIEU = -1, EXTERIEUR_AU_BORD = -2 };
  70         int Donne_Directions(int num_maille,const NOEUD &n,int etat_face[NBR_FACES_MAX]);
  71         // Méthode interne de localisation
  72         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);
  73         void Cree_Mapping(int flag_convexe=0);                                             // SUPPOSE NON CONVEXE PAR DEFAUT
  74         void Cree_Mapping(NUAGENOEUD * nf, int flag_convexe=0);                            // SUPPOSE NON CONVEXE PAR DEFAUT
  75         inline int operator[](int i) const {return resultat_mapping[i];}                   // Renvoie la valeur mappé, si le mapping a été fait, sinon, n'importe quoi
  76         inline vector<int> & Get_Mapping() {return resultat_mapping;}                        // Renvoie le vector contenant le mapping
  77         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
  78         inline int Exist_dTree() const {return (my_dTree);}                                // Teste si le dTree existe
  79         void affiche()
  80                 {
  81                 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;
  82                 }
  83 };
  84
  85 //////////////////////////////////////////////////////////////////
  86 ///                                                            ///
  87 ///                            CODE                            ///
  88 ///                                                            ///
  89 //////////////////////////////////////////////////////////////////
  90
  91 _TEMPLATE_ void _MAPPING_::Cree_Mapping(NUAGENOEUD * nf, int flag_convexe)
  92         {
  93         noeuds_front=nf;
  94         Cree_Mapping(flag_convexe);
  95         }
  96
  97 _TEMPLATE_ void _MAPPING_::Cree_Mapping(int flag_convexe)
  98         {
  99
 100         if (resultat_mapping.size()==0)
 101                 {
 102                 if (!noeuds_front)
 103                          {
 104                          cerr<<"Mapping : Pas de noeuds à mapper !"<<endl;
 105                          exit(-1);
 106                          }
 107
 108                 int i;
 109                 int nbr_noeuds=noeuds_front->SIZE();
 110                 int num_maille_depart;
 111                 int nma=0;
 112                 resultat_mapping     = vector<int>(nbr_noeuds,MED_UNDEFINED);
 113                 point_le_plus_proche = vector<int>(nbr_noeuds,MED_UNDEFINED);
 114
 115                 // noeuds_back->affiche();
 116
 117                 for (i=0;i<nbr_noeuds;i++)
 118                         {
 119                         point_le_plus_proche[i]=my_dTree->trouve_plus_proche_point((*noeuds_front)[i]);
 120                         num_maille_depart=maillage_back->DONNE_PREMIERE_MAILLE_CONTENANT(point_le_plus_proche[i]);
 121                         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);
 122                         }
 123                 }
 124
 125         else
 126                 {
 127                 cout<<"Le mapping semble déja existé, interrogation sur l'existant"<<endl;
 128                 }
 129
 130         }
 131 _TEMPLATE_ _MAPPING_::Mapping(MAILLAGE * mb,NUAGENOEUD * nb,NUAGENOEUD * nf):maillage_back(mb),noeuds_back(nb),noeuds_front(nf),my_dTree(NULL)
 132         {
 133
 134         if (!maillage_back)
 135                 {
 136                 cerr<<"Mapping : constructeur appelé avec Maillage Vide"<<endl;
 137                 exit(-1);
 138                 }
 139
 140         if (!noeuds_back)
 141                 {
 142                 cerr<<"Mapping : constructeur appelé avec Nuage Noeuds Vide"<<endl;
 143                 exit(-1);
 144                 }
 145
 146         mailles_back=maillage_back->DONNE_POINTEUR_NUAGEMAILLE();
 147
 148         CB=new Coordonnees_Barycentriques<NUAGEMAILLE,NUAGENOEUD,NOEUD,DIMENSION>(mailles_back,noeuds_back);
 149
 150         // TEST REDONDANCE
 151         /*
 152         int nnb=noeuds_back->SIZE();
 153         if (nnb<20000)
 154                 {
 155                 cout<<"MAPPING : VERIFICATION REDONDANCES DANS NUAGE NOEUD BACK"<<endl;
 156                 noeuds_back->affiche();
 157                 int i,j;
 158                 vector<int> redondance(nnb,0);
 159                 for (i=0;i<nnb;i++)
 160                         {
 161                         for (j=i+1;j<nnb;j++) if ((*noeuds_back)[i]==(*noeuds_back)[j])
 162                                 {
 163                                 cerr<<"Le Noeud "<<j<<" est identique au Noeud "<<i<<endl;
 164                                 exit(-1);
 165                                 }
 166                         }
 167                 cout<<"FIN TEST VERIFICATION REDONDANCES"<<endl;
 168                 }
 169         // FIN TEST */
 170
 171         my_dTree=new dTree<NOEUD,NUAGENOEUD,DIMENSION>(noeuds_back);
 172
 173         }
 174 // Renvoie :
 175 //     1 si le point est intérieur
 176 //    -1 si le point est extérieur à la maille via uniquement des faces qui ne sont pas au bord
 177 //    -2 si le point est extérieur à la maille par au moins une face de bord
 178 // Et modifie etat_face de telle sorte que :
 179 // etat_face[i] = -1 s'il n'existe pas de voisin via la face i
 180 // etat_face[i] =  0 si le point est intérieur via la face i et que le voisin i existe
 181 // etat_face[i] =  1 si le point est extérieur via la face i et que le voisin i existe
 182 _TEMPLATE_ int _MAPPING_::Donne_Directions(int num_maille,const NOEUD &n,int etat_face[NBR_FACES_MAX])
 183         {
 184         vector<double> ef=CB->Donne_Pseudo_Coord_Baryc(num_maille,n);
 185         int etat_int=VRAI;
 186         int etat_ext_bord=FAUX;
 187         int tf,tv,tb;
 188         int nbr_faces=(*mailles_back)[num_maille].DONNE_NBR_FACES();
 189         for (int i=0;i<nbr_faces;i++)
 190                 {
 191                 tf=(ef[i]<0);
 192                 tv=(maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,i)==MED_UNDEFINED);
 193                 tb=(maillage_back->EST_AU_BORD_FACE_DE_MAILLE(num_maille,i));
 194                 if (tf) // extérieur
 195                         {
 196                         etat_int=FAUX;
 197                         if (tb) etat_ext_bord=VRAI;
 198                         }
 199                 if (tv) etat_face[i]=-1; // ya pas de voisin
 200                 else
 201                         {
 202                         if (tf) etat_face[i]=1;
 203                         else etat_face[i]=0;
 204                         }
 205                 }
 206         if (etat_int) return INTERIEUR;
 207         if (etat_ext_bord) return EXTERIEUR_AU_BORD;
 208         return EXTERIEUR_AU_MILIEU;
 209         }
 210 _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)
 211         {
 212
 213         int etat_face[NBR_FACES_MAX];
 214         int i,tmp,nbr_rnd;
 215         int indirection[NBR_FACES_MAX];
 216         int ind_reel;
 217         int num_reel;
 218         int new_num=MED_UNDEFINED;
 219
 220         int test=Donne_Directions(num_maille,n,etat_face);
 221
 222         int nbr_faces=maillage_back->DONNE_NBR_FACES_MAILLE(num_maille);
 223
 224         if ( test != INTERIEUR ) { // EAP, for PAL11458
 225           // check neighbors
 226           int etat_face_for_check[NBR_FACES_MAX];
 227           for (i=0;i<nbr_faces;i++) {
 228             int num_neighbor=maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,i);
 229             if ( num_neighbor != MED_UNDEFINED &&
 230                  Donne_Directions(num_neighbor,n,etat_face_for_check) == INTERIEUR )
 231               return num_neighbor;
 232             indirection[i]=i;
 233           }
 234         }
 235
 236         nbr_mailles_examinees=0;
 237
 238         while (nbr_mailles_examinees<max_loop)
 239                 {
 240                 if (test==INTERIEUR)
 241                         {
 242                         return num_maille;
 243                         }
 244                 if ((test==EXTERIEUR_AU_BORD)&&(flag_convexe))
 245                         {
 246                         return 2*MED_UNDEFINED;
 247                         }
 248                 nbr_mailles_examinees++;
 249                 for (i=0;i<nbr_faces;i++)
 250                         {
 251                         tmp=indirection[i];
 252                         nbr_rnd=rand()%nbr_faces;
 253                         indirection[i]=indirection[nbr_rnd];
 254                         indirection[nbr_rnd]=tmp;
 255                         }
 256                 for (i=0;(i<nbr_faces)&&(new_num==MED_UNDEFINED);i++)
 257                         {
 258                         ind_reel=indirection[i];
 259                         num_reel=maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,ind_reel);
 260                         if ((etat_face[ind_reel]==1)&&(num_reel!=num_maille_interdit))
 261                                 {
 262                                 new_num=num_reel;
 263                                 }
 264                         }
 265                 for (i=0;(i<nbr_faces)&&(new_num==MED_UNDEFINED);i++)
 266                         {
 267                         ind_reel=indirection[i];
 268                         num_reel=maillage_back->DONNE_VOISIN_DE_MAILLE(num_maille,ind_reel);
 269                         if ((etat_face[ind_reel]==0)&&(num_reel!=num_maille_interdit))
 270                                 {
 271                                 new_num=num_reel;
 272                                 }
 273                         }
 274                 if (new_num==MED_UNDEFINED)
 275                         {
 276                         new_num=num_maille_interdit;
 277                         }
 278                 num_maille_interdit=num_maille;
 279                 num_maille=new_num;
 280                 new_num=MED_UNDEFINED;
 281                 test=Donne_Directions(num_maille,n,etat_face);
 282                 }
 283         return MED_UNDEFINED;
 284         }
 285
 286 #undef _TEMPLATE_
 287 #undef _MAPPING_
 288
 289 #endif