Back to index

salome-med  6.5.0
MEDMEMTest_MeshAndMeshing.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDMEMTest.hxx"
00021 #include "MEDMEM_Meshing.hxx"
00022 #include "MEDMEM_Group.hxx"
00023 #include "MEDMEM_define.hxx"
00024 #include "MEDMEM_MedMeshDriver.hxx"
00025 #include "MEDMEM_Field.hxx"
00026 #include "MEDMEM_Grid.hxx"
00027 
00028 #include <sstream>
00029 #include <cmath>
00030 
00031 #include <cppunit/Message.h>
00032 #include <cppunit/TestAssert.h>
00033 
00034 // use this define to enable lines, execution of which leads to Segmentation Fault
00035 //#define ENABLE_FAULTS
00036 
00037 // use this define to enable CPPUNIT asserts and fails, showing bugs
00038 //#define ENABLE_FORCED_FAILURES
00039 
00040 using namespace std;
00041 using namespace MEDMEM;
00042 using namespace MED_EN;
00043 
00044 static void addMedFacesGroup (MESHING& meshing, int nFaces, const int *groupValue,
00045                               string groupName, const MED_EN::medGeometryElement *mytypes,
00046                               const int *index, const int *myNumberOfElements, int nbOfGeomTypes)
00047 {
00048   GROUP *faces=new GROUP;
00049   faces->setName(groupName);
00050   faces->setMesh(&meshing);
00051   faces->setEntity(MED_EN::MED_FACE);
00052   faces->setNumberOfGeometricType(nbOfGeomTypes);
00053   faces->setGeometricType(mytypes);
00054   faces->setNumberOfElements(myNumberOfElements);
00055   faces->setNumber(index, groupValue);
00056   meshing.addGroup(*faces);
00057   faces->removeReference();
00058 }
00059 
00197 void MEDMEMTest::testMeshAndMeshing()
00198 {
00199   string filename = getResourceFile("pointe.med");
00200   string meshname = "maa1";
00201   string filenameout21        = makeTmpFile("myMeshWrite4_pointe21.med");
00202   string filename_profiles_wr = makeTmpFile("myMedProfilesFieldfile.med");
00203   string filenameout          = makeTmpFile("out.med");
00204 
00205   // To remove tmp files from disk
00206   MEDMEMTest_TmpFilesRemover aRemover;
00207   aRemover.Register(filenameout21);
00208   aRemover.Register(filename_profiles_wr);
00209   aRemover.Register(filenameout);
00210 
00212   // TEST 1 //
00214 
00215   MESH * myMesh= new MESH();
00216   myMesh->setName("FIRST_MESH");
00217   CPPUNIT_ASSERT(myMesh != NULL);
00218 
00219   //test operator <<
00220   {
00221     ostringstream out;
00222     CPPUNIT_ASSERT_NO_THROW(out << *myMesh << endl);
00223   }
00224 
00225   //test operator =
00226   MESH *myMesh1 =new MESH( *myMesh);
00227 
00228   //deepCompare
00229   bool isEqual = false;
00230   CPPUNIT_ASSERT_NO_THROW(isEqual = myMesh1->deepCompare(*myMesh));
00231   CPPUNIT_ASSERT(isEqual);
00232   myMesh1->removeReference();
00233 
00234   //ensure it imposible to compare meshes
00235   MESH *myMeshPointer =  myMesh;
00236   //test operator ==
00237   CPPUNIT_ASSERT(*myMeshPointer == *myMesh);
00238 
00239   myMesh->removeReference();
00240 
00241   //set a MESH object
00242   MESHING *myMeshing=new MESHING;
00243   myMeshing->setName("meshing");
00244   // define coordinates
00245 
00246   int SpaceDimension = 3;
00247   int NumberOfNodes = 19;
00248   double Coordinates[57] = 
00249     {
00250       0.0, 0.0, 0.0,
00251       0.0, 0.0, 1.0,
00252       2.0, 0.0, 1.0,
00253       0.0, 2.0, 1.0,
00254       -2.0, 0.0, 1.0,
00255       0.0, -2.0, 1.0,
00256       1.0, 1.0, 2.0,
00257       -1.0, 1.0, 2.0,
00258       -1.0, -1.0, 2.0,
00259       1.0, -1.0, 2.0,
00260       1.0, 1.0, 3.0,
00261       -1.0, 1.0, 3.0,
00262       -1.0, -1.0, 3.0,
00263       1.0, -1.0, 3.0,
00264       1.0, 1.0, 4.0,
00265       -1.0, 1.0, 4.0,
00266       -1.0, -1.0, 4.0,
00267       1.0, -1.0, 4.0,
00268       0.0, 0.0, 5.0
00269     };
00270   try
00271     {
00272       myMeshing->setCoordinates(SpaceDimension,NumberOfNodes,Coordinates,"CARTESIAN",MED_FULL_INTERLACE);
00273     }
00274   catch (const std::exception &e)
00275     {
00276       CPPUNIT_FAIL(e.what());
00277     }
00278   catch (...)
00279     {
00280       CPPUNIT_FAIL("Unknown exception");
00281     }
00282 
00283   string Names[3] = 
00284     {
00285       "X","Y","Z"
00286     };
00287   try
00288     {
00289       myMeshing->setCoordinatesNames(Names);
00290     }
00291   catch (const std::exception &e)
00292     {
00293       CPPUNIT_FAIL(e.what());
00294     }
00295   catch (...)
00296     {
00297       CPPUNIT_FAIL("Unknown exception");
00298     }
00299 
00300   string Units[3] = 
00301     {
00302       "cm","cm","cm"
00303     };
00304   try
00305     {
00306       myMeshing->setCoordinatesUnits(Units);
00307     }
00308   catch (const std::exception &e)
00309     {
00310       CPPUNIT_FAIL(e.what());
00311     }
00312   catch (...)
00313     {
00314       CPPUNIT_FAIL("Unknown exception");
00315     }
00316 
00317   // define conectivities
00318 
00319   // cell part
00320 
00321   const int NumberOfTypes = 3;
00322   medGeometryElement Types[NumberOfTypes] = 
00323     {
00324       MED_TETRA4,MED_PYRA5,MED_HEXA8
00325     };
00326   const int NumberOfElements[NumberOfTypes] = 
00327     {
00328       12,2,2
00329     };
00330 
00331   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfTypes(NumberOfTypes,MED_CELL));
00332 
00333   CPPUNIT_ASSERT_NO_THROW(myMeshing->setTypes(Types,MED_CELL));
00334 
00335   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfElements(NumberOfElements,MED_CELL));
00336 
00337   const int sizeTetra = 12*4;
00338   int ConnectivityTetra[sizeTetra]=
00339     {
00340       1,2,3,6,
00341       1,2,4,3,
00342       1,2,5,4,
00343       1,2,6,5,
00344       2,7,4,3,
00345       2,8,5,4,
00346       2,9,6,5,
00347       2,10,3,6,
00348       2,7,3,10,
00349       2,8,4,7,
00350       2,9,5,8,
00351       2,10,6,9
00352     };
00353 
00354   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity(MED_CELL,MED_TETRA4,ConnectivityTetra));
00355 
00356   int ConnectivityPyra[2*5]=
00357     {
00358       7,8,9,10,2,
00359       15,18,17,16,19
00360     };
00361 
00362   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity(MED_CELL,MED_PYRA5,ConnectivityPyra));
00363 
00364   int ConnectivityHexa[2*8]=
00365     {
00366       11,12,13,14,7,8,9,10,
00367       15,16,17,18,11,12,13,14
00368     };
00369 
00370   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity(MED_CELL,MED_HEXA8,ConnectivityHexa));
00371 
00372   // face part
00373   const int NumberOfFacesTypes = 2;
00374   medGeometryElement FacesTypes[NumberOfFacesTypes] = 
00375     {
00376       MED_TRIA3,MED_QUAD4
00377     };
00378   const int NumberOfFacesElements[NumberOfFacesTypes] = 
00379     {
00380       4,4
00381     };
00382 
00383   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfTypes(NumberOfFacesTypes,MED_FACE));
00384   CPPUNIT_ASSERT_NO_THROW(myMeshing->setTypes(FacesTypes,MED_FACE));
00385   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfElements(NumberOfFacesElements,MED_FACE));
00386   const int nbTria = 4;
00387   const int sizeTria = nbTria*3;
00388   int ConnectivityTria[sizeTria]=
00389     {
00390       1,4,3,
00391       1,5,4,
00392       1,6,5,
00393       1,3,6
00394     };
00395 
00396   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity(MED_FACE,MED_TRIA3,ConnectivityTria));
00397   const int nbQua = 4;
00398   int ConnectivityQua[nbQua*4]=
00399     {
00400       7,8,9,10,
00401       11,12,13,14,
00402       11,7,8,12,
00403       12,8,9,13
00404     };
00405 
00406   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity(MED_FACE,MED_QUAD4,ConnectivityQua));
00407 
00408   int meshDimension = SpaceDimension; // because there 3D cells in the mesh
00409   //CPPUNIT_ASSERT_NO_THROW(myMeshing->setMeshDimension(meshDimension));
00410 
00411   // edge part
00412 
00413   // not yet implemented : if set, results are unpredictable.
00414 
00415   // Some groups :
00416 
00417   // Node :
00418   {
00419     GROUP *myGroup=new GROUP;
00420     myGroup->setName("SomeNodes");
00421     myGroup->setMesh(myMeshing);
00422     myGroup->setEntity(MED_NODE);
00423     myGroup->setNumberOfGeometricType(1);
00424     medGeometryElement myTypes[1] = 
00425       {
00426         MED_NONE
00427       };
00428     myGroup->setGeometricType(myTypes);
00429     const int myNumberOfElements[1] = 
00430       {
00431         4
00432       };
00433     myGroup->setNumberOfElements(myNumberOfElements);
00434     const int index[1+1] = 
00435       {
00436         1,5
00437       };
00438     const int value[4]= 
00439       {
00440         1,4,5,7
00441       };
00442     myGroup->setNumber(index,value);
00443     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00444     myGroup->removeReference();
00445   }
00446   {
00447     GROUP *myGroup=new GROUP;
00448     myGroup->setName("OtherNodes");
00449     myGroup->setMesh(myMeshing);
00450     myGroup->setEntity(MED_NODE);
00451     myGroup->setNumberOfGeometricType(1);
00452     medGeometryElement myTypes[1] = 
00453       {
00454         MED_NONE
00455       };
00456     myGroup->setGeometricType(myTypes);
00457     const int myNumberOfElements[1] = 
00458       {
00459         3
00460       };
00461     myGroup->setNumberOfElements(myNumberOfElements);
00462     const int index[1+1] = 
00463       {
00464         1,4
00465       };
00466     const int value[3]= 
00467       {
00468         2,3,6
00469       };
00470     myGroup->setNumber(index,value);
00471     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00472     myGroup->removeReference();
00473   }
00474 
00475   // Cell :
00476   {
00477     GROUP *myGroup=new GROUP;
00478     myGroup->setName("SomeCells");
00479     myGroup->setMesh(myMeshing);
00480     myGroup->setEntity(MED_CELL);
00481     myGroup->setNumberOfGeometricType(3);
00482     medGeometryElement myTypes[3] = 
00483       {
00484         MED_TETRA4,MED_PYRA5,MED_HEXA8
00485       };
00486     myGroup->setGeometricType(myTypes);
00487     const int myNumberOfElements[3] = 
00488       {
00489         4,1,2
00490       };
00491     myGroup->setNumberOfElements(myNumberOfElements);
00492     const int index[3+1] = 
00493       {
00494         1,5,6,8
00495       };
00496     const int value[4+1+2]=
00497       {
00498         2,7,8,12,
00499         13,
00500         15,16
00501       };
00502     myGroup->setNumber(index,value);
00503     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00504     myGroup->removeReference();
00505   }
00506   {
00507     GROUP *myGroup=new GROUP;
00508     myGroup->setName("OtherCells");
00509     myGroup->setMesh(myMeshing);
00510     myGroup->setEntity(MED_CELL);
00511     myGroup->setNumberOfGeometricType(2);
00512     medGeometryElement myTypes[] = 
00513       {
00514         MED_TETRA4,MED_PYRA5
00515       };
00516     myGroup->setGeometricType(myTypes);
00517     const int myNumberOfElements[] = 
00518       {
00519         4,1
00520       };
00521     myGroup->setNumberOfElements(myNumberOfElements);
00522     const int index[2+1] = 
00523       {
00524         1,5,6
00525       };
00526     const int value[4+1]=
00527       {
00528         3,4,5,9,
00529         14
00530       };
00531     myGroup->setNumber(index,value);
00532     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00533     myGroup->removeReference();
00534   }
00535 
00536   // Face :
00537   {
00538     GROUP *myGroup=new GROUP;
00539     myGroup->setName("SomeFaces");
00540     myGroup->setMesh(myMeshing);
00541     myGroup->setEntity(MED_FACE);
00542     myGroup->setNumberOfGeometricType(2);
00543     medGeometryElement myTypes[2] = 
00544       {
00545         MED_TRIA3,MED_QUAD4
00546       };
00547     myGroup->setGeometricType(myTypes);
00548     const int myNumberOfElements[2] = 
00549       {
00550         2,3
00551       };
00552     myGroup->setNumberOfElements(myNumberOfElements);
00553     const int index[2+1] = 
00554       {
00555         1,3,6
00556       };
00557     const int value[2+3]=
00558       {
00559         2,4,
00560         5,6,8
00561       };
00562     myGroup->setNumber(index,value);
00563     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00564     myGroup->removeReference();
00565   }
00566   {
00567     GROUP *myGroup=new GROUP;
00568     myGroup->setName("OtherFaces");
00569     myGroup->setMesh(myMeshing);
00570     myGroup->setEntity(MED_FACE);
00571     myGroup->setNumberOfGeometricType(1);
00572     medGeometryElement myTypes[1] = 
00573       {
00574         MED_TRIA3
00575       };
00576     myGroup->setGeometricType(myTypes);
00577     const int myNumberOfElements[1] = 
00578       {
00579         2
00580       };
00581     myGroup->setNumberOfElements(myNumberOfElements);
00582     const int index[1+1] = 
00583       {
00584         1,3
00585       };
00586     const int value[2]=
00587       {
00588         1,3
00589       };
00590     myGroup->setNumber(index,value);
00591     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00592     myGroup->removeReference();
00593   }
00594 
00595   //test Mesh(MESH &m)
00596   {
00597     MESH * myMesh2 = new MESH( *myMeshing );
00598     CPPUNIT_ASSERT(myMesh2->deepCompare(*myMeshing));
00599     myMeshing->removeReference();
00600     //cout<<*myMesh2<<endl;
00601     ostringstream os;
00602     os << * myMesh2;
00603     CPPUNIT_ASSERT(os.str() != "");
00604 
00605     CPPUNIT_ASSERT_EQUAL(myMesh2->getName(),(string)"meshing");
00606     CPPUNIT_ASSERT((myMesh2->getDescription()).size() == 0);
00607     myMesh2->setDescription("This class contains all information related to a 'meshing' mesh ");
00608     CPPUNIT_ASSERT((myMesh2->getDescription()).size() != 0);
00609 
00610     CPPUNIT_ASSERT(myMesh2->getSpaceDimension() == SpaceDimension);
00611     CPPUNIT_ASSERT(myMesh2->getMeshDimension() == meshDimension);
00612     CPPUNIT_ASSERT(myMesh2->getNumberOfNodes() == NumberOfNodes);
00613 
00614     const COORDINATE* coord = myMesh2->getCoordinateptr();
00615     try
00616       {
00617         CPPUNIT_ASSERT(myMesh2->getCoordinatesSystem() != "catresian");
00618       }
00619     catch (const std::exception &e)
00620       {
00621         CPPUNIT_FAIL(e.what());
00622       }
00623     catch (...)
00624       {
00625         CPPUNIT_FAIL("Unknown exception");
00626       }
00627 
00628     const string * units;
00629     try
00630       {
00631         units = myMesh2->getCoordinatesUnits();
00632         for (int axe = 0; axe < SpaceDimension; axe++) 
00633           {
00634             string verif = coord->getCoordinateUnit(axe+1);
00635             CPPUNIT_ASSERT(verif == units[axe]);
00636           }
00637       }
00638     catch (const std::exception &e)
00639       {
00640         CPPUNIT_FAIL(e.what());
00641       }
00642     catch (...)
00643       {
00644         CPPUNIT_FAIL("Unknown exception");
00645       }
00646 
00647     const string * noms;
00648     try
00649       {
00650         noms = myMesh2->getCoordinatesNames();
00651         for (int axe = 0; axe < SpaceDimension; axe++) 
00652           {
00653             string verif = coord->getCoordinateName(axe+1);
00654             CPPUNIT_ASSERT(verif == noms[axe]);
00655           }
00656       }
00657     catch (const std::exception &e)
00658       {
00659         CPPUNIT_FAIL(e.what());
00660       }
00661     catch (...)
00662       {
00663         CPPUNIT_FAIL("Unknown exception");
00664       }
00665 
00666     try
00667       {
00668         const double * coor2 = myMesh2->getCoordinates(MED_FULL_INTERLACE);
00669 
00670         for (int axe = 0; axe < SpaceDimension; axe++) 
00671           {
00672             try
00673               {
00674                 for (int num = 0; num < NumberOfNodes; num++) 
00675                   {
00676                     try
00677                       {
00678                         const double d = myMesh2->getCoordinate(num + 1, axe + 1);
00679                         CPPUNIT_ASSERT(fabs(d - coor2[(num * SpaceDimension)+axe]) < 0.001);
00680                       }
00681                     catch (const std::exception &e)
00682                       {
00683                         CPPUNIT_FAIL(e.what());
00684                       }
00685                     catch (...)
00686                       {
00687                         CPPUNIT_FAIL("Unknown exception");
00688                       }
00689                   }
00690               }
00691             catch (const std::exception &e)
00692               {
00693                 CPPUNIT_FAIL(e.what());
00694               }
00695             catch (...)
00696               {
00697                 CPPUNIT_FAIL("Unknown exception");
00698               }
00699           }
00700       }
00701     catch (const std::exception &e)
00702       {
00703         CPPUNIT_FAIL(e.what());
00704       }
00705     catch (...)
00706       {
00707         CPPUNIT_FAIL("Unknown exception");
00708       }
00709 
00710     const CONNECTIVITY * myConnectivity = myMesh2->getConnectivityptr();
00711 
00712     // MED_EN::MED_CELL
00713     MED_EN::medEntityMesh entity = myConnectivity->getEntity();
00714     CPPUNIT_ASSERT_EQUAL(MED_CELL, entity);
00715 
00716     int typesNb;
00717     CPPUNIT_ASSERT_NO_THROW(typesNb= myConnectivity->getNumberOfTypes(entity));
00718     CPPUNIT_ASSERT_EQUAL(NumberOfTypes, typesNb);
00719 
00720     const MED_EN::medGeometryElement * Types1;
00721     CPPUNIT_ASSERT_NO_THROW(Types1 = myMesh2->getTypes(entity));
00722 
00723     medConnectivity myMedConnect;
00724     bool existConnect = false;
00725     if (myMesh2->existConnectivity(MED_NODAL, entity))
00726       {
00727         existConnect = true;
00728         myMedConnect = MED_NODAL;
00729       }
00730     else if(myMesh2->existConnectivity(MED_DESCENDING, entity))
00731       {
00732         existConnect = true;
00733         myMedConnect = MED_DESCENDING;
00734       }
00735 
00736     for(int t = 0; t < NumberOfTypes; t++ )
00737       {
00738         CPPUNIT_ASSERT_EQUAL(Types1[t], Types[t]);
00739         int NumberOfElements1 = 0;
00740         CPPUNIT_ASSERT_NO_THROW(NumberOfElements1 = myMesh2->getNumberOfElements(entity, Types1[t]));
00741         CPPUNIT_ASSERT_EQUAL(NumberOfElements1, NumberOfElements[t]);
00742         if(existConnect)
00743           {
00744             ostringstream out;
00745             const int * connectivity;
00746             const int * connectivity_index;
00747             CPPUNIT_ASSERT_NO_THROW(connectivity = myMesh2->getConnectivity (myMedConnect, entity, Types1[t]));
00748             connectivity_index = myMesh2->getConnectivityIndex(myMedConnect, entity);
00749             for (int j = 0; j < NumberOfElements1; j++) 
00750               {
00751                 out<<"!!!!!!!!!!!!!!!"<<endl;
00752                 for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
00753                   out << connectivity[k-1] << " ";
00754                 out << endl;
00755               }
00756           }
00757 
00758         const CELLMODEL* myCellModel = myMesh2->getCellsTypes(entity);
00759         string* TypeNames;
00760         CPPUNIT_ASSERT_NO_THROW(TypeNames = myMesh2->getCellTypeNames(entity));
00761 
00762         for(int k = 0; k < NumberOfTypes; k++ )
00763           {
00764             CPPUNIT_ASSERT_EQUAL(TypeNames[k], myCellModel[k].getName());
00765           }
00766         delete [] TypeNames;
00767 
00768         const int* myGlobalNbIdx;
00769         CPPUNIT_ASSERT_NO_THROW(myGlobalNbIdx = myMesh2->getGlobalNumberingIndex(MED_FACE));
00770         for(int i = 0; i <= NumberOfFacesTypes; i++)
00771           {
00772             if(i == NumberOfFacesTypes)
00773               {
00774                 CPPUNIT_ASSERT_EQUAL(myGlobalNbIdx[i],nbTria+nbQua+1);
00775                 CPPUNIT_ASSERT_THROW(myMesh2->getElementType(MED_FACE, myGlobalNbIdx[i]), MEDEXCEPTION);
00776                 break;
00777               }
00778             //cout<<"Global number of first element of each geom type : "<<myGlobalNbIdx[i]<<endl;
00779           }
00780         //cout<<"Global number of first element of each geom type : "<<myGlobalNbIdx[i]<<endl;
00781       }
00782     {
00783       const int * ReverseNodalConnectivity;
00784 
00785       ((CONNECTIVITY*)myConnectivity)->setNumberOfNodes(NumberOfNodes);
00786 
00787       CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivity = myMesh2->getReverseConnectivity(MED_NODAL, entity));
00788       CPPUNIT_ASSERT_NO_THROW(myMesh2->getReverseConnectivityLength(MED_NODAL, entity));
00789       const int * ReverseNodalConnectivityIndex = myMesh2->getReverseConnectivityIndex(MED_NODAL, entity);
00790       const int ReverseIdxLength = myMesh2->getReverseConnectivityIndexLength(MED_NODAL, entity);
00791       CPPUNIT_ASSERT(ReverseIdxLength == NumberOfNodes+1);
00792       ostringstream out;
00793       for (int i = 0; i < NumberOfNodes; i++) 
00794         {
00795           out << "Node "<< i+1 << " : ";
00796           for (int j = ReverseNodalConnectivityIndex[i]; j < ReverseNodalConnectivityIndex[i+1]; j++)
00797             out << ReverseNodalConnectivity[j-1] << " ";
00798           out << endl;
00799         }
00800 
00801       // Show Descending Connectivity
00802       int NumberOfElements1;
00803       const int * connectivity;
00804       const int * connectivity_index;
00805       myMesh2->calculateConnectivity( MED_DESCENDING, entity);
00806       try
00807         {
00808           NumberOfElements1 = myMesh2->getNumberOfElements(entity, MED_ALL_ELEMENTS);
00809           connectivity = myMesh2->getConnectivity( MED_DESCENDING, entity, MED_ALL_ELEMENTS);
00810           connectivity_index = myMesh2->getConnectivityIndex(MED_DESCENDING, entity);
00811         }
00812       catch (MEDEXCEPTION m) 
00813         {
00814           CPPUNIT_FAIL(m.what());
00815         }
00816 
00817       for (int j = 0; j < NumberOfElements1; j++) 
00818         {
00819           out << "Element " << j+1 << " : ";
00820           for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
00821             out << connectivity[k-1] << " ";
00822           out << endl;
00823         }
00824 
00825       // getElementNumber
00826       if (myMesh2->existConnectivity(MED_NODAL, MED_FACE)) 
00827         {
00828           int myTr[3] = 
00829             {
00830               1,5,4
00831             };
00832           CPPUNIT_ASSERT_NO_THROW(myMesh2->getElementNumber(MED_NODAL,MED_FACE,MED_TRIA3,myTr));
00833         }
00834     }
00835 
00836     //test family and group
00837     int NumberOfGroups;
00838     CPPUNIT_ASSERT_THROW(myMesh2->getNumberOfGroups(MED_ALL_ENTITIES), MEDEXCEPTION);
00839     CPPUNIT_ASSERT_NO_THROW(NumberOfGroups = myMesh2->getNumberOfGroups(MED_CELL));
00840     CPPUNIT_ASSERT_EQUAL(NumberOfGroups, 2);
00841     vector<GROUP*> groups;
00842     CPPUNIT_ASSERT_NO_THROW(groups = myMesh2->getGroups(MED_CELL));
00843     CPPUNIT_ASSERT(groups.size() != 0);
00844     for(int nb = 1; nb <= NumberOfGroups; nb++ )
00845       {
00846         const GROUP* group;
00847         CPPUNIT_ASSERT_NO_THROW(group = myMesh2->getGroup(MED_CELL, nb));
00848         CPPUNIT_ASSERT_EQUAL(group->getName(), groups[nb-1]->getName());
00849       }
00850 
00851     int NumberOfFamilies;
00852     CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies = myMesh2->getNumberOfFamilies(MED_CELL));
00853     CPPUNIT_ASSERT_MESSAGE("Current mesh hasn't Families", NumberOfFamilies == 0);
00854 
00855     //create families - it's not possible to create, becase not all entities are defined
00856     // EAP: the problem has been fixed for ENSIGHT Industrialization project
00857     CPPUNIT_ASSERT_NO_THROW( myMesh2->createFamilies() );
00858 
00859     /*CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies = myMesh2->getNumberOfFamilies(MED_CELL));
00860       CPPUNIT_ASSERT( NumberOfFamilies != 0);*/
00861 
00862     myMesh2->removeReference();
00863   }
00864 
00866   // TEST 2: Polygon and Polyhedron(only NODAL connectivity)  //
00868 
00869   double CoordinatesPoly[57] = 
00870     {
00871       2.0, 3.0, 2.0,
00872       3.0, 2.0, 2.0,
00873       4.0, 1.0, 2.0,
00874       2.0, 0.0, 2.0,
00875       0.0, 1.0, 2.0,
00876       1.0, 2.0, 2.0,
00877       2.0, 3.0, 1.0,
00878       3.0, 2.0, 0.0,
00879       4.0, 1.0, 0.0,
00880       2.0, 0.0, 0.0,
00881       0.0, 1.0, 0.0,
00882       1.0, 2.0, 0.0,
00883       5.0, 3.0, 2.0,
00884       7.0, 2.0, 2.0,
00885       6.0, 0.0, 2.0,
00886       6.0, 3.0, 0.0,
00887       7.0, 2.0, 0.0,
00888       6.0, 0.0, -1.0,
00889       5.0, 1.0, -3.0
00890     };
00891 
00892   const int REFnodalConnOfFaces[91] = 
00893     {
00894       1, 2, 3, 4, 5, 6, -1,// Polyhedron 1
00895       1, 7, 8, 2,       -1,
00896       2, 8, 9, 3,       -1,
00897       4, 3, 9, 10,      -1,
00898       5, 4, 10, 11,     -1,
00899       6, 5, 11, 12,     -1,
00900       1, 6, 12, 7,      -1,
00901       7, 12, 8, 10,     -1,
00902       9, 8, 12, 11,     
00903 
00904       13, 14, 15, 3, 2, -1,// Polyhedron 2
00905       13, 2, 8, 16,     -1,
00906       14, 13, 16, 17,   -1,
00907       15, 14, 17, 15,   -1,
00908       17, 18, 15,       -1,
00909       18, 9, 3,         -1,
00910       15, 9, 2,         -1,
00911       3, 9, 8,          -1,
00912       8, 9, 17, 16,     -1,
00913       9, 18, 17
00914     };
00915   const int NumberOfFaces = 19;
00916   const int NumberOfPolyhedron = 2;
00917   const int nbOfPolygons = 2;
00918   const int REFpolyIndex[NumberOfPolyhedron+1] = 
00919     {
00920       1,47,92
00921     };
00922 
00923   double PolygonCoordinates[27] = 
00924     {
00925       2.0, 3.0, 12.0,
00926       3.0, 2.0, 12.0,
00927       4.0, 1.0, 12.0,
00928       2.0, 0.0, 12.0,
00929       0.0, 1.0, 12.0,
00930       1.0, 2.0, 12.0,
00931       5.0, 3.0, 12.0,
00932       7.0, 2.0, 12.0,
00933       6.0, 0.0, 12.0
00934     };
00935 
00936   const int REFpolygonFaces[11] = 
00937     {
00938       1, 2, 3, 4, 5, 6, // Polygon 1
00939       7, 8, 9, 3, 2     // Polygon 2
00940     };
00941 
00942   const int REFpolygonIndex[nbOfPolygons+1] = 
00943     {
00944       1, 7, 12
00945     };
00946 
00947   MESHING *myMeshingPoly=new MESHING;
00948   myMeshingPoly->setName("meshingpoly");
00949 
00950   const int NbOfTypes = 2;
00951   medGeometryElement TypesPoly[NbOfTypes] = 
00952     {
00953       MED_TETRA4, MED_POLYHEDRA
00954     };
00955   const int NbOfElements[NbOfTypes] = 
00956     {
00957       1,2
00958     };
00959 
00960   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setNumberOfTypes(NbOfTypes, MED_CELL));
00961 
00962   CPPUNIT_ASSERT_NO_THROW
00963     (myMeshingPoly->setCoordinates(SpaceDimension, NumberOfNodes, CoordinatesPoly,
00964                                    "CARTESIAN", MED_FULL_INTERLACE));
00965 
00966   //CPPUNIT_ASSERT_NO_THROW( myMeshingPoly->setSpaceDimension(SpaceDimension));
00967 
00968   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setTypes(TypesPoly, MED_CELL));
00969   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setNumberOfElements(NbOfElements, MED_CELL));
00970 
00971   string Unit ="cm";
00972   for(int i = 0; i < SpaceDimension; i++ )
00973     {
00974       CPPUNIT_ASSERT_NO_THROW( myMeshingPoly->setCoordinateName(Names[i],i) );
00975       CPPUNIT_ASSERT_NO_THROW( myMeshingPoly->setCoordinateUnit(Unit, i) );
00976     }
00977 
00978   int ConnectivityTetraPoly[4*1]=
00979     {
00980       17, 9, 18, 19
00981     };
00982 
00983   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setConnectivity(MED_CELL, MED_TETRA4,ConnectivityTetraPoly));
00984 
00985   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setConnectivity(MED_CELL, MED_POLYHEDRA,REFnodalConnOfFaces,REFpolyIndex));
00986 
00987   bool PolyConn = false;
00988   CPPUNIT_ASSERT_NO_THROW(PolyConn = myMeshingPoly->existConnectivity(MED_NODAL, MED_CELL));
00989   CPPUNIT_ASSERT(PolyConn);
00990   {
00991     CPPUNIT_ASSERT_EQUAL(NumberOfPolyhedron,
00992                          myMeshingPoly->getNumberOfElements(MED_CELL,MED_POLYHEDRA));
00993     CPPUNIT_ASSERT_NO_THROW( myMeshingPoly->calculateConnectivity (MED_NODAL,MED_FACE));
00994     CPPUNIT_ASSERT_EQUAL(NumberOfFaces-1, myMeshingPoly->getNumberOfElements(MED_FACE,MED_POLYGON)); // -1: one face is shared with tetra
00995     CPPUNIT_ASSERT_EQUAL(91,myMeshingPoly->getConnectivityLength(MED_NODAL,MED_CELL,MED_POLYHEDRA));
00996     const int * PolyConn;
00997     const int * PolyIdx;
00998     CPPUNIT_ASSERT_NO_THROW(PolyConn = myMeshingPoly->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS));
00999     CPPUNIT_ASSERT_NO_THROW(PolyIdx = myMeshingPoly->getConnectivityIndex(MED_NODAL,MED_CELL));
01000     for(int i = NbOfElements[0], iRef=0; i<NbOfElements[0]+NumberOfPolyhedron; i++)
01001       {
01002         int NodeIdxBegin = PolyIdx[i];
01003         int NodeIdxEnd = PolyIdx[i+1];
01004         for(int k = NodeIdxBegin; k < NodeIdxEnd; k++)
01005           CPPUNIT_ASSERT_EQUAL(REFnodalConnOfFaces[iRef++], PolyConn[k-1]);
01006       }
01007   }
01008 
01009     MESHING *myPolygonMeshing=new MESHING;
01010     myPolygonMeshing->setName("PolygonMeshing");
01011 
01012     medGeometryElement PolygonTypes[NbOfTypes] = 
01013       {
01014         MED_TRIA3,MED_POLYGON
01015       };
01016     const int PolygonNumberOfElements[NbOfTypes] = 
01017       {
01018         2,nbOfPolygons
01019       };
01020 
01021     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setNumberOfTypes(NbOfTypes, MED_CELL));
01022 
01023     try
01024     {
01025       myPolygonMeshing->setCoordinates(SpaceDimension, NumberOfNodes, PolygonCoordinates,
01026                                        "CARTESIAN", MED_FULL_INTERLACE);
01027     }
01028     catch (const std::exception &e)
01029     {
01030       CPPUNIT_FAIL(e.what());
01031     }
01032     catch (...)
01033     {
01034       CPPUNIT_FAIL("Unknown exception");
01035     }
01036 
01037     CPPUNIT_ASSERT_EQUAL(SpaceDimension, myPolygonMeshing->getSpaceDimension());
01038     CPPUNIT_ASSERT_EQUAL(NumberOfNodes, myPolygonMeshing->getNumberOfNodes());
01039 
01040     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setTypes(PolygonTypes, MED_CELL));
01041     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setNumberOfElements(PolygonNumberOfElements, MED_CELL));
01042     CPPUNIT_ASSERT_EQUAL(2, myPolygonMeshing->getMeshDimension());
01043 
01044     try
01045     {
01046       myPolygonMeshing->setCoordinatesNames(Names);
01047     }
01048     catch (const std::exception &e)
01049     {
01050       CPPUNIT_FAIL(e.what());
01051     }
01052     catch (...)
01053     {
01054       CPPUNIT_FAIL("Unknown exception");
01055     }
01056 
01057     try
01058     {
01059       myPolygonMeshing->setCoordinatesUnits(Units);
01060     }
01061     catch (const std::exception &e)
01062     {
01063       CPPUNIT_FAIL(e.what());
01064     }
01065     catch (...)
01066     {
01067       CPPUNIT_FAIL("Unknown exception");
01068     }
01069 
01070     const int sizeTri = 3*2;
01071     int ConnectivityTri[sizeTri]=
01072       {
01073         1, 7, 2, 3, 9, 4
01074       };
01075 
01076     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setConnectivity(MED_CELL, MED_TRIA3,ConnectivityTri));
01077     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setConnectivity(MED_CELL, MED_POLYGON,REFpolygonFaces, REFpolygonIndex));
01078 
01079     bool PolygonConn = false;
01080     CPPUNIT_ASSERT_NO_THROW(PolygonConn = myPolygonMeshing->existConnectivity(MED_NODAL, MED_CELL));
01081     if(PolygonConn)
01082     {
01083       int Polytypes;
01084       CPPUNIT_ASSERT_NO_THROW(Polytypes = myPolygonMeshing->getNumberOfTypes(MED_CELL));
01085       CPPUNIT_ASSERT_EQUAL(NbOfTypes,Polytypes);
01086 
01087       const MED_EN::medGeometryElement * PolyTypes;
01088       CPPUNIT_ASSERT_NO_THROW(PolyTypes = myPolygonMeshing->getTypes(MED_CELL));
01089       CPPUNIT_ASSERT_EQUAL(PolyTypes[NbOfTypes-1],MED_POLYGON);
01090 
01091       for(int t = 0; t < Polytypes; t++)
01092       {
01093         CPPUNIT_ASSERT_NO_THROW( myPolygonMeshing->getNumberOfElements(MED_CELL, PolyTypes[t]));
01094       }
01095       medGeometryElement geomPolyElem;
01096       CPPUNIT_ASSERT_NO_THROW(geomPolyElem = myPolygonMeshing->getElementType(MED_CELL, 1));
01097       CPPUNIT_ASSERT_EQUAL(MED_TRIA3,geomPolyElem);
01098 
01099       CPPUNIT_ASSERT_EQUAL(myPolygonMeshing->getNumberOfElements(MED_CELL,MED_POLYGON),nbOfPolygons);
01100       CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->getConnectivityLength(MED_NODAL,MED_CELL,MED_POLYGON));
01101       myPolygonMeshing->removeReference();
01102       const int * PolygonConn;
01103       CPPUNIT_ASSERT_THROW(PolygonConn = myMeshingPoly->getConnectivity(MED_NODAL,MED_CELL,MED_POLYGON),MEDEXCEPTION);
01104     }
01105     myMeshingPoly->removeReference();
01107     // TEST : SUPPORT* sup = myMeshPointe->getSupportOnAll()) //
01109     {
01110       MESH * myMeshPointe = new MESH(MED_DRIVER, filename, meshname);
01111       const SUPPORT* sup = myMeshPointe->getSupportOnAll(MED_CELL);
01112       CPPUNIT_ASSERT( sup->isOnAllElements() );
01113       CPPUNIT_ASSERT_EQUAL( myMeshPointe->getNumberOfTypes( sup->getEntity() ),
01114                             sup->getNumberOfTypes());
01115       CPPUNIT_ASSERT( sup->getNumber( MED_ALL_ELEMENTS ));
01116       myMeshPointe->removeReference();
01117     }
01118 
01120     // TEST 3: test MESH on  MEDMEMTest::createTestMesh()//
01122 
01123     MESH* myMesh3 = MEDMEMTest_createTestMesh();
01124 
01125     int MeshDim  = myMesh3->getMeshDimension();
01126     medEntityMesh constituentEntity;
01127     if (MeshDim==3) 
01128     {
01129       constituentEntity = MED_CELL;
01130     }
01131     if (MeshDim==2) 
01132     {
01133       constituentEntity = MED_FACE;
01134     }
01135     if (MeshDim==1) 
01136     {
01137       constituentEntity = MED_EDGE;
01138     }
01139 
01140     int SpaceDim = myMesh3->getSpaceDimension();
01141 
01142     // Show Reverse Nodal Connectivity
01143     const int* ReverseNodalConnectivity;
01144     const int* ReverseNodalConnectivityIndex;
01145     int ReverseLength;
01146     int ReverseIdxLength;
01147 
01148     CONNECTIVITY* myConnectivity3 = (CONNECTIVITY*)myMesh3->getConnectivityptr();
01149     myConnectivity3->setNumberOfNodes(myMesh3->getNumberOfNodes());
01150 
01151     CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivity= myMesh3->getReverseConnectivity(MED_NODAL, MED_CELL));
01152     CPPUNIT_ASSERT_NO_THROW(ReverseLength = myMesh3->getReverseConnectivityLength(MED_NODAL, MED_CELL));
01153     CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivityIndex = myMesh3->getReverseConnectivityIndex(MED_NODAL, MED_CELL));
01154     CPPUNIT_ASSERT_NO_THROW(ReverseIdxLength = myMesh3->getReverseConnectivityIndexLength(MED_NODAL, MED_CELL));
01155     CPPUNIT_ASSERT(ReverseIdxLength == myMesh3->getNumberOfNodes()+1);
01156 
01157     ostringstream out;
01158     for (int i = 0; i < myMesh3->getNumberOfNodes(); i++) 
01159     {
01160       out << "Node "<< i+1 << " : ";
01161       for (int j = ReverseNodalConnectivityIndex[i]; j < ReverseNodalConnectivityIndex[i+1]; j++)
01162         out << ReverseNodalConnectivity[j-1] << " ";
01163       out << endl;
01164     }
01165 
01166     // Show Descending Connectivity
01167     int NumberOfElements1;
01168     const int * connectivity;
01169     const int * connectivity_index;
01170     myMesh3->calculateConnectivity(MED_DESCENDING, MED_EN::MED_CELL);
01171     try
01172     {
01173       NumberOfElements1 = myMesh3->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
01174       connectivity = myMesh3->getConnectivity( MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS);
01175       connectivity_index = myMesh3->getConnectivityIndex(MED_DESCENDING, MED_CELL);
01176     }
01177     catch (MEDEXCEPTION m) 
01178     {
01179       CPPUNIT_FAIL(m.what());
01180     }
01181 
01182     for (int j = 0; j < NumberOfElements1; j++) 
01183     {
01184       out << "Element " << j+1 << " : ";
01185       for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
01186         out << connectivity[k-1] << " ";
01187       out << endl;
01188     }
01189 
01190     //test 3D mesh
01191     for(int ind = SpaceDim; ind > 1; ind-- )
01192     {
01193       int NumberOfElem = myMesh3->getNumberOfElements (constituentEntity,MED_ALL_ELEMENTS);
01194       if(NumberOfElem < 1) continue;
01195 
01196       const SUPPORT * sup = myMesh3->getSupportOnAll( constituentEntity );
01197 
01198       if (ind == 2)
01199       {
01200         // test of normal(for 1d or 2d elements)
01201         FIELD<double>* normal;
01202         CPPUNIT_ASSERT_NO_THROW(normal = myMesh3->getNormal(sup));
01203 
01204         const int nbNormVals = 47*3;
01205         double refNormals[nbNormVals] = 
01206           {
01207             0, 1, 0    ,// #1 
01208             -1, 0, 0   ,// #2 
01209             -0, 0, 2   ,// #3 
01210             1, -1, -2  ,// #4 
01211             -1, 0, 0   ,// #5 
01212             0, 0, 2    ,// #6 
01213             1, 1, -2   ,// #7 
01214             0, -1, 0   ,// #8 
01215             0, -0, 2   ,// #9 
01216             -1, 1, -2  , // #10
01217             0, 0, 2    , // #11
01218             -1, -1, -2 , // #12
01219             -1, 0, 1   , // #13
01220             0, -1, 1   , // #14
01221             1, 1, 0    , // #15
01222             0, -1, 1   , // #16
01223             1, -0, 1   , // #17
01224             -1, 1, 0   , // #18
01225             1, 0, 1    , // #19
01226             0, 1, 1    , // #20
01227             -1, -1, 0  , // #21
01228             -0, 1, 1   , // #22
01229             -1, 0, 1   , // #23
01230             1, -1, 0   , // #24
01231             -1, 0, 1   , // #25
01232             1, 0, 1    , // #26
01233             0, -1, 1   , // #27
01234             -0, 1, 1   , // #28
01235             1, 0, 1    , // #29
01236             -1, 0, 1   , // #30
01237             0, 1, 1    , // #31
01238             0, -1, 1   , // #32
01239             1, 0, 1    , // #33
01240             0, -1, 1   , // #34
01241             -1, 0, 1   , // #35
01242             -0, 1, 1   , // #36
01243             0, 0, 4    , // #37
01244             -0, -0, -4 , // #38
01245             0, 0, 4    , // #39
01246             0, 2, 0    , // #40
01247             -2, -0, -0 , // #41
01248             0, -2, 0   , // #42
01249             2, -0, 0   , // #43
01250             0, 2, 0    , // #44
01251             -2, -0, -0 , // #45
01252             0, -2, 0   , // #46
01253             2, -0, 0     // #47
01254           };
01255         for ( int i = 0; i < nbNormVals; ++i )
01256           CPPUNIT_ASSERT_DOUBLES_EQUAL( refNormals[i], normal->getValue()[i], 1e-6);
01257         //       double normal_square, norm;
01258         //       double maxnorm=0.;
01259         //       double minnorm=0.;
01260         //       double tmp_value;
01261         //       for (int i = 1; i<=NumberOfElem; i++) {
01262         //         normal_square = 0.;
01263         //         cout << "Normal " << i << " ";
01264         //         for (int j=1; j<=SpaceDim; j++) {
01265         //           tmp_value = normal->getValueIJ(i,j);
01266         //           normal_square += tmp_value*tmp_value;
01267         //           cout << tmp_value << " ";
01268         //         }
01269         //         norm = sqrt(normal_square);
01270         //         maxnorm = dmax(maxnorm,norm);
01271         //         minnorm = dmin(minnorm,norm);
01272         //         cout << ", Norm = " << norm << endl;
01273         //       }
01274         //       cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
01275         delete normal;
01276 
01277         // test of area(for 2d elements)
01278         FIELD<double>* area;
01279         CPPUNIT_ASSERT_NO_THROW(area = myMesh3->getArea(sup));
01280 
01281         //       double maxarea,minarea,areatot;
01282         //       maxarea = 0.;
01283         //       minarea = 0.;
01284         //       areatot = 0.0;
01285         //       for (int i = 1; i<=NumberOfElem;i++)
01286         //       {
01287         //         cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
01288         //         maxarea = dmax(maxarea,area->getValueIJ(i,1));
01289         //         minarea = dmin(minarea,area->getValueIJ(i,1));
01290         //         areatot = areatot + area->getValueIJ(i,1);
01291         //       }
01292 
01293         //       cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
01294         //       cout << "Support Area " << areatot << endl;
01295         const int nbAreas = 47;
01296         double refArea [nbAreas] = 
01297           {
01298             1, 1, 2, 2.44949, 1, 2, 2.44949, 1, 2,  2.44949,  2,  2.44949,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  1.41421,  4,  4,  4,  2,  2,  2,  2,  2,  2,  2,  2 
01299           };
01300         for ( int i = 0; i < nbAreas; ++i )
01301           CPPUNIT_ASSERT_DOUBLES_EQUAL( refArea[i], area->getValue()[i], 1e-5);
01302         delete area;
01303       }
01304 
01305       // test of barycenter(for 3d and 2d elements)
01306       FIELD<double>* barycenter;
01307       CPPUNIT_ASSERT_NO_THROW(barycenter = myMesh3->getBarycenter(sup));
01308 
01309       CPPUNIT_ASSERT_NO_THROW(NumberOfElem = myMesh3->getNumberOfElements(constituentEntity,MED_ALL_ELEMENTS));
01310       if ( ind == 3 )
01311       {
01312         double refBC[16*3] = 
01313           {
01314             0.5, -0.5, 0.75    ,// #1 
01315             0.5, 0.5, 0.75     ,// #2 
01316             -0.5, 0.5, 0.75    ,// #3 
01317             -0.5, -0.5, 0.75   ,// #4 
01318             0.75, 0.75, 1.25   ,// #5 
01319             -0.75, 0.75, 1.25  ,// #6 
01320             -0.75, -0.75, 1.25 ,// #7 
01321             0.75, -0.75, 1.25  ,// #8 
01322             1, 0, 1.5          ,// #9 
01323             0, 1, 1.5          , // #10
01324             -1, 0, 1.5         , // #11
01325             0, -1, 1.5         , // #12
01326             0, 0, 1.8          , // #13
01327             0, 0, 4.2          , // #14
01328             0, 0, 2.5          , // #15
01329             0, 0, 3.5          , // #16
01330           };
01331         for ( int i = 0; i < 16*3; ++i )
01332           CPPUNIT_ASSERT_DOUBLES_EQUAL( refBC[i], barycenter->getValue()[i], 1e-6);
01333       }
01334       if ( ind == 2 )
01335       {
01336         double refBC[47*3] = 
01337           {
01338             0.666667, 0, 0.666667          ,// #1 
01339             0, -0.666667, 0.666667         ,// #2 
01340             0.666667, -0.666667, 1         ,// #3 
01341             0.666667, -0.666667, 0.666667  ,// #4 
01342             0, 0.666667, 0.666667          ,// #5 
01343             0.666667, 0.666667, 1          ,// #6 
01344             0.666667, 0.666667, 0.666667   ,// #7 
01345             -0.666667, 0, 0.666667         ,// #8 
01346             -0.666667, 0.666667, 1         ,// #9 
01347             -0.666667, 0.666667, 0.666667 , // #10
01348             -0.666667, -0.666667, 1       , // #11
01349             -0.666667, -0.666667, 0.666667, // #12
01350             0.333333, 1, 1.33333          , // #13
01351             1, 0.333333, 1.33333          , // #14
01352             1, 1, 1.33333                 , // #15
01353             -1, 0.333333, 1.33333         , // #16
01354             -0.333333, 1, 1.33333         , // #17
01355             -1, 1, 1.33333                , // #18
01356             -0.333333, -1, 1.33333        , // #19
01357             -1, -0.333333, 1.33333        , // #20
01358             -1, -1, 1.33333               , // #21
01359             1, -0.333333, 1.33333         , // #22
01360             0.333333, -1, 1.33333         , // #23
01361             1, -1, 1.33333                , // #24
01362             0.666667, 0, 1.66667          , // #25
01363             1.33333, 0, 1.66667           , // #26
01364             0, 0.666667, 1.66667          , // #27
01365             0, 1.33333, 1.66667           , // #28
01366             -0.666667, 0, 1.66667         , // #29
01367             -1.33333, 0, 1.66667          , // #30
01368             0, -0.666667, 1.66667         , // #31
01369             0, -1.33333, 1.66667          , // #32
01370             0.666667, 0, 4.33333          , // #33
01371             0, -0.666667, 4.33333         , // #34
01372             -0.666667, 0, 4.33333         , // #35
01373             0, 0.666667, 4.33333          , // #36
01374             0, 0, 2                       , // #37
01375             0, 0, 4                       , // #38
01376             0, 0, 3                       , // #39
01377             0, 1, 2.5                     , // #40
01378             -1, 0, 2.5                    , // #41
01379             0, -1, 2.5                    , // #42
01380             1, 0, 2.5                     , // #43
01381             0, 1, 3.5                     , // #44
01382             -1, 0, 3.5                    , // #45
01383             0, -1, 3.5                    , // #46
01384             1, 0, 3.5                     , // #47
01385           };
01386         for ( int i = 0; i < 47*3; ++i )
01387         {
01388           CPPUNIT_ASSERT_DOUBLES_EQUAL( refBC[i], barycenter->getValue()[i], 1e-5);
01389         }
01390       }
01391       //     for (int i = 1; i<=NumberOfElem;i++)
01392       //     {
01393       //     cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
01394       //     }
01395 
01396       delete barycenter;
01397 
01398       // test of volume(for 3d elements)
01399       if (ind == 3)
01400       {
01401         FIELD<double>* volume;
01402         CPPUNIT_ASSERT_NO_THROW(volume= myMesh3->getVolume(sup));
01403 
01404         //       double maxvol,minvol,voltot;
01405         //       maxvol = 0.;
01406         //       minvol = 0.;
01407         //       voltot = 0.0;
01408         //       for (int i = 1; i<=NumberOfElem;i++)
01409         //       {
01410         //         cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
01411         //         maxvol = dmax(maxvol,volume->getValueIJ(i,1));
01412         //         minvol = dmin(minvol,volume->getValueIJ(i,1));
01413         //         voltot = voltot + volume->getValueIJ(i,1);
01414         //       }
01415 
01416         //       cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
01417         //       cout << "Support Volume " << voltot << endl;
01418 
01419         double refVol[16] = 
01420           {
01421             0.666667, 0.666667, 0.666667, 0.666667, 0.666667, 0.666667, 0.666667, 0.666667, 0.666667,  0.666667,  0.666667,  0.666667,  1.333333,  1.333333,  4,  4 
01422           };
01423         for ( int i = 0; i < 16; ++i )
01424           CPPUNIT_ASSERT_DOUBLES_EQUAL( refVol[i], volume->getValue()[i], 1e-6);
01425 
01426         delete volume;
01427 
01428         // test of skin
01429         SUPPORT *skin;
01430         CPPUNIT_ASSERT_NO_THROW(skin = myMesh3->getSkin(sup));
01431 
01432         //test mergeSupports and intersectSupports. vactor contains only 1 elements
01433         vector<SUPPORT *> myVectSup;
01434         myVectSup.push_back(skin);
01435 
01436         //method return a copy of skin object
01437         SUPPORT *copyMergeSkin;
01438         CPPUNIT_ASSERT_NO_THROW(copyMergeSkin = myMesh3->mergeSupports(myVectSup));
01439         try
01440         {
01441           CPPUNIT_ASSERT(copyMergeSkin->deepCompare(*skin));
01442         }
01443         catch (const std::exception &e)
01444         {
01445           CPPUNIT_FAIL(e.what());
01446         }
01447         catch (...)
01448         {
01449           CPPUNIT_FAIL("Unknown exception");
01450         }
01451 
01452         //method return a copy of skin object
01453         SUPPORT *copyIntersectSkin;
01454         CPPUNIT_ASSERT_NO_THROW(copyIntersectSkin = myMesh3->intersectSupports(myVectSup));
01455         try
01456         {
01457           CPPUNIT_ASSERT(copyIntersectSkin->deepCompare(*skin));
01458         }
01459         catch (const std::exception &e)
01460         {
01461           CPPUNIT_FAIL(e.what());
01462         }
01463         catch (...)
01464         {
01465           CPPUNIT_FAIL("Unknown exception");
01466         }
01467 
01468         skin->removeReference();
01469         copyMergeSkin->removeReference();
01470         copyIntersectSkin->removeReference();
01471       }
01472       constituentEntity++;
01473     }
01474 
01475     // 0020911: [CEA 413] getBarycenter on polygons
01476     {
01477       MESHING* polygonMesh = new MESHING();
01478       MEDMEM::AutoDeref derefMesh( polygonMesh );
01479 
01480       polygonMesh->setName("polygonMesh");
01481 
01482       const int spaceDim = 2, nbNodes = 6;
01483       const double coords[spaceDim*nbNodes] = 
01484         {
01485           0,0, 1,0, 2,0,  0,1, 1,1, 2,1 
01486         };
01487       polygonMesh->setCoordinates( spaceDim, nbNodes, coords, "CART", MED_EN::MED_FULL_INTERLACE);
01488 
01489       const MED_EN::medGeometryElement type = MED_EN::MED_POLYGON;
01490       polygonMesh->setNumberOfTypes(1, MED_EN::MED_CELL);
01491       polygonMesh->setTypes( &type, MED_EN::MED_CELL );
01492 
01493       const int nbPolygons = 2;
01494       const int conn[nbPolygons*4] = 
01495         {
01496           1,2,5,4, 2,3,6,5 
01497         };
01498       const int index[nbPolygons+1] = 
01499         {
01500           1,5,9
01501         };
01502       polygonMesh->setNumberOfElements( &nbPolygons, MED_EN::MED_CELL );
01503       polygonMesh->setConnectivity( MED_EN::MED_CELL, type, conn, index );
01504 
01505       FIELD<double>* barycenter;
01506       const SUPPORT* sup = polygonMesh->getSupportOnAll(MED_CELL);
01507       CPPUNIT_ASSERT_NO_THROW(barycenter = polygonMesh->getBarycenter(sup));
01508       CPPUNIT_ASSERT_EQUAL( 2, barycenter->getNumberOfValues() );
01509       CPPUNIT_ASSERT_EQUAL( 2, barycenter->getNumberOfComponents() );
01510       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, barycenter->getValueIJ(1,1), 1e-10);
01511       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, barycenter->getValueIJ(1,2), 1e-10);
01512       CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.5, barycenter->getValueIJ(2,1), 1e-10);
01513       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, barycenter->getValueIJ(2,2), 1e-10);
01514       barycenter->removeReference();
01515     }
01516 
01517 
01518     // Testing length and normal vectors on 1d elements
01519     {
01520       // coordinates
01521       int NumberOfNodes3 = 4;
01522 
01523       string Names3[3] = 
01524         {
01525           "X","Y","Z"
01526         };
01527       string Units3[3] = 
01528         {
01529           "cm","cm","cm"
01530         };
01531 
01532       double Coordinates3[4*2] = 
01533         {
01534           0.0,  0.0,  // n1
01535           1.0,  1.0,  // n2
01536           0.0,  1.0,  // n3
01537           1.0,  0.0
01538         }; // n4
01539 
01540       const int NumberOfEdgeTypes = 1;
01541       MED_EN::medGeometryElement EdgeTypes[NumberOfEdgeTypes] = 
01542         {
01543           MED_SEG2
01544         };
01545       const int NumberOfEdges[NumberOfEdgeTypes] = 
01546         {
01547           4
01548         };
01549       int ConnectivityEdge[4*2] = 
01550         {
01551           1,2, 2,3, 3,4, 4,1
01552         };
01553 
01554       // CREATE THE MESH
01555       MEDMEM::MESHING* myMeshing3 = new MEDMEM::MESHING;
01556       myMeshing3->setName("meshing3");
01557       myMeshing3->setCoordinates(/*SpaceDimension*/2, NumberOfNodes3, Coordinates3,
01558                                  "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
01559       myMeshing3->setCoordinatesNames(Names3);
01560       myMeshing3->setCoordinatesUnits(Units3);
01561 
01562       // define connectivities
01563       //      cell part
01564       const int NumberOfTypes3 = 1;
01565       medGeometryElement Types3[NumberOfTypes3] = 
01566         {
01567           MED_QUAD4
01568         };
01569       const int NumberOfElements3[NumberOfTypes3] = 
01570         {
01571           1
01572         };
01573 
01574       myMeshing3->setNumberOfTypes(NumberOfTypes3,MED_CELL);
01575       myMeshing3->setTypes(Types3,MED_CELL);
01576       myMeshing3->setNumberOfElements(NumberOfElements3,MED_CELL);
01577 
01578       int Connectivityquad[1*4] = 
01579         {
01580           1,2,3,4
01581         };
01582 
01583       myMeshing3->setConnectivity(MED_CELL,MED_QUAD4,Connectivityquad);
01584 
01585       myMeshing3->setNumberOfTypes(NumberOfEdgeTypes, MED_EDGE);
01586       myMeshing3->setTypes(EdgeTypes, MED_EDGE);
01587       myMeshing3->setNumberOfElements(NumberOfEdges, MED_EDGE);
01588 
01589       myMeshing3->setConnectivity( MED_EDGE, MED_SEG2,ConnectivityEdge);
01590 
01591       //test 2D mesh
01592       //int NumberOfElem = myMeshing3->getNumberOfElements (MED_EDGE, MED_ALL_ELEMENTS);
01593 
01594       const SUPPORT * sup = myMeshing3->getSupportOnAll( MED_EDGE );
01595 
01596       // test of normal(for 1d or 2d elements)
01597       FIELD<double>* normal;
01598       CPPUNIT_ASSERT_NO_THROW(normal = myMeshing3->getNormal(sup));
01599 
01600       //     double normal_square, norm;
01601       //     double maxnorm=0.;
01602       //     double minnorm=0.;
01603       //     double tmp_value;
01604       //     for (int i = 1; i<=NumberOfElem; i++) {
01605       //       normal_square = 0.;
01606       //       cout << "Normal " << i << " ";
01607       //       for (int j=1; j<=/*SpaceDimension*/2; j++) {
01608       //         tmp_value = normal->getValueIJ(i,j);
01609       //         normal_square += tmp_value*tmp_value;
01610       //         cout << tmp_value << " ";
01611       //       }
01612       //       norm = sqrt(normal_square);
01613       //       maxnorm = dmax(maxnorm,norm);
01614       //       minnorm = dmin(minnorm,norm);
01615       //       cout << ", Norm = " << norm << endl;
01616       //     }
01617       //     cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
01618       {
01619         double refNormals[8] = 
01620           {
01621             -1, 1 ,
01622             -0, -1,
01623             1, 1 ,
01624             -0, -1
01625           };
01626         for ( int i = 0; i < 8; ++i )
01627           CPPUNIT_ASSERT_DOUBLES_EQUAL( refNormals[i], normal->getValue()[i], 1e-6);
01628       }
01629       // test of length(for 1d elements)
01630       FIELD<double>* length;
01631       CPPUNIT_ASSERT_NO_THROW(length = myMeshing3->getLength(sup));
01632 
01633       //     double length_value,maxlength,minlength;
01634       //     maxlength = 0;
01635       //     minlength = 0;
01636       //     for (int i = 1; i<=NumberOfElem;i++) {
01637       //       length_value = length->getValueIJ(i,1);
01638       //       cout << "Length " << i << " " << length_value << endl;
01639       //       maxlength = dmax(maxlength,length_value);
01640       //       minlength = dmin(minlength,length_value);
01641       //     }
01642       //     cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
01643       double refLen[4] = 
01644         {
01645           1.41421,1,1.41421,1
01646         };
01647       for ( int i = 0; i < 4; ++i )
01648         CPPUNIT_ASSERT_DOUBLES_EQUAL( refLen[i], length->getValue()[i], 1e-5);
01649 
01650       vector< FIELD<double> *> myVectField1;
01651       myVectField1.push_back(normal);
01652       myVectField1.push_back(length);
01653       CPPUNIT_ASSERT_NO_THROW(myMeshing3->mergeFields(myVectField1)->removeReference());
01654 
01655       normal->removeReference();
01656       length->removeReference();
01657       {
01658         vector<SUPPORT *> myVectSupEmpty;
01659         CPPUNIT_ASSERT_THROW(myMesh3->mergeSupports(myVectSupEmpty), MEDEXCEPTION);
01660       }
01661 
01662       // test mergeFields method: Fields have the same value type
01663       //intersectSupports and mergeSupports methods
01664       {
01665         SUPPORT * sup1 = new SUPPORT;
01666         sup1->setMesh( myMeshing3 );
01667         sup1->setEntity( MED_EDGE );
01668         SUPPORT * sup2 = new SUPPORT;
01669         sup2->setMesh( myMeshing3 );
01670         sup2->setEntity( MED_EDGE );
01671         MED_EN::medGeometryElement gtEdges[1] = 
01672           {
01673             MED_SEG2
01674           };
01675         int nbEdges1[1] = 
01676           {
01677             1
01678           };
01679         int edges1[1] = 
01680           {
01681             1
01682           };
01683         int nbEdges2[1] = 
01684           {
01685             2
01686           };
01687         int edges2[2] = 
01688           {
01689             2,3
01690           };
01691         sup1->setpartial("description 1", 1, 1, gtEdges, nbEdges1, edges1);
01692         sup2->setpartial("description 1", 1, 2, gtEdges, nbEdges2, edges2);
01693 
01694         vector<SUPPORT *> myVectSup3;
01695         myVectSup3.push_back(sup1);
01696         myVectSup3.push_back(sup2);
01697         //method return a MergeSup on the union of all SUPPORTs in Supports.
01698         SUPPORT *MergeSup;
01699         CPPUNIT_ASSERT_NO_THROW(MergeSup = myMesh3->mergeSupports(myVectSup3));
01700         {
01701           ostringstream out;
01702           out << *MergeSup << endl;
01703         }
01704         MergeSup->removeReference();
01705 
01706         //method return a intersection of all SUPPORTs in IntersectSup
01707         SUPPORT *IntersectSup;
01708         CPPUNIT_ASSERT_NO_THROW(IntersectSup = myMesh3->intersectSupports(myVectSup3));
01709         {
01710           ostringstream out;
01711           if (IntersectSup != NULL) out<< *IntersectSup <<endl;
01712         }
01713         IntersectSup->removeReference();
01714 
01715         FIELD<double> * length1 = myMeshing3->getLength(sup1);
01716         FIELD<double> * length2 = myMeshing3->getLength(sup2);
01717 
01718         vector< FIELD<double> *> myVect12;
01719         myVect12.push_back(length1);
01720         myVect12.push_back(length2);
01721 
01722         FIELD<double> * length12;
01723         CPPUNIT_ASSERT_NO_THROW(length12 = myMeshing3->mergeFields(myVect12));
01724         length12->removeReference();
01725 
01726         sup1->removeReference();
01727         sup2->removeReference();
01728         length1->removeReference();
01729         length2->removeReference();
01730       }
01731       myMeshing3->removeReference();
01732     }
01733     myMesh3->removeReference();
01734 
01736     // TEST 4: test MESH constructed from file pointe.med //
01738     MESH * myMesh4 = new MESH();
01739     myMesh4->setName(meshname);
01740     MED_MESH_RDONLY_DRIVER myMeshDriver (filename, myMesh4);
01741     myMeshDriver.setMeshName(meshname);
01742 
01743     //Mesh has no driver->segmentation violation
01744     CPPUNIT_ASSERT_THROW(myMesh4->read(), MEDEXCEPTION);
01745 
01746     //Add an existing MESH driver.
01747     int myDriver4;
01748     CPPUNIT_ASSERT_NO_THROW(myDriver4 = myMesh4->addDriver(myMeshDriver));
01749 
01750     //read all objects in the file
01751     CPPUNIT_ASSERT_NO_THROW(myMesh4->read(myDriver4));
01752 
01753     if (myMesh4->getIsAGrid()) 
01754     {
01755       GRID* myGrid = dynamic_cast<GRID*>(myMesh4);
01756       CPPUNIT_ASSERT(myGrid);
01757     }
01758 
01759     //myDriver4->DRONLY->can't write
01760     CPPUNIT_ASSERT_THROW(myMesh4->write(myDriver4), MEDEXCEPTION);
01761 
01762     // add new driver
01763     int idMeshV21;
01764     CPPUNIT_ASSERT_NO_THROW(idMeshV21 = myMesh4->addDriver(MED_DRIVER,filenameout21));
01765 
01766     //Write all the content of the MESH using driver referenced by the integer handler index.
01767     CPPUNIT_ASSERT_NO_THROW(myMesh4->write(idMeshV21));
01768 
01769     // remove driver from mesh
01770     CPPUNIT_ASSERT_NO_THROW(myMesh4->rmDriver(myDriver4));
01771 
01772     // ensure exception is raised on second attempt to remove driver
01773     CPPUNIT_ASSERT_THROW(myMesh4->rmDriver(myDriver4),MEDEXCEPTION);
01774 
01775     // Create a MESH object using a MESH driver of type MED_DRIVER associated with file fileName.
01776     MESH* myMesh5;
01777     CPPUNIT_ASSERT_NO_THROW(myMesh5 = new MESH(MED_DRIVER, filename, meshname));
01778     if(myMesh5->getIsAGrid())
01779     {
01780       GRID* myGrid = dynamic_cast<GRID*>(myMesh4);
01781       CPPUNIT_ASSERT(myGrid);
01782     }
01783 
01784     //ensure two meshes constructed from one file in two different ways are equal
01785     CPPUNIT_ASSERT(myMesh5->deepCompare(*myMesh4));
01786 
01787     // test other variants of read() and write()
01788     {
01789       const string otherName1("otherName1"), otherName2("otherName2");
01790       MESH mesh1, mesh2, mesh3;
01791 
01792       //myMesh5 -> filenameout21
01793       // GMESH::write(driverTypes driverType, const string& filename,const string& meshname)
01794       myMesh5->write( MED_DRIVER, filenameout21, otherName1);
01795       CPPUNIT_ASSERT_THROW( myMesh5->write(myMeshDriver), MEDEXCEPTION); // write with RDONLY driver
01796 
01797       //filenameout21 -> mesh1
01798       // GMESH::read(driverTypes driverType, const string& filename, const string& meshname);
01799       CPPUNIT_ASSERT_THROW( mesh1.read(VTK_DRIVER,filenameout21,otherName1), MEDEXCEPTION);
01800       CPPUNIT_ASSERT_THROW( mesh1.read(MED_DRIVER, filenameout21,otherName2), MEDEXCEPTION);
01801       mesh1.read(MED_DRIVER,filenameout21,otherName1);
01802       CPPUNIT_ASSERT(myMesh5->deepCompare(mesh1));
01803 
01804       MED_MESH_RDWR_DRIVER driver;
01805       driver.setFileName( filenameout21 );
01806       driver.setMeshName( otherName2 );
01807 
01808       // mesh1 -> filenameout21
01809       // GMESH::write( const GENDRIVER& )
01810       mesh1.write( driver );
01811 
01812       // filenameout21 -> mesh2
01813       // GMESH::read( const GENDRIVER& )
01814       mesh2.read( driver );
01815       CPPUNIT_ASSERT(myMesh5->deepCompare(mesh2));
01816       // check that GMESH::write() clears the file by default
01817       CPPUNIT_ASSERT_THROW( mesh3.read( MED_DRIVER, filenameout21,otherName1), MEDEXCEPTION);
01818       // write by adding a mesh
01819       myMesh5->write( MED_DRIVER, filenameout21, otherName1, MED_EN::RDWR );
01820       CPPUNIT_ASSERT_NO_THROW( mesh3.read( MED_DRIVER, filenameout21,otherName1));
01821     }
01822 
01823     int myDriver6;
01824     MESH* myMesh6 = new MESH();
01825     myDriver6 = myMesh6->addDriver(MED_DRIVER, filename, meshname, RDONLY);
01826 
01827     myMesh6->read(myDriver6);
01828 
01829     //ensure two meshes constracted from one file in two different ways are equal
01830     CPPUNIT_ASSERT(myMesh6->deepCompare(*myMesh4));
01831 
01832     //test FAMILY
01833     int NumberOfFamilies4;
01834     CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies4 = myMesh6->getNumberOfFamilies(MED_CELL));
01835     CPPUNIT_ASSERT_MESSAGE("Current mesh hasn't Families", NumberOfFamilies4 != 0);
01836 
01837     vector<FAMILY*> families4;
01838     CPPUNIT_ASSERT_NO_THROW(families4 = myMesh6->getFamilies(MED_CELL));
01839     CPPUNIT_ASSERT((int)families4.size() == NumberOfFamilies4);
01840     for(int nb = 1; nb <= NumberOfFamilies4; nb++ )
01841     {
01842       const FAMILY* family;
01843       CPPUNIT_ASSERT_NO_THROW(family = myMesh6->getFamily(MED_CELL, nb));
01844       CPPUNIT_ASSERT_EQUAL(family->getName(), families4[nb-1]->getName());
01845     }
01846 
01847     //get support which reference all elements on the boundary of mesh.
01848     SUPPORT*myBndSup;
01849     CPPUNIT_ASSERT_THROW(myMesh6->getBoundaryElements(MED_CELL), MEDEXCEPTION);
01850     //get only face in 3D.
01851     CPPUNIT_ASSERT_NO_THROW(myBndSup = myMesh6->getBoundaryElements(MED_FACE));
01852 
01853     //test buildSupportOnElementsFromElementList and buildSupportOnNodeFromElementList
01854     const int * myConnectivityValue6;
01855     CPPUNIT_ASSERT_NO_THROW(myConnectivityValue6 = myMesh6->getReverseConnectivity(MED_DESCENDING));
01856     const int * myConnectivityIndex6;
01857     CPPUNIT_ASSERT_NO_THROW(myConnectivityIndex6 = myMesh6->getReverseConnectivityIndex(MED_DESCENDING));
01858     int numberOfElem6;
01859     CPPUNIT_ASSERT_NO_THROW(numberOfElem6 = myMesh6->getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS));
01860     list<int> myElementsList6;
01861 
01862     for (int i=0; i<numberOfElem6; i++)
01863       if (myConnectivityValue6[myConnectivityIndex6[i]] == 0) 
01864       {
01865         myElementsList6.push_back(i+1);
01866       }
01867 
01868     SUPPORT * mySupportOnNode;
01869     SUPPORT * mySupportOnElem;
01870     CPPUNIT_ASSERT_NO_THROW(mySupportOnElem = myMesh6->buildSupportOnElementsFromElementList(myElementsList6,MED_FACE));
01871     CPPUNIT_ASSERT(mySupportOnElem->deepCompare(*myBndSup));
01872     myBndSup->removeReference();
01873     CPPUNIT_ASSERT_EQUAL(MED_FACE, mySupportOnElem->getEntity());
01874 
01875     list<int>::const_iterator iteronelem = myElementsList6.begin();
01876     for (int i = 1; i <= 3; i++, iteronelem++) 
01877     {
01878       CPPUNIT_ASSERT_EQUAL(i, mySupportOnElem->getValIndFromGlobalNumber(*iteronelem));
01879     }
01880     mySupportOnElem->removeReference();
01881     CPPUNIT_ASSERT_NO_THROW(mySupportOnNode = myMesh6->buildSupportOnNodeFromElementList(myElementsList6,MED_FACE));
01882     SUPPORT *suppp=myMesh6->getBoundaryElements(MED_NODE);
01883     CPPUNIT_ASSERT(mySupportOnNode->deepCompare( *(suppp)));
01884     suppp->removeReference();
01885     mySupportOnNode->removeReference();
01886 
01887     //sets mesh fields to initial values
01888     myMesh6->init();
01889 
01890     //ensure two meshes constracted from one file in two different ways are equal
01891     CPPUNIT_ASSERT(!myMesh6->deepCompare(*myMesh4));
01892 
01893     //ensure mesh is empty
01894     CPPUNIT_ASSERT(myMesh6->getSpaceDimension() == MED_INVALID);
01895     CPPUNIT_ASSERT(myMesh6->getNumberOfNodes() == MED_INVALID);
01896     CPPUNIT_ASSERT(myMesh6->getCoordinateptr() == NULL);
01897 
01898     myMesh4->removeReference();
01899     myMesh5->removeReference();
01900     myMesh6->removeReference();
01901 
01902     MESH* myMesh7 = MEDMEMTest_createTestMesh();
01903     vector< vector<double> > myBndBox;
01904     try
01905     {
01906       myBndBox = myMesh7->getBoundingBox();
01907     }
01908     catch (const std::exception &e)
01909     {
01910       CPPUNIT_FAIL(e.what());
01911     }
01912     catch (...)
01913     {
01914       CPPUNIT_FAIL("Unknown exception");
01915     }
01916 
01917     //cout<<"Bounding box for createTestMesh()"<<endl;
01918     double refBox[6] = 
01919       {
01920         -2, -2, 0, 2, 2, 5 
01921       };
01922     for(int i = 0, iB = 0; i < myBndBox.size(); i++)
01923     {
01924       for(int j = 0; j < myBndBox[i].size(); j++)
01925         CPPUNIT_ASSERT_DOUBLES_EQUAL( refBox[iB++], myBndBox[i][j], 1e-6);
01926       //cout<<" "<< myBndBox[i][j]<<" ";
01927       //cout<<endl;
01928     }
01929 
01930     double CoorPoint[3] = 
01931       {
01932         0.0,  0.0, 1.0
01933       }; //n2
01934     int idxElem;
01935     try
01936     {
01937       idxElem = myMesh7->getElementContainingPoint(CoorPoint);
01938     }
01939     catch (const std::exception &e)
01940     {
01941       CPPUNIT_FAIL(e.what());
01942     }
01943     catch (...)
01944     {
01945       CPPUNIT_FAIL("Unknown exception");
01946     }
01947     CPPUNIT_ASSERT(idxElem != -1);
01948 
01949     double CoorNoPoint[3] = 
01950       {
01951         5.0,  0.0, -5.0
01952       }; //there is no such point
01953     int idxNoElem;
01954     try
01955     {
01956       idxNoElem = myMesh7->getElementContainingPoint(CoorNoPoint);
01957     }
01958     catch (const std::exception &e)
01959     {
01960       CPPUNIT_FAIL(e.what());
01961     }
01962     catch (...)
01963     {
01964       CPPUNIT_FAIL("Unknown exception");
01965     }
01966     myMesh7->removeReference();
01967     CPPUNIT_ASSERT(idxNoElem == -1);
01968 
01970     // TEST 5: test desactivateFacesComputation() method //
01971     //         of driver: NPAL17670                      //
01973     double coords[54] = 
01974       {
01975         -0.215040, -0.107520, +0.000000,
01976         +0.000000, -0.107520, +0.000000, 
01977         +0.000000, +0.107520, +0.000000, 
01978         -0.215040, +0.107520, +0.000000, 
01979         +0.215040, -0.107520, +0.000000, 
01980         +0.215040, +0.107520, +0.000000, 
01981         -0.215040, -0.107520, +1.500000, 
01982         -0.215040, -0.107520, +4.080623, 
01983         +0.000000, -0.107520, +1.500000, 
01984         +0.000000, -0.107520, +4.080623, 
01985         +0.000000, +0.107520, +1.500000, 
01986         +0.000000, +0.107520, +4.080623, 
01987         -0.215040, +0.107520, +1.500000, 
01988         -0.215040, +0.107520, +4.080623, 
01989         +0.215040, -0.107520, +1.500000, 
01990         +0.215040, -0.107520, +4.080623, 
01991         +0.215040, +0.107520, +1.500000, 
01992         +0.215040, +0.107520, +4.080623 
01993       };
01994 
01995     int connQuad4[] = 
01996       {
01997         2 ,         5   ,      15   ,       9 ,
01998         10,         16  ,       18  ,       12, 
01999         11 ,        13  ,       14  ,       12, 
02000         7   ,       9   ,      11   ,      13, 
02001         3   ,       4   ,      13   ,      11, 
02002         1   ,       2   ,       9   ,       7, 
02003         1   ,       2   ,       3   ,       4, 
02004         15   ,      17  ,       18  ,       16 ,
02005         5    ,      6   ,      17   ,      15, 
02006         9    ,     15   ,      17   ,      11 ,
02007         13    ,      7  ,        8  ,       14, 
02008         4     ,     1   ,       7   ,      13, 
02009         9     ,    11   ,      12   ,      10, 
02010         8     ,    10   ,      12   ,      14, 
02011         2     ,     5   ,       6   ,       3, 
02012         17    ,     11,         12  ,       18 ,
02013         2     ,     3  ,       11   ,       9, 
02014         6     ,     3  ,       11   ,      17,
02015         7     ,     9  ,       10   ,       8, 
02016         9     ,    15  ,       16   ,      10 
02017       };
02018 
02019     int connHexa8[] = 
02020       {
02021         3       ,   2      ,    1   ,       4    ,     11     ,     9    ,      7    ,     13, 
02022         17     ,    15     ,     9    ,     11   ,      18    ,     16   ,      10   ,      12, 
02023         11    ,      9     ,     7    ,     13   ,      12    ,     10   ,       8   ,       14 ,
02024         6 ,         5      ,    2     ,     3    ,     17     ,    15    ,      9    ,     11
02025       };
02026 
02027     int bottom[2] = 
02028       {
02029         7,15
02030       };
02031     MED_EN::medGeometryElement bottomTypes[1] = 
02032       {
02033         MED_EN::MED_QUAD4
02034       };
02035     int bottomIndex[2] = 
02036       {
02037         1,3
02038       };
02039     int bottomNbOfElts[1] = 
02040       {
02041         2
02042       };
02043 
02044     MESHING* meshing = new MESHING();
02045     meshing->setName( "TESTMESH" );
02046     const int nFaces=20;
02047     const int nNodes=18;
02048     meshing->setCoordinates(3, nNodes, coords, "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
02049     string coordname[3] = 
02050       {
02051         "x", "y", "z" 
02052       };
02053     meshing->setCoordinatesNames(coordname);
02054     string coordunit[3] = 
02055       {
02056         "m", "m", "m" 
02057       };
02058     meshing->setCoordinatesUnits(coordunit);
02059     //Cell connectivity info for classical elts
02060     const MED_EN::medGeometryElement classicalTypesCell[1]=
02061       {
02062         MED_EN::MED_HEXA8
02063       };
02064     const int nbOfCellElts[1]=
02065       {
02066         4
02067       };
02068     meshing->setNumberOfTypes(1,MED_EN::MED_CELL);
02069     meshing->setTypes(classicalTypesCell,MED_EN::MED_CELL);
02070     meshing->setNumberOfElements(nbOfCellElts,MED_EN::MED_CELL);
02071     //Face connectivity info for classical elts
02072     const MED_EN::medGeometryElement classicalTypesFace[1]=
02073       {
02074         MED_EN::MED_QUAD4
02075       };
02076     const int nbOfFaceElts[1]=
02077       {
02078         nFaces
02079       };
02080     meshing->setNumberOfTypes(1,MED_EN::MED_FACE);
02081     meshing->setTypes(classicalTypesFace,MED_EN::MED_FACE);
02082     meshing->setNumberOfElements(nbOfFaceElts,MED_EN::MED_FACE);
02083     //All cell conn
02084     meshing->setConnectivity(MED_EN::MED_CELL,MED_EN::MED_HEXA8,connHexa8);
02085     //All face conn
02086     meshing->setConnectivity(MED_EN::MED_FACE,MED_EN::MED_QUAD4,connQuad4);
02087     //Adding some groups on faces
02088     addMedFacesGroup( *meshing, 4,  bottom, "BottomFaces",bottomTypes,bottomIndex,bottomNbOfElts,1) ;
02089     //addMedFacesGroupAll( *meshing, "AllFaces");
02090     //writing...
02091     int id=meshing->addDriver(MED_DRIVER,filenameout,meshing->getName());
02092     meshing->write(id);
02093     // Field writing
02094     const SUPPORT *sup=meshing->getSupportOnAll( MED_FACE );
02095     FIELD<double> * field = new FIELD<double>(sup, 1);
02096     field->setName("temperature");
02097     field->setComponentName(1,"T"); field->setMEDComponentUnit(1,"K");
02098     double *tab=(double *)field->getValue();
02099     for(int i=0;i<nFaces;i++)
02100       tab[i]=i*(1.22);
02101     field->setIterationNumber(0);
02102     field->setOrderNumber(-1);
02103     field->setTime(12.);
02104     id=field->addDriver(MED_DRIVER,filenameout,field->getName());
02105     field->write(id);
02106     field->removeReference();
02107     meshing->removeReference();
02108     //
02109     MESH mesh;
02110     MED_MESH_RDONLY_DRIVER drv(filenameout,&mesh);
02111     drv.desactivateFacesComputation();
02112     drv.setMeshName("TESTMESH");
02113     mesh.addDriver(drv);
02114     mesh.read();
02115     const int *conn=mesh.getConnectivity(MED_NODAL,MED_FACE,MED_ALL_ELEMENTS);
02116     for (int j = 0; j < nFaces; j++) 
02117     {
02118       for (int k = 0; k < 4; k++)
02119         CPPUNIT_ASSERT_EQUAL(conn[4*j+k], connQuad4[4*j+k]);
02120     }
02121     FIELD<double> f;
02122     f.addDriver(MED_DRIVER,filenameout,"temperature");
02123     f.setIterationNumber(0);
02124     f.setOrderNumber(-1);
02125     CPPUNIT_ASSERT_NO_THROW( f.read() );
02126 
02128     // TEST 6: Test Reading of a Field with given Mesh.     //
02129     // Group from the Mesh must be taken for Field Support. //
02131     {
02132       // mesh creation
02133       MESHING* mesh_prof = new MESHING();
02134       mesh_prof->setName("TESTMESH");
02135       mesh_prof->setCoordinates(3, nNodes, coords, "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
02136       mesh_prof->setCoordinatesNames(coordname);
02137       mesh_prof->setCoordinatesUnits(coordunit);
02138 
02139       //Cell connectivity info for classical elts
02140       //mesh_prof->setNumberOfTypes(1,MED_EN::MED_CELL);
02141       //mesh_prof->setTypes(classicalTypesCell,MED_EN::MED_CELL);
02142       //mesh_prof->setNumberOfElements(nbOfCellElts,MED_EN::MED_CELL);
02143       //mesh_prof->setConnectivity(connHexa8,MED_EN::MED_CELL,MED_EN::MED_HEXA8);
02144 
02145       //Face connectivity info for classical elts
02146       //mesh_prof->setNumberOfTypes(1,MED_EN::MED_FACE);
02147       //mesh_prof->setTypes(classicalTypesFace,MED_EN::MED_FACE);
02148       //mesh_prof->setNumberOfElements(nbOfFaceElts,MED_EN::MED_FACE);
02149       //mesh_prof->setConnectivity(connQuad4,MED_EN::MED_FACE,MED_EN::MED_QUAD4);
02150       mesh_prof->setNumberOfTypes(1,MED_EN::MED_CELL);
02151       mesh_prof->setTypes(classicalTypesFace,MED_EN::MED_CELL);
02152       mesh_prof->setNumberOfElements(nbOfFaceElts,MED_EN::MED_CELL);
02153       mesh_prof->setConnectivity(MED_EN::MED_CELL,MED_EN::MED_QUAD4,connQuad4);
02154 
02155       //Adding some groups on faces
02156       GROUP *faces_prof=new GROUP;
02157       faces_prof->setName("BottomFaces");
02158       faces_prof->setMesh(mesh_prof);
02159       //faces_prof->setEntity(MED_EN::MED_FACE);
02160       faces_prof->setEntity(MED_EN::MED_CELL);
02161       faces_prof->setNumberOfGeometricType(1);
02162       faces_prof->setGeometricType(bottomTypes);
02163       faces_prof->setNumberOfElements(bottomNbOfElts);
02164       faces_prof->setNumber(bottomIndex, bottom);
02165       mesh_prof->addGroup(*faces_prof);
02166 
02167       // Field creation
02168       FIELD<double> * field_prof = new FIELD<double>(faces_prof, 1);
02169       faces_prof->removeReference();
02170       field_prof->setName("temperature");
02171       field_prof->setComponentName(1,"T");
02172       field_prof->setMEDComponentUnit(1,"K");
02173       double *tab = (double *)field_prof->getValue();
02174       for (int i = 0; i < 2; i++)
02175         tab[i] = i*(1.22);
02176       field_prof->setTime(12.);
02177 
02178       // Writing...
02179       int id_prof = mesh_prof->addDriver(MED_DRIVER, filename_profiles_wr, mesh_prof->getName());
02180       mesh_prof->write(id_prof);
02181       mesh_prof->rmDriver(id_prof);
02182 
02183       // Field writing
02184       id_prof = field_prof->addDriver(MED_DRIVER, filename_profiles_wr, field_prof->getName());
02185       field_prof->write(id_prof);
02186 
02187       field_prof->removeReference();
02188       mesh_prof->removeReference();
02189 
02190       // Reading...
02191       MESH *mesh_rd=new MESH(MED_DRIVER, filename_profiles_wr, "TESTMESH");
02192       FIELD<double> *field_rd=new FIELD<double>(MED_DRIVER, filename_profiles_wr, "temperature", -1, -1, mesh_rd);
02193       //const vector <GROUP*> groups_rd = mesh_rd.getGroups(MED_EN::MED_FACE);
02194       vector <GROUP*> groups_rd = mesh_rd->getGroups(MED_EN::MED_CELL);
02195       vector <SUPPORT*> supports_rd(groups_rd.size());//Because virtual inheritance
02196       int iii=0;
02197       for(vector <GROUP*>::iterator iter=groups_rd.begin();iter!=groups_rd.end();iter++,iii++)
02198         supports_rd[iii]=dynamic_cast<SUPPORT*>(*iter);
02199       CPPUNIT_ASSERT(find(supports_rd.begin(),supports_rd.end(),field_rd->getSupport()) != supports_rd.end());
02200       field_rd->removeReference();
02201       mesh_rd->removeReference();
02202     }
02203 }
02204