src/MEDMEM/med_test.cxx

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