src/MEDMEM/med_test.cxx

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