src/HOMARD_I/HomardMedCommun.cxx

   1 // Copyright (C) 2011-2015  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "HomardMedCommun.h"
  21
  22 #include <iostream>
  23 #include <cstdlib>
  24 #include <cmath>
  25 #include <cstring>
  26 #include <algorithm>
  27
  28 extern "C"
  29 {
  30 #include <med.h>
  31 }
  32
  33 #include "utilities.h"
  34 // =======================================================================
  35 int MEDFileExist( const char * aFile )
  36 // Retourne 1 si le fichier existe, 0 sinon
  37 // =======================================================================
  38 {
  39   int existe ;
  40   med_idt medIdt = MEDfileOpen(aFile,MED_ACC_RDONLY);
  41   if ( medIdt < 0 ) { existe = 0 ; }
  42   else              { MEDfileClose(medIdt);
  43                       existe = 1 ; }
  44   return existe ;
  45 }
  46 // =======================================================================
  47 std::set<std::string> GetListeGroupesInMedFile(const char * aFile)
  48 // =======================================================================
  49 {
  50   std::set<std::string> ListeGroupes;
  51   med_err erreur = 0 ;
  52   med_int medIdt ;
  53   while ( erreur == 0 )
  54   {
  55     //  Ouverture du fichier
  56     medIdt = MEDfileOpen(aFile,MED_ACC_RDONLY);
  57     if ( medIdt < 0 )
  58     {
  59       erreur = 1 ;
  60       break ;
  61     }
  62     // Caracteristiques du maillage
  63     char meshname[MED_NAME_SIZE+1];
  64     med_int spacedim,meshdim;
  65     med_mesh_type meshtype;
  66     char descriptionription[MED_COMMENT_SIZE+1];
  67     char dtunit[MED_SNAME_SIZE+1];
  68     med_sorting_type sortingtype;
  69     med_int nstep;
  70     med_axis_type axistype;
  71     int naxis = MEDmeshnAxis(medIdt,1);
  72     char *axisname=new char[naxis*MED_SNAME_SIZE+1];
  73     char *axisunit=new char[naxis*MED_SNAME_SIZE+1];
  74     erreur = MEDmeshInfo(medIdt,
  75                             1,
  76                             meshname,
  77                             &spacedim,
  78                             &meshdim,
  79                             &meshtype,
  80                             descriptionription,
  81                             dtunit,
  82                             &sortingtype,
  83                             &nstep,
  84                             &axistype,
  85                             axisname,
  86                             axisunit);
  87     delete[] axisname ;
  88     delete[] axisunit ;
  89     if ( erreur < 0 ) { break ; }
  90     // Nombre de familles
  91     med_int nfam ;
  92     nfam = MEDnFamily(medIdt,meshname) ;
  93     if ( nfam < 0 )
  94     {
  95       erreur = 2 ;
  96       break ;
  97     }
  98   // Lecture des caracteristiques des familles
  99     for (int i=0;i<nfam;i++)
 100     {
 101 //       Lecture du nombre de groupes
 102       med_int ngro = MEDnFamilyGroup(medIdt,meshname,i+1);
 103       if ( ngro < 0 )
 104       {
 105         erreur = 3 ;
 106         break ;
 107       }
 108 //       Lecture de la famille
 109       else if ( ngro > 0 )
 110       {
 111         char familyname[MED_NAME_SIZE+1];
 112         med_int numfam;
 113         char* gro = (char*) malloc(MED_LNAME_SIZE*ngro+1);
 114         erreur = MEDfamilyInfo(medIdt,
 115                                meshname,
 116                                i+1,
 117                                familyname,
 118                                &numfam,
 119                                gro);
 120         if ( erreur < 0 )
 121         {
 122           free(gro);
 123           break ;
 124         }
 125         // Lecture des groupes pour une famille de mailles
 126         if ( numfam < 0)
 127         {
 128           for (int j=0;j<ngro;j++)
 129           {
 130             char str2[MED_LNAME_SIZE+1];
 131             strncpy(str2,gro+j*MED_LNAME_SIZE,MED_LNAME_SIZE);
 132             str2[MED_LNAME_SIZE] = '\0';
 133             ListeGroupes.insert(std::string(str2));
 134           }
 135         }
 136         free(gro);
 137       }
 138     }
 139     break ;
 140   }
 141   // Fermeture du fichier
 142   if ( medIdt > 0 ) MEDfileClose(medIdt);
 143
 144   return ListeGroupes;
 145 }
 146
 147 // =======================================================================
 148 std::vector<double> GetBoundingBoxInMedFile(const char * aFile)
 149 // =======================================================================
 150 {
 151 // Le vecteur en retour contiendra les informations suivantes :
 152 // en position 0 et 1 Xmin, Xmax et en position 2 Dx si < 0  2D
 153 // en position 3 et 4 Ymin, Ymax et en position 5 Dy si < 0  2D
 154 // en position 6 et 7 Zmin, Zmax et en position 8 Dz si < 0  2D
 155 //  9 distance max dans le maillage
 156
 157   std::vector<double> LesExtremes;
 158   med_err erreur = 0 ;
 159   med_int medIdt ;
 160   while ( erreur == 0 )
 161   {
 162     //  Ouverture du fichier
 163     medIdt = MEDfileOpen(aFile,MED_ACC_RDONLY);
 164     if ( medIdt < 0 )
 165     {
 166       erreur = 1 ;
 167       break ;
 168     }
 169     //Nombre de maillage : on ne peut en lire qu'un seul
 170     med_int numberOfMeshes = MEDnMesh(medIdt) ;
 171     if (numberOfMeshes != 1 )
 172     {
 173       erreur = 2 ;
 174       break ;
 175     }
 176     // Caracteristiques du maillage
 177     char meshname[MED_NAME_SIZE+1];
 178     med_int spacedim,meshdim;
 179     med_mesh_type meshtype;
 180     char descriptionription[MED_COMMENT_SIZE+1];
 181     char dtunit[MED_SNAME_SIZE+1];
 182     med_sorting_type sortingtype;
 183     med_int nstep;
 184     med_axis_type axistype;
 185     int naxis = MEDmeshnAxis(medIdt,1);
 186     char *axisname=new char[naxis*MED_SNAME_SIZE+1];
 187     char *axisunit=new char[naxis*MED_SNAME_SIZE+1];
 188     erreur = MEDmeshInfo(medIdt,
 189                             1,
 190                             meshname,
 191                             &spacedim,
 192                             &meshdim,
 193                             &meshtype,
 194                             descriptionription,
 195                             dtunit,
 196                             &sortingtype,
 197                             &nstep,
 198                             &axistype,
 199                             axisname,
 200                             axisunit);
 201     delete[] axisname ;
 202     delete[] axisunit ;
 203     if ( erreur < 0 ) { break ; }
 204
 205     // Nombre de noeuds
 206     med_bool chgt,trsf;
 207     med_int nnoe  = MEDmeshnEntity(medIdt,
 208                               meshname,
 209                               MED_NO_DT,
 210                               MED_NO_IT,
 211                               MED_NODE,
 212                               MED_NO_GEOTYPE,
 213                               MED_COORDINATE,
 214                               MED_NO_CMODE,
 215                               &chgt,
 216                               &trsf);
 217     if ( nnoe < 0 )
 218     {
 219       erreur =  4;
 220       break ;
 221     }
 222
 223     // Les coordonnees
 224     med_float* coo    = (med_float*) malloc(sizeof(med_float)*nnoe*spacedim);
 225
 226     erreur = MEDmeshNodeCoordinateRd(medIdt,
 227                                       meshname,
 228                                       MED_NO_DT,
 229                                       MED_NO_IT,
 230                                       MED_NO_INTERLACE,
 231                                       coo);
 232     if ( erreur < 0 )
 233     {
 234       free(coo) ;
 235       break ;
 236     }
 237
 238     // Calcul des extremes
 239     med_float xmin,xmax,ymin,ymax,zmin,zmax;
 240
 241     xmin=coo[0];
 242     xmax=coo[0];
 243     for (int i=1;i<nnoe;i++)
 244     {
 245         xmin = std::min(xmin,coo[i]);
 246         xmax = std::max(xmax,coo[i]);
 247     }
 248   //
 249     if (spacedim > 1)
 250     {
 251         ymin=coo[nnoe]; ymax=coo[nnoe];
 252         for (int i=nnoe+1;i<2*nnoe;i++)
 253         {
 254             ymin = std::min(ymin,coo[i]);
 255             ymax = std::max(ymax,coo[i]);
 256         }
 257     }
 258     else
 259     {
 260         ymin=0;
 261         ymax=0;
 262         zmin=0;
 263         zmax=0;
 264     }
 265 //
 266     if (spacedim > 2)
 267     {
 268         zmin=coo[2*nnoe]; zmax=coo[2*nnoe];
 269         for (int i=2*nnoe+1;i<3*nnoe;i++)
 270         {
 271             zmin = std::min(zmin,coo[i]);
 272             zmax = std::max(zmax,coo[i]);
 273         }
 274     }
 275     else
 276     {
 277         zmin=0;
 278         zmax=0;
 279     }
 280
 281     MESSAGE( "_______________________________________");
 282     MESSAGE( "xmin : " << xmin << " xmax : " << xmax );
 283     MESSAGE( "ymin : " << ymin << " ymax : " << ymax );
 284     MESSAGE( "zmin : " << zmin << " zmax : " << zmax );
 285     MESSAGE( "_______________________________________" );
 286     double epsilon = 1.e-6 ;
 287     LesExtremes.push_back(xmin);
 288     LesExtremes.push_back(xmax);
 289     LesExtremes.push_back(0);
 290     LesExtremes.push_back(ymin);
 291     LesExtremes.push_back(ymax);
 292     LesExtremes.push_back(0);
 293     LesExtremes.push_back(zmin);
 294     LesExtremes.push_back(zmax);
 295     LesExtremes.push_back(0);
 296
 297
 298    double max1=std::max ( LesExtremes[1] - LesExtremes[0] , LesExtremes[4] - LesExtremes[3] ) ;
 299    double max2=std::max ( max1 , LesExtremes[7] - LesExtremes[6] ) ;
 300    LesExtremes.push_back(max2);
 301
 302 // LesExtremes[0] = Xmini du maillage
 303 // LesExtremes[1] = Xmaxi du maillage
 304 // LesExtremes[2] = increment de progression en X
 305 // LesExtremes[3,4,5] : idem pour Y
 306 // LesExtremes[6,7,8] : idem pour Z
 307 // LesExtremes[9] = ecart maximal entre coordonnees
 308 // On fait un traitement pour dans le cas d'une coordonnee constante
 309 // inhiber ce cas en mettant un increment negatif
 310 //
 311     double diff = LesExtremes[1] - LesExtremes[0];
 312     if ( fabs(diff) > epsilon*max2 ) { LesExtremes[2] = diff/100.; }
 313     else                             { LesExtremes[2] = -1. ; }
 314
 315     diff = LesExtremes[4] - LesExtremes[3];
 316     if ( fabs(diff) > epsilon*max2 ) { LesExtremes[5]=diff/100.; }
 317     else                             { LesExtremes[5] = -1. ; }
 318
 319     diff = LesExtremes[7] - LesExtremes[6];
 320     if ( fabs(diff) > epsilon*max2 ) { LesExtremes[8]=diff/100.; }
 321     else                             { LesExtremes[8] = -1. ;  }
 322
 323     MESSAGE ( "_______________________________________" );
 324     MESSAGE ( "xmin : " << LesExtremes[0] << " xmax : " << LesExtremes[1] << " xincr : " << LesExtremes[2] );
 325     MESSAGE ( "ymin : " << LesExtremes[3] << " ymax : " << LesExtremes[4] << " yincr : " << LesExtremes[5] );
 326     MESSAGE ( "zmin : " << LesExtremes[6] << " zmax : " << LesExtremes[7] << " zincr : " << LesExtremes[8] );
 327     MESSAGE ( "dmax : " << LesExtremes[9] );
 328     MESSAGE ( "_______________________________________" );
 329
 330     free(coo) ;
 331     break ;
 332   }
 333   // Fermeture du fichier
 334   if ( medIdt > 0 ) MEDfileClose(medIdt);
 335
 336    return  LesExtremes;
 337 }
 338