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