src/MEDMEM/test_MEDMEM_poly3D.cxx

   1 // This program test the implementation of the class MED_MESH_RDONLY_DRIVER for Polygons and Polyhedron. //
   2
   3 #include "MEDMEM_define.hxx"
   4 #include "MEDMEM_Mesh.hxx"
   5 #include "MEDMEM_Support.hxx"
   6 #include "MEDMEM_Family.hxx"
   7 #include "MEDMEM_Field.hxx"
   8 #include "MEDMEM_ModulusArray.hxx"
   9 #include "MEDMEM_MedMeshDriver.hxx"
  10 #include <vector>
  11
  12 #define MESHNAME "poly3D"
  13
  14 using namespace std;
  15 using namespace MEDMEM;
  16 using namespace MED_EN;
  17
  18 #define DIM_OF_FIELD 3
  19
  20
  21 class SupportTester {
  22 private:
  23   const int *_tabOfNodes;
  24   vector<int> _eltsActiveYet;
  25   vector<int> _lgthOfEltsActiveYet;
  26 public:
  27   SupportTester(const int *tabOfNodes, int nbOfElts, int nbOfNodesPerElt);
  28   SupportTester(const int *tabOfNodes, int nbOfElts, const int *nbOfNodesPerElt);
  29   bool isIncludedAndNotAlreadyConsumed(const int *tabOfNodesOfTheElementToTest);
  30   bool areAllEltsConsumed();
  31 private:
  32   static bool areEquivalent(const int *nodes1, const int *nodes2, int nbOfNodes);
  33 };
  34
  35 SupportTester::SupportTester(const int *tabOfNodes, int nbOfElts, int nbOfNodesPerElt):_tabOfNodes(tabOfNodes)
  36 {
  37   for(int i=0;i<nbOfElts;i++)
  38     {
  39       _eltsActiveYet.push_back(i*nbOfNodesPerElt);
  40       _lgthOfEltsActiveYet.push_back(nbOfNodesPerElt);
  41     }
  42 }
  43
  44 SupportTester::SupportTester(const int *tabOfNodes, int nbOfElts, const int *nbOfNodesPerElt):_tabOfNodes(tabOfNodes)
  45 {
  46   int offset=0;
  47   for(int i=0;i<nbOfElts;i++)
  48     {
  49       _eltsActiveYet.push_back(offset);
  50       _lgthOfEltsActiveYet.push_back(nbOfNodesPerElt[i]);
  51     }
  52 }
  53
  54 bool SupportTester::isIncludedAndNotAlreadyConsumed(const int *tabOfNodesOfTheElementToTest)
  55 {
  56   vector<int>::iterator iter2=_lgthOfEltsActiveYet.begin();
  57   for(vector<int>::iterator iter=_eltsActiveYet.begin();iter!=_eltsActiveYet.end();iter++)
  58     if(areEquivalent(_tabOfNodes+(*iter),tabOfNodesOfTheElementToTest,*iter2))
  59       {
  60         _eltsActiveYet.erase(iter);
  61         _lgthOfEltsActiveYet.erase(iter2);
  62         return true;
  63       }
  64     else
  65       {
  66         iter2++;
  67       }
  68   return false;
  69 }
  70
  71 bool SupportTester::areAllEltsConsumed()
  72 {
  73   return _eltsActiveYet.size()==0;
  74 }
  75
  76 bool SupportTester::areEquivalent(const int *nodes1, const int *nodes2, int nbOfNodes)
  77 {
  78   MEDMODULUSARRAY arr1(nbOfNodes,nodes1);
  79   MEDMODULUSARRAY arr2(nbOfNodes,nodes2);
  80   return arr1.compare(arr2)!=0;
  81 }
  82
  83 int main (int argc, char ** argv)
  84 {
  85   if (argc<2) // after 2, ignored !
  86     {
  87       cerr << "Usage : " << argv[0] << " poly3D.med typically in ../../share/salome/resources" << endl << endl;
  88       exit(-1);
  89     }
  90   int nbOfPtsForTest=0;
  91   int i;
  92   string filename = argv[1];
  93   MESH * myMesh = new MESH;
  94   myMesh->setName(MESHNAME);
  95   MED_MESH_RDONLY_DRIVER myMeshReadDriver(filename,myMesh);
  96   myMeshReadDriver.setMeshName(MESHNAME);
  97   myMeshReadDriver.open();
  98   myMeshReadDriver.read();
  99   myMeshReadDriver.close();
 100   //Test 1 : test if connectivity of poly3D mesh is OK
 101   if(myMesh->getMeshDimension()==3 && myMesh->getSpaceDimension()==3)
 102     nbOfPtsForTest++;
 103   if(myMesh->getNumberOfElementsWithPoly(MED_EN::MED_CELL,MED_EN::MED_TETRA4)==1 && myMesh->getNumberOfElementsWithPoly(MED_EN::MED_CELL,MED_EN::MED_POLYHEDRA)==2)
 104     nbOfPtsForTest++;
 105   const int REFnodalConnForTetra[4]={17, 9, 18, 19};
 106   const int *connForTetraToTest=myMesh->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_TETRA4);
 107   const int *connIndForTetraToTest=myMesh->getConnectivityIndex(MED_NODAL,MED_CELL);
 108   for(i=connIndForTetraToTest[0]-1;i<connIndForTetraToTest[1]-1;i++)
 109     if(connForTetraToTest[i]==REFnodalConnForTetra[i])
 110       nbOfPtsForTest++;
 111   //6
 112   const int *globIndex=myMesh->getPolyhedronIndex(MED_NODAL);
 113   const int *nodalConnOfFaces=myMesh->getPolyhedronConnectivity(MED_NODAL);
 114   const int *facesIndex=myMesh->getPolyhedronFacesIndex();
 115   if(globIndex[1]-globIndex[0]==9 && globIndex[2]-globIndex[1]==10)// resp 9 faces and 10 faces are in polyh 1 and 2.
 116     nbOfPtsForTest++;
 117   //7
 118   const int REFfacesIndex[20]={1, 7, 11, 15, 19, 23, 27, 31, 34, 39, 44, 48, 52, 55, 58, 61, 64, 68, 72, 75};
 119   const int REFnodalConnOfFaces[74]={1, 2, 3, 4, 5, 6, // Polyhedron 1
 120                                      1, 7, 8, 2,
 121                                      2, 8, 9, 3,
 122                                      4, 3, 9, 10,
 123                                      5, 4, 10, 11,
 124                                      6, 5, 11, 12,
 125                                      1, 6, 12, 7,
 126                                      7, 12, 8, 10,
 127                                      9, 8, 12, 11,
 128
 129                                      13, 14, 15, 3, 2, // Polyhedron 2
 130                                      13, 2, 8, 16,
 131                                      14, 13, 16, 17,
 132                                      15, 14, 17, 15,
 133                                      17, 18, 15,
 134                                      18, 9, 3,
 135                                      15, 9, 2,
 136                                      3, 9, 8,
 137                                      8, 9, 17, 16,
 138                                      9, 18, 17 };
 139   for(i=0;i<20;i++)
 140     if(REFfacesIndex[i]==facesIndex[i])
 141       nbOfPtsForTest++;
 142   for(i=0;i<74;i++)
 143     if(REFnodalConnOfFaces[i]==nodalConnOfFaces[i])
 144       nbOfPtsForTest++;
 145   if(nbOfPtsForTest!=101)
 146     {
 147       cout << "TEST1 K0 ! : Invalid Globaldata in memory..." << endl;
 148       return 1;
 149     }
 150   // TEST 2 : FAMILY
 151   nbOfPtsForTest=0;
 152   vector<FAMILY*> families=myMesh->getFamilies(MED_FACE);
 153   if(families.size()==3)
 154     nbOfPtsForTest++;
 155   vector<FAMILY *>::iterator iter=families.begin();
 156   FAMILY *fam1=*(iter++);
 157   FAMILY *fam2=*(iter++);
 158   FAMILY *fam3=*(iter);
 159   const int *nbs;
 160   // family 1
 161   if(fam1->getNumberOfTypes()==1 && fam1->getTypes()[0]==MED_POLYGON && fam1->getNumberOfElements(MED_ALL_ELEMENTS)==3)
 162     nbOfPtsForTest++;
 163   nbs=fam1->getNumber(MED_ALL_ELEMENTS);
 164   const int REFTabForPolyg[16]={1, 2, 3, 4, 5, 6, 10, 9, 8, 12, 11, 13, 14, 15, 3, 2};
 165   const int REFTabForPolygLgth[3]={6,5,5};
 166   SupportTester test1(REFTabForPolyg,3,REFTabForPolygLgth);
 167   for(i=0;i<3;i++)
 168     {
 169       int lgth;
 170       const int *conn=((CONNECTIVITY *)myMesh->getConnectivityptr())->getConnectivityOfAnElementWithPoly(MED_NODAL,MED_FACE,nbs[i],lgth);
 171       if(test1.isIncludedAndNotAlreadyConsumed(conn))
 172         nbOfPtsForTest++;
 173     }
 174   if(test1.areAllEltsConsumed())
 175     nbOfPtsForTest++;
 176   // family 2
 177   if(fam2->getNumberOfElements(MED_ALL_ELEMENTS)==8)
 178     nbOfPtsForTest++;
 179   nbs=fam2->getNumber(MED_ALL_ELEMENTS);
 180   const int REFTabForQuad[32]={1, 7, 8, 2, 2, 8, 9, 3, 4, 3, 9, 10, 5, 4, 10, 11, 6, 5, 11, 12, 1, 6, 12, 7, 14, 13, 16, 17, 8, 9, 17, 16};
 181   SupportTester test2(REFTabForQuad,8,4);
 182   for(i=0;i<8;i++)
 183     {
 184       int lgth;
 185       const int *conn=((CONNECTIVITY *)myMesh->getConnectivityptr())->getConnectivityOfAnElementWithPoly(MED_NODAL,MED_FACE,nbs[i],lgth);
 186       if(test2.isIncludedAndNotAlreadyConsumed(conn))
 187         nbOfPtsForTest++;
 188     }
 189   if(test2.areAllEltsConsumed())
 190     nbOfPtsForTest++;
 191   // family 3
 192   if(fam3->getNumberOfElements(MED_ALL_ELEMENTS)==6)
 193     nbOfPtsForTest++;
 194   nbs=fam3->getNumber(MED_ALL_ELEMENTS);
 195   const int REFTabForTria[18]={7, 12, 8, 15, 14, 17, 15, 17, 18, 15, 18, 9, 3, 15, 9, 18, 17, 9};
 196   SupportTester test3(REFTabForTria,6,3);
 197   for(i=0;i<6;i++)
 198     {
 199       int lgth;
 200       const int *conn=((CONNECTIVITY *)myMesh->getConnectivityptr())->getConnectivityOfAnElementWithPoly(MED_NODAL,MED_FACE,nbs[i],lgth);
 201       if(test3.isIncludedAndNotAlreadyConsumed(conn))
 202         nbOfPtsForTest++;
 203     }
 204   if(test3.areAllEltsConsumed())
 205     nbOfPtsForTest++;
 206   if(nbOfPtsForTest!=21)
 207     {
 208       cout << "TEST2 K0 ! : Invalid data in memory for families !!!" << endl;
 209       return 1;
 210     }
 211   // TEST 3 : volumes, aeras, barycenter
 212   nbOfPtsForTest=0;
 213   SUPPORT *supOnCell=new SUPPORT(myMesh);
 214   FIELD<double>* vol1=myMesh->getVolume(supOnCell);
 215   int lgth=vol1->getValueLength(MED_FULL_INTERLACE);
 216   const double *vals=vol1->getValue(MED_FULL_INTERLACE);
 217   if(lgth==3)
 218     nbOfPtsForTest++;
 219   const double REFVolOfPolyHedron[3]={2.333333333333333,11.66666666666666,13.83224131414673};
 220   for(i=0;i<3;i++)
 221     if(fabs(REFVolOfPolyHedron[i]-vals[i])<1e-12)
 222       nbOfPtsForTest++;
 223   delete vol1;
 224   FIELD<double>* bary=myMesh->getBarycenter(supOnCell);
 225   lgth=bary->getValueLength(MED_FULL_INTERLACE);
 226   vals=bary->getValue(MED_FULL_INTERLACE);
 227   if(lgth==9)
 228     nbOfPtsForTest++;
 229   const double REFBaryOfPolyHedron[9]= {5.5, 1, -1, 2, 1.5, 1.0833333333333333, 5.1, 1.6, 0.9};
 230   for(i=0;i<9;i++)
 231     if(fabs(REFBaryOfPolyHedron[i]-vals[i])<1e-12)
 232       nbOfPtsForTest++;
 233   delete bary;
 234   delete supOnCell;
 235   //area
 236   vol1=myMesh->getArea(fam1);
 237   lgth=vol1->getValueLength(MED_FULL_INTERLACE);
 238   vals=vol1->getValue(MED_FULL_INTERLACE);
 239   if(lgth==3)
 240     nbOfPtsForTest++;
 241   const double REFAreaForPolyg[3]={6,5,6.5};
 242   for(i=0;i<3;i++)
 243     if(fabs(REFAreaForPolyg[i]-vals[i])<1e-12)
 244       nbOfPtsForTest++;
 245
 246   vol1=myMesh->getArea(fam2);
 247   lgth=vol1->getValueLength(MED_FULL_INTERLACE);
 248   vals=vol1->getValue(MED_FULL_INTERLACE);
 249   if(lgth==8)
 250     nbOfPtsForTest++;
 251   const double REFAreaForQuad[8]={2.1213203435596424, 2.8284271247461903, 4.4721359549995796, 4.4721359549995796,
 252                                   2.8284271247461903, 2.1213203435596428, 3.6798724963767362, 4};
 253   for(i=0;i<8;i++)
 254     if(fabs(REFAreaForQuad[i]-vals[i])<1e-12)
 255       nbOfPtsForTest++;
 256   delete vol1;
 257
 258   vol1=myMesh->getArea(fam3);
 259   lgth=vol1->getValueLength(MED_FULL_INTERLACE);
 260   vals=vol1->getValue(MED_FULL_INTERLACE);
 261   if(lgth==6)
 262     nbOfPtsForTest++;
 263   const double REFAreaForTri[6]={2.9580398915498081, 1.4142135623730951, 2.2360679774997898,
 264                                  3.3541019662496847, 3.3541019662496847, 2.2360679774997898};
 265   for(i=0;i<6;i++)
 266     if(fabs(REFAreaForTri[i]-vals[i])<1e-12)
 267       nbOfPtsForTest++;
 268   delete vol1;
 269   if(nbOfPtsForTest!=34)
 270     {
 271       cout << "TEST3 K0 ! : Error in caluclation of basic properties !!!" << endl;
 272       return 1;
 273     }
 274   // TEST 4 -- CHECK FOR Reverse descending using getBoundaryElements.
 275   nbOfPtsForTest=0;
 276   SUPPORT *bound=myMesh->getBoundaryElements(MED_NODE);
 277   if(bound->getNumberOfElements(MED_ALL_ELEMENTS)==19)
 278     nbOfPtsForTest++;
 279   if(bound->isOnAllElements())
 280       nbOfPtsForTest++;
 281   if(nbOfPtsForTest!=2)
 282     {
 283       cout << "TEST4 K0 ! : Error in getBoundaryElements probably due to Reverse descending !!!" << endl;
 284       return 1;
 285     }
 286   delete bound;
 287   ///////////
 288   cout << "ALL TESTS OK !!!" << endl;
 289   delete myMesh;
 290   return 0;
 291 }