src/MEDMEM/med_test.cxx

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