Back to index

salome-med  6.5.0
med_test.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 #include "MEDMEM_Exception.hxx"
00024 #include "MEDMEM_Mesh.hxx"
00025 #include "MEDMEM_Family.hxx"
00026 #include "MEDMEM_Group.hxx"
00027 
00028 #include "MEDMEM_Support.hxx"
00029 #include "MEDMEM_Field.hxx"
00030 #include "MEDMEM_MedMeshDriver.hxx"
00031 #include "MEDMEM_MedFieldDriver.hxx"
00032 #include "MEDMEM_define.hxx"
00033 
00034 #include<string>
00035 
00036 #include <math.h>
00037 #include <stdlib.h>
00038 
00039 using namespace std;
00040 using namespace MEDMEM;
00041 using namespace MED_EN;
00042 
00043 static double dmax(double x, double y) { return (x>y)?x:y;}
00044 
00045 static double dmin(double x, double y) { return (x>y)?y:x;}
00046 
00047 static double infty = 1.e20;
00048 
00049 static void affiche_support(const SUPPORT * mySupport) 
00050 {
00051   cout << "  - Name : "<<mySupport->getName().c_str()<<endl ;
00052   cout << "  - Description : "<<mySupport->getDescription().c_str()<<endl ;
00053   cout << "  - Entity : "<<mySupport->getEntity()<<endl ;
00054   cout << "  - Entities list : "<<endl ;
00055   if (!(mySupport->isOnAllElements())) {
00056     int NumberOfTypes = mySupport->getNumberOfTypes() ;
00057     cout<<"  - NumberOfTypes : "<<NumberOfTypes<<endl;
00058     const medGeometryElement * Types = mySupport->getTypes() ;
00059     for (int j=0;j<NumberOfTypes;j++) {
00060       cout << "    * Type "<<Types[j]<<" : " ;
00061       int NumberOfElements = mySupport->getNumberOfElements(Types[j]) ;
00062       const int * Number = mySupport->getNumber(Types[j]) ;
00063       for (int k=0; k<NumberOfElements;k++)
00064         cout << Number[k] << " ";
00065       cout << endl ;
00066     }
00067   } else
00068     cout << "    Is on all entities !"<< endl;
00069 }
00070 
00071 
00072 static void affiche_famille(MESH *myMesh,medEntityMesh Entity) 
00073 {
00074   int NumberOfFamilies = myMesh->getNumberOfFamilies(Entity) ;
00075   cout << "NumberOfFamilies : "<<NumberOfFamilies<<endl;
00076   for (int i=1; i<NumberOfFamilies+1;i++) {
00077     const FAMILY* myFamily = myMesh->getFamily(Entity,i);
00078     affiche_support(myFamily);
00079     cout << "  - Identifier : "<<myFamily->getIdentifier()<<endl ;
00080     int NumberOfAttributes = myFamily->getNumberOfAttributes() ;
00081     cout << "  - Attributes ("<<NumberOfAttributes<<") :"<<endl;
00082     for (int j=1;j<NumberOfAttributes+1;j++)
00083       cout << "    * "<<myFamily->getAttributeIdentifier(j)<<" : "<<myFamily->getAttributeValue(j)<<", "<<myFamily->getAttributeDescription(j).c_str()<<endl ;
00084     int NumberOfGroups = myFamily->getNumberOfGroups() ;
00085     cout << "  - Groups ("<<NumberOfGroups<<") :"<<endl;
00086     for (int j=1;j<NumberOfGroups+1;j++)
00087       cout << "    * "<<myFamily->getGroupName(j).c_str()<<endl ;
00088   }
00089 }
00090 
00091 static void affiche_groupe(MESH *myMesh,medEntityMesh Entity) 
00092 {
00093   int NumberOfGroups = myMesh->getNumberOfGroups(Entity) ;
00094   cout << "NumberOfGroups : "<<NumberOfGroups<<endl;
00095   for (int i=1; i<NumberOfGroups+1;i++) {
00096     const GROUP* myGroup = myMesh->getGroup(Entity,i);
00097     affiche_support(myGroup);
00098     int NumberOfFamillies = myGroup->getNumberOfFamilies() ;
00099     cout << "  - Families ("<<NumberOfFamillies<<") :"<<endl;
00100     for (int j=1;j<NumberOfFamillies+1;j++)
00101       cout << "    * "<<myGroup->getFamily(j)->getName().c_str()<<endl ;
00102   }
00103 }
00104 
00105 int main (int argc, char ** argv) {
00106 
00107   if ((argc !=3) && (argc != 4)) {
00108     cerr << "Usage : " << argv[0] 
00109          << " filename meshname [fieldname]" << endl << endl;
00110     exit(-1);
00111   }
00112 
00113   string filename = argv[1] ;
00114   string meshname = argv[2] ;
00115 
00116   MESH * myMesh= new MESH(MED_DRIVER,filename,meshname) ;
00117   
00118   int SpaceDimension = myMesh->getSpaceDimension() ;
00119   int MeshDimension  = myMesh->getMeshDimension() ;
00120   int NumberOfNodes  = myMesh->getNumberOfNodes() ;
00121 
00122   cout << "Space Dimension : " << SpaceDimension << endl << endl ; 
00123 
00124   cout << "Mesh Dimension : " << MeshDimension << endl << endl ; 
00125 
00126   const double * Coordinates = myMesh->getCoordinates(MED_FULL_INTERLACE) ;
00127 
00128   cout << "Show Nodes Coordinates : " << endl ;
00129 
00130   cout << "Name :" << endl ;
00131   const string * CoordinatesNames = myMesh->getCoordinatesNames() ;
00132   for(int i=0; i<SpaceDimension ; i++) {
00133     cout << " - " << CoordinatesNames[i] << endl ;
00134   }
00135   cout << "Unit :" << endl ;
00136   const string * CoordinatesUnits = myMesh->getCoordinatesUnits() ;
00137   for(int i=0; i<SpaceDimension ; i++) {
00138     cout << " - " << CoordinatesUnits[i] << endl ;
00139   }
00140   for(int i=0; i<NumberOfNodes ; i++) {
00141     cout << "Nodes " << i+1 << " : " ;
00142     for (int j=0; j<SpaceDimension ; j++)
00143       cout << Coordinates[i*SpaceDimension+j] << " " ;
00144     cout << endl ;
00145   }
00146 
00147   int NumberOfTypes = myMesh->getNumberOfTypes(MED_CELL) ;
00148   const medGeometryElement  * Types = myMesh->getTypes(MED_CELL) ;
00149 
00150   cout << "Show Connectivity (Nodal) :" << endl ;
00151   for (int i=0; i<NumberOfTypes; i++) {
00152     cout << "For type " << Types[i] << " : " << endl ;
00153     int NumberOfElements = myMesh->getNumberOfElements(MED_CELL,Types[i]);
00154     const int * connectivity =  myMesh->getConnectivity(MED_NODAL,MED_CELL,Types[i]);
00155     int NomberOfNodesPerCell = Types[i]%100 ;
00156     for (int j=0;j<NumberOfElements;j++){
00157       cout << "Element "<< j+1 <<" : " ;
00158       for (int k=0;k<NomberOfNodesPerCell;k++)
00159         cout << connectivity[j*NomberOfNodesPerCell+k]<<" ";
00160       cout << endl ;
00161     }
00162   }
00163 
00164   cout << "Show Family :"<<endl ;
00165   affiche_famille(myMesh,MED_NODE);
00166   affiche_famille(myMesh,MED_CELL);
00167   affiche_famille(myMesh,MED_FACE);
00168   affiche_famille(myMesh,MED_EDGE);
00169 
00170   cout << "Show Group :"<<endl ;
00171   affiche_groupe(myMesh,MED_NODE);
00172   affiche_groupe(myMesh,MED_CELL);
00173   affiche_groupe(myMesh,MED_FACE);
00174   affiche_groupe(myMesh,MED_EDGE);
00175 
00176   cout << "Show Reverse Nodal Connectivity :" << endl ;
00177   const int * ReverseNodalConnectivity = myMesh->getReverseConnectivity(MED_NODAL) ;
00178   const int * ReverseNodalConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_NODAL) ;
00179   for (int i=0; i<NumberOfNodes; i++) {
00180     cout << "Node "<<i+1<<" : " ;
00181     for (int j=ReverseNodalConnectivityIndex[i];j<ReverseNodalConnectivityIndex[i+1];j++)
00182       cout << ReverseNodalConnectivity[j-1] << " " ;
00183     cout << endl ;
00184   }
00185 
00186   cout << "Show Connectivity (Descending) :" << endl ;
00187   int NumberOfElements ;
00188   const int * connectivity ;
00189   const int * connectivity_index ;
00190   myMesh->calculateConnectivity(MED_DESCENDING,MED_CELL);
00191   try {
00192     NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
00193     connectivity =  myMesh->getConnectivity(MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS);
00194     connectivity_index =  myMesh->getConnectivityIndex(MED_DESCENDING,MED_CELL);
00195   }
00196   catch (MEDEXCEPTION& m) {
00197     cout << m.what() << endl ;
00198     exit (-1) ;
00199   }
00200   for (int j=0;j<NumberOfElements;j++) {
00201     cout << "Element "<<j+1<<" : " ;
00202     for (int k=connectivity_index[j];k<connectivity_index[j+1];k++)
00203       cout << connectivity[k-1]<<" ";
00204     cout << endl ;
00205   }
00206 
00207   cout << "Show Reverse Descending Connectivity :" << endl ;
00208   const int * ReverseDescendingConnectivity = myMesh->getReverseConnectivity(MED_DESCENDING) ;
00209   const int * ReverseDescendingConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_DESCENDING) ;
00210 
00211   int NumberOfConstituents  = 0;
00212   string constituent ;
00213   medEntityMesh constituentEntity ;
00214 
00215   if (MeshDimension==3) {
00216     constituent = "Face" ;
00217     constituentEntity = MED_FACE ;
00218   }
00219 
00220   if (MeshDimension==2) {
00221     constituent = "Edge" ;
00222     constituentEntity = MED_EDGE ;
00223   }
00224 
00225   if (MeshDimension==1) {
00226     INFOS_MED("ERROR : MeshDimension = 1 !");
00227     INFOS_MED("We could not see Reverse Descending Connectivity.") ;
00228   } else {
00229     NumberOfConstituents = myMesh->getNumberOfElements (constituentEntity,MED_ALL_ELEMENTS);
00230     for (int i=0; i<NumberOfConstituents; i++) {
00231       cout << constituent <<i+1<<" : " ;
00232       for (int j=ReverseDescendingConnectivityIndex[i];j<ReverseDescendingConnectivityIndex[i+1];j++)
00233         cout << ReverseDescendingConnectivity[j-1] << " " ;
00234       cout << endl ;
00235     }
00236   }
00237   cout << "Show "<<constituent<<" Connectivity (Nodal) :" << endl ;
00238   const int * face_connectivity =  myMesh->getConnectivity(MED_NODAL,constituentEntity,MED_ALL_ELEMENTS);
00239   const int * face_connectivity_index =  myMesh->getConnectivityIndex(MED_NODAL,constituentEntity);
00240   for (int i=0; i<NumberOfConstituents; i++) {
00241     cout << constituent <<i+1<<" : " ;
00242     for (int j=face_connectivity_index[i];j<face_connectivity_index[i+1];j++)
00243       cout << face_connectivity[j-1]<<" ";
00244     cout << endl ;
00245   }
00246 
00247   /* test of normal, area, volume, barycenter */
00248 
00249   const SUPPORT* support1 = myMesh->getSupportOnAll(constituentEntity);
00250   cout << "Building of the Support on all cells dimensionned (Meshdim-1) of the mesh :"<< endl ;
00251   cout << "Face in 3D or Edge in 2D" << endl;
00252 
00253   cout << "Getting the normal of each face of this support !" << endl ;
00254 
00255   FIELD<double>* normal = myMesh->getNormal(support1);
00256 
00257   double normal_square, norm ;
00258   double maxnorm=-infty;
00259   double minnorm=infty;
00260   double tmp_value ;
00261   for (int i = 1; i<=NumberOfConstituents;i++)
00262     {
00263       normal_square = 0. ;
00264       cout << "Normal " << i << " " ; 
00265       for (int j=1; j<=SpaceDimension; j++)
00266         {
00267           tmp_value = normal->getValueIJ(i,j) ;
00268           normal_square += tmp_value*tmp_value ;
00269           cout << tmp_value << " " ;
00270         }
00271       norm = sqrt(normal_square);
00272       maxnorm = dmax(maxnorm,norm);
00273       minnorm = dmin(minnorm,norm);
00274       cout << ", Norm = " << norm << endl;
00275     }
00276   cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
00277 
00278   if(normal)
00279     normal->removeReference() ;
00280 
00281   if (SpaceDimension == 2)
00282     {
00283       cout << "Getting the length of each edge !" << endl ;
00284 
00285       FIELD<double>* length = myMesh->getLength(support1);
00286 
00287       double length_value,maxlength,minlength;
00288       maxlength = -infty;
00289       minlength = infty;
00290       for (int i = 1; i<=NumberOfConstituents;i++)
00291         {
00292           length_value = length->getValueIJ(i,1) ;
00293           cout << "Length " << i << " " << length_value << endl;
00294           maxlength = dmax(maxlength,length_value);
00295           minlength = dmin(minlength,length_value);
00296         }
00297       cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
00298       if(length)
00299         length->removeReference();
00300     }
00301 
00302   cout << "Building of the Support on all space-dimensionned cells of the mesh :"<< endl ;
00303   const SUPPORT * support = myMesh->getSupportOnAll( MED_CELL );
00304 
00305   cout << "Getting the barycenter of each element of this support !" << endl ;
00306 
00307   FIELD<double>* barycenter = myMesh->getBarycenter(support);
00308   NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
00309 
00310   for (int i = 1; i<=NumberOfElements;i++)
00311     {
00312       if (SpaceDimension == 3)
00313         cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
00314 
00315       if (SpaceDimension == 2)
00316         cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << endl;
00317     }
00318   if(barycenter)
00319     barycenter->removeReference();
00320 
00321   if (SpaceDimension == 3)
00322     {
00323       cout << "Getting the Volume of each element of this support which is a 3D one !" << endl;
00324 
00325       FIELD<double>* volume = myMesh->getVolume(support);
00326 
00327       double maxvol,minvol,voltot;
00328       maxvol = -infty;
00329       minvol = infty;
00330       voltot = 0.0;
00331       for (int i = 1; i<=NumberOfElements;i++)
00332         {
00333           cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
00334           maxvol = dmax(maxvol,volume->getValueIJ(i,1));
00335           minvol = dmin(minvol,volume->getValueIJ(i,1));
00336           voltot = voltot + volume->getValueIJ(i,1);
00337         }
00338 
00339       cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
00340       cout << "Support Volume " << voltot << endl;
00341       if(volume)
00342         volume->removeReference() ;
00343     }
00344   else if (SpaceDimension == 2)
00345     {
00346       cout << "Getting the Area of each element of this support which is a 2D one !" << endl;
00347 
00348       FIELD<double>* area = myMesh->getArea(support);
00349 
00350       double maxarea,minarea,areatot;
00351       maxarea = -infty;
00352       minarea = infty;
00353       areatot = 0.0;
00354       for (int i = 1; i<=NumberOfElements;i++)
00355         {
00356           cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
00357           maxarea = dmax(maxarea,area->getValueIJ(i,1));
00358           minarea = dmin(minarea,area->getValueIJ(i,1));
00359           areatot = areatot + area->getValueIJ(i,1);
00360         }
00361 
00362       cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
00363       cout << "Support Area " << areatot << endl;
00364       if(area)
00365         area->removeReference();
00366     }
00367 
00368   if (argc < 4) return 0;
00369 
00370   // read field :
00371 
00372   if (argc != 4) exit(0) ;
00373   // else we have a field !
00374 
00375   string fieldname = argv[3];
00376 
00377   const SUPPORT * mySupport = myMesh->getSupportOnAll(MED_CELL);
00378   FIELD<double> * myField= new FIELD<double>() ;
00379 
00380   myField->setName(fieldname);
00381   myField->setSupport(mySupport);
00382   MED_FIELD_RDONLY_DRIVER<double> myFieldDriver(filename,myField) ;
00383   myFieldDriver.setFieldName(fieldname);
00384   myFieldDriver.open() ;
00385 
00386   try
00387     {
00388       myFieldDriver.read() ;
00389     }
00390   catch (...)
00391     {
00392       mySupport = myMesh->getSupportOnAll(MED_NODE);
00393       myField->setSupport(mySupport);
00394       try
00395         {
00396           myFieldDriver.read() ;
00397         }
00398       catch (...)
00399         {
00400           cout << "Field " << fieldname << " not found !!!" << endl ;
00401           exit (-1) ;
00402         }
00403     }
00404 
00405   myFieldDriver.close() ;
00406 
00407   cout << "Field "<< myField->getName() << " : " <<myField->getDescription() <<  endl ;
00408   int NumberOfComponents = myField->getNumberOfComponents() ;
00409   cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
00410   for (int i=1; i<NumberOfComponents+1; i++)
00411     {
00412       cout << "  - composante "<<i<<" :"<<endl ;
00413       cout << "      - nom         : "<<myField->getComponentName(i)<< endl;
00414       cout << "      - description : "<<myField->getComponentDescription(i) << endl;
00415       cout << "      - units       : "<<myField->getMEDComponentUnit(i) << endl;
00416     }
00417   cout << "- iteration :" << endl ;
00418   cout << "    - numero : " << myField->getIterationNumber()<< endl  ;
00419   cout << "    - ordre  : " << myField->getOrderNumber()<< endl  ;
00420   cout << "    - temps  : " << myField->getTime()<< endl  ;
00421 
00422   cout << "- Valeurs :"<<endl;
00423   int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
00424   MEDMEM_Array<double> * myvalue = myField->getArrayNoGauss();
00425   const double * value ;
00426   for (int i=1; i<NumberOf+1; i++)
00427     {
00428       value = myvalue->getRow(i) ;
00429       for (int j=0; j<NumberOfComponents; j++)
00430         cout << value[j]<< " ";
00431       cout<<endl;
00432     }
00433   cout<<endl;
00434 
00435   myField->removeReference();
00436   myMesh->removeReference();
00437 
00438   return 0;
00439 }