src/HOMARD_I/HomardMedCommun.cxx

   1 // Copyright (C) 2011-2012  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.
   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_idt medIdt = MEDfileOpen(aFile,MED_ACC_RDONLY);
  52   if ( medIdt < 0 ) { return ListeGroupes; };
  53
  54   char meshname[MED_NAME_SIZE+1];
  55   med_int spacedim,meshdim;
  56   med_mesh_type meshtype;
  57   char descriptionription[MED_COMMENT_SIZE+1];
  58   char dtunit[MED_SNAME_SIZE+1];
  59   med_sorting_type sortingtype;
  60   med_int nstep;
  61   med_axis_type axistype;
  62   int naxis = MEDmeshnAxis(medIdt,1);
  63   char *axisname=new char[naxis*MED_SNAME_SIZE+1];
  64   char *axisunit=new char[naxis*MED_SNAME_SIZE+1];
  65   med_err aRet = MEDmeshInfo(medIdt,
  66                           1,
  67                           meshname,
  68                           &spacedim,
  69                           &meshdim,
  70                           &meshtype,
  71                           descriptionription,
  72                           dtunit,
  73                           &sortingtype,
  74                           &nstep,
  75                           &axistype,
  76                           axisname,
  77                           axisunit);
  78    if ( aRet < 0 ) { return ListeGroupes; };
  79
  80    med_int nfam, ngro, natt;
  81    if ((nfam = MEDnFamily(medIdt,meshname)) < 0) { return ListeGroupes; };
  82
  83   char familyname[MED_NAME_SIZE+1];
  84   med_int numfam;
  85   for (int i=0;i<nfam;i++)
  86   {
  87     if ((ngro = MEDnFamilyGroup(medIdt,meshname,i+1)) < 0)
  88     {
  89       // GERALD -- QMESSAGE BOX
  90       std::cerr << " Error : Families are unreadable" << std::endl;
  91       std::cerr << "Pb avec la famille : " << i+1 << std::endl;
  92       break;
  93     }
  94     if (ngro == 0) continue;
  95
  96     if ((natt = MEDnFamily23Attribute(medIdt,meshname,i+1)) < 0)
  97     {
  98       // GERALD -- QMESSAGE BOX
  99       std::cerr << " Error : Families are unreadable" << std::endl;
 100       std::cerr << "Pb avec la famille : " << i+1 << std::endl;
 101       break;
 102     }
 103
 104     med_int* attide = (med_int*) malloc(sizeof(med_int)*natt);
 105     med_int* attval = (med_int*) malloc(sizeof(med_int)*natt);
 106     char*    attdes = (char *)   malloc(MED_COMMENT_SIZE*natt+1);
 107     char*    gro    = (char*)    malloc(MED_LNAME_SIZE*ngro+1);
 108
 109     med_err aRet = MEDfamily23Info(medIdt,
 110                                 meshname,
 111                                 i+1,
 112                                 familyname,
 113                                 attide,
 114                                 attval,
 115                                 attdes,
 116                                 &numfam,
 117                                 gro);
 118
 119     if (aRet < 0)
 120     {
 121       // GERALD -- QMESSAGE BOX
 122       std::cerr << " Error : Families are unreadable" << std::endl;
 123       std::cerr << "Pb avec la famille : " << i+1 << std::endl;
 124         break;
 125     }
 126     free(attide);
 127     free(attval);
 128     free(attdes);
 129     if ((numfam )> 0) { continue;} // On ne garde que les familles d elts
 130
 131     for (int j=0;j<ngro;j++)
 132     {
 133           char str2[MED_LNAME_SIZE+1];
 134           strncpy(str2,gro+j*MED_LNAME_SIZE,MED_LNAME_SIZE);
 135           str2[MED_LNAME_SIZE] = '\0';
 136           ListeGroupes.insert(std::string(str2));
 137     }
 138     free(gro);
 139   }
 140   MEDfileClose(medIdt);
 141   return ListeGroupes;
 142 }
 143
 144 // =======================================================================
 145 std::vector<double> GetBoundingBoxInMedFile(const char * aFile)
 146 // =======================================================================
 147 {
 148 // Le vecteur en retour contiendra les informations suivantes :
 149 // en position 0 et 1 Xmin, Xmax et en position 2 Dx si < 0  2D
 150 // en position 3 et 4 Ymin, Ymax et en position 5 Dy si < 0  2D
 151 // en position 6 et 7 Zmin, Zmax et en position 8 Dz si < 0  2D
 152 //  9 distance max dans le maillage
 153
 154    std::vector<double> LesExtremes;
 155
 156    // Ouverture du Fichier Med
 157    med_idt medIdt = MEDfileOpen(aFile,MED_ACC_RDONLY);
 158    if (medIdt <0)
 159    {
 160           // GERALD -- QMESSAGE BOX
 161           std::cerr << "Error : mesh is unreadable" << std::endl;
 162           return LesExtremes;
 163    }
 164
 165                                 // Le fichier Med est lisible
 166     // Boucle sur les noms de maillage
 167    med_int numberOfMeshes = MEDnMesh(medIdt) ;
 168    if (numberOfMeshes != 1 )
 169    {
 170           // GERALD -- QMESSAGE BOX
 171           std::cerr << "Error : file contains more than one mesh" << std::endl;
 172           return LesExtremes;
 173    }
 174
 175   char meshname[MED_NAME_SIZE+1];
 176   med_int spacedim,meshdim;
 177   med_mesh_type meshtype;
 178   char descriptionription[MED_COMMENT_SIZE+1];
 179   char dtunit[MED_SNAME_SIZE+1];
 180   med_sorting_type sortingtype;
 181   med_int nstep;
 182   med_axis_type axistype;
 183   int naxis = MEDmeshnAxis(medIdt,1);
 184   char *axisname=new char[naxis*MED_SNAME_SIZE+1];
 185   char *axisunit=new char[naxis*MED_SNAME_SIZE+1];
 186   med_err aRet = MEDmeshInfo(medIdt,
 187                           1,
 188                           meshname,
 189                           &spacedim,
 190                           &meshdim,
 191                           &meshtype,
 192                           descriptionription,
 193                           dtunit,
 194                           &sortingtype,
 195                           &nstep,
 196                           &axistype,
 197                           axisname,
 198                           axisunit);
 199
 200    if (aRet < 0)
 201    {
 202           // GERALD -- QMESSAGE BOX
 203           std::cerr << "Error : mesh is unreadable" << std::endl;
 204           return LesExtremes;
 205    }
 206
 207   med_bool chgt,trsf;
 208   med_int nnoe  = MEDmeshnEntity(medIdt,
 209                             meshname,
 210                             MED_NO_DT,
 211                             MED_NO_IT,
 212                             MED_NODE,
 213                             MED_NO_GEOTYPE,
 214                             MED_COORDINATE,
 215                             MED_NO_CMODE,
 216                             &chgt,
 217                             &trsf);
 218    if ( nnoe < 0)
 219    {
 220           // GERALD -- QMESSAGE BOX
 221           std::cerr << "Error : mesh is unreadable" << std::endl;
 222           return LesExtremes;
 223    }
 224
 225   med_float* coo    = (med_float*) malloc(sizeof(med_float)*nnoe*spacedim);
 226
 227   aRet = MEDmeshNodeCoordinateRd(medIdt,
 228                                       meshname,
 229                                       MED_NO_DT,
 230                                       MED_NO_IT,
 231                                       MED_NO_INTERLACE,
 232                                       coo);
 233    if ( aRet < 0)
 234    {
 235           // GERALD -- QMESSAGE BOX
 236           std::cerr << "Error : mesh coordinates are unreadable" << std::endl;
 237           return LesExtremes;
 238    }
 239
 240    med_float xmin,xmax,ymin,ymax,zmin,zmax;
 241
 242    xmin=coo[0];
 243    xmax=coo[0];
 244    for (int i=1;i<nnoe;i++)
 245    {
 246       xmin = std::min(xmin,coo[i]);
 247       xmax = std::max(xmax,coo[i]);
 248    }
 249 //
 250    if (spacedim > 1)
 251    {
 252        ymin=coo[nnoe]; ymax=coo[nnoe];
 253        for (int i=nnoe+1;i<2*nnoe;i++)
 254        {
 255            ymin = std::min(ymin,coo[i]);
 256            ymax = std::max(ymax,coo[i]);
 257        }
 258    }
 259    else
 260    {
 261        ymin=0;
 262        ymax=0;
 263        zmin=0;
 264        zmax=0;
 265    }
 266 //
 267    if (spacedim > 2)
 268    {
 269        zmin=coo[2*nnoe]; zmax=coo[2*nnoe];
 270        for (int i=2*nnoe+1;i<3*nnoe;i++)
 271        {
 272            zmin = std::min(zmin,coo[i]);
 273            zmax = std::max(zmax,coo[i]);
 274        }
 275    }
 276    else
 277    {
 278        zmin=0;
 279        zmax=0;
 280    }
 281    MEDfileClose(medIdt);
 282
 283    MESSAGE( "_______________________________________");
 284    MESSAGE( "xmin : " << xmin << " xmax : " << xmax );
 285    MESSAGE( "ymin : " << ymin << " ymax : " << ymax );
 286    MESSAGE( "zmin : " << zmin << " zmax : " << zmax );
 287    MESSAGE( "_______________________________________" );
 288    double epsilon = 1.e-6 ;
 289    LesExtremes.push_back(xmin);
 290    LesExtremes.push_back(xmax);
 291    LesExtremes.push_back(0);
 292    LesExtremes.push_back(ymin);
 293    LesExtremes.push_back(ymax);
 294    LesExtremes.push_back(0);
 295    LesExtremes.push_back(zmin);
 296    LesExtremes.push_back(zmax);
 297    LesExtremes.push_back(0);
 298
 299
 300    double max1=std::max ( LesExtremes[1] - LesExtremes[0] , LesExtremes[4] - LesExtremes[3] ) ;
 301    double max2=std::max ( max1 , LesExtremes[7] - LesExtremes[6] ) ;
 302    LesExtremes.push_back(max2);
 303
 304 // LesExtremes[0] = Xmini du maillage
 305 // LesExtremes[1] = Xmaxi du maillage
 306 // LesExtremes[2] = increment de progression en X
 307 // LesExtremes[3,4,5] : idem pour Y
 308 // LesExtremes[6,7,8] : idem pour Z
 309 // LesExtremes[9] = ecart maximal entre coordonnees
 310 // On fait un traitement pour dans le cas d'une coordonnee constante
 311 // inhiber ce cas en mettant un increment negatif
 312 //
 313    double diff = LesExtremes[1] - LesExtremes[0];
 314    if (fabs(diff)  > epsilon*max2)
 315    {
 316       LesExtremes[2] = diff/100.;
 317    }
 318    else
 319    {
 320       LesExtremes[2] = -1. ;
 321    }
 322
 323    diff = LesExtremes[4] - LesExtremes[3];
 324    if (fabs(diff)  > epsilon*max2)
 325    {
 326       LesExtremes[5]=diff/100.;
 327    }
 328    else
 329    {
 330       LesExtremes[5] = -1. ;
 331    }
 332
 333    diff = LesExtremes[7] - LesExtremes[6];
 334    if (fabs(diff)  > epsilon*max2)
 335    {
 336       LesExtremes[8]=diff/100.;
 337    }
 338    else
 339    {
 340       LesExtremes[8] = -1. ;
 341    }
 342
 343    MESSAGE ( "_______________________________________" );
 344    MESSAGE ( "xmin : " << LesExtremes[0] << " xmax : " << LesExtremes[1] << " xincr : " << LesExtremes[2] );
 345    MESSAGE ( "ymin : " << LesExtremes[3] << " ymax : " << LesExtremes[4] << " yincr : " << LesExtremes[5] );
 346    MESSAGE ( "zmin : " << LesExtremes[6] << " zmax : " << LesExtremes[7] << " zincr : " << LesExtremes[8] );
 347    MESSAGE ( "dmax : " << LesExtremes[9] );
 348    MESSAGE ( "_______________________________________" );
 349
 350    return  LesExtremes;
 351 }
 352