src/MEDMEM/MEDMEM_DriverTools.cxx

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