Back to index

salome-med  6.5.0
MEDMEMTest_Field_fault.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_Mesh.hxx"
00025 #include "MEDMEM_Group.hxx"
00026 #include "MEDMEM_Support.hxx"
00027 #include "MEDMEM_VtkMeshDriver.hxx"
00028 #include "MEDMEM_MedMeshDriver.hxx"
00029 
00030 #include <cppunit/TestAssert.h>
00031 
00032 #include <sstream>
00033 #include <cmath>
00034 
00035 // use this define to enable lines, execution of which leads to Segmentation Fault
00036 //#define ENABLE_FAULTS
00037 
00038 // use this define to enable CPPUNIT asserts and fails, showing bugs
00039 //#define ENABLE_FORCED_FAILURES
00040 
00041 using namespace std;
00042 using namespace MEDMEM;
00043 
00044 // #14,15: MEDMEMTest_Field.cxx
00045 // Check methods from MEDMEM_Field.hxx, MEDMEM_FieldConvert.hxx
00046 
00231  static void compareField_(const FIELD_ * theField_1, const FIELD_ * theField_2, bool isFIELD, bool isValue)
00232 {
00233   // name, description, support
00234   CPPUNIT_ASSERT_EQUAL(theField_1->getName(), theField_2->getName());
00235   CPPUNIT_ASSERT_EQUAL(theField_1->getDescription(), theField_2->getDescription());
00236   CPPUNIT_ASSERT_EQUAL(theField_1->getSupport(), theField_2->getSupport());
00237 
00238   // components information
00239   int aNbComps = theField_1->getNumberOfComponents();
00240   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_2->getNumberOfComponents());
00241 
00242   for (int i = 1; i <= aNbComps; i++) 
00243     {
00244       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentName(i), theField_2->getComponentName(i));
00245       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentDescription(i), theField_2->getComponentDescription(i));
00246       CPPUNIT_ASSERT_EQUAL(theField_1->getMEDComponentUnit(i), theField_2->getMEDComponentUnit(i));
00247     }
00248 
00249   // iteration information
00250   CPPUNIT_ASSERT_EQUAL(theField_1->getIterationNumber(), theField_2->getIterationNumber());
00251   CPPUNIT_ASSERT_EQUAL(theField_1->getOrderNumber(), theField_2->getOrderNumber());
00252   CPPUNIT_ASSERT_DOUBLES_EQUAL(theField_1->getTime(), theField_2->getTime(), 0.0000001);
00253 
00254   // Value
00255   int nbOfValues = theField_1->getNumberOfValues();
00256   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_2->getNumberOfValues());
00257 
00258   if (isFIELD) 
00259     {
00260       // Value type and Interlacing type
00261       CPPUNIT_ASSERT_EQUAL(theField_1->getValueType(), theField_2->getValueType());
00262       CPPUNIT_ASSERT_EQUAL(theField_1->getInterlacingType(), theField_2->getInterlacingType());
00263 
00264       // Gauss Presence
00265       if (isValue) 
00266         {
00267           CPPUNIT_ASSERT_EQUAL(theField_1->getGaussPresence(), theField_2->getGaussPresence());
00268         }
00269       else
00270         {
00271           CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
00272           CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
00273         }
00274     }
00275   else
00276     {
00277       CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
00278       CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
00279     }
00280 }
00281 
00282 static void checkField_(FIELD_ * theField_, const SUPPORT * theSupport,
00283                         MED_EN::med_type_champ theValueType,
00284                         MED_EN::medModeSwitch theInterlace)
00285 {
00286   // name
00287   const string aFieldName = "a_name_of_a_field";
00288   theField_->setName(aFieldName);
00289   CPPUNIT_ASSERT_EQUAL(aFieldName, theField_->getName());
00290 
00291   // description
00292   const string aFieldDescr = "a_description_of_a_field";
00293   theField_->setDescription(aFieldDescr);
00294   CPPUNIT_ASSERT_EQUAL(aFieldDescr, theField_->getDescription());
00295 
00296   // support
00297   theField_->setSupport(theSupport);
00298   CPPUNIT_ASSERT(theField_->getSupport() == theSupport);
00299 
00300   // components information
00301   int aNbComps = 3;
00302 
00303   string aCompsNames[3] = 
00304     {
00305       "Vx", "Vy", "Vz"
00306     };
00307   string aCompsDescs[3] = 
00308     {
00309       "vitesse selon x", "vitesse selon y", "vitesse selon z"
00310     };
00311   string aCompsUnits[3] = 
00312     {
00313       "m.s-1", "m.s-1", "m.s-1"
00314     };
00315 
00316   theField_->setNumberOfComponents(aNbComps);
00317   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_->getNumberOfComponents());
00318 
00319   theField_->setComponentsNames(aCompsNames);
00320 
00321   //#ifdef ENABLE_FAULTS
00322   try
00323     {
00324       theField_->setNumberOfComponents(7);
00325       // Segmentation fault here because array of components names is not resized
00326       for (int i = 1; i <= 7; i++) 
00327         {
00328           theField_->setComponentName(i, "AnyComponent");
00329         }
00330     }
00331   catch (MEDEXCEPTION& ex) 
00332     {
00333       // Ok, it is good to have MEDEXCEPTION here
00334     }
00335   catch (...) 
00336     {
00337       CPPUNIT_FAIL("Unknown exception cought");
00338     }
00339   // restore components names
00340   theField_->setNumberOfComponents(aNbComps);
00341   theField_->setComponentsNames(aCompsNames);
00342   //#endif
00343   //#ifdef ENABLE_FORCED_FAILURES
00344   //  CPPUNIT_FAIL("FIELD_::_componentsNames bad management");
00345   //#endif
00346 
00347   theField_->setComponentsDescriptions(aCompsDescs);
00348   theField_->setMEDComponentsUnits(aCompsUnits);
00349 
00350   const string * aCompsNamesBack = theField_->getComponentsNames();
00351   const string * aCompsDescsBack = theField_->getComponentsDescriptions();
00352   const string * aCompsUnitsBack = theField_->getMEDComponentsUnits();
00353   for (int i = 1; i <= aNbComps; i++) 
00354     {
00355       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], theField_->getComponentName(i));
00356       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], aCompsNames[i-1]);
00357 
00358       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], theField_->getComponentDescription(i));
00359       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], aCompsDescs[i-1]);
00360 
00361       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], theField_->getMEDComponentUnit(i));
00362       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], aCompsUnits[i-1]);
00363     }
00364 
00365   const string aCompName2 ("Name of second component");
00366   const string aCompDesc2 ("Description of second component");
00367   const string aCompUnit2 ("Unit of second MED component");
00368 
00369   theField_->setComponentName(2, aCompName2);
00370   theField_->setComponentDescription(2, aCompDesc2);
00371   theField_->setMEDComponentUnit(2, aCompUnit2);
00372 
00373   const string * aCompsNamesBack2 = theField_->getComponentsNames();
00374   const string * aCompsDescsBack2 = theField_->getComponentsDescriptions();
00375   const string * aCompsUnitsBack2 = theField_->getMEDComponentsUnits();
00376 
00377   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], theField_->getComponentName(2));
00378   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], aCompName2);
00379 
00380   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], theField_->getComponentDescription(2));
00381   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], aCompDesc2);
00382 
00383   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], theField_->getMEDComponentUnit(2));
00384   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], aCompUnit2);
00385 
00386   //#ifdef ENABLE_FAULTS
00387   // (BUG) No index checking
00388   CPPUNIT_ASSERT_THROW(theField_->setComponentName(0, "str"), MEDEXCEPTION);
00389   CPPUNIT_ASSERT_THROW(theField_->setComponentName(aNbComps + 1, "str"), MEDEXCEPTION);
00390   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(0, "str"), MEDEXCEPTION);
00391   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(aNbComps + 1, "str"), MEDEXCEPTION);
00392   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(0, "str"), MEDEXCEPTION);
00393   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(aNbComps + 1, "str"), MEDEXCEPTION);
00394   //#endif
00395   //#ifdef ENABLE_FORCED_FAILURES
00396   //  CPPUNIT_FAIL("FIELD::setComponentXXX() does not check component index");
00397   //#endif
00398 
00399   // iteration information
00400   int anIterNumber = 10; // set value to MED_NOPDT if undefined (default)
00401   theField_->setIterationNumber(anIterNumber);
00402   CPPUNIT_ASSERT_EQUAL(anIterNumber, theField_->getIterationNumber());
00403 
00404   int anOrderNumber = 1; // set value to MED_NONOR if undefined (default)
00405   theField_->setOrderNumber(anOrderNumber);
00406   CPPUNIT_ASSERT_EQUAL(anOrderNumber, theField_->getOrderNumber());
00407 
00408   double aTime = 3.435678; // in second
00409   theField_->setTime(aTime);
00410   CPPUNIT_ASSERT_DOUBLES_EQUAL(aTime, theField_->getTime(), 0.0000001);
00411 
00412   // Value
00413   int nbOfValues = 10;
00414   // dangerous method, because it does not reallocate values array
00415   theField_->setNumberOfValues(nbOfValues);
00416   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_->getNumberOfValues());
00417 
00418   // Value type and Interlacing type
00419   CPPUNIT_ASSERT_EQUAL(theValueType, theField_->getValueType());
00420   CPPUNIT_ASSERT_EQUAL(theInterlace, theField_->getInterlacingType());
00421 }
00422 
00423 template<class T, class INTERLACING_TAG>
00424 void compareField(const FIELD<T, INTERLACING_TAG> * theField_1,
00425                   const FIELD<T, INTERLACING_TAG> * theField_2, bool isValue)
00426 {
00427   // compare FIELD_ part
00428   compareField_(theField_1, theField_2, /*isFIELD = */true, isValue);
00429 
00430   // compare FIELD part
00431   // TO DO
00432 }
00433 
00434 template<class T, class INTERLACING_TAG>
00435 void checkField (FIELD<T, INTERLACING_TAG> * theField, const SUPPORT * theSupport)
00436 {
00437   // check FIELD_ part
00438   MED_EN::med_type_champ aValueType = SET_VALUE_TYPE<T>::_valueType;
00439   MED_EN::medModeSwitch  anInterlace = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
00440   checkField_(theField, theSupport, aValueType, anInterlace);
00441 
00442   // check FIELD part
00443 
00444   // filling by support charackteristics (NOT IMPLEMENTED METHODS!!!):
00445   // value type must be MED_REEL64 (i.e. a FIELD<double>) for these methods,
00446   // nb. of components must be equal 1 (for Volume, Area, Length) or
00447   // space dimension (for Normal, Barycenter, )
00448   {
00449     const GMESH* aMesh = theSupport->getMesh();
00450     int spaceDim = 3;
00451     if (aMesh) spaceDim = aMesh->getSpaceDimension();
00452     theField->deallocValue();
00453     theField->allocValue(/*NumberOfComponents = */spaceDim + 1);
00454 
00455     //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
00456     // getVolume() etc. does nothing
00457     //     CPPUNIT_ASSERT_THROW(theField->getVolume(), MEDEXCEPTION);
00458     //     CPPUNIT_ASSERT_THROW(theField->getArea(), MEDEXCEPTION);
00459     //     CPPUNIT_ASSERT_THROW(theField->getLength(), MEDEXCEPTION);
00460     //     if (aMesh) {
00461     //       CPPUNIT_ASSERT_THROW(theField->getNormal(), MEDEXCEPTION);
00462     //       CPPUNIT_ASSERT_THROW(theField->getBarycenter(), MEDEXCEPTION);
00463     //     }
00464 
00465     theField->deallocValue();
00466     theField->allocValue(/*NumberOfComponents = */1);
00467     //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
00468     // getVolume() etc. does nothing
00469     //     if (aValueType == MED_EN::MED_REEL64) {
00470     //       CPPUNIT_ASSERT_NO_THROW(theField->getVolume());
00471     //       CPPUNIT_ASSERT_NO_THROW(theField->getArea());
00472     //       CPPUNIT_ASSERT_NO_THROW(theField->getLength());
00473     //     }
00474     //     else {
00475     //       CPPUNIT_ASSERT_THROW(theField->getVolume(), MEDEXCEPTION);
00476     //       CPPUNIT_ASSERT_THROW(theField->getArea(), MEDEXCEPTION);
00477     //       CPPUNIT_ASSERT_THROW(theField->getLength(), MEDEXCEPTION);
00478     //     }
00479 
00480     if (aMesh) 
00481       {
00482         theField->deallocValue();
00483         theField->allocValue(/*NumberOfComponents = */spaceDim);
00484         //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
00485         // getVolume() etc. does nothing
00486         //       if (aValueType == MED_EN::MED_REEL64) {
00487         //         CPPUNIT_ASSERT_NO_THROW(theField->getNormal());
00488         //         CPPUNIT_ASSERT_NO_THROW(theField->getBarycenter());
00489         //       }
00490         //       else {
00491         //         CPPUNIT_ASSERT_THROW(theField->getNormal(), MEDEXCEPTION);
00492         //         CPPUNIT_ASSERT_THROW(theField->getBarycenter(), MEDEXCEPTION);
00493         //       }
00494       }
00495   }
00496 
00497   // values
00498   theField->deallocValue();
00499   theField->allocValue(/*NumberOfComponents = */2);
00500   int nbElemSupport = theSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
00501   CPPUNIT_ASSERT_EQUAL(nbElemSupport, theField->getNumberOfValues());
00502 
00503   //#ifdef ENABLE_FAULTS
00504   // (BUG) FIELD::deallocValue() does not nullify _value pointer,
00505   // that is why there can be failures in other methods
00506   // (even if simply call deallocValue() two times)
00507   theField->deallocValue();
00508   theField->getGaussPresence();
00509   //#endif
00510   //#ifdef ENABLE_FORCED_FAILURES
00511   //  CPPUNIT_FAIL("FIELD::deallocValue() does not nullify _value pointer");
00512   //#endif
00513 
00514   // copy constructor
00515   FIELD<T, INTERLACING_TAG> aField_copy1 (*theField);
00516   //compareField(theField, &aField_copy1, /*isValue = */false);
00517   compareField(theField, &aField_copy1, /*isValue = */true);
00518 
00519   // operator=
00520   //#ifdef ENABLE_FAULTS
00521   // (BUG) This fails (Segmentation fault) if not set:
00522   // _componentsNames or _componentsDescriptions, or _componentsUnits, or _MEDComponentsUnits
00523   FIELD<T, INTERLACING_TAG> aField_copy2;
00524   aField_copy2 = *theField;
00525   //compareField(theField, &aField_copy2, /*isValue = */false);
00526   compareField(theField, &aField_copy2, /*isValue = */true);
00527   //#endif
00528   //#ifdef ENABLE_FORCED_FAILURES
00529   //  CPPUNIT_FAIL("FIELD_::operator=() fails if _componentsUnits is not set");
00530   //#endif
00531 }
00532 
00533 template<class T>
00534 FIELD<T> * createFieldOnGroup(MESH* theMesh, const GROUP* theGroup,
00535                               const string theName, const string theDescr)
00536 {
00537   FIELD<T> * aFieldOnGroup = new FIELD<T> (theGroup, /*NumberOfComponents = */2);
00538 
00539   aFieldOnGroup->setName(theName);
00540   aFieldOnGroup->setDescription(theDescr);
00541 
00542   string aCompsNames[2] = 
00543     {
00544       "Pos", "Neg"
00545     };
00546   string aCompsDescs[2] = 
00547     {
00548       "+", "-"
00549     };
00550   string aCompsUnits[2] = 
00551     {
00552       "unit1", "unit2"
00553     };
00554 
00555   aFieldOnGroup->setComponentsNames(aCompsNames);
00556   aFieldOnGroup->setComponentsDescriptions(aCompsDescs);
00557   aFieldOnGroup->setMEDComponentsUnits(aCompsUnits);
00558 
00559   return aFieldOnGroup;
00560 }
00561 
00562 double plus13 (double val);
00563 double plus13 (double val)
00564 {
00565   return val + 13;
00566 }
00567 
00568 // function to calculate field values from coordinates of an element
00569 // typedef void (*myFuncType)(const double * temp, T* output);
00570 // size of temp array = space dim = 3
00571 // size of output array = nb. comps = 2
00572 static void proj2d (const double * temp, double* output)
00573 {
00574   // dimetric projection with coefficients:
00575   // 1.0 along Oy and Oz, 0.5 along Ox
00576   //
00577   //    ^ z (y_)
00578   //    |
00579   //    |
00580   //    .----> y (x_)
00581   //   /
00582   //  L x
00583   //
00584   // x_ = y - x * sqrt(2.) / 4.
00585   // y_ = z - x * sqrt(2.) / 4.
00586 
00587   double dx = temp[0] * std::sqrt(2.) / 4.;
00588   output[0] = temp[1] - dx;
00589   output[1] = temp[2] - dx;
00590 }
00591 
00592 static void testDrivers()
00593 {
00594   string filename_rd                  = getResourceFile("pointe.med");
00595   string filename_wr                  = makeTmpFile("myMedFieldfile.med", filename_rd);
00596   string filename_support_wr          = makeTmpFile("myMedSupportFiledfile.med");
00597   string filename22_rd                = getResourceFile("pointe.med");
00598   string filenamevtk_wr               = makeTmpFile("myMedFieldfile22.vtk");
00599 
00600   string fieldname_celldouble_rd      = "fieldcelldoublescalar";
00601   string fieldname_celldouble_wr      = fieldname_celldouble_rd + "_cpy";
00602   string fieldname_nodeint_rd         = "fieldnodeint";
00603   string fieldname_nodeint_wr         = fieldname_nodeint_rd + "_cpy";
00604   string fieldname_nodeint_wr1        = fieldname_nodeint_rd + "_cpy1";
00605   string meshname                     = "maa1";
00606 
00607   // To remove tmp files from disk
00608   MEDMEMTest_TmpFilesRemover aRemover;
00609   aRemover.Register(filename_wr);
00610   aRemover.Register(filenamevtk_wr);
00611   aRemover.Register(filename_support_wr);
00612 
00613   FIELD<int> *aInvalidField=new FIELD<int>();
00614   //must throw becase only VTK_DRIVER or MED_DRIVER may be specified as driverType for FIELD
00615   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(NO_DRIVER, filename_rd, fieldname_nodeint_rd)),
00616                        MEDEXCEPTION);
00617   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(GIBI_DRIVER, filename_rd, fieldname_nodeint_rd)),
00618                        MEDEXCEPTION);
00619   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(PORFLOW_DRIVER, filename_rd, fieldname_nodeint_rd)),
00620                        MEDEXCEPTION);
00621   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(ASCII_DRIVER, filename_rd, fieldname_nodeint_rd)),
00622                        MEDEXCEPTION);
00623 
00625   //TestRead Part//
00627   FIELD<double> *aField_1 = NULL;
00628   CPPUNIT_ASSERT_NO_THROW(aField_1 = new FIELD<double>(MED_DRIVER, filename_rd, fieldname_celldouble_rd));
00629 
00630   //Test read(int index) method
00631   int IdDriver_rd = aField_1->addDriver(MED_DRIVER,filename_rd,fieldname_celldouble_rd);
00632   //#ifdef ENABLE_FORCED_FAILURES
00633   // (BUG) Cannot open file, but file exist
00634   CPPUNIT_ASSERT_NO_THROW(aField_1->read(IdDriver_rd));
00635   //#endif
00636 
00637   //Test read(GENDRIVER & genDriver) method
00638   //Creation a Driver
00639   MED_FIELD_RDONLY_DRIVER<int> *aMedRdFieldDriver_1 = new MED_FIELD_RDONLY_DRIVER<int>();
00640   //Creation a Field
00641   FIELD<int> *aField_2 = new FIELD<int>();
00642   aField_2->setName(fieldname_nodeint_rd);
00643   aField_2->addDriver(*aMedRdFieldDriver_1);
00644   aField_2->read(*aMedRdFieldDriver_1);
00645 
00647   //Test Write Part//
00649   int IdDriver;
00650   MESH *aMesh = new MESH(MED_DRIVER,filename_rd,meshname);
00651   const SUPPORT *aSupport = aMesh->getSupportOnAll( MED_CELL );
00652   FIELD<int> *aFieldSupport;
00653   //#ifdef ENABLE_FORCED_FAILURES  
00654   CPPUNIT_ASSERT_NO_THROW(aFieldSupport = 
00655                           new FIELD<int>(aSupport, MED_DRIVER,filename_support_wr,fieldname_nodeint_rd));
00656   //(BUG) Can not open file
00657   MED_FIELD_WRONLY_DRIVER<int> * aFieldWrDriver = 
00658     new MED_FIELD_WRONLY_DRIVER<int>(filename_support_wr,aFieldSupport);
00659   aFieldWrDriver->setFieldName(aFieldSupport->getName() + "_copy");
00660   CPPUNIT_ASSERT_NO_THROW(IdDriver= aFieldSupport->addDriver(*aFieldWrDriver));
00661   CPPUNIT_ASSERT_NO_THROW(aFieldSupport->write(IdDriver));
00662   aFieldSupport->removeReference();
00663   delete aFieldWrDriver;
00664   //#endif
00665 
00666   //Create fileds
00667   FIELD<double> * aField_3 = new FIELD<double>();
00668   MED_FIELD_RDONLY_DRIVER<double> *aMedRdFieldDriver_2 =
00669     new MED_FIELD_RDONLY_DRIVER<double>(filename_rd, aField_3);
00670   aMedRdFieldDriver_2->open();
00671   aMedRdFieldDriver_2->setFieldName(fieldname_celldouble_rd);
00672   aMedRdFieldDriver_2->read();
00673   aMedRdFieldDriver_2->close();
00674 
00675   //Test write(int index) method
00676   //Add drivers to FIELDs
00677   int IdDriver1 = -1;
00678   try
00679     {
00680       IdDriver1 = aField_3->addDriver(MED_DRIVER,filename_wr,fieldname_celldouble_wr);
00681     }
00682   catch(MEDEXCEPTION &e)
00683     {
00684       e.what();
00685     }
00686   catch( ... )
00687     {
00688       CPPUNIT_FAIL("Unknown exception");
00689     }
00690   //Trying call write(int index) method with incorrect index
00691   //#ifdef ENABLE_FAULTS
00692   CPPUNIT_ASSERT_THROW(aField_3->write(IdDriver1+1),MEDEXCEPTION);
00693   // => Segmentation fault
00694   //#endif
00695 
00696   //Write field to file
00697   //#ifdef ENABLE_FAULTS
00698   try
00699     {
00700       aField_3->write(IdDriver1);
00701       // => Segmentation fault
00702     }
00703   catch(MEDEXCEPTION &e)
00704     {
00705       e.what();
00706     }
00707   catch( ... )
00708     {
00709       CPPUNIT_FAIL("Unknown exception");
00710     }
00711   //#endif
00712 
00713   CPPUNIT_ASSERT_NO_THROW(aField_3->rmDriver(IdDriver1));
00714 
00715   //Test write(const GENDRIVER &);
00716   //Create a driver
00717   MED_FIELD_WRONLY_DRIVER<int> *aMedWrFieldDriver = new MED_FIELD_WRONLY_DRIVER<int>();
00718   aMedWrFieldDriver->setFileName(filename_wr);
00719   aField_2->setName(fieldname_nodeint_wr1);
00720   //Add driver to a field
00721   aField_2->addDriver(*aMedWrFieldDriver);
00722 
00723   try
00724     {
00725       aField_2->write(*aMedWrFieldDriver);
00726     }
00727   catch(MEDEXCEPTION &e)
00728     {
00729       e.what();
00730     }
00731   catch( ... )
00732     {
00733       CPPUNIT_FAIL("Unknown exception");
00734     }
00735 
00736   //Test writeAppend(int index) method
00737   //Create a vtk file
00738   MESH * aMesh_1 = new MESH;
00739   MED_MESH_RDONLY_DRIVER *aMedMeshRdDriver22 = new MED_MESH_RDONLY_DRIVER(filename22_rd, aMesh_1);
00740   aMedMeshRdDriver22->open();
00741   aMedMeshRdDriver22->setMeshName(meshname);
00742   aMedMeshRdDriver22->read();
00743   aMedMeshRdDriver22->close();
00744   VTK_MESH_DRIVER *aVtkDriver = new VTK_MESH_DRIVER(filenamevtk_wr, aMesh_1);
00745   aVtkDriver->open();
00746   aVtkDriver->write();
00747   aVtkDriver->close();
00748 
00749   //Create a field
00750   FIELD<int> * aField_4 = new FIELD<int>();
00751   MED_FIELD_RDONLY_DRIVER<int> *aMedRdFieldDriver22 =
00752     new MED_FIELD_RDONLY_DRIVER<int>(filename22_rd, aField_2);
00753   aMedRdFieldDriver22->open();
00754   aMedRdFieldDriver22->setFieldName(fieldname_nodeint_rd);
00755   aMedRdFieldDriver22->read();
00756   aMedRdFieldDriver22->close();
00757 
00758   //Add Driver to a field
00759   int IdDriver2;
00760   try
00761     {
00762       IdDriver2 = aField_4->addDriver(VTK_DRIVER, filenamevtk_wr ,fieldname_nodeint_wr);
00763     }
00764   catch(MEDEXCEPTION &e)
00765     {
00766       e.what();
00767     }
00768   catch( ... )
00769     {
00770       CPPUNIT_FAIL("Unknown exception");
00771     }
00772   //#ifdef ENABLE_FAULTS
00773   //Trying call writeAppend() method with incorrect index
00774   CPPUNIT_ASSERT_THROW(aField_4->writeAppend(IdDriver2+1,fieldname_nodeint_wr),MEDEXCEPTION);
00775   // => Segmentation fault
00776   //#endif
00777 
00778   //#ifdef ENABLE_FAULTS
00779   // (BUG) => Segmentation fault
00780   CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(IdDriver2, fieldname_nodeint_wr));
00781   //#endif
00782 
00783   //Test writeAppend(const GENDRIVER &) method
00784   aField_4->setName(fieldname_nodeint_wr1);
00785 
00786   //Add driver to a field
00787   //#ifdef ENABLE_FAULTS
00788   //Create a driver
00789   VTK_FIELD_DRIVER<int> *aVtkFieldDriver = new VTK_FIELD_DRIVER<int>(filenamevtk_wr, aField_4);
00790   CPPUNIT_ASSERT_NO_THROW(aField_4->addDriver(*aVtkFieldDriver));
00791   //(BUG) => Segmentation fault after addDriver(const GENDRIVER &)
00792   CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(*aVtkFieldDriver));
00793   //#endif
00794 
00795 
00796   //Delete objects
00797   aField_1->removeReference();
00798   delete aMedRdFieldDriver_1;
00799   aField_2->removeReference();
00800   aField_3->removeReference();
00801   delete aMedRdFieldDriver_2;
00802   aField_4->removeReference();
00803   delete aMedMeshRdDriver22;
00804   delete aMedWrFieldDriver;
00805   delete aVtkDriver;
00806   aMesh->removeReference();
00807   aMesh_1->removeReference();
00808   delete aMedRdFieldDriver22;
00809 }
00810 
00811 static void MEDMEMTest_testField()
00812 {
00813   SUPPORT *anEmptySupport=new SUPPORT;
00815   // TEST 1: FIELD_ //
00817   FIELD_ *aField_=new FIELD_ ;
00818 
00819   // check set/get methods
00820   MED_EN::med_type_champ aValueType = MED_EN::MED_UNDEFINED_TYPE;
00821   MED_EN::medModeSwitch  anInterlace = MED_EN::MED_UNDEFINED_INTERLACE;
00822   checkField_(aField_, anEmptySupport, aValueType, anInterlace);
00823 
00824   // copy constructor
00825   // This fails (Segmentation fault) if not set:
00826   // _componentsNames or _componentsDescriptions, or _MEDComponentsUnits
00827   FIELD_ *aField_copy1=new FIELD_(*aField_);
00828   compareField_(aField_, aField_copy1, /*isFIELD = */false, /*isValue = */false);
00829 
00830   // operator=
00831   //#ifdef ENABLE_FAULTS
00832   // (BUG) This fails (Segmentation fault) if not set:
00833   // _componentsNames or _componentsDescriptions, or _componentsUnits, or _MEDComponentsUnits
00834   // (BUG) Code duplication with copyGlobalInfo(), called from copy constructor
00835   FIELD_ *aField_copy2=new FIELD_;
00836   *aField_copy2 = *aField_;
00837   compareField_(aField_, aField_copy2, /*isFIELD = */false, /*isValue = */false);
00838   //#endif
00839   //#ifdef ENABLE_FORCED_FAILURES
00840   //  CPPUNIT_FAIL("FIELD_::operator=() fails if _componentsUnits is not set");
00841   //#endif
00842 
00843   // following test is commented since method
00844   // setTotalNumberOfElements() is removed.
00845   /*
00846   // construction on a given support
00847   {
00848   anEmptySupport.setTotalNumberOfElements(11);
00849   // CASE1:
00850   FIELD_ aField_case1 (&anEmptySupport, 10);
00851   // CASE2:
00852   FIELD_ aField_case2;
00853   aField_case2.setSupport(&anEmptySupport);
00854   aField_case2.setNumberOfComponents(10);
00855 
00856   #ifdef ENABLE_FORCED_FAILURES
00857   CPPUNIT_ASSERT_EQUAL_MESSAGE("No correspondance between CASE1 and CASE2",
00858   aField_case1.getNumberOfValues(),
00859   aField_case2.getNumberOfValues());
00860   #endif
00861   }
00862   */
00863 
00865   // TEST 2: FIELD<int> //
00867   FIELD<int> *aFieldInt=new FIELD<int>();
00868   checkField(aFieldInt, anEmptySupport);
00869 
00871   // TEST 3: FIELD<double, NoInterlace> //
00873   MESH * aMesh = MEDMEMTest_createTestMesh();
00874   const GROUP* aGroup = aMesh->getGroup(MED_EN::MED_FACE, 1);
00875 
00876   FIELD<double, NoInterlace> aFieldDouble;
00877   checkField(&aFieldDouble, aGroup);
00878 
00880   // TEST 4: FIELD<double, FullInterlace> //
00882   FIELD<double> * aFieldOnGroup1 = createFieldOnGroup<double>(aMesh, aGroup, "Linear", "N");
00883   FIELD<double> * aFieldOnGroup2 = createFieldOnGroup<double>(aMesh, aGroup, "Quadratic", "N**2");
00884 
00885   int nbVals = aFieldOnGroup1->getNumberOfValues();
00886   CPPUNIT_ASSERT(nbVals);
00887 
00888   // numbers of elements in group,
00889   // they are needed in method FIELD::setValueIJ()
00890   const int *anElems = aGroup->getnumber()->getValue();
00891   double eucl1 = 0., eucl2 = 0.;
00892 
00893   for (int i = 1; i <= nbVals; i++) 
00894     {
00895       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, (double)i);
00896       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, (double)(-i));
00897 
00898       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, (double)i*i);
00899       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, (double)(-i*i));
00900 
00901       eucl1 += 2. * i * i;
00902       eucl2 += 2. * i * i * i * i;
00903     }
00904 
00905   // out of bound (inexisting 33-th component of last element)
00906   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->setValueIJ(anElems[nbVals-1], 33, 10.), MEDEXCEPTION);
00907 
00908   // normMax
00909   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals, aFieldOnGroup1->normMax(), 0.000001);
00910   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals*nbVals, aFieldOnGroup2->normMax(), 0.000001);
00911 
00912   // norm2
00913   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl1), aFieldOnGroup1->norm2(), 0.000001); // 10.4881
00914   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl2), aFieldOnGroup2->norm2(), 0.000001); // 44.2493
00915 
00916   // check getXXX methods
00917   CPPUNIT_ASSERT(!aFieldOnGroup1->getGaussPresence());
00918   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getArrayGauss(), MEDEXCEPTION);
00919   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * anArrayNoGauss =
00920     aFieldOnGroup1->getArrayNoGauss();
00921 
00922   MEDMEM_Array_ * aMEDMEM_Array_ = aFieldOnGroup1->getArray();
00923   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * aMEDMEM_Array_conv =
00924     static_cast<MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array *>(aMEDMEM_Array_);
00925 
00926   const double * aValues = aFieldOnGroup1->getValue();
00927 
00928   // out of range
00929   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(3), MEDEXCEPTION);
00930   // cannot get column in FullInterlace
00931   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(1), MEDEXCEPTION);
00932 
00933   for (int i = 1; i <= nbVals; i++) 
00934     {
00935       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aFieldOnGroup1->getValueIJK(anElems[i-1], 1, 1), 0.000001);
00936       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aFieldOnGroup1->getValueIJK(anElems[i-1], 2, 1), 0.000001);
00937 
00938       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aValues[(i-1)*2 + 0], 0.000001);
00939       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aValues[(i-1)*2 + 1], 0.000001);
00940 
00941       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , anArrayNoGauss->getIJ(i, 1), 0.000001);
00942       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), anArrayNoGauss->getIJ(i, 2), 0.000001);
00943 
00944       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aMEDMEM_Array_conv->getIJ(i, 1), 0.000001);
00945       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aMEDMEM_Array_conv->getIJ(i, 2), 0.000001);
00946 
00947       const double* row_i = aFieldOnGroup1->getRow(anElems[i-1]);
00948       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , row_i[0], 0.000001);
00949       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), row_i[1], 0.000001);
00950 
00951       double vals_i [2];
00952       aFieldOnGroup1->getValueOnElement(anElems[i-1], vals_i);
00953       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , vals_i[0], 0.000001);
00954       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), vals_i[1], 0.000001);
00955     }
00956 
00957   // modify all values of aFieldOnGroup2 by formula a*x + b (a = 2, b = 3)
00958   aFieldOnGroup2->applyLin(2., 3.);
00959   for (int i = 1; i <= nbVals; i++) 
00960     {
00961       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. + 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
00962       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. - 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
00963     }
00964 
00965   // apply function plus13() to aFieldOnGroup1
00966   aFieldOnGroup1->applyFunc<plus13>();
00967   for (int i = 1; i <= nbVals; i++) 
00968     {
00969       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
00970       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
00971     }
00972 
00973   // scalarProduct
00974   FIELD<double, FullInterlace> * aScalarProduct =
00975     FIELD<double, FullInterlace>::scalarProduct(*aFieldOnGroup1, *aFieldOnGroup2, /*deepCheck = */true);
00976   CPPUNIT_ASSERT_EQUAL(nbVals, aScalarProduct->getNumberOfValues());
00977   CPPUNIT_ASSERT_EQUAL(1, aScalarProduct->getNumberOfComponents());
00978   for (int i = 1; i <= nbVals; i++) 
00979     {
00980       CPPUNIT_ASSERT_DOUBLES_EQUAL(78. + 4.*i*i*i, //(3. + 2.*i*i)*(13 + i) + (3. - 2.*i*i)*(13 - i)
00981                                    aScalarProduct->getValueIJ(anElems[i-1], 1), 0.000001);
00982     }
00983 
00984   // fillFromAnalytic
00985   aFieldOnGroup2->fillFromAnalytic(proj2d);
00986 
00987   double bary [3];
00988   double outp [2];
00989   const SUPPORT * aSupp = aFieldOnGroup2->getSupport();
00990   FIELD<double,FullInterlace> * barycenter = aMesh->getBarycenter(aSupp);
00991   for (int i = 1; i <= nbVals; i++) 
00992     {
00993       bary[0] = barycenter->getValueIJ(anElems[i-1], 1);
00994       bary[1] = barycenter->getValueIJ(anElems[i-1], 2);
00995       bary[2] = barycenter->getValueIJ(anElems[i-1], 3);
00996 
00997       proj2d(bary, outp);
00998 
00999       //cout << "barycenter (" << bary[0] << ", " << bary[1] << ", " << bary[2] << ")" << endl;
01000       //cout << "proj2d     (" << outp[0] << ", " << outp[1] << ")" << endl;
01001 
01002       //bary (-0.666667,  0.666667, 0.666667) -> outp ( 0.902369, 0.902369)
01003       //bary ( 0.666667, -0.666667, 0.666667) -> outp (-0.902369, 0.430964)
01004       //bary ( 0.      ,  0.      , 2.      ) -> outp ( 0.      , 2.      )
01005       //bary ( 0.      ,  0.      , 3.      ) -> outp ( 0.      , 3.      )
01006       //bary (-1.      ,  0.      , 2.5     ) -> outp ( 0.353553, 2.85355 )
01007 
01008       //#ifdef ENABLE_FORCED_FAILURES
01009       // (BUG) in FIELD::fillFromAnalytic() in case of support, different from nodes:
01010       //       barycenterField in FullInterlace, but values extracted like from NoInterlace
01011       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[0], aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
01012       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[1], aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
01013       //#endif
01014 
01015       // currently it gives values, that are wrong:
01016       //aFieldOnGroup2 row1 ( 0.902369,  0.235702)
01017       //aFieldOnGroup2 row2 (-0.235702,  2.7643  )
01018       //aFieldOnGroup2 row3 (-0.235702, -1.2357  )
01019       //aFieldOnGroup2 row4 ( 1.7643  , -0.235702)
01020       //aFieldOnGroup2 row5 ( 0.235702,  2.7357  )
01021     }
01022 
01023   // info about support (Group1)
01024   CPPUNIT_ASSERT(!aFieldOnGroup1->isOnAllElements()); // because we build Group1 so
01025   int nbTypes = aFieldOnGroup1->getNumberOfGeometricTypes();
01026   //CPPUNIT_ASSERT(nbTypes);
01027   CPPUNIT_ASSERT_EQUAL(2, nbTypes);
01028   const int * nbElemsInEachType = aFieldOnGroup1->getNumberOfElements();
01029   const MED_EN::medGeometryElement * aGeomTypes = aFieldOnGroup1->getGeometricTypes();
01030 
01031   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_TRIA3, aGeomTypes[0]);
01032   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_QUAD4, aGeomTypes[1]);
01033 
01034   // GAUSS
01035 
01036   // now we have no gauss localization in aFieldOnGroup1
01037   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
01038   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
01039   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
01040   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
01041 
01042   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3), MEDEXCEPTION);
01043   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3), MEDEXCEPTION);
01044 
01045   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
01046 
01047   // set a gauss localization into aFieldOnGroup1
01048   double cooRef[6] = 
01049     {
01050       1.,1., 2.,4., 3.,9.
01051     }; // xy xy xy
01052   double cooGauss[10] = 
01053     {
01054       7.,7., 6.,6., 5.,5., 4.,3., 2.,1.
01055     }; // x1,y1  x2,y2  x3,y3  x4,y4  x5,y5
01056   double wg[5] = 
01057     {
01058       1., 2., 3., 4., 5.
01059     };
01060   GAUSS_LOCALIZATION<> gl1 ("GL1", MED_EN::MED_TRIA3, /*nGauss*/5, cooRef, cooGauss, wg);
01061 
01062   aFieldOnGroup1->setGaussLocalization(MED_EN::MED_TRIA3, gl1);
01063 
01064   // now we have a gauss localization for MED_TRIA3 type
01065   CPPUNIT_ASSERT_EQUAL(5, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
01066   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
01067   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
01068   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
01069 
01070   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_QUAD4), MEDEXCEPTION);
01071   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_QUAD4), MEDEXCEPTION);
01072 
01073   GAUSS_LOCALIZATION<> gl1Back = aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3);
01074   const GAUSS_LOCALIZATION<> * gl1BackPtr = aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3);
01075 
01076   CPPUNIT_ASSERT(gl1 == gl1Back);
01077   CPPUNIT_ASSERT(gl1 == *gl1BackPtr);
01078 
01079   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
01080 
01081   // sub-support of Group1 on one (first) geometric type
01082   SUPPORT * aSubSupport1 = new SUPPORT;
01083   aSubSupport1->setMesh( aMesh );
01084   aSubSupport1->setName( "Sub-Support 1 of Group1" );
01085   aSubSupport1->setEntity( MED_EN::MED_FACE );
01086 
01087   int nbTypes1 = 1;
01088   int nbElemsInEachType1[1];
01089   nbElemsInEachType1[0] = nbElemsInEachType[0];
01090   int nbElems1 = nbElemsInEachType1[0];
01091   MED_EN::medGeometryElement aGeomTypes1[1];
01092   aGeomTypes1[0] = aGeomTypes[0];
01093   int * anElems1 = new int[nbElems1];
01094   for (int i = 0; i < nbElems1; i++) 
01095     {
01096       anElems1[i] = anElems[i];
01097     }
01098 
01099   aSubSupport1->setpartial("Support for sub-field 1 on one type of elements",
01100                            nbTypes1, nbElems1, aGeomTypes1, nbElemsInEachType1, anElems1);
01101 
01102   //cout << "aSubSupport1:" << endl;
01103   //cout << *aSubSupport1 << endl;
01104 
01105   // extract sub-field on aSubSupport1
01106   FIELD<double, FullInterlace> * aSubField1 = aFieldOnGroup1->extract(aSubSupport1);
01107   CPPUNIT_ASSERT_EQUAL(nbElems1 * /*NumberOfComponents = */2, aSubField1->getValueLength());
01108 
01109   // aSubField1:
01110   // elt\comp |  1 |  2
01111   //--------------------
01112   //  1       | 14 | 12
01113   //  2       | 15 | 11
01114 
01115   // check normL2() and normL1()
01116   FIELD<double>* anAreaField = aMesh->getArea(aSubSupport1);
01117   double area1 = anAreaField->getValueIJ(anElems1[0], 1);
01118   double area2 = anAreaField->getValueIJ(anElems1[1], 1);
01119   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area1, 0.00001);
01120   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area2, 0.00001);
01121 
01122   CPPUNIT_ASSERT_DOUBLES_EQUAL(210.5, aSubField1->normL2(1), 0.00001); // (14*14 + 15*15)/2
01123   //#ifdef ENABLE_FORCED_FAILURES
01124   // (BUG) FIELD::normL2(int component, const FIELD * p_field_volume):
01125   //       component is not taken into account
01126   CPPUNIT_ASSERT_DOUBLES_EQUAL(132.5, aSubField1->normL2(2), 0.00001); // (12*12 + 11*11)/2
01127   //#endif
01128   CPPUNIT_ASSERT_DOUBLES_EQUAL(343.0, aSubField1->normL2() , 0.00001); // 210.5 + 132.5
01129 
01130   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.5, aSubField1->normL1(1), 0.00001); // (14 + 15)/2
01131   CPPUNIT_ASSERT_DOUBLES_EQUAL(11.5, aSubField1->normL1(2), 0.00001); // (12 + 11)/2
01132   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1() , 0.00001); // 14.5 + 11.5
01133 
01134   double aNewArea [2] = 
01135     {
01136       1., 0.
01137     }; // only first element will be taken into account
01138   anAreaField->setValue(aNewArea);
01139 
01140   CPPUNIT_ASSERT_DOUBLES_EQUAL(196.0, aSubField1->normL2(1, anAreaField), 0.00001); // 14*14/1
01141   //#ifdef ENABLE_FORCED_FAILURES
01142   // (BUG) FIELD::normL2(int component, const FIELD * p_field_volume):
01143   //       component is not taken into account
01144   CPPUNIT_ASSERT_DOUBLES_EQUAL(144.0, aSubField1->normL2(2, anAreaField), 0.00001); // 12*12/1
01145   //#endif
01146   CPPUNIT_ASSERT_DOUBLES_EQUAL(340.0, aSubField1->normL2(anAreaField) , 0.00001); // 196 + 144
01147 
01148   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.0, aSubField1->normL1(1, anAreaField), 0.00001); // 14/1
01149   CPPUNIT_ASSERT_DOUBLES_EQUAL(12.0, aSubField1->normL1(2, anAreaField), 0.00001); // 12/1
01150   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1(anAreaField) , 0.00001); // 14 + 12
01151 
01152   // applyPow
01153   aSubField1->applyPow(2.);
01154   CPPUNIT_ASSERT_DOUBLES_EQUAL(196., aSubField1->getValueIJ(anElems1[0], 1), 0.000001); // 14*14
01155   CPPUNIT_ASSERT_DOUBLES_EQUAL(144., aSubField1->getValueIJ(anElems1[0], 2), 0.000001); // 12*12
01156   CPPUNIT_ASSERT_DOUBLES_EQUAL(225., aSubField1->getValueIJ(anElems1[1], 1), 0.000001); // 15*15
01157   CPPUNIT_ASSERT_DOUBLES_EQUAL(121., aSubField1->getValueIJ(anElems1[1], 2), 0.000001); // 11*11
01158 
01159   // setArray (NoGauss)
01160   MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array * aNewArrayNoGauss =
01161     new MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array(/*dim*/2, /*nbelem*/2);
01162   aNewArrayNoGauss->setIJ(1, 1, 4.);
01163   aNewArrayNoGauss->setIJ(1, 2, 2.);
01164   aNewArrayNoGauss->setIJ(2, 1, 5.);
01165   aNewArrayNoGauss->setIJ(2, 2, 1.);
01166   aSubField1->setArray(aNewArrayNoGauss);
01167   // no need to delete aNewArrayNoGauss, because it will be deleted
01168   // in destructor or in deallocValue() method of aSubField1
01169 
01170   CPPUNIT_ASSERT_DOUBLES_EQUAL(4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01171   CPPUNIT_ASSERT_DOUBLES_EQUAL(2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01172   CPPUNIT_ASSERT_DOUBLES_EQUAL(5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01173   CPPUNIT_ASSERT_DOUBLES_EQUAL(1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01174 
01175   // setRow
01176   double row[2] = 
01177     {
01178       -1., -3.
01179     };
01180   aSubField1->setRow(anElems1[0], row);
01181   CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01182   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01183   CPPUNIT_ASSERT_DOUBLES_EQUAL( 5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01184   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01185   // out of range
01186   CPPUNIT_ASSERT_THROW(aSubField1->setRow(3, row), MEDEXCEPTION);
01187 
01188   // setColumn
01189   double col[2] = 
01190     {
01191       -7., -9.
01192     };
01193   aSubField1->setColumn(1, col);
01194   CPPUNIT_ASSERT_DOUBLES_EQUAL(-7., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01195   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01196   //#ifdef ENABLE_FORCED_FAILURES
01197   // (BUG) in MEDMEM_Array::setColumn()
01198   CPPUNIT_ASSERT_DOUBLES_EQUAL(-9., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01199   //#endif
01200   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01201   // out of range
01202   CPPUNIT_ASSERT_THROW(aSubField1->setColumn(3, col), MEDEXCEPTION);
01203 
01204   // setArray (Gauss)
01205   {
01206     int nbelgeoc[2] = 
01207       {
01208         1, 3
01209       }; // 3 - 1 = two elements for the first (and the only) type
01210     int nbgaussgeo[2] = 
01211       {
01212         -1, 1
01213       }; // one gauss point per each element
01214     MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array * aNewArrayGauss =
01215       new MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array
01216       (/*dim*/2, /*nbelem*/2, /*nbtypegeo*/1, /*nbelgeoc*/nbelgeoc, /*nbgaussgeo*/nbgaussgeo);
01217 
01218     //#ifdef ENABLE_FAULTS
01219     aNewArrayGauss->setIJ(1, 1, -4.);
01220     aNewArrayGauss->setIJ(1, 2, -2.);
01221     aNewArrayGauss->setIJ(2, 1, -5.);
01222     aNewArrayGauss->setIJ(2, 2, -1.);
01223     //#endif
01224     //#ifdef ENABLE_FORCED_FAILURES
01225     // ? (BUG) in FullInterlaceGaussPolicy::getIndex(int i,int j)
01226     // FullInterlaceGaussPolicy::getIndex(2,2) returns 4!!!
01227     //    CPPUNIT_FAIL("? Bug in FullInterlaceGaussPolicy::getIndex(int i,int j) ?");
01228     //#endif
01229 
01230     aNewArrayGauss->setIJK(1, 1, 1, -4.);
01231     aNewArrayGauss->setIJK(1, 2, 1, -2.);
01232     aNewArrayGauss->setIJK(2, 1, 1, -5.);
01233     aNewArrayGauss->setIJK(2, 2, 1, -1.);
01234 
01235     aSubField1->setArray(aNewArrayGauss);
01236     // no need to delete aNewArrayGauss, because it will be deleted
01237     // in destructor or in deallocValue() method of aSubField1
01238 
01239     //#ifdef ENABLE_FAULTS
01240     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
01241     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
01242     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
01243     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
01244     //#endif
01245     //#ifdef ENABLE_FORCED_FAILURES
01246     // ? (BUG) in FullInterlaceGaussPolicy::getIndex(int i,int j)
01247     // Must be   : return _G[i-1]-1 + (j-1);
01248     // Instead of: return _G[i-1]-1 + (j-1)*_dim;
01249     //    CPPUNIT_FAIL("? Bug in FullInterlaceGaussPolicy::getIndex(int i,int j) ?");
01250     //#endif
01251 
01252     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJK(anElems1[0], 1, 1), 0.000001);
01253     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJK(anElems1[0], 2, 1), 0.000001);
01254     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJK(anElems1[1], 1, 1), 0.000001);
01255     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJK(anElems1[1], 2, 1), 0.000001);
01256   }
01257 
01258   // alloc/dealloc; compatibility of new size with support
01259   try
01260     {
01261       aSubField1->deallocValue();
01262       aSubField1->allocValue(/*NumberOfComponents*/2, /*LengthValue*/5);
01263       //#ifdef ENABLE_FAULTS
01264       // (BUG) No compatibility between Support and allocated value
01265       aSubField1->normL1();
01266       //#endif
01267       //#ifdef ENABLE_FORCED_FAILURES
01268       //    CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
01269       //#endif
01270     }
01271   catch (MEDEXCEPTION & ex) 
01272     {
01273       // normal behaviour
01274     }
01275   catch (...) 
01276     {
01277       CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
01278     }
01279 
01280   // check that aFieldOnGroup1 is not changed after aSubField1 modifications
01281   for (int i = 1; i <= nbVals; i++) 
01282     {
01283       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01284       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01285     }
01286 
01287   // reset aFieldOnGroup2 values for simple control of operators results
01288   for (int i = 1; i <= nbVals; i++) 
01289     {
01290       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, i*i);
01291       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, -i*i);
01292     }
01293 
01294   int len = aFieldOnGroup1->getValueLength();
01295   const double * val1 = aFieldOnGroup1->getValue();
01296   const double * val2 = aFieldOnGroup2->getValue();
01297   const double * val_res;
01298 
01299   // operators and add, sub, mul, div
01300 
01301   // +
01302   FIELD<double> *aSum = *aFieldOnGroup1 + *aFieldOnGroup2;
01303   aSum->setName(aFieldOnGroup1->getName());
01304   aSum->setDescription(aFieldOnGroup1->getDescription());
01305   compareField_(aFieldOnGroup1, aSum, true, true);
01306   val_res = aSum->getValue();
01307   for (int i = 0; i < len; i++) 
01308     {
01309       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
01310     }
01311   aSum->removeReference();
01312 
01313   // -
01314   FIELD<double> *aDifference = *aFieldOnGroup1 - *aFieldOnGroup2;
01315   aDifference->setName(aFieldOnGroup1->getName());
01316   aDifference->setDescription(aFieldOnGroup1->getDescription());
01317   compareField_(aFieldOnGroup1, aDifference, true, true);
01318   val_res = aDifference->getValue();
01319   for (int i = 0; i < len; i++) 
01320     {
01321       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
01322     }
01323   aDifference->removeReference();
01324 
01325   // - (unary)
01326   FIELD<double> *aNegative = - *aFieldOnGroup1;
01327   aNegative->setName(aFieldOnGroup1->getName());
01328   aNegative->setDescription(aFieldOnGroup1->getDescription());
01329   compareField_(aFieldOnGroup1, aNegative, true, true);
01330   val_res = aNegative->getValue();
01331   for (int i = 0; i < len; i++) 
01332     {
01333       CPPUNIT_ASSERT_DOUBLES_EQUAL(- val1[i], val_res[i], 0.000001);
01334     }
01335   aNegative->removeReference();
01336 
01337   // *
01338   FIELD<double> *aProduct = (*aFieldOnGroup1) * (*aFieldOnGroup2);
01339   aProduct->setName(aFieldOnGroup1->getName());
01340   aProduct->setDescription(aFieldOnGroup1->getDescription());
01341   compareField_(aFieldOnGroup1, aProduct, true, true);
01342   val_res = aProduct->getValue();
01343   for (int i = 0; i < len; i++) 
01344     {
01345       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
01346     }
01347   aProduct->removeReference();
01348   // /
01349   FIELD<double> *aQuotient = *aFieldOnGroup1 / *aFieldOnGroup2;
01350   aQuotient->setName(aFieldOnGroup1->getName());
01351   aQuotient->setDescription(aFieldOnGroup1->getDescription());
01352   compareField_(aFieldOnGroup1, aQuotient, true, true);
01353   val_res = aQuotient->getValue();
01354   for (int i = 0; i < len; i++) 
01355     {
01356       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
01357     }
01358   aQuotient->removeReference();
01359 
01360   double val22 = aFieldOnGroup2->getValueIJ(anElems[2], 2);
01361   aFieldOnGroup2->setValueIJ(anElems[2], 2, 0.);
01362 
01363   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup2, MEDEXCEPTION);
01364   //#ifdef ENABLE_FORCED_FAILURES
01365   // (BUG) is it up to user to control validity of data to avoid division on zero?
01366   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup2, MEDEXCEPTION);
01367   //#endif
01368   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01369   CPPUNIT_ASSERT_THROW(FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01370 
01371   // restore value
01372   aFieldOnGroup2->setValueIJ(anElems[2], 2, val22);
01373 
01374   // restore values
01375   for (int i = 1; i <= nbVals; i++) 
01376     {
01377       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, 13 + i);
01378       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, 13 - i);
01379     }
01380 
01381   // static methods
01382   FIELD<double> * aPtr;
01383 
01384   // add
01385   aPtr = FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup2);
01386   aPtr->setName(aFieldOnGroup1->getName());
01387   aPtr->setDescription(aFieldOnGroup1->getDescription());
01388   compareField_(aFieldOnGroup1, aPtr, true, true);
01389   val_res = aPtr->getValue();
01390   for (int i = 0; i < len; i++) 
01391     {
01392       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
01393     }
01394   delete aPtr;
01395 
01396   // sub
01397   aPtr = FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2);
01398   aPtr->setName(aFieldOnGroup1->getName());
01399   aPtr->setDescription(aFieldOnGroup1->getDescription());
01400   compareField_(aFieldOnGroup1, aPtr, true, true);
01401   val_res = aPtr->getValue();
01402   for (int i = 0; i < len; i++) 
01403     {
01404       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
01405     }
01406   delete aPtr;
01407 
01408   // mul
01409   aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
01410   aPtr->setName(aFieldOnGroup1->getName());
01411   aPtr->setDescription(aFieldOnGroup1->getDescription());
01412   compareField_(aFieldOnGroup1, aPtr, true, true);
01413   val_res = aPtr->getValue();
01414   for (int i = 0; i < len; i++) 
01415     {
01416       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
01417     }
01418   delete aPtr;
01419 
01420   // div
01421   aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
01422   aPtr->setName(aFieldOnGroup1->getName());
01423   aPtr->setDescription(aFieldOnGroup1->getDescription());
01424   compareField_(aFieldOnGroup1, aPtr, true, true);
01425   val_res = aPtr->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   delete aPtr;
01431 
01432   // addDeep
01433   aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01434   aPtr->setName(aFieldOnGroup1->getName());
01435   aPtr->setDescription(aFieldOnGroup1->getDescription());
01436   compareField_(aFieldOnGroup1, aPtr, true, true);
01437   val_res = aPtr->getValue();
01438   for (int i = 0; i < len; i++) 
01439     {
01440       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
01441     }
01442   delete aPtr;
01443 
01444   // subDeep
01445   aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01446   aPtr->setName(aFieldOnGroup1->getName());
01447   aPtr->setDescription(aFieldOnGroup1->getDescription());
01448   compareField_(aFieldOnGroup1, aPtr, true, true);
01449   val_res = aPtr->getValue();
01450   for (int i = 0; i < len; i++) 
01451     {
01452       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
01453     }
01454   delete aPtr;
01455 
01456   // mulDeep
01457   aPtr = FIELD<double>::mulDeep(*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   delete aPtr;
01467 
01468   // divDeep
01469   aPtr = FIELD<double>::divDeep(*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   delete aPtr;
01479 
01480   // +=
01481   *aFieldOnGroup1 += *aFieldOnGroup2;
01482   for (int i = 1; i <= nbVals; i++) 
01483     {
01484       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i + i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01485       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i - i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01486     }
01487 
01488   // -=
01489   *aFieldOnGroup1 -= *aFieldOnGroup2;
01490   for (int i = 1; i <= nbVals; i++) 
01491     {
01492       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01493       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01494     }
01495 
01496   // *=
01497   *aFieldOnGroup1 *= *aFieldOnGroup2;
01498   for (int i = 1; i <= nbVals; i++) 
01499     {
01500       CPPUNIT_ASSERT_DOUBLES_EQUAL( (13 + i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01501       CPPUNIT_ASSERT_DOUBLES_EQUAL(-(13 - i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01502     }
01503 
01504   // /=
01505   *aFieldOnGroup1 /= *aFieldOnGroup2;
01506   for (int i = 1; i <= nbVals; i++) 
01507     {
01508       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
01509       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
01510     }
01511 
01512   // check case of different operands: support
01513   MESH * aMeshOneMore = MEDMEMTest_createTestMesh();
01514   const GROUP* aGroupOneMore = aMeshOneMore->getGroup(MED_EN::MED_FACE, 1);
01515   FIELD<double> * aFieldOnGroup3 =
01516     createFieldOnGroup<double>(aMeshOneMore, aGroupOneMore, "Test_Diff_Mesh", "test");
01517   for (int i = 1; i <= nbVals; i++) 
01518     {
01519       aFieldOnGroup3->setValueIJ(anElems[i-1], 1, 2*i);
01520       aFieldOnGroup3->setValueIJ(anElems[i-1], 2, 3*i);
01521     }
01522   const double * val3 = aFieldOnGroup3->getValue();
01523 
01524   //CPPUNIT_ASSERT_NO_THROW();
01525   try
01526     {
01527       aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01528       aPtr->setName(aFieldOnGroup1->getName());
01529       aPtr->setDescription(aFieldOnGroup1->getDescription());
01530       compareField_(aFieldOnGroup1, aPtr, true, true);
01531       val_res = aPtr->getValue();
01532       for (int i = 0; i < len; i++) 
01533         {
01534           CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val3[i], val_res[i], 0.000001);
01535         }
01536       delete aPtr;
01537 
01538       aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01539       delete aPtr;
01540       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01541       delete aPtr;
01542       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup3);
01543       delete aPtr;
01544     }
01545   catch (MEDEXCEPTION & ex) 
01546     {
01547       CPPUNIT_FAIL(ex.what());
01548     }
01549   catch (...) 
01550     {
01551       CPPUNIT_FAIL("Unknown exception in FIELD::xxxDeep()");
01552     }
01553 
01554   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01555   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01556   CPPUNIT_ASSERT_THROW(FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01557   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
01558 
01559   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup3, MEDEXCEPTION);
01560   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup3, MEDEXCEPTION);
01561   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 * *aFieldOnGroup3, MEDEXCEPTION);
01562   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup3, MEDEXCEPTION);
01563 
01564   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup3, MEDEXCEPTION);
01565   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup3, MEDEXCEPTION);
01566   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 *= *aFieldOnGroup3, MEDEXCEPTION);
01567   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup3, MEDEXCEPTION);
01568 
01569   // check case of different operands: MEDComponentsUnits
01570   aFieldOnGroup1->setMEDComponentUnit(1, "unit3");
01571 
01572   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup2, MEDEXCEPTION);
01573   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup2, MEDEXCEPTION);
01574   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup2, MEDEXCEPTION);
01575   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup2, MEDEXCEPTION);
01576   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01577   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01578   CPPUNIT_ASSERT_THROW(FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01579   CPPUNIT_ASSERT_THROW(FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
01580 
01581   //CPPUNIT_ASSERT_NO_THROW();
01582   try
01583     {
01584       aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
01585       delete aPtr;
01586       aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
01587       delete aPtr;
01588       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01589       delete aPtr;
01590       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2);
01591       delete aPtr;
01592 
01593       *aFieldOnGroup1 *= *aFieldOnGroup2;
01594       *aFieldOnGroup1 /= *aFieldOnGroup2;
01595 
01596       FIELD<double> *aPr = *aFieldOnGroup1 * *aFieldOnGroup2;
01597       FIELD<double> *aQu = *aFieldOnGroup1 / *aFieldOnGroup2;
01598       aPr->removeReference();
01599       aQu->removeReference();
01600     }
01601   catch (MEDEXCEPTION & ex) 
01602     {
01603       CPPUNIT_FAIL(ex.what());
01604     }
01605   catch (...) 
01606     {
01607       CPPUNIT_FAIL("Unknown exception");
01608     }
01609 
01610   // restore MED units
01611   aFieldOnGroup1->setMEDComponentUnit(1, "unit1");
01612 
01613   // check case of different operands: numberOfComponents
01614   //#ifdef ENABLE_FAULTS
01615   // (BUG) Cannot allocate value of higher dimension because of _componentsTypes reinitialization
01616   // Must be MEDEXCEPTION instead. And on attempt to change nb.components must be the same behaviour.
01617   aFieldOnGroup1->deallocValue();
01618   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->allocValue(/*dim*/5), MEDEXCEPTION);
01619   //#endif
01620 
01621   aSubSupport1->removeReference();
01622   delete [] anElems1;
01623 
01624   delete aScalarProduct;
01625   delete aSubField1;
01626   delete anAreaField;
01627   delete barycenter;
01628   delete aFieldOnGroup1;
01629   delete aFieldOnGroup2;
01630   delete aFieldOnGroup3;
01631 
01632   delete aMesh;
01633   delete aMeshOneMore;
01634 
01636   // TEST 5: Drivers //
01638   testDrivers();
01639 }
01640 
01641 int main (int argc, char** argv)
01642 {
01643   MEDMEMTest_testField();
01644 }