Back to index

salome-med  6.5.0
test_MEDMEM_poly3D.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 #include "MEDMEM_define.hxx"
00023 #include "MEDMEM_Mesh.hxx"
00024 #include "MEDMEM_Support.hxx"
00025 #include "MEDMEM_Family.hxx"
00026 #include "MEDMEM_Field.hxx"
00027 #include "MEDMEM_ModulusArray.hxx"
00028 #include "MEDMEM_MedMeshDriver.hxx"
00029 #include <vector>
00030 
00031 #define MESHNAME "poly3D"
00032 
00033 using namespace std;
00034 using namespace MEDMEM;
00035 using namespace MED_EN;
00036 
00037 #define DIM_OF_FIELD 3
00038 
00039 
00040 class SupportTester 
00041 {
00042 private:
00043   const int *_tabOfNodes;
00044   vector<int> _eltsActiveYet;
00045   vector<int> _lgthOfEltsActiveYet;
00046 public:
00047   SupportTester(const int *tabOfNodes, int nbOfElts, int nbOfNodesPerElt);
00048   SupportTester(const int *tabOfNodes, int nbOfElts, const int *nbOfNodesPerElt);
00049   bool isIncludedAndNotAlreadyConsumed(const int *tabOfNodesOfTheElementToTest);
00050   bool areAllEltsConsumed();
00051 private:
00052   static bool areEquivalent(const int *nodes1, const int *nodes2, int nbOfNodes);
00053 };
00054 
00055 SupportTester::SupportTester(const int *tabOfNodes, int nbOfElts, int nbOfNodesPerElt):_tabOfNodes(tabOfNodes)
00056 {
00057   for(int i=0;i<nbOfElts;i++)
00058     {
00059       _eltsActiveYet.push_back(i*nbOfNodesPerElt);
00060       _lgthOfEltsActiveYet.push_back(nbOfNodesPerElt);
00061     }
00062 }
00063 
00064 SupportTester::SupportTester(const int *tabOfNodes, int nbOfElts, const int *nbOfNodesPerElt):_tabOfNodes(tabOfNodes)
00065 {
00066   int offset=0;
00067   for(int i=0;i<nbOfElts;i++)
00068     {
00069       _eltsActiveYet.push_back(offset);
00070       _lgthOfEltsActiveYet.push_back(nbOfNodesPerElt[i]);
00071     }
00072 }
00073 
00074 bool SupportTester::isIncludedAndNotAlreadyConsumed(const int *tabOfNodesOfTheElementToTest)
00075 {
00076   vector<int>::iterator iter2=_lgthOfEltsActiveYet.begin();
00077   for(vector<int>::iterator iter=_eltsActiveYet.begin();iter!=_eltsActiveYet.end();iter++)
00078     if(areEquivalent(_tabOfNodes+(*iter),tabOfNodesOfTheElementToTest,*iter2))
00079       {
00080         _eltsActiveYet.erase(iter);
00081         _lgthOfEltsActiveYet.erase(iter2);
00082         return true;
00083       }
00084     else
00085       {
00086         iter2++;
00087       }
00088   return false;
00089 }
00090 
00091 bool SupportTester::areAllEltsConsumed()
00092 {
00093   return _eltsActiveYet.size()==0;
00094 }
00095 
00096 bool SupportTester::areEquivalent(const int *nodes1, const int *nodes2, int nbOfNodes)
00097 {
00098   MEDMODULUSARRAY arr1(nbOfNodes,nodes1);
00099   MEDMODULUSARRAY arr2(nbOfNodes,nodes2);
00100   return arr1.compare(arr2)!=0;
00101 }
00102 
00103 int main (int argc, char ** argv)
00104 {
00105   if (argc<2) // after 2, ignored !
00106     {
00107       cerr << "Usage : " << argv[0] << " poly3D.med typically in ../../share/salome/resources/med" << endl << endl;
00108       exit(-1);
00109     }
00110   int nbOfPtsForTest=0;
00111   int i;
00112   string filename = argv[1];
00113   MESH * myMesh = new MESH;
00114   myMesh->setName(MESHNAME);
00115   MED_MESH_RDONLY_DRIVER myMeshReadDriver(filename,myMesh);
00116   myMeshReadDriver.setMeshName(MESHNAME);
00117   myMeshReadDriver.open();
00118   myMeshReadDriver.read();
00119   myMeshReadDriver.close();
00120   //Test 1 : test if connectivity of poly3D mesh is OK
00121   if(myMesh->getMeshDimension()==3 && myMesh->getSpaceDimension()==3)
00122     nbOfPtsForTest++;
00123   if(myMesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_TETRA4)==1 && myMesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_POLYHEDRA)==2)
00124     nbOfPtsForTest++;
00125   const int REFnodalConnForTetra[4]=
00126     {
00127       17, 9, 18, 19
00128     };
00129   const int *connForTetraToTest=myMesh->getConnectivity(MED_NODAL,MED_CELL,MED_TETRA4);
00130   const int *connIndForTetraToTest=myMesh->getConnectivityIndex(MED_NODAL,MED_CELL);
00131   for(i=connIndForTetraToTest[0]-1;i<connIndForTetraToTest[1]-1;i++)
00132     if(connForTetraToTest[i]==REFnodalConnForTetra[i])
00133       nbOfPtsForTest++;
00134   //6
00135   const int *globIndex=connIndForTetraToTest + 1; // skip 1 tetra
00136   const int *nodalConnOfFaces=myMesh->getConnectivity(MED_NODAL,MED_CELL,MED_POLYHEDRA);
00137   if(globIndex[1]-globIndex[0]==46 && globIndex[2]-globIndex[1]==45)// resp 46 nodes and 45 nodes are in polyh 1 and 2.
00138     nbOfPtsForTest++;
00139   //7
00140   const int REFnodalConnOfFaces[91]=
00141     {
00142       1, 2, 3, 4, 5, 6, -1,// Polyhedron 1
00143       1, 7, 8, 2,       -1,
00144       2, 8, 9, 3,       -1,
00145       4, 3, 9, 10,      -1,
00146       5, 4, 10, 11,     -1,
00147       6, 5, 11, 12,     -1,
00148       1, 6, 12, 7,      -1,
00149       7, 12, 8,         -1,
00150       10, 9, 8, 12, 11,     
00151 
00152       13, 14, 15, 3, 2, -1,// Polyhedron 2
00153       13, 2, 8, 16,     -1,
00154       14, 13, 16, 17,   -1,
00155       15, 14, 17,       -1,
00156       15, 17, 18,       -1,
00157       15, 18, 9,        -1,
00158       3, 15, 9,         -1,
00159       2, 3, 9, 8,       -1,
00160       8, 9, 17, 16,     -1,
00161       9, 18, 17 
00162     };
00163   for(i=0;i<74;i++)
00164     if(REFnodalConnOfFaces[i]==nodalConnOfFaces[i])
00165       nbOfPtsForTest++;
00166   if(nbOfPtsForTest!=7+74)
00167     {
00168       cout << "TEST1 K0 ! : Invalid Globaldata in memory..." << endl;
00169       return 1;
00170     }
00171   // TEST 2 : FAMILY 
00172   nbOfPtsForTest=0;
00173   vector<FAMILY*> families=myMesh->getFamilies(MED_FACE);
00174   if(families.size()==3)
00175     nbOfPtsForTest++;
00176   vector<FAMILY *>::iterator iter=families.begin();
00177   FAMILY *fam1=*(iter++);
00178   FAMILY *fam2=*(iter++);
00179   FAMILY *fam3=*(iter);
00180   const int *nbs;
00181   // family 1
00182   if(fam1->getNumberOfTypes()==1 && fam1->getTypes()[0]==MED_POLYGON && fam1->getNumberOfElements(MED_ALL_ELEMENTS)==3)
00183     nbOfPtsForTest++;
00184   nbs=fam1->getNumber(MED_ALL_ELEMENTS);
00185   const int REFTabForPolyg[16]=
00186     {
00187       1, 2, 3, 4, 5, 6, 10, 9, 8, 12, 11, 13, 14, 15, 3, 2
00188     };
00189   const int REFTabForPolygLgth[3]=
00190     {
00191       6,5,5
00192     };
00193   SupportTester test1(REFTabForPolyg,3,REFTabForPolygLgth);
00194   for(i=0;i<3;i++)
00195     {
00196       int lgth;
00197       const int *conn=((CONNECTIVITY *)myMesh->getConnectivityptr())->getConnectivityOfAnElement(MED_NODAL,MED_FACE,nbs[i],lgth);
00198       if(test1.isIncludedAndNotAlreadyConsumed(conn))
00199         nbOfPtsForTest++;
00200     }
00201   if(test1.areAllEltsConsumed())
00202     nbOfPtsForTest++;
00203   // family 2
00204   if(fam2->getNumberOfElements(MED_ALL_ELEMENTS)==8)
00205     nbOfPtsForTest++;
00206   nbs=fam2->getNumber(MED_ALL_ELEMENTS);
00207   const int REFTabForQuad[32]=
00208     {
00209       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
00210     };
00211   SupportTester test2(REFTabForQuad,8,4);
00212   for(i=0;i<8;i++)
00213     {
00214       int lgth;
00215       const int *conn=((CONNECTIVITY *)myMesh->getConnectivityptr())->getConnectivityOfAnElement(MED_NODAL,MED_FACE,nbs[i],lgth);
00216       if(test2.isIncludedAndNotAlreadyConsumed(conn))
00217         nbOfPtsForTest++;
00218     }
00219   if(test2.areAllEltsConsumed())
00220     nbOfPtsForTest++;
00221   // family 3
00222   if(fam3->getNumberOfElements(MED_ALL_ELEMENTS)==6)
00223     nbOfPtsForTest++;
00224   nbs=fam3->getNumber(MED_ALL_ELEMENTS);
00225   const int REFTabForTria[18]=
00226     {
00227       7, 12, 8, 15, 14, 17, 15, 17, 18, 15, 18, 9, 3, 15, 9, 18, 17, 9
00228     };
00229   SupportTester test3(REFTabForTria,6,3);
00230   for(i=0;i<6;i++)
00231     {
00232       int lgth;
00233       const int *conn=((CONNECTIVITY *)myMesh->getConnectivityptr())->getConnectivityOfAnElement(MED_NODAL,MED_FACE,nbs[i],lgth);
00234       if(test3.isIncludedAndNotAlreadyConsumed(conn))
00235         nbOfPtsForTest++;
00236     }
00237   if(test3.areAllEltsConsumed())
00238     nbOfPtsForTest++;
00239   if(nbOfPtsForTest!=21)
00240     {
00241       cout << "TEST2 K0 ! : Invalid data in memory for families !!!" << endl;
00242       return 1;
00243     }
00244   // TEST 3 : volumes, areas, barycenter
00245   nbOfPtsForTest=0;
00246   const SUPPORT *supOnCell=myMesh->getSupportOnAll(MED_CELL);
00247   FIELD<double>* vol1=myMesh->getVolume(supOnCell, false);
00248   int lgth=vol1->getValueLength();
00249   const double *vals=vol1->getValue();
00250   if(lgth==3)
00251     nbOfPtsForTest++;
00252   const double REFVolOfPolyHedron[3]=
00253     {
00254       2.333333333333333,-11.66666666666666,-13.83224131414673
00255     };
00256   for(i=0;i<3;i++)
00257     if(fabs(REFVolOfPolyHedron[i]-vals[i])<1e-12)
00258       nbOfPtsForTest++;
00259   vol1->removeReference();
00260   vol1=myMesh->getVolume(supOnCell, true);
00261   lgth=vol1->getValueLength();
00262   vals=vol1->getValue();
00263   if(lgth==3)
00264     nbOfPtsForTest++;
00265   for(i=0;i<3;i++)
00266     if(fabs(fabs(REFVolOfPolyHedron[i])-vals[i])<1e-12)
00267       nbOfPtsForTest++;
00268   vol1->removeReference();
00269   FIELD<double>* bary=myMesh->getBarycenter(supOnCell);
00270   lgth=bary->getValueLength();
00271   vals=bary->getValue();
00272   if(lgth==9)
00273     nbOfPtsForTest++;
00274   const double REFBaryOfPolyHedron[9]= 
00275     {
00276       5.5, 1, -1, 2, 1.5, 1.0833333333333333, 5.1, 1.6, 0.9
00277     };
00278   for(i=0;i<9;i++)
00279     if(fabs(REFBaryOfPolyHedron[i]-vals[i])<1e-12)
00280       nbOfPtsForTest++;
00281   bary->removeReference();
00282   //area
00283   vol1=myMesh->getArea(fam1);
00284   lgth=vol1->getValueLength();
00285   vals=vol1->getValue();
00286   if(lgth==3)
00287     nbOfPtsForTest++;
00288   const double REFAreaForPolyg[3]=
00289     {
00290       6,5,6.5
00291     };
00292   for(i=0;i<3;i++)
00293     if(fabs(REFAreaForPolyg[i]-vals[i])<1e-12)
00294       nbOfPtsForTest++;
00295   vol1->removeReference();
00296 
00297   vol1=myMesh->getArea(fam2);
00298   lgth=vol1->getValueLength();
00299   vals=vol1->getValue();
00300   if(lgth==8)
00301     nbOfPtsForTest++;
00302   const double REFAreaForQuad[8]=
00303     {
00304       2.1213203435596424, 2.8284271247461903, 4.4721359549995796, 4.4721359549995796, 
00305       2.8284271247461903, 2.1213203435596428, 3.6798724963767362, 4
00306     };
00307   for(i=0;i<8;i++)
00308     if(fabs(REFAreaForQuad[i]-vals[i])<1e-12)
00309       nbOfPtsForTest++;
00310   vol1->removeReference();
00311 
00312   vol1=myMesh->getArea(fam3);
00313   lgth=vol1->getValueLength();
00314   vals=vol1->getValue();
00315   if(lgth==6)
00316     nbOfPtsForTest++;
00317   const double REFAreaForTri[6]=
00318     {
00319       2.9580398915498081, 1.4142135623730951, 2.2360679774997898, 
00320       3.3541019662496847, 3.3541019662496847, 2.2360679774997898
00321     };
00322   for(i=0;i<6;i++)
00323     if(fabs(REFAreaForTri[i]-vals[i])<1e-12)
00324       nbOfPtsForTest++;
00325   vol1->removeReference();
00326   if(nbOfPtsForTest!=38)
00327     {
00328       cout << "TEST3 K0 ! : Error in calculation of basic properties !!!" << endl;
00329       return 1;
00330     }
00331   // TEST 4 -- CHECK FOR Reverse descending using getBoundaryElements.
00332   nbOfPtsForTest=0;
00333   SUPPORT *bound=myMesh->getBoundaryElements(MED_NODE);
00334   if(bound->getNumberOfElements(MED_ALL_ELEMENTS)==19)
00335     nbOfPtsForTest++;
00336   if(bound->isOnAllElements())
00337     nbOfPtsForTest++;
00338   if(nbOfPtsForTest!=2)
00339     {
00340       cout << "TEST4 K0 ! : Error in getBoundaryElements probably due to Reverse descending !!!" << endl;
00341       return 1;
00342     }
00343   bound->removeReference();
00345   cout << "ALL TESTS OK !!!" << endl;
00346   myMesh->removeReference();
00347   return 0;
00348 }