src/MEDMEMBinTest/med_test.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "MEDMEM_Exception.hxx"
  24 #include "MEDMEM_Mesh.hxx"
  25 #include "MEDMEM_Family.hxx"
  26 #include "MEDMEM_Group.hxx"
  27
  28 #include "MEDMEM_Support.hxx"
  29 #include "MEDMEM_Field.hxx"
  30 #include "MEDMEM_MedMeshDriver.hxx"
  31 #include "MEDMEM_MedFieldDriver.hxx"
  32 #include "MEDMEM_define.hxx"
  33
  34 #include<string>
  35
  36 #include <math.h>
  37 #include <stdlib.h>
  38
  39 using namespace std;
  40 using namespace MEDMEM;
  41 using namespace MED_EN;
  42
  43 static double dmax(double x, double y) { return (x>y)?x:y;}
  44
  45 static double dmin(double x, double y) { return (x>y)?y:x;}
  46
  47 static double infty = 1.e20;
  48
  49 static void affiche_support(const SUPPORT * mySupport)
  50 {
  51   cout << "  - Name : "<<mySupport->getName().c_str()<<endl ;
  52   cout << "  - Description : "<<mySupport->getDescription().c_str()<<endl ;
  53   cout << "  - Entity : "<<mySupport->getEntity()<<endl ;
  54   cout << "  - Entities list : "<<endl ;
  55   if (!(mySupport->isOnAllElements())) {
  56     int NumberOfTypes = mySupport->getNumberOfTypes() ;
  57     cout<<"  - NumberOfTypes : "<<NumberOfTypes<<endl;
  58     const medGeometryElement * Types = mySupport->getTypes() ;
  59     for (int j=0;j<NumberOfTypes;j++) {
  60       cout << "    * Type "<<Types[j]<<" : " ;
  61       int NumberOfElements = mySupport->getNumberOfElements(Types[j]) ;
  62       const int * Number = mySupport->getNumber(Types[j]) ;
  63       for (int k=0; k<NumberOfElements;k++)
  64         cout << Number[k] << " ";
  65       cout << endl ;
  66     }
  67   } else
  68     cout << "    Is on all entities !"<< endl;
  69 }
  70
  71
  72 static void affiche_famille(MESH *myMesh,medEntityMesh Entity)
  73 {
  74   int NumberOfFamilies = myMesh->getNumberOfFamilies(Entity) ;
  75   cout << "NumberOfFamilies : "<<NumberOfFamilies<<endl;
  76   for (int i=1; i<NumberOfFamilies+1;i++) {
  77     const FAMILY* myFamily = myMesh->getFamily(Entity,i);
  78     affiche_support(myFamily);
  79     cout << "  - Identifier : "<<myFamily->getIdentifier()<<endl ;
  80     int NumberOfAttributes = myFamily->getNumberOfAttributes() ;
  81     cout << "  - Attributes ("<<NumberOfAttributes<<") :"<<endl;
  82     for (int j=1;j<NumberOfAttributes+1;j++)
  83       cout << "    * "<<myFamily->getAttributeIdentifier(j)<<" : "<<myFamily->getAttributeValue(j)<<", "<<myFamily->getAttributeDescription(j).c_str()<<endl ;
  84     int NumberOfGroups = myFamily->getNumberOfGroups() ;
  85     cout << "  - Groups ("<<NumberOfGroups<<") :"<<endl;
  86     for (int j=1;j<NumberOfGroups+1;j++)
  87       cout << "    * "<<myFamily->getGroupName(j).c_str()<<endl ;
  88   }
  89 }
  90
  91 static void affiche_groupe(MESH *myMesh,medEntityMesh Entity)
  92 {
  93   int NumberOfGroups = myMesh->getNumberOfGroups(Entity) ;
  94   cout << "NumberOfGroups : "<<NumberOfGroups<<endl;
  95   for (int i=1; i<NumberOfGroups+1;i++) {
  96     const GROUP* myGroup = myMesh->getGroup(Entity,i);
  97     affiche_support(myGroup);
  98     int NumberOfFamillies = myGroup->getNumberOfFamilies() ;
  99     cout << "  - Families ("<<NumberOfFamillies<<") :"<<endl;
 100     for (int j=1;j<NumberOfFamillies+1;j++)
 101       cout << "    * "<<myGroup->getFamily(j)->getName().c_str()<<endl ;
 102   }
 103 }
 104
 105 int main (int argc, char ** argv) {
 106
 107   if ((argc !=3) && (argc != 4)) {
 108     cerr << "Usage : " << argv[0]
 109          << " filename meshname [fieldname]" << endl << endl;
 110     exit(-1);
 111   }
 112
 113   string filename = argv[1] ;
 114   string meshname = argv[2] ;
 115
 116   MESH * myMesh= new MESH(MED_DRIVER,filename,meshname) ;
 117
 118   int SpaceDimension = myMesh->getSpaceDimension() ;
 119   int MeshDimension  = myMesh->getMeshDimension() ;
 120   int NumberOfNodes  = myMesh->getNumberOfNodes() ;
 121
 122   cout << "Space Dimension : " << SpaceDimension << endl << endl ;
 123
 124   cout << "Mesh Dimension : " << MeshDimension << endl << endl ;
 125
 126   const double * Coordinates = myMesh->getCoordinates(MED_FULL_INTERLACE) ;
 127
 128   cout << "Show Nodes Coordinates : " << endl ;
 129
 130   cout << "Name :" << endl ;
 131   const string * CoordinatesNames = myMesh->getCoordinatesNames() ;
 132   for(int i=0; i<SpaceDimension ; i++) {
 133     cout << " - " << CoordinatesNames[i] << endl ;
 134   }
 135   cout << "Unit :" << endl ;
 136   const string * CoordinatesUnits = myMesh->getCoordinatesUnits() ;
 137   for(int i=0; i<SpaceDimension ; i++) {
 138     cout << " - " << CoordinatesUnits[i] << endl ;
 139   }
 140   for(int i=0; i<NumberOfNodes ; i++) {
 141     cout << "Nodes " << i+1 << " : " ;
 142     for (int j=0; j<SpaceDimension ; j++)
 143       cout << Coordinates[i*SpaceDimension+j] << " " ;
 144     cout << endl ;
 145   }
 146
 147   int NumberOfTypes = myMesh->getNumberOfTypes(MED_CELL) ;
 148   const medGeometryElement  * Types = myMesh->getTypes(MED_CELL) ;
 149
 150   cout << "Show Connectivity (Nodal) :" << endl ;
 151   for (int i=0; i<NumberOfTypes; i++) {
 152     cout << "For type " << Types[i] << " : " << endl ;
 153     int NumberOfElements = myMesh->getNumberOfElements(MED_CELL,Types[i]);
 154     const int * connectivity =  myMesh->getConnectivity(MED_NODAL,MED_CELL,Types[i]);
 155     int NomberOfNodesPerCell = Types[i]%100 ;
 156     for (int j=0;j<NumberOfElements;j++){
 157       cout << "Element "<< j+1 <<" : " ;
 158       for (int k=0;k<NomberOfNodesPerCell;k++)
 159         cout << connectivity[j*NomberOfNodesPerCell+k]<<" ";
 160       cout << endl ;
 161     }
 162   }
 163
 164   cout << "Show Family :"<<endl ;
 165   affiche_famille(myMesh,MED_NODE);
 166   affiche_famille(myMesh,MED_CELL);
 167   affiche_famille(myMesh,MED_FACE);
 168   affiche_famille(myMesh,MED_EDGE);
 169
 170   cout << "Show Group :"<<endl ;
 171   affiche_groupe(myMesh,MED_NODE);
 172   affiche_groupe(myMesh,MED_CELL);
 173   affiche_groupe(myMesh,MED_FACE);
 174   affiche_groupe(myMesh,MED_EDGE);
 175
 176   cout << "Show Reverse Nodal Connectivity :" << endl ;
 177   const int * ReverseNodalConnectivity = myMesh->getReverseConnectivity(MED_NODAL) ;
 178   const int * ReverseNodalConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_NODAL) ;
 179   for (int i=0; i<NumberOfNodes; i++) {
 180     cout << "Node "<<i+1<<" : " ;
 181     for (int j=ReverseNodalConnectivityIndex[i];j<ReverseNodalConnectivityIndex[i+1];j++)
 182       cout << ReverseNodalConnectivity[j-1] << " " ;
 183     cout << endl ;
 184   }
 185
 186   cout << "Show Connectivity (Descending) :" << endl ;
 187   int NumberOfElements ;
 188   const int * connectivity ;
 189   const int * connectivity_index ;
 190   myMesh->calculateConnectivity(MED_DESCENDING,MED_CELL);
 191   try {
 192     NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
 193     connectivity =  myMesh->getConnectivity(MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS);
 194     connectivity_index =  myMesh->getConnectivityIndex(MED_DESCENDING,MED_CELL);
 195   }
 196   catch (MEDEXCEPTION& m) {
 197     cout << m.what() << endl ;
 198     exit (-1) ;
 199   }
 200   for (int j=0;j<NumberOfElements;j++) {
 201     cout << "Element "<<j+1<<" : " ;
 202     for (int k=connectivity_index[j];k<connectivity_index[j+1];k++)
 203       cout << connectivity[k-1]<<" ";
 204     cout << endl ;
 205   }
 206
 207   cout << "Show Reverse Descending Connectivity :" << endl ;
 208   const int * ReverseDescendingConnectivity = myMesh->getReverseConnectivity(MED_DESCENDING) ;
 209   const int * ReverseDescendingConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_DESCENDING) ;
 210
 211   int NumberOfConstituents  = 0;
 212   string constituent ;
 213   medEntityMesh constituentEntity ;
 214
 215   if (MeshDimension==3) {
 216     constituent = "Face" ;
 217     constituentEntity = MED_FACE ;
 218   }
 219
 220   if (MeshDimension==2) {
 221     constituent = "Edge" ;
 222     constituentEntity = MED_EDGE ;
 223   }
 224
 225   if (MeshDimension==1) {
 226     INFOS_MED("ERROR : MeshDimension = 1 !");
 227     INFOS_MED("We could not see Reverse Descending Connectivity.") ;
 228   } else {
 229     NumberOfConstituents = myMesh->getNumberOfElements (constituentEntity,MED_ALL_ELEMENTS);
 230     for (int i=0; i<NumberOfConstituents; i++) {
 231       cout << constituent <<i+1<<" : " ;
 232       for (int j=ReverseDescendingConnectivityIndex[i];j<ReverseDescendingConnectivityIndex[i+1];j++)
 233         cout << ReverseDescendingConnectivity[j-1] << " " ;
 234       cout << endl ;
 235     }
 236   }
 237   cout << "Show "<<constituent<<" Connectivity (Nodal) :" << endl ;
 238   const int * face_connectivity =  myMesh->getConnectivity(MED_NODAL,constituentEntity,MED_ALL_ELEMENTS);
 239   const int * face_connectivity_index =  myMesh->getConnectivityIndex(MED_NODAL,constituentEntity);
 240   for (int i=0; i<NumberOfConstituents; i++) {
 241     cout << constituent <<i+1<<" : " ;
 242     for (int j=face_connectivity_index[i];j<face_connectivity_index[i+1];j++)
 243       cout << face_connectivity[j-1]<<" ";
 244     cout << endl ;
 245   }
 246
 247   /* test of normal, area, volume, barycenter */
 248
 249   const SUPPORT* support1 = myMesh->getSupportOnAll(constituentEntity);
 250   cout << "Building of the Support on all cells dimensionned (Meshdim-1) of the mesh :"<< endl ;
 251   cout << "Face in 3D or Edge in 2D" << endl;
 252
 253   cout << "Getting the normal of each face of this support !" << endl ;
 254
 255   FIELD<double>* normal = myMesh->getNormal(support1);
 256
 257   double normal_square, norm ;
 258   double maxnorm=-infty;
 259   double minnorm=infty;
 260   double tmp_value ;
 261   for (int i = 1; i<=NumberOfConstituents;i++)
 262     {
 263       normal_square = 0. ;
 264       cout << "Normal " << i << " " ;
 265       for (int j=1; j<=SpaceDimension; j++)
 266         {
 267           tmp_value = normal->getValueIJ(i,j) ;
 268           normal_square += tmp_value*tmp_value ;
 269           cout << tmp_value << " " ;
 270         }
 271       norm = sqrt(normal_square);
 272       maxnorm = dmax(maxnorm,norm);
 273       minnorm = dmin(minnorm,norm);
 274       cout << ", Norm = " << norm << endl;
 275     }
 276   cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
 277
 278   if(normal)
 279     normal->removeReference() ;
 280
 281   if (SpaceDimension == 2)
 282     {
 283       cout << "Getting the length of each edge !" << endl ;
 284
 285       FIELD<double>* length = myMesh->getLength(support1);
 286
 287       double length_value,maxlength,minlength;
 288       maxlength = -infty;
 289       minlength = infty;
 290       for (int i = 1; i<=NumberOfConstituents;i++)
 291         {
 292           length_value = length->getValueIJ(i,1) ;
 293           cout << "Length " << i << " " << length_value << endl;
 294           maxlength = dmax(maxlength,length_value);
 295           minlength = dmin(minlength,length_value);
 296         }
 297       cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
 298       if(length)
 299         length->removeReference();
 300     }
 301
 302   cout << "Building of the Support on all space-dimensionned cells of the mesh :"<< endl ;
 303   const SUPPORT * support = myMesh->getSupportOnAll( MED_CELL );
 304
 305   cout << "Getting the barycenter of each element of this support !" << endl ;
 306
 307   FIELD<double>* barycenter = myMesh->getBarycenter(support);
 308   NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
 309
 310   for (int i = 1; i<=NumberOfElements;i++)
 311     {
 312       if (SpaceDimension == 3)
 313         cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
 314
 315       if (SpaceDimension == 2)
 316         cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << endl;
 317     }
 318   if(barycenter)
 319     barycenter->removeReference();
 320
 321   if (SpaceDimension == 3)
 322     {
 323       cout << "Getting the Volume of each element of this support which is a 3D one !" << endl;
 324
 325       FIELD<double>* volume = myMesh->getVolume(support);
 326
 327       double maxvol,minvol,voltot;
 328       maxvol = -infty;
 329       minvol = infty;
 330       voltot = 0.0;
 331       for (int i = 1; i<=NumberOfElements;i++)
 332         {
 333           cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
 334           maxvol = dmax(maxvol,volume->getValueIJ(i,1));
 335           minvol = dmin(minvol,volume->getValueIJ(i,1));
 336           voltot = voltot + volume->getValueIJ(i,1);
 337         }
 338
 339       cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
 340       cout << "Support Volume " << voltot << endl;
 341       if(volume)
 342         volume->removeReference() ;
 343     }
 344   else if (SpaceDimension == 2)
 345     {
 346       cout << "Getting the Area of each element of this support which is a 2D one !" << endl;
 347
 348       FIELD<double>* area = myMesh->getArea(support);
 349
 350       double maxarea,minarea,areatot;
 351       maxarea = -infty;
 352       minarea = infty;
 353       areatot = 0.0;
 354       for (int i = 1; i<=NumberOfElements;i++)
 355         {
 356           cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
 357           maxarea = dmax(maxarea,area->getValueIJ(i,1));
 358           minarea = dmin(minarea,area->getValueIJ(i,1));
 359           areatot = areatot + area->getValueIJ(i,1);
 360         }
 361
 362       cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
 363       cout << "Support Area " << areatot << endl;
 364       if(area)
 365         area->removeReference();
 366     }
 367
 368   if (argc < 4) return 0;
 369
 370   // read field :
 371
 372   if (argc != 4) exit(0) ;
 373   // else we have a field !
 374
 375   string fieldname = argv[3];
 376
 377   const SUPPORT * mySupport = myMesh->getSupportOnAll(MED_CELL);
 378   FIELD<double> * myField= new FIELD<double>() ;
 379
 380   myField->setName(fieldname);
 381   myField->setSupport(mySupport);
 382   MED_FIELD_RDONLY_DRIVER<double> myFieldDriver(filename,myField) ;
 383   myFieldDriver.setFieldName(fieldname);
 384   myFieldDriver.open() ;
 385
 386   try
 387     {
 388       myFieldDriver.read() ;
 389     }
 390   catch (...)
 391     {
 392       mySupport = myMesh->getSupportOnAll(MED_NODE);
 393       myField->setSupport(mySupport);
 394       try
 395         {
 396           myFieldDriver.read() ;
 397         }
 398       catch (...)
 399         {
 400           cout << "Field " << fieldname << " not found !!!" << endl ;
 401           exit (-1) ;
 402         }
 403     }
 404
 405   myFieldDriver.close() ;
 406
 407   cout << "Field "<< myField->getName() << " : " <<myField->getDescription() <<  endl ;
 408   int NumberOfComponents = myField->getNumberOfComponents() ;
 409   cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
 410   for (int i=1; i<NumberOfComponents+1; i++)
 411     {
 412       cout << "  - composante "<<i<<" :"<<endl ;
 413       cout << "      - nom         : "<<myField->getComponentName(i)<< endl;
 414       cout << "      - description : "<<myField->getComponentDescription(i) << endl;
 415       cout << "      - units       : "<<myField->getMEDComponentUnit(i) << endl;
 416     }
 417   cout << "- iteration :" << endl ;
 418   cout << "    - numero : " << myField->getIterationNumber()<< endl  ;
 419   cout << "    - ordre  : " << myField->getOrderNumber()<< endl  ;
 420   cout << "    - temps  : " << myField->getTime()<< endl  ;
 421
 422   cout << "- Valeurs :"<<endl;
 423   int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
 424   MEDMEM_Array<double> * myvalue = myField->getArrayNoGauss();
 425   const double * value ;
 426   for (int i=1; i<NumberOf+1; i++)
 427     {
 428       value = myvalue->getRow(i) ;
 429       for (int j=0; j<NumberOfComponents; j++)
 430         cout << value[j]<< " ";
 431       cout<<endl;
 432     }
 433   cout<<endl;
 434
 435   myField->removeReference();
 436   myMesh->removeReference();
 437
 438   return 0;
 439 }