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(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     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       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     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     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   string * CoordinatesNames = myMesh->getCoordinatesNames() ;
 121   for(int i=0; i<SpaceDimension ; i++) {
 122     cout << " - " << CoordinatesNames[i] << endl ;
 123   }
 124   cout << "Unit :" << endl ;
 125   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   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     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   int * ReverseNodalConnectivity = myMesh->getReverseConnectivity(MED_NODAL) ;
 167   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   int * connectivity ;
 178   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   int * ReverseDescendingConnectivity = myMesh->getReverseConnectivity(MED_DESCENDING) ;
 198   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     MESSAGE("ERROR : MeshDimension = 1 !");
 216     MESSAGE("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   int * face_connectivity =  myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,constituentEntity,MED_ALL_ELEMENTS);
 228   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   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
 274   if (SpaceDimension == 2)
 275     {
 276       cout << "Getting the length of each edge !" << endl ;
 277
 278       length = myMesh->getLength(support1);
 279
 280       double length_value,maxlength,minlength;
 281       maxlength = -infty;
 282       minlength = infty;
 283       for (int i = 1; i<=NumberOfConstituents;i++)
 284         {
 285           length_value = length->getValueIJ(i,1) ;
 286           cout << "Length " << i << " " << length_value << endl;
 287           maxlength = dmax(maxlength,length_value);
 288           minlength = dmin(minlength,length_value);
 289         }
 290       cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
 291     }
 292
 293   cout << "Building of the Support on all space-dimensionned cells of the mesh :"<< endl ;
 294   SUPPORT * support = new SUPPORT(myMesh);
 295
 296   cout << "Getting the barycenter of each element of this support !" << endl ;
 297
 298   FIELD<double>* barycenter = new FIELD<double>::FIELD();
 299
 300   barycenter = myMesh->getBarycenter(support);
 301   NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
 302
 303   for (int i = 1; i<=NumberOfElements;i++)
 304     {
 305       if (SpaceDimension == 3)
 306         cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
 307
 308       if (SpaceDimension == 2)
 309         cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << endl;
 310     }
 311
 312   FIELD<double>* volume = new FIELD<double>::FIELD();
 313   FIELD<double>* area = new FIELD<double>::FIELD();
 314   volume = NULL;
 315   area = NULL;
 316
 317   if (SpaceDimension == 3)
 318     {
 319       cout << "Getting the Volume of each element of this support which is a 3D one !" << endl;
 320
 321       volume = myMesh->getVolume(support);
 322
 323       double maxvol,minvol,voltot;
 324       maxvol = -infty;
 325       minvol = infty;
 326       voltot = 0.0;
 327       for (int i = 1; i<=NumberOfElements;i++)
 328         {
 329           cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
 330           maxvol = dmax(maxvol,volume->getValueIJ(i,1));
 331           minvol = dmin(minvol,volume->getValueIJ(i,1));
 332           voltot = voltot + volume->getValueIJ(i,1);
 333         }
 334
 335       cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
 336       cout << "Support Volume " << voltot << endl;
 337     }
 338   else if (SpaceDimension == 2)
 339     {
 340       cout << "Getting the Area of each element of this support which is a 2D one !" << endl;
 341
 342       area = myMesh->getArea(support);
 343
 344       //    cout << "nb of comp "<< area->getNumberOfComponents() << " length " << area->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS) << endl;
 345
 346       double maxarea,minarea,areatot;
 347       maxarea = -infty;
 348       minarea = infty;
 349       areatot = 0.0;
 350       for (int i = 1; i<=NumberOfElements;i++)
 351         {
 352           cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
 353           maxarea = dmax(maxarea,area->getValueIJ(i,1));
 354           minarea = dmin(minarea,area->getValueIJ(i,1));
 355           areatot = areatot + area->getValueIJ(i,1);
 356         }
 357
 358       cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
 359       cout << "Support Area " << areatot << endl;
 360     }
 361
 362   if (barycenter != NULL) delete barycenter;
 363   if (volume != NULL ) delete volume;
 364   if (area != NULL ) delete area;
 365
 366
 367   if (argc < 4) return 0;
 368
 369   // read field :
 370
 371   if (argc != 4) exit(0) ;
 372   // else we have a field !
 373
 374   string fieldname = argv[3];
 375
 376   //  SUPPORT * mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
 377   SUPPORT * mySupport = new SUPPORT(myMesh,"On_all_cell",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     myFieldDriver.read() ;
 388   } catch (...) {
 389     delete mySupport ;
 390     mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
 391     myField->setSupport(mySupport);
 392     try {
 393       myFieldDriver.read() ;
 394     } catch (...) {
 395       cout << "Field " << fieldname << " not found !!!" << endl ;
 396       exit (-1) ;
 397     }
 398   }
 399
 400   myFieldDriver.close() ;
 401
 402   cout << "Field "<< myField->getName() << " : " <<myField->getDescription() <<  endl ;
 403   int NumberOfComponents = myField->getNumberOfComponents() ;
 404   cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
 405   for (int i=1; i<NumberOfComponents+1; i++) {
 406     cout << "  - composante "<<i<<" :"<<endl ;
 407     cout << "      - nom         : "<<myField->getComponentName(i)<< endl;
 408     cout << "      - description : "<<myField->getComponentDescription(i) << endl;
 409     cout << "      - units       : "<<myField->getMEDComponentUnit(i) << endl;
 410   }
 411   cout << "- iteration :" << endl ;
 412   cout << "    - numero : " << myField->getIterationNumber()<< endl  ;
 413   cout << "    - ordre  : " << myField->getOrderNumber()<< endl  ;
 414   cout << "    - temps  : " << myField->getTime()<< endl  ;
 415
 416   cout << "- Valeurs :"<<endl;
 417   int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
 418   //    for (int i=1; i<NumberOfComponents+1; i++) {
 419   //      double * value = myField->getValueI(MED_NO_INTERLACE,i) ;
 420   //      for (int j=0; j<NumberOf; j++)
 421   //        cout << value[j]<< " ";
 422   //      cout<<endl;
 423   //    }
 424   for (int i=1; i<NumberOf+1; i++) {
 425     double * value = myField->getValueI(MED_FULL_INTERLACE,i) ;
 426     for (int j=0; j<NumberOfComponents; j++)
 427       cout << value[j]<< " ";
 428     cout<<endl;
 429   }
 430
 431
 432   return 0;
 433 }