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