Back to index

salome-med  6.5.0
MEDMEMTest_MeshAndMeshing_fault.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 <cppunit/Message.h>
00029 #include <cppunit/TestAssert.h>
00030 
00031 #include <sstream>
00032 #include <cmath>
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 double dmax(double x, double y) 
00045 {
00046 return (x>y)?x:y;
00047 }
00048 
00049 double dmin(double x, double y) 
00050 {
00051 return (x>y)?y:x;
00052 }
00053 
00222 void MEDMEMTest_testMeshAndMeshing()
00223 {
00224   string filename = getResourceFile("pointe.med");
00225   string meshname = "maa1";
00226   string filenameout21 = makeTmpFile("myMeshWrite4_pointe21.med");
00227 
00228   // To remove tmp files from disk
00229   MEDMEMTest_TmpFilesRemover aRemover;
00230   aRemover.Register(filenameout21);
00231 
00233   // TEST 1 //
00235 
00236   MESH * myMesh= new MESH();
00237   myMesh->setName("FIRST_MESH");
00238   CPPUNIT_ASSERT(myMesh != NULL);
00239 
00240   //test operator <<
00241   //#ifdef ENABLE_FAULTS
00242   //CPPUNIT_ASSERT_NO_THROW(cout << *myMesh << endl);
00243   //#endif
00244   //#ifdef ENABLE_FORCED_FAILURES
00245   //  CPPUNIT_FAIL("ERROR: operator << : if mesh is empty then attempt"
00246   //               " to get values from null object causes error");
00247   //#endif
00248 
00249   //test operator =
00250   MESH myMesh1 = *myMesh;
00251 
00252   //deepCompare
00253   //#ifdef ENABLE_FAULTS
00254   bool isEqual = false;
00255   CPPUNIT_ASSERT_NO_THROW(isEqual = myMesh1.deepCompare(*myMesh));
00256   CPPUNIT_ASSERT(isEqual);
00257   //#endif
00258   //#ifdef ENABLE_FORCED_FAILURES
00259   //  CPPUNIT_FAIL("ERROR: deepCompare(...) fails if mesh is empty");
00260   //#endif
00261 
00262   //ensure it imposible to compare meshes
00263   MESH *myMeshPointer =  myMesh;
00264   //test operator ==
00265   CPPUNIT_ASSERT(*myMeshPointer == *myMesh);
00266 
00267   delete myMesh;
00268 
00269   //set a MESH object
00270   MESHING *myMeshing=new MESHING;
00271   myMeshing->setName("meshing");
00272   // define coordinates
00273 
00274   int SpaceDimension = 3;
00275   int NumberOfNodes = 19;
00276   double Coordinates[57] = 
00277     {
00278       0.0, 0.0, 0.0,
00279       0.0, 0.0, 1.0,
00280       2.0, 0.0, 1.0,
00281       0.0, 2.0, 1.0,
00282       -2.0, 0.0, 1.0,
00283       0.0, -2.0, 1.0,
00284       1.0, 1.0, 2.0,
00285       -1.0, 1.0, 2.0,
00286       -1.0, -1.0, 2.0,
00287       1.0, -1.0, 2.0,
00288       1.0, 1.0, 3.0,
00289       -1.0, 1.0, 3.0,
00290       -1.0, -1.0, 3.0,
00291       1.0, -1.0, 3.0,
00292       1.0, 1.0, 4.0,
00293       -1.0, 1.0, 4.0,
00294       -1.0, -1.0, 4.0,
00295       1.0, -1.0, 4.0,
00296       0.0, 0.0, 5.0
00297     };
00298   try
00299     {
00300       myMeshing->setCoordinates(SpaceDimension,NumberOfNodes,Coordinates,"CARTESIAN",MED_FULL_INTERLACE);
00301     }
00302   catch (const std::exception &e)
00303     {
00304       CPPUNIT_FAIL(e.what());
00305     }
00306   catch (...)
00307     {
00308       CPPUNIT_FAIL("Unknown exception");
00309     }
00310 
00311   string Names[3] = 
00312     {
00313       "X","Y","Z"
00314     };
00315   try
00316     {
00317       myMeshing->setCoordinatesNames(Names);
00318     }
00319   catch (const std::exception &e)
00320     {
00321       CPPUNIT_FAIL(e.what());
00322     }
00323   catch (...)
00324     {
00325       CPPUNIT_FAIL("Unknown exception");
00326     }
00327 
00328   string Units[3] = 
00329     {
00330       "cm","cm","cm"
00331     };
00332   try
00333     {
00334       myMeshing->setCoordinatesUnits(Units);
00335     }
00336   catch (const std::exception &e)
00337     {
00338       CPPUNIT_FAIL(e.what());
00339     }
00340   catch (...)
00341     {
00342       CPPUNIT_FAIL("Unknown exception");
00343     }
00344 
00345   // define conectivities
00346 
00347   // cell part
00348 
00349   const int NumberOfTypes = 3;
00350   medGeometryElement Types[NumberOfTypes] = 
00351     {
00352       MED_TETRA4,MED_PYRA5,MED_HEXA8
00353     };
00354   const int NumberOfElements[NumberOfTypes] = 
00355     {
00356       12,2,2
00357     };
00358 
00359   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfTypes(NumberOfTypes,MED_CELL));
00360 
00361   CPPUNIT_ASSERT_NO_THROW(myMeshing->setTypes(Types,MED_CELL));
00362 
00363   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfElements(NumberOfElements,MED_CELL));
00364 
00365   const int sizeTetra = 12*4;
00366   int ConnectivityTetra[sizeTetra]=
00367     {
00368       1,2,3,6,
00369       1,2,4,3,
00370       1,2,5,4,
00371       1,2,6,5,
00372       2,7,4,3,
00373       2,8,5,4,
00374       2,9,6,5,
00375       2,10,3,6,
00376       2,7,3,10,
00377       2,8,4,7,
00378       2,9,5,8,
00379       2,10,6,9
00380     };
00381 
00382   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_CELL, MED_TETRA4, ConnectivityTetra ));
00383 
00384   int ConnectivityPyra[2*5]=
00385     {
00386       7,8,9,10,2,
00387       15,18,17,16,19
00388     };
00389 
00390   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_CELL, MED_PYRA5, ConnectivityPyra ));
00391 
00392   int ConnectivityHexa[2*8]=
00393     {
00394       11,12,13,14,7,8,9,10,
00395       15,16,17,18,11,12,13,14
00396     };
00397 
00398   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_CELL, MED_HEXA8, ConnectivityHexa ));
00399 
00400   // face part
00401   const int NumberOfFacesTypes = 2;
00402   medGeometryElement FacesTypes[NumberOfFacesTypes] = 
00403     {
00404       MED_TRIA3,MED_QUAD4
00405     };
00406   const int NumberOfFacesElements[NumberOfFacesTypes] = 
00407     {
00408       4,4
00409     };
00410 
00411   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfTypes(NumberOfFacesTypes,MED_FACE));
00412   CPPUNIT_ASSERT_NO_THROW(myMeshing->setTypes(FacesTypes,MED_FACE));
00413   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfElements(NumberOfFacesElements,MED_FACE));
00414   const int nbTria = 4;
00415   const int sizeTria = nbTria*3;
00416   int ConnectivityTria[sizeTria]=
00417     {
00418       1,4,3,
00419       1,5,4,
00420       1,6,5,
00421       1,3,6
00422     };
00423 
00424   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_FACE, MED_TRIA3, ConnectivityTria ));
00425   const int nbQua = 4;
00426   int ConnectivityQua[nbQua*4]=
00427     {
00428       7,8,9,10,
00429       11,12,13,14,
00430       11,7,8,12,
00431       12,8,9,13
00432     };
00433 
00434   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_FACE, MED_QUAD4, ConnectivityQua ));
00435 
00436   int meshDimension = SpaceDimension; // because there 3D cells in the mesh
00437 
00438   // edge part
00439 
00440   // not yet implemented : if set, results are unpredictable.
00441 
00442   // Some groups :
00443 
00444   // Node :
00445   {
00446     GROUP *myGroup=new GROUP;
00447     myGroup->setName("SomeNodes");
00448     myGroup->setMesh(myMeshing);
00449     myGroup->setEntity(MED_NODE);
00450     myGroup->setNumberOfGeometricType(1);
00451     medGeometryElement myTypes[1] = 
00452       {
00453         MED_NONE
00454       };
00455     myGroup->setGeometricType(myTypes);
00456     const int myNumberOfElements[1] = 
00457       {
00458         4
00459       };
00460     myGroup->setNumberOfElements(myNumberOfElements);
00461     const int index[1+1] = 
00462       {
00463         1,5
00464       };
00465     const int value[4]= 
00466       {
00467         1,4,5,7
00468       };
00469     myGroup->setNumber(index,value);
00470     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00471     myGroup->removeReference();
00472   }
00473   {
00474     GROUP *myGroup=new GROUP;
00475     myGroup->setName("OtherNodes");
00476     myGroup->setMesh(myMeshing);
00477     myGroup->setEntity(MED_NODE);
00478     myGroup->setNumberOfGeometricType(1);
00479     medGeometryElement myTypes[1] = 
00480       {
00481         MED_NONE
00482       };
00483     myGroup->setGeometricType(myTypes);
00484     const int myNumberOfElements[1] = 
00485       {
00486         3
00487       };
00488     myGroup->setNumberOfElements(myNumberOfElements);
00489     const int index[1+1] = 
00490       {
00491         1,4
00492       };
00493     const int value[3]= 
00494       {
00495         2,3,6
00496       };
00497     myGroup->setNumber(index,value);
00498     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00499     myGroup->removeReference();
00500   }
00501 
00502   // Cell :
00503   {
00504     GROUP *myGroup=new GROUP;
00505     myGroup->setName("SomeCells");
00506     myGroup->setMesh(myMeshing);
00507     myGroup->setEntity(MED_CELL);
00508     myGroup->setNumberOfGeometricType(3);
00509     medGeometryElement myTypes[3] = 
00510       {
00511         MED_TETRA4,MED_PYRA5,MED_HEXA8
00512       };
00513     myGroup->setGeometricType(myTypes);
00514     const int myNumberOfElements[3] = 
00515       {
00516         4,1,2
00517       };
00518     myGroup->setNumberOfElements(myNumberOfElements);
00519     const int index[3+1] = 
00520       {
00521         1,5,6,8
00522       };
00523     const int value[4+1+2]=
00524       {
00525         2,7,8,12,
00526         13,
00527         15,16
00528       };
00529     myGroup->setNumber(index,value);
00530     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00531     myGroup->removeReference();
00532   }
00533   {
00534     GROUP *myGroup=new GROUP;
00535     myGroup->setName("OtherCells");
00536     myGroup->setMesh(myMeshing);
00537     myGroup->setEntity(MED_CELL);
00538     myGroup->setNumberOfGeometricType(2);
00539     medGeometryElement myTypes[] = 
00540       {
00541         MED_TETRA4,MED_PYRA5
00542       };
00543     myGroup->setGeometricType(myTypes);
00544     const int myNumberOfElements[] = 
00545       {
00546         4,1
00547       };
00548     myGroup->setNumberOfElements(myNumberOfElements);
00549     const int index[2+1] = 
00550       {
00551         1,5,6
00552       };
00553     const int value[4+1]=
00554       {
00555         3,4,5,9,
00556         14
00557       };
00558     myGroup->setNumber(index,value);
00559     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00560     myGroup->removeReference();
00561   }
00562 
00563   // Face :
00564   {
00565     GROUP *myGroup=new GROUP;
00566     myGroup->setName("SomeFaces");
00567     myGroup->setMesh(myMeshing);
00568     myGroup->setEntity(MED_FACE);
00569     myGroup->setNumberOfGeometricType(2);
00570     medGeometryElement myTypes[2] = 
00571       {
00572         MED_TRIA3,MED_QUAD4
00573       };
00574     myGroup->setGeometricType(myTypes);
00575     const int myNumberOfElements[2] = 
00576       {
00577         2,3
00578       };
00579     myGroup->setNumberOfElements(myNumberOfElements);
00580     const int index[2+1] = 
00581       {
00582         1,3,6
00583       };
00584     const int value[2+3]=
00585       {
00586         2,4,
00587         5,6,8
00588       };
00589     myGroup->setNumber(index,value);
00590     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00591     myGroup->removeReference();
00592   }
00593   {
00594     GROUP *myGroup=new GROUP;
00595     myGroup->setName("OtherFaces");
00596     myGroup->setMesh(myMeshing);
00597     myGroup->setEntity(MED_FACE);
00598     myGroup->setNumberOfGeometricType(1);
00599     medGeometryElement myTypes[1] = 
00600       {
00601         MED_TRIA3
00602       };
00603     myGroup->setGeometricType(myTypes);
00604     const int myNumberOfElements[1] = 
00605       {
00606         2
00607       };
00608     myGroup->setNumberOfElements(myNumberOfElements);
00609     const int index[1+1] = 
00610       {
00611         1,3
00612       };
00613     const int value[2]=
00614       {
00615         1,3
00616       };
00617     myGroup->setNumber(index,value);
00618     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
00619     myGroup->removeReference();
00620   }
00621 
00622   //test Mesh(MESH &m)
00623   {
00624     MESH * myMesh2 = new MESH( *myMeshing );
00625     CPPUNIT_ASSERT(myMesh2->deepCompare(*myMeshing));
00626 
00627     cout<<*myMesh2<<endl;
00628     ostringstream os;
00629     os << * myMesh2;
00630     CPPUNIT_ASSERT(os.str() != "");
00631 
00632     CPPUNIT_ASSERT_EQUAL(myMesh2->getName(),(string)"meshing");
00633     CPPUNIT_ASSERT((myMesh2->getDescription()).size() == 0);
00634     myMesh2->setDescription("This class contains all information related to a 'meshing' mesh ");
00635     CPPUNIT_ASSERT((myMesh2->getDescription()).size() != 0);
00636 
00637     CPPUNIT_ASSERT(myMesh2->getSpaceDimension() == SpaceDimension);
00638     CPPUNIT_ASSERT(myMesh2->getMeshDimension() == meshDimension);
00639     CPPUNIT_ASSERT(myMesh2->getNumberOfNodes() == NumberOfNodes);
00640 
00641     const COORDINATE* coord = myMesh2->getCoordinateptr();
00642     try
00643       {
00644         CPPUNIT_ASSERT(myMesh2->getCoordinatesSystem() != "catresian");
00645       }
00646     catch (const std::exception &e)
00647       {
00648         CPPUNIT_FAIL(e.what());
00649       }
00650     catch (...)
00651       {
00652         CPPUNIT_FAIL("Unknown exception");
00653       }
00654 
00655     const string * units;
00656     try
00657       {
00658         units = myMesh2->getCoordinatesUnits();
00659         for (int axe = 0; axe < SpaceDimension; axe++) 
00660           {
00661             string verif = coord->getCoordinateUnit(axe+1);
00662             CPPUNIT_ASSERT(verif == units[axe]);
00663           }
00664       }
00665     catch (const std::exception &e)
00666       {
00667         CPPUNIT_FAIL(e.what());
00668       }
00669     catch (...)
00670       {
00671         CPPUNIT_FAIL("Unknown exception");
00672       }
00673 
00674     const string * noms;
00675     try
00676       {
00677         noms = myMesh2->getCoordinatesNames();
00678         for (int axe = 0; axe < SpaceDimension; axe++) 
00679           {
00680             string verif = coord->getCoordinateName(axe+1);
00681             CPPUNIT_ASSERT(verif == noms[axe]);
00682           }
00683       }
00684     catch (const std::exception &e)
00685       {
00686         CPPUNIT_FAIL(e.what());
00687       }
00688     catch (...)
00689       {
00690         CPPUNIT_FAIL("Unknown exception");
00691       }
00692 
00693     try
00694       {
00695         const double * coor2 = myMesh2->getCoordinates(MED_FULL_INTERLACE);
00696 
00697         for (int axe = 0; axe < SpaceDimension; axe++) 
00698           {
00699             try
00700               {
00701                 for (int num = 0; num < NumberOfNodes; num++) 
00702                   {
00703                     try
00704                       {
00705                         const double d = myMesh2->getCoordinate(num + 1, axe + 1);
00706                         CPPUNIT_ASSERT(fabs(d - coor2[(num * SpaceDimension)+axe]) < 0.001);
00707                       }
00708                     catch (const std::exception &e)
00709                       {
00710                         CPPUNIT_FAIL(e.what());
00711                       }
00712                     catch (...)
00713                       {
00714                         CPPUNIT_FAIL("Unknown exception");
00715                       }
00716                   }
00717               }
00718             catch (const std::exception &e)
00719               {
00720                 CPPUNIT_FAIL(e.what());
00721               }
00722             catch (...)
00723               {
00724                 CPPUNIT_FAIL("Unknown exception");
00725               }
00726           }
00727       }
00728     catch (const std::exception &e)
00729       {
00730         CPPUNIT_FAIL(e.what());
00731       }
00732     catch (...)
00733       {
00734         CPPUNIT_FAIL("Unknown exception");
00735       }
00736 
00737     const CONNECTIVITY * myConnectivity = myMesh2->getConnectivityptr();
00738 
00739     // MED_EN::MED_CELL
00740     MED_EN::medEntityMesh entity = myConnectivity->getEntity();
00741     CPPUNIT_ASSERT_EQUAL(MED_CELL, entity);
00742 
00743     int typesNb;
00744     CPPUNIT_ASSERT_NO_THROW(typesNb= myConnectivity->getNumberOfTypes(entity));
00745     CPPUNIT_ASSERT_EQUAL(NumberOfTypes, typesNb);
00746 
00747     const MED_EN::medGeometryElement * Types1;
00748     CPPUNIT_ASSERT_NO_THROW(Types1 = myMesh2->getTypes(entity));
00749 
00750     medConnectivity myMedConnect;
00751     bool existConnect = false;
00752     if (myMesh2->existConnectivity(MED_NODAL, entity))
00753       {
00754         existConnect = true;
00755         myMedConnect = MED_NODAL;
00756       }
00757     else if(myMesh2->existConnectivity(MED_DESCENDING, entity))
00758       {
00759         existConnect = true;
00760         myMedConnect = MED_DESCENDING;
00761       }
00762 
00763     for(int t = 0; t < NumberOfTypes; t++ )
00764       {
00765         CPPUNIT_ASSERT_EQUAL(Types1[t], Types[t]);
00766         int NumberOfElements1 = 0;
00767         CPPUNIT_ASSERT_NO_THROW(NumberOfElements1 = myMesh2->getNumberOfElements(entity, Types1[t]));
00768         CPPUNIT_ASSERT_EQUAL(NumberOfElements1, NumberOfElements[t]);
00769         if(existConnect)
00770           {
00771             const int * connectivity;
00772             const int * connectivity_index;
00773             CPPUNIT_ASSERT_NO_THROW(connectivity = myMesh2->getConnectivity
00774                                     (myMedConnect, entity, Types1[t]));
00775             connectivity_index = myMesh2->getConnectivityIndex(myMedConnect, entity);
00776             for (int j = 0; j < NumberOfElements1; j++) 
00777               {
00778                 cout<<"!!!!!!!!!!!!!!!"<<endl;
00779                 for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
00780                   cout << connectivity[k-1] << " ";
00781                 cout << endl;
00782               }
00783           }
00784       }
00785 
00786     const CELLMODEL* myCellModel = myMesh2->getCellsTypes(entity);
00787     string* TypeNames;
00788     CPPUNIT_ASSERT_NO_THROW(TypeNames = myMesh2->getCellTypeNames(entity));
00789 
00790     for(int k = 0; k < NumberOfTypes; k++ )
00791       {
00792         CPPUNIT_ASSERT_EQUAL(TypeNames[k], myCellModel[k].getName());
00793       }
00794     delete [] TypeNames;
00795 
00796     const int* myGlobalNbIdx;
00797     CPPUNIT_ASSERT_NO_THROW(myGlobalNbIdx = myMesh2->getGlobalNumberingIndex(MED_FACE));
00798     for(int i = 0; i <= NumberOfFacesTypes; i++)
00799       {
00800         if(i == NumberOfFacesTypes)
00801           {
00802             CPPUNIT_ASSERT_EQUAL(myGlobalNbIdx[i],nbTria+nbQua+1);
00803             CPPUNIT_ASSERT_THROW(myMesh2->getElementType(MED_FACE, myGlobalNbIdx[i]), MEDEXCEPTION);
00804             break;
00805           }
00806         cout<<"Global number of first element of each geom type : "<<myGlobalNbIdx[i]<<endl;
00807       }
00808     {
00809       const int * ReverseNodalConnectivity;
00810 
00811       // Show Reverse Nodal Connectivity
00812       //#ifndef ENABLE_FAULTS
00813       // (BUG) CONNECTIVITY::_numberOfNodes is not set
00814       ((CONNECTIVITY*)myConnectivity)->setNumberOfNodes(NumberOfNodes);
00815       //#endif
00816       //#ifdef ENABLE_FORCED_FAILURES
00817       //      CPPUNIT_FAIL("ERROR in CONNECTIVITY::calculateReverseNodalConnectivity()"
00818       //                   " because myMesh2->_connectivity->_numberOfNodes is not set");
00819       //#endif
00820 
00821       CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivity = myMesh2->getReverseConnectivity(MED_NODAL, entity));
00822       CPPUNIT_ASSERT_NO_THROW(myMesh2->getReverseConnectivityLength(MED_NODAL, entity));
00823       const int * ReverseNodalConnectivityIndex = myMesh2->getReverseConnectivityIndex(MED_NODAL, entity);
00824       const int ReverseIdxLength = myMesh2->getReverseConnectivityIndexLength(MED_NODAL, entity);
00825       CPPUNIT_ASSERT(ReverseIdxLength == NumberOfNodes+1);
00826       for (int i = 0; i < NumberOfNodes; i++) 
00827         {
00828           cout << "Node "<< i+1 << " : ";
00829           for (int j = ReverseNodalConnectivityIndex[i]; j < ReverseNodalConnectivityIndex[i+1]; j++)
00830             cout << ReverseNodalConnectivity[j-1] << " ";
00831           cout << endl;
00832         }
00833 
00834       // Show Descending Connectivity
00835       int NumberOfElements1;
00836       const int * connectivity;
00837       const int * connectivity_index;
00838       myMesh2->calculateConnectivity(MED_DESCENDING, entity);
00839       try
00840         {
00841           NumberOfElements1 = myMesh2->getNumberOfElements(entity, MED_ALL_ELEMENTS);
00842           connectivity = myMesh2->getConnectivity( MED_DESCENDING, entity, MED_ALL_ELEMENTS);
00843           connectivity_index = myMesh2->getConnectivityIndex(MED_DESCENDING, entity);
00844         }
00845       catch (MEDEXCEPTION m) 
00846         {
00847           CPPUNIT_FAIL(m.what());
00848         }
00849 
00850       for (int j = 0; j < NumberOfElements1; j++) 
00851         {
00852           cout << "Element " << j+1 << " : ";
00853           for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
00854             cout << connectivity[k-1] << " ";
00855           cout << endl;
00856         }
00857 
00858       // getElementNumber
00859       if (myMesh2->existConnectivity(MED_NODAL, MED_FACE)) 
00860         {
00861           int myTr[3] = 
00862             {
00863               1,5,4
00864             };
00865           CPPUNIT_ASSERT_NO_THROW(myMesh2->getElementNumber(MED_NODAL,MED_FACE,MED_TRIA3,myTr));
00866         }
00867     }
00868 
00869     //test family and group
00870     int NumberOfGroups;
00871     CPPUNIT_ASSERT_THROW(myMesh2->getNumberOfGroups(MED_ALL_ENTITIES), MEDEXCEPTION);
00872     CPPUNIT_ASSERT_NO_THROW(NumberOfGroups = myMesh2->getNumberOfGroups(MED_CELL));
00873     CPPUNIT_ASSERT_EQUAL(NumberOfGroups, 2);
00874     vector<GROUP*> groups;
00875     CPPUNIT_ASSERT_NO_THROW(groups = myMesh2->getGroups(MED_CELL));
00876     CPPUNIT_ASSERT(groups.size() != 0);
00877     for(int nb = 1; nb <= NumberOfGroups; nb++ )
00878       {
00879         const GROUP* group;
00880         CPPUNIT_ASSERT_NO_THROW(group = myMesh2->getGroup(MED_CELL, nb));
00881         CPPUNIT_ASSERT_EQUAL(group->getName(), groups[nb-1]->getName());
00882       }
00883 
00884     int NumberOfFamilies;
00885     CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies = myMesh2->getNumberOfFamilies(MED_CELL));
00886     CPPUNIT_ASSERT_MESSAGE("Current mesh hasn't Families", NumberOfFamilies == 0);
00887 
00888     //create families - it's not possible to create, becase not all entities are defined
00889     CPPUNIT_ASSERT_THROW( myMesh2->createFamilies(),MEDEXCEPTION);
00890 
00891     /*CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies = myMesh2->getNumberOfFamilies(MED_CELL));
00892       CPPUNIT_ASSERT( NumberOfFamilies != 0);*/
00893 
00894     delete myMesh2;
00895   }
00896 
00898   // TEST 2: Polygon and Polyhedron(only NODAL connectivity)  //
00900 
00901   double CoordinatesPoly[57] = 
00902     {
00903       2.0, 3.0, 2.0,
00904       3.0, 2.0, 2.0,
00905       4.0, 1.0, 2.0,
00906       2.0, 0.0, 2.0,
00907       0.0, 1.0, 2.0,
00908       1.0, 2.0, 2.0,
00909       2.0, 3.0, 1.0,
00910       3.0, 2.0, 0.0,
00911       4.0, 1.0, 0.0,
00912       2.0, 0.0, 0.0,
00913       0.0, 1.0, 0.0,
00914       1.0, 2.0, 0.0,
00915       5.0, 3.0, 2.0,
00916       7.0, 2.0, 2.0,
00917       6.0, 0.0, 2.0,
00918       6.0, 3.0, 0.0,
00919       7.0, 2.0, 0.0,
00920       6.0, 0.0, -1.0,
00921       5.0, 1.0, -3.0
00922     };
00923 
00924   const int REFnodalConnOfFaces[91] = 
00925     {
00926       1, 2, 3, 4, 5, 6, -1,// Polyhedron 1
00927       1, 7, 8, 2,       -1,
00928       2, 8, 9, 3,       -1,
00929       4, 3, 9, 10,      -1,
00930       5, 4, 10, 11,     -1,
00931       6, 5, 11, 12,     -1,
00932       1, 6, 12, 7,      -1,
00933       7, 12, 8, 10,     -1,
00934       9, 8, 12, 11,     
00935 
00936       13, 14, 15, 3, 2, -1,// Polyhedron 2
00937       13, 2, 8, 16,     -1,
00938       14, 13, 16, 17,   -1,
00939       15, 14, 17, 15,   -1,
00940       17, 18, 15,       -1,
00941       18, 9, 3,         -1,
00942       15, 9, 2,         -1,
00943       3, 9, 8,          -1,
00944       8, 9, 17, 16,     -1,
00945       9, 18, 17
00946     };
00947   const int NumberOfFaces = 19;
00948   const int NumberOfPolyhedron = 2;
00949   const int nbOfPolygons = 2;
00950 
00951   const int REFpolyIndex[NumberOfPolyhedron+1] = 
00952     {
00953       1,47,92
00954     };
00955 
00956   double PolygonCoordinates[27] = 
00957     {
00958       2.0, 3.0, 12.0,
00959       3.0, 2.0, 12.0,
00960       4.0, 1.0, 12.0,
00961       2.0, 0.0, 12.0,
00962       0.0, 1.0, 12.0,
00963       1.0, 2.0, 12.0,
00964       5.0, 3.0, 12.0,
00965       7.0, 2.0, 12.0,
00966       6.0, 0.0, 12.0
00967     };
00968 
00969   const int REFpolygonFaces[11] = 
00970     {
00971       1, 2, 3, 4, 5, 6, // Polygon 1
00972       7, 8, 9, 3, 2
00973     }; // Polygon 2
00974 
00975   const int REFpolygonIndex[nbOfPolygons+1] = 
00976     {
00977       1, 7, 12
00978     };
00979 
00980   MESHING *myMeshingPoly=new MESHING;
00981   myMeshingPoly->setName("meshingpoly");
00982 
00983   const int NbOfTypes = 1;
00984   medGeometryElement TypesPoly[NbOfTypes] = 
00985     {
00986       MED_TETRA4
00987     };
00988   const int NbOfElements[NbOfTypes] = 
00989     {
00990       1
00991     };
00992 
00993   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setNumberOfTypes(NbOfTypes, MED_CELL));
00994 
00995   try
00996     {
00997       myMeshingPoly->setCoordinates(SpaceDimension, NumberOfNodes, CoordinatesPoly,
00998                                     "CARTESIAN", MED_FULL_INTERLACE);
00999     }
01000   catch (const std::exception &e)
01001     {
01002       CPPUNIT_FAIL(e.what());
01003     }
01004   catch (...)
01005     {
01006       CPPUNIT_FAIL("Unknown exception");
01007     }
01008 
01009   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setTypes(TypesPoly, MED_CELL));
01010   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setNumberOfElements(NbOfElements, MED_CELL));
01011 
01012   string Unit ="cm";
01013   for(int i = 0; i < SpaceDimension; i++ )
01014     {
01015       try
01016         {
01017           myMeshingPoly->setCoordinateName(Names[i],i);
01018         }
01019       catch (const std::exception &e)
01020         {
01021           CPPUNIT_FAIL(e.what());
01022         }
01023       catch (...)
01024         {
01025           CPPUNIT_FAIL("Unknown exception");
01026         }
01027       try
01028         {
01029           myMeshingPoly->setCoordinateUnit(Unit, i);
01030         }
01031       catch (const std::exception &e)
01032         {
01033           CPPUNIT_FAIL(e.what());
01034         }
01035       catch (...)
01036         {
01037           CPPUNIT_FAIL("Unknown exception");
01038         }
01039     }
01040 
01041   int ConnectivityTetraPoly[4*1]=
01042     {
01043       17, 9, 18, 19
01044     };
01045 
01046   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setConnectivity( MED_CELL, MED_TETRA4, ConnectivityTetraPoly ));
01047 
01048   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setConnectivity(MED_CELL, MED_POLYHEDRA,REFnodalConnOfFaces,REFpolyIndex));
01049 
01050   bool PolyConn = false;
01051   CPPUNIT_ASSERT_NO_THROW(PolyConn = myMeshingPoly->existConnectivity(MED_NODAL, MED_CELL));
01052   CPPUNIT_ASSERT(PolyConn);
01053     {
01054     CPPUNIT_ASSERT_EQUAL(NumberOfPolyhedron,
01055                          myMeshingPoly->getNumberOfElements(MED_CELL,MED_POLYHEDRA));
01056     CPPUNIT_ASSERT_NO_THROW( myMeshingPoly->calculateConnectivity (MED_NODAL,MED_FACE));
01057     CPPUNIT_ASSERT_EQUAL(NumberOfFaces-1, myMeshingPoly->getNumberOfElements(MED_FACE,MED_POLYGON)); // -1: one face is shared with tetra
01058     CPPUNIT_ASSERT_EQUAL(91,myMeshingPoly->getConnectivityLength(MED_NODAL,MED_CELL,MED_POLYHEDRA));
01059     const int * PolyConn;
01060     const int * PolyIdx;
01061     CPPUNIT_ASSERT_NO_THROW(PolyConn = myMeshingPoly->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS));
01062     CPPUNIT_ASSERT_NO_THROW(PolyIdx = myMeshingPoly->getConnectivityIndex(MED_NODAL,MED_CELL));
01063     for(int i = NbOfElements[0], iRef=0; i<NbOfElements[0]+NumberOfPolyhedron; i++)
01064       {
01065         int NodeIdxBegin = PolyIdx[i];
01066         int NodeIdxEnd = PolyIdx[i+1];
01067         for(int k = NodeIdxBegin; k < NodeIdxEnd; k++)
01068           CPPUNIT_ASSERT_EQUAL(REFnodalConnOfFaces[iRef++], PolyConn[k-1]);
01069       }
01070   }
01071 
01072   MESHING *myPolygonMeshing=new MESHING;
01073   myPolygonMeshing->setName("PolygonMeshing");
01074 
01075   medGeometryElement PolygonTypes[NbOfTypes] = 
01076     {
01077       MED_TRIA3
01078     };
01079   const int PolygonNumberOfElements[NbOfTypes] = 
01080     {
01081       2
01082     };
01083 
01084   CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setNumberOfTypes(NbOfTypes, MED_CELL));
01085 
01086   try
01087     {
01088       myPolygonMeshing->setCoordinates(SpaceDimension, NumberOfNodes, PolygonCoordinates,
01089                                        "CARTESIAN", MED_FULL_INTERLACE);
01090     }
01091   catch (const std::exception &e)
01092     {
01093       CPPUNIT_FAIL(e.what());
01094     }
01095   catch (...)
01096     {
01097       CPPUNIT_FAIL("Unknown exception");
01098     }
01099 
01100     CPPUNIT_ASSERT_EQUAL(SpaceDimension, myPolygonMeshing->getSpaceDimension());
01101     CPPUNIT_ASSERT_EQUAL(NumberOfNodes, myPolygonMeshing->getNumberOfNodes());
01102 
01103     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setTypes(PolygonTypes, MED_CELL));
01104     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setNumberOfElements(PolygonNumberOfElements, MED_CELL));
01105     CPPUNIT_ASSERT_EQUAL(2, myPolygonMeshing->getMeshDimension());
01106 
01107   try
01108     {
01109       myPolygonMeshing->setCoordinatesNames(Names);
01110     }
01111   catch (const std::exception &e)
01112     {
01113       CPPUNIT_FAIL(e.what());
01114     }
01115   catch (...)
01116     {
01117       CPPUNIT_FAIL("Unknown exception");
01118     }
01119 
01120   try
01121     {
01122       myPolygonMeshing->setCoordinatesUnits(Units);
01123     }
01124   catch (const std::exception &e)
01125     {
01126       CPPUNIT_FAIL(e.what());
01127     }
01128   catch (...)
01129     {
01130       CPPUNIT_FAIL("Unknown exception");
01131     }
01132 
01133   const int sizeTri = 3*2;
01134   int ConnectivityTri[sizeTri]=
01135     {
01136       1, 7, 2, 3, 9, 4
01137     };
01138 
01139   CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setConnectivity( MED_CELL, MED_TRIA3, ConnectivityTri ));
01140   CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setConnectivity( MED_CELL, MED_POLYGON, REFpolygonFaces, REFpolygonIndex ));
01141 
01142   bool PolygonConn = false;
01143   CPPUNIT_ASSERT_NO_THROW(PolygonConn = myPolygonMeshing->existConnectivity(MED_NODAL, MED_CELL));
01144   if(PolygonConn)
01145     {
01146       int Polytypes;
01147       CPPUNIT_ASSERT_NO_THROW(Polytypes = myPolygonMeshing->getNumberOfTypes(MED_CELL));
01148       CPPUNIT_ASSERT_EQUAL(NbOfTypes,Polytypes);
01149 
01150       const MED_EN::medGeometryElement * PolyTypes;
01151       CPPUNIT_ASSERT_NO_THROW(PolyTypes = myPolygonMeshing->getTypes(MED_CELL));
01152       CPPUNIT_ASSERT_EQUAL(PolyTypes[NbOfTypes-1],MED_POLYGON);
01153 
01154       for(int t = 0; t < Polytypes; t++)
01155         {
01156           CPPUNIT_ASSERT_NO_THROW( myPolygonMeshing->getNumberOfElements(MED_CELL, PolyTypes[t]));
01157         }
01158 
01159       medGeometryElement geomPolyElem;
01160       CPPUNIT_ASSERT_NO_THROW(geomPolyElem = myPolygonMeshing->getElementType(MED_CELL, 1));
01161       CPPUNIT_ASSERT_EQUAL(geomPolyElem, MED_TRIA3);
01162 
01163       CPPUNIT_ASSERT_EQUAL(myPolygonMeshing->getNumberOfElements(MED_CELL,MED_POLYGON),nbOfPolygons);
01164       CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->getConnectivityLength(MED_NODAL,MED_CELL,MED_POLYGON));
01165       myPolygonMeshing->removeReference();
01166       const int * PolygonConn;
01167       CPPUNIT_ASSERT_THROW(PolygonConn = myMeshingPoly->getConnectivity(MED_NODAL,MED_CELL,MED_POLYGON),MEDEXCEPTION);
01168     }
01169 
01171   // TEST : SUPPORT* sup = myMeshPointe->getSupportOnAll()) //
01173 
01174   //#ifdef ENABLE_FAULTS
01175   {
01176     MESH * myMeshPointe = new MESH(MED_DRIVER, filename, meshname);
01177     const SUPPORT* sup = myMeshPointe->getSupportOnAll(MED_CELL);
01178     CPPUNIT_ASSERT( sup->isOnAllElements() );
01179     CPPUNIT_ASSERT_EQUAL( myMeshPointe->getNumberOfTypes( sup->getEntity() ),
01180                           sup->getNumberOfTypes());
01181     CPPUNIT_ASSERT( sup->getNumber( MED_ALL_ELEMENTS ));
01182     myMeshPointe->removeReference();
01183   }
01184   //#endif
01185   //#ifdef ENABLE_FORCED_FAILURES
01186   //  CPPUNIT_FAIL("ERROR: can not create SUPPORT on mesh, read from pointe.med");
01187   //#endif
01188 
01190   // TEST 3: test MESH on  MEDMEMTest::createTestMesh()//
01192 
01193   MESH* myMesh3 = MEDMEMTest_createTestMesh();
01194 
01195   int MeshDim  = myMesh3->getMeshDimension();
01196   medEntityMesh constituentEntity;
01197   if (MeshDim==3) 
01198     {
01199       constituentEntity = MED_CELL;
01200     }
01201   if (MeshDim==2) 
01202     {
01203       constituentEntity = MED_FACE;
01204     }
01205   if (MeshDim==1) 
01206     {
01207       constituentEntity = MED_EDGE;
01208     }
01209 
01210   int SpaceDim = myMesh3->getSpaceDimension();
01211 
01212   // Show Reverse Nodal Connectivity
01213   const int* ReverseNodalConnectivity;
01214   const int* ReverseNodalConnectivityIndex;
01215   int ReverseLength;
01216   int ReverseIdxLength;
01217 
01218   CONNECTIVITY* myConnectivity3 = (CONNECTIVITY*)myMesh3->getConnectivityptr();
01219   myConnectivity3->setNumberOfNodes(myMesh3->getNumberOfNodes());
01220 
01221   CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivity= myMesh3->getReverseConnectivity(MED_NODAL, MED_CELL));
01222   CPPUNIT_ASSERT_NO_THROW(ReverseLength = myMesh3->getReverseConnectivityLength(MED_NODAL, MED_CELL));
01223   CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivityIndex = myMesh3->getReverseConnectivityIndex(MED_NODAL, MED_CELL));
01224   CPPUNIT_ASSERT_NO_THROW(ReverseIdxLength = myMesh3->getReverseConnectivityIndexLength(MED_NODAL, MED_CELL));
01225   CPPUNIT_ASSERT(ReverseIdxLength == myMesh3->getNumberOfNodes()+1);
01226 
01227   for (int i = 0; i < myMesh3->getNumberOfNodes(); i++) 
01228     {
01229       cout << "Node "<< i+1 << " : ";
01230       for (int j = ReverseNodalConnectivityIndex[i]; j < ReverseNodalConnectivityIndex[i+1]; j++)
01231         cout << ReverseNodalConnectivity[j-1] << " ";
01232       cout << endl;
01233     }
01234 
01235   // Show Descending Connectivity
01236   int NumberOfElements1;
01237   const int * connectivity;
01238   const int * connectivity_index;
01239   myMesh3->calculateConnectivity(MED_DESCENDING, MED_EN::MED_CELL);
01240   try
01241     {
01242       NumberOfElements1 = myMesh3->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
01243       connectivity = myMesh3->getConnectivity( MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS);
01244       connectivity_index = myMesh3->getConnectivityIndex(MED_DESCENDING, MED_CELL);
01245     }
01246   catch (MEDEXCEPTION m) 
01247     {
01248       CPPUNIT_FAIL(m.what());
01249     }
01250 
01251   for (int j = 0; j < NumberOfElements1; j++) 
01252     {
01253       cout << "Element " << j+1 << " : ";
01254       for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
01255         cout << connectivity[k-1] << " ";
01256       cout << endl;
01257     }
01258 
01259   //test 3D mesh
01260   for(int ind = SpaceDim; ind > 1; ind-- )
01261     {
01262       int NumberOfElem = myMesh3->getNumberOfElements (constituentEntity,MED_ALL_ELEMENTS);
01263       if(NumberOfElem < 1) continue;
01264 
01265       const SUPPORT * sup = myMesh3->getSupportOnAll( constituentEntity );
01266 
01267       if (ind == 2)
01268         {
01269           // test of normal(for 1d or 2d elements)
01270           FIELD<double>* normal;
01271           CPPUNIT_ASSERT_NO_THROW(normal = myMesh3->getNormal(sup));
01272 
01273           double normal_square, norm;
01274           double maxnorm=0.;
01275           double minnorm=0.;
01276           double tmp_value;
01277           for (int i = 1; i<=NumberOfElem; i++) 
01278             {
01279               normal_square = 0.;
01280               cout << "Normal " << i << " ";
01281               for (int j=1; j<=SpaceDim; j++) 
01282                 {
01283                   tmp_value = normal->getValueIJ(i,j);
01284                   normal_square += tmp_value*tmp_value;
01285                   cout << tmp_value << " ";
01286                 }
01287               norm = sqrt(normal_square);
01288               maxnorm = dmax(maxnorm,norm);
01289               minnorm = dmin(minnorm,norm);
01290               cout << ", Norm = " << norm << endl;
01291             }
01292           cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
01293           delete normal;
01294 
01295           // test of area(for 2d elements)
01296           FIELD<double>* area;
01297           CPPUNIT_ASSERT_NO_THROW(area = myMesh3->getArea(sup));
01298 
01299           double maxarea,minarea,areatot;
01300           maxarea = 0.;
01301           minarea = 0.;
01302           areatot = 0.0;
01303           for (int i = 1; i<=NumberOfElem;i++)
01304             {
01305               cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
01306               maxarea = dmax(maxarea,area->getValueIJ(i,1));
01307               minarea = dmin(minarea,area->getValueIJ(i,1));
01308               areatot = areatot + area->getValueIJ(i,1);
01309             }
01310 
01311           cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
01312           cout << "Support Area " << areatot << endl;
01313 
01314           delete area;
01315         }
01316 
01317       // test of barycenter(for 3d and 2d elements)
01318       FIELD<double>* barycenter;
01319       CPPUNIT_ASSERT_NO_THROW(barycenter = myMesh3->getBarycenter(sup));
01320 
01321       CPPUNIT_ASSERT_NO_THROW(NumberOfElem = myMesh3->getNumberOfElements(constituentEntity,MED_ALL_ELEMENTS));
01322 
01323       for (int i = 1; i<=NumberOfElem;i++)
01324         {
01325           if (ind == 3)
01326             cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
01327 
01328           if (ind == 2)
01329             cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << endl;
01330         }
01331       delete barycenter;
01332 
01333       // test of volume(for 3d elements)
01334       if (ind == 3)
01335         {
01336           FIELD<double>* volume;
01337           CPPUNIT_ASSERT_NO_THROW(volume= myMesh3->getVolume(sup));
01338 
01339           double maxvol,minvol,voltot;
01340           maxvol = 0.;
01341           minvol = 0.;
01342           voltot = 0.0;
01343           for (int i = 1; i<=NumberOfElem;i++)
01344             {
01345               cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
01346               maxvol = dmax(maxvol,volume->getValueIJ(i,1));
01347               minvol = dmin(minvol,volume->getValueIJ(i,1));
01348               voltot = voltot + volume->getValueIJ(i,1);
01349             }
01350 
01351           cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
01352           cout << "Support Volume " << voltot << endl;
01353 
01354           delete volume;
01355 
01356           // test of skin
01357           SUPPORT *skin;
01358           CPPUNIT_ASSERT_NO_THROW(skin = myMesh3->getSkin(sup));
01359 
01360           //test mergeSupports and intersectSupports. vactor contains only 1 elements
01361           vector<SUPPORT *> myVectSup;
01362           myVectSup.push_back(skin);
01363 
01364           //method return a copy of skin object
01365           SUPPORT *copyMergeSkin;
01366           CPPUNIT_ASSERT_NO_THROW(copyMergeSkin = myMesh3->mergeSupports(myVectSup));
01367           try
01368             {
01369               CPPUNIT_ASSERT(copyMergeSkin->deepCompare(*skin));
01370             }
01371           catch (const std::exception &e)
01372             {
01373               CPPUNIT_FAIL(e.what());
01374             }
01375           catch (...)
01376             {
01377               CPPUNIT_FAIL("Unknown exception");
01378             }
01379 
01380           //method return a copy of skin object
01381           SUPPORT *copyIntersectSkin;
01382           CPPUNIT_ASSERT_NO_THROW(copyIntersectSkin = myMesh3->intersectSupports(myVectSup));
01383           try
01384             {
01385               CPPUNIT_ASSERT(copyIntersectSkin->deepCompare(*skin));
01386             }
01387           catch (const std::exception &e)
01388             {
01389               CPPUNIT_FAIL(e.what());
01390             }
01391           catch (...)
01392             {
01393               CPPUNIT_FAIL("Unknown exception");
01394             }
01395 
01396           skin->removeReference();
01397           copyMergeSkin->removeReference();
01398           copyIntersectSkin->removeReference();
01399         }
01400 
01401       constituentEntity++;
01402     }
01403 
01404 
01405   // Testing length and normal vectors on 1d elements
01406   {
01407     // coordinates
01408     int NumberOfNodes3 = 4;
01409 
01410     string Names3[3] = 
01411       {
01412         "X","Y","Z"
01413       };
01414     string Units3[3] = 
01415       {
01416         "cm","cm","cm"
01417       };
01418 
01419     double Coordinates3[4*2] = 
01420       {
01421         0.0,  0.0,  // n1
01422         1.0,  1.0,  // n2
01423         0.0,  1.0,  // n3
01424         1.0,  0.0
01425       }; // n4
01426 
01427     const int NumberOfEdgeTypes = 1;
01428     MED_EN::medGeometryElement EdgeTypes[NumberOfEdgeTypes] = 
01429       {
01430         MED_SEG2
01431       };
01432     const int NumberOfEdges[NumberOfEdgeTypes] = 
01433       {
01434         4
01435       };
01436     int ConnectivityEdge[4*2] = 
01437       {
01438         1,2, 2,3, 3,4, 4,1
01439       };
01440 
01441     // CREATE THE MESH
01442     MEDMEM::MESHING* myMeshing3 = new MEDMEM::MESHING;
01443     myMeshing3->setName("meshing3");
01444     myMeshing3->setCoordinates(/*SpaceDimension*/2, NumberOfNodes3, Coordinates3,
01445                                "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
01446     myMeshing3->setCoordinatesNames(Names3);
01447     myMeshing3->setCoordinatesUnits(Units3);
01448 
01449     // define connectivities
01450     //      cell part
01451     const int NumberOfTypes3 = 1;
01452     medGeometryElement Types3[NumberOfTypes3] = 
01453       {
01454         MED_QUAD4
01455       };
01456     const int NumberOfElements3[NumberOfTypes3] = 
01457       {
01458         1
01459       };
01460 
01461     myMeshing3->setNumberOfTypes(NumberOfTypes3,MED_CELL);
01462     myMeshing3->setTypes(Types3,MED_CELL);
01463     myMeshing3->setNumberOfElements(NumberOfElements3,MED_CELL);
01464 
01465     int Connectivityquad[1*4] = 
01466       {
01467         1,2,3,4
01468       };
01469 
01470     myMeshing3->setConnectivity( MED_CELL, MED_QUAD4, Connectivityquad );
01471 
01472     myMeshing3->setNumberOfTypes(NumberOfEdgeTypes, MED_EDGE);
01473     myMeshing3->setTypes(EdgeTypes, MED_EDGE);
01474     myMeshing3->setNumberOfElements(NumberOfEdges, MED_EDGE);
01475 
01476     myMeshing3->setConnectivity( MED_EDGE, MED_SEG2, ConnectivityEdge );
01477 
01478     //test 2D mesh
01479     int NumberOfElem = myMeshing3->getNumberOfElements (MED_EDGE, MED_ALL_ELEMENTS);
01480 
01481     const SUPPORT * sup = myMeshing3->getSupportOnAll( MED_EDGE );
01482 
01483     // test of normal(for 1d or 2d elements)
01484     FIELD<double>* normal;
01485     CPPUNIT_ASSERT_NO_THROW(normal = myMeshing3->getNormal(sup));
01486 
01487     double normal_square, norm;
01488     double maxnorm=0.;
01489     double minnorm=0.;
01490     double tmp_value;
01491     for (int i = 1; i<=NumberOfElem; i++) 
01492       {
01493         normal_square = 0.;
01494         cout << "Normal " << i << " ";
01495         for (int j=1; j<=/*SpaceDimension*/2; j++) 
01496           {
01497             tmp_value = normal->getValueIJ(i,j);
01498             normal_square += tmp_value*tmp_value;
01499             cout << tmp_value << " ";
01500           }
01501         norm = sqrt(normal_square);
01502         maxnorm = dmax(maxnorm,norm);
01503         minnorm = dmin(minnorm,norm);
01504         cout << ", Norm = " << norm << endl;
01505       }
01506     cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
01507 
01508     // test of length(for 1d elements)
01509     FIELD<double>* length;
01510     CPPUNIT_ASSERT_NO_THROW(length = myMeshing3->getLength(sup));
01511 
01512     double length_value,maxlength,minlength;
01513     maxlength = 0;
01514     minlength = 0;
01515     for (int i = 1; i<=NumberOfElem;i++) 
01516       {
01517         length_value = length->getValueIJ(i,1);
01518         cout << "Length " << i << " " << length_value << endl;
01519         maxlength = dmax(maxlength,length_value);
01520         minlength = dmin(minlength,length_value);
01521       }
01522     cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
01523 
01524     vector< FIELD<double> *> myVectField1;
01525     myVectField1.push_back(normal);
01526     myVectField1.push_back(length);
01527     CPPUNIT_ASSERT_NO_THROW(myMeshing3->mergeFields(myVectField1));
01528 
01529     delete normal;
01530     delete length;
01531 
01532     //#ifdef ENABLE_FAULTS
01533     {
01534       // (BUG) Segmentation fault if vector is empty
01535       vector<SUPPORT *> myVectSupEmpty;
01536       CPPUNIT_ASSERT_THROW(myMesh3->mergeSupports(myVectSupEmpty), MEDEXCEPTION);
01537     }
01538     //#endif
01539 
01540     // test mergeFields method: Fields have the same value type
01541     //intersectSupports and mergeSupports methods
01542     {
01543       SUPPORT * sup1 = new SUPPORT;
01544       sup1->setMesh( myMeshing3 );
01545       sup1->setEntity( MED_EDGE );
01546       SUPPORT * sup2 = new SUPPORT;
01547       sup2->setMesh( myMeshing3 );
01548       sup2->setEntity( MED_EDGE );
01549       MED_EN::medGeometryElement gtEdges[1] = 
01550         {
01551           MED_SEG2
01552         };
01553       int nbEdges1[1] = 
01554         {
01555           1
01556         };
01557       int edges1[1] = 
01558         {
01559           1
01560         };
01561       int nbEdges2[1] = 
01562         {
01563           2
01564         };
01565       int edges2[2] = 
01566         {
01567           2,3
01568         };
01569       sup1->setpartial("description 1", 1, 1, gtEdges, nbEdges1, edges1);
01570       sup2->setpartial("description 1", 1, 2, gtEdges, nbEdges2, edges2);
01571 
01572       vector<SUPPORT *> myVectSup3;
01573       myVectSup3.push_back(sup1);
01574       myVectSup3.push_back(sup2);
01575       //method return a MergeSup on the union of all SUPPORTs in Supports.
01576       SUPPORT *MergeSup;
01577       CPPUNIT_ASSERT_NO_THROW(MergeSup = myMesh3->mergeSupports(myVectSup3));
01578       cout << *MergeSup << endl;
01579       MergeSup->removeReference();
01580 
01581       //method return a intersection of all SUPPORTs in IntersectSup
01582       SUPPORT *IntersectSup;
01583       CPPUNIT_ASSERT_NO_THROW(IntersectSup = myMesh3->intersectSupports(myVectSup3));
01584       if (IntersectSup != NULL) cout<< *IntersectSup <<endl;
01585       IntersectSup->removeReference();
01586 
01587       FIELD<double> * length1 = myMeshing3->getLength(sup1);
01588       FIELD<double> * length2 = myMeshing3->getLength(sup2);
01589 
01590       vector< FIELD<double> *> myVect12;
01591       myVect12.push_back(length1);
01592       myVect12.push_back(length2);
01593 
01594       FIELD<double> * length12;
01595       CPPUNIT_ASSERT_NO_THROW(length12 = myMeshing3->mergeFields(myVect12));
01596       delete length12;
01597 
01598       sup1->removeReference();
01599       sup2->removeReference();
01600       delete length1;
01601       delete length2;
01602     }
01603   }
01604 
01606   // TEST 4: test MESH constructed from file pointe.med //
01608   MESH * myMesh4 = new MESH();
01609   myMesh4->setName(meshname);
01610   MED_MESH_RDONLY_DRIVER myMeshDriver (filename, myMesh4);
01611   myMeshDriver.setMeshName(meshname);
01612 
01613   //Mesh has no driver->segmentation violation
01614   //CPPUNIT_ASSERT_THROW(myMesh4->read(), MEDEXCEPTION);
01615 
01616   //Add an existing MESH driver.
01617   int myDriver4;
01618   CPPUNIT_ASSERT_NO_THROW(myDriver4 = myMesh4->addDriver(myMeshDriver));
01619 
01620   //read all objects in the file
01621   CPPUNIT_ASSERT_NO_THROW(myMesh4->read(myDriver4));
01622 
01623   if (myMesh4->getIsAGrid()) 
01624     {
01625       GRID* myGrid = dynamic_cast<GRID*>(myMesh4);
01626       CPPUNIT_ASSERT(myGrid);
01627     }
01628 
01629   //myDriver4->DRONLY->can't write
01630   CPPUNIT_ASSERT_THROW(myMesh4->write(myDriver4), MEDEXCEPTION);
01631 
01632   // add new driver
01633   int idMeshV21;
01634   CPPUNIT_ASSERT_NO_THROW(idMeshV21 = myMesh4->addDriver(MED_DRIVER,filenameout21));
01635 
01636   //Write all the content of the MESH using driver referenced by the integer handler index.
01637   CPPUNIT_ASSERT_NO_THROW(myMesh4->write(idMeshV21));
01638 
01639   // remove driver from mesh
01640   CPPUNIT_ASSERT_NO_THROW(myMesh4->rmDriver(myDriver4));
01641   //#ifdef ENABLE_FORCED_FAILURES
01642   //  CPPUNIT_FAIL("ERROR: driver with index idMedV21 has not been removed");
01643   //#endif
01644   // ensure exception is raised on second attempt to remove driver
01645   //CPPUNIT_ASSERT_THROW(myMesh4->rmDriver(myDriver4),MEDEXCEPTION);
01646 
01647   // Create a MESH object using a MESH driver of type MED_DRIVER associated with file fileName.
01648   MESH* myMesh5;
01649   CPPUNIT_ASSERT_NO_THROW(myMesh5 = new MESH(MED_DRIVER, filename, meshname));
01650   if(myMesh5->getIsAGrid())
01651     {
01652       GRID* myGrid = dynamic_cast<GRID*>(myMesh4);
01653       CPPUNIT_ASSERT(myGrid);
01654     }
01655 
01656   //ensure two meshes constracted from one file in two different ways are equal
01657   CPPUNIT_ASSERT(myMesh5->deepCompare(*myMesh4));
01658 
01659   int myDriver6;
01660   MESH* myMesh6 = new MESH();
01661   try
01662     {
01663       myDriver6 = myMesh6->addDriver(MED_DRIVER, filename, meshname, RDONLY);
01664     }
01665   catch (const std::exception &e)
01666     {
01667       CPPUNIT_FAIL(e.what());
01668     }
01669   catch (...)
01670     {
01671       CPPUNIT_FAIL("Unknown exception");
01672     }
01673 
01674   try
01675     {
01676       myMesh6->read(myDriver6);
01677     }
01678   catch (const std::exception &e)
01679     {
01680       CPPUNIT_FAIL(e.what());
01681     }
01682   catch (...)
01683     {
01684       CPPUNIT_FAIL("Unknown exception");
01685     }
01686 
01687   //ensure two meshes constracted from one file in two different ways are equal
01688   CPPUNIT_ASSERT(myMesh6->deepCompare(*myMesh4));
01689 
01690   //test FAMILY
01691   int NumberOfFamilies4;
01692   CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies4 = myMesh6->getNumberOfFamilies(MED_CELL));
01693   CPPUNIT_ASSERT_MESSAGE("Current mesh hasn't Families", NumberOfFamilies4 != 0);
01694 
01695   vector<FAMILY*> families4;
01696   CPPUNIT_ASSERT_NO_THROW(families4 = myMesh6->getFamilies(MED_CELL));
01697   CPPUNIT_ASSERT(families4.size() == NumberOfFamilies4);
01698   for(int nb = 1; nb <= NumberOfFamilies4; nb++ )
01699     {
01700       const FAMILY* family;
01701       CPPUNIT_ASSERT_NO_THROW(family = myMesh6->getFamily(MED_CELL, nb));
01702       CPPUNIT_ASSERT_EQUAL(family->getName(), families4[nb-1]->getName());
01703     }
01704 
01705   //get support which reference all elements on the boundary of mesh.
01706   SUPPORT*myBndSup;
01707   CPPUNIT_ASSERT_THROW(myMesh6->getBoundaryElements(MED_CELL), MEDEXCEPTION);
01708   //get only face in 3D.
01709   CPPUNIT_ASSERT_NO_THROW(myBndSup = myMesh6->getBoundaryElements(MED_FACE));
01710 
01711   //test buildSupportOnElementsFromElementList and buildSupportOnNodeFromElementList
01712   const int * myConnectivityValue6;
01713   CPPUNIT_ASSERT_NO_THROW(myConnectivityValue6 = myMesh6->getReverseConnectivity(MED_DESCENDING));
01714   const int * myConnectivityIndex6;
01715   CPPUNIT_ASSERT_NO_THROW(myConnectivityIndex6 = myMesh6->getReverseConnectivityIndex(MED_DESCENDING));
01716   int numberOfElem6;
01717   CPPUNIT_ASSERT_NO_THROW(numberOfElem6 = myMesh6->getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS));
01718   list<int> myElementsList6;
01719 
01720   for (int i=0; i<numberOfElem6; i++)
01721     if (myConnectivityValue6[myConnectivityIndex6[i]] == 0) 
01722       {
01723         myElementsList6.push_back(i+1);
01724       }
01725 
01726   SUPPORT * mySupportOnNode;
01727   SUPPORT * mySupportOnElem;
01728   CPPUNIT_ASSERT_NO_THROW(mySupportOnElem = myMesh6->buildSupportOnElementsFromElementList(myElementsList6,MED_FACE));
01729   CPPUNIT_ASSERT(mySupportOnElem->deepCompare(*myBndSup));
01730   CPPUNIT_ASSERT_EQUAL(MED_FACE, mySupportOnElem->getEntity());
01731 
01732   list<int>::const_iterator iteronelem = myElementsList6.begin();
01733   for (int i = 1; i <= 3; i++, iteronelem++) 
01734     {
01735       CPPUNIT_ASSERT_EQUAL(i, mySupportOnElem->getValIndFromGlobalNumber(*iteronelem));
01736     }
01737 
01738   CPPUNIT_ASSERT_NO_THROW(mySupportOnNode = myMesh6->buildSupportOnNodeFromElementList(myElementsList6,MED_FACE));
01739   CPPUNIT_ASSERT(mySupportOnNode->deepCompare( *(myMesh6->getBoundaryElements(MED_NODE))));
01740 
01741   //sets mesh fields to initial values
01742   try
01743     {
01744       myMesh6->init();
01745     }
01746   catch (const std::exception &e)
01747     {
01748       CPPUNIT_FAIL(e.what());
01749     }
01750   catch (...)
01751     {
01752       CPPUNIT_FAIL("Unknown exception");
01753     }
01754 
01755   //ensure two meshes constracted from one file in two different ways are equal
01756   CPPUNIT_ASSERT(!myMesh6->deepCompare(*myMesh4));
01757 
01758   //ensure mesh is empty
01759   CPPUNIT_ASSERT(myMesh6->getSpaceDimension() == MED_INVALID);
01760   CPPUNIT_ASSERT(myMesh6->getNumberOfNodes() == MED_INVALID);
01761   CPPUNIT_ASSERT(myMesh6->getCoordinateptr() == NULL);
01762 
01763   delete myMesh4;
01764   delete myMesh5;
01765   delete myMesh6;
01766 
01767   MESH* myMesh7 = MEDMEMTest_createTestMesh();
01768   vector< vector<double> > myBndBox;
01769   try
01770     {
01771       myBndBox = myMesh7->getBoundingBox();
01772     }
01773   catch (const std::exception &e)
01774     {
01775       CPPUNIT_FAIL(e.what());
01776     }
01777   catch (...)
01778     {
01779       CPPUNIT_FAIL("Unknown exception");
01780     }
01781 
01782   cout<<"Bounding box for createTestMesh()"<<endl;
01783   for(int i = 0; i < myBndBox.size(); i++)
01784     {
01785       for(int j = 0; j < myBndBox[i].size(); j++)
01786         cout<<" "<< myBndBox[i][j]<<" ";
01787       cout<<endl;
01788     }
01789 
01790   double CoorPoint[3] = 
01791     {
01792       0.0,  0.0, 1.0
01793     }; //n2
01794   int idxElem;
01795   try
01796     {
01797       idxElem = myMesh7->getElementContainingPoint(CoorPoint);
01798     }
01799   catch (const std::exception &e)
01800     {
01801       CPPUNIT_FAIL(e.what());
01802     }
01803   catch (...)
01804     {
01805       CPPUNIT_FAIL("Unknown exception");
01806     }
01807   CPPUNIT_ASSERT(idxElem != -1);
01808 
01809   double CoorNoPoint[3] = 
01810     {
01811       5.0,  0.0, -5.0
01812     }; //there is no such point
01813   int idxNoElem;
01814   try
01815     {
01816       idxNoElem = myMesh7->getElementContainingPoint(CoorNoPoint);
01817     }
01818   catch (const std::exception &e)
01819     {
01820       CPPUNIT_FAIL(e.what());
01821     }
01822   catch (...)
01823     {
01824       CPPUNIT_FAIL("Unknown exception");
01825     }
01826   CPPUNIT_ASSERT(idxNoElem == -1);
01827 }
01828 
01829 int main (int argc, char** argv)
01830 {
01831   MEDMEMTest_testMeshAndMeshing();
01832 }