src/MEDMEMBinTest/test_MEDMEM_poly3D.cxx

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