src/MEDMEM/MEDMEM_DriverTools.cxx

   1 #include "MEDMEM_DriverTools.hxx"
   2 #include "MEDMEM_STRING.hxx"
   3 #include "MEDMEM_Exception.hxx"
   4 #include "MEDMEM_Mesh.hxx"
   5 #include "MEDMEM_Group.hxx"
   6 #include <iomanip>
   7 #include <algorithm>
   8
   9 using namespace std;
  10 using namespace MEDMEM;
  11 using namespace MED_EN;
  12
  13 // Cet opérateur permet d'ordonner les mailles dans un set suivant l'ordre requis par MED
  14 bool _maille::operator < (const _maille& ma) const
  15 {
  16     // si le type géométrique differe, la comparaison est basée dessus
  17     // sinon on se base sur une comparaison des numéros de sommets
  18     if(geometricType==ma.geometricType)
  19     {
  20         // construction de deux vecteur temporaire contenant les numeros de sommets
  21         // pour faire le tri et les comparaisons
  22         size_t l=sommets.size();
  23         std::vector<int> v1(l);
  24         std::vector<int> v2(l);
  25         for (unsigned int i=0; i!=l; ++i)
  26           {
  27             v1[i]=sommets[i]->second.number;
  28             v2[i]=ma.sommets[i]->second.number;
  29           }
  30         std::sort(v1.begin(), v1.end());
  31         std::sort(v2.begin(), v2.end());
  32         for(std::vector<int>::const_iterator i1=v1.begin(), i2=v2.begin(); i1!=v1.end(); ++i1, ++i2)
  33             if(*i1 != *i2)
  34                 return *i1 < *i2;
  35         return false; // cas d'égalité
  36     }
  37     else
  38         return geometricType<ma.geometricType;
  39 };
  40
  41
  42 // retourne l'entité d'une maille en fonction de la dimension du maillage.
  43 MED_EN::medEntityMesh _maille::getEntity(const int meshDimension) const throw (MEDEXCEPTION)
  44 {
  45   const char * LOC = "_maille::getEntity(const int meshDimension)" ;
  46   BEGIN_OF(LOC);
  47
  48   int mailleDimension = this->dimension();
  49   medEntityMesh entity;
  50   if (mailleDimension == meshDimension)
  51     entity = MED_CELL;
  52   else
  53     switch (mailleDimension)
  54       {
  55       case 0 :
  56         entity = MED_NODE;
  57         break;
  58       case 1 :
  59         entity = MED_EDGE;
  60         break;
  61       case 2 :
  62         entity = MED_FACE;
  63         break;
  64       default :
  65         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Impossible de determiner l'entite de la maille."));
  66       }
  67 return entity;
  68
  69 END_OF(LOC);
  70 };
  71
  72 std::ostream& MEDMEM::operator << (std::ostream& os, const _maille& ma)
  73 {
  74     os << "maille " << ma.ordre << " (" << ma.geometricType << ") : < ";
  75     std::vector< std::map<int,_noeud>::iterator >::const_iterator i=ma.sommets.begin();
  76     os << (*i++)->second.number ;
  77     for( ; i!=ma.sommets.end(); ++i)
  78         os << ", " << (*i)->second.number;
  79     os << " >";
  80     return os;
  81 }
  82
  83 std::ostream& MEDMEM::operator << (std::ostream& os, const _groupe& gr)
  84 {
  85     os << "--- Groupe " << gr.nom << " --- " << std::endl ;
  86     os << " -> liste des sous-groupes : ";
  87     for( std::vector<int>::const_iterator i=gr.groupes.begin(); i!=gr.groupes.end(); ++i)
  88         os << *i << " ";
  89     os << std::endl << " -> liste des mailles : " << std::endl;
  90     std::set< std::set< _maille, std::less< _maille >,
  91     std::allocator< _maille > >::iterator,
  92     _mailleIteratorCompare,
  93     std::allocator< std::set< _maille, std::less< _maille >,
  94     std::allocator< _maille > >::iterator > >::const_iterator i ;
  95     for( i=gr.mailles.begin(); i!=gr.mailles.end(); i++)
  96         os << "    " << *(*i) << std::endl;
  97     return os;
  98 }
  99
 100 std::ostream& MEDMEM::operator << (std::ostream& os, const _noeud& no)
 101 {
 102     os << "noeud " << no.number << " : < ";
 103     std::vector<double>::const_iterator i=no.coord.begin();
 104     os << *i++ ;
 105     for( ; i!=no.coord.end(); ++i)
 106         os << ", " << *i;
 107     os << " >";
 108     return os;
 109 }
 110
 111 std::ostream& MEDMEM::operator << (std::ostream& os, const _intermediateMED& mi)
 112 {
 113     os << "Set des mailles : " << std::endl;
 114     for( std::set<_maille>::const_iterator i=mi.maillage.begin(); i!=mi.maillage.end(); ++i)
 115         os << *i << std::endl;
 116
 117     os << std::endl << "Vector des groupes : " << std::endl;
 118    // for( std::vector<_groupe>::const_iterator i=mi.groupes.begin(); i!=mi.groupes.end(); ++i)
 119     for (unsigned int i=0; i!=mi.groupes.size(); i++)
 120       os << std::setw(3) << i << " " << mi.groupes[i] << std::endl;
 121
 122     os << std::endl << "map des noeuds : " << std::endl;
 123     for( std::map<int,_noeud>::const_iterator i=mi.points.begin(); i!=mi.points.end(); ++i)
 124         os << i->second << std::endl;
 125     return os;
 126 }
 127
 128 void _intermediateMED::numerotationMaillage()
 129 {
 130     // numerotation des mailles par entité
 131     int i_maille=0;
 132     std::set<_maille>::iterator i=maillage.begin();
 133     int dimension=i->dimension();
 134     for( ; i!=maillage.end(); ++i)
 135     {
 136         if ( dimension!=i->dimension() ) // on change d'entite
 137         {
 138             dimension=i->dimension();
 139             i_maille=0;
 140         }
 141         (*i).ordre=++i_maille;
 142     }
 143 }
 144
 145 void _intermediateMED::numerotationPoints()
 146 {
 147     // Fonction de renumerotation des noeuds (necessaire quand il y a des trous dans la numerotation.
 148     int i_noeud=0;
 149     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i)
 150         i->second.number = ++i_noeud ;
 151 }
 152
 153
 154 /*!
 155  * \if developper
 156  * create a MED COORDINATE from the intermediate structure.
 157  * \endif
 158  */
 159 COORDINATE *
 160 _intermediateMED::getCoordinate(const string & coordinateSystem)
 161 {
 162     const medModeSwitch mode=MED_FULL_INTERLACE;
 163     int spaceDimension=points.begin()->second.coord.size();
 164     int numberOfNodes=points.size();
 165
 166
 167     // creation du tableau des coordonnees en mode MED_FULL_INTERLACE
 168     double * coord = new double[spaceDimension*numberOfNodes];
 169     int k=0;
 170     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i, ++k)
 171         std::copy(i->second.coord.begin(), i->second.coord.end(), coord+k*spaceDimension);
 172
 173     // creation de l'objet COORDINATE
 174     COORDINATE * coordinate = new COORDINATE(spaceDimension, numberOfNodes, mode);
 175     coordinate->setCoordinates(mode,coord);
 176     delete [] coord;
 177     coordinate->setCoordinatesSystem(coordinateSystem);
 178     return coordinate;
 179 }
 180
 181
 182   /*!
 183    * \if developper
 184    * create a MED CONNECTIVITY from the intermediate structure.
 185    * for memory optimisation, clear the intermediate structure (the "maillage" set )
 186    * \endif
 187    */
 188 CONNECTIVITY *
 189 _intermediateMED::getConnectivity()
 190 {
 191     const char * LOC = "_intermediateMED::getConnectivity() : ";
 192     BEGIN_OF(LOC);
 193     int numberOfNodes=points.size(); // number of nodes in the mesh
 194     int numberOfTypes=0;
 195     medEntityMesh entity;
 196     medGeometryElement * types=NULL; // pointeurs pour allouer les structures MED
 197     int * count=NULL;
 198     int * connectivity=NULL;
 199     CONNECTIVITY *Connectivity, *Constituent;
 200
 201     medGeometryElement type=0; // variables de travail
 202     int nbtype=0;
 203     int dimension=0;
 204     bool first = true;
 205
 206     std::vector<medGeometryElement> vtype; // tableau de travail : stockage des types geometriques pour UNE entite
 207     std::vector<int> vcount; // tableau de travail : nombre d'elements pour chaque type geometrique de vtype
 208
 209     std::set<_maille>::const_iterator i=maillage.begin(); // iterateurs sur les mailles
 210     std::set<_maille>::const_iterator j=maillage.begin();
 211
 212     // renumerote les points de 1 a n (pour le cas ou certains points ne sont pas presents dans le maillage d'origine)
 213     numerotationPoints();
 214
 215     do
 216     {
 217         // boucle sur les entites de la structure maillage :
 218         //   - on parcourt une premiere fois avec i les mailles d'une entite pour récupérer
 219         //     des infos sur les types geometriques, leur nombre et le nombre d'elements.
 220         //   - on alloue la connectivite
 221         //   - on parcourt une deuxieme fois avec j pour lire les noeuds.
 222
 223
 224         type=i->geometricType; // init boucle for
 225         dimension=i->dimension();
 226         nbtype=0;
 227         vtype.push_back(type);
 228         // Boucle sur i de parcours des mailles d'une entite
 229         // Une entite se termine lorsqu'on atteint la fin de maillage ou lorsque la dimension des mailles change
 230         for( ; i!=maillage.end() && dimension==i->dimension() ; ++i)
 231         {
 232             if (type != i->geometricType) // si changement de type geometrique
 233             {
 234                 vcount.push_back(nbtype); // stocke le nombre d'elements de type nbtype
 235                 nbtype=0;
 236                 type=i->geometricType;
 237                 vtype.push_back(type); // stocke le nouveau type geometrique rencontre
 238             }
 239
 240             ++nbtype;
 241         }
 242         vcount.push_back(nbtype); // n'a pas été stocké dans la boucle
 243
 244
 245         // Pour l'entite qu'on vient de parcourir, allocation des tableau et creation de la connectivite
 246 //      cout << "Dimension = " << dimension << endl;
 247 //      cout << "Nombre de type geometriques : " << vtype.size() << endl;
 248 //      for (unsigned k=0; k!=vtype.size(); ++k )
 249 //          cout << "  -> " << vtype[k] << " : " << vcount[k] << endl;
 250
 251         numberOfTypes=vtype.size(); // nombre de types de l'entite
 252
 253         if ( i==maillage.end() ) // cas de l'entite de plus haut niveau
 254             entity=MED_CELL;
 255         else if (dimension==2 )
 256             entity=MED_FACE;
 257         else if (dimension==1 )
 258             entity=MED_EDGE;
 259         else if (dimension==0 )
 260             entity=MED_NODE;
 261
 262         Connectivity = new CONNECTIVITY(numberOfTypes,entity);
 263         Connectivity->setEntityDimension(dimension);
 264         Connectivity->setNumberOfNodes(numberOfNodes);
 265
 266         types = new medGeometryElement[numberOfTypes];
 267         std::copy(vtype.begin(),vtype.end(),types);
 268         Connectivity->setGeometricTypes(types,entity);
 269         delete [] types;
 270
 271         count = new int[numberOfTypes+1];
 272         count[0]=1;
 273         for (unsigned int k=0; k!=vcount.size(); ++k )
 274           count[k+1]=count[k]+vcount[k];
 275         Connectivity->setCount (count, entity);
 276         delete [] count;
 277
 278         for (int k=0; k!=numberOfTypes; ++k )
 279           {
 280             // pour chaque type géometrique k, copie des sommets dans connectivity et set dans Connectivity
 281             int nbSommetsParMaille = j->sommets.size();
 282             int nbSommets = vcount[k] * j->sommets.size();
 283             connectivity = new int[ nbSommets ];
 284             for (int l=0; l!=vcount[k]; ++l, ++j)
 285             {
 286                 for ( unsigned n=0; n != j->sommets.size(); ++n)
 287                     connectivity[nbSommetsParMaille*l+n] = j->sommets[n]->second.number;
 288                 //CCRT maillage.erase(j);    ; // dangereux, mais optimise la mémoire consommée!
 289             }
 290
 291             Connectivity->setNodal  (connectivity, entity, vtype[k]);
 292             delete [] connectivity;
 293           }
 294
 295         if (i!=j)
 296             throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Erreur de lecture des mailles ."));
 297
 298         if ( ! first)
 299             Connectivity->setConstituent (Constituent);
 300         else
 301             first = false;
 302         Constituent = Connectivity; // stocke Connectivity pour utilisation dans setConstituent lors de la boucle suivante
 303
 304         vtype.clear();
 305         vcount.clear();
 306
 307     }
 308     while ( i != maillage.end() );
 309
 310     END_OF(LOC);
 311     return Connectivity;
 312 }
 313
 314
 315   /*!
 316    * \if developper
 317    * fill the arguments vector of groups from the intermediate structure.
 318    * This function must be called before getConnectivity()
 319    * \endif
 320    */
 321 void
 322 _intermediateMED::getGroups(std::vector<GROUP *> & _groupCell, std::vector<GROUP *> & _groupFace, std::vector<GROUP *> & _groupEdge, std::vector<GROUP *> & _groupNode, MESH * _ptrMesh)
 323 {
 324     const char * LOC = "_intermediateMED::getGroups() : ";
 325     BEGIN_OF(LOC);
 326     if (maillage.size() == 0)
 327         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Erreur : il n'y a plus de mailles (appeler getConnectivity apres getGroups!)."));
 328
 329     int dimension_maillage=maillage.rbegin()->dimension();
 330
 331     numerotationMaillage(); // Renumerotation des mailles par entite
 332
 333     for (unsigned int i=0; i!=this->groupes.size(); ++i)
 334     {
 335         // si le groupe est vide, ou si le groupe n'est pas nommé : on passe au suivant
 336         if (groupes[i].mailles.empty() || groupes[i].nom.empty())
 337             continue;
 338
 339         int nb_geometric_types=1;
 340         std::set< std::set< _maille, std::less< _maille >,
 341         std::allocator< _maille > >::iterator,
 342         _mailleIteratorCompare,
 343         std::allocator< std::set< _maille, std::less<_maille >,
 344         std::allocator< _maille > >::iterator > >::iterator j ;
 345         j=groupes[i].mailles.begin();
 346         // initialise groupe_entity a l'entite de la premiere maille du groupe
 347         medEntityMesh groupe_entity = (**j).getEntity(dimension_maillage);
 348         medGeometryElement geometrictype=(**j).geometricType;
 349
 350         //Parcours des mailles (a partir de la deuxieme) pour compter les types geometriques
 351         for ( ++j ; j!=groupes[i].mailles.end(); ++j )
 352         {
 353             //Compte nombre de types geometriques
 354             if ( (**j).geometricType != geometrictype ) // si on change de type geometrique
 355             {
 356                 nb_geometric_types++;
 357                 geometrictype=(**j).geometricType;
 358             }
 359
 360             //Test si groupe valide : le groupe doit pointer vers des entites de meme dimension
 361             if ((**j).dimension() != dimension_maillage)
 362                 continue;
 363         }
 364
 365         // le groupe est valide -> on le traite
 366         MED_EN::medGeometryElement * tab_types_geometriques = new MED_EN::medGeometryElement[nb_geometric_types];
 367         int * tab_index_types_geometriques = new int[nb_geometric_types+1];
 368         int * tab_numeros_elements = new int[groupes[i].mailles.size()];
 369         int * tab_nombres_elements = new int[nb_geometric_types];
 370
 371         //Remplit tableaux entree des methodes set
 372         int indice_mailles=0;
 373         j=groupes[i].mailles.begin();
 374         geometrictype=(**j).geometricType;
 375         tab_index_types_geometriques[0]=1;
 376         int indice_types_geometriques=1;
 377         tab_types_geometriques[0]=geometrictype;
 378         //parcours des mailles du groupe
 379         for (  ; j!=groupes[i].mailles.end(); ++j , ++indice_mailles)
 380         {
 381             tab_numeros_elements[indice_mailles]=((**j).ordre);
 382             if ((**j).geometricType != geometrictype)
 383             {
 384                 tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 385                 geometrictype=(**j).geometricType;
 386                 tab_types_geometriques[indice_types_geometriques]=geometrictype;
 387                 ++indice_types_geometriques;
 388             }
 389         }
 390         tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
 391         for (int k=0; k != nb_geometric_types; ++k)
 392         {
 393             tab_nombres_elements[k] = tab_index_types_geometriques[k+1]-tab_index_types_geometriques[k];
 394         }
 395
 396         //Determination type entite du groupe
 397         vector <GROUP *> * vect_group;
 398         switch ( groupe_entity )
 399         {
 400             case MED_CELL :
 401                 vect_group= & _groupCell;
 402                 break;
 403             case MED_FACE :
 404                 vect_group= & _groupFace;
 405                 break;
 406             case MED_EDGE :
 407                 vect_group= & _groupEdge;
 408                 break;
 409             case MED_NODE :
 410                 vect_group= & _groupNode;
 411                 break;
 412         }
 413         //Creation nouveau groupe MED
 414         GROUP * new_group = new GROUP();
 415         //Appel methodes set
 416         new_group->setTotalNumberOfElements(groupes[i].mailles.size());
 417         new_group->setName(groupes[i].nom);
 418         new_group->setMesh(_ptrMesh);
 419         new_group->setNumberOfGeometricType(nb_geometric_types);
 420         new_group->setGeometricType(tab_types_geometriques);
 421         new_group->setNumberOfElements(tab_nombres_elements);
 422         new_group->setNumber(tab_index_types_geometriques,tab_numeros_elements);
 423         new_group->setEntity(groupe_entity);
 424         new_group->setAll(groupes[i].mailles.size() == maillage.size());
 425         vect_group->push_back(new_group);
 426         delete [] tab_types_geometriques;
 427         delete [] tab_index_types_geometriques;
 428         delete [] tab_numeros_elements;
 429         delete [] tab_nombres_elements;
 430     }
 431
 432     END_OF(LOC);
 433 }
 434
 435 /////