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