Back to index

salome-med  6.5.0
MEDMEMTest_Field.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 #include "MEDMEMTest.hxx"
00021 
00022 #include "MEDMEM_FieldConvert.hxx"
00023 #include "MEDMEM_Field.hxx"
00024 #include "MEDMEM_MedMeshDriver.hxx"
00025 #include "MEDMEM_Grid.hxx"
00026 #include "MEDMEM_Group.hxx"
00027 #include "MEDMEM_Mesh.hxx"
00028 #include "MEDMEM_Meshing.hxx"
00029 #include "MEDMEM_Support.hxx"
00030 #include "MEDMEM_VtkMeshDriver.hxx"
00031 
00032 
00033 #include <cppunit/TestAssert.h>
00034 #include <sstream>
00035 #include <cmath>
00036 
00037 // use this define to enable lines, execution of which leads to Segmentation Fault
00038 //#define ENABLE_FAULTS
00039 
00040 // use this define to enable CPPUNIT asserts and fails, showing bugs
00041 //#define ENABLE_FORCED_FAILURES
00042 
00043 using namespace std;
00044 using namespace MEDMEM;
00045 
00046 // #14,15: MEDMEMTest_Field.cxx
00047 // Check methods from MEDMEM_Field.hxx, MEDMEM_FieldConvert.hxx
00048 
00238  static void compareField_(const FIELD_ * theField_1, const FIELD_ * theField_2, bool isFIELD, bool isValue)
00239 {
00240   // name, description, support
00241   CPPUNIT_ASSERT_EQUAL(theField_1->getName(), theField_2->getName());
00242   CPPUNIT_ASSERT_EQUAL(theField_1->getDescription(), theField_2->getDescription());
00243   CPPUNIT_ASSERT_EQUAL(theField_1->getSupport(), theField_2->getSupport());
00244 
00245   // components information
00246   int aNbComps = theField_1->getNumberOfComponents();
00247   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_2->getNumberOfComponents());
00248 
00249   for (int i = 1; i <= aNbComps; i++)
00250     {
00251       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentName(i), theField_2->getComponentName(i));
00252       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentDescription(i), theField_2->getComponentDescription(i));
00253       CPPUNIT_ASSERT_EQUAL(theField_1->getMEDComponentUnit(i), theField_2->getMEDComponentUnit(i));
00254     }
00255 
00256   // iteration information
00257   CPPUNIT_ASSERT_EQUAL(theField_1->getIterationNumber(), theField_2->getIterationNumber());
00258   CPPUNIT_ASSERT_EQUAL(theField_1->getOrderNumber(), theField_2->getOrderNumber());
00259   CPPUNIT_ASSERT_DOUBLES_EQUAL(theField_1->getTime(), theField_2->getTime(), 0.0000001);
00260 
00261   // Value
00262   int nbOfValues = theField_1->getNumberOfValues();
00263   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_2->getNumberOfValues());
00264 
00265   if (isFIELD)
00266     {
00267       // Value type and Interlacing type
00268       CPPUNIT_ASSERT_EQUAL(theField_1->getValueType(), theField_2->getValueType());
00269       CPPUNIT_ASSERT_EQUAL(theField_1->getInterlacingType(), theField_2->getInterlacingType());
00270 
00271       // Gauss Presence
00272       if (isValue)
00273         {
00274           CPPUNIT_ASSERT_EQUAL(theField_1->getGaussPresence(), theField_2->getGaussPresence());
00275         }
00276       else
00277         {
00278           CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
00279           CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
00280         }
00281     }
00282   else
00283     {
00284       CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
00285       CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
00286     }
00287 }
00288 
00289 static void checkField_(FIELD_ * theField_, const SUPPORT * theSupport,
00290                         MED_EN::med_type_champ theValueType,
00291                         MED_EN::medModeSwitch theInterlace)
00292 {
00293   // name
00294   const string aFieldName = "a_name_of_a_field";
00295   theField_->setName(aFieldName);
00296   CPPUNIT_ASSERT_EQUAL(aFieldName, theField_->getName());
00297 
00298   // description
00299   const string aFieldDescr = "a_description_of_a_field";
00300   theField_->setDescription(aFieldDescr);
00301   CPPUNIT_ASSERT_EQUAL(aFieldDescr, theField_->getDescription());
00302 
00303   // support
00304   theField_->setSupport(theSupport);
00305   CPPUNIT_ASSERT(theField_->getSupport() == theSupport);
00306 
00307   // components information
00308   int aNbComps = 3;
00309 
00310   string aCompsNames[3] =
00311     {
00312       "Vx", "Vy", "Vz"
00313     };
00314   string aCompsDescs[3] =
00315     {
00316       "vitesse selon x", "vitesse selon y", "vitesse selon z"
00317     };
00318   string aCompsUnits[3] =
00319     {
00320       "m.s-1", "m.s-1", "m.s-1"
00321     };
00322 
00323   theField_->setNumberOfComponents(aNbComps);
00324   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_->getNumberOfComponents());
00325 
00326   theField_->setComponentsNames(aCompsNames);
00327 
00328   //#ifdef ENABLE_FAULTS
00329   try
00330     {
00331       theField_->setNumberOfComponents(7);
00332       // Segmentation fault here because array of components names is not resized
00333       for (int i = 1; i <= 7; i++)
00334         {
00335           theField_->setComponentName(i, "AnyComponent");
00336         }
00337     }
00338   catch (MEDEXCEPTION& ex)
00339     {
00340       // Ok, it is good to have MEDEXCEPTION here
00341     }
00342   catch (...)
00343     {
00344       CPPUNIT_FAIL("Unknown exception cought");
00345     }
00346   // restore components names
00347   theField_->setNumberOfComponents(aNbComps);
00348   theField_->setComponentsNames(aCompsNames);
00349 
00350   theField_->setComponentsDescriptions(aCompsDescs);
00351   theField_->setMEDComponentsUnits(aCompsUnits);
00352 
00353   const string * aCompsNamesBack = theField_->getComponentsNames();
00354   const string * aCompsDescsBack = theField_->getComponentsDescriptions();
00355   const string * aCompsUnitsBack = theField_->getMEDComponentsUnits();
00356   for (int i = 1; i <= aNbComps; i++)
00357     {
00358       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], theField_->getComponentName(i));
00359       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], aCompsNames[i-1]);
00360 
00361       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], theField_->getComponentDescription(i));
00362       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], aCompsDescs[i-1]);
00363 
00364       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], theField_->getMEDComponentUnit(i));
00365       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], aCompsUnits[i-1]);
00366     }
00367 
00368   const string aCompName2 ("Name of second component");
00369   const string aCompDesc2 ("Description of second component");
00370   const string aCompUnit2 ("Unit of second MED component");
00371 
00372   theField_->setComponentName(2, aCompName2);
00373   theField_->setComponentDescription(2, aCompDesc2);
00374   theField_->setMEDComponentUnit(2, aCompUnit2);
00375 
00376   const string * aCompsNamesBack2 = theField_->getComponentsNames();
00377   const string * aCompsDescsBack2 = theField_->getComponentsDescriptions();
00378   const string * aCompsUnitsBack2 = theField_->getMEDComponentsUnits();
00379 
00380   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], theField_->getComponentName(2));
00381   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], aCompName2);
00382 
00383   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], theField_->getComponentDescription(2));
00384   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], aCompDesc2);
00385 
00386   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], theField_->getMEDComponentUnit(2));
00387   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], aCompUnit2);
00388 
00389   CPPUNIT_ASSERT_THROW(theField_->setComponentName(0, "str"), MEDEXCEPTION);
00390   CPPUNIT_ASSERT_THROW(theField_->setComponentName(aNbComps + 1, "str"), MEDEXCEPTION);
00391   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(0, "str"), MEDEXCEPTION);
00392   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(aNbComps + 1, "str"), MEDEXCEPTION);
00393   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(0, "str"), MEDEXCEPTION);
00394   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(aNbComps + 1, "str"), MEDEXCEPTION);
00395 
00396   // iteration information
00397   int anIterNumber = 10; // set value to MED_NOPDT if undefined (default)
00398   theField_->setIterationNumber(anIterNumber);
00399   CPPUNIT_ASSERT_EQUAL(anIterNumber, theField_->getIterationNumber());
00400 
00401   int anOrderNumber = 1; // set value to MED_NONOR if undefined (default)
00402   theField_->setOrderNumber(anOrderNumber);
00403   CPPUNIT_ASSERT_EQUAL(anOrderNumber, theField_->getOrderNumber());
00404 
00405   double aTime = 3.435678; // in second
00406   theField_->setTime(aTime);
00407   CPPUNIT_ASSERT_DOUBLES_EQUAL(aTime, theField_->getTime(), 0.0000001);
00408 
00409   // Value
00410   int nbOfValues = 10;
00411   // dangerous method, because it does not reallocate values array
00412   theField_->setNumberOfValues(nbOfValues);
00413   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_->getNumberOfValues());
00414 
00415   // Value type and Interlacing type
00416   CPPUNIT_ASSERT_EQUAL(theValueType, theField_->getValueType());
00417   CPPUNIT_ASSERT_EQUAL(theInterlace, theField_->getInterlacingType());
00418 }
00419 
00420 template<class T, class INTERLACING_TAG>
00421 void compareField(const FIELD<T, INTERLACING_TAG> * theField_1,
00422                   const FIELD<T, INTERLACING_TAG> * theField_2, bool isValue)
00423 {
00424   // compare FIELD_ part
00425   compareField_(theField_1, theField_2, /*isFIELD = */true, isValue);
00426 
00427   // compare FIELD part
00428   // TO DO
00429 }
00430 
00431 template<class T, class INTERLACING_TAG>
00432 void checkField (FIELD<T, INTERLACING_TAG> * theField, const SUPPORT * theSupport)
00433 {
00434   // check FIELD_ part
00435   MED_EN::med_type_champ aValueType = SET_VALUE_TYPE<T>::_valueType;
00436   MED_EN::medModeSwitch  anInterlace = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
00437   checkField_(theField, theSupport, aValueType, anInterlace);
00438 
00439   // check FIELD part
00440 
00441   // filling by support charackteristics (NOT IMPLEMENTED METHODS!!!):
00442   // value type must be MED_REEL64 (i.e. a FIELD<double>) for these methods,
00443   // nb. of components must be equal 1 (for Volume, Area, Length) or
00444   // space dimension (for Normal, Barycenter, )
00445   {
00446     const GMESH* aMesh = theSupport->getMesh();
00447     int spaceDim = 3;
00448     if (aMesh) spaceDim = aMesh->getSpaceDimension();
00449     theField->deallocValue();
00450     theField->allocValue(/*NumberOfComponents = */spaceDim + 1);
00451 
00452     theField->deallocValue();
00453     theField->allocValue(/*NumberOfComponents = */1);
00454 
00455     if (aMesh) 
00456       {
00457         theField->deallocValue();
00458         theField->allocValue(/*NumberOfComponents = */spaceDim);
00459       }
00460 
00461     if (aMesh)
00462       {
00463         theField->deallocValue();
00464         theField->allocValue(/*NumberOfComponents = */spaceDim);
00465         //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
00466         // getVolume() etc. does nothing
00467         //       if (aValueType == MED_EN::MED_REEL64)
00468         {
00469           //         CPPUNIT_ASSERT_NO_THROW(theField->getNormal());
00470           //         CPPUNIT_ASSERT_NO_THROW(theField->getBarycenter());
00471           //
00472         }
00473         //       else
00474         {
00475           //         CPPUNIT_ASSERT_THROW(theField->getNormal(), MEDEXCEPTION);
00476           //         CPPUNIT_ASSERT_THROW(theField->getBarycenter(), MEDEXCEPTION);
00477           //
00478         }
00479       }
00480   }
00481 
00482   // values
00483   theField->deallocValue();
00484   theField->allocValue(/*NumberOfComponents = */2);
00485   int nbElemSupport = theSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
00486   CPPUNIT_ASSERT_EQUAL(nbElemSupport, theField->getNumberOfValues());
00487 
00488   theField->deallocValue();
00489   CPPUNIT_ASSERT_THROW(theField->getGaussPresence(), MEDEXCEPTION);
00490 
00491   // copy constructor
00492   FIELD<T, INTERLACING_TAG> *aField_copy1= new FIELD<T, INTERLACING_TAG>(*theField);
00493   compareField(theField, aField_copy1, /*isValue = */false);
00494   aField_copy1->removeReference();
00495 
00496   // operator=
00497   FIELD<T, INTERLACING_TAG> *aField_copy2=new FIELD<T, INTERLACING_TAG>();
00498   *aField_copy2 = *theField;
00499   compareField(theField, aField_copy2, /*isValue = */false);
00500   aField_copy2->removeReference();
00501 }
00502 
00503 template<class T>
00504 FIELD<T> * createFieldOnGroup(MESH* theMesh, const GROUP* theGroup,
00505                               const string theName, const string theDescr)
00506 {
00507   FIELD<T> * aFieldOnGroup = new FIELD<T>(theGroup, /*NumberOfComponents = */2);
00508 
00509   aFieldOnGroup->setName(theName);
00510   aFieldOnGroup->setDescription(theDescr);
00511 
00512   string aCompsNames[2] =
00513     {
00514       "Pos", "Neg"
00515     };
00516   string aCompsDescs[2] =
00517     {
00518       "+", "-"
00519     };
00520   string aCompsUnits[2] =
00521     {
00522       "unit1", "unit2"
00523     };
00524 
00525   aFieldOnGroup->setComponentsNames(aCompsNames);
00526   aFieldOnGroup->setComponentsDescriptions(aCompsDescs);
00527   aFieldOnGroup->setMEDComponentsUnits(aCompsUnits);
00528 
00529   return aFieldOnGroup;
00530 }
00531 
00532 double plus13 (double val);
00533 double plus13 (double val)
00534 {
00535   return val + 13;
00536 }
00537 
00538 // function to calculate field values from coordinates of an element
00539 // typedef void (*myFuncType)(const double * temp, T* output);
00540 // size of temp array = space dim = 3
00541 // size of output array = nb. comps = 2
00542 static void proj2d (const double * temp, double* output)
00543 {
00544   // dimetric projection with coefficients:
00545   // 1.0 along Oy and Oz, 0.5 along Ox
00546   //
00547   //    ^ z (y_)
00548   //    |
00549   //    |
00550   //    .----> y (x_)
00551   //   /
00552   //  L x
00553   //
00554   // x_ = y - x * sqrt(2.) / 4.
00555   // y_ = z - x * sqrt(2.) / 4.
00556 
00557   double dx = temp[0] * std::sqrt(2.) / 4.;
00558   output[0] = temp[1] - dx;
00559   output[1] = temp[2] - dx;
00560 }
00561 
00562 static void testDrivers()
00563 {
00564   string filename_rd                  = getResourceFile("pointe.med");
00565   string filename_wr                  = makeTmpFile("myMedFieldfile.med", filename_rd);
00566   string filename_support_wr          = makeTmpFile("myMedSupportFiledfile.med");
00567   string filename22_rd                = getResourceFile("pointe.med");
00568   string filenamevtk_wr               = makeTmpFile("myMedFieldfile22.vtk");
00569 
00570   string fieldname_celldouble_rd      = "fieldcelldoublescalar";
00571   string fieldname_celldouble_wr      = fieldname_celldouble_rd + "_cpy";
00572   string fieldname_nodeint_rd         = "fieldnodeint";
00573   string fieldname_nodeint_wr         = fieldname_nodeint_rd + "_cpy";
00574   string fieldname_nodeint_wr1        = fieldname_nodeint_rd + "_cpy1";
00575   string meshname                     = "maa1";
00576 
00577   // To remove tmp files from disk
00578   MEDMEMTest_TmpFilesRemover aRemover;
00579   aRemover.Register(filename_wr);
00580   aRemover.Register(filenamevtk_wr);
00581   aRemover.Register(filename_support_wr);
00582 
00583   FIELD<int> *aInvalidField=new FIELD<int>();
00584   //must throw becase only VTK_DRIVER or MED_DRIVER may be specified as driverType for FIELD
00585   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(NO_DRIVER, filename_rd, fieldname_nodeint_rd)),
00586                        MEDEXCEPTION);
00587   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(GIBI_DRIVER, filename_rd, fieldname_nodeint_rd)),
00588                        MEDEXCEPTION);
00589   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(PORFLOW_DRIVER, filename_rd, fieldname_nodeint_rd)),
00590                        MEDEXCEPTION);
00591   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(ASCII_DRIVER, filename_rd, fieldname_nodeint_rd)),
00592                        MEDEXCEPTION);
00593   aInvalidField->removeReference();
00595   //TestRead Part//
00597   FIELD<double> *aField_1 = NULL;
00598   CPPUNIT_ASSERT_NO_THROW(aField_1 = new FIELD<double>(MED_DRIVER, filename_rd,
00599                                                        fieldname_celldouble_rd));
00600 
00601   //Test read(int index) method
00602   int IdDriver_rd = aField_1->addDriver(MED_DRIVER,filename_rd,fieldname_celldouble_rd,MED_EN::RDONLY);
00603   // TODO: throw if read for the second time
00604   // (BUG) Cannot open file, but file exist
00605   // EAP: no more pb with opening the file for the second time with "weaker" mode,
00606   // but why to re-read the field?
00607   //CPPUNIT_ASSERT_THROW(aField_1->read(IdDriver_rd),MEDEXCEPTION);
00608   CPPUNIT_ASSERT_NO_THROW(aField_1->read(IdDriver_rd));
00609 
00610   //Test read(GENDRIVER & genDriver) method
00611   FIELD<int> *aField_2 = new FIELD<int>();
00612   aField_2->setName(fieldname_nodeint_rd);
00613   {
00614     MED_FIELD_RDONLY_DRIVER<int> aMedRdFieldDriver;
00615     aMedRdFieldDriver.setFileName(filename_rd);
00616     aField_2->read( aMedRdFieldDriver );
00617   }
00618   //Test read(driverTypes driverType, const std::string & fileName);
00619   FIELD<double> * aField_3 = new FIELD<double>();
00620   aField_3->setName(fieldname_celldouble_rd);
00621   aField_3->read( MED_DRIVER, filename_rd);
00622 
00624   //Test Write Part//
00626   //int IdDriver;
00627   MESH *aMesh = new MESH(MED_DRIVER,filename_rd,meshname);
00628   const SUPPORT *aSupport = aMesh->getSupportOnAll(MED_CELL);
00629   FIELD<int> *aFieldSupport;
00630   CPPUNIT_ASSERT_THROW(aFieldSupport = 
00631                        new FIELD<int>(aSupport, MED_DRIVER,filename_rd,
00632                                       fieldname_nodeint_rd), MEDMEM::MEDEXCEPTION);
00633   aSupport = aMesh->getSupportOnAll(MED_NODE);
00634   CPPUNIT_ASSERT_NO_THROW(aFieldSupport = 
00635                           new FIELD<int>(aSupport, MED_DRIVER, filename_rd,
00636                                          fieldname_nodeint_rd));
00637   aFieldSupport->removeReference();
00638 
00639   //Test write(int index) method
00640   // Add drivers to FIELDs
00641   int IdDriver1 = aField_3->addDriver(MED_DRIVER,filename_wr,fieldname_celldouble_wr);
00642 
00643   //Trying call write(int index) method with incorrect index
00644   CPPUNIT_ASSERT_THROW(aField_3->write(IdDriver1+1),MEDEXCEPTION);
00645 
00646   //Write field to file
00647   aField_3->write(IdDriver1);
00648 
00649   CPPUNIT_ASSERT_NO_THROW(aField_3->rmDriver(IdDriver1));
00650 
00651   //Test write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
00652   {
00653     MED_FIELD_WRONLY_DRIVER<int> aMedWrFieldDriver;
00654     aMedWrFieldDriver.setFileName(filename_wr);
00655 
00656     aField_3->setName(fieldname_nodeint_wr);
00657     aField_3->write(aMedWrFieldDriver);
00658     FIELD<double> aField_3_RD;
00659     aField_3_RD.setName(fieldname_nodeint_wr);
00660     aField_3_RD.read(MED_DRIVER,filename_wr);
00661   }
00662 
00663   // Test write(driverTypes driverType, const std::string& filename
00664   //            MED_EN::med_mode_acces medMode=MED_EN::RDWR)
00665   aField_3->setName(fieldname_nodeint_wr1);
00666   // wrong mode
00667   CPPUNIT_ASSERT_THROW(aField_3->write(MED_DRIVER,filename_wr,MED_EN::RDONLY), MEDEXCEPTION);
00668   aField_3->write(MED_DRIVER,filename_wr);
00669   {
00670     FIELD<double> aField_3_RD;
00671     aField_3_RD.setName(fieldname_nodeint_wr1);
00672     aField_3_RD.read(MED_DRIVER,filename_wr);
00673   }
00674   aField_3->setName("fieldname_nodeint_wr1");
00675   aField_3->write(MED_DRIVER,filename_wr, MED_EN::WRONLY);// overwrite filename_wr
00676   {
00677     FIELD<double> aField_3_RD;
00678     aField_3_RD.setName(fieldname_nodeint_wr1);
00679     CPPUNIT_ASSERT_THROW(aField_3_RD.read(MED_DRIVER,filename_wr),MEDEXCEPTION);
00680   }
00681   //Test writeAppend(int index) method
00682   //Create a vtk file
00683   MESH * aMesh_1 = new MESH(MED_DRIVER,filename22_rd, meshname);
00684   aMesh_1->write(VTK_DRIVER, filenamevtk_wr);
00685   //Create a field
00686   FIELD<int> * aField_4 =
00687     new FIELD<int>(MED_DRIVER,filename22_rd,fieldname_nodeint_rd,-1,-1,aMesh_1);
00688   //Add Driver to a field
00689   int IdDriver2 = aField_4->addDriver(VTK_DRIVER, filenamevtk_wr ,fieldname_nodeint_wr);
00690   //Trying call writeAppend() method with incorrect index
00691   CPPUNIT_ASSERT_THROW(aField_4->writeAppend(IdDriver2+1,fieldname_nodeint_wr),MEDEXCEPTION);
00692 
00693   CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(IdDriver2, fieldname_nodeint_wr));
00694 
00695   //Test writeAppend(const GENDRIVER &) method
00696   aField_4->setName(fieldname_nodeint_wr1);
00697   {
00698     VTK_FIELD_DRIVER<int> aVtkFieldDriver(filenamevtk_wr, aField_4);
00699     CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(aVtkFieldDriver));
00700   }
00701 
00702   //Delete objects
00703   aField_1->removeReference();
00704   aField_2->removeReference();
00705   aField_3->removeReference();
00706   aField_4->removeReference();
00707   aMesh->removeReference();
00708   aMesh_1->removeReference();
00709 }
00710 
00711 void MEDMEMTest::testField()
00712 {
00713   SUPPORT *anEmptySupport=new SUPPORT;
00715   // TEST 1: FIELD_ //
00717   FIELD_ *aField_=new FIELD_;
00718 
00719   // check set/get methods
00720   MED_EN::med_type_champ aValueType = MED_EN::MED_UNDEFINED_TYPE;
00721   MED_EN::medModeSwitch  anInterlace = MED_EN::MED_UNDEFINED_INTERLACE;
00722   checkField_(aField_, anEmptySupport, aValueType, anInterlace);
00723 
00724   // copy constructor
00725   // This fails (Segmentation fault) if not set:
00726   // _componentsNames or _componentsDescriptions, or _MEDComponentsUnits
00727   FIELD_ *aField_copy1=new FIELD_(*aField_);
00728   compareField_(aField_, aField_copy1, /*isFIELD = */false, /*isValue = */false);
00729   aField_copy1->removeReference();
00730   // operator=
00731   FIELD_ *aField_copy2=new FIELD_;
00732   *aField_copy2 = *aField_;
00733   compareField_(aField_, aField_copy2,/*isFIELD = */false, /*isValue = */false);
00734   aField_copy2->removeReference();
00735   aField_->removeReference();
00736 
00738   // TEST 2: FIELD<int> //
00740   FIELD<int> *aFieldInt=new FIELD<int>();
00741   checkField(aFieldInt, anEmptySupport);
00742   aFieldInt->removeReference();
00743   anEmptySupport->removeReference();
00745   // TEST 3: FIELD<double, NoInterlace> //
00747   MESH * aMesh = MEDMEMTest_createTestMesh();
00748   const GROUP* aGroup = aMesh->getGroup(MED_EN::MED_FACE, 1);
00749 
00750   FIELD<double, NoInterlace> *aFieldDouble=new FIELD<double, NoInterlace>();
00751   checkField(aFieldDouble, aGroup);
00752   aFieldDouble->removeReference();
00754   // TEST 4: FIELD<double, FullInterlace> //
00756   FIELD<double> * aFieldOnGroup1 = createFieldOnGroup<double>(aMesh, aGroup, "Linear", "N");
00757   FIELD<double> * aFieldOnGroup2 = createFieldOnGroup<double>(aMesh, aGroup, "Quadratic", "N**2");
00758 
00759   int nbVals = aFieldOnGroup1->getNumberOfValues();
00760   CPPUNIT_ASSERT(nbVals);
00761 
00762   // numbers of elements in group,
00763   // they are needed in method FIELD::setValueIJ()
00764   const int *anElems = aGroup->getnumber()->getValue();
00765   double eucl1 = 0., eucl2 = 0.;
00766 
00767   for (int i = 1; i <= nbVals; i++)
00768     {
00769       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, (double)i);
00770       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, (double)(-i));
00771 
00772       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, (double)i*i);
00773       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, (double)(-i*i));
00774 
00775       eucl1 += 2. * i * i;
00776       eucl2 += 2. * i * i * i * i;
00777     }
00778 
00779   // out of bound (inexisting 33-th component of last element)
00780   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->setValueIJ(anElems[nbVals-1], 33, 10.), MEDEXCEPTION);
00781 
00782   // normMax
00783   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals, aFieldOnGroup1->normMax(), 0.000001);
00784   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals*nbVals, aFieldOnGroup2->normMax(), 0.000001);
00785 
00786   // norm2
00787   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl1), aFieldOnGroup1->norm2(), 0.000001); // 10.4881
00788   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl2), aFieldOnGroup2->norm2(), 0.000001); // 44.2493
00789 
00790   // check getXXX methods
00791   CPPUNIT_ASSERT(!aFieldOnGroup1->getGaussPresence());
00792   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getArrayGauss(), MEDEXCEPTION);
00793   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * anArrayNoGauss =
00794     aFieldOnGroup1->getArrayNoGauss();
00795 
00796   MEDMEM_Array_ * aMEDMEM_Array_ = aFieldOnGroup1->getArray();
00797   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * aMEDMEM_Array_conv =
00798     static_cast<MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array *>(aMEDMEM_Array_);
00799 
00800   const double * aValues = aFieldOnGroup1->getValue();
00801 
00802   // out of range
00803   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(3), MEDEXCEPTION);
00804   // cannot get column in FullInterlace
00805   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(1), MEDEXCEPTION);
00806 
00807   for (int i = 1; i <= nbVals; i++)
00808     {
00809       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aFieldOnGroup1->getValueIJK(anElems[i-1], 1, 1), 0.000001);
00810       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aFieldOnGroup1->getValueIJK(anElems[i-1], 2, 1), 0.000001);
00811 
00812       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aValues[(i-1)*2 + 0], 0.000001);
00813       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aValues[(i-1)*2 + 1], 0.000001);
00814 
00815       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , anArrayNoGauss->getIJ(i, 1), 0.000001);
00816       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), anArrayNoGauss->getIJ(i, 2), 0.000001);
00817 
00818       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aMEDMEM_Array_conv->getIJ(i, 1), 0.000001);
00819       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aMEDMEM_Array_conv->getIJ(i, 2), 0.000001);
00820 
00821       const double* row_i = aFieldOnGroup1->getRow(anElems[i-1]);
00822       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , row_i[0], 0.000001);
00823       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), row_i[1], 0.000001);
00824 
00825       double vals_i [2];
00826       aFieldOnGroup1->getValueOnElement(anElems[i-1], vals_i);
00827       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , vals_i[0], 0.000001);
00828       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), vals_i[1], 0.000001);
00829     }
00830   {
00831     // Test getValueOnPoints()
00832     // -----------------------
00833     const string file = getResourceFile("cube_hexa8.med");
00834     MESH mesh(MED_DRIVER, file, "CUBE_EN_HEXA8");
00835     double value[6];
00836     {
00837       FIELD<double> fieldcelldouble(MED_DRIVER, file, "fieldcelldouble", -1,-1, &mesh);
00838       fieldcelldouble.getSupport()->setMesh( &mesh );
00839       double point[6] =
00840         {
00841           0.25, 0.75, 0.25, // the 3-d cell
00842           1.0,  1.0,  1.0
00843         };// the 8-th cell
00844       fieldcelldouble.getValueOnPoint(point, value);
00845       const double tol = 1e-20;
00846       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[0], tol );
00847       CPPUNIT_ASSERT_DOUBLES_EQUAL( 3.0, value[1], tol );
00848       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, value[2], tol );
00849       fieldcelldouble.getValueOnPoints(2, point, value);
00850       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[0], tol );
00851       CPPUNIT_ASSERT_DOUBLES_EQUAL( 3.0, value[1], tol );
00852       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, value[2], tol );
00853       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[3], tol );
00854       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[4], tol );
00855       CPPUNIT_ASSERT_DOUBLES_EQUAL( 3.0, value[5], tol );
00856     }
00857     {
00858       FIELD<double> fieldnodedouble(MED_DRIVER, file, "fieldnodedouble", -1,-1, &mesh); // t=0.0
00859       fieldnodedouble.getSupport()->setMesh( &mesh );
00860       double point[6] =
00861         {
00862           1.0, 1.0, 0.0, // the 9-th node
00863           1.0, 0.0, 1.0
00864         };// the 21-st node
00865       fieldnodedouble.getValueOnPoint(point, value);
00866       const double tol = 1e-20;
00867       CPPUNIT_ASSERT_DOUBLES_EQUAL( 5.0, value[0], tol );
00868       fieldnodedouble.getValueOnPoints(2, point, value);
00869       CPPUNIT_ASSERT_DOUBLES_EQUAL( 5.0, value[0], tol );
00870       CPPUNIT_ASSERT_DOUBLES_EQUAL( 4.0, value[1], tol );
00871 
00872       point[0] = 1.1; // point out of mesh
00873       CPPUNIT_ASSERT_THROW(fieldnodedouble.getValueOnPoint(point, value), MEDEXCEPTION);
00874     }
00875   }
00876 
00877   // modify all values of aFieldOnGroup2 by formula a*x + b (a = 2, b = 3)
00878   aFieldOnGroup2->applyLin(2., 3.);
00879   for (int i = 1; i <= nbVals; i++)
00880     {
00881       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. + 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
00882       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. - 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
00883     }
00884 
00885   // apply function plus13() to aFieldOnGroup1
00886   aFieldOnGroup1->applyFunc<plus13>();
00887   for (int i = 1; i <= nbVals; i++)
00888     {
00889       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
00890       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
00891     }
00892 
00893   // scalarProduct
00894   FIELD<double, FullInterlace> * aScalarProduct =
00895     FIELD<double, FullInterlace>::scalarProduct(*aFieldOnGroup1, *aFieldOnGroup2, /*deepCheck = */true);
00896   CPPUNIT_ASSERT_EQUAL(nbVals, aScalarProduct->getNumberOfValues());
00897   CPPUNIT_ASSERT_EQUAL(1, aScalarProduct->getNumberOfComponents());
00898   for (int i = 1; i <= nbVals; i++)
00899     {
00900       CPPUNIT_ASSERT_DOUBLES_EQUAL(78. + 4.*i*i*i, //(3. + 2.*i*i)*(13 + i) + (3. - 2.*i*i)*(13 - i)
00901                                    aScalarProduct->getValueIJ(anElems[i-1], 1), 0.000001);
00902     }
00903 
00904   // fillFromAnalytic
00905   aFieldOnGroup2->fillFromAnalytic(proj2d);
00906 
00907   double bary [3];
00908   double outp [2];
00909   const SUPPORT * aSupp = aFieldOnGroup2->getSupport();
00910   FIELD<double,FullInterlace> * barycenter = aMesh->getBarycenter(aSupp);
00911   for (int i = 1; i <= nbVals; i++)
00912     {
00913       bary[0] = barycenter->getValueIJ(anElems[i-1], 1);
00914       bary[1] = barycenter->getValueIJ(anElems[i-1], 2);
00915       bary[2] = barycenter->getValueIJ(anElems[i-1], 3);
00916 
00917       proj2d(bary, outp);
00918 
00919       //cout << "barycenter (" << bary[0] << ", " << bary[1] << ", " << bary[2] << ")" << endl;
00920       //cout << "proj2d     (" << outp[0] << ", " << outp[1] << ")" << endl;
00921 
00922       //bary (-0.666667,  0.666667, 0.666667) -> outp ( 0.902369, 0.902369)
00923       //bary ( 0.666667, -0.666667, 0.666667) -> outp (-0.902369, 0.430964)
00924       //bary ( 0.      ,  0.      , 2.      ) -> outp ( 0.      , 2.      )
00925       //bary ( 0.      ,  0.      , 3.      ) -> outp ( 0.      , 3.      )
00926       //bary (-1.      ,  0.      , 2.5     ) -> outp ( 0.353553, 2.85355 )
00927 
00928       //#ifdef ENABLE_FORCED_FAILURES
00929       // (BUG) in FIELD::fillFromAnalytic() in case of support, different from nodes:
00930       //       barycenterField in FullInterlace, but values extracted like from NoInterlace
00931       // TODO: fix MEDMEM_Field:3483 - code should depend on interlace
00932       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[0], aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
00933       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[1], aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
00934       //#endif
00935 
00936       // currently it gives values, that are wrong:
00937       //aFieldOnGroup2 row1 ( 0.902369,  0.235702)
00938       //aFieldOnGroup2 row2 (-0.235702,  2.7643  )
00939       //aFieldOnGroup2 row3 (-0.235702, -1.2357  )
00940       //aFieldOnGroup2 row4 ( 1.7643  , -0.235702)
00941       //aFieldOnGroup2 row5 ( 0.235702,  2.7357  )
00942     }
00943 
00944   // info about support (Group1)
00945   CPPUNIT_ASSERT(!aFieldOnGroup1->isOnAllElements()); // because we build Group1 so
00946   int nbTypes = aFieldOnGroup1->getNumberOfGeometricTypes();
00947   //CPPUNIT_ASSERT(nbTypes);
00948   CPPUNIT_ASSERT_EQUAL(2, nbTypes);
00949   const int * nbElemsInEachType = aFieldOnGroup1->getNumberOfElements();
00950   const MED_EN::medGeometryElement * aGeomTypes = aFieldOnGroup1->getGeometricTypes();
00951 
00952   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_TRIA3, aGeomTypes[0]);
00953   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_QUAD4, aGeomTypes[1]);
00954 
00955   // GAUSS
00956 
00957   // now we have no gauss localization in aFieldOnGroup1
00958   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
00959   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
00960   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
00961   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
00962 
00963   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3), MEDEXCEPTION);
00964   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3), MEDEXCEPTION);
00965 
00966   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
00967 
00968   // set a gauss localization into aFieldOnGroup1
00969   double cooRef[6] = 
00970     {
00971       1.,1., 2.,4., 3.,9.
00972     }; // xy xy xy
00973   double cooGauss[10] = 
00974     {
00975       7.,7., 6.,6., 5.,5., 4.,3., 2.,1.
00976     }; // x1,y1  x2,y2  x3,y3  x4,y4  x5,y5
00977   double wg[5] = 
00978     {
00979       1., 2., 3., 4., 5.
00980     };
00981   GAUSS_LOCALIZATION<> gl1 ("GL1", MED_EN::MED_TRIA3, /*nGauss*/5, cooRef, cooGauss, wg);
00982 
00983   aFieldOnGroup1->setGaussLocalization(MED_EN::MED_TRIA3, gl1);
00984 
00985   // now we have a gauss localization for MED_TRIA3 type
00986   CPPUNIT_ASSERT_EQUAL(5, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
00987   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
00988   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
00989   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
00990 
00991   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_QUAD4), MEDEXCEPTION);
00992   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_QUAD4), MEDEXCEPTION);
00993 
00994   GAUSS_LOCALIZATION<> gl1Back = aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3);
00995   const GAUSS_LOCALIZATION<> * gl1BackPtr = aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3);
00996 
00997   CPPUNIT_ASSERT(gl1 == gl1Back);
00998   CPPUNIT_ASSERT(gl1 == *gl1BackPtr);
00999 
01000   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
01001 
01002   // sub-support of Group1 on one (first) geometric type
01003   SUPPORT * aSubSupport1 = new SUPPORT();
01004   aSubSupport1->setMesh( aMesh );
01005   aSubSupport1->setName( "Sub-Support 1 of Group1" );
01006   aSubSupport1->setEntity( MED_EN::MED_FACE );
01007 
01008   int nbTypes1 = 1;
01009   int nbElemsInEachType1[1];
01010   nbElemsInEachType1[0] = nbElemsInEachType[0];
01011   int nbElems1 = nbElemsInEachType1[0];
01012   MED_EN::medGeometryElement aGeomTypes1[1];
01013   aGeomTypes1[0] = aGeomTypes[0];
01014   int * anElems1 = new int[nbElems1];
01015   for (int i = 0; i < nbElems1; i++)
01016     {
01017       anElems1[i] = anElems[i];
01018     }
01019 
01020   aSubSupport1->setpartial("Support for sub-field 1 on one type of elements",
01021                            nbTypes1, nbElems1, aGeomTypes1, nbElemsInEachType1, anElems1);
01022 
01023   // extract sub-field on aSubSupport1
01024   FIELD<double, FullInterlace> * aSubField1 = aFieldOnGroup1->extract(aSubSupport1);
01025   CPPUNIT_ASSERT_EQUAL(nbElems1 * /*NumberOfComponents = */2, aSubField1->getValueLength());
01026 
01027   // aSubField1:
01028   // elt\comp |  1 |  2
01029   //--------------------
01030   //  1       | 14 | 12
01031   //  2       | 15 | 11
01032 
01033   // check normL2() and normL1()
01034   FIELD<double>* anAreaField = aMesh->getArea(aSubSupport1);
01035   double area1 = anAreaField->getValueIJ(anElems1[0], 1);
01036   double area2 = anAreaField->getValueIJ(anElems1[1], 1);
01037   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area1, 0.00001);
01038   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area2, 0.00001);
01039 
01040   CPPUNIT_ASSERT_DOUBLES_EQUAL(210.5, aSubField1->normL2(1), 0.00001); // (14*14 + 15*15)/2
01041   CPPUNIT_ASSERT_DOUBLES_EQUAL(132.5, aSubField1->normL2(2), 0.00001); // (12*12 + 11*11)/2
01042   CPPUNIT_ASSERT_DOUBLES_EQUAL(343.0, aSubField1->normL2() , 0.00001); // 210.5 + 132.5
01043 
01044   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.5, aSubField1->normL1(1), 0.00001); // (14 + 15)/2
01045   CPPUNIT_ASSERT_DOUBLES_EQUAL(11.5, aSubField1->normL1(2), 0.00001); // (12 + 11)/2
01046   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1() , 0.00001); // 14.5 + 11.5
01047 
01048   double aNewArea [2] =
01049     {
01050       1., 0.
01051     }; // only first element will be taken into account
01052   anAreaField->setValue(aNewArea);
01053 
01054   CPPUNIT_ASSERT_DOUBLES_EQUAL(196.0, aSubField1->normL2(1, anAreaField), 0.00001); // 14*14/1
01055   CPPUNIT_ASSERT_DOUBLES_EQUAL(144.0, aSubField1->normL2(2, anAreaField), 0.00001); // 12*12/1
01056   CPPUNIT_ASSERT_DOUBLES_EQUAL(340.0, aSubField1->normL2(anAreaField) , 0.00001); // 196 + 144
01057 
01058   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.0, aSubField1->normL1(1, anAreaField), 0.00001); // 14/1
01059   CPPUNIT_ASSERT_DOUBLES_EQUAL(12.0, aSubField1->normL1(2, anAreaField), 0.00001); // 12/1
01060   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1(anAreaField) , 0.00001); // 14 + 12
01061 
01062   // check normL2() on nodal field (issue 0020120)
01063   {
01064     // read nodal field from pointe.med
01065     string filename  = getResourceFile("pointe.med");
01066     string fieldname = "fieldnodedouble";
01067     string meshname  = "maa1";
01068     FIELD<double> *nodalField=new FIELD<double>( MED_DRIVER, filename, fieldname);
01069     MESH *mesh=new MESH( MED_DRIVER, filename, meshname);
01070     nodalField->getSupport()->setMesh( mesh );
01071 
01072     // make a field on the nodes of first cell only
01073     SUPPORT *oneCellNodesSup=new SUPPORT;
01074     oneCellNodesSup->setMesh(mesh);
01075     oneCellNodesSup->setName("Sub-Support of nodes of 1 cell");
01076     oneCellNodesSup->setEntity(MED_NODE);
01077     int NumberOfElements[] =
01078       {
01079         mesh->getTypes(MED_CELL)[0]%100
01080       };
01081     medGeometryElement GeometricType[] =
01082       {
01083         MED_POINT1
01084       };
01085     oneCellNodesSup->setpartial("Support for sub-field of one cell nodes",
01086                                 /*NumberOfGeometricType=*/1,
01087                                 /*TotalNumberOfElements=*/ *NumberOfElements,
01088                                 GeometricType,
01089                                 NumberOfElements,
01090                                 /*NumberValue=*/ mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS ));
01091     FIELD<double, FullInterlace> * oneCellNodesField = nodalField->extract( oneCellNodesSup );
01092     oneCellNodesSup->removeReference();
01093     // compute normL2 by avarage nodal value on the cell
01094 
01095     const SUPPORT *allCellsSupport=mesh->getSupportOnAll( MED_CELL );
01096     FIELD<double>* volumeField = mesh->getVolume(allCellsSupport);
01097     // mdump output:
01098     // - Mailles de type MED_TETRA4 : 
01099     //  [     1 ] :          1          2          3          6
01100     // - Valeurs :
01101     // | 1.000000 | 1.000000 | 1.000000 | 2.000000 | 2.000000 | 2.000000 | ...
01102     const double cellVal = ( 1.000000 + 1.000000 + 1.000000 + 2.000000 ) / 4;
01103     const double vol = abs( volumeField->getValueIJ( 1, 1 ));
01104     const double norm = cellVal*cellVal*vol/vol; // v*v*vol/totVol;
01105     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(), 0.000001);
01106     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(1), 0.000001);
01107     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(volumeField), 0.000001);
01108     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(1, volumeField), 0.000001);
01109 
01110     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL1(), MEDEXCEPTION);
01111     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL1(1), MEDEXCEPTION);
01112     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(nodalField), MEDEXCEPTION);
01113     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(anAreaField), MEDEXCEPTION);
01114     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(aSubField1), MEDEXCEPTION);
01115     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(aFieldOnGroup1), MEDEXCEPTION);
01116     nodalField->removeReference();
01117     volumeField->removeReference();
01118     oneCellNodesField->removeReference();
01119     mesh->removeReference();
01120   }
01121 
01122   // double integral(SUPPORT*) - mantis issue 0020460:
01123   // [CEA 352] Integral calculus on all field or on a field subarea (groupe or family)
01124   {
01125     // make 2D grid 2x2 with steps along axis 1.0, 2.0
01126     const int dim = 2;
01127     double coord[3] =
01128       {
01129         0.0, 1.0, 3.0
01130       };
01131     vector< vector<double> > xyz_array( dim, vector<double>( coord, coord + 3 ));
01132     vector<string> coord_name( dim, "coord_name");
01133     vector<string> coord_unit( dim, "m");
01134     MESH *mesh= const_cast<MESH*>( GRID( xyz_array, coord_name, coord_unit ).convertInMESH());
01135 
01136     // make supports on the grid
01137     const SUPPORT *supOnAll=mesh->getSupportOnAll(MED_CELL);
01138     SUPPORT *sup123=new SUPPORT;
01139     sup123->setMesh(mesh);
01140     sup123->setName("sup123");
01141     int nbEl123[] =
01142       {
01143         3
01144       };
01145     int elems123[] =
01146       {
01147         1,2,3
01148       };
01149     SUPPORT *sup12=new SUPPORT;
01150     sup12->setMesh(mesh);
01151     sup12->setName("sup12");
01152     int nbEl12 [] =
01153       {
01154         2
01155       };
01156     int elems12 [] =
01157       {
01158         1,2
01159       };
01160     SUPPORT *sup34=new SUPPORT;
01161     sup34->setMesh(mesh);
01162     sup34->setName("sup34");
01163     int nbEl34 [] =
01164       {
01165         2
01166       };
01167     int elems34 [] =
01168       {
01169         3,4
01170       };
01171     const int nbGeomTypes = 1;
01172     const medGeometryElement * geomType = mesh->getTypes(MED_EN::MED_CELL);
01173     mesh->removeReference();
01174     sup123->setpartial("test", nbGeomTypes, *nbEl123, geomType, nbEl123, elems123 );
01175     sup12->setpartial("test", nbGeomTypes, *nbEl12 , geomType, nbEl12 , elems12 );
01176     sup34->setpartial("test", nbGeomTypes, *nbEl34 , geomType, nbEl34 , elems34 );
01177 
01178     // make vectorial fields with values of i-th elem { i, i*10, i*100 }
01179     const int nbComp = 3, nbElems = 4;
01180     const int mult[nbComp] =
01181       {
01182         1, 10, 100
01183       };
01184     FIELD<int,NoInterlaceByType> *fAllNoTy=new FIELD<int,NoInterlaceByType>(supOnAll, nbComp), *f12NoTy=new FIELD<int,NoInterlaceByType>(sup12, nbComp);
01185     FIELD<int,NoInterlace>       *fAllNo=new FIELD<int,NoInterlace>(supOnAll, nbComp), *f12No=new FIELD<int,NoInterlace>(sup12, nbComp);
01186     FIELD<int,FullInterlace>     *fAllFull=new FIELD<int,FullInterlace>(supOnAll, nbComp), *f12Full=new FIELD<int,FullInterlace>(sup12, nbComp);
01187     int i, j;
01188     for ( i = 1; i <= nbElems; ++i )
01189       for ( j = 1; j <= nbComp; ++j )
01190         {
01191           fAllFull->setValueIJ( i, j, i * mult[j-1]);
01192           fAllNoTy->setValueIJ( i, j, i * mult[j-1]);
01193           fAllNo  ->setValueIJ( i, j, i * mult[j-1]);
01194           if ( i < 3 )
01195             {
01196               f12Full->setValueIJ( i, j, i * mult[j-1]);
01197               f12NoTy->setValueIJ( i, j, i * mult[j-1]);
01198               f12No  ->setValueIJ( i, j, i * mult[j-1]);
01199             }
01200         }
01201     // Test
01202     double preci = 1e-18, integral;
01203     // Integral = SUM( area * (i*1 + i*10 + i*100)) == 111 * SUM( area * i )
01204     // elem area: { 1, 2, 2, 4 }
01205     integral = 111*( 1*1 + 2*2 + 2*3 + 4*4 );
01206     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNoTy->integral(), preci );
01207     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNo  ->integral(), preci );
01208     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllFull->integral(), preci );
01209     integral = 111*( 1*1 + 2*2 );
01210     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(), preci );
01211     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(), preci );
01212     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(), preci );
01213 
01214     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNoTy->integral(sup12), preci );
01215     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNo  ->integral(sup12), preci );
01216     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllFull->integral(sup12), preci );
01217 
01218     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(sup12), preci );
01219     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(sup12), preci );
01220     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(sup12), preci );
01221     sup12->removeReference();
01222     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(sup123), preci );
01223     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(sup123), preci );
01224     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(sup123), preci );
01225     integral = 111*( 1*1 + 2*2 + 2*3 );
01226     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNoTy->integral(sup123), preci );
01227     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNo  ->integral(sup123), preci );
01228     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllFull->integral(sup123), preci );
01229     fAllNoTy->removeReference();
01230     fAllNo->removeReference();
01231     sup123->removeReference();
01232     fAllFull->removeReference();
01233     integral = 0;
01234     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(sup34), preci );
01235     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(sup34), preci );
01236     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(sup34), preci );
01237     sup34->removeReference();
01238     f12NoTy->removeReference();
01239     f12No->removeReference();
01240     f12Full->removeReference();
01241   }
01242 
01243   // applyPow
01244   aSubField1->applyPow(2.);
01245   CPPUNIT_ASSERT_DOUBLES_EQUAL(196., aSubField1->getValueIJ(anElems1[0], 1), 0.000001); // 14*14
01246   CPPUNIT_ASSERT_DOUBLES_EQUAL(144., aSubField1->getValueIJ(anElems1[0], 2), 0.000001); // 12*12
01247   CPPUNIT_ASSERT_DOUBLES_EQUAL(225., aSubField1->getValueIJ(anElems1[1], 1), 0.000001); // 15*15
01248   CPPUNIT_ASSERT_DOUBLES_EQUAL(121., aSubField1->getValueIJ(anElems1[1], 2), 0.000001); // 11*11
01249 
01250   // setArray (NoGauss)
01251   MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array * aNewArrayNoGauss =
01252     new MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array(/*dim*/2, /*nbelem*/2);
01253   aNewArrayNoGauss->setIJ(1, 1, 4.);
01254   aNewArrayNoGauss->setIJ(1, 2, 2.);
01255   aNewArrayNoGauss->setIJ(2, 1, 5.);
01256   aNewArrayNoGauss->setIJ(2, 2, 1.);
01257   aSubField1->setArray(aNewArrayNoGauss);
01258   // no need to delete aNewArrayNoGauss, because it will be deleted
01259   // in destructor or in deallocValue() method of aSubField1
01260 
01261   CPPUNIT_ASSERT_DOUBLES_EQUAL(4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01262   CPPUNIT_ASSERT_DOUBLES_EQUAL(2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01263   CPPUNIT_ASSERT_DOUBLES_EQUAL(5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01264   CPPUNIT_ASSERT_DOUBLES_EQUAL(1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01265 
01266   // setRow
01267   double row[2] =
01268     {
01269       -1., -3.
01270     };
01271   aSubField1->setRow(anElems1[0], row);
01272   CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01273   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01274   CPPUNIT_ASSERT_DOUBLES_EQUAL( 5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01275   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01276   // out of range
01277   CPPUNIT_ASSERT_THROW(aSubField1->setRow(3, row), MEDEXCEPTION);
01278 
01279   // setColumn
01280   double col[2] =
01281     {
01282       -7., -9.
01283     };
01284   aSubField1->setColumn(1, col);
01285   CPPUNIT_ASSERT_DOUBLES_EQUAL(-7., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01286   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01287   CPPUNIT_ASSERT_DOUBLES_EQUAL(-9., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01288   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01289   // out of range
01290   CPPUNIT_ASSERT_THROW(aSubField1->setColumn(3, col), MEDEXCEPTION);
01291 
01292   // setArray (Gauss)
01293   {
01294     int nbelgeoc[2] =
01295       {
01296         1, 3
01297       }; // 3 - 1 = two elements for the first (and the only) type
01298     int nbgaussgeo[2] =
01299       {
01300         0, 1
01301       }; // one gauss point per each element
01302     MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array * aNewArrayGauss =
01303       new MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array
01304       (/*dim*/2, /*nbelem*/2, /*nbtypegeo*/1, /*nbelgeoc*/nbelgeoc, /*nbgaussgeo*/nbgaussgeo);
01305 
01306     aNewArrayGauss->setIJ(1, 1, -4.);
01307     aNewArrayGauss->setIJ(1, 2, -2.);
01308     aNewArrayGauss->setIJ(2, 1, -5.);
01309     aNewArrayGauss->setIJ(2, 2, -1.);
01310 
01311     aNewArrayGauss->setIJK(1, 1, 1, -4.);
01312     aNewArrayGauss->setIJK(1, 2, 1, -2.);
01313     aNewArrayGauss->setIJK(2, 1, 1, -5.);
01314     aNewArrayGauss->setIJK(2, 2, 1, -1.);
01315 
01316     aSubField1->setArray(aNewArrayGauss);
01317     // no need to delete aNewArrayGauss, because it will be deleted
01318     // in destructor or in deallocValue() method of aSubField1
01319 
01320     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01321     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01322     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01323     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01324 
01325     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJK(anElems1[0], 1, 1), 0.000001);
01326     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJK(anElems1[0], 2, 1), 0.000001);
01327     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJK(anElems1[1], 1, 1), 0.000001);
01328     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJK(anElems1[1], 2, 1), 0.000001);
01329   }
01330 
01331   // alloc/dealloc; compatibility of new size with support
01332   try
01333     {
01334       aSubField1->deallocValue();
01335       aSubField1->allocValue(/*NumberOfComponents*/2, /*LengthValue*/5);
01336       //#ifdef ENABLE_FAULTS
01337       // (BUG) No compatibility between Support and allocated value
01338       //aSubField1->normL1();
01339       CPPUNIT_ASSERT_THROW(aSubField1->normL1(),MEDEXCEPTION);
01340       //#endif
01341       //#ifdef ENABLE_FORCED_FAILURES
01342       // TODO: FIX to throw an EXCEPTION
01343       //    CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
01344       //#endif
01345     }
01346   catch (MEDEXCEPTION & ex)
01347     {
01348       // normal behaviour
01349     }
01350   catch (...)
01351     {
01352       CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
01353     }
01354 
01355   // check that aFieldOnGroup1 is not changed after aSubField1 modifications
01356   for (int i = 1; i <= nbVals; i++)
01357     {
01358       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01359       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01360     }
01361 
01362   // reset aFieldOnGroup2 values for simple control of operators results
01363   for (int i = 1; i <= nbVals; i++)
01364     {
01365       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, i*i);
01366       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, -i*i);
01367     }
01368 
01369   int len = aFieldOnGroup1->getValueLength();
01370   const double * val1 = aFieldOnGroup1->getValue();
01371   const double * val2 = aFieldOnGroup2->getValue();
01372   const double * val_res;
01373 
01374   // operators and add, sub, mul, div
01375 
01376   // +
01377   FIELD<double> *aSum = *aFieldOnGroup1 + *aFieldOnGroup2;
01378   aSum->setName(aFieldOnGroup1->getName());
01379   aSum->setDescription(aFieldOnGroup1->getDescription());
01380   compareField_(aFieldOnGroup1, aSum, true, true);
01381   val_res = aSum->getValue();
01382   for (int i = 0; i < len; i++)
01383     {
01384       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
01385     }
01386   aSum->removeReference();
01387   // -
01388   FIELD<double> *aDifference = *aFieldOnGroup1 - *aFieldOnGroup2;
01389   aDifference->setName(aFieldOnGroup1->getName());
01390   aDifference->setDescription(aFieldOnGroup1->getDescription());
01391   compareField_(aFieldOnGroup1, aDifference, true, true);
01392   val_res = aDifference->getValue();
01393   for (int i = 0; i < len; i++)
01394     {
01395       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
01396     }
01397   aDifference->removeReference();
01398   // - (unary)
01399   FIELD<double> *aNegative = - *aFieldOnGroup1;
01400   aNegative->setName(aFieldOnGroup1->getName());
01401   aNegative->setDescription(aFieldOnGroup1->getDescription());
01402   compareField_(aFieldOnGroup1, aNegative, true, true);
01403   val_res = aNegative->getValue();
01404   for (int i = 0; i < len; i++)
01405     {
01406       CPPUNIT_ASSERT_DOUBLES_EQUAL(- val1[i], val_res[i], 0.000001);
01407     }
01408   aNegative->removeReference();
01409   // *
01410   FIELD<double> *aProduct = (*aFieldOnGroup1) * (*aFieldOnGroup2);
01411   aProduct->setName(aFieldOnGroup1->getName());
01412   aProduct->setDescription(aFieldOnGroup1->getDescription());
01413   compareField_(aFieldOnGroup1, aProduct, true, true);
01414   val_res = aProduct->getValue();
01415   for (int i = 0; i < len; i++)
01416     {
01417       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
01418     }
01419   aProduct->removeReference();
01420   // /
01421   FIELD<double> *aQuotient = *aFieldOnGroup1 / *aFieldOnGroup2;
01422   aQuotient->setName(aFieldOnGroup1->getName());
01423   aQuotient->setDescription(aFieldOnGroup1->getDescription());
01424   compareField_(aFieldOnGroup1, aQuotient, true, true);
01425   val_res = aQuotient->getValue();
01426   for (int i = 0; i < len; i++)
01427     {
01428       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
01429     }
01430   aQuotient->removeReference();
01431   double val22 = aFieldOnGroup2->getValueIJ(anElems[2], 2);
01432   aFieldOnGroup2->setValueIJ(anElems[2], 2, 0.);
01433 
01434   CPPUNIT_ASSERT_THROW((*aFieldOnGroup1 / *aFieldOnGroup2), MEDEXCEPTION);
01435   //#ifdef ENABLE_FORCED_FAILURES
01436   // (BUG) is it up to user to control validity of data to avoid division on zero?
01437   // YES: USER SHOULD CARE OF IT
01438   //CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup2, MEDEXCEPTION);
01439   //#endif
01440   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2)->removeReference(), MEDEXCEPTION);
01441   CPPUNIT_ASSERT_THROW(FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2)->removeReference(), MEDEXCEPTION);
01442 
01443   // restore value
01444   aFieldOnGroup2->setValueIJ(anElems[2], 2, val22);
01445 
01446   // restore values
01447   for (int i = 1; i <= nbVals; i++)
01448     {
01449       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, 13 + i);
01450       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, 13 - i);
01451     }
01452 
01453   // static methods
01454   FIELD<double> * aPtr;
01455 
01456   // add
01457   aPtr = FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup2);
01458   aPtr->setName(aFieldOnGroup1->getName());
01459   aPtr->setDescription(aFieldOnGroup1->getDescription());
01460   compareField_(aFieldOnGroup1, aPtr, true, true);
01461   val_res = aPtr->getValue();
01462   for (int i = 0; i < len; i++)
01463     {
01464       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
01465     }
01466   aPtr->removeReference();
01467 
01468   // sub
01469   aPtr = FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2);
01470   aPtr->setName(aFieldOnGroup1->getName());
01471   aPtr->setDescription(aFieldOnGroup1->getDescription());
01472   compareField_(aFieldOnGroup1, aPtr, true, true);
01473   val_res = aPtr->getValue();
01474   for (int i = 0; i < len; i++)
01475     {
01476       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
01477     }
01478   aPtr->removeReference();
01479 
01480   // mul
01481   aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
01482   aPtr->setName(aFieldOnGroup1->getName());
01483   aPtr->setDescription(aFieldOnGroup1->getDescription());
01484   compareField_(aFieldOnGroup1, aPtr, true, true);
01485   val_res = aPtr->getValue();
01486   for (int i = 0; i < len; i++)
01487     {
01488       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
01489     }
01490   aPtr->removeReference();
01491 
01492   // div
01493   aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
01494   aPtr->setName(aFieldOnGroup1->getName());
01495   aPtr->setDescription(aFieldOnGroup1->getDescription());
01496   compareField_(aFieldOnGroup1, aPtr, true, true);
01497   val_res = aPtr->getValue();
01498   for (int i = 0; i < len; i++)
01499     {
01500       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
01501     }
01502   aPtr->removeReference();
01503 
01504   // addDeep
01505   aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01506   aPtr->setName(aFieldOnGroup1->getName());
01507   aPtr->setDescription(aFieldOnGroup1->getDescription());
01508   compareField_(aFieldOnGroup1, aPtr, true, true);
01509   val_res = aPtr->getValue();
01510   for (int i = 0; i < len; i++)
01511     {
01512       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
01513     }
01514   aPtr->removeReference();
01515 
01516   // subDeep
01517   aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01518   aPtr->setName(aFieldOnGroup1->getName());
01519   aPtr->setDescription(aFieldOnGroup1->getDescription());
01520   compareField_(aFieldOnGroup1, aPtr, true, true);
01521   val_res = aPtr->getValue();
01522   for (int i = 0; i < len; i++)
01523     {
01524       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
01525     }
01526   aPtr->removeReference();
01527 
01528   // mulDeep
01529   aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01530   aPtr->setName(aFieldOnGroup1->getName());
01531   aPtr->setDescription(aFieldOnGroup1->getDescription());
01532   compareField_(aFieldOnGroup1, aPtr, true, true);
01533   val_res = aPtr->getValue();
01534   for (int i = 0; i < len; i++)
01535     {
01536       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
01537     }
01538   aPtr->removeReference();
01539 
01540   // divDeep
01541   aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01542   aPtr->setName(aFieldOnGroup1->getName());
01543   aPtr->setDescription(aFieldOnGroup1->getDescription());
01544   compareField_(aFieldOnGroup1, aPtr, true, true);
01545   val_res = aPtr->getValue();
01546   for (int i = 0; i < len; i++)
01547     {
01548       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
01549     }
01550   aPtr->removeReference();
01551 
01552   // +=
01553   *aFieldOnGroup1 += *aFieldOnGroup2;
01554   for (int i = 1; i <= nbVals; i++)
01555     {
01556       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i + i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01557       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i - i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01558     }
01559 
01560   // -=
01561   *aFieldOnGroup1 -= *aFieldOnGroup2;
01562   for (int i = 1; i <= nbVals; i++)
01563     {
01564       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01565       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01566     }
01567 
01568   // *=
01569   *aFieldOnGroup1 *= *aFieldOnGroup2;
01570   for (int i = 1; i <= nbVals; i++)
01571     {
01572       CPPUNIT_ASSERT_DOUBLES_EQUAL( (13 + i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01573       CPPUNIT_ASSERT_DOUBLES_EQUAL(-(13 - i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01574     }
01575 
01576   // /=
01577   *aFieldOnGroup1 /= *aFieldOnGroup2;
01578   for (int i = 1; i <= nbVals; i++)
01579     {
01580       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01581       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01582     }
01583 
01584   // check case of different operands: support
01585   MESH * aMeshOneMore = MEDMEMTest_createTestMesh();
01586   const GROUP* aGroupOneMore = aMeshOneMore->getGroup(MED_EN::MED_FACE, 1);
01587   FIELD<double> * aFieldOnGroup3 =
01588     createFieldOnGroup<double>(aMeshOneMore, aGroupOneMore, "Test_Diff_Mesh", "test");
01589   for (int i = 1; i <= nbVals; i++)
01590     {
01591       aFieldOnGroup3->setValueIJ(anElems[i-1], 1, 2*i);
01592       aFieldOnGroup3->setValueIJ(anElems[i-1], 2, 3*i);
01593     }
01594   const double * val3 = aFieldOnGroup3->getValue();
01595 
01596   //CPPUNIT_ASSERT_NO_THROW();
01597   try
01598     {
01599       aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01600       aPtr->setName(aFieldOnGroup1->getName());
01601       aPtr->setDescription(aFieldOnGroup1->getDescription());
01602       compareField_(aFieldOnGroup1, aPtr, true, true);
01603       val_res = aPtr->getValue();
01604       for (int i = 0; i < len; i++)
01605         {
01606           CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val3[i], val_res[i], 0.000001);
01607         }
01608       aPtr->removeReference();
01609 
01610       aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01611       aPtr->removeReference();
01612       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01613       aPtr->removeReference();
01614       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01615       aPtr->removeReference();
01616     }
01617   catch (MEDEXCEPTION & ex)
01618     {
01619       CPPUNIT_FAIL(ex.what());
01620     }
01621   catch (...)
01622     {
01623       CPPUNIT_FAIL("Unknown exception in FIELD::xxxDeep()");
01624     }
01625 
01626   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01627   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01628   CPPUNIT_ASSERT_THROW(FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01629   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01630 
01631   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup3, MEDEXCEPTION);
01632   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup3, MEDEXCEPTION);
01633   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 * *aFieldOnGroup3, MEDEXCEPTION);
01634   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup3, MEDEXCEPTION);
01635 
01636   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup3, MEDEXCEPTION);
01637   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup3, MEDEXCEPTION);
01638   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 *= *aFieldOnGroup3, MEDEXCEPTION);
01639   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup3, MEDEXCEPTION);
01640 
01641   // check case of different operands: MEDComponentsUnits
01642   aFieldOnGroup1->setMEDComponentUnit(1, "unit3");
01643 
01644   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup2, MEDEXCEPTION);
01645   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup2, MEDEXCEPTION);
01646   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup2, MEDEXCEPTION);
01647   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup2, MEDEXCEPTION);
01648   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01649   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01650   CPPUNIT_ASSERT_THROW(FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01651   CPPUNIT_ASSERT_THROW(FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01652 
01653   //CPPUNIT_ASSERT_NO_THROW();
01654   try
01655     {
01656       aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
01657       aPtr->removeReference();
01658       aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
01659       aPtr->removeReference();
01660       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01661       aPtr->removeReference();
01662       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01663       aPtr->removeReference();
01664 
01665       *aFieldOnGroup1 *= *aFieldOnGroup2;
01666       *aFieldOnGroup1 /= *aFieldOnGroup2;
01667 
01668       FIELD<double> *aPr = *aFieldOnGroup1 * *aFieldOnGroup2;
01669       FIELD<double> *aQu = *aFieldOnGroup1 / *aFieldOnGroup2;
01670       aPr->removeReference();
01671       aQu->removeReference();
01672     }
01673   catch (MEDEXCEPTION & ex)
01674     {
01675       CPPUNIT_FAIL(ex.what());
01676     }
01677   catch (...)
01678     {
01679       CPPUNIT_FAIL("Unknown exception");
01680     }
01681 
01682   // restore MED units
01683   aFieldOnGroup1->setMEDComponentUnit(1, "unit1");
01684 
01685   // check case of different operands: valueType
01686   //FIELD<int> * aFieldOnGroup4 =
01687   //  createFieldOnGroup<int>(aMeshOneMore, aGroupOneMore, "Test_Diff_Mesh", "test");
01688   //
01689   //CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup4, *aFieldOnGroup3), MEDEXCEPTION);
01690   //CPPUNIT_ASSERT_THROW(*aFieldOnGroup4 - *aFieldOnGroup3, MEDEXCEPTION);
01691   //CPPUNIT_ASSERT_THROW(*aFieldOnGroup4 *= *aFieldOnGroup3, MEDEXCEPTION);
01692   //delete aFieldOnGroup4;
01693 
01694   // check case of different operands: numberOfComponents
01695   //#ifdef ENABLE_FAULTS
01696   // (BUG) Cannot allocate value of higher dimension because of _componentsTypes reinitialization
01697   aFieldOnGroup1->deallocValue();
01698   //CPPUNIT_ASSERT_THROW(aFieldOnGroup1->allocValue(/*dim*/5), MEDEXCEPTION);
01699   CPPUNIT_ASSERT_NO_THROW(aFieldOnGroup1->allocValue(/*dim*/5));
01700   //#endif
01701   //#ifdef ENABLE_FORCED_FAILURES
01702   // YES THERE MUST BE AN EXCEPTION
01703   //   CPPUNIT_FAIL("Segmentation fault on attempt to allocate value of higher dimension."
01704   //                " Must be MEDEXCEPTION instead. And on attempt to change nb.components"
01705   //                " must be the same behaviour.");
01706   //#endif
01707   aFieldOnGroup1->setNumberOfComponents(5);
01708 
01709   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01710   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 * *aFieldOnGroup2, MEDEXCEPTION);
01711   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup2, MEDEXCEPTION);
01712 
01713   // check case of different operands: numberOfValues
01714   aFieldOnGroup1->deallocValue();
01715   aFieldOnGroup1->allocValue(2, nbVals + 1);
01716   // be carefull: aFieldOnGroup1 reallocated and contains random values
01717 
01718   CPPUNIT_ASSERT_THROW(FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01719   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup2, MEDEXCEPTION);
01720   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup2, MEDEXCEPTION);
01721 
01722   // restore aFieldOnGroup1
01723   aFieldOnGroup1->deallocValue();
01724   aFieldOnGroup1->allocValue(2, nbVals);
01725   // be carefull: aFieldOnGroup1 reallocated and contains random values
01726 
01727   aSubSupport1->removeReference();
01728   delete [] anElems1;
01729 
01730   aScalarProduct->removeReference();
01731   aSubField1->removeReference();
01732   anAreaField->removeReference();
01733   barycenter->removeReference();
01734   aFieldOnGroup1->removeReference();
01735   aFieldOnGroup2->removeReference();
01736   aFieldOnGroup3->removeReference();
01737 
01738   aMesh->removeReference();
01739   aMeshOneMore->removeReference();
01740 
01742   // TEST 5: Drivers //
01744   testDrivers();
01745 }
01746 
01752 void MEDMEMTest::testFieldConvert()
01753 {
01754   // create an empty integer field 2x10
01755   FIELD<int, FullInterlace> * aField_FING = new FIELD<int, FullInterlace>();
01756 
01757   aField_FING->setName("Field_FING");
01758   aField_FING->setDescription("Field full interlace no gauss");
01759 
01760   aField_FING->setNumberOfComponents(2);
01761   aField_FING->setNumberOfValues(10);
01762 
01763   string aCompsNames[2] =
01764     {
01765       "Pos", "Neg"
01766     };
01767   string aCompsDescs[2] =
01768     {
01769       "+", "-"
01770     };
01771   string aMEDCompsUnits[2] =
01772     {
01773       "unit1", "unit2"
01774     };
01775   UNIT   aCompsUnits[2];
01776 
01777   aCompsUnits[0] = UNIT("u1", "descr1");
01778   aCompsUnits[1] = UNIT("u2", "descr2");
01779 
01780   aField_FING->setComponentsNames(aCompsNames);
01781   aField_FING->setComponentsDescriptions(aCompsDescs);
01782   aField_FING->setMEDComponentsUnits(aMEDCompsUnits);
01783   aField_FING->setComponentsUnits(aCompsUnits);
01784 
01785   // check UNITs (for testField())
01786   const UNIT * aCompsUnitsBack = aField_FING->getComponentsUnits();
01787   CPPUNIT_ASSERT(aCompsUnits[0].getName() == aCompsUnitsBack[0].getName());
01788   CPPUNIT_ASSERT(aCompsUnits[1].getName() == aCompsUnitsBack[1].getName());
01789 
01790   const UNIT * aCompUnitBack1 = aField_FING->getComponentUnit(1);
01791   const UNIT * aCompUnitBack2 = aField_FING->getComponentUnit(2);
01792   CPPUNIT_ASSERT(aCompsUnits[0].getName() == aCompUnitBack1->getName());
01793   CPPUNIT_ASSERT(aCompsUnits[1].getName() == aCompUnitBack2->getName());
01794 
01795   // create one more field by copy
01796   FIELD<int, FullInterlace> * aField_FIGG = new FIELD<int, FullInterlace>(*aField_FING);
01797 
01798   // values
01799   int values_FING[20] =
01800     {
01801       7,- 7, 14,-14, 21,-21, 28,-28, 35,-35,
01802       42,-42, 49,-49, 56,-56, 63,-63, 70,-70
01803     };
01804 
01806   // TEST 1: NoGauss //
01808 
01809   MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array * anArray_FING =
01810     new MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array
01811     (values_FING, /*dim*/2, /*nbelem*/10, /*shallowCopy*/false, /*ownershipOfValues*/false);
01812   aField_FING->setArray(anArray_FING);
01813   // no need to delete anArray_FING, because it will be deleted in destructor of aField_FING
01814 
01815   // 1. FullInterlace -> NoInterlace
01816   FIELD<int, NoInterlace> * aField_NING = FieldConvert(*aField_FING);
01817   const int * values_NING = aField_NING->getValue();
01818 
01819   for (int i = 0; i < 10; i++)
01820     {
01821       for (int j = 0; j < 2; j++)
01822         {
01823           CPPUNIT_ASSERT_EQUAL(values_FING[2*i + j], values_NING[10*j + i]);
01824         }
01825     }
01826 
01827   // 2. NoInterlace -> FullInterlace
01828   FIELD<int, FullInterlace> * aField_FING_conv = FieldConvert(*aField_NING);
01829   const int * values_FING_conv = aField_FING_conv->getValue();
01830 
01831   for (int i = 0; i < 10; i++)
01832     {
01833       for (int j = 0; j < 2; j++)
01834         {
01835           CPPUNIT_ASSERT_EQUAL(values_FING_conv[2*i + j], values_FING[2*i + j]);
01836           CPPUNIT_ASSERT_EQUAL(values_FING_conv[2*i + j], values_NING[10*j + i]);
01837         }
01838     }
01839 
01840   aField_FING->removeReference();
01841   aField_NING->removeReference();
01842   aField_FING_conv->removeReference();
01843 
01845   // TEST 2: Gauss //
01847   int nbelgeoc[2] =
01848     {
01849       1, 11
01850     };
01851   int nbgaussgeo[2] =
01852     {
01853       -1, 1
01854     };
01855   MEDMEM_ArrayInterface<int,FullInterlace,Gauss>::Array * anArray_FIGG =
01856     new MEDMEM_ArrayInterface<int,FullInterlace,Gauss>::Array
01857     (values_FING, /*dim*/2, /*nbelem*/10, /*nbtypegeo*/1, /*nbelgeoc*/nbelgeoc,
01858      /*nbgaussgeo*/nbgaussgeo, /*shallowCopy*/false, /*ownershipOfValues*/false);
01859   aField_FIGG->setArray(anArray_FIGG);
01860   // no need to delete anArray_FIGG, because it will be deleted in destructor of aField_FIGG
01861 
01862   // 1. FullInterlace -> NoInterlace
01863   FIELD<int, NoInterlace> * aField_NIGG = FieldConvert(*aField_FIGG);
01864   const int * values_NIGG = aField_NIGG->getValue();
01865 
01866   for (int i = 0; i < 10; i++)
01867     {
01868       for (int j = 0; j < 2; j++)
01869         {
01870           CPPUNIT_ASSERT_EQUAL(values_FING[2*i + j], values_NIGG[10*j + i]);
01871         }
01872     }
01873 
01874   // 2. NoInterlace -> FullInterlace
01875   FIELD<int, FullInterlace> * aField_FIGG_conv = FieldConvert(*aField_NIGG);
01876   const int * values_FIGG_conv = aField_FIGG_conv->getValue();
01877 
01878   for (int i = 0; i < 10; i++)
01879     {
01880       for (int j = 0; j < 2; j++)
01881         {
01882           CPPUNIT_ASSERT_EQUAL(values_FIGG_conv[2*i + j], values_FING[2*i + j]);
01883           CPPUNIT_ASSERT_EQUAL(values_FIGG_conv[2*i + j], values_NIGG[10*j + i]);
01884         }
01885     }
01886 
01887   aField_FIGG->removeReference();
01888   aField_NIGG->removeReference();
01889   aField_FIGG_conv->removeReference();
01890   {
01891     // create an empty integer field 2x10
01892     FIELD<int, FullInterlace> * aField = new FIELD<int, FullInterlace>();
01893 
01894     aField->setName("aField");
01895     aField->setDescription("Field full interlace no gauss");
01896 
01897     aField->setNumberOfComponents(2);
01898     aField->setNumberOfValues(10);
01899 
01900     aField->setComponentsNames(aCompsNames);
01901     aField->setComponentsDescriptions(aCompsDescs);
01902     aField->setMEDComponentsUnits(aMEDCompsUnits);
01903 
01904     MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array * anArray =
01905       new MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array
01906       (values_FING, /*dim*/2, /*nbelem*/10, /*shallowCopy*/false, /*ownershipOfValues*/false);
01907     aField->setArray(anArray);
01908     // no need to delete anArray, because it will be deleted in destructor of aField
01909 
01910     FIELD<int, NoInterlace> * aField_conv = FieldConvert(*aField);
01911     aField->removeReference();
01912     CPPUNIT_ASSERT(aField_conv);
01913     aField_conv->removeReference();
01914   }
01915 }
01916 
01917 //================================================================================
01923 //================================================================================
01924 
01925 void MEDMEMTest::testReadFieldOnNodesAndCells()
01926 {
01927   const string outfile = makeTmpFile("field_on_nodes_and_cells.med");
01928   const string fieldName = "field on NODEs and CELLs";
01929 
01930   MEDMEMTest_TmpFilesRemover aTmpFilesRemover;
01931   aTmpFilesRemover.Register( outfile );
01932 
01933   // write file with a field on NODEs and CELLs
01934 
01935   using namespace MED_EN;
01936 
01937   MEDMEM::MESH *mesh=MEDMEMTest_createTestMesh();
01938   int drv = mesh->addDriver( MED_DRIVER, outfile, mesh->getName());
01939   mesh->write(drv);
01940 
01941   const SUPPORT *supportOnCells=mesh->getSupportOnAll( MED_CELL );
01942   const SUPPORT *supportOnNodes=mesh->getSupportOnAll( MED_NODE );
01943   supportOnCells->addReference();
01944   supportOnNodes->addReference();
01945   mesh->removeReference();
01946   int numberOfCells = supportOnCells->getNumberOfElements(MED_ALL_ELEMENTS);
01947   int numberOfNodes = supportOnNodes->getNumberOfElements(MED_ALL_ELEMENTS);
01948 
01949   PointerOf<double> cellValues( numberOfCells ), nodeValues( numberOfNodes );
01950   for ( int i = 0; i < numberOfCells; ++i ) cellValues[i] = i;
01951   for ( int i = 0; i < numberOfNodes; ++i ) nodeValues[i] = -i;
01952 
01953   FIELD<double> *wrFieldOnCells=new FIELD<double>(supportOnCells,1);
01954   wrFieldOnCells->setName(fieldName);
01955   wrFieldOnCells->setComponentName(1,"Vx");
01956   wrFieldOnCells->setComponentDescription(1,"comp1");
01957   wrFieldOnCells->setMEDComponentUnit(1,"unit1");
01958   wrFieldOnCells->setValue( cellValues );
01959   drv = wrFieldOnCells->addDriver(MED_DRIVER, outfile, fieldName);
01960   wrFieldOnCells->write( drv );
01961   wrFieldOnCells->removeReference();
01962 
01963   FIELD<double> *wrFieldOnNodes=new FIELD<double>(supportOnNodes,1);
01964   wrFieldOnNodes->setName(fieldName);
01965   wrFieldOnNodes->setComponentName(1,"Vx");
01966   wrFieldOnNodes->setComponentDescription(1,"comp1");
01967   wrFieldOnNodes->setMEDComponentUnit(1,"unit1");
01968   wrFieldOnNodes->setValue( nodeValues );
01969   drv = wrFieldOnNodes->addDriver(MED_DRIVER, outfile, fieldName);
01970   wrFieldOnNodes->write( drv );
01971   wrFieldOnNodes->removeReference();
01972 
01973   //  READ FIELDS BACK
01974 
01975   //  field on CELLs
01976   FIELD<double> *cellField=new FIELD<double>(supportOnCells, MED_DRIVER, outfile, fieldName, -1, -1 );
01977   CPPUNIT_ASSERT_EQUAL( MED_CELL, cellField->getSupport()->getEntity());
01978   CPPUNIT_ASSERT_DOUBLES_EQUAL( numberOfCells-1, cellField->getValueIJ( numberOfCells, 1 ),1e-20);
01979   cellField->removeReference();
01980 
01981   //  field on NODEs
01982   FIELD<double> *nodeField=new FIELD<double>(supportOnNodes, MED_DRIVER, outfile, fieldName, -1, -1 );
01983   CPPUNIT_ASSERT_EQUAL( MED_NODE, nodeField->getSupport()->getEntity());
01984   CPPUNIT_ASSERT_DOUBLES_EQUAL( -(numberOfNodes-1), nodeField->getValueIJ( numberOfNodes, 1 ),1e-20);
01985   nodeField->removeReference();
01986 
01987   supportOnCells->removeReference();
01988   supportOnNodes->removeReference();
01989 }
01990 
01991 
01992 //================================================================================
01998 //================================================================================
01999 void MEDMEMTest::testGetGaussPointsCoordinates() 
02000 {
02001   const int spaceDimension = 3;
02002   const int numberOfNodes = 56;
02003 
02004   string names[3] =
02005     {
02006       "X","Y","Z"
02007     };
02008   string units[3] =
02009     {
02010       "cm","cm","cm"
02011     };
02012 
02013   const int numberOfCellTypes = 14;
02014 
02015   double coordinates [ spaceDimension*numberOfNodes ] =
02016     {
02017       0.0,     0.0,      0.0,   //N1
02018       0.5,     0.5,      0.5,   //N2
02019       1.0,     1.0,      1.0,   //N3
02020 
02021       1.0,     1.0,      0.0,   //N4
02022       2.0,     2.5,      0.0,   //N5
02023       6.0,     1.5,      0.0,   //N6
02024       1.0,     2.0,      0.0,   //N7
02025       4.5,     2.5,      0.0,   //N8
02026       4.0,     0.5,      0.0,   //N9
02027 
02028       0.0,     4.0,      0.0,   //N10
02029       4.0,     4.0,      0.0,   //N11
02030       4.0,     0.0,      0.0,   //N12
02031       0.0,     2.0,      0.0,   //N13
02032       2.0,     4.0,      0.0,   //N14
02033       4.0,     2.0,      0.0,   //N15
02034       2.0,     0.0,      0.0,   //N16
02035 
02036       0.0,     6.0,      0.0,   //N17
02037       3.0,     3.0,      0.0,   //N18
02038       1.3,     3.0,      3.0,   //N19
02039 
02040       0.0,     3.0,      0.0,   //N20
02041       1.5,     4.5,      0.0,   //N21
02042       1.5,     1.5,      0.0,   //N22
02043       0.65,    1.5,      1.5,   //N23
02044       0.65,    4.5,      1.5,   //N24
02045       2.15,    3.0,      1.5,   //N25
02046 
02047       2.0,     2.0,      2.0,   //N26
02048       3.0,     1.0,      1.0,   //N27
02049       3.0,     3.0,      1.0,   //N28
02050       1.0,     3.0,      1.0,   //N29
02051       1.0,     1.0,      1.0,   //N30
02052 
02053       0.0,     3.0,      0.0,   //N31 
02054       2.0,     0.0,      0.0,   //N32 
02055       0.0,     0.0,      6.0,   //N33 
02056       0.0,     3.0,      6.0,   //N34 
02057       3.0,     0.0,      6.0,   //N35 
02058 
02059       0.0,     1.5,      0.0,   //N36 
02060       1.5,     1.5,      0.0,   //N37 
02061       1.5,     0.0,      0.0,   //N38 
02062       0.0,     1.5,      6.0,   //N39 
02063       1.5,     1.5,      6.0,   //N40 
02064       1.5,     0.0,      6.0,   //N41 
02065       0.0,     0.0,      3.0,   //N42
02066       0.0,     3.0,      3.0,   //N43
02067       3.0,     0.0,      3.0,   //N44
02068 
02069       0.0,     0.0,      4.0,   //N45
02070       0.0,     4.0,      4.0,   //N46
02071       4.0,     4.0,      4.0,   //N47
02072       4.0,     0.0,      4.0,   //N48
02073       0.0,     2.0,      4.0,   //N49 
02074       2.0,     4.0,      4.0,   //N50
02075       4.0,     2.0,      4.0,   //N51
02076       2.0,     0.0,      4.0,   //N52
02077       0.0,     0.0,      2.0,   //N53
02078       0.0,     4.0,      2.0,   //N54
02079       4.0,     4.0,      2.0,   //N55
02080       4.0,     0.0,      2.0,   //N56
02081     };
02082 
02083   MED_EN::medGeometryElement cellTypes[ numberOfCellTypes ] =
02084     {
02085       MED_EN::MED_SEG2,
02086       MED_EN::MED_SEG3,
02087       MED_EN::MED_TRIA3,
02088       MED_EN::MED_TRIA6,
02089       MED_EN::MED_QUAD4,
02090       MED_EN::MED_QUAD8,
02091       MED_EN::MED_TETRA4,
02092       MED_EN::MED_TETRA10,
02093       MED_EN::MED_PYRA5,
02094       MED_EN::MED_PYRA13,
02095       MED_EN::MED_PENTA6,
02096       MED_EN::MED_PENTA15,
02097       MED_EN::MED_HEXA8,
02098       MED_EN::MED_HEXA20
02099     };
02100 
02101   const int numberOfCells[numberOfCellTypes] =
02102     {
02103       1, 1,                   //1D
02104       1, 1, 1, 1,             //2D 
02105       1, 1, 1, 1, 1, 1, 1, 1  //3D
02106     };
02107 
02108   //Seg2 Connectivity
02109   int seg2C [ 2 ] =
02110     {
02111       1,3
02112     };
02113 
02114   //Seg3 Connectivity
02115   int seg3C [ 3 ] =
02116     {
02117       1,3,2
02118     };
02119 
02120   //Tria3 Connectivity
02121   int tria3C [ 3 ] =
02122     {
02123       4,5,6
02124     };
02125 
02126   //Tria6 Connectivity
02127   int tria6C [ 6 ] =
02128     {
02129       4,5,6,7,8,9
02130     };
02131 
02132   //Quad4 Connectivity
02133   int quad4C [4] =
02134     {
02135       1, 10, 11, 12
02136     };
02137 
02138   //Quad8 Connectivity
02139   int quad8C [8] =
02140     {
02141       1, 10, 11, 12, 13, 14, 15, 16
02142     };
02143 
02144   //Tetra4 Connectivity
02145   int tetra4C [4] =
02146     {
02147       1, 17, 18, 19
02148     };
02149 
02150   //Tetra10 Connectivity
02151   int tetra10C [10] =
02152     {
02153       1, 17, 18, 19, 20, 21, 22, 23, 24, 25
02154     };
02155 
02156   //Pyra13 Connectivity
02157   int pyra5C [5] =
02158     {
02159       1, 12, 11, 10, 26
02160     };
02161 
02162   //Pyra13 Connectivity
02163   int pyra13C [13] =
02164     {
02165       1, 12, 11, 10, 26, 16, 15, 14, 13, 27, 28, 29, 30
02166     };
02167 
02168   //Penta6 Connectivity
02169   int penta6C [6] =
02170     {
02171       1, 31, 32, 33, 34, 35
02172     };
02173 
02174   //Penta6 Connectivity
02175   int penta15C [15] =
02176     {
02177       1, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44
02178     };
02179 
02180   //Hexa8 Connectivity 
02181   int hexa8C[8] =
02182     {
02183       1, 10, 11, 12, 45, 46, 47, 48
02184     };
02185 
02186   //Hexa20 Connectivity 
02187   int hexa20C[20] =
02188     {
02189       1, 10, 11, 12, 45, 46, 47, 48, 13, 14, 15, 16, 49, 50, 51, 52, 53, 54, 55, 56
02190     };
02191 
02192 
02193 
02194   MEDMEM::MESHING* mesh = new MEDMEM::MESHING;
02195   mesh->setCoordinates(spaceDimension, numberOfNodes, coordinates,
02196                        "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
02197   mesh->setCoordinatesNames(names);
02198   mesh->setCoordinatesUnits(units);
02199 
02200   //connectivities
02201   mesh->setNumberOfTypes(numberOfCellTypes, MED_EN::MED_CELL);
02202   mesh->setTypes(cellTypes, MED_EN::MED_CELL);
02203   mesh->setNumberOfElements(numberOfCells, MED_EN::MED_CELL);
02204 
02205   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_SEG2, seg2C );
02206   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_SEG3, seg3C );
02207   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TRIA3, tria3C );
02208   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TRIA6, tria6C );
02209   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_QUAD4, quad4C );
02210   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_QUAD8, quad8C );
02211   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TETRA4, tetra4C );
02212   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TETRA10, tetra10C );
02213   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PYRA5, pyra5C );
02214   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PYRA13, pyra13C );
02215   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PENTA6, penta6C );
02216   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PENTA15, penta15C );
02217   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_HEXA8, hexa8C );
02218   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_HEXA20, hexa20C );
02219 
02220 
02221   //Support definition 
02222   const SUPPORT* sup = mesh->getSupportOnAll( MED_EN::MED_CELL );
02223 
02224   //Create test field
02225   FIELD<int>* aField = new FIELD<int>(sup,1);
02226 
02227   //Gauss Localization definition:
02228   double v[1] =
02229     {
02230       0.3
02231     };
02232   double v_2[2] =
02233     {
02234       0.3, 0.3
02235     };
02236   double v_3[3] =
02237     {
02238       0.3, 0.3, 0.3
02239     };
02240   double v_4[4] =
02241     {
02242       0.3, 0.3, 0.3, 0.3
02243     };
02244 
02245   //              --------- Seg2 localization ---------
02246   //              Nb of the gauss points = 1 
02247   double seg2CooRef[2] =
02248     {
02249       -1.0 , 1.0
02250     };
02251   double seg2CooGauss[1] =
02252     {
02253       0.2
02254     };
02255 
02256   GAUSS_LOCALIZATION<>* aseg2L = new GAUSS_LOCALIZATION<>("Seg2 Gauss localization",
02257                                                           MED_EN::MED_SEG2,
02258                                                           1,
02259                                                           seg2CooRef,
02260                                                           seg2CooGauss,
02261                                                           v);
02262 
02263   //              --------- Seg3 localization ---------
02264   //              Nb of the gauss points = 1
02265   double seg3CooRef[3] =
02266     {
02267       -1.0, 1.0, 0.0
02268     };
02269   double seg3CooGauss[1] =
02270     {
02271       0.2
02272     };
02273 
02274   GAUSS_LOCALIZATION<>* aseg3L = new GAUSS_LOCALIZATION<>("Seg3 Gauss localization",
02275                                                           MED_EN::MED_SEG3,
02276                                                           1,
02277                                                           seg3CooRef,
02278                                                           seg3CooGauss,
02279                                                           v);
02280   //              --------- Tria3 localization ---------
02281   //              Nb of the gauss points = 2
02282   double tria3CooRef[6] =
02283     {
02284       0.0, 0.0, 1.0 , 0.0, 0.0, 1.0
02285     };
02286 
02287   double tria3CooGauss[4] =
02288     {
02289       0.1, 0.8, 0.2, 0.7
02290     };
02291 
02292   GAUSS_LOCALIZATION<>* atria3L = new GAUSS_LOCALIZATION<>("Tria3 Gauss localization",
02293                                                            MED_EN::MED_TRIA3,
02294                                                            2,
02295                                                            tria3CooRef,
02296                                                            tria3CooGauss,
02297                                                            v_2);
02298 
02299   //              --------- Tria6 localization ---------
02300   //              Nb of the gauss points = 3
02301   double tria6CooRef[12] =
02302     {
02303       0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5
02304     };
02305 
02306   double tria6CooGauss[6] =
02307     {
02308       0.3, 0.2, 0.2, 0.1, 0.2, 0.4
02309     };
02310 
02311   GAUSS_LOCALIZATION<>* atria6L = new GAUSS_LOCALIZATION<>("Tria6 Gauss localization",
02312                                                            MED_EN::MED_TRIA6,
02313                                                            3,
02314                                                            tria6CooRef,
02315                                                            tria6CooGauss,
02316                                                            v_3);
02317   //              --------- Quad4 localization ---------
02318   //              Nb of the gauss points = 4
02319   double quad4CooRef[8] =
02320     {
02321       -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0
02322     };
02323 
02324   double quad4CooGauss[8] =
02325     { //Size is type/10*NbGauss = (204/100)*4 = 8
02326       0.3, 0.2, 0.2, 0.1, 0.2, 0.4, 0.15, 0.27
02327     };
02328 
02329   GAUSS_LOCALIZATION<>* aquad4L = new GAUSS_LOCALIZATION<>("Quad8 Gauss localization",
02330                                                            MED_EN::MED_QUAD4,
02331                                                            4,
02332                                                            quad4CooRef,
02333                                                            quad4CooGauss,
02334                                                            v_4);
02335   //              --------- Quad8 localization ---------
02336   //              Nb of the gauss points = 4
02337   double quad8CooRef[16] =
02338     {
02339       -1.0, -1.0,
02340       1.0, -1.0,
02341       1.0,  1.0,
02342       -1.0,  1.0,
02343       0.0, -1.0,
02344       1.0,  0.0,
02345       0.0,  1.0,
02346       -1.0,  0.0
02347     };
02348 
02349   double quad8CooGauss[8] =
02350     {
02351       0.34, 0.16, 0.21, 0.3, 0.23, 0.4, 0.14, 0.37
02352     };
02353 
02354   GAUSS_LOCALIZATION<>* aquad8L = new GAUSS_LOCALIZATION<>("Quad8 Gauss localization",
02355                                                            MED_EN::MED_QUAD8,
02356                                                            4,
02357                                                            quad8CooRef,
02358                                                            quad8CooGauss,
02359                                                            v_4);
02360 
02361   //              --------- Tetra4 localization
02362   //              Nb of the gauss points = 1
02363   double tetra4CooRef[12] =
02364     {
02365       0.0,  1.0,  0.0,
02366       0.0,  0.0,  1.0, 
02367       0.0,  0.0,  0.0,
02368       1.0,  0.0,  0.0
02369     };
02370 
02371   double tetra4CooGauss[3] =
02372     {
02373       0.34, 0.16, 0.21
02374     };
02375 
02376   GAUSS_LOCALIZATION<>* atetra4L = new GAUSS_LOCALIZATION<>("Tetra4 Gauss localization",
02377                                                             MED_EN::MED_TETRA4,
02378                                                             1,
02379                                                             tetra4CooRef,
02380                                                             tetra4CooGauss,
02381                                                             v);
02382 
02383   //              --------- Tetra10 localization
02384   //              Nb of the gauss points = 1
02385   double tetra10CooRef[30] =
02386     {
02387       0.0, 1.0, 0.0,
02388       0.0, 0.0, 0.0,
02389       0.0, 0.0, 1.0,
02390       1.0, 0.0, 0.0,
02391       0.0, 0.5, 0.0,
02392       0.0, 0.0, 0.5,
02393       0.0, 0.5, 0.5,
02394       0.5, 0.5, 0.0,
02395       0.5, 0.0, 0.0,
02396       0.5, 0.0, 0.5,
02397     };
02398 
02399   double tetra10CooGauss[3] =
02400     {
02401       0.2, 0.3, 0.1
02402     };
02403 
02404   GAUSS_LOCALIZATION<>* atetra10L = new GAUSS_LOCALIZATION<>("Tetra10 Gauss localization",
02405                                                              MED_EN::MED_TETRA10,
02406                                                              1,
02407                                                              tetra10CooRef,
02408                                                              tetra10CooGauss,
02409                                                              v);
02410   //              --------- Pyra5 localization
02411   //              Nb of the gauss points = 1
02412   double pyra5CooRef[15] =
02413     {
02414       1.0,  0.0,  0.0,
02415       0.0,  1.0,  0.0,
02416       -1.0,  0.0,  0.0,
02417       0.0, -1.0,  0.0,
02418       0.0,  0.0,  1.0
02419     };
02420 
02421   double pyra5CooGauss[3] =
02422     {
02423       0.2, 0.3, 0.1
02424     };
02425 
02426   GAUSS_LOCALIZATION<>* apyra5L = new GAUSS_LOCALIZATION<>("Pyra5 Gauss localization",
02427                                                            MED_EN::MED_PYRA5,
02428                                                            1,
02429                                                            pyra5CooRef,
02430                                                            pyra5CooGauss,
02431                                                            v);
02432 
02433   //              --------- Pyra13 localization
02434   //              Nb of the gauss points = 1
02435   double pyra13CooRef[39] =
02436     {
02437       1.0,  0.0, 0.0,
02438       0.0,  1.0, 0.0,
02439       -1.0,  0.0, 0.0,
02440       0.0, -1.0, 0.0,
02441       0.0,  0.0, 1.0,
02442       0.5,  0.5, 0.0,
02443       -0.5,  0.5, 0.0,
02444       -0.5, -0.5, 0.0,
02445       0.5, -0.5, 0.0,
02446       0.5,  0.0, 0.5,
02447       0.0,  0.5, 0.5,
02448       -0.5,  0.0, 0.5,
02449       0.0, -0.5, 0.5
02450     };
02451 
02452   double pyra13CooGauss[3] =
02453     {
02454       0.1, 0.2, 0.7
02455     };
02456 
02457   GAUSS_LOCALIZATION<>* apyra13L = new GAUSS_LOCALIZATION<>("Pyra13 Gauss localization",
02458                                                             MED_EN::MED_PYRA13,
02459                                                             1,
02460                                                             pyra13CooRef,
02461                                                             pyra13CooGauss,
02462                                                             v);
02463   //              --------- Penta6 localization
02464   //              Nb of the gauss points = 1
02465   double penta6CooRef[18] =
02466     {
02467       -1.0,   1.0,   0.0,
02468       -1.0,  -0.0,   1.0,
02469       -1.0,   0.0,   0.0,
02470       1.0,   1.0,   0.0,
02471       1.0,   0.0,   1.0,
02472       1.0,   0.0,   0.0
02473     };
02474 
02475   double penta6CooGauss[3] =
02476     {
02477       0.2, 0.3, 0.1
02478     };
02479 
02480   GAUSS_LOCALIZATION<>* apenta6L = new GAUSS_LOCALIZATION<>("Penta6 Gauss localization",
02481                                                             MED_EN::MED_PENTA6,
02482                                                             1,
02483                                                             penta6CooRef,
02484                                                             penta6CooGauss,
02485                                                             v);
02486 
02487   //              --------- Penta15 localization
02488   //              Nb of the gauss points = 1
02489   double penta15CooRef[45] =
02490     {
02491       -1.0,  1.0,   0.0,
02492       -1.0,  0.0,   1.0,
02493       -1.0,  0.0,   0.0,
02494       1.0,  1.0,   0.0,
02495       1.0,  0.0,   1.0,
02496       1.0,  0.0,   0.0,
02497       -1.0,  0.5,   0.5,
02498       -1.0,  0.0,   0.5,
02499       -1.0,  0.5,   0.0,
02500       0.0,  1.0,   0.0,
02501       0.0,  0.0,   1.0,
02502       0.0,  0.0,   0.0,
02503       1.0,  0.5,   0.5,
02504       1.0,  0.0,   0.5,
02505       1.0,  0.5,   0.0
02506     };
02507 
02508   double penta15CooGauss[3] =
02509     {
02510       0.2, 0.3, 0.15
02511     };
02512 
02513   GAUSS_LOCALIZATION<>* apenta15L = new GAUSS_LOCALIZATION<>("Penta15 Gauss localization",
02514                                                              MED_EN::MED_PENTA15,
02515                                                              1,
02516                                                              penta15CooRef,
02517                                                              penta15CooGauss,
02518                                                              v);
02519 
02520   //              --------- Hexa8 localization
02521   //              Nb of the gauss points = 1
02522   double hexa8CooRef [24] =
02523     {
02524       -1.0,  -1.0,  -1.0,
02525       1.0,  -1.0,  -1.0,
02526       1.0,   1.0,  -1.0,
02527       -1.0,   1.0,  -1.0,
02528       -1.0,  -1.0,   1.0,
02529       1.0,  -1.0,   1.0,
02530       1.0,   1.0,   1.0,
02531       -1.0,   1.0,   1.0
02532     };
02533 
02534   double hexa8CooGauss[3] =
02535     {
02536       0.2, 0.3, 0.15
02537     };
02538 
02539   GAUSS_LOCALIZATION<>* ahexa8L = new GAUSS_LOCALIZATION<>("Hexa8 Gauss localization",
02540                                                            MED_EN::MED_HEXA8,
02541                                                            1,
02542                                                            hexa8CooRef,
02543                                                            hexa8CooGauss,
02544                                                            v);
02545 
02546   //              --------- Hexa20 localization
02547   //              Nb of the gauss points = 1
02548   double hexa20CooRef[60] =
02549     {
02550       -1.0,  -1.0,  -1.0,
02551       1.0,  -1.0,  -1.0,
02552       1.0,   1.0,  -1.0,
02553       -1.0,   1.0,  -1.0,
02554       -1.0,  -1.0,   1.0,
02555       1.0,  -1.0,   1.0,
02556       1.0,   1.0,   1.0,
02557       -1.0,   1.0,   1.0,
02558       0.0,  -1.0,  -1.0,
02559       1.0,   0.0,  -1.0,
02560       0.0,   1.0,  -1.0,
02561       -1.0,   0.0,  -1.0,
02562       -1.0,  -1.0,   0.0,
02563       1.0,  -1.0,   0.0,
02564       1.0,   1.0,   0.0,
02565       -1.0,   1.0,   0.0,
02566       0.0,  -1.0,   1.0,
02567       1.0,   0.0,   1.0,
02568       0.0,   1.0,   1.0,
02569       -1.0,   0.0,   1.0
02570     };
02571 
02572   double hexa20CooGauss[3] =
02573     {
02574       0.11, 0.3, 0.55
02575     };
02576 
02577   GAUSS_LOCALIZATION<>* ahexa20L = new GAUSS_LOCALIZATION<>("Hexa20 Gauss localization",
02578                                                             MED_EN::MED_HEXA20,
02579                                                             1,
02580                                                             hexa20CooRef,
02581                                                             hexa20CooGauss,
02582                                                             v);
02583 
02584 
02585 
02586   aField->setGaussLocalization(MED_EN::MED_SEG2, aseg2L);
02587   aField->setGaussLocalization(MED_EN::MED_SEG3, aseg3L);
02588   aField->setGaussLocalization(MED_EN::MED_TRIA3, atria3L);
02589   aField->setGaussLocalization(MED_EN::MED_TRIA6, atria6L);
02590   aField->setGaussLocalization(MED_EN::MED_QUAD4, aquad4L);
02591   aField->setGaussLocalization(MED_EN::MED_QUAD8, aquad8L);
02592   aField->setGaussLocalization(MED_EN::MED_TETRA4, atetra4L);
02593   aField->setGaussLocalization(MED_EN::MED_TETRA10, atetra10L);
02594   aField->setGaussLocalization(MED_EN::MED_PYRA5, apyra5L);
02595   aField->setGaussLocalization(MED_EN::MED_PYRA13, apyra13L);
02596   aField->setGaussLocalization(MED_EN::MED_PENTA6, apenta6L);
02597   aField->setGaussLocalization(MED_EN::MED_PENTA15, apenta15L);
02598   aField->setGaussLocalization(MED_EN::MED_HEXA8, ahexa8L);
02599   aField->setGaussLocalization(MED_EN::MED_HEXA20, ahexa20L);
02600 
02601   FIELD<double>* aGaussCoords = aField->getGaussPointsCoordinates();
02602 
02603 
02604   //Coordinates to check result
02605   double seg2Coord[3] =
02606     {
02607       0.6,0.6,0.6
02608     };
02609   double seg3Coord[3] =
02610     {
02611       0.6,0.6,0.6
02612     };
02613 
02614   double tria3Coord[3*2]  =
02615     {
02616       5.1, 1.55, 0.0,   //First GP Coordinates 
02617       4.7, 1.65, 0.0    //Second GP Coordinates
02618     };
02619 
02620   double tria6Coord[3*3] =
02621     {
02622       2.32, 1.52, 0.0, //First GP Coordinates  
02623       1.6 , 1.32, 0.0, //Second GP Coordinates  
02624       3.52, 1.26, 0.0  //Third GP Coordinates
02625     };
02626 
02627   double quad4Coord[4*3] =
02628     {
02629       2.6, 1.6,  0.0,
02630       2.4, 1.8,  0.0,
02631       2.4, 1.2,  0.0,
02632       2.3, 1.46, 0.0
02633     };
02634 
02635   double quad8Coord[4*3] =
02636     {
02637       2.32, 2.68,  0.0,
02638       2.6,  2.42,  0.0,
02639       2.8,  2.46,  0.0,
02640       2.74, 2.28,  0.0
02641     };
02642 
02643   double tetra4Coord[3] =
02644     {
02645       1.312, 3.15, 1.02
02646     };
02647 
02648   double tetra10Coord[3] =
02649     {
02650       0.56, 3.3, 0.6
02651     };
02652 
02653   double pyra5Coord [3]=
02654     {
02655       2.18, 1.1, 0.2
02656     };
02657 
02658   double pyra13Coord [3] =
02659     {
02660       1.18, 1.54, 0.98
02661     };
02662 
02663   double penta6Coord [3] =
02664     {
02665       1.56, 0.3, 3.6
02666     };
02667 
02668   double penta15Coord [3] =
02669     {
02670       1.613, 0.801, 4.374
02671     };
02672 
02673   double hexa8Coord [3] =
02674     {
02675       2.6, 2.4, 2.3
02676     };
02677 
02678   double hexa20Coord [3] =
02679     {
02680       2.31232, 2.3934, 1.55326
02681     };
02682 
02683   //Check result of the calculation
02684   int ElemId = 1;
02685   double EPS = 0.000001;
02686   int idx = 0;
02687 
02688   //Seg2:
02689   for(int j = 1; j<=3;j++)
02690     {
02691       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02692                                    seg2Coord[j-1], EPS);
02693     }
02694 
02695   //Seg3:
02696   ElemId++;
02697   for(int j = 1; j<=3;j++)
02698     {
02699       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02700                                    seg3Coord[j-1], EPS);
02701     }
02702 
02703   //Tria3
02704   ElemId++;
02705   for(int k = 1; k <=2 ; k++)
02706     {
02707       for(int j = 1; j<=3;j++)
02708         {
02709           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
02710                                        tria3Coord[idx++], EPS);
02711         }
02712     }
02713 
02714   //Tria6
02715   ElemId++;
02716   idx=0;
02717   for(int k = 1; k <=3 ; k++)
02718     {
02719       for(int j = 1; j<=3;j++)
02720         {
02721           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
02722                                        tria6Coord[idx++], EPS);
02723         }
02724     }
02725 
02726   //Quad4
02727   ElemId++;
02728   idx=0;
02729   for(int k = 1; k <=4 ; k++)
02730     {
02731       for(int j = 1; j<=3;j++)
02732         {
02733           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
02734                                        quad4Coord[idx++], EPS);
02735         }
02736     }
02737 
02738   //Quad8
02739   ElemId++;
02740   idx=0;
02741   for(int k = 1; k <=4 ; k++)
02742     {
02743       for(int j = 1; j<=3;j++)
02744         {
02745           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
02746                                        quad8Coord[idx++], EPS);
02747         }
02748     }
02749 
02750   //Tetra4
02751   ElemId++;
02752   for(int j = 1; j<=3;j++)
02753     {
02754       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02755                                    tetra4Coord[j-1], EPS);
02756     }
02757 
02758   //Tetra10
02759   ElemId++;
02760   for(int j = 1; j<=3;j++)
02761     {
02762       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02763                                    tetra10Coord[j-1], EPS);
02764     }
02765 
02766   //Pyra5
02767   ElemId++;
02768   for(int j = 1; j<=3;j++)
02769     {
02770       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02771                                    pyra5Coord[j-1], EPS);
02772     }
02773 
02774   //Penta15
02775   ElemId++;
02776   for(int j = 1; j<=3;j++)
02777     {
02778       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02779                                    pyra13Coord[j-1], EPS);
02780     }
02781 
02782   //Penta6
02783   ElemId++;
02784   for(int j = 1; j<=3;j++)
02785     {
02786       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02787                                    penta6Coord[j-1], EPS);
02788     }
02789 
02790   //Penta15
02791   ElemId++;
02792   for(int j = 1; j<=3;j++)
02793     {
02794       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02795                                    penta15Coord[j-1], EPS);
02796     }
02797 
02798   //Hexa8
02799   ElemId++;
02800   for(int j = 1; j<=3;j++)
02801     {
02802       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02803                                    hexa8Coord[j-1], EPS);
02804     }
02805 
02806   //Hexa20
02807   ElemId++;
02808   for(int j = 1; j<=3;j++)
02809     {
02810       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
02811                                    hexa20Coord[j-1], 0.0001);
02812     }
02813 
02814   aGaussCoords->removeReference();
02815   aField->removeReference();
02816   mesh->removeReference();
02817 }