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