7 #include "MEDMEM_Exception.hxx"
8 #include "MEDMEM_Mesh.hxx"
9 #include "MEDMEM_Family.hxx"
10 #include "MEDMEM_Group.hxx"
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"
19 double dmax(double x, double y) { return (x>y)?x:y;}
21 double dmin(double x, double y) { return (x>y)?y:x;}
25 void affiche_support(const SUPPORT * mySupport)
27 cout << " - Name : "<<mySupport->getName().c_str()<<endl ;
28 cout << " - Description : "<<mySupport->getDescription().c_str()<<endl ;
29 cout << " - Entity : "<<mySupport->getEntity()<<endl ;
30 cout << " - Entities list : "<<endl ;
31 if (!(mySupport->isOnAllElements())) {
32 int NumberOfTypes = mySupport->getNumberOfTypes() ;
33 cout<<" - NumberOfTypes : "<<NumberOfTypes<<endl;
34 const medGeometryElement * Types = mySupport->getTypes() ;
35 for (int j=0;j<NumberOfTypes;j++) {
36 cout << " * Type "<<Types[j]<<" : " ;
37 int NumberOfElements = mySupport->getNumberOfElements(Types[j]) ;
38 const int * Number = mySupport->getNumber(Types[j]) ;
39 for (int k=0; k<NumberOfElements;k++)
40 cout << Number[k] << " ";
44 cout << " Is on all entities !"<< endl;
48 void affiche_famille(MESH *myMesh,medEntityMesh Entity)
50 int NumberOfFamilies = myMesh->getNumberOfFamilies(Entity) ;
51 cout << "NumberOfFamilies : "<<NumberOfFamilies<<endl;
52 for (int i=1; i<NumberOfFamilies+1;i++) {
53 const FAMILY* myFamily = myMesh->getFamily(Entity,i);
54 affiche_support(myFamily);
55 cout << " - Identifier : "<<myFamily->getIdentifier()<<endl ;
56 int NumberOfAttributes = myFamily->getNumberOfAttributes() ;
57 cout << " - Attributes ("<<NumberOfAttributes<<") :"<<endl;
58 for (int j=1;j<NumberOfAttributes+1;j++)
59 cout << " * "<<myFamily->getAttributeIdentifier(j)<<" : "<<myFamily->getAttributeValue(j)<<", "<<myFamily->getAttributeDescription(j).c_str()<<endl ;
60 int NumberOfGroups = myFamily->getNumberOfGroups() ;
61 cout << " - Groups ("<<NumberOfGroups<<") :"<<endl;
62 for (int j=1;j<NumberOfGroups+1;j++)
63 cout << " * "<<myFamily->getGroupName(j).c_str()<<endl ;
67 void affiche_groupe(MESH *myMesh,medEntityMesh Entity)
69 int NumberOfGroups = myMesh->getNumberOfGroups(Entity) ;
70 cout << "NumberOfGroups : "<<NumberOfGroups<<endl;
71 for (int i=1; i<NumberOfGroups+1;i++) {
72 const GROUP* myGroup = myMesh->getGroup(Entity,i);
73 affiche_support(myGroup);
74 int NumberOfFamillies = myGroup->getNumberOfFamilies() ;
75 cout << " - Families ("<<NumberOfFamillies<<") :"<<endl;
76 for (int j=1;j<NumberOfFamillies+1;j++)
77 cout << " * "<<myGroup->getFamily(j)->getName().c_str()<<endl ;
81 int main (int argc, char ** argv) {
85 if ((argc !=3) && (argc != 4)) {
86 cerr << "Usage : " << argv[0]
87 << " filename meshname [fieldname]" << endl << endl;
91 string filename = argv[1] ;
92 string meshname = argv[2] ;
94 // MESH * myMesh= new MESH(MED_DRIVER,filename,meshname) ;
95 MESH * myMesh= new MESH() ;
96 myMesh->setName(meshname);
97 MED_MESH_RDONLY_DRIVER myMeshDriver(filename,myMesh) ;
98 myMeshDriver.setMeshName(meshname);
100 myMeshDriver.read() ;
101 myMeshDriver.close() ;
103 // int drv = myMesh->addDriver(MED_DRIVER,"sortie.med",meshname);
104 // myMesh->write(drv);
107 int SpaceDimension = myMesh->getSpaceDimension() ;
108 int MeshDimension = myMesh->getMeshDimension() ;
109 int NumberOfNodes = myMesh->getNumberOfNodes() ;
111 cout << "Space Dimension : " << SpaceDimension << endl << endl ;
113 cout << "Mesh Dimension : " << MeshDimension << endl << endl ;
115 const double * Coordinates = myMesh->getCoordinates(MED_FULL_INTERLACE) ;
117 cout << "Show Nodes Coordinates : " << endl ;
119 cout << "Name :" << endl ;
120 const string * CoordinatesNames = myMesh->getCoordinatesNames() ;
121 for(int i=0; i<SpaceDimension ; i++) {
122 cout << " - " << CoordinatesNames[i] << endl ;
124 cout << "Unit :" << endl ;
125 const string * CoordinatesUnits = myMesh->getCoordinatesUnits() ;
126 for(int i=0; i<SpaceDimension ; i++) {
127 cout << " - " << CoordinatesUnits[i] << endl ;
129 for(int i=0; i<NumberOfNodes ; i++) {
130 cout << "Nodes " << i+1 << " : " ;
131 for (int j=0; j<SpaceDimension ; j++)
132 cout << Coordinates[i*SpaceDimension+j] << " " ;
136 int NumberOfTypes = myMesh->getNumberOfTypes(MED_CELL) ;
137 const medGeometryElement * Types = myMesh->getTypes(MED_CELL) ;
139 cout << "Show Connectivity (Nodal) :" << endl ;
140 for (int i=0; i<NumberOfTypes; i++) {
141 cout << "For type " << Types[i] << " : " << endl ;
142 int NumberOfElements = myMesh->getNumberOfElements(MED_CELL,Types[i]);
143 const int * connectivity = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,Types[i]);
144 int NomberOfNodesPerCell = Types[i]%100 ;
145 for (int j=0;j<NumberOfElements;j++){
146 cout << "Element "<< j+1 <<" : " ;
147 for (int k=0;k<NomberOfNodesPerCell;k++)
148 cout << connectivity[j*NomberOfNodesPerCell+k]<<" ";
153 cout << "Show Family :"<<endl ;
154 affiche_famille(myMesh,MED_NODE);
155 affiche_famille(myMesh,MED_CELL);
156 affiche_famille(myMesh,MED_FACE);
157 affiche_famille(myMesh,MED_EDGE);
159 cout << "Show Group :"<<endl ;
160 affiche_groupe(myMesh,MED_NODE);
161 affiche_groupe(myMesh,MED_CELL);
162 affiche_groupe(myMesh,MED_FACE);
163 affiche_groupe(myMesh,MED_EDGE);
165 cout << "Show Reverse Nodal Connectivity :" << endl ;
166 const int * ReverseNodalConnectivity = myMesh->getReverseConnectivity(MED_NODAL) ;
167 const int * ReverseNodalConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_NODAL) ;
168 for (int i=0; i<NumberOfNodes; i++) {
169 cout << "Node "<<i+1<<" : " ;
170 for (int j=ReverseNodalConnectivityIndex[i];j<ReverseNodalConnectivityIndex[i+1];j++)
171 cout << ReverseNodalConnectivity[j-1] << " " ;
175 cout << "Show Connectivity (Descending) :" << endl ;
176 int NumberOfElements ;
177 const int * connectivity ;
178 const int * connectivity_index ;
179 myMesh->calculateConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL);
181 NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
182 connectivity = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS);
183 connectivity_index = myMesh->getConnectivityIndex(MED_DESCENDING,MED_CELL);
185 catch (MEDEXCEPTION m) {
186 cout << m.what() << endl ;
189 for (int j=0;j<NumberOfElements;j++) {
190 cout << "Element "<<j+1<<" : " ;
191 for (int k=connectivity_index[j];k<connectivity_index[j+1];k++)
192 cout << connectivity[k-1]<<" ";
196 cout << "Show Reverse Descending Connectivity :" << endl ;
197 const int * ReverseDescendingConnectivity = myMesh->getReverseConnectivity(MED_DESCENDING) ;
198 const int * ReverseDescendingConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_DESCENDING) ;
200 int NumberOfConstituents = 0;
202 medEntityMesh constituentEntity ;
204 if (MeshDimension==3) {
205 constituent = "Face" ;
206 constituentEntity = MED_FACE ;
209 if (MeshDimension==2) {
210 constituent = "Edge" ;
211 constituentEntity = MED_EDGE ;
214 if (MeshDimension==1) {
215 INFOS("ERROR : MeshDimension = 1 !");
216 INFOS("We could not see Reverse Descending Connectivity.") ;
218 NumberOfConstituents = myMesh->getNumberOfElements (constituentEntity,MED_ALL_ELEMENTS);
219 for (int i=0; i<NumberOfConstituents; i++) {
220 cout << constituent <<i+1<<" : " ;
221 for (int j=ReverseDescendingConnectivityIndex[i];j<ReverseDescendingConnectivityIndex[i+1];j++)
222 cout << ReverseDescendingConnectivity[j-1] << " " ;
226 cout << "Show "<<constituent<<" Connectivity (Nodal) :" << endl ;
227 const int * face_connectivity = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,constituentEntity,MED_ALL_ELEMENTS);
228 const int * face_connectivity_index = myMesh->getConnectivityIndex(MED_NODAL,constituentEntity);
229 for (int i=0; i<NumberOfConstituents; i++) {
230 cout << constituent <<i+1<<" : " ;
231 for (int j=face_connectivity_index[i];j<face_connectivity_index[i+1];j++)
232 cout << face_connectivity[j-1]<<" ";
236 /* test of normal, area, volume, barycenter */
238 SUPPORT * support1 = (SUPPORT*) NULL;
240 //FIELD<double>* normal = new FIELD<double>::FIELD();
241 //FIELD<double>* length = new FIELD<double>::FIELD();
244 string support_name = "Support on all " ;
245 support_name+=constituent;
246 support1 = new SUPPORT(myMesh,support_name,constituentEntity);
247 cout << "Building of the Support on all cells dimensionned (Meshdim-1) of the mesh :"<< endl ;
248 cout << "Face in 3D or Edge in 2D" << endl;
250 cout << "Getting the normal of each face of this support !" << endl ;
252 FIELD<double>* normal = myMesh->getNormal(support1);
254 double normal_square, norm ;
255 double maxnorm=-infty;
256 double minnorm=infty;
258 for (int i = 1; i<=NumberOfConstituents;i++) {
260 cout << "Normal " << i << " " ;
261 for (int j=1; j<=SpaceDimension; j++) {
262 tmp_value = normal->getValueIJ(i,j) ;
263 normal_square += tmp_value*tmp_value ;
264 cout << tmp_value << " " ;
266 norm = sqrt(normal_square);
267 maxnorm = dmax(maxnorm,norm);
268 minnorm = dmin(minnorm,norm);
269 cout << ", Norm = " << norm << endl;
271 cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
275 if (SpaceDimension == 2)
277 cout << "Getting the length of each edge !" << endl ;
279 FIELD<double>* length = myMesh->getLength(support1);
281 double length_value,maxlength,minlength;
284 for (int i = 1; i<=NumberOfConstituents;i++)
286 length_value = length->getValueIJ(i,1) ;
287 cout << "Length " << i << " " << length_value << endl;
288 maxlength = dmax(maxlength,length_value);
289 minlength = dmin(minlength,length_value);
291 cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
298 cout << "Building of the Support on all space-dimensionned cells of the mesh :"<< endl ;
299 SUPPORT * support = new SUPPORT(myMesh);
301 cout << "Getting the barycenter of each element of this support !" << endl ;
303 //FIELD<double>* barycenter = new FIELD<double>::FIELD();
305 FIELD<double>* barycenter = myMesh->getBarycenter(support);
306 NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
308 for (int i = 1; i<=NumberOfElements;i++)
310 if (SpaceDimension == 3)
311 cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
313 if (SpaceDimension == 2)
314 cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << endl;
319 //FIELD<double>* volume = new FIELD<double>::FIELD();
320 //FIELD<double>* area = new FIELD<double>::FIELD();
324 if (SpaceDimension == 3)
326 cout << "Getting the Volume of each element of this support which is a 3D one !" << endl;
328 FIELD<double>* volume = myMesh->getVolume(support);
330 double maxvol,minvol,voltot;
334 for (int i = 1; i<=NumberOfElements;i++)
336 cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
337 maxvol = dmax(maxvol,volume->getValueIJ(i,1));
338 minvol = dmin(minvol,volume->getValueIJ(i,1));
339 voltot = voltot + volume->getValueIJ(i,1);
342 cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
343 cout << "Support Volume " << voltot << endl;
347 else if (SpaceDimension == 2)
349 cout << "Getting the Area of each element of this support which is a 2D one !" << endl;
351 FIELD<double>* area = myMesh->getArea(support);
353 // cout << "nb of comp "<< area->getNumberOfComponents() << " length " << area->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS) << endl;
355 double maxarea,minarea,areatot;
359 for (int i = 1; i<=NumberOfElements;i++)
361 cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
362 maxarea = dmax(maxarea,area->getValueIJ(i,1));
363 minarea = dmin(minarea,area->getValueIJ(i,1));
364 areatot = areatot + area->getValueIJ(i,1);
367 cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
368 cout << "Support Area " << areatot << endl;
375 //if (barycenter != NULL) delete barycenter;
376 //if (volume != NULL ) delete volume;
377 //if (area != NULL ) delete area;
380 if (argc < 4) return 0;
384 if (argc != 4) exit(0) ;
385 // else we have a field !
387 string fieldname = argv[3];
389 // SUPPORT * mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
390 SUPPORT * mySupport = new SUPPORT(myMesh,"On_all_cell",MED_CELL);
391 FIELD<double> * myField= new FIELD<double>() ;
393 myField->setName(fieldname);
394 myField->setSupport(mySupport);
395 MED_FIELD_RDONLY_DRIVER<double> myFieldDriver(filename,myField) ;
396 myFieldDriver.setFieldName(fieldname);
397 myFieldDriver.open() ;
400 myFieldDriver.read() ;
403 mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
404 myField->setSupport(mySupport);
406 myFieldDriver.read() ;
408 cout << "Field " << fieldname << " not found !!!" << endl ;
413 myFieldDriver.close() ;
415 cout << "Field "<< myField->getName() << " : " <<myField->getDescription() << endl ;
416 int NumberOfComponents = myField->getNumberOfComponents() ;
417 cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
418 for (int i=1; i<NumberOfComponents+1; i++) {
419 cout << " - composante "<<i<<" :"<<endl ;
420 cout << " - nom : "<<myField->getComponentName(i)<< endl;
421 cout << " - description : "<<myField->getComponentDescription(i) << endl;
422 cout << " - units : "<<myField->getMEDComponentUnit(i) << endl;
424 cout << "- iteration :" << endl ;
425 cout << " - numero : " << myField->getIterationNumber()<< endl ;
426 cout << " - ordre : " << myField->getOrderNumber()<< endl ;
427 cout << " - temps : " << myField->getTime()<< endl ;
429 cout << "- Valeurs :"<<endl;
430 int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
431 // for (int i=1; i<NumberOfComponents+1; i++) {
432 // double * value = myField->getValueI(MED_NO_INTERLACE,i) ;
433 // for (int j=0; j<NumberOf; j++)
434 // cout << value[j]<< " ";
437 MEDARRAY<double> * myvalue = myField->getvalue();
438 const double * value ;
439 for (int i=1; i<NumberOf+1; i++) {
440 value = myvalue->getRow(i) ;
441 for (int j=0; j<NumberOfComponents; j++)
442 cout << value[j]<< " ";