src/MEDMEM/test_MEDMEM_poly3D.cxx

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