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