1 // MED MEDMEM : MED files in memory
3 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
24 // File : med_test.cxx
33 #include "MEDMEM_Exception.hxx"
34 #include "MEDMEM_Mesh.hxx"
35 #include "MEDMEM_Family.hxx"
36 #include "MEDMEM_Group.hxx"
38 #include "MEDMEM_MedMeshDriver.hxx"
39 #include "MEDMEM_MedFieldDriver.hxx"
40 #include "MEDMEM_Support.hxx"
41 #include "MEDMEM_Field.hxx"
42 #include "MEDMEM_define.hxx"
45 double dmax(double x, double y) { return (x>y)?x:y;}
47 double dmin(double x, double y) { return (x>y)?y:x;}
51 void affiche_support(const SUPPORT * mySupport)
53 cout << " - Name : "<<mySupport->getName().c_str()<<endl ;
54 cout << " - Description : "<<mySupport->getDescription().c_str()<<endl ;
55 cout << " - Entity : "<<mySupport->getEntity()<<endl ;
56 cout << " - Entities list : "<<endl ;
57 if (!(mySupport->isOnAllElements())) {
58 int NumberOfTypes = mySupport->getNumberOfTypes() ;
59 cout<<" - NumberOfTypes : "<<NumberOfTypes<<endl;
60 const medGeometryElement * Types = mySupport->getTypes() ;
61 for (int j=0;j<NumberOfTypes;j++) {
62 cout << " * Type "<<Types[j]<<" : " ;
63 int NumberOfElements = mySupport->getNumberOfElements(Types[j]) ;
64 const int * Number = mySupport->getNumber(Types[j]) ;
65 for (int k=0; k<NumberOfElements;k++)
66 cout << Number[k] << " ";
70 cout << " Is on all entities !"<< endl;
74 void affiche_famille(MESH *myMesh,medEntityMesh Entity)
76 int NumberOfFamilies = myMesh->getNumberOfFamilies(Entity) ;
77 cout << "NumberOfFamilies : "<<NumberOfFamilies<<endl;
78 for (int i=1; i<NumberOfFamilies+1;i++) {
79 const FAMILY* myFamily = myMesh->getFamily(Entity,i);
80 affiche_support(myFamily);
81 cout << " - Identifier : "<<myFamily->getIdentifier()<<endl ;
82 int NumberOfAttributes = myFamily->getNumberOfAttributes() ;
83 cout << " - Attributes ("<<NumberOfAttributes<<") :"<<endl;
84 for (int j=1;j<NumberOfAttributes+1;j++)
85 cout << " * "<<myFamily->getAttributeIdentifier(j)<<" : "<<myFamily->getAttributeValue(j)<<", "<<myFamily->getAttributeDescription(j).c_str()<<endl ;
86 int NumberOfGroups = myFamily->getNumberOfGroups() ;
87 cout << " - Groups ("<<NumberOfGroups<<") :"<<endl;
88 for (int j=1;j<NumberOfGroups+1;j++)
89 cout << " * "<<myFamily->getGroupName(j).c_str()<<endl ;
93 void affiche_groupe(MESH *myMesh,medEntityMesh Entity)
95 int NumberOfGroups = myMesh->getNumberOfGroups(Entity) ;
96 cout << "NumberOfGroups : "<<NumberOfGroups<<endl;
97 for (int i=1; i<NumberOfGroups+1;i++) {
98 const GROUP* myGroup = myMesh->getGroup(Entity,i);
99 affiche_support(myGroup);
100 int NumberOfFamillies = myGroup->getNumberOfFamilies() ;
101 cout << " - Families ("<<NumberOfFamillies<<") :"<<endl;
102 for (int j=1;j<NumberOfFamillies+1;j++)
103 cout << " * "<<myGroup->getFamily(j)->getName().c_str()<<endl ;
107 int main (int argc, char ** argv) {
111 if ((argc !=3) && (argc != 4)) {
112 cerr << "Usage : " << argv[0]
113 << " filename meshname [fieldname]" << endl << endl;
117 string filename = argv[1] ;
118 string meshname = argv[2] ;
120 // MESH * myMesh= new MESH(MED_DRIVER,filename,meshname) ;
121 MESH * myMesh= new MESH() ;
122 myMesh->setName(meshname);
123 MED_MESH_RDONLY_DRIVER myMeshDriver(filename,myMesh) ;
124 myMeshDriver.setMeshName(meshname);
125 myMeshDriver.open() ;
126 myMeshDriver.read() ;
127 myMeshDriver.close() ;
129 // int drv = myMesh->addDriver(MED_DRIVER,"sortie.med",meshname);
130 // myMesh->write(drv);
133 int SpaceDimension = myMesh->getSpaceDimension() ;
134 int MeshDimension = myMesh->getMeshDimension() ;
135 int NumberOfNodes = myMesh->getNumberOfNodes() ;
137 cout << "Space Dimension : " << SpaceDimension << endl << endl ;
139 cout << "Mesh Dimension : " << MeshDimension << endl << endl ;
141 const double * Coordinates = myMesh->getCoordinates(MED_FULL_INTERLACE) ;
143 cout << "Show Nodes Coordinates : " << endl ;
145 cout << "Name :" << endl ;
146 const string * CoordinatesNames = myMesh->getCoordinatesNames() ;
147 for(int i=0; i<SpaceDimension ; i++) {
148 cout << " - " << CoordinatesNames[i] << endl ;
150 cout << "Unit :" << endl ;
151 const string * CoordinatesUnits = myMesh->getCoordinatesUnits() ;
152 for(int i=0; i<SpaceDimension ; i++) {
153 cout << " - " << CoordinatesUnits[i] << endl ;
155 for(int i=0; i<NumberOfNodes ; i++) {
156 cout << "Nodes " << i+1 << " : " ;
157 for (int j=0; j<SpaceDimension ; j++)
158 cout << Coordinates[i*SpaceDimension+j] << " " ;
162 int NumberOfTypes = myMesh->getNumberOfTypes(MED_CELL) ;
163 const medGeometryElement * Types = myMesh->getTypes(MED_CELL) ;
165 cout << "Show Connectivity (Nodal) :" << endl ;
166 for (int i=0; i<NumberOfTypes; i++) {
167 cout << "For type " << Types[i] << " : " << endl ;
168 int NumberOfElements = myMesh->getNumberOfElements(MED_CELL,Types[i]);
169 const int * connectivity = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,Types[i]);
170 int NomberOfNodesPerCell = Types[i]%100 ;
171 for (int j=0;j<NumberOfElements;j++){
172 cout << "Element "<< j+1 <<" : " ;
173 for (int k=0;k<NomberOfNodesPerCell;k++)
174 cout << connectivity[j*NomberOfNodesPerCell+k]<<" ";
179 cout << "Show Family :"<<endl ;
180 affiche_famille(myMesh,MED_NODE);
181 affiche_famille(myMesh,MED_CELL);
182 affiche_famille(myMesh,MED_FACE);
183 affiche_famille(myMesh,MED_EDGE);
185 cout << "Show Group :"<<endl ;
186 affiche_groupe(myMesh,MED_NODE);
187 affiche_groupe(myMesh,MED_CELL);
188 affiche_groupe(myMesh,MED_FACE);
189 affiche_groupe(myMesh,MED_EDGE);
191 cout << "Show Reverse Nodal Connectivity :" << endl ;
192 const int * ReverseNodalConnectivity = myMesh->getReverseConnectivity(MED_NODAL) ;
193 const int * ReverseNodalConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_NODAL) ;
194 for (int i=0; i<NumberOfNodes; i++) {
195 cout << "Node "<<i+1<<" : " ;
196 for (int j=ReverseNodalConnectivityIndex[i];j<ReverseNodalConnectivityIndex[i+1];j++)
197 cout << ReverseNodalConnectivity[j-1] << " " ;
201 cout << "Show Connectivity (Descending) :" << endl ;
202 int NumberOfElements ;
203 const int * connectivity ;
204 const int * connectivity_index ;
205 myMesh->calculateConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL);
207 NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
208 connectivity = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_DESCENDING,MED_CELL,MED_ALL_ELEMENTS);
209 connectivity_index = myMesh->getConnectivityIndex(MED_DESCENDING,MED_CELL);
211 catch (MEDEXCEPTION m) {
212 cout << m.what() << endl ;
215 for (int j=0;j<NumberOfElements;j++) {
216 cout << "Element "<<j+1<<" : " ;
217 for (int k=connectivity_index[j];k<connectivity_index[j+1];k++)
218 cout << connectivity[k-1]<<" ";
222 cout << "Show Reverse Descending Connectivity :" << endl ;
223 const int * ReverseDescendingConnectivity = myMesh->getReverseConnectivity(MED_DESCENDING) ;
224 const int * ReverseDescendingConnectivityIndex = myMesh->getReverseConnectivityIndex(MED_DESCENDING) ;
226 int NumberOfConstituents = 0;
228 medEntityMesh constituentEntity ;
230 if (MeshDimension==3) {
231 constituent = "Face" ;
232 constituentEntity = MED_FACE ;
235 if (MeshDimension==2) {
236 constituent = "Edge" ;
237 constituentEntity = MED_EDGE ;
240 if (MeshDimension==1) {
241 INFOS("ERROR : MeshDimension = 1 !");
242 INFOS("We could not see Reverse Descending Connectivity.") ;
244 NumberOfConstituents = myMesh->getNumberOfElements (constituentEntity,MED_ALL_ELEMENTS);
245 for (int i=0; i<NumberOfConstituents; i++) {
246 cout << constituent <<i+1<<" : " ;
247 for (int j=ReverseDescendingConnectivityIndex[i];j<ReverseDescendingConnectivityIndex[i+1];j++)
248 cout << ReverseDescendingConnectivity[j-1] << " " ;
252 cout << "Show "<<constituent<<" Connectivity (Nodal) :" << endl ;
253 const int * face_connectivity = myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,constituentEntity,MED_ALL_ELEMENTS);
254 const int * face_connectivity_index = myMesh->getConnectivityIndex(MED_NODAL,constituentEntity);
255 for (int i=0; i<NumberOfConstituents; i++) {
256 cout << constituent <<i+1<<" : " ;
257 for (int j=face_connectivity_index[i];j<face_connectivity_index[i+1];j++)
258 cout << face_connectivity[j-1]<<" ";
262 /* test of normal, area, volume, barycenter */
264 SUPPORT * support1 = (SUPPORT*) NULL;
266 //FIELD<double>* normal = new FIELD<double>::FIELD();
267 //FIELD<double>* length = new FIELD<double>::FIELD();
270 string support_name = "Support on all " ;
271 support_name+=constituent;
272 support1 = new SUPPORT(myMesh,support_name,constituentEntity);
273 cout << "Building of the Support on all cells dimensionned (Meshdim-1) of the mesh :"<< endl ;
274 cout << "Face in 3D or Edge in 2D" << endl;
276 cout << "Getting the normal of each face of this support !" << endl ;
278 FIELD<double>* normal = myMesh->getNormal(support1);
280 double normal_square, norm ;
281 double maxnorm=-infty;
282 double minnorm=infty;
284 for (int i = 1; i<=NumberOfConstituents;i++) {
286 cout << "Normal " << i << " " ;
287 for (int j=1; j<=SpaceDimension; j++) {
288 tmp_value = normal->getValueIJ(i,j) ;
289 normal_square += tmp_value*tmp_value ;
290 cout << tmp_value << " " ;
292 norm = sqrt(normal_square);
293 maxnorm = dmax(maxnorm,norm);
294 minnorm = dmin(minnorm,norm);
295 cout << ", Norm = " << norm << endl;
297 cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
301 if (SpaceDimension == 2)
303 cout << "Getting the length of each edge !" << endl ;
305 FIELD<double>* length = myMesh->getLength(support1);
307 double length_value,maxlength,minlength;
310 for (int i = 1; i<=NumberOfConstituents;i++)
312 length_value = length->getValueIJ(i,1) ;
313 cout << "Length " << i << " " << length_value << endl;
314 maxlength = dmax(maxlength,length_value);
315 minlength = dmin(minlength,length_value);
317 cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
324 cout << "Building of the Support on all space-dimensionned cells of the mesh :"<< endl ;
325 SUPPORT * support = new SUPPORT(myMesh);
327 cout << "Getting the barycenter of each element of this support !" << endl ;
329 //FIELD<double>* barycenter = new FIELD<double>::FIELD();
331 FIELD<double>* barycenter = myMesh->getBarycenter(support);
332 NumberOfElements = myMesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
334 for (int i = 1; i<=NumberOfElements;i++)
336 if (SpaceDimension == 3)
337 cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
339 if (SpaceDimension == 2)
340 cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << endl;
345 //FIELD<double>* volume = new FIELD<double>::FIELD();
346 //FIELD<double>* area = new FIELD<double>::FIELD();
350 if (SpaceDimension == 3)
352 cout << "Getting the Volume of each element of this support which is a 3D one !" << endl;
354 FIELD<double>* volume = myMesh->getVolume(support);
356 double maxvol,minvol,voltot;
360 for (int i = 1; i<=NumberOfElements;i++)
362 cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
363 maxvol = dmax(maxvol,volume->getValueIJ(i,1));
364 minvol = dmin(minvol,volume->getValueIJ(i,1));
365 voltot = voltot + volume->getValueIJ(i,1);
368 cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
369 cout << "Support Volume " << voltot << endl;
373 else if (SpaceDimension == 2)
375 cout << "Getting the Area of each element of this support which is a 2D one !" << endl;
377 FIELD<double>* area = myMesh->getArea(support);
379 // cout << "nb of comp "<< area->getNumberOfComponents() << " length " << area->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS) << endl;
381 double maxarea,minarea,areatot;
385 for (int i = 1; i<=NumberOfElements;i++)
387 cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
388 maxarea = dmax(maxarea,area->getValueIJ(i,1));
389 minarea = dmin(minarea,area->getValueIJ(i,1));
390 areatot = areatot + area->getValueIJ(i,1);
393 cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
394 cout << "Support Area " << areatot << endl;
401 //if (barycenter != NULL) delete barycenter;
402 //if (volume != NULL ) delete volume;
403 //if (area != NULL ) delete area;
406 if (argc < 4) return 0;
410 if (argc != 4) exit(0) ;
411 // else we have a field !
413 string fieldname = argv[3];
415 // SUPPORT * mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
416 SUPPORT * mySupport = new SUPPORT(myMesh,"On_all_cell",MED_CELL);
417 FIELD<double> * myField= new FIELD<double>() ;
419 myField->setName(fieldname);
420 myField->setSupport(mySupport);
421 MED_FIELD_RDONLY_DRIVER<double> myFieldDriver(filename,myField) ;
422 myFieldDriver.setFieldName(fieldname);
423 myFieldDriver.open() ;
426 myFieldDriver.read() ;
429 mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
430 myField->setSupport(mySupport);
432 myFieldDriver.read() ;
434 cout << "Field " << fieldname << " not found !!!" << endl ;
439 myFieldDriver.close() ;
441 cout << "Field "<< myField->getName() << " : " <<myField->getDescription() << endl ;
442 int NumberOfComponents = myField->getNumberOfComponents() ;
443 cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
444 for (int i=1; i<NumberOfComponents+1; i++) {
445 cout << " - composante "<<i<<" :"<<endl ;
446 cout << " - nom : "<<myField->getComponentName(i)<< endl;
447 cout << " - description : "<<myField->getComponentDescription(i) << endl;
448 cout << " - units : "<<myField->getMEDComponentUnit(i) << endl;
450 cout << "- iteration :" << endl ;
451 cout << " - numero : " << myField->getIterationNumber()<< endl ;
452 cout << " - ordre : " << myField->getOrderNumber()<< endl ;
453 cout << " - temps : " << myField->getTime()<< endl ;
455 cout << "- Valeurs :"<<endl;
456 int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
457 // for (int i=1; i<NumberOfComponents+1; i++) {
458 // double * value = myField->getValueI(MED_NO_INTERLACE,i) ;
459 // for (int j=0; j<NumberOf; j++)
460 // cout << value[j]<< " ";
463 MEDARRAY<double> * myvalue = myField->getvalue();
464 const double * value ;
465 for (int i=1; i<NumberOf+1; i++) {
466 value = myvalue->getRow(i) ;
467 for (int j=0; j<NumberOfComponents; j++)
468 cout << value[j]<< " ";