Back to index

salome-med  6.5.0
MEDMEMTest_Extractor.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 // File      : MEDMEMTest_Extractor.cxx
00021 // Created   : Mon Dec 22 11:12:57 2008
00022 // Author    : Edward AGAPOV (eap)
00023 //
00024 #include "MEDMEMTest.hxx"
00025 
00026 #include "MEDMEM_Extractor.hxx"
00027 #include "MEDMEM_Meshing.hxx"
00028 
00029 #include <cppunit/TestAssert.h>
00030 
00031 using namespace std;
00032 using namespace MEDMEM;
00033 
00034 //================================================================================
00038 //================================================================================
00039 
00040 static void test_extractLine( Extractor* extractor,
00041                               const double* coords, const double* direction,
00042                               int nbSegments,
00043                               const char* name,
00044                               const string& result_file)
00045 {
00046   FIELD<double>* resField =0;
00047   CPPUNIT_ASSERT_NO_THROW( resField = extractor->extractLine(coords,direction));
00048   CPPUNIT_ASSERT( bool( resField ) == bool( nbSegments > 0 ));
00049 
00050   // store extracted field
00051   if ( resField )
00052     {
00053       const GMESH* mesh = resField->getSupport()->getMesh();
00054       mesh->write( MED_DRIVER, result_file, name );
00055       resField->write(MED_DRIVER, result_file );
00056 
00057       CPPUNIT_ASSERT_EQUAL( nbSegments, resField->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS));
00058       CPPUNIT_ASSERT_EQUAL( nbSegments+1, resField->getSupport()->getMesh()->getNumberOfNodes());
00059       resField->removeReference();
00060     }
00061 }
00062 
00063 void MEDMEMTest::testExtractor()
00064 {
00065   // =======================
00066   // TEST 3D->2D extraction
00067   // =======================
00068 
00069   // read a mesh from a MED file
00070   string filename  = getResourceFile("pointe.med");
00071   string meshname  = "maa1";
00072   string fieldname = "fieldcelldoublescalar";
00073 
00074   string result_file = makeTmpFile("extracted_2D.med", filename);
00075 
00076   // To remove tmp files from disk
00077   MEDMEMTest_TmpFilesRemover aRemover;
00078   aRemover.Register(result_file);
00079 
00080   // field with no mesh
00081   FIELD<double> * aField1 = new FIELD<double> (MED_DRIVER, filename, fieldname);
00082   CPPUNIT_ASSERT_THROW( new Extractor(*aField1), MEDEXCEPTION);
00083 
00084   const SUPPORT * aSupport = aField1->getSupport();
00085   MESH * aMesh = new MESH(MED_DRIVER, filename, aSupport->getMeshName());
00086   aSupport->setMesh(aMesh);
00087 
00088   Extractor* extractor = 0;
00089   CPPUNIT_ASSERT_NO_THROW( extractor = new Extractor(*aField1));
00090 
00091   FIELD<double>* resField = 0;
00092   double coords [3] = 
00093     {
00094       0,0,0
00095     };
00096   {
00097     // bad normal
00098     double normal [3] = 
00099       {
00100         0,0,0
00101       };
00102     CPPUNIT_ASSERT_THROW( extractor->extractPlane(coords,normal), MEDEXCEPTION);
00103   }
00104   double normal [3] = 
00105     {
00106       10,0,0
00107     };
00108   {
00109     // no intersection
00110     double coords [3] = 
00111       {
00112         10,0,0
00113       };
00114     CPPUNIT_ASSERT_NO_THROW( resField = extractor->extractPlane(coords,normal));
00115     CPPUNIT_ASSERT( !bool( resField ));
00116   }
00117   CPPUNIT_ASSERT_NO_THROW( resField = extractor->extractPlane(coords,normal ));
00118 
00119   // store extracted mesh
00120   const GMESH* mesh = resField->getSupport()->getMesh();
00121   mesh->write( MED_DRIVER, result_file );
00122   CPPUNIT_ASSERT_EQUAL( 2, mesh->getNumberOfTypes(MED_CELL));
00123   CPPUNIT_ASSERT_EQUAL( 6, mesh->getNumberOfElements(MED_CELL, MED_TRIA3));
00124   CPPUNIT_ASSERT_EQUAL( 2, mesh->getNumberOfElements(MED_CELL, MED_QUAD4));
00125 
00126   // store extracted field
00127   resField->write(MED_DRIVER, result_file);
00128 
00129   delete extractor; extractor=0;
00130   aMesh->removeReference(); aMesh=0;
00131   resField->removeReference(); resField=0;
00132   aField1->removeReference(); aField1=0;
00133 
00134   // ===================================================================================
00135   // TEST 2D->1D extraction in 3d space
00136   // ===================================
00137 
00138   // create 3d shell
00139   int SpaceDimension = 3 ;
00140   int NumberOfNodes = 8 ;
00141   double Coordinates[24] = 
00142     {
00143       0.0, 0.0, 200,
00144       0.0, 0.0, 0.0,
00145       0.0, 200, 200,
00146       0.0, 200, 0.0,
00147       200, 0.0, 200,
00148       200, 0.0, 0.0,
00149       200, 200, 200,
00150       200, 200, 0.0
00151     };
00152   MESHING* myMeshing  = new MESHING();
00153   myMeshing->setName("shell") ;
00154 
00155   myMeshing->setCoordinates(SpaceDimension,NumberOfNodes,Coordinates,
00156                             "CARTESIAN",MED_FULL_INTERLACE);
00157   string Names[3] = 
00158     {
00159       "X","Y","Z"
00160     } ;
00161   myMeshing->setCoordinatesNames(Names);
00162   string Units[3] = 
00163     {
00164       "cm","cm","cm"
00165     } ;
00166   myMeshing->setCoordinatesUnits(Units) ;
00167 
00168   // define conectivities
00169 
00170   const int NumberOfTypes = 1;
00171   medGeometryElement Types[NumberOfTypes] = 
00172     {
00173       MED_TRIA3
00174     };
00175   const int NumberOfElements[NumberOfTypes] = 
00176     {
00177       12
00178     };
00179 
00180   myMeshing->setNumberOfTypes(NumberOfTypes,MED_CELL);
00181   myMeshing->setTypes(Types,MED_CELL);
00182   myMeshing->setNumberOfElements(NumberOfElements,MED_CELL);
00183 
00184   const int sizeTria = 12*3 ;
00185   int ConnectivityTria[sizeTria]= 
00186     {
00187       2, 5, 1, 
00188       4, 1, 3, 
00189       3, 1, 7, 
00190       4, 8, 2, 
00191       1, 5, 7, 
00192       6, 5, 2, 
00193       2, 8, 6, 
00194       8, 7, 5, 
00195       6, 8, 5, 
00196       2, 1, 4, 
00197       8, 4, 7, 
00198       4, 3, 7 
00199     };
00200   myMeshing->setConnectivity(MED_CELL,MED_TRIA3,ConnectivityTria);
00201 
00202   // store input mesh
00203 
00204   fieldname   = "doubleOnTria";
00205   result_file = makeTmpFile("extracted_1D.med");
00206 
00207   // To remove tmp files from disk
00208   aRemover.Register(result_file);
00209 
00210   CPPUNIT_ASSERT_NO_THROW( myMeshing->write( MED_DRIVER, result_file ));
00211 
00212   // Make input field
00213 
00214   aSupport = myMeshing->getSupportOnAll( MED_CELL );
00215   FIELD<double>* inField = new FIELD<double>( aSupport, 1 );
00216   inField->setName(fieldname);
00217   string str = "";
00218   inField->setComponentsNames( &str );
00219   inField->setDescription( str );
00220   inField->setComponentsDescriptions( &str );
00221   inField->setMEDComponentsUnits( &str );
00222 
00223   vector<double> value( NumberOfElements[0] );
00224   for ( unsigned i = 0; i < value.size(); ++i )
00225     value[i] = double ( i % 10 );
00226   inField->setValue( &value[0] );
00227 
00228   // store input field
00229   int drv = inField->addDriver(MED_DRIVER, result_file, fieldname);
00230   inField->write(drv);
00231 
00232   // Extraction
00233   {
00234     double coords [3] = 
00235       {
00236         0,0,0
00237       };
00238     double normal [3] = 
00239       {
00240         -2, -2, 1 
00241       };
00242     CPPUNIT_ASSERT_NO_THROW( extractor = new Extractor(*inField));
00243     CPPUNIT_ASSERT_NO_THROW( resField = extractor->extractPlane(coords,normal));
00244     CPPUNIT_ASSERT( bool( resField ));
00245   }
00246 
00247   // store extracted mesh
00248   mesh = resField->getSupport()->getMesh();
00249   mesh->write( MED_DRIVER, result_file );
00250   //   CPPUNIT_ASSERT_EQUAL( 2, mesh->getNumberOfTypes(MED_CELL));
00251   //   CPPUNIT_ASSERT_EQUAL( 6, mesh->getNumberOfElements(MED_CELL, MED_TRIA3));
00252   //   CPPUNIT_ASSERT_EQUAL( 2, mesh->getNumberOfElements(MED_CELL, MED_QUAD4));
00253 
00254   // store extracted field
00255   drv = resField->addDriver(MED_DRIVER, result_file, resField->getName());
00256   resField->write( drv );
00257 
00258   aSupport=0;
00259   myMeshing->removeReference(); myMeshing=0;
00260   inField->removeReference(); inField=0;
00261   delete extractor; extractor=0;
00262   resField->removeReference(); resField=0;
00263 
00264   // ======================================================================================
00265   // TEST 3D->1D extraction
00266   // =======================
00267 
00268   // create 3d 2x2x2 hexa box
00269   myMeshing  = new MESHING();
00270   {
00271     SpaceDimension = 3 ;
00272     NumberOfNodes = 27 ;
00273     double Coordinates[27*3] = 
00274       {
00275         0.0, 0.0, 200,
00276         0.0, 0.0, 0.0,
00277         0.0, 200, 200,
00278         0.0, 200, 0.0,
00279         200, 0.0, 200,
00280         200, 0.0, 0.0,
00281         200, 200, 200,
00282         200, 200, 0.0,
00283         0.0, 0.0, 100,
00284         0.0, 100, 200,
00285         0.0, 200, 100,
00286         0.0, 100, 0.0,
00287         200, 0.0, 100,
00288         200, 100, 200,
00289         200, 200, 100,
00290         200, 100, 0.0,
00291         100, 0.0, 0.0,
00292         100, 0.0, 200,
00293         100, 200, 0.0,
00294         100, 200, 200,
00295         0.0, 100, 100,
00296         200, 100, 100,
00297         100, 0.0, 100,
00298         100, 200, 100,
00299         100, 100, 0.0,
00300         100, 100, 200,
00301         100, 100, 100,
00302       };
00303     myMeshing->setName("box") ;
00304 
00305     myMeshing->setCoordinates(SpaceDimension,NumberOfNodes,Coordinates,
00306                               "CARTESIAN",MED_FULL_INTERLACE);
00307     string Names[3] = 
00308       {
00309         "X","Y","Z"
00310       } ;
00311     myMeshing->setCoordinatesNames(Names);
00312     string Units[3] = 
00313       {
00314         "cm","cm","cm"
00315       } ;
00316     myMeshing->setCoordinatesUnits(Units) ;
00317 
00318     // define conectivities
00319 
00320     const int NumberOfTypes = 1;
00321     medGeometryElement Types[NumberOfTypes] = 
00322       {
00323         MED_HEXA8
00324       };
00325     const int NumberOfElements[NumberOfTypes] = 
00326       {
00327         8
00328       };
00329 
00330     myMeshing->setNumberOfTypes(NumberOfTypes,MED_CELL);
00331     myMeshing->setTypes(Types,MED_CELL);
00332     myMeshing->setNumberOfElements(NumberOfElements,MED_CELL);
00333 
00334     int ConnectivityHex[8*8]= 
00335       {
00336         9, 23, 18, 1, 21, 27, 26, 10, 
00337         25, 16, 22, 27, 19, 8, 15, 24, 
00338         27, 22, 14, 26, 24, 15, 7, 20, 
00339         17, 6, 13, 23, 25, 16, 22, 27, 
00340         21, 27, 26, 10, 11, 24, 20, 3, 
00341         2, 17, 23, 9, 12, 25, 27, 21, 
00342         23, 13, 5, 18, 27, 22, 14, 26, 
00343         12, 25, 27, 21, 4, 19, 24, 11 
00344       };
00345     myMeshing->setConnectivity(MED_CELL,MED_HEXA8,ConnectivityHex);
00346   }
00347   // store input mesh
00348 
00349   fieldname   = "doubleOnHex";
00350   result_file = makeTmpFile("extracted3D_1D.med");
00351 
00352   // To remove tmp files from disk
00353   aRemover.Register(result_file);
00354 
00355   drv = myMeshing->addDriver( MED_DRIVER, result_file, myMeshing->getName() );
00356   CPPUNIT_ASSERT_NO_THROW( myMeshing->write(drv) );
00357 
00358   // Make input field
00359 
00360   aSupport = myMeshing->getSupportOnAll( MED_CELL );
00361   inField = new FIELD<double>( aSupport, 1 );
00362   inField->setName(fieldname);
00363   inField->setComponentsNames( &str );
00364   inField->setDescription( str );
00365   inField->setComponentsDescriptions( &str );
00366   inField->setMEDComponentsUnits( &str );
00367 
00368   value.resize( NumberOfElements[0] );
00369   for ( unsigned i = 0; i < value.size(); ++i )
00370     value[i] = double (i+1);
00371   inField->setValue( &value[0] );
00372   // store input field
00373   drv = inField->addDriver(MED_DRIVER, result_file, fieldname);
00374   inField->write(drv);
00375 
00376   // Extraction
00377   CPPUNIT_ASSERT_NO_THROW( extractor = new Extractor(*inField));
00378   {
00379     // Test corner-corner intersection
00380     double coords [3] = 
00381       {
00382         1, 1, 1 
00383       };
00384     double direct [3] = 
00385       {
00386         2, 2, 2 
00387       };
00388     test_extractLine( extractor, coords, direct, 2, "corner-corner", result_file );
00389   }
00390   {
00391     // Test corner-face intersection
00392     double coords [3] = 
00393       {
00394         0,0,0
00395       };
00396     double direct [3] = 
00397       {
00398         2, 1, 1 
00399       };
00400     test_extractLine( extractor, coords, direct, 2, "corner-face", result_file );
00401   }
00402   {
00403     // Test on-face intersection
00404     double coords [3] = 
00405       {
00406         -2, 0,-1 
00407       };
00408     double direct [3] = 
00409       {
00410         2, 0, 1 
00411       };
00412     test_extractLine( extractor, coords, direct, 2, "on-face", result_file );
00413   }
00414   {
00415     // Test between-cells intersection
00416     double coords [3] = 
00417       {
00418         100,0,0
00419       };
00420     double direct [3] = 
00421       {
00422         0, 2, 2 
00423       };
00424     test_extractLine( extractor, coords, direct, 2, "between-cells", result_file );
00425   }
00426   {
00427     // Test between-cells-entrance intersection
00428     double coords [3] = 
00429       {
00430         100,0,0
00431       };
00432     double direct [3] = 
00433       {
00434         1, 2, 2 
00435       };
00436     test_extractLine( extractor, coords, direct, 2, "between-cells-entrance", result_file );
00437   }
00438   {
00439     // Test edge-entrance intersection
00440     double coords [3] = 
00441       {
00442         100,0,50
00443       };
00444     double direct [3] = 
00445       {
00446         1, 2, 2 
00447       };
00448     test_extractLine( extractor, coords, direct, 3, "edge-entrance", result_file );
00449   }
00450   {
00451     // Test touch intersection - expect no result
00452     double coords [3] = 
00453       {
00454         0,0,0
00455       };
00456     double direct [3] = 
00457       {
00458         0, 2, -2 
00459       };
00460     test_extractLine( extractor, coords, direct, 0, "touch", result_file );
00461   }
00462   {
00463     // Test face-face intersection
00464     double coords [3] = 
00465       {
00466         50,50,0
00467       };
00468     double direct [3] = 
00469       {
00470         2, 2, 0 
00471       };
00472     test_extractLine( extractor, coords, direct, 2, "corner-corner-on", result_file );
00473   }
00474   {
00475     // Test face-face intersection
00476     double coords [3] = 
00477       {
00478         50,50,0
00479       };
00480     double direct [3] = 
00481       {
00482         2, 2, 2 
00483       };
00484     test_extractLine( extractor, coords, direct, 3, "face-face", result_file );
00485   }
00486   {
00487     // Test external edge intersection
00488     double coords [3] = 
00489       {
00490         0, 0,200 
00491       };
00492     double direct [3] = 
00493       {
00494         -1, 0, 0 
00495       };
00496     test_extractLine( extractor, coords, direct, 2, "external edge", result_file );
00497   }
00498   {
00499     // Test internal edge intersection
00500     double coords [3] = 
00501       {
00502         100,0,100
00503       };
00504     double direct [3] = 
00505       {
00506         0, -2, 0 
00507       };
00508     test_extractLine( extractor, coords, direct, 2, "internal edge", result_file );
00509   }
00510 
00511   inField->setSupport(0);
00512   aSupport=0;
00513   delete extractor; extractor=0;
00514   myMeshing->removeReference(); myMeshing=0;
00515 
00516   // ======================================================================================
00517   // TEST 3D->1D extraction on a large model
00518   // =======================================
00519 
00520   // read a mesh
00521   filename = getResourceFile("geomMesh22.med");
00522   meshname = "GeomMesh";
00523   aMesh = new MESH(MED_DRIVER, filename, meshname);
00524   aSupport = aMesh->getSupportOnAll( MED_CELL );
00525 
00526   // make a field
00527   int nbValues = aSupport->getNumberOfElements(MED_ALL_ELEMENTS);
00528   inField->setSupport( aSupport );
00529   inField->allocValue( 1, nbValues );
00530   {
00531     double* value = const_cast<double*>( inField->getValue() );
00532     for ( int i = 0; i < nbValues; ++i )
00533       value[i] = double ( i % 7 );
00534   }
00535   // extract a field
00536   CPPUNIT_ASSERT_NO_THROW( extractor = new Extractor(*inField));
00537   {
00538     double coords [3] = 
00539       {
00540         20,0,10
00541       };
00542     double direct [3] = 
00543       {
00544         1, 1,1.5 
00545       };
00546     CPPUNIT_ASSERT_NO_THROW( resField = extractor->extractLine(coords,direct));
00547   }
00548   CPPUNIT_ASSERT( resField );
00549   CPPUNIT_ASSERT_EQUAL( 31, resField->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS));
00550   CPPUNIT_ASSERT_EQUAL( 33, resField->getSupport()->getMesh()->getNumberOfNodes());
00551   delete extractor;
00552   resField->removeReference(); resField=0;
00553   aMesh->removeReference(); aMesh=0;
00554   aSupport=0;
00555   inField->removeReference(); inField=0;
00556 }
00557