src/MEDMEM/MEDMEM_Field.hxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 /*
  23  File Field.hxx
  24 */
  25
  26 #ifndef FIELD_HXX
  27 #define FIELD_HXX
  28
  29 #include "MEDMEM.hxx"
  30
  31 #include <vector>
  32 #include <map>
  33 #include <algorithm>
  34 #include <math.h>
  35 #include <cmath>
  36 #include <float.h>
  37
  38 #include "MEDMEM_Utilities.hxx"
  39 #include "MEDMEM_Exception.hxx"
  40 #include "MEDMEM_define.hxx"
  41 #include "MEDMEM_Support.hxx"
  42 #include "MEDMEM_Unit.hxx"
  43 #include "MEDMEM_nArray.hxx"
  44 #include "MEDMEM_GenDriver.hxx"
  45 #include "MEDMEM_ArrayInterface.hxx"
  46 #include "MEDMEM_SetInterlacingType.hxx"
  47 #include "MEDMEM_FieldForward.hxx"
  48 #include "MEDMEM_GaussLocalization.hxx"
  49
  50
  51
  52 namespace MEDMEM {
  53
  54   template<class T>
  55   struct MinMax {
  56   };
  57
  58   template<>
  59   struct MinMax<double> {
  60     static double getMin() { return DBL_MIN; }
  61     static double getMax() { return DBL_MAX; }
  62   };
  63
  64   template<>
  65   struct MinMax<int> {
  66     static int getMin() { return INT_MIN; }
  67     static int getMax() { return INT_MAX; }
  68   };
  69
  70   template < typename T > struct SET_VALUE_TYPE {
  71     static const MED_EN::med_type_champ _valueType = MED_EN::MED_UNDEFINED_TYPE;};
  72   template < > struct SET_VALUE_TYPE<double> {
  73     static const MED_EN::med_type_champ _valueType = MED_EN::MED_REEL64; };
  74   template < > struct SET_VALUE_TYPE<int> {
  75     static const MED_EN::med_type_champ _valueType = MED_EN::MED_INT32; };
  76
  77         /*!\defgroup FIELD_io Reading and writing files
  78
  79 Fields can be read or written to/from MED files.
  80
  81 \par Reading fields
  82
  83  For reading a field a typical use consists in :
  84 - reading the mesh associated on which the field lies
  85 - retrieve the support on which the field will be defined
  86 - read the field, specifying its time step and order number
  87
  88 As an example :
  89 \verbatim
  90 //reading mesh from file
  91 MESH mesh(MED_DRIVER, "file.med", "my_Mesh");
  92 //retrieving group in the mesh structure
  93 GROUP* group= mesh->getGroup("myGroup");
  94 //reading the field from the file
  95 FIELD<double> field(group,MED_DRIVER,"file.med","my_Field",1,1);
  96 \endverbatim
  97
  98 If the field is defined on all elements, one could have :
  99 \verbatim
 100 //reading mesh from file
 101 MESH mesh(MED_DRIVER, "file.med", "my_Mesh");
 102 //creating a support on all faces
 103 SUPPORT support (mesh,"mySupport",MED_FACE);
 104 //reading the field from the file
 105 FIELD<double> field(&support,MED_DRIVER,"file.med","my_FieldOnFaces",1,1);
 106 \endverbatim
 107
 108 It is also possible to read a field without specifying its support. In this case, the field constructor
 109 creates a support with no link to the initial mesh:
 110 \verbatim
 111 FIELD<double> field(MED_DRIVER, "file.med", "myField",1,1);
 112 SUPPORT* support= field->getSupport();
 113 \endverbatim
 114
 115 See also \ref FIELD_constructors
 116
 117 \par Writing fields
 118
 119 When it comes to write fields, it is necessary to use addDriver and then write.
 120 A typical use will be :
 121
 122 \verbatim
 123 mesh.addDriver(MED_DRIVER, "myResultFile.med", "myMesh");
 124 mesh.write();
 125 field.addDriver(MED_DRIVER, "myResultFile.med, "myField");
 126 field.write();
 127 \endverbatim
 128
 129 \defgroup FIELD_constructors
 130
 131 The different field constructors correspond to the two main
 132 ways a field is used :
 133 - either it is read from a file to be consulted,
 134 - or it can be created from scratch with a link to a support on which the values will be built.
 135
 136 \defgroup FIELD_algo Numerical operations on fields
 137 This section groups together the different operators that enable the user to
 138 treat the FIELD objects as high-level numerical arrays, giving operators for
 139 numerical treatment (overloading of basic operators, algorithms, etc...)
 140
 141 \defgroup FIELD_getset Basic Get/Set operations
 142
 143 This sections groups together the basic operations
 144 that describe access to all the elements constitutive of the description of the field :
 145 - name (compulsory),
 146 - time iteration number(compulsory),
 147 - inner loop iteration number(compulsory),
 148 - time(compulsory),
 149 - description(optional),
 150 - number of components(compulsory),
 151 - components names(optional),
 152 - components description(optional).
 153
 154 Some of these items are compulsory because they are essential to the field in order to define
 155 its structure or to be identified inside a MED file during the write process. The other ones
 156 are there for additional information and can be overlooked if not necessary.
 157
 158 When creating a field by reading a file, all the parameters are set according to the file
 159 data and can be consulted via the get methods. When creating a file from scratch, the
 160  name and number of components are set by the constructor, but the other items have to be
 161  set via the setXXX methods.
 162
 163 \defgroup FIELD_gauss Gauss points
 164
 165 In MED, it is possible to declare a Gauss model that
 166 describes the location of Gauss points in a reference cell.
 167 This Gauss model definition is contained in the
 168 \a GAUSS_LOCALIZATION class. A \a GAUSS_LOCALIZATION object
 169 is associated to a field and to a type.
 170
 171 It is not permitted to define a Gauss model in a polygonal
 172 or polyhedric element.
 173
 174 The Gauss model can be :
 175 - loaded from a MED file,
 176 - written to a MED file,
 177 - used to define a FIELD with multiple localizations per element.
 178
 179 \section gauss_constructors Constructing a Gauss Model
 180
 181 A Gauss model can be constructed with the following constructor :
 182 \param locName defines a name associated with the gauss model
 183 \param typeGeo names the type to which the Gauss model is assocaited
 184 \param nGauss defines the number of Gauss points
 185 \param cooRef defines an array giving the coordinates of the nodes of the reference element (dimension : spaceDimension * number of nodes for type \a typeGeo)
 186 \param cooGauss defines an array giving the coordinates of the nodes of the Gauss points (dimension : spaceDimension * \a nGauss
 187 )
 188 \param wg weights associated with each Gauss point (dimension : \a nGauss)
 189
 190 Example : in 2D, a Gauss model definition for a triangle
 191 would be written as :
 192
 193 \code
 194 string locname("gauss model");
 195 double cooRef[6] ={0.0, 0.0, 1.0, 0.0, 0.0, 1.0};
 196 double cooGauss[6]={0.2, 0.2, 0.8, 0.1, 0.1, 0.8};
 197 double wg[3]={0.3334, 0.3334, 0.3334};
 198 GAUSS_LOCALIZATION model(locname,
 199                          MED_EN::MED_TRIA3,
 200                          3,
 201                          cooRef,
 202                          cooGauss,
 203                          wg);
 204 \endcode
 205
 206 */
 207
 208
 209 /*!
 210
 211   This class contains all the informations related with a template class FIELD :
 212   - Components descriptions
 213   - Time step description
 214   - Location of the values (a SUPPORT class)
 215
 216 */
 217 class MEDMEM_EXPORT FIELD_    // GENERIC POINTER TO a template <class T, class INTERLACING_TAG> class FIELD
 218 {
 219 protected:
 220
 221   bool            _isRead ;
 222   bool            _isMinMax;
 223
 224   /*!
 225     \if developper
 226     Field name.
 227     \endif
 228   */
 229   string          _name ;
 230   /*!
 231     \if developper
 232     Field description.
 233     \endif
 234   */
 235   string          _description ;
 236   /*!
 237     \if developper
 238     Pointer to the support the field deals with.
 239     \endif
 240   */
 241   const SUPPORT * _support ;
 242
 243   /*!
 244     \if developper
 245     Number of field's components.
 246     \endif
 247   */
 248   int             _numberOfComponents ;
 249   /*!
 250     \if developper
 251     Number of field's values.
 252     doesn't take care of _numberOfComponents
 253     and number of Gauss points.
 254     \endif
 255   */
 256   int             _numberOfValues ;
 257
 258   /*!
 259     \if developper
 260     Array of size _numberOfComponents. \n
 261     (constant, scalar, vector, tensor)\n
 262     We could use an array of integer to store
 263     numbers of values: \n
 264     - 1 for scalar,\n
 265     - space dimension for vector,\n
 266     - space dimension square for tensor.\n
 267     So numbers of values per entities would be
 268     sum of _componentsTypes array.
 269
 270     Not implemented yet! All type are scalar !
 271     \endif
 272   */
 273   //int *           _componentsTypes ;
 274   vector<int>     _componentsTypes ;
 275   /*!
 276     \if developper
 277     Array of size _numberOfComponents
 278     storing components names if any.
 279     \endif
 280   */
 281   //string *      _componentsNames;
 282   vector<string>  _componentsNames;
 283   /*!
 284     \if developper
 285     Array of size _numberOfComponents
 286     storing components descriptions if any.
 287     \endif
 288   */
 289   //string *        _componentsDescriptions;
 290   vector<string>  _componentsDescriptions;
 291   /*!
 292     \if developper
 293     Array of size _numberOfComponents
 294     storing components units if any.
 295     \endif
 296   */
 297   //UNIT *          _componentsUnits;
 298   vector<UNIT>    _componentsUnits;
 299   /*!
 300     \if developper
 301     Array of size _numberOfComponents
 302     storing components units if any.
 303     \endif
 304   */
 305   //string *        _MEDComponentsUnits;
 306   vector<string>  _MEDComponentsUnits;
 307   /*!
 308     \if developper
 309     Iteration number of the field.
 310     \endif
 311   */
 312   int             _iterationNumber ;
 313   /*!
 314     \if developper
 315     Time of the field.
 316     \endif
 317   */
 318   double          _time;
 319   /*!
 320     \if developper
 321     Order number of the field.
 322     \endif
 323   */
 324   int             _orderNumber ;
 325   /*!
 326     \if developper
 327     At the initialization step of the field using the constructors; this attribute,
 328     the value type (integer or real) , is set automatically. There is a get method
 329     but not a set method for this attribute.
 330     \endif
 331   */
 332   MED_EN::med_type_champ _valueType ;
 333   /*!
 334     \if developper
 335     At the initialization step of the field using the constructors; this attribute,
 336     the interlacing type (full interlace or no interlace field value storage), is set
 337     automatically. There is a get method but not a set method for this attribute.
 338     \endif
 339   */
 340   MED_EN::medModeSwitch _interlacingType;
 341
 342   vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
 343   static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
 344   static void _deepCheckFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
 345   void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL) const  throw (MEDEXCEPTION);
 346   FIELD<double>* _getFieldSize() const;
 347
 348 public:
 349
 350   friend class MED_MED_RDONLY_DRIVER21;
 351   friend class MED_MED_WRONLY_DRIVER21;
 352   friend class MED_MED_RDWR_DRIVER21;
 353   friend class MED_MED_RDONLY_DRIVER22;
 354   friend class MED_MED_WRONLY_DRIVER22;
 355   friend class MED_MED_RDWR_DRIVER22;
 356   friend class VTK_MED_DRIVER;
 357
 358   /*!
 359     Constructor.
 360   */
 361   FIELD_ ();
 362   /*! \ifnot MEDMEM_ug
 363     Constructor.
 364 \endif
 365   */
 366   FIELD_(const SUPPORT * Support, const int NumberOfComponents);
 367   /*!  \ifnot MEDMEM_ug
 368     Copy constructor.
 369 \endif
 370   */
 371   FIELD_(const FIELD_ &m);
 372
 373   /*!
 374     Destructor.
 375   */
 376   virtual ~FIELD_();
 377
 378  FIELD_& operator=(const FIELD_ &m);
 379
 380   virtual  void     rmDriver(int index=0);
 381
 382   /*! \if MEDMEM_ug
 383     \addtogroup FIELD_io
 384     @{
 385     \endif
 386   */
 387
 388   /*! Creates a driver for reading/writing fields in a file.
 389     \param driverType specifies the file type (MED_DRIVER, VTK_DRIVER)
 390     \param fileName name of the output file
 391     \param driverFieldName name of the field
 392     \param access specifies whether the file is opened for read, write or both.
 393   */
 394
 395   virtual   int     addDriver(driverTypes driverType,
 396                               const string & fileName="Default File Name.med",
 397                               const string & driverFieldName="Default Field Nam",
 398                               MED_EN::med_mode_acces access=MED_EN::RDWR) ;
 399
 400   virtual  int      addDriver( GENDRIVER & driver);
 401   virtual  void     read (const GENDRIVER &);
 402   virtual  void     read(int index=0);
 403   virtual  void     openAppend( void );
 404   virtual  void     write(const GENDRIVER &);
 405
 406   /*! Triggers the writing of the field with respect to the driver handle
 407     \a index given by \a addDriver(...) method. */
 408   virtual  void     write(int index=0, const string & driverName="");
 409   /*!\if MEDMEM_ug @} \endif */
 410
 411   virtual  void     writeAppend(const GENDRIVER &);
 412   virtual  void     writeAppend(int index=0, const string & driverName="");
 413
 414   inline void     setName(const string Name);
 415   inline string   getName() const;
 416   inline void     setDescription(const string Description);
 417   inline string   getDescription() const;
 418   inline const SUPPORT * getSupport() const;
 419   inline void     setSupport(const SUPPORT * support);
 420   /*inline*/ void     setNumberOfComponents(const int NumberOfComponents);
 421   inline int      getNumberOfComponents() const;
 422   inline void     setNumberOfValues(const int NumberOfValues);
 423   inline int      getNumberOfValues() const;
 424   //    inline void     setComponentType(int *ComponentType);
 425   //    inline int *    getComponentType() const;
 426   //    inline int      getComponentTypeI(int i) const;
 427   inline void     setComponentsNames(const string * ComponentsNames);
 428   inline void     setComponentName(int i, const string ComponentName);
 429   inline const string * getComponentsNames() const;
 430   inline string   getComponentName(int i) const;
 431   inline void     setComponentsDescriptions(const string * ComponentsDescriptions);
 432   inline void     setComponentDescription(int i, const string ComponentDescription);
 433   inline const string * getComponentsDescriptions() const;
 434   inline string   getComponentDescription(int i) const;
 435
 436   // provisoire : en attendant de regler le probleme des unites !
 437   inline void     setComponentsUnits(const UNIT * ComponentsUnits);
 438   inline const UNIT *   getComponentsUnits() const;
 439   inline const UNIT *   getComponentUnit(int i) const;
 440   inline void     setMEDComponentsUnits(const string * MEDComponentsUnits);
 441   inline void     setMEDComponentUnit(int i, const string MEDComponentUnit);
 442   inline const string * getMEDComponentsUnits() const;
 443   inline string   getMEDComponentUnit(int i) const;
 444
 445   inline void     setIterationNumber(int IterationNumber);
 446   inline int      getIterationNumber() const;
 447   inline void     setTime(double Time);
 448   inline double   getTime() const;
 449   inline void     setOrderNumber(int OrderNumber);
 450   inline int      getOrderNumber() const;
 451
 452   inline MED_EN::med_type_champ getValueType () const;
 453   inline MED_EN::medModeSwitch  getInterlacingType() const;
 454   virtual inline bool getGaussPresence() const throw (MEDEXCEPTION);
 455 protected:
 456   void copyGlobalInfo(const FIELD_& m);
 457 };
 458
 459 // ---------------------------------
 460 // Implemented Methods : constructor
 461 // ---------------------------------
 462
 463 // -----------------
 464 // Methodes Inline
 465 // -----------------
 466 /*! \if MEDMEM_ug
 467 \addtogroup FIELD_getset
 468 @{
 469 \endif
 470 */
 471 /*!
 472   Sets FIELD name. The length should not exceed MED_TAILLE_NOM
 473 as defined in Med (i.e. 32 characters).
 474 */
 475 inline void FIELD_::setName(const string Name)
 476 {
 477   _name=Name;
 478 }
 479 /*!
 480   Gets FIELD name.
 481 */
 482 inline string FIELD_::getName() const
 483 {
 484   return _name;
 485 }
 486 /*!
 487   Sets FIELD description. The length should not exceed MED_TAILLE_DESC as defined in Med (i.e. 200 characters).
 488 */
 489 inline void FIELD_::setDescription(const string Description)
 490 {
 491   _description=Description;
 492 }
 493 /*!
 494   Gets FIELD description.
 495 */
 496 inline string FIELD_::getDescription() const
 497 {
 498   return _description;
 499 }
 500 /*!
 501   Sets FIELD number of components.
 502 */
 503 inline void FIELD_::setNumberOfComponents(const int NumberOfComponents)
 504 {
 505   _numberOfComponents=NumberOfComponents;
 506   _componentsTypes.resize(_numberOfComponents);
 507   _componentsNames.resize(_numberOfComponents);
 508   _componentsDescriptions.resize(_numberOfComponents);
 509   _componentsUnits.resize(_numberOfComponents);
 510   _MEDComponentsUnits.resize(_numberOfComponents);
 511 }
 512 /*!
 513   Gets FIELD number of components.
 514 */
 515 inline int FIELD_::getNumberOfComponents() const
 516 {
 517   return _numberOfComponents ;
 518 }
 519 /*!
 520   Sets FIELD number of values.
 521
 522   It must be the same than in the associated SUPPORT object.
 523 */
 524 inline void FIELD_::setNumberOfValues(const int NumberOfValues)
 525 {
 526   _numberOfValues=NumberOfValues;
 527 }
 528 /*!
 529   Gets FIELD number of value.
 530 */
 531 inline int FIELD_::getNumberOfValues() const
 532 {
 533   return _numberOfValues ;
 534 }
 535
 536 //  inline void FIELD_::setComponentType(int *ComponentType)
 537 //  {
 538 //    _componentsTypes=ComponentType ;
 539 //  }
 540 //  inline int * FIELD_::getComponentType() const
 541 //  {
 542 //    return _componentsTypes ;
 543 //  }
 544 //  inline int FIELD_::getComponentTypeI(int i) const
 545 //  {
 546 //    return _componentsTypes[i-1] ;
 547 //  }
 548
 549 /*!
 550   Sets FIELD components names.
 551
 552   Duplicates the ComponentsNames string array to put components names in
 553   FIELD. ComponentsNames size must be equal to number of components.
 554 */
 555 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
 556 {
 557   //if (NULL == _componentsNames)
 558   //  _componentsNames = new string[_numberOfComponents] ;
 559   _componentsNames.resize(_numberOfComponents);
 560   for (int i=0; i<_numberOfComponents; i++)
 561     _componentsNames[i]=ComponentsNames[i] ;
 562 }
 563 /*! \ifnot MEDMEM_ug
 564   Sets FIELD i^th component name.
 565
 566   i must be >=1 and <= number of components.
 567 \endif
 568 */
 569 inline void FIELD_::setComponentName(int i, const string ComponentName)
 570 {
 571   const char * LOC = " FIELD_::setComponentName() : ";
 572   BEGIN_OF_MED(LOC);
 573   if( i<1 || i>_numberOfComponents )
 574     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 575
 576   _componentsNames[i-1]=ComponentName ;
 577 }
 578 /*!
 579   Gets a reference to the string array which contain the components names.
 580
 581   This Array size is equal to number of components
 582 */
 583 inline const string * FIELD_::getComponentsNames() const
 584 {
 585   return &(_componentsNames[0]) ;
 586 }
 587 /*!\ifnot MEDMEM_ug
 588   Gets the name of the i^th component.
 589 \endif
 590 */
 591 inline string FIELD_::getComponentName(int i) const
 592 {
 593   const char * LOC = " FIELD_::getComponentName() : ";
 594   BEGIN_OF_MED(LOC);
 595   if( i<1 || i>_numberOfComponents )
 596     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 597
 598   return _componentsNames[i-1] ;
 599 }
 600 /*!
 601   Sets FIELD components descriptions.
 602
 603   Duplicates the ComponentsDescriptions string array to put components
 604   descriptions in FIELD.
 605   ComponentsDescriptions size must be equal to number of components.
 606 */
 607 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
 608 {
 609   //if (NULL == _componentsDescriptions)
 610   //  _componentsDescriptions = new string[_numberOfComponents] ;
 611   _componentsDescriptions.resize(_numberOfComponents);
 612   for (int i=0; i<_numberOfComponents; i++)
 613     _componentsDescriptions[i]=ComponentsDescriptions[i] ;
 614 }
 615 /*!\ifnot MEDMEM_ug
 616   Sets FIELD i^th component description.
 617
 618   i must be >=1 and <= number of components.
 619 \endif
 620 */
 621 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
 622 {
 623   const char * LOC = " FIELD_::setComponentDescription() : ";
 624   BEGIN_OF_MED(LOC);
 625   if( i<1 || i>_numberOfComponents )
 626     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 627
 628   _componentsDescriptions[i-1]=ComponentDescription ;
 629 }
 630 /*!
 631   Gets a reference to the string array which contain the components descriptions.
 632
 633   This Array size is equal to number of components
 634 */
 635 inline const string * FIELD_::getComponentsDescriptions() const
 636 {
 637   return &(_componentsDescriptions[0]);
 638 }
 639 /*!\ifnot MEDMEM_ug
 640   Gets the description of the i^th component.
 641 \endif
 642 */
 643 inline string FIELD_::getComponentDescription(int i) const
 644 {
 645   const char * LOC = " FIELD_::setComponentDescription() : ";
 646   BEGIN_OF_MED(LOC);
 647   if( i<1 || i>_numberOfComponents )
 648     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 649
 650   return _componentsDescriptions[i-1];
 651 }
 652
 653 /*!\ifnot MEDMEM_ug
 654   \todo
 655   Sets FIELD components UNIT.
 656
 657   Duplicates the ComponentsUnits UNIT array to put components
 658   units in FIELD.
 659   ComponentsUnits size must be equal to number of components.
 660 \endif
 661 */
 662 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
 663 {
 664   //if (NULL == _componentsUnits)
 665   //  _componentsUnits = new UNIT[_numberOfComponents] ;
 666   _componentsUnits.resize(_numberOfComponents);
 667   for (int i=0; i<_numberOfComponents; i++)
 668     _componentsUnits[i]=ComponentsUnits[i] ;
 669 }
 670 /*!\ifnot MEDMEM_ug
 671   Gets a reference to the UNIT array which contain the components units.
 672
 673   This array size is equal to number of components
 674 \endif
 675 */
 676 inline const UNIT * FIELD_::getComponentsUnits() const
 677 {
 678   return &(_componentsUnits[0]);
 679 }
 680 /*!\ifnot MEDMEM_ug
 681   Gets the UNIT of the i^th component.
 682 \endif
 683 */
 684 inline const UNIT * FIELD_::getComponentUnit(int i) const
 685 {
 686   const char * LOC = " FIELD_::getComponentUnit() : ";
 687   BEGIN_OF_MED(LOC);
 688   if( i<1 || i>_numberOfComponents )
 689     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 690
 691   return &_componentsUnits[i-1] ;
 692 }
 693 /*!
 694   Sets FIELD components unit.
 695
 696   Duplicates the MEDComponentsUnits string array to put components
 697   units in FIELD.
 698   MEDComponentsUnits size must be equal to number of components.
 699
 700 */
 701 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
 702 {
 703   //if (NULL == _MEDComponentsUnits)
 704   //  _MEDComponentsUnits = new string[_numberOfComponents] ;
 705   _MEDComponentsUnits.resize(_numberOfComponents);
 706   for (int i=0; i<_numberOfComponents; i++)
 707     _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
 708 }
 709 /*!\ifnot MEDMEM_ug
 710   Sets FIELD i^th component unit.
 711
 712   i must be >=1 and <= number of components.
 713 \endif
 714 */
 715 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
 716 {
 717   const char * LOC = " FIELD_::setMEDComponentUnit() : ";
 718   BEGIN_OF_MED(LOC);
 719   if( i<1 || i>_numberOfComponents )
 720     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 721
 722   _MEDComponentsUnits[i-1]=MEDComponentUnit ;
 723 }
 724 /*!
 725   Gets a reference to the string array which contain the components units.
 726
 727   This array size is equal to number of components
 728 */
 729 inline const string * FIELD_::getMEDComponentsUnits() const
 730 {
 731   return &(_MEDComponentsUnits[0]);
 732 }
 733 /*! \ifnot MEDMEM_ug
 734   Gets the string for unit of the i^th component.
 735 \endif
 736 */
 737 inline string FIELD_::getMEDComponentUnit(int i) const
 738 {
 739   const char * LOC = " FIELD_::getMEDComponentUnit() : ";
 740   BEGIN_OF_MED(LOC);
 741   if( i<1 || i>_numberOfComponents )
 742     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 743
 744   return _MEDComponentsUnits[i-1] ;
 745 }
 746 /*!
 747   Sets the iteration number where FIELD has been calculated.
 748 */
 749 inline void FIELD_::setIterationNumber(int IterationNumber)
 750 {
 751   _iterationNumber=IterationNumber;
 752 }
 753 /*!
 754   Gets the iteration number where FIELD has been calculated.
 755 */
 756 inline int FIELD_::getIterationNumber() const
 757 {
 758   return _iterationNumber ;
 759 }
 760 /*!
 761   Sets the time when FIELD has been calculated.
 762 */
 763 inline void FIELD_::setTime(double Time)
 764 {
 765   _time=Time ;
 766 }
 767 /*!
 768   Gets the time when FIELD has been calculated.
 769 */
 770 inline double FIELD_::getTime() const
 771 {
 772   return _time ;
 773 }
 774 /*!
 775   Sets the order number where FIELD has been calculated.
 776
 777   It corresponds to internal iteration during one time step.
 778 */
 779 inline void FIELD_::setOrderNumber(int OrderNumber)
 780 {
 781   _orderNumber=OrderNumber ;
 782 }
 783 /*!
 784   Gets the order number where FIELD has been calculated.
 785 */
 786 inline int FIELD_::getOrderNumber() const
 787 {
 788   return _orderNumber ;
 789 }
 790 /*!
 791   Gets a reference to the SUPPORT object associated to FIELD.
 792 */
 793 inline  const SUPPORT * FIELD_::getSupport() const
 794 {
 795   return _support ;
 796 }
 797 /*!
 798   Sets the reference to the SUPPORT object associated to FIELD.
 799
 800   Reference is not duplicate, so it must not be deleted.
 801 */
 802 inline void FIELD_::setSupport(const SUPPORT * support)
 803 {
 804   //A.G. Addings for RC
 805   if(_support)
 806     _support->removeReference();
 807   _support = support ;
 808   if(_support)
 809     _support->addReference();
 810 }
 811 /*!
 812   Gets the FIELD med value type (MED_INT32 or MED_REEL64).
 813 */
 814 inline MED_EN::med_type_champ FIELD_::getValueType () const
 815 {
 816   return _valueType ;
 817 }
 818
 819 /*!
 820   Gets the FIELD med interlacing type (MED_FULL_INTERLACE or MED_NO_INTERLACE).
 821 */
 822   inline MED_EN::medModeSwitch FIELD_::getInterlacingType () const
 823 {
 824   return _interlacingType ;
 825 }
 826  /*!\if MEDMEM_ug @} \endif*/
 827
 828 /*!\if MEDMEM_ug
 829 \addtogroup FIELD_gauss
 830 @{
 831 \endif */
 832
 833 /*!
 834  Determines whether the field stores several Gauss points per element.
 835 */
 836   inline bool  FIELD_::getGaussPresence() const throw (MEDEXCEPTION)
 837 {
 838   const char * LOC = "FIELD_::getGaussPresence() : ";
 839   throw MEDEXCEPTION(STRING(LOC) << " This FIELD_ doesn't rely on a FIELD<T>" );
 840 }
 841
 842   /*!\if MEDMEM_ug @} \endif*/
 843
 844 } //End namespace MEDMEM
 845
 846 /////////////////////////
 847 // END OF CLASS FIELD_ //
 848 /////////////////////////
 849
 850 /*!
 851
 852   This template class contains informations related with a FIELD :
 853   - Values of the field, their type (real or integer), the storage mode (full interlace,
 854     no interlace or no interlace by type).
 855
 856 */
 857
 858
 859 namespace MEDMEM {
 860
 861   template<class T2> class MED_FIELD_RDONLY_DRIVER21;
 862   template<class T2> class MED_FIELD_WRONLY_DRIVER21;
 863   template<class T2> class MED_FIELD_RDONLY_DRIVER22;
 864   template<class T2> class MED_FIELD_WRONLY_DRIVER22;
 865   template<class T2> class VTK_FIELD_DRIVER;
 866
 867
 868   template <class T,
 869             class INTERLACING_TAG
 870             > class FIELD : public FIELD_
 871 {
 872 protected:
 873
 874   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array   ArrayNoGauss;
 875   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array     ArrayGauss;
 876   typedef typename MEDMEM_ArrayInterface<T,NoInterlace,NoGauss>::Array       ArrayNo;
 877   typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array     ArrayFull;
 878   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,NoGauss>::Array ArrayNoByType;
 879   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,Gauss>::Array   ArrayNoByTypeGauss;
 880   typedef MEDMEM_Array_ Array;
 881   typedef T ElementType;
 882   typedef INTERLACING_TAG InterlacingTag;
 883   typedef map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> locMap;
 884
 885   // array of value of type T
 886   Array *_value ;
 887
 888   // MESH, to be used for field reading from a file (if desired to link
 889   // to existing support instead of new support creation for the field)
 890   MESH* _mesh;
 891
 892   // extrema values
 893   T _vmin;
 894   T _vmax;
 895
 896   map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> _gaussModel; //A changer quand les drivers seront template de l'entrelacement
 897
 898   static T _scalarForPow;
 899   static T pow(T x);
 900
 901 private:
 902   void _operation(const FIELD& m,const FIELD& n, const char* Op);
 903   void _operationInitialize(const FIELD& m,const FIELD& n, const char* Op);
 904   void _add_in_place(const FIELD& m,const FIELD& n);
 905   void _sub_in_place(const FIELD& m,const FIELD& n);
 906   void _mul_in_place(const FIELD& m,const FIELD& n);
 907   void _div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION);
 908   //setValueType() ;
 909
 910 public:
 911   FIELD();
 912   FIELD(const FIELD &m);
 913   FIELD(const SUPPORT * Support, const int NumberOfComponents) throw (MEDEXCEPTION);
 914   FIELD(driverTypes driverType,
 915         const string & fileName, const string & fieldDriverName,
 916         const int iterationNumber=-1, const int orderNumber=-1,
 917         MESH* mesh = 0)
 918     throw (MEDEXCEPTION);
 919   FIELD(const SUPPORT * Support, driverTypes driverType,
 920         const string & fileName="", const string & fieldName="",
 921         const int iterationNumber = -1, const int orderNumber = -1)
 922     throw (MEDEXCEPTION);
 923   ~FIELD();
 924
 925   FIELD & operator=(const FIELD &m);
 926         FIELD & operator=(T value);
 927   const FIELD operator+(const FIELD& m) const;
 928   const FIELD operator-(const FIELD& m) const;
 929   const FIELD operator*(const FIELD& m) const;
 930   const FIELD operator/(const FIELD& m) const;
 931   const FIELD operator-() const;
 932   FIELD& operator+=(const FIELD& m);
 933   FIELD& operator-=(const FIELD& m);
 934   FIELD& operator*=(const FIELD& m);
 935   FIELD& operator/=(const FIELD& m);
 936
 937   void          applyLin(T a, T b, int icomp);
 938   static FIELD* add(const FIELD& m, const FIELD& n);
 939   static FIELD* addDeep(const FIELD& m, const FIELD& n);
 940   static FIELD* sub(const FIELD& m, const FIELD& n);
 941   static FIELD* subDeep(const FIELD& m, const FIELD& n);
 942   static FIELD* mul(const FIELD& m, const FIELD& n);
 943   static FIELD* mulDeep(const FIELD& m, const FIELD& n);
 944   static FIELD* div(const FIELD& m, const FIELD& n);
 945   static FIELD* divDeep(const FIELD& m, const FIELD& n);
 946   double normMax() const throw (MEDEXCEPTION);
 947
 948   //------- TDG and BS addings
 949
 950   void getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION);
 951   vector<int> getHistogram(int &nbint) throw (MEDEXCEPTION);
 952   FIELD<double>* buildGradient() const throw (MEDEXCEPTION);
 953   FIELD<double>* buildNorm2Field() const throw (MEDEXCEPTION);
 954
 955   //-------------------
 956
 957   double norm2() const throw (MEDEXCEPTION);
 958   void   applyLin(T a, T b);
 959   template <T T_function(T)> void applyFunc();
 960   void applyPow(T scalar);
 961   static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
 962   double normL2(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 963   double normL2(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 964   double normL1(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 965   double normL1(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 966   FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
 967
 968   friend class MED_FIELD_RDONLY_DRIVER21<T>;
 969   friend class MED_FIELD_WRONLY_DRIVER21<T>;
 970   friend class MED_FIELD_RDONLY_DRIVER22<T>;
 971   friend class MED_FIELD_WRONLY_DRIVER22<T>;
 972   friend class VTK_FIELD_DRIVER<T>;
 973   //friend class MED_FIELD_RDWR_DRIVER  <T>;
 974
 975   void init ();
 976   void rmDriver(int index=0);
 977   int  addDriver(driverTypes driverType,
 978                  const string & fileName="Default File Name.med",
 979                  const string & driverFieldName="Default Field Name",
 980                  MED_EN::med_mode_acces access=MED_EN::RDWR) ;
 981
 982   int  addDriver(GENDRIVER & driver);
 983
 984   void allocValue(const int NumberOfComponents);
 985   void allocValue(const int NumberOfComponents, const int LengthValue);
 986
 987   void deallocValue();
 988
 989   inline void read(int index=0);
 990   inline void read(const GENDRIVER & genDriver);
 991   inline void write(int index=0, const string & driverName = "");
 992   inline void write(const GENDRIVER &);
 993
 994   inline void writeAppend(int index=0, const string & driverName = "");
 995   inline void writeAppend(const GENDRIVER &);
 996
 997   inline MEDMEM_Array_  * getArray()        const throw (MEDEXCEPTION);
 998   inline ArrayGauss     * getArrayGauss()   const throw (MEDEXCEPTION);
 999   inline ArrayNoGauss   * getArrayNoGauss() const throw (MEDEXCEPTION);
1000   inline bool            getGaussPresence() const throw (MEDEXCEPTION);
1001
1002   inline int          getValueLength() const throw (MEDEXCEPTION);
1003
1004   /*! \if MEDMEM_ug
1005 \addtogroup FIELD_value
1006 @{
1007 \endif  */
1008   /*! Returns a pointer to the value array.*/
1009   inline const T*     getValue()       const throw (MEDEXCEPTION);
1010   inline const T*     getRow(int i)    const throw (MEDEXCEPTION);
1011   inline const T*     getColumn(int j) const throw (MEDEXCEPTION);
1012 /*!
1013   Returns the value of \f$ i^{th} \f$ element and \f$ j^{th}\f$ component.
1014   This method only works with fields having no particular Gauss point
1015 definition (i.e., fields having one value per element).
1016  This method makes the retrieval of the value independent from the
1017   interlacing pattern, but it is slower than the complete retrieval
1018   obtained by the \b getValue() method.
1019 */
1020
1021   inline T            getValueIJ(int i,int j) const throw (MEDEXCEPTION);
1022
1023 /*!
1024   Returns the \f$ j^{th}\f$  component of \f$ k^{th}\f$  Gauss points of \f$ i^{th}\f$  value.
1025   This method is compatible with elements having more than one Gauss point.
1026   This method makes the retrieval of the value independent from the
1027   interlacing pattern, but it is slower than the complete retrieval
1028   obtained by the \b getValue() method.
1029 */
1030   inline T            getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION);
1031
1032   inline int          getValueByTypeLength(int t)                const throw (MEDEXCEPTION);
1033   inline const T*     getValueByType(int t)                      const throw (MEDEXCEPTION);
1034   inline T            getValueIJByType(int i,int j,int t)        const throw (MEDEXCEPTION);
1035   inline T            getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION);
1036
1037   /*!
1038                 The following example describes the creation of a FIELD.
1039
1040                 \example FIELDcreate.cxx
1041
1042  \if MEDMEM_ug @} \endif */
1043
1044   bool                getValueOnElement(int eltIdInSup,T* retValues) const throw (MEDEXCEPTION);
1045
1046   const int   getNumberOfGeometricTypes() const throw (MEDEXCEPTION);
1047   const GAUSS_LOCALIZATION<INTERLACING_TAG> & getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1048   const GAUSS_LOCALIZATION<INTERLACING_TAG> * getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1049   const GAUSS_LOCALIZATION_* getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1050   void setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG> & gaussloc);
1051   void setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc);
1052   const int * getNumberOfGaussPoints() const throw (MEDEXCEPTION);
1053   const int   getNumberOfGaussPoints( MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1054   const int   getNbGaussI(int i)          const throw (MEDEXCEPTION);
1055   const int * getNumberOfElements()       const throw (MEDEXCEPTION);
1056   const MED_EN::medGeometryElement  * getGeometricTypes()  const throw (MEDEXCEPTION);
1057   bool        isOnAllElements()           const throw (MEDEXCEPTION);
1058
1059   inline void setArray(MEDMEM_Array_ *value) throw (MEDEXCEPTION);
1060
1061   /*! \if MEDMEM_ug
1062  \addtogroup FIELD_value
1063 @{
1064 \endif
1065   */
1066 /*!
1067 This method makes it possible to have the field pointing to
1068 an existing value array. The ordering of the elements in the value array must
1069 conform to the MEDMEM ordering (I,K,J) : the outer loop is on the elements,
1070 the intermediate loop is on the Gauss points, the inner loop is on
1071 the components.
1072 */
1073   inline void setValue( T* value) throw (MEDEXCEPTION);
1074   inline void setRow( int i, T* value) throw (MEDEXCEPTION);
1075   inline void setColumn( int i, T* value) throw (MEDEXCEPTION);
1076 /*!
1077   Sets the value of \f$ i^{th} \f$ element and \f$ j^{th}\f$ component with \a value.
1078 */
1079   inline void setValueIJ(int i, int j, T value) throw (MEDEXCEPTION);
1080   /*! \if MEDMEM_ug @} \endif */
1081   inline void setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION);
1082   inline void setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION);
1083   inline void setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION);
1084
1085   /*!
1086     This fonction feeds the FIELD<double> private attributs _value with the
1087     volume of each cells belonging to the argument Support. The field has to be
1088     initialised via the constructor FIELD<double>(const SUPPORT * , const int )
1089     with Support as SUPPORT argument, 1 has the number of components, and Support
1090     has be a SUPPORT on 3D cells. This initialisation could be done by the empty
1091     constructor followed by a setSupport and setNumberOfComponents call.
1092    */
1093   //void getVolume() const throw (MEDEXCEPTION) ;
1094   /*!
1095     This fonction feeds the FIELD<double> private attributs _value with the
1096     area of each cells (or faces) belonging to the attribut _support. The field
1097     has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
1098     const int ) with 1 has the number of components, and _support has be a
1099     SUPPORT on 2D cells or 3D faces. This initialisation could be done by the
1100     empty constructor followed by a setSupport and setNumberOfComponents call.
1101    */
1102   //void getArea() const throw (MEDEXCEPTION) ;
1103   /*!
1104     This fonction feeds the FIELD<double> private attributs _value with the
1105     length of each segments belonging to the attribut _support. The field has
1106     to be initialised via the constructor FIELD<double>(const SUPPORT * ,
1107     const int ) with 1 has the number of components, and _support has be a
1108     SUPPORT on 3D edges or 2D faces. This initialisation could be done by the
1109     empty constructor followed by a setSupport and setNumberOfComponents call.
1110    */
1111   //void getLength() const throw (MEDEXCEPTION) ;
1112   /*!
1113     This fonction feeds the FIELD<double> private attributs _value with the
1114     normal vector of each faces belonging to the attribut _support. The field
1115     has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
1116     const int ) with the space dimension has the number of components, and
1117     _support has be a SUPPORT on 3D or 2D faces. This initialisation could be done
1118     by the empty constructor followed by a setSupport and setNumberOfComponents
1119     call.
1120    */
1121   //void getNormal() const throw (MEDEXCEPTION) ;
1122   /*!
1123     This fonction feeds the FIELD<double> private attributs _value with the
1124     barycenter of each faces or cells or edges belonging to the attribut _support.
1125     The field has to be initialised via the constructor
1126     FIELD<double>(const SUPPORT * ,const int ) with the space dimension has the
1127     number of components, and _support has be a SUPPORT on 3D cells or 2D faces.
1128     This initialisation could be done by the empty constructor followed by a
1129     setSupport and setNumberOfComponents call.
1130    */
1131   //void getBarycenter() const throw (MEDEXCEPTION) ;
1132
1133   typedef void (*myFuncType)(const double *,T*);
1134   void fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION);
1135   typedef void (*myFuncType2)(const T *,T*);
1136   FIELD<T,INTERLACING_TAG> *execFunc(int nbOfComponents, myFuncType2 f) throw (MEDEXCEPTION);
1137 };
1138 }
1139
1140 #include "MEDMEM_DriverFactory.hxx"
1141
1142 namespace MEDMEM {
1143
1144 template <class T,class INTERLACING_TAG> T FIELD<T, INTERLACING_TAG>::_scalarForPow=1;
1145
1146 // --------------------
1147 // Implemented Methods
1148 // --------------------
1149
1150 /*!
1151   Constructor with no parameter, most of the attribut members are set to NULL.
1152 */
1153 template <class T, class INTERLACING_TAG>
1154 FIELD<T, INTERLACING_TAG>::FIELD():FIELD_()
1155 {
1156   MESSAGE_MED("Constructeur FIELD sans parametre");
1157
1158   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
1159   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE);
1160   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
1161
1162   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
1163   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE);
1164   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
1165
1166   _value = ( ArrayNoGauss * ) NULL;
1167 }
1168
1169         /*!
1170 \addtogroup FIELD_constructors FIELD<T> constructors
1171 @{
1172         */
1173
1174 /*!
1175   Constructor that allocates the value array with the dimensions provided by
1176 \a NumberOfComponents and the dimension of \a Support. The value array is
1177  allocated but not initialized.
1178 This constructor does not allow the creation of fields with Gauss points.
1179 \param Support support on which the field lies
1180 \param NumberOfComponents number of components of the variable stored. For instance,
1181 it will be 3 for a (vx,vy,vz) vector.
1182
1183 \code
1184 FIELD<double> field (support, 3);
1185 int nbelem = support->getNumberOfElements(MED_ALL_ELEMENTS);
1186 for (int i=1; i<=nbelem; i++)
1187    for (j=1; j<=3;j++)
1188        field->setValueIJ(i,j,0.0);
1189 \endcode
1190 */
1191 template <class T, class INTERLACING_TAG>
1192 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
1193                                  const int NumberOfComponents) throw (MEDEXCEPTION) :
1194   FIELD_(Support, NumberOfComponents),_value(NULL)
1195 {
1196   const char* LOC = "FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)";
1197   BEGIN_OF_MED(LOC);
1198   SCRUTE_MED(this);
1199
1200   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
1201   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
1202   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
1203
1204   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
1205   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
1206   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
1207
1208   try {
1209     // becarefull about the numbre of gauss point
1210     _numberOfValues = Support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
1211   }
1212 #if defined(_DEBUG_) || defined(_DEBUG)
1213   catch (MEDEXCEPTION &ex) {
1214 #else
1215   catch (MEDEXCEPTION ) {
1216 #endif
1217     MESSAGE_MED("No value defined ! ("<<ex.what()<<")");
1218   }
1219   MESSAGE_MED("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
1220   if (0<_numberOfValues) {
1221     _value = new ArrayNoGauss (_numberOfComponents,_numberOfValues);
1222     _isRead = true ;
1223   }
1224
1225   END_OF_MED(LOC);
1226 }
1227         /*!
1228 @}
1229         */
1230 /*!
1231   \if developper
1232   \endif
1233 */
1234 template <class T, class INTERLACING_TAG> void FIELD<T, INTERLACING_TAG>::init ()
1235 {
1236 }
1237
1238 /*!
1239   Copy constructor.
1240 */
1241 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::FIELD(const FIELD & m):
1242   FIELD_((FIELD_) m)
1243 {
1244   MESSAGE_MED("Constructeur FIELD de recopie");
1245
1246   // RECOPIE PROFONDE <> de l'operateur= Rmq from EF
1247   if (m._value != NULL)
1248     {
1249       if ( m.getGaussPresence() )
1250         _value = new ArrayGauss( *(static_cast< ArrayGauss * > (m._value) ) ,false);
1251       else
1252         _value = new ArrayNoGauss( *(static_cast< ArrayNoGauss * > (m._value)) ,false);
1253     }
1254   else
1255     _value = (ArrayNoGauss *) NULL;
1256   locMap::const_iterator it;
1257   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
1258     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
1259       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
1260                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
1261                                               );
1262
1263   _valueType       = m._valueType;
1264   _interlacingType = m._interlacingType;
1265   //drivers = m._drivers;
1266 }
1267
1268 /*!
1269   \if developper
1270   Not implemented.
1271   \endif
1272 */
1273 template <class T, class INTERLACING_TAG>
1274 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(const FIELD &m)
1275 {
1276   MESSAGE_MED("Appel de FIELD<T>::operator=") ;
1277   if ( this == &m) return *this;
1278
1279   // copy values array
1280   FIELD_::operator=(m);  // Driver are ignored & ?copie su pointeur de Support?
1281
1282   _value           = m._value; //PROBLEME RECOPIE DES POINTEURS PAS COHERENT AVEC LE
1283                                //CONSTRUCTEUR PAR RECOPIE
1284                                //CF :Commentaire dans MEDMEM_Array
1285   locMap::const_iterator it;
1286   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
1287     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
1288       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
1289                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
1290                                               );
1291
1292   _valueType       = m._valueType;
1293   _interlacingType = m._interlacingType;
1294
1295   return *this;
1296 }
1297
1298 /*!
1299         Initializes all the field values to \a value
1300 */
1301 template <class T, class INTERLACING_TAG>
1302 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(T value)
1303 {
1304   MESSAGE_MED("Appel de FIELD<T>::operator= T") ;
1305         int size=getNumberOfComponents()*getNumberOfValues();
1306         T* ptr= const_cast<T*>( getValue());
1307         for (int i=0; i< size; i++)
1308                 {*ptr++=value;}
1309
1310   return *this;
1311 }
1312
1313 /*!\addtogroup FIELD_algo
1314 @{
1315 */
1316
1317 /*!
1318      Overload addition operator.
1319      This operation is authorized only for compatible fields that have the same support.
1320      The compatibility checking includes equality tests of the folowing data members:\n
1321      - _support
1322      - _numberOfComponents
1323      - _numberOfValues
1324      - _componentsTypes
1325      - _MEDComponentsUnits.
1326
1327      The data members of the returned field are initialized, based on the first field, except for the name,
1328      which is the combination of the two field's names and the operator.
1329      Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> \n
1330      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1331      When using python, this operator calls the copy constructor in any case.
1332      The user has to be aware that when using operator + in associatives expressions like
1333      <tt> a = b + c + d +e; </tt> \n
1334      no optimisation is performed : the evaluation of last expression requires the construction of
1335      3 temporary fields.
1336 */
1337 template <class T, class INTERLACING_TAG>
1338 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator+(const FIELD & m) const
1339 {
1340   const char* LOC = "FIELD<T>::operator+(const FIELD & m)";
1341   BEGIN_OF_MED(LOC);
1342     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1343
1344     // Creation of the result - memory is allocated by FIELD constructor
1345     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1346     //result._operation(*this,m,mode,"+"); // perform Atribute's initialization & addition
1347     result._operationInitialize(*this,m,"+"); // perform Atribute's initialization
1348     result._add_in_place(*this,m); // perform addition
1349
1350   END_OF_MED(LOC);
1351     return result;
1352 }
1353
1354 /*!  Overloaded Operator +=
1355  *   Operations are directly performed in the first field's array.
1356  *   This operation is authorized only for compatible fields that have the same support.
1357  */
1358 template <class T, class INTERLACING_TAG>
1359 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator+=(const FIELD & m)
1360 {
1361   const char* LOC = "FIELD<T>::operator+=(const FIELD & m)";
1362   BEGIN_OF_MED(LOC);
1363     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1364
1365     const T* value1=m.getValue(); // get pointers to the values we are adding
1366
1367     // get a non const pointer to the inside array of values and perform operation
1368     T * value=const_cast<T *> (getValue());
1369     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1370     const T* endV=value+size; // pointer to the end of value
1371     for(;value!=endV; value1++,value++)
1372         *value += *value1;
1373   END_OF_MED(LOC);
1374     return *this;
1375 }
1376
1377
1378 /*! Addition of fields. Static member function.
1379  *  The function return a pointer to a new created field that holds the addition.
1380  *  Data members are checked for compatibility and initialized.
1381  *  The user is in charge of memory deallocation.
1382  */
1383 template <class T, class INTERLACING_TAG>
1384 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::add(const FIELD& m, const FIELD& n)
1385 {
1386   const char* LOC = "FIELD<T>::add(const FIELD & m, const FIELD& n)";
1387   BEGIN_OF_MED(LOC);
1388     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1389
1390     // Creation of a new field
1391     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1392                                                                       m.getNumberOfComponents());
1393     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1394     result->_add_in_place(m,n); // perform addition
1395
1396   END_OF_MED(LOC);
1397     return result;
1398 }
1399
1400 /*! Same as add method except that field check is deeper.
1401  */
1402 template <class T, class INTERLACING_TAG>
1403 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::addDeep(const FIELD& m, const FIELD& n)
1404 {
1405   const char* LOC = "FIELD<T>::addDeep(const FIELD & m, const FIELD& n)";
1406   BEGIN_OF_MED(LOC);
1407     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1408
1409     // Creation of a new field
1410     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1411                                                                       m.getNumberOfComponents());
1412     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1413     result->_add_in_place(m,n); // perform addition
1414
1415   END_OF_MED(LOC);
1416     return result;
1417 }
1418
1419 /*!
1420      Overload substraction operator.
1421      This operation is authorized only for compatible fields that have the same support.
1422      The compatibility checking includes equality tests of the folowing data members:\n
1423      - _support
1424      - _numberOfComponents
1425      - _numberOfValues
1426      - _componentsTypes
1427      - _MEDComponentsUnits.
1428
1429      The data members of the returned field are initialized, based on the first field, except for the name,
1430      which is the combination of the two field's names and the operator.
1431      Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> \n
1432      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1433      When using python, this operator calls the copy constructor in any case.
1434      The user has to be aware that when using operator - in associatives expressions like
1435      <tt> a = b - c - d -e; </tt> \n
1436      no optimisation is performed : the evaluation of last expression requires the construction of
1437      3 temporary fields.
1438 */
1439 template <class T, class INTERLACING_TAG>
1440 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator-(const FIELD & m) const
1441 {
1442   const char* LOC = "FIELD<T>::operator-(const FIELD & m)";
1443   BEGIN_OF_MED(LOC);
1444     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1445
1446     // Creation of the result - memory is allocated by FIELD constructor
1447     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1448     //result._operation(*this,m,mode,"-"); // perform Atribute's initialization & substraction
1449     result._operationInitialize(*this,m,"-"); // perform Atribute's initialization
1450     result._sub_in_place(*this,m); // perform substracion
1451
1452   END_OF_MED(LOC);
1453     return result;
1454 }
1455
1456 template <class T, class INTERLACING_TAG>
1457 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator-() const
1458 {
1459   const char* LOC = "FIELD<T>::operator-()";
1460   BEGIN_OF_MED(LOC);
1461
1462     // Creation of the result - memory is allocated by FIELD constructor
1463     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1464     // Atribute's initialization
1465     result.setName("- "+getName());
1466     result.setComponentsNames(getComponentsNames());
1467     // not yet implemented    setComponentType(getComponentType());
1468     result.setComponentsDescriptions(getComponentsDescriptions());
1469     result.setMEDComponentsUnits(getMEDComponentsUnits());
1470     result.setComponentsUnits(getComponentsUnits());
1471     result.setIterationNumber(getIterationNumber());
1472     result.setTime(getTime());
1473     result.setOrderNumber(getOrderNumber());
1474
1475     const T* value1=getValue();
1476     // get a non const pointer to the inside array of values and perform operation
1477     T * value=const_cast<T *> (result.getValue());
1478     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1479     const T* endV=value+size; // pointer to the end of value
1480
1481     for(;value!=endV; value1++,value++)
1482         *value = -(*value1);
1483   END_OF_MED(LOC);
1484     return result;
1485 }
1486
1487 /*!  Overloaded Operator -=
1488  *   Operations are directly performed in the first field's array.
1489  *   This operation is authorized only for compatible fields that have the same support.
1490  */
1491 template <class T, class INTERLACING_TAG>
1492 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator-=(const FIELD & m)
1493 {
1494   const char* LOC = "FIELD<T>::operator-=(const FIELD & m)";
1495   BEGIN_OF_MED(LOC);
1496     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1497
1498     const T* value1=m.getValue();
1499
1500     // get a non const pointer to the inside array of values and perform operation
1501     T * value=const_cast<T *> (getValue());
1502     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1503     const T* endV=value+size; // pointer to the end of value
1504
1505     for(;value!=endV; value1++,value++)
1506         *value -= *value1;
1507
1508   END_OF_MED(LOC);
1509     return *this;
1510 }
1511
1512
1513
1514 /*!  Apply to a given field component the linear function x -> ax+b.
1515  *   calculation is done "in place".
1516  */
1517 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b, int icomp)
1518 {
1519     // get a non const pointer to the inside array of values and perform operation in place
1520     T * value=const_cast<T *> (getValue());
1521
1522     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1523
1524     if (size>0) // for a negative size, there is nothing to do
1525     {
1526                         value+=icomp-1;
1527                         const T* lastvalue=value+size; // pointer to the end of value
1528
1529                         for(;value!=lastvalue; value+=getNumberOfComponents()) // apply linear transformation
1530                                 *value = a*(*value)+b;
1531     }
1532 }
1533
1534 /*! Substraction of fields. Static member function.
1535  *  The function return a pointer to a new created field that holds the substraction.
1536  *  Data members are checked for compatibility and initialized.
1537  *  The user is in charge of memory deallocation.
1538  */
1539 template <class T, class INTERLACING_TAG>
1540 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::sub(const FIELD& m, const FIELD& n)
1541 {
1542   const char* LOC = "FIELD<T>::sub(const FIELD & m, const FIELD& n)";
1543   BEGIN_OF_MED(LOC);
1544     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1545
1546     // Creation of a new field
1547     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1548                                                                       m.getNumberOfComponents());
1549     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1550     result->_sub_in_place(m,n); // perform substraction
1551
1552   END_OF_MED(LOC);
1553     return result;
1554 }
1555
1556 /*! Same as sub method except that field check is deeper.
1557  */
1558 template <class T, class INTERLACING_TAG>
1559 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::subDeep(const FIELD& m, const FIELD& n)
1560 {
1561   const char* LOC = "FIELD<T>::subDeep(const FIELD & m, const FIELD& n)";
1562   BEGIN_OF_MED(LOC);
1563     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1564
1565     // Creation of a new field
1566     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1567                                                                       m.getNumberOfComponents());
1568     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1569     result->_sub_in_place(m,n); // perform substraction
1570
1571   END_OF_MED(LOC);
1572     return result;
1573 }
1574
1575 /*!
1576      Overload multiplication operator.
1577      This operation is authorized only for compatible fields that have the same support.
1578      The compatibility checking includes equality tests of the folowing data members:\n
1579      - _support
1580      - _numberOfComponents
1581      - _numberOfValues
1582      - _componentsTypes
1583      - _MEDComponentsUnits.
1584
1585      The data members of the returned field are initialized, based on the first field, except for the name,
1586      which is the combination of the two field's names and the operator.
1587      Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> \n
1588      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1589      When using python, this operator calls the copy constructor in any case.
1590      The user has to be aware that when using operator * in associatives expressions like
1591      <tt> a = b * c * d *e; </tt> \n
1592      no optimisation is performed : the evaluation of last expression requires the construction of
1593      3 temporary fields.
1594 */
1595 template <class T, class INTERLACING_TAG>
1596 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator*(const FIELD & m) const
1597 {
1598   const char* LOC = "FIELD<T>::operator*(const FIELD & m)";
1599   BEGIN_OF_MED(LOC);
1600     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1601
1602     // Creation of the result - memory is allocated by FIELD constructor
1603     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1604     //result._operation(*this,m,mode,"*"); // perform Atribute's initialization & multiplication
1605     result._operationInitialize(*this,m,"*"); // perform Atribute's initialization
1606     result._mul_in_place(*this,m); // perform multiplication
1607
1608   END_OF_MED(LOC);
1609     return result;
1610 }
1611
1612 /*!  Overloaded Operator *=
1613  *   Operations are directly performed in the first field's array.
1614  *   This operation is authorized only for compatible fields that have the same support.
1615  */
1616 template <class T, class INTERLACING_TAG>
1617 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator*=(const FIELD & m)
1618 {
1619   const char* LOC = "FIELD<T>::operator*=(const FIELD & m)";
1620   BEGIN_OF_MED(LOC);
1621     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1622
1623     const T* value1=m.getValue();
1624
1625     // get a non const pointer to the inside array of values and perform operation
1626     T * value=const_cast<T *> (getValue());
1627     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1628     const T* endV=value+size; // pointer to the end of value
1629
1630     for(;value!=endV; value1++,value++)
1631         *value *= *value1;
1632
1633   END_OF_MED(LOC);
1634     return *this;
1635 }
1636
1637
1638 /*! Multiplication of fields. Static member function.
1639  *  The function return a pointer to a new created field that holds the multiplication.
1640  *  Data members are checked for compatibility and initialized.
1641  *  The user is in charge of memory deallocation.
1642  */
1643 template <class T, class INTERLACING_TAG>
1644 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mul(const FIELD& m, const FIELD& n)
1645 {
1646   const char* LOC = "FIELD<T>::mul(const FIELD & m, const FIELD& n)";
1647   BEGIN_OF_MED(LOC);
1648     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1649
1650     // Creation of a new field
1651     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1652                                                                       m.getNumberOfComponents());
1653     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1654     result->_mul_in_place(m,n); // perform multiplication
1655
1656   END_OF_MED(LOC);
1657     return result;
1658 }
1659
1660 /*! Same as mul method except that field check is deeper.
1661  */
1662 template <class T, class INTERLACING_TAG>
1663 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mulDeep(const FIELD& m, const FIELD& n)
1664 {
1665   const char* LOC = "FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)";
1666   BEGIN_OF_MED(LOC);
1667     FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1668
1669     // Creation of a new field
1670     FIELD<T, INTERLACING_TAG>* result = new FIELD<T,INTERLACING_TAG>(m.getSupport(),
1671                                                                      m.getNumberOfComponents());
1672     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1673     result->_mul_in_place(m,n); // perform multiplication
1674
1675   END_OF_MED(LOC);
1676     return result;
1677 }
1678
1679 /*!
1680      Overload division operator.
1681      This operation is authorized only for compatible fields that have the same support.
1682      The compatibility checking includes equality tests of the folowing data members:\n
1683      - _support
1684      - _numberOfComponents
1685      - _numberOfValues
1686      - _componentsTypes
1687      - _MEDComponentsUnits.
1688
1689      The data members of the returned field are initialized, based on the first field, except for the name,
1690      which is the combination of the two field's names and the operator.
1691      Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> \n
1692      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1693      When using python, this operator calls the copy constructor in any case.
1694      The user has to be aware that when using operator / in associatives expressions like
1695      <tt> a = b / c / d /e; </tt> \n
1696      no optimisation is performed : the evaluation of last expression requires the construction of
1697      3 temporary fields.
1698 */
1699 template <class T, class INTERLACING_TAG>
1700 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator/(const FIELD & m) const
1701 {
1702   const char* LOC = "FIELD<T>::operator/(const FIELD & m)";
1703   BEGIN_OF_MED(LOC);
1704     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1705
1706     // Creation of the result - memory is allocated by FIELD constructor
1707     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1708     //result._operation(*this,m,mode,"/"); // perform Atribute's initialization & division
1709     result._operationInitialize(*this,m,"/"); // perform Atribute's initialization
1710     result._div_in_place(*this,m); // perform division
1711
1712   END_OF_MED(LOC);
1713     return result;
1714 }
1715
1716
1717 /*!  Overloaded Operator /=
1718  *   Operations are directly performed in the first field's array.
1719  *   This operation is authorized only for compatible fields that have the same support.
1720  */
1721 template <class T, class INTERLACING_TAG>
1722 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator/=(const FIELD & m)
1723 {
1724   const char* LOC = "FIELD<T>::operator/=(const FIELD & m)";
1725   BEGIN_OF_MED(LOC);
1726     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1727
1728     const T* value1=m.getValue(); // get pointers to the values we are adding
1729
1730     // get a non const pointer to the inside array of values and perform operation
1731     T * value=const_cast<T *> (getValue());
1732     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1733     const T* endV=value+size; // pointer to the end of value
1734
1735     for(;value!=endV; value1++,value++)
1736         *value /= *value1;
1737
1738   END_OF_MED(LOC);
1739     return *this;
1740 }
1741
1742
1743 /*! Division of fields. Static member function.
1744  *  The function return a pointer to a new created field that holds the division.
1745  *  Data members are checked for compatibility and initialized.
1746  *  The user is in charge of memory deallocation.
1747  */
1748 template <class T, class INTERLACING_TAG>
1749 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::div(const FIELD& m, const FIELD& n)
1750 {
1751   const char* LOC = "FIELD<T>::div(const FIELD & m, const FIELD& n)";
1752   BEGIN_OF_MED(LOC);
1753     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1754
1755     // Creation of a new field
1756     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1757                                                                       m.getNumberOfComponents());
1758     result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1759     result->_div_in_place(m,n); // perform division
1760
1761   END_OF_MED(LOC);
1762     return result;
1763 }
1764
1765 /*! Same as div method except that field check is deeper.
1766  */
1767 template <class T,class INTERLACING_TAG>
1768 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::divDeep(const FIELD& m, const FIELD& n)
1769 {
1770   const char* LOC = "FIELD<T>::divDeep(const FIELD & m, const FIELD& n)";
1771   BEGIN_OF_MED(LOC);
1772   FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1773
1774   // Creation of a new field
1775   FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1776                                                                     m.getNumberOfComponents());
1777   result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1778   result->_div_in_place(m,n); // perform division
1779
1780   END_OF_MED(LOC);
1781   return result;
1782 }
1783
1784
1785 /*!
1786 @}
1787 */
1788
1789 /*!
1790   \if developper
1791   This internal method initialize the members of a new field created to hold the result of the operation Op .
1792   Initialization is based on the first field, except for the name, which is the combination of the two field's names
1793   and the operator.
1794   \endif
1795 */
1796 template <class T, class INTERLACING_TAG>
1797 void FIELD<T, INTERLACING_TAG>::_operationInitialize(const FIELD& m,const FIELD& n, const char* Op)
1798 {
1799     MESSAGE_MED("Appel methode interne " << Op);
1800
1801     // Atribute's initialization (copy of the first field's attributes)
1802     // Other data members (_support, _numberOfValues) are initialized in the field's constr.
1803     setName(m.getName()+" "+Op+" "+n.getName());
1804     setComponentsNames(m.getComponentsNames());
1805     // not yet implemented    setComponentType(m.getComponentType());
1806     setComponentsDescriptions(m.getComponentsDescriptions());
1807     setMEDComponentsUnits(m.getMEDComponentsUnits());
1808
1809     // The following data member may differ from field m to n.
1810     // The initialization is done based on the first field.
1811
1812     //if(m.getComponentsUnits() != NULL)
1813     setComponentsUnits(m.getComponentsUnits());
1814     //else
1815     //  _componentsUnits = (UNIT *) NULL;
1816
1817     setIterationNumber(m.getIterationNumber());
1818     setTime(m.getTime());
1819     setOrderNumber(m.getOrderNumber());
1820 }
1821
1822
1823 /*!
1824   \if developper
1825   Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
1826   This method is applied to a just created field with medModeSwitch mode.
1827   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1828   it doesn't exist!
1829   \endif
1830 */
1831 template <class T, class INTERLACING_TAG>
1832 void FIELD<T, INTERLACING_TAG>::_add_in_place(const FIELD& m,const FIELD& n)
1833 {
1834     // get pointers to the values we are adding
1835     const T* value1=m.getValue();
1836     const T* value2=n.getValue();
1837     // get a non const pointer to the inside array of values and perform operation
1838     T * value=const_cast<T *> (getValue());
1839
1840     const int size=getNumberOfValues()*getNumberOfComponents();
1841     SCRUTE_MED(size);
1842     const T* endV1=value1+size;
1843     for(;value1!=endV1; value1++,value2++,value++)
1844         *value=(*value1)+(*value2);
1845 }
1846
1847 /*!
1848   \if developper
1849   Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
1850   This method is applied to a just created field with medModeSwitch mode.
1851   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1852   it doesn't exist!
1853   \endif
1854 */
1855 template <class T, class INTERLACING_TAG>
1856 void FIELD<T, INTERLACING_TAG>::_sub_in_place(const FIELD& m,const FIELD& n)
1857 {
1858     // get pointers to the values we are adding
1859     const T* value1=m.getValue();
1860     const T* value2=n.getValue();
1861     // get a non const pointer to the inside array of values and perform operation
1862     T * value=const_cast<T *> (getValue());
1863
1864     const int size=getNumberOfValues()*getNumberOfComponents();
1865     SCRUTE_MED(size);
1866     const T* endV1=value1+size;
1867     for(;value1!=endV1; value1++,value2++,value++)
1868         *value=(*value1)-(*value2);
1869 }
1870
1871 /*!
1872   \if developper
1873   Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
1874   This method is applied to a just created field with medModeSwitch mode.
1875   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1876   it doesn't exist!
1877   \endif
1878 */
1879 template <class T, class INTERLACING_TAG>
1880 void FIELD<T, INTERLACING_TAG>::_mul_in_place(const FIELD& m,const FIELD& n)
1881 {
1882     // get pointers to the values we are adding
1883     const T* value1=m.getValue();
1884     const T* value2=n.getValue();
1885     // get a non const pointer to the inside array of values and perform operation
1886     T * value=const_cast<T *> (getValue());
1887
1888     const int size=getNumberOfValues()*getNumberOfComponents();
1889     SCRUTE_MED(size);
1890     const T* endV1=value1+size;
1891     for(;value1!=endV1; value1++,value2++,value++)
1892         *value=(*value1)*(*value2);
1893 }
1894
1895 /*!
1896   \if developper
1897   Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
1898   This method is applied to a just created field with medModeSwitch mode.
1899   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1900   it doesn't exist!
1901   \endif
1902 */
1903 template <class T, class INTERLACING_TAG>
1904 void FIELD<T, INTERLACING_TAG>::_div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION)
1905 {
1906     // get pointers to the values we are adding
1907     const T* value1=m.getValue();
1908     const T* value2=n.getValue();
1909     // get a non const pointer to the inside array of values and perform operation
1910     T * value=const_cast<T *> (getValue());
1911
1912     const int size=getNumberOfValues()*getNumberOfComponents();
1913     SCRUTE_MED(size);
1914     const T* endV1=value1+size;
1915     for(;value1!=endV1; value1++,value2++,value++){
1916       if ( *value2 == 0 ) { // FAIRE PLUTOT UN TRY CATCH Rmq from EF
1917           string diagnosis;
1918           diagnosis="FIELD<T,INTERLACING_TAG>::_div_in_place(...) : Divide by zero !";
1919           throw MEDEXCEPTION(diagnosis.c_str());
1920         }
1921         *value=(*value1)/(*value2);
1922     }
1923 }
1924
1925 /*!
1926 \addtogroup FIELD_algo
1927 @{
1928 */
1929
1930 /*!  Return maximum of all absolute values contained in the array (all elements and all components are browsed).
1931  */
1932 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::normMax() const throw (MEDEXCEPTION)
1933 {
1934     const T* value=getValue(); // get pointer to the values
1935     const int size=getNumberOfValues()*getNumberOfComponents();
1936     if (size <= 0) // Size of array has to be strictly positive
1937     {
1938         string diagnosis;
1939         diagnosis="FIELD<T,INTERLACIN_TAG>::normMax() : cannot compute the norm of "+getName()+
1940             " : it size is non positive!";
1941         throw MEDEXCEPTION(diagnosis.c_str());
1942     }
1943     const T* lastvalue=value+size; // get pointer just after last value
1944     const T* pMax=value; // pointer to the max value
1945     const T* pMin=value; // pointer to the min value
1946
1947     // get pointers to the max & min value of array
1948     while ( ++value != lastvalue )
1949     {
1950         if ( *pMin > *value )
1951             pMin=value;
1952         if ( *pMax < *value )
1953             pMax=value;
1954     }
1955
1956     T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
1957     T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
1958
1959     return Max>Min ? static_cast<double>(Max) : static_cast<double>(Min);
1960 }
1961
1962 /*!  Return Euclidian norm for all elements of the array.
1963  */
1964 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::norm2() const throw (MEDEXCEPTION)
1965 {
1966     const T* value=this->getValue(); // get const pointer to the values
1967     const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
1968     if (size <= 0) // Size of array has to be strictly positive
1969     {
1970         string diagnosis;
1971         diagnosis="FIELD<T,INTERLACIN_TAG>::norm2() : cannot compute the norm of "+getName()+
1972             " : it size is non positive!";
1973         throw MEDEXCEPTION(diagnosis.c_str());
1974     }
1975     const T* lastvalue=value+size; // point just after last value
1976
1977     T result((T)0); // init
1978     for( ; value!=lastvalue ; ++value)
1979         result += (*value) * (*value);
1980
1981     return std::sqrt(static_cast<double> (result));
1982 }
1983
1984
1985 //------------- TDG and BS addings
1986
1987 /*!  Return Extrema of field
1988  */
1989  template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION)
1990 {
1991   const T* value=getValue(); // get pointer to the values
1992   const int size=getNumberOfValues()*getNumberOfComponents();
1993   const T* lastvalue=value+size; // point just after last value
1994
1995   if (size <= 0){ // Size of array has to be strictly positive
1996
1997     string diagnosis;
1998     diagnosis="FIELD<T,INTERLACIN_TAG>::getMinMax() : cannot compute the extremums of "+getName()+
1999       " : its size is non positive!";
2000     throw MEDEXCEPTION(diagnosis.c_str());
2001   }
2002
2003   if (!_isMinMax){
2004     vmax=MinMax<T>::getMin(); // init a max value
2005     vmin=MinMax<T>::getMax(); // init a min value
2006
2007     for( ; value!=lastvalue ; ++value){
2008       if ( vmin > *value )
2009         vmin=*value;
2010       if ( vmax < *value )
2011         vmax=*value;
2012     }
2013     _isMinMax=true;
2014     _vmin=vmin;
2015     _vmax=vmax;
2016   }
2017   else{
2018     vmin = _vmin;
2019     vmax = _vmax;
2020   }
2021
2022 }
2023
2024 /*!  Return Histogram of field
2025  */
2026  template <class T, class INTERLACIN_TAG> vector<int> FIELD<T, INTERLACIN_TAG>::getHistogram(int &nbint) throw (MEDEXCEPTION)
2027 {
2028   const T* value=getValue(); // get pointer to the values
2029   const int size=getNumberOfValues()*getNumberOfComponents();
2030   const T* lastvalue=value+size; // point just after last value
2031
2032   if (size <= 0){ // Size of array has to be strictly positive
2033
2034     string diagnosis;
2035     diagnosis="FIELD<T,INTERLACIN_TAG>::getHistogram() : cannot compute the histogram of "+getName()+
2036       " : it size is non positive!";
2037     throw MEDEXCEPTION(diagnosis.c_str());
2038   }
2039   //    return static_cast<ArrayGauss *>(_value)->getIJ(valIndex,j) ;
2040
2041   vector<int> Histogram(nbint) ;
2042   T vmin,vmax;
2043   int j;
2044
2045   for( j=0 ; j!=nbint ; j++) Histogram[j]=0 ;
2046
2047   getMinMax(vmin,vmax);
2048   for( ; value!=lastvalue ; ++value){
2049     if(*value==vmax) j = nbint-1;
2050     else j = (int)(((double)nbint * (*value-vmin))/(vmax-vmin));
2051     Histogram[j]+=1 ;
2052   }
2053
2054   return Histogram ;
2055
2056 }
2057
2058 /*!  Return vectorial gradient field
2059  */
2060 template <class T, class INTERLACIN_TAG>
2061 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildGradient() const throw (MEDEXCEPTION)
2062 {
2063   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildGradient() : ";
2064   BEGIN_OF_MED(LOC);
2065
2066   // space dimension of input mesh
2067   int spaceDim = getSupport()->getMesh()->getSpaceDimension();
2068   double *x = new double[spaceDim];
2069
2070   FIELD<double, FullInterlace>* Gradient =
2071     new FIELD<double, FullInterlace>(getSupport(),spaceDim);
2072
2073   string name("gradient of ");
2074   name += getName();
2075   Gradient->setName(name);
2076   string descr("gradient of ");
2077   descr += getDescription();
2078   Gradient->setDescription(descr);
2079
2080   if( _numberOfComponents > 1 ){
2081     delete Gradient;
2082     delete [] x;
2083     throw MEDEXCEPTION("gradient calculation only on scalar field");
2084   }
2085
2086   for(int i=1;i<=spaceDim;i++){
2087     string nameC("gradient of ");
2088     nameC += getName();
2089     Gradient->setComponentName(i,nameC);
2090     Gradient->setComponentDescription(i,"gradient");
2091     string MEDComponentUnit = getMEDComponentUnit(1)+getSupport()->getMesh()->getCoordinatesUnits()[i-1];
2092     Gradient->setMEDComponentUnit(i,MEDComponentUnit);
2093   }
2094
2095   Gradient->setIterationNumber(getIterationNumber());
2096   Gradient->setOrderNumber(getOrderNumber());
2097   Gradient->setTime(getTime());
2098
2099   // typ of entity on what is field
2100   MED_EN::medEntityMesh typ = getSupport()->getEntity();
2101
2102   const int *C;
2103   const int *iC;
2104   const int *revC;
2105   const int *indC;
2106   const double *coord;
2107   int NumberOf;
2108
2109   switch (typ) {
2110   case MED_CELL:
2111   case MED_FACE:
2112   case MED_EDGE:
2113     {
2114       // read connectivity array to have the list of nodes contained by an element
2115       C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,typ,MED_ALL_ELEMENTS);
2116       iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,typ);
2117       // calculate reverse connectivity to have the list of elements which contains node i
2118       revC = getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
2119       indC = getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
2120       // coordinates of each node
2121       coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
2122       // number of elements
2123       NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2124       // barycenter field of elements
2125       FIELD<double, FullInterlace>* barycenter = getSupport()->getMesh()->getBarycenter(getSupport());
2126
2127       // calculate gradient vector for each element i
2128       for (int i = 1; i < NumberOf + 1; i++) {
2129
2130         // listElements contains elements which contains a node of element i
2131         set <int> listElements;
2132         set <int>::iterator elemIt;
2133         listElements.clear();
2134
2135         // for each node j of element i
2136         for (int ij = iC[i-1]; ij < iC[i]; ij++) {
2137           int j = C[ij-1];
2138           for (int k = indC[j-1]; k < indC[j]; k++) {
2139             // c element contains node j
2140             int c = revC[k-1];
2141             // we put the elements in set
2142             if (c != i)
2143               listElements.insert(c);
2144           }
2145         }
2146         // coordinates of barycentre of element i in space of dimension spaceDim
2147         for (int j = 0; j < spaceDim; j++)
2148           x[j] = barycenter->getValueIJ(i,j+1);
2149
2150         for (int j = 0; j < spaceDim; j++) {
2151           // value of field of element i
2152           double val = getValueIJ(i,1);
2153           double grad = 0.;
2154           // calculate gradient for each neighbor element
2155           for (elemIt = listElements.begin(); elemIt != listElements.end(); elemIt++) {
2156             int elem = *elemIt;
2157             double d2 = 0.;
2158             for (int l = 0; l < spaceDim; l++) {
2159               // coordinate of barycenter of element elem
2160               double xx = barycenter->getValueIJ(elem, l+1);
2161               d2 += (x[l]-xx) * (x[l]-xx);
2162             }
2163             grad += (barycenter->getValueIJ(elem,j+1)-x[j])*(getValueIJ(elem,1)-val)/sqrt(d2);
2164           }
2165           if (listElements.size() != 0) grad /= listElements.size();
2166           Gradient->setValueIJ(i,j+1,grad);
2167         }
2168       }
2169       delete barycenter;
2170     }
2171     break;
2172   case MED_NODE:
2173     // read connectivity array to have the list of nodes contained by an element
2174     C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
2175     iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,MED_CELL);
2176     // calculate reverse connectivity to have the list of elements which contains node i
2177     revC=getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
2178     indC=getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
2179     // coordinates of each node
2180     coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
2181
2182     // calculate gradient for each node
2183     NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2184     for (int i=1; i<NumberOf+1; i++){
2185       // listNodes contains nodes neigbor of node i
2186       set <int> listNodes;
2187       set <int>::iterator nodeIt ;
2188       listNodes.clear();
2189       for(int j=indC[i-1];j<indC[i];j++){
2190         // c element contains node i
2191         int c=revC[j-1];
2192         // we put the nodes of c element in set
2193         for(int k=iC[c-1];k<iC[c];k++)
2194           if(C[k-1] != i)
2195             listNodes.insert(C[k-1]);
2196       }
2197       // coordinates of node i in space of dimension spaceDim
2198       for(int j=0;j<spaceDim;j++)
2199         x[j] = coord[(i-1)*spaceDim+j];
2200
2201       for(int j=0;j<spaceDim;j++){
2202         // value of field
2203         double val = getValueIJ(i,1);
2204         double grad = 0.;
2205         // calculate gradient for each neighbor node
2206         for(nodeIt=listNodes.begin();nodeIt!=listNodes.end();nodeIt++){
2207           int node = *nodeIt;
2208           double d2 = 0.;
2209           for(int l=0;l<spaceDim;l++){
2210             double xx = coord[(node-1)*spaceDim+l];
2211             d2 += (x[l]-xx) * (x[l]-xx);
2212           }
2213           grad += (coord[(node-1)*spaceDim+j]-x[j])*(getValueIJ(node,1)-val)/sqrt(d2);
2214         }
2215         if(listNodes.size() != 0) grad /= listNodes.size();
2216         Gradient->setValueIJ(i,j+1,grad);
2217       }
2218     }
2219     break;
2220   case MED_ALL_ENTITIES:
2221     delete [] x;
2222     delete Gradient;
2223     throw MEDEXCEPTION("gradient calculation not yet implemented on all elements");
2224     break;
2225   }
2226
2227   delete [] x;
2228
2229   END_OF_MED(LOC);
2230   return Gradient;
2231 }
2232
2233 /*!  Return scalar norm2 field
2234  */
2235 template <class T, class INTERLACIN_TAG>
2236 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildNorm2Field() const throw (MEDEXCEPTION)
2237 {
2238   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildNorm2Field() : ";
2239   BEGIN_OF_MED(LOC);
2240
2241   FIELD<double, FullInterlace>* Norm2Field =
2242     new FIELD<double, FullInterlace>(getSupport(),1);
2243
2244   string name("norm2 of ");
2245   name += getName();
2246   Norm2Field->setName(name);
2247   string descr("norm2 of ");
2248   descr += getDescription();
2249   Norm2Field->setDescription(descr);
2250
2251   string nameC("norm2 of ");
2252   nameC += getName();
2253   Norm2Field->setComponentName(1,nameC);
2254   Norm2Field->setComponentDescription(1,"norm2");
2255   string MEDComponentUnit = getMEDComponentUnit(1);
2256   Norm2Field->setMEDComponentUnit(1,MEDComponentUnit);
2257
2258   Norm2Field->setIterationNumber(getIterationNumber());
2259   Norm2Field->setOrderNumber(getOrderNumber());
2260   Norm2Field->setTime(getTime());
2261
2262   // calculate nom2 for each element
2263   int NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2264   for (int i=1; i<NumberOf+1; i++){
2265     double norm2 = 0.;
2266     for(int j=1;j<=getNumberOfComponents();j++)
2267       norm2 += getValueIJ(i,j)*getValueIJ(i,j);
2268     Norm2Field->setValueIJ(i,1,sqrt(norm2));
2269   }
2270
2271   END_OF_MED(LOC);
2272   return Norm2Field;
2273
2274 }
2275
2276 /*!  Apply to each (scalar) field component the template parameter T_function,
2277  *   which is a pointer to function.
2278  *   Since the pointer is known at compile time, the function is inlined into the inner loop!
2279  *   calculation is done "in place".
2280  *   Use examples :
2281  *
2282  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
2283  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
2284  */
2285 template <class T, class INTERLACIN_TAG> template <T T_function(T)>
2286 void FIELD<T, INTERLACIN_TAG>::applyFunc()
2287 {
2288   // get a non const pointer to the inside array of values and perform operation
2289     T * value=const_cast<T *> (getValue());
2290     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
2291
2292     if (size>0) // for a negative size, there is nothing to do
2293     {
2294         const T* lastvalue=value+size; // pointer to the end of value
2295         for(;value!=lastvalue; ++value) // apply linear transformation
2296             *value = T_function(*value);
2297     }
2298 }
2299
2300 template <class T, class INTERLACIN_TAG> T FIELD<T, INTERLACIN_TAG>::pow(T x)
2301 {
2302   return (T)::pow((double)x,FIELD<T, INTERLACIN_TAG>::_scalarForPow);
2303 }
2304
2305 /*!  Apply to each (scalar) field component the math function pow.
2306  *   calculation is done "in place".
2307  *   Use examples :
2308  *
2309  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
2310  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
2311  */
2312 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyPow(T scalar)
2313 {
2314   FIELD<T, INTERLACIN_TAG>::_scalarForPow=scalar;
2315   applyFunc<FIELD<T, INTERLACIN_TAG>::pow>();
2316 }
2317
2318 /*!  Apply to each (scalar) field component the linear function x -> ax+b.
2319  *   calculation is done "in place".
2320  */
2321 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b)
2322 {
2323     // get a non const pointer to the inside array of values and perform operation in place
2324     T * value=const_cast<T *> (getValue());
2325     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
2326
2327     if (size>0) // for a negative size, there is nothing to do
2328     {
2329         const T* lastvalue=value+size; // pointer to the end of value
2330         for(;value!=lastvalue; ++value) // apply linear transformation
2331             *value = a*(*value)+b;
2332     }
2333 }
2334
2335
2336 /*!
2337  *   Return a pointer to a new field that holds the scalar product. Static member function.
2338  *   This operation is authorized only for compatible fields that have the same support.
2339  *   The compatibility checking includes equality tests of the folowing data members:\n
2340  *   - _support
2341  *   - _numberOfComponents
2342  *   - _numberOfValues
2343  *   - _componentsTypes
2344  *   - _MEDComponentsUnits.
2345  *   Data members are initialized.
2346  *   The new field point to the same support and has one component.
2347  *   Each value of it is the scalar product of the two argument's fields.
2348  *   The user is in charge of memory deallocation.
2349  */
2350 template <class T, class INTERLACING_TAG>
2351 FIELD<T, INTERLACING_TAG>*
2352 FIELD<T, INTERLACING_TAG>::scalarProduct(const FIELD & m, const FIELD & n, bool deepCheck)
2353 {
2354     if(!deepCheck)
2355       FIELD_::_checkFieldCompatibility( m, n, false); // may throw exception
2356     else
2357       FIELD_::_deepCheckFieldCompatibility(m, n, false);
2358
2359     // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
2360     // result type imply INTERLACING_TAG=FullInterlace for m & n
2361
2362     const int numberOfElements=m.getNumberOfValues(); // strictly positive
2363     const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
2364
2365     // Creation & init of a the result field on the same support, with one component
2366     // You have to be careful about the interlacing mode, because in the computation step,
2367     // it seems to assume the the interlacing mode is the FullInterlacing
2368
2369     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),1);
2370     result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
2371     result->setIterationNumber(m.getIterationNumber());
2372     result->setTime(m.getTime());
2373     result->setOrderNumber(m.getOrderNumber());
2374
2375     const T* value1=m.getValue(); // get const pointer to the values
2376     const T* value2=n.getValue(); // get const pointer to the values
2377     // get a non const pointer to the inside array of values and perform operation
2378     T * value=const_cast<T *> (result->getValue());
2379
2380     const T* lastvalue=value+numberOfElements; // pointing just after last value of result
2381     for ( ; value!=lastvalue ; ++value ) // loop on all elements
2382     {
2383         *value=(T)0; // initialize value
2384         const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
2385         for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
2386             *value += (*value1) * (*value2);
2387     }
2388     return result;
2389 }
2390
2391 /*!  Return L2 Norm  of the field's component.
2392  *   Cannot be applied to a field with a support on nodes.
2393  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2394  */
2395 template <class T, class INTERLACING_TAG>
2396 double FIELD<T, INTERLACING_TAG>::normL2(int component,
2397                                          const FIELD<double, FullInterlace> * p_field_volume) const
2398 {
2399     _checkNormCompatibility(p_field_volume); // may throw exception
2400     if ( component<1 || component>getNumberOfComponents() )
2401         throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
2402
2403     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2404     if(!p_field_volume) // if the user don't supply the volume
2405         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
2406
2407     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2408     const double* vol=p_field_size->getValue();
2409     // Il n'est vraiment pas optimal de mixer des champs dans des modes d'entrelacement
2410     // different juste pour le calcul
2411
2412
2413     double integrale=0.0;
2414     double totVol=0.0;
2415
2416     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2417       const T* value = getValue();
2418       value = value + (component-1) * getNumberOfValues();
2419       const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
2420       for (; value!=lastvalue ; ++value ,++vol) {
2421         integrale += static_cast<double>((*value) * (*value)) * (*vol);
2422         totVol+=*vol;
2423       }
2424     }
2425     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2426       ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2427       //for (int i=1; i <= anArray->getNbElem() ; i++ ) {
2428       //  for (int j=1; j<= anArray->getDim(); j++, ++vol ) {
2429       //    integrale += static_cast<double>( anArray->getIJ(i,j) * anArray->getIJ(i,j) * (*vol) );
2430       //    totVol+=*vol;
2431       //  }
2432       //}
2433       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2434         integrale += static_cast<double>( anArray->getIJ(i,component) * anArray->getIJ(i,component) * (*vol) );
2435         totVol+=*vol;
2436       }
2437       //delete anArray;
2438     }
2439     else { // FULL_INTERLACE
2440       ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2441       //for (int i=1; i <= anArray->getNbElem() ; i++ ) {
2442       //  for (int j=1; j<= anArray->getDim(); j++, ++vol ) {
2443       //for (int j=1; j<= anArray->getDim(); j++ ) {
2444       //  for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2445       //    integrale += static_cast<double>( anArray->getIJ(i,j) * anArray->getIJ(i,j) * (*vol) );
2446       //    totVol+=*vol;
2447       //  }
2448       //}
2449       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2450         integrale += static_cast<double>( anArray->getIJ(i,component) * anArray->getIJ(i,component) * (*vol) );
2451         totVol+=*vol;
2452       }
2453       //delete anArray;
2454     }
2455
2456     //const T * value     = NULL;
2457     //ArrayNo * myArray   = NULL;
2458     //if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2459     //  value = getValue();
2460     //else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2461     //  myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2462     //  value   = myArray->getPtr();
2463     //}
2464     //else {
2465     //  myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2466     //  value   = myArray->getPtr();
2467     //}
2468
2469     //value = value + (component-1) * getNumberOfValues();
2470     //const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
2471
2472     //double integrale=0.0;
2473     //double totVol=0.0;
2474     //for (; value!=lastvalue ; ++value ,++vol)
2475     //{
2476     //integrale += static_cast<double>((*value) * (*value)) * (*vol);
2477     //totVol+=*vol;
2478     //}
2479
2480     if(!p_field_volume) // if the user didn't supply the volume
2481         delete p_field_size; // delete temporary volume field
2482     //if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2483     if( totVol <= 0)
2484         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2485
2486     return integrale/totVol;
2487 }
2488
2489 /*!  Return L2 Norm  of the field.
2490  *   Cannot be applied to a field with a support on nodes.
2491  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2492  */
2493 template <class T, class INTERLACING_TAG>
2494 double FIELD<T, INTERLACING_TAG>::normL2(const FIELD<double, FullInterlace> * p_field_volume) const
2495 {
2496     _checkNormCompatibility(p_field_volume); // may throw exception
2497     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2498     if(!p_field_volume) // if the user don't supply the volume
2499         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃ©mentation dans mesh]
2500
2501     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2502     const double* vol=p_field_size->getValue();
2503     const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
2504
2505
2506     double integrale=0.0;
2507     double totVol=0.0;
2508     const double* p_vol=vol;
2509     for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2510         totVol+=*p_vol;
2511
2512     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2513       const T* value = getValue();
2514       for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
2515         for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
2516           integrale += static_cast<double>((*value) * (*value)) * (*p_vol);
2517         }
2518       }
2519     }
2520     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2521       ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2522       for (int j=1; j<=anArray->getDim(); j++) {
2523         int i = 1;
2524         for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2525           integrale += static_cast<double>( anArray->getIJ(i,j) * anArray->getIJ(i,j) * (*p_vol) );
2526         }
2527       }
2528       //delete anArray;
2529     }
2530     else { // FULL_INTERLACE
2531       ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2532       for (int j=1; j<=anArray->getDim(); j++) {
2533         int i = 1;
2534         for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2535           integrale += static_cast<double>( anArray->getIJ(i,j) * anArray->getIJ(i,j) * (*p_vol) );
2536         }
2537       }
2538       //delete anArray;
2539     }
2540
2541
2542     //const T * value     = NULL;
2543     //ArrayNo * myArray   = NULL;
2544     //if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2545     //  value = getValue();
2546     //else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ){
2547     //  myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2548     //  value   = myArray->getPtr();
2549     //}
2550     //else {
2551     //  myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2552     //  value   = myArray->getPtr();
2553     //}
2554
2555     //double totVol=0.0;
2556     //const double* p_vol=vol;
2557     //for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2558     //totVol+=*p_vol;
2559
2560     //double integrale=0.0;
2561     //for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
2562     //for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
2563     //    integrale += static_cast<double>((*value) * (*value)) * (*p_vol);
2564
2565     if(!p_field_volume) // if the user didn't supply the volume
2566         delete p_field_size; // delete temporary volume field
2567     //if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2568     if( totVol <= 0)
2569         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2570
2571     return integrale/totVol;
2572 }
2573
2574 /*!  Return L1 Norm  of the field's component.
2575  *   Cannot be applied to a field with a support on nodes.
2576  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2577  */
2578 template <class T, class INTERLACING_TAG>
2579 double FIELD<T, INTERLACING_TAG>::normL1(int component,
2580                                          const FIELD<double, FullInterlace > * p_field_volume) const
2581 {
2582     _checkNormCompatibility(p_field_volume); // may throw exception
2583     if ( component<1 || component>getNumberOfComponents() )
2584         throw MEDEXCEPTION(STRING("FIELD<T,INTERLACING_TAG>::normL2() : The component argument should be between 1 and the number of components"));
2585
2586     const FIELD<double,FullInterlace> * p_field_size=p_field_volume;
2587     if(!p_field_volume) // if the user don't supply the volume
2588         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃ©mentation dans mesh]
2589
2590     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2591     const double* vol = p_field_size->getValue();
2592
2593     double integrale=0.0;
2594     double totVol=0.0;
2595
2596     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2597       const T* value = getValue();
2598       const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
2599       for (; value!=lastvalue ; ++value ,++vol) {
2600         integrale += std::fabs( static_cast<double>(*value) ) * (*vol);
2601         totVol+=*vol;
2602       }
2603     }
2604     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2605       ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2606       //for (int i=1; i <= anArray->getNbElem() ; i++ ) {
2607       //  for (int j=1; j<= anArray->getDim(); j++, ++vol ) {
2608       //    integrale += std::abs(static_cast<double>( anArray->getIJ(i,j)) ) * (*vol);
2609       //    totVol+=*vol;
2610       //  }
2611       //}
2612       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2613         integrale += std::abs(static_cast<double>( anArray->getIJ(i,component)) ) * (*vol);
2614         totVol+=*vol;
2615       }
2616       //delete anArray;
2617     }
2618     else { // FULL_INTERLACE
2619       ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2620       //for (int i=1; i <= anArray->getNbElem() ; i++ ) {
2621       //  for (int j=1; j<= anArray->getDim(); j++, ++vol ) {
2622       //    integrale += std::abs(static_cast<double>( anArray->getIJ(i,j)) ) * (*vol);
2623       //    totVol+=*vol;
2624       //  }
2625       //}
2626       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2627         integrale += std::abs(static_cast<double>( anArray->getIJ(i,component)) ) * (*vol);
2628         totVol+=*vol;
2629       }
2630       //delete anArray;
2631     }
2632
2633     //const T * value     = NULL;
2634     //ArrayNo * myArray   = NULL;
2635     //if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2636     //  value = getColumn(component);
2637     //else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2638     //  myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2639     //  value   = myArray->getColumn(component);
2640     //}
2641     //else {
2642     //  myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2643     //  value   = myArray->getColumn(component);
2644     //}
2645
2646     //const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
2647
2648     //double integrale=0.0;
2649     //double totVol=0.0;
2650     //for (; value!=lastvalue ; ++value ,++vol)
2651     //{
2652     //integrale += std::abs( static_cast<double>(*value) ) * (*vol);
2653     //totVol+=*vol;
2654     //}
2655
2656     if(!p_field_volume) // if the user didn't supply the volume
2657         delete p_field_size; // delete temporary volume field
2658     //if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2659     if( totVol <= 0)
2660         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2661
2662     return integrale/totVol;
2663 }
2664
2665 /*!  Return L1 Norm  of the field.
2666  *   Cannot be applied to a field with a support on nodes.
2667  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2668  */
2669 template <class T, class INTERLACING_TAG>
2670 double FIELD<T, INTERLACING_TAG>::normL1(const FIELD<double, FullInterlace> * p_field_volume) const
2671 {
2672   _checkNormCompatibility(p_field_volume); // may throw exception
2673   const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2674   if(!p_field_volume) // if the user don't supply the volume
2675     p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃÂ©mentation dans mesh]
2676
2677   // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2678   const double* vol = p_field_size->getValue();
2679   const double* lastvol = vol+getNumberOfValues(); // pointing just after the end of vol
2680
2681   double integrale=0.0;
2682   double totVol=0.0;
2683   const double* p_vol=vol;
2684   for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2685     totVol+=*p_vol;
2686
2687   if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2688     const T* value = getValue();
2689     for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
2690       for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
2691         integrale += std::abs( static_cast<double>(*value) ) * (*p_vol);
2692       }
2693     }
2694   }
2695   else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2696     ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2697     for (int j=1; j<=anArray->getDim(); j++) {
2698       int i = 1;
2699       for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2700         integrale += std::abs(static_cast<double>(anArray->getIJ(i,j))) * (*p_vol);
2701       }
2702     }
2703     //delete anArray;
2704   }
2705   else { // FULL_INTERLACE
2706     ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2707     for (int j=1; j<=anArray->getDim(); j++) {
2708       int i = 1;
2709       for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2710         integrale += std::abs(static_cast<double>(anArray->getIJ(i,j))) * (*p_vol);
2711       }
2712     }
2713     //delete anArray;
2714   }
2715
2716
2717   //const T * value     = NULL;
2718   //ArrayNo * myArray   = NULL;
2719   //if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2720   //  value = getValue();
2721   //else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2722   //  myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2723   //  value   = myArray->getPtr();
2724   //}
2725   //else {
2726   //  myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2727   //  value   = myArray->getPtr();
2728   //}
2729
2730   //double totVol=0.0;
2731   //const double* p_vol=vol;
2732   //for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2733   //  totVol+=*p_vol;
2734
2735   //double integrale=0.0;
2736   //for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
2737   //  for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
2738   //    integrale += std::abs( static_cast<double>(*value) ) * (*p_vol);
2739
2740   if(!p_field_volume) // if the user didn't supply the volume
2741     delete p_field_size; // delete temporary volume field
2742   //if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2743   if( totVol <= 0)
2744     throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2745
2746   return integrale/totVol;
2747 }
2748
2749 /*! Return a new field (to deallocate with delete) lying on subSupport that is included by
2750  *   this->_support with corresponding values extracting from this->_value.
2751  */
2752 template <class T, class INTERLACING_TAG>
2753 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
2754 {
2755   if(!subSupport->belongsTo(*_support))
2756     throw MEDEXCEPTION("FIELD<T>::extract : subSupport not included in this->_support !");
2757   if(_support->isOnAllElements() && subSupport->isOnAllElements())
2758     return new FIELD<T, INTERLACING_TAG>(*this);
2759
2760   FIELD<T, INTERLACING_TAG> *ret = new FIELD<T, INTERLACING_TAG>(subSupport,
2761                                                                  _numberOfComponents);
2762
2763   if(!ret->_value)
2764     throw MEDEXCEPTION("FIELD<T>::extract : unvalid support detected !");
2765
2766   T* valuesToSet=(T*)ret->getValue();
2767
2768   int nbOfEltsSub=subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2769   const int *eltsSub=subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS);
2770   T* tempVals=new T[_numberOfComponents];
2771   for(int i=0;i<nbOfEltsSub;i++)
2772     {
2773       if(!getValueOnElement(eltsSub[i],tempVals))
2774         throw MEDEXCEPTION("Problem in belongsTo function !!!");
2775       for(int j=0;j<_numberOfComponents;j++)
2776         valuesToSet[i*_numberOfComponents+j]=tempVals[j];
2777     }
2778   delete [] tempVals;
2779
2780   ret->copyGlobalInfo(*this);
2781   return ret;
2782 }
2783 /*!
2784 @}
2785 */
2786
2787 /*!
2788         \addtogroup FIELD_io
2789         @{
2790 */
2791 /*!
2792   Constructor with parameters; the object is set via a file and its associated
2793   driver. For the moment only the MED_DRIVER is considered and if the last two
2794   argument (iterationNumber and orderNumber) are not set; their default value
2795   is -1. If the field fieldDriverName with the iteration number
2796   iterationNumber and the order number orderNumber does not exist in the file
2797   fieldDriverName; the constructor raises an exception.
2798 */
2799 template <class T, class INTERLACING_TAG>
2800 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
2801                                  driverTypes driverType,
2802                                  const string & fileName/*=""*/,
2803                                  const string & fieldDriverName/*=""*/,
2804                                  const int iterationNumber,
2805                                  const int orderNumber) throw (MEDEXCEPTION)
2806 {
2807   const char* LOC = "template <class T> FIELD<T>::FIELD(const SUPPORT * Support, driverTypes driverType, const string & fileName=\"\", const string & fieldName=\"\", const int iterationNumber=-1, const int orderNumber=-1) : ";
2808   BEGIN_OF_MED(LOC);
2809
2810   int current;
2811
2812   init();
2813
2814   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2815   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2816   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
2817
2818   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2819   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2820   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2821
2822   _support = Support;
2823   //A.G. Addings for RC
2824   if(_support)
2825     _support->addReference();
2826   // OCC 10/03/2006 -- According to the rules defined with help of
2827   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2828   // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2829   // FIELD template - MSVC++ 2003 compiler generated an error here.
2830   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2831   _value = NULL;
2832
2833   _iterationNumber = iterationNumber;
2834   _time = 0.0;
2835   _orderNumber = orderNumber;
2836
2837   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
2838
2839   _drivers[current]->open();
2840   _drivers[current]->read();
2841   _drivers[current]->close();
2842
2843   END_OF_MED(LOC);
2844 }
2845
2846 /*!
2847   If the mesh argument is not initialized or passed NULL,
2848   this constructor, at least, allows to create a FIELD without creating any
2849   SUPPORT then without having to load a MESH object, a support is created. It
2850   provides the meshName related mesh but doesn't not set a mesh in the created
2851   support.
2852   If the passed mesh contains corresponding support, this support will be used
2853   for the field. This support will be found in mesh by name of one of profiles,
2854   on which the FIELD lays in MED-file. This has sense for the case, then MED-file
2855   was created by MEDMEM, and so name of profile contains name of corresponding support.
2856 */
2857 template <class T, class INTERLACING_TAG>
2858 FIELD<T,INTERLACING_TAG>::FIELD(driverTypes driverType,
2859                                 const string & fileName,
2860                                 const string & fieldDriverName,
2861                                 const int iterationNumber,
2862                                 const int orderNumber,
2863                                 MESH* mesh)
2864   throw (MEDEXCEPTION) :FIELD_()
2865 {
2866   int current;
2867   const char* LOC = "FIELD<T,INTERLACING_TAG>::FIELD( driverTypes driverType, const string & fileName, string & fieldDriverName, int iterationNumber, int orderNumber) : ";
2868   BEGIN_OF_MED(LOC);
2869
2870   init();
2871
2872   _mesh = mesh;
2873
2874   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2875   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2876   FIELD_::_valueType = SET_VALUE_TYPE<T>::_valueType;
2877
2878   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2879   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2880   FIELD_::_interlacingType = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2881
2882   _support = (SUPPORT *) NULL;
2883   // OCC 10/03/2006 -- According to the rules defined with help of
2884   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2885   // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2886   // FIELD template - MSVC++ 2003 compiler generated an error here.
2887   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2888   _value = NULL;
2889
2890   _iterationNumber = iterationNumber;
2891   _time = 0.0;
2892   _orderNumber = orderNumber;
2893
2894   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
2895
2896   _drivers[current]->open();
2897   _drivers[current]->read();
2898   _drivers[current]->close();
2899
2900   END_OF_MED(LOC);
2901 }
2902 /*!
2903 @}
2904 */
2905
2906 /*!
2907   Destructor.
2908 */
2909 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::~FIELD()
2910 {
2911   const char* LOC = " Destructeur FIELD<T, INTERLACING_TAG>::~FIELD()";
2912   BEGIN_OF_MED(LOC);
2913   SCRUTE_MED(this);
2914   if (_value) delete _value;
2915   locMap::const_iterator it;
2916   for ( it = _gaussModel.begin();it != _gaussModel.end(); it++ )
2917     delete (*it).second;
2918
2919   END_OF_MED(LOC);
2920 }
2921
2922 /*!
2923
2924 */
2925 template <class T, class INTERLACING_TAG>
2926 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)
2927 {
2928   const char* LOC = "FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)";
2929   BEGIN_OF_MED(LOC);
2930
2931   _numberOfComponents = NumberOfComponents ;
2932   //if (_componentsTypes == NULL)
2933   //  _componentsTypes = new int[NumberOfComponents] ;
2934   _componentsTypes.resize(NumberOfComponents);
2935   //if (_componentsNames == NULL)
2936   //  _componentsNames = new string[NumberOfComponents];
2937   _componentsNames.resize(NumberOfComponents);
2938   //if (_componentsDescriptions == NULL)
2939   //  _componentsDescriptions = new string[NumberOfComponents];
2940   _componentsDescriptions.resize(NumberOfComponents);
2941   //if (_componentsUnits == NULL)
2942   //  _componentsUnits = new UNIT[NumberOfComponents];
2943   _componentsUnits.resize(NumberOfComponents);
2944   //if (_MEDComponentsUnits == NULL)
2945   //  _MEDComponentsUnits = new string[NumberOfComponents];
2946   _MEDComponentsUnits.resize(NumberOfComponents);
2947   for (int i=0;i<NumberOfComponents;i++) {
2948     _componentsTypes[i] = 0 ;
2949   }
2950
2951   try {
2952     // becarefull about the number of gauss point
2953     _numberOfValues = _support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2954     MESSAGE_MED(PREFIX_MED <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
2955
2956     //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
2957     _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
2958
2959     _isRead = true ;
2960   }
2961   catch (MEDEXCEPTION &ex) {
2962     MESSAGE_MED("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
2963     // OCC 10/03/2006 -- According to the rules defined with help of
2964     // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2965     // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2966     // FIELD template - MSVC++ 2003 compiler generated an error here.
2967     // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2968     _value = NULL;
2969   }
2970
2971   SCRUTE_MED(_value);
2972   END_OF_MED(LOC);
2973 }
2974
2975 /*!
2976
2977 */
2978 template <class T, class INTERLACING_TAG>
2979 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents,
2980                                            const int LengthValue)
2981 {
2982   const char* LOC = "void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)";
2983   BEGIN_OF_MED(LOC);
2984
2985   _numberOfComponents = NumberOfComponents ;
2986   //if (_componentsTypes == NULL)
2987   //  _componentsTypes = new int[NumberOfComponents] ;
2988   //if (_componentsNames == NULL)
2989   //  _componentsNames = new string[NumberOfComponents];
2990   //if (_componentsDescriptions == NULL)
2991   //  _componentsDescriptions = new string[NumberOfComponents];
2992   //if (_componentsUnits == NULL)
2993   //  _componentsUnits = new UNIT[NumberOfComponents];
2994   //if (_MEDComponentsUnits == NULL)
2995   //  _MEDComponentsUnits = new string[NumberOfComponents];
2996   _componentsTypes.resize(NumberOfComponents);
2997   _componentsNames.resize(NumberOfComponents);
2998   _componentsDescriptions.resize(NumberOfComponents);
2999   _componentsUnits.resize(NumberOfComponents);
3000   _MEDComponentsUnits.resize(NumberOfComponents);
3001   for (int i=0;i<NumberOfComponents;i++) {
3002     _componentsTypes[i] = 0 ;
3003   }
3004
3005   MESSAGE_MED("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
3006   _numberOfValues = LengthValue ;
3007
3008   //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
3009   _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
3010
3011   _isRead = true ;
3012
3013   SCRUTE_MED(_value);
3014   END_OF_MED(LOC);
3015 }
3016
3017 /*!
3018
3019 */
3020 template <class T, class INTERLACING_TAG>
3021 void FIELD<T, INTERLACING_TAG>::deallocValue()
3022 {
3023   const char* LOC = "void FIELD<T, INTERLACING_TAG>::deallocValue()";
3024   BEGIN_OF_MED(LOC);
3025   _numberOfValues = 0 ;
3026   _numberOfComponents = 0 ;
3027   if (_value != NULL) {
3028     delete _value;
3029     _value = NULL;
3030   }
3031
3032   END_OF_MED(LOC);
3033 }
3034
3035
3036
3037 /*!\if MEDMEM_ug
3038
3039 \addtogroup FIELD_io
3040 @{
3041 \endif */
3042 // -----------------
3043 // Methodes Inline
3044 // -----------------
3045
3046 /*!
3047   Creates the specified driver and return its index reference to path to
3048   read or write methods.
3049
3050 \param driverType specifies the file type (MED_DRIVER or VTK_DRIVER)
3051 \param fileName name of the output file
3052 \param driverName name of the field
3053 \param access access type (read, write or both)
3054
3055 */
3056
3057 template <class T, class INTERLACING_TAG>
3058 int FIELD<T, INTERLACING_TAG>::addDriver(driverTypes driverType,
3059                                          const string & fileName/*="Default File Name.med"*/,
3060                                          const string & driverName/*="Default Field Name"*/,
3061                                          MED_EN::med_mode_acces access)
3062 {
3063   //jfa tmp (as last argument has no default value):const char * LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Field Name\",MED_EN::med_mode_acces access) : ";
3064
3065   GENDRIVER * driver;
3066
3067   const char* LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName,const string & driverName,MED_EN::med_mode_acces access) :";
3068   BEGIN_OF_MED(LOC);
3069
3070   SCRUTE_MED(driverType);
3071
3072   driver = DRIVERFACTORY::buildDriverForField(driverType,fileName,this,access);
3073
3074   _drivers.push_back(driver);
3075
3076   int current = _drivers.size()-1;
3077
3078   _drivers[current]->setFieldName(driverName);
3079
3080   END_OF_MED(LOC);
3081
3082   return current;
3083 }
3084 /*! \if MEDMEM_ug @} \endif */
3085
3086 /*!
3087   Duplicates the given driver and return its index reference to path to
3088   read or write methods.
3089 */
3090 template <class T, class INTERLACING_TAG>
3091 inline int FIELD<T, INTERLACING_TAG>::addDriver (GENDRIVER & driver )
3092 {
3093   int current;
3094
3095   const char* LOC = "FIELD<T, INTERLACING_TAG>::addDriver(GENDRIVER &) : ";
3096   BEGIN_OF_MED(LOC);
3097
3098   // duplicate driver to delete it with destructor !
3099   //GENDRIVER * newDriver = driver.copy() ;
3100
3101   // for FIELD->read( genDriver ) if FIELD was not passed to genDriver
3102   GENDRIVER * newDriver =
3103     DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
3104                                        driver.getFileName(), this,
3105                                        driver.getAccessMode());
3106   _drivers.push_back(newDriver);
3107
3108   current = _drivers.size()-1;
3109   SCRUTE_MED(current);
3110   driver.setId(current);
3111
3112   newDriver->merge( driver );
3113   newDriver->setId(current);
3114
3115   return current ;
3116 }
3117
3118 /*!
3119   Remove the driver referenced by its index.
3120 */
3121 template <class T, class INTERLACING_TAG>
3122 void FIELD<T, INTERLACING_TAG>::rmDriver (int index/*=0*/)
3123 {
3124   const char * LOC = "FIELD<T, INTERLACING_TAG>::rmDriver (int index=0): ";
3125   BEGIN_OF_MED(LOC);
3126
3127   if ( index>=0 && index<_drivers.size() && _drivers[index] ) {
3128     //_drivers.erase(&_drivers[index]);
3129     // why not ????
3130     MESSAGE_MED ("detruire");
3131   }
3132   else
3133     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3134                                      << "The <index given is invalid, index must be between  0 and  |"
3135                                      << _drivers.size()
3136                                      )
3137                           );
3138
3139   END_OF_MED(LOC);
3140 }
3141
3142 /*!
3143   Read FIELD in the file specified in the driver given by its index.
3144 */
3145 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(int index/*=0*/)
3146 {
3147   const char * LOC = "FIELD<T, INTERLACING_TAG>::read(int index=0) : ";
3148   BEGIN_OF_MED(LOC);
3149
3150   if ( index>=0 && index<_drivers.size() && _drivers[index] ) {
3151     _drivers[index]->open();
3152     _drivers[index]->read();
3153     _drivers[index]->close();
3154   }
3155   else
3156     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3157                                      << "The index given is invalid, index must be between  0 and |"
3158                                      << _drivers.size()
3159                                      )
3160                           );
3161   END_OF_MED(LOC);
3162 }
3163
3164 /*! \if MEDMEM_ug
3165 \addtogroup FIELD_io
3166 @{
3167 \endif */
3168
3169 /*!
3170   Writes FIELD in the file specified by the driver handle \a index.
3171
3172 Example :
3173 \verbatim
3174 //...
3175 // Attaching the friver to file "output.med", meshname "Mesh"
3176 int driver_handle = mesh.addDriver(MED_DRIVER, "output.med", "Mesh");
3177 // Writing the content of mesh to the file
3178 mesh.write(driver_handle);
3179 \endverbatim
3180 */
3181 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(int index/*=0*/, const string & driverName /*= ""*/)
3182 {
3183   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
3184   BEGIN_OF_MED(LOC);
3185
3186   if( index>=0 && index<_drivers.size() && _drivers[index] ) {
3187     _drivers[index]->open();
3188     if (driverName != "") _drivers[index]->setFieldName(driverName);
3189     _drivers[index]->write();
3190     _drivers[index]->close();
3191   }
3192   else
3193     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3194                                      << "The index given is invalid, index must be between  0 and |"
3195                                      << _drivers.size()
3196                                      )
3197                           );
3198   END_OF_MED(LOC);
3199 }
3200 /*! \if MEDMEM_ug @} \endif */
3201 /*!
3202   Write FIELD in the file specified in the driver given by its index. Use this
3203   method for ASCII drivers (e.g. VTK_DRIVER)
3204 */
3205 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
3206 {
3207   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
3208   BEGIN_OF_MED(LOC);
3209
3210   if( index>=0 && index<_drivers.size() && _drivers[index] ) {
3211     _drivers[index]->openAppend();
3212     if (driverName != "") _drivers[index]->setFieldName(driverName);
3213     _drivers[index]->writeAppend();
3214     _drivers[index]->close();
3215   }
3216   else
3217     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3218                                      << "The index given is invalid, index must be between  0 and |"
3219                                      << _drivers.size()
3220                                      )
3221                           );
3222   END_OF_MED(LOC);
3223 }
3224
3225 /*!
3226   \internal
3227   Write FIELD with the driver which is equal to the given driver.
3228
3229   Use by MED object.
3230 */
3231 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(const GENDRIVER & genDriver)
3232 {
3233   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
3234   BEGIN_OF_MED(LOC);
3235
3236   for (unsigned int index=0; index < _drivers.size(); index++ )
3237     if ( *_drivers[index] == genDriver ) {
3238       _drivers[index]->open();
3239       _drivers[index]->write();
3240       _drivers[index]->close();
3241     }
3242
3243   END_OF_MED(LOC);
3244
3245 }
3246
3247 /*!
3248   \internal
3249   Write FIELD with the driver which is equal to the given driver.
3250
3251   Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
3252 */
3253 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(const GENDRIVER & genDriver)
3254 {
3255   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
3256   BEGIN_OF_MED(LOC);
3257
3258   for (unsigned int index=0; index < _drivers.size(); index++ )
3259     if ( *_drivers[index] == genDriver ) {
3260       _drivers[index]->openAppend();
3261       _drivers[index]->writeAppend();
3262       _drivers[index]->close();
3263     }
3264
3265   END_OF_MED(LOC);
3266
3267 }
3268
3269 /*!
3270   \internal
3271   Read FIELD with the driver which is equal to the given driver.
3272
3273   Use by MED object.
3274 */
3275 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(const GENDRIVER & genDriver)
3276 {
3277   const char* LOC = " FIELD<T, INTERLACING_TAG>::read(const GENDRIVER &) : ";
3278   BEGIN_OF_MED(LOC);
3279
3280   for (unsigned int index=0; index < _drivers.size(); index++ )
3281     if ( *_drivers[index] == genDriver ) {
3282       _drivers[index]->open();
3283       _drivers[index]->read();
3284       _drivers[index]->close();
3285     }
3286
3287   END_OF_MED(LOC);
3288
3289 }
3290
3291 /*!
3292   Fills in already allocated retValues array the values related to eltIdInSup.
3293   If the element does not exist in this->_support false is returned, true otherwise.
3294  */
3295 template <class T, class INTERLACING_TAG>
3296 bool FIELD<T, INTERLACING_TAG>::getValueOnElement(int eltIdInSup,T* retValues)
3297   const throw (MEDEXCEPTION)
3298 {
3299
3300   if(eltIdInSup<1)
3301     return false;
3302   if(_support->isOnAllElements())
3303     {
3304       int nbOfEltsThis=_support->getMesh()->getNumberOfElements(_support->getEntity(),MED_EN::MED_ALL_ELEMENTS);
3305       if(eltIdInSup>nbOfEltsThis)
3306         return false;
3307       const T* valsThis=getValue();
3308       for(int j=0;j<_numberOfComponents;j++)
3309         retValues[j]=valsThis[(eltIdInSup-1)*_numberOfComponents+j];
3310       return true;
3311     }
3312   else
3313     {
3314       int nbOfEltsThis=_support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
3315       const int *eltsThis=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
3316       int iThis;
3317       bool found=false;
3318       for(iThis=0;iThis<nbOfEltsThis && !found;)
3319         if(eltsThis[iThis]==eltIdInSup)
3320           found=true;
3321         else
3322           iThis++;
3323       if(!found)
3324         return false;
3325       const T* valsThis=getValue();
3326       for(int j=0;j<_numberOfComponents;j++)
3327         retValues[j]=valsThis[iThis*_numberOfComponents+j];
3328       return true;
3329     }
3330 }
3331
3332 /*!
3333   \if developper
3334   Destroy the MEDARRAY<T> in FIELD and put the new one without copy.
3335   \endif
3336 */
3337 template <class T, class INTERLACING_TAG>
3338 inline void FIELD<T, INTERLACING_TAG>::setArray(MEDMEM_Array_ * Value)
3339   throw (MEDEXCEPTION)
3340 {
3341   if (NULL != _value) delete _value ;
3342   _value=Value ;
3343 }
3344
3345 /*!
3346   \if developper
3347   Return a reference to  the MEDARRAY<T, INTERLACING_TAG> in FIELD.
3348   \endif
3349 */
3350 template <class T, class INTERLACING_TAG>
3351 inline MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() const throw (MEDEXCEPTION)
3352 {
3353   const char* LOC = "MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() : ";
3354   BEGIN_OF_MED(LOC);
3355   END_OF_MED(LOC);
3356   return _value ;
3357 }
3358 template <class T,class INTERLACING_TAG>  inline
3359 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array *
3360 FIELD<T, INTERLACING_TAG>::getArrayGauss() const throw (MEDEXCEPTION)
3361 {
3362   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayGauss() : ";
3363   BEGIN_OF_MED(LOC);
3364
3365   if ( getGaussPresence() )
3366     return static_cast<ArrayGauss *> (_value);
3367   else
3368     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
3369                                  "The field has no Gauss Point"));
3370
3371   END_OF_MED(LOC);
3372
3373 }
3374
3375 template <class T,class INTERLACING_TAG>  inline
3376 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array *
3377 FIELD<T, INTERLACING_TAG>::getArrayNoGauss() const throw (MEDEXCEPTION)
3378 {
3379   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayNoGauss() : ";
3380   BEGIN_OF_MED(LOC);
3381
3382   if ( ! getGaussPresence() )
3383     return static_cast < ArrayNoGauss * > (_value);
3384   else
3385     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
3386                                  "The field has Gauss Point"));
3387
3388   END_OF_MED(LOC);
3389 }
3390
3391
3392 template <class T,class INTERLACING_TAG> inline bool
3393 FIELD<T, INTERLACING_TAG>::getGaussPresence() const throw (MEDEXCEPTION)
3394 {
3395   const char * LOC = "FIELD<T, INTERLACING_TAG>::getGaussPresence() const :";
3396   //BEGIN_OF_MED(LOC);
3397
3398   if (_value != NULL)
3399     return _value->getGaussPresence();
3400   else
3401     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't call getGaussPresence on a null _value"));
3402
3403   //END_OF_MED();
3404 }
3405
3406 /*!
3407   Return the actual length of the reference to values array returned by getValue.
3408   Take care of number of components and number of Gauss points by geometric type
3409 */
3410 template <class T, class INTERLACING_TAG>
3411 inline int FIELD<T, INTERLACING_TAG>::getValueLength() const
3412   throw (MEDEXCEPTION)
3413 {
3414   if ( getGaussPresence() )
3415     return static_cast<ArrayGauss *>(_value)->getArraySize() ;
3416   else
3417     return static_cast<ArrayNoGauss *>(_value)->getArraySize() ;
3418 }
3419
3420
3421 /*! \if MEDMEM_ug
3422 \defgroup FIELD_value Field values
3423
3424 These methods are provided for accessing the values of a field.
3425 There are two ways to do so : one consists in using accessors
3426 that retrieve elements or group of elements from the entire field. Typical use is
3427 \verbatim
3428 FIELD field(MED_DRIVER, "result.med","Pressure");
3429 double P0=field.getValueIJ(1,1);
3430 \endverbatim
3431
3432 Another way is to retrieve the pointer to the array that contains the
3433 variable values. In this case, the user should be aware of the interlacing mode
3434 so that no mistakes are made when retrieving the values.
3435
3436 \verbatim
3437 FIELD field(MED_DRIVER, "result.med","Pressure");
3438 double* ptrP=field.getValue();
3439 double P0=ptrP[0];
3440 \endverbatim
3441
3442 @{
3443 \endif
3444 */
3445 /*!
3446   Returns a reference to values array to read them.
3447 */
3448 template <class T, class INTERLACIN_TAG>
3449 inline const T* FIELD<T, INTERLACIN_TAG>::getValue() const throw (MEDEXCEPTION)
3450 {
3451   const char* LOC = "FIELD<T, INTERLACING_TAG>::getValue() : ";
3452   BEGIN_OF_MED(LOC);
3453   if ( getGaussPresence() )
3454     return static_cast<ArrayGauss *>(_value)->getPtr() ;
3455   else
3456     return static_cast<ArrayNoGauss *>(_value)->getPtr() ;
3457 }
3458 /*!
3459   Returns a reference to \f$ i^{th} \f$ row
3460   of FIELD values array.
3461   If a faster accessor is intended you may use getArray() once,
3462   then MEDMEM_Array accessors.
3463   Be careful if field support is not on all elements getRow
3464   use support->getValIndFromGlobalNumber(i).
3465 */
3466 template <class T,class INTERLACING_TAG> inline
3467 const T*
3468 FIELD<T,INTERLACING_TAG>::getRow(int i) const throw (MEDEXCEPTION)
3469 {
3470   const char * LOC = "FIELD<T,INTERLACING_TAG>::getRow(int i) : ";
3471   //BEGIN_OF_MED(LOC);
3472
3473   int valIndex=-1;
3474   if (_support)
3475     valIndex = _support->getValIndFromGlobalNumber(i);
3476   else
3477     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
3478
3479   //cout << endl << "getRow Valindex : " << valIndex << endl;
3480   if ( getGaussPresence() )
3481     return static_cast<ArrayGauss *>(_value)->getRow(valIndex) ;
3482   else
3483     return static_cast<ArrayNoGauss *>(_value)->getRow(valIndex) ;
3484   //END_OF_MED();
3485 }
3486
3487 /*!
3488   Returns a reference to $j^{th}$ column
3489   of FIELD values array.
3490 */
3491 template <class T,class INTERLACING_TAG> inline const T*
3492 FIELD<T,INTERLACING_TAG>::getColumn(int j) const throw (MEDEXCEPTION)
3493 {
3494   //BEGIN_OF_MED(LOC);
3495   if ( getGaussPresence() )
3496     return static_cast<ArrayGauss *>(_value)->getColumn(j) ;
3497   else
3498     return static_cast<ArrayNoGauss *>(_value)->getColumn(j) ;
3499 }
3500
3501 /*!
3502   Returns the value of $i^{th}$ element and $j^{th}$ component.
3503   This method only works with fields having no particular Gauss point
3504 definition (i.e., fields having one value per element).
3505  This method makes the retrieval of the value independent from the
3506   interlacing pattern, but it is slower than the complete retrieval
3507   obtained by the \b getValue() method.
3508 */
3509 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJ(int i,int j) const throw (MEDEXCEPTION)
3510 {
3511   const char * LOC = "getValueIJ(..)";
3512   //BEGIN_OF_MED(LOC);
3513   int valIndex=-1;
3514   if (_support)
3515     valIndex = _support->getValIndFromGlobalNumber(i);
3516   else
3517     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3518
3519   if ( getGaussPresence() )
3520     return static_cast<ArrayGauss *>(_value)->getIJ(valIndex,j) ;
3521   else
3522     return static_cast<ArrayNoGauss *>(_value)->getIJ(valIndex,j) ;
3523 }
3524
3525 /*!
3526   Returns the $j^{th}$ component of $k^{th}$ Gauss points of $i^{th}$ value.
3527   This method is compatible with elements having more than one Gauss point.
3528   This method makes the retrieval of the value independent from the
3529   interlacing pattern, but it is slower than the complete retrieval
3530   obtained by the \b getValue() method.
3531 */
3532 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION)
3533 {
3534   const char * LOC = "getValueIJK(..)";
3535   //BEGIN_OF_MED(LOC);
3536   int valIndex=-1;
3537   if (_support)
3538     valIndex = _support->getValIndFromGlobalNumber(i);
3539   else
3540     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3541
3542   if ( getGaussPresence() )
3543     return static_cast<ArrayGauss *>(_value)->getIJK(valIndex,j,k) ;
3544   else
3545     return static_cast<ArrayNoGauss *>(_value)->getIJK(valIndex,j,k) ;
3546 }
3547 /*! \if MEDMEM_ug @} \endif */
3548
3549 /*!
3550   Return number of values of a geomertic type in NoInterlaceByType mode
3551 */
3552 template <class T, class INTERLACIN_TAG>
3553 inline int FIELD<T, INTERLACIN_TAG>::getValueByTypeLength(int t) const throw (MEDEXCEPTION)
3554 {
3555   const char * LOC ="getValueByTypeLength() : ";
3556   //BEGIN_OF_MED(LOC);
3557   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3558     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3559
3560   if ( getGaussPresence() ) {
3561     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
3562     if (  t < 1 || t > array->getNbGeoType() )
3563       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
3564     return array->getLengthOfType( t );
3565   }
3566   else {
3567     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
3568     if (  t < 1 || t > array->getNbGeoType() )
3569       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
3570     return array->getLengthOfType( t );
3571   }
3572 }
3573
3574 /*!
3575   Return a reference to values array to read them.
3576 */
3577 template <class T, class INTERLACIN_TAG>
3578 inline const T* FIELD<T, INTERLACIN_TAG>::getValueByType(int t) const throw (MEDEXCEPTION)
3579 {
3580   const char * LOC ="getValueByType() : ";
3581   //BEGIN_OF_MED(LOC);
3582   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3583     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3584
3585   if ( getGaussPresence() ) {
3586     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
3587     return array->getPtr() + array->getIndex( t );
3588   }
3589   else {
3590     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
3591     return array->getPtr() + array->getIndex( t );
3592   }
3593 }
3594
3595 /*!
3596   Return the value of i^{th} element in indicated type t and j^{th} component.
3597 */
3598 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJByType(int i,int j, int t) const throw (MEDEXCEPTION)
3599 {
3600   const char * LOC = "getValueIJByType(..)";
3601   //BEGIN_OF_MED(LOC);
3602   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3603     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3604
3605   if ( getGaussPresence() )
3606     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJByType(i,j,t) ;
3607   else
3608     return static_cast<ArrayNoByType *>(_value)->getIJByType(i,j,t) ;
3609 }
3610
3611 /*!
3612   Return the j^{th} component of k^{th} gauss points of i^{th} value with type t.
3613 */
3614 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION)
3615 {
3616   const char * LOC = "getValueIJKByType(..)";
3617   //BEGIN_OF_MED(LOC);
3618   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3619     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3620
3621   if ( getGaussPresence() )
3622     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJKByType(i,j,k,t) ;
3623   else
3624     return static_cast<ArrayNoByType      *>(_value)->getIJKByType(i,j,k,t) ;
3625 }
3626
3627
3628 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGeometricTypes() const throw (MEDEXCEPTION)
3629 {
3630   const char * LOC = "getNumberOfGeometricTypes(..)";
3631   BEGIN_OF_MED(LOC);
3632   if (_support)
3633     return _support->getNumberOfTypes();
3634   else
3635     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3636   END_OF_MED(LOC);
3637 }
3638
3639 /*! \if MEDMEM_ug
3640 \addtogroup FIELD_gauss
3641 @{
3642 \endif
3643  */
3644
3645 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> &
3646 FIELD<T,INTERLACING_TAG>::getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
3647 {
3648   const char * LOC ="getGaussLocalization(MED_EN::medGeometryElement geomElement) : ";
3649   const GAUSS_LOCALIZATION_ * locPtr=0;
3650
3651   locMap::const_iterator it;
3652   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3653         locPtr = (*it).second;
3654         return  *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
3655   }
3656   else
3657     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
3658
3659 }
3660
3661 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> *
3662 FIELD<T,INTERLACING_TAG>::getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
3663 {
3664   const char * LOC ="getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) : ";
3665   const GAUSS_LOCALIZATION_ * locPtr=0;
3666
3667   locMap::const_iterator it;
3668   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3669         locPtr = (*it).second;
3670         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
3671   }
3672   else
3673     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
3674
3675 }
3676
3677 /*!
3678  * \brief Return GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
3679  */
3680 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION_ *
3681 FIELD<T,INTERLACING_TAG>::getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
3682 {
3683   const char * LOC ="getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) : ";
3684
3685   locMap::const_iterator it;
3686   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3687         return (*it).second;
3688   }
3689   else
3690     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type: "<< geomElement ));
3691
3692 }
3693
3694 /*!
3695  * \brief Take onership of GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
3696  */
3697 template <class T,class INTERLACING_TAG> void
3698 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc)
3699 {
3700   locMap::iterator it = _gaussModel.find(geomElement);
3701   if ( it != _gaussModel.end() ) {
3702     delete it->second;
3703     it->second = gaussloc;
3704   }
3705   else {
3706     _gaussModel[ geomElement ] = gaussloc;
3707   }
3708 }
3709
3710
3711 template <class T,class INTERLACING_TAG> void
3712 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG>& gaussloc)
3713 {
3714   locMap::iterator it = _gaussModel.find(geomElement);
3715   if ( it != _gaussModel.end() ) {
3716     delete it->second;
3717     it->second = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
3718   }
3719   else {
3720     _gaussModel[ geomElement ] = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
3721   }
3722 }
3723
3724 /*!
3725   Returns number of Gauss points for this medGeometryElement.
3726
3727   Note :
3728   if there is no GAUSS_LOCALIZATION having this medGeometryElement but
3729   the medGeometryElement exist in the SUPPORT, getNumberOfGaussPoints
3730   return 1
3731 */
3732 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) const
3733   throw (MEDEXCEPTION)
3734 {
3735   const char * LOC ="getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
3736   const GAUSS_LOCALIZATION_ * locPtr=0;
3737
3738   locMap::const_iterator it;
3739   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3740         locPtr = (*it).second;
3741         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr)->getNbGauss();
3742   }
3743   else
3744     if (_support)
3745       try {
3746         if ( _support->getNumberOfElements(geomElement) ) return 1;
3747       } catch ( MEDEXCEPTION & ex) {
3748         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<< "GeometricType not found !" )) ;
3749       }
3750     else
3751       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3752
3753   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Should never execute this!" ));
3754
3755 }
3756
3757 /*!
3758   Returns number of Gauss points for each geometric type.
3759
3760   Note :
3761   if there is no gauss points whatever the geometric type is
3762   it returns an exception. (renvoyer un tableau de 1 ?)
3763 */
3764 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints() const
3765   throw (MEDEXCEPTION)
3766 {
3767   const char * LOC ="const int * getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
3768
3769   if (_value)
3770     if ( getGaussPresence() ) {
3771       return static_cast<ArrayGauss *>(_value)->getNbGaussGeo()+1;
3772     } else
3773       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"value hasn't Gauss points " ));
3774
3775     else
3776       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Value not defined" ));
3777 }
3778
3779 /*!
3780   Returns number of Gauss points for element n°i.
3781   The i index is a global index (take care of previous element
3782   on different geometric type).
3783 */
3784 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNbGaussI(int i) const throw (MEDEXCEPTION)
3785 {
3786   const char * LOC = "getNbGaussI(..)";
3787 //   BEGIN_OF_MED(LOC);
3788
3789   int valIndex=-1;
3790   if (_support)
3791     valIndex = _support->getValIndFromGlobalNumber(i);
3792   else
3793     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3794
3795   if (_value)
3796    if ( getGaussPresence() )
3797      return static_cast<ArrayGauss *>(_value)->getNbGauss(valIndex) ;
3798    else
3799      return static_cast<ArrayNoGauss *>(_value)->getNbGauss(valIndex) ;
3800  else
3801    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_value not defined" ));
3802 //   END_OF_MED();
3803 }
3804 /*!
3805 @}
3806  */
3807
3808
3809 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfElements() const throw (MEDEXCEPTION)
3810 {
3811   const char * LOC = "getNumberOfElements(..)";
3812   BEGIN_OF_MED(LOC);
3813   if (_support)
3814     return _support->getNumberOfElements();
3815   else
3816     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3817   END_OF_MED(LOC);
3818 }
3819
3820 template <class T,class INTERLACING_TAG> const MED_EN::medGeometryElement  * FIELD<T,INTERLACING_TAG>::getGeometricTypes()  const throw (MEDEXCEPTION)
3821 {
3822   const char * LOC = "getGeometricTypes(..)";
3823   BEGIN_OF_MED(LOC);
3824   if (_support)
3825     return _support->getTypes();
3826   else
3827     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3828   END_OF_MED(LOC);
3829 }
3830 template <class T,class INTERLACING_TAG> bool  FIELD<T,INTERLACING_TAG>::isOnAllElements() const throw (MEDEXCEPTION)
3831 {
3832   const char * LOC = "isOnAllElements(..)";
3833   BEGIN_OF_MED(LOC);
3834   if (_support)
3835     return _support->isOnAllElements();
3836   else
3837     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3838   END_OF_MED(LOC);
3839 }
3840
3841
3842 /*! \if MEDMEM_ug
3843 \addtogroup FIELD_value
3844 @{
3845 \endif */
3846
3847 /*!
3848   Copy new values array in FIELD according to the given mode.
3849
3850   Array must have right size. If not results are unpredicable.
3851   In MED_FULL_INTERLACE mode, values are stored elementwise in X1,Y1,Z1,X2,Y2,Z2.. order.
3852 In MED_NO_INTERLACE mode, values are stored componentwise in X1,X2,X3,...,Y1,Y2,Y3,... order.
3853 */
3854 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValue( T* value) throw (MEDEXCEPTION)
3855 {
3856   if ( getGaussPresence() )
3857     static_cast<ArrayGauss *>(_value)->setPtr(value) ;
3858   else
3859     static_cast<ArrayNoGauss *>(_value)->setPtr(value) ;
3860 }
3861
3862 /*!
3863   Update values array in the j^{th} row of FIELD values array with the given ones and
3864   according to specified mode.
3865 */
3866 template <class T,class INTERLACING_TAG>
3867 inline void FIELD<T,INTERLACING_TAG>::setRow( int i, T* value) throw (MEDEXCEPTION)
3868 {
3869   const char * LOC = "FIELD<T,INTERLACING_TAG>::setRow(int i, T* value) : ";
3870   int valIndex=i;
3871   if (_support)
3872     valIndex = _support->getValIndFromGlobalNumber(i);
3873   else
3874     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
3875
3876   if ( getGaussPresence() )
3877     static_cast<ArrayGauss *>(_value)->setRow(valIndex, value) ;
3878   else
3879     static_cast<ArrayNoGauss *>(_value)->setRow(valIndex, value) ;
3880 }
3881
3882 /*!
3883   Update values array in the $j^{th}$ column of FIELD values array with the given ones and
3884   according to specified mode.
3885 */
3886 template <class T,class INTERLACING_TAG>
3887 inline void FIELD<T,INTERLACING_TAG>::setColumn( int j, T* value) throw (MEDEXCEPTION)
3888 {
3889   if ( getGaussPresence() )
3890     static_cast<ArrayGauss *>(_value)->setColumn(j, value) ;
3891   else
3892     static_cast<ArrayNoGauss *>(_value)->setColumn(j, value) ;
3893 }
3894
3895 /*!
3896   Sets the value of i^{th} element and j^{th} component with the given one.
3897 */
3898 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) throw (MEDEXCEPTION)
3899 {
3900   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
3901   int valIndex=-1;
3902   if (_support)
3903     valIndex = _support->getValIndFromGlobalNumber(i);
3904   else
3905     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
3906
3907   if ( getGaussPresence() )
3908     static_cast<ArrayGauss *>(_value)->setIJ(valIndex,j,value) ;
3909   else
3910     static_cast<ArrayNoGauss *>(_value)->setIJ(valIndex,j,value) ;
3911 }
3912
3913 /*!
3914   Set the value of i^{th} element, j^{th} component and k^{th} gauss point with the given one.
3915 */
3916 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION)
3917 {
3918   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
3919   int valIndex=-1;
3920   if (_support)
3921     valIndex = _support->getValIndFromGlobalNumber(i);
3922   else
3923     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
3924
3925   if ( getGaussPresence() )
3926     static_cast<ArrayGauss *>(_value)->setIJK(valIndex,j,k,value) ;
3927   else
3928     static_cast<ArrayNoGauss *>(_value)->setIJK(valIndex,j,k,value) ;
3929 }
3930
3931 /*!
3932   Set the value of i^{th} element and j^{th} component with the given one.
3933 */
3934 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION)
3935 {
3936   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) : ";
3937   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3938     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3939
3940   if ( getGaussPresence() )
3941     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJByType(i,j,t,value) ;
3942   else
3943     return static_cast<ArrayNoByType *>(_value)->setIJByType(i,j,t,value) ;
3944 }
3945
3946 /*!
3947   Set the value of component of k^{th} gauss points of i^{th} element and j^{th} component with the given one.
3948 */
3949 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION)
3950 {
3951   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int t, int k, T value) : ";
3952   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3953     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3954
3955   if ( getGaussPresence() )
3956     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJKByType(i,j,k,t,value) ;
3957   else
3958     return static_cast<ArrayNoByType *>(_value)->setIJKByType(i,j,k,t,value) ;
3959 }
3960 /*! \if MEDMEM_ug @} \endif */
3961
3962 /*
3963   METHODS
3964 */
3965
3966 /*!
3967   Fill values array with volume values.
3968 */
3969 // template <class T, class INTERLACING_TAG>
3970 // void FIELD<T, INTERLACING_TAG>::getVolume() const throw (MEDEXCEPTION)
3971 // {
3972 //   const char * LOC = "FIELD<double>::getVolume() const : ";
3973 //   BEGIN_OF_MED(LOC);
3974
3975 //   // The field has to be initilised by a non empty support and a
3976 //   // number of components = 1 and its value type has to be set to MED_REEL64
3977 //   // (ie a FIELD<double>)
3978
3979 //   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
3980 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to 1 and a value type set to MED_REEL64"));
3981
3982 //   END_OF_MED(LOC);
3983 // }
3984
3985 // /*!
3986 //   Fill values array with area values.
3987 // */
3988 // template <class T, class INTERLACING_TAG>
3989 // void FIELD<T, INTERLACING_TAG>::getArea() const throw (MEDEXCEPTION)
3990 // {
3991 //   const char * LOC = "FIELD<double>::getArea() const : ";
3992 //   BEGIN_OF_MED(LOC);
3993
3994 //   // The field has to be initilised by a non empty support and a
3995 //   // number of components = 1 and its value type has to be set to MED_REEL64
3996 //   // (ie a FIELD<double>)
3997
3998 //   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
3999 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to 1 and a value type set to MED_REEL64"));
4000
4001 //   END_OF_MED(LOC);
4002 // }
4003
4004 // /*!
4005 //   Fill values array with length values.
4006 // */
4007 // template <class T, class INTERLACING_TAG>
4008 // void FIELD<T, INTERLACING_TAG>::getLength() const throw (MEDEXCEPTION)
4009 // {
4010 //   const char * LOC = "FIELD<double>::getLength() const : ";
4011 //   BEGIN_OF_MED(LOC);
4012
4013 //   // The field has to be initilised by a non empty support and a
4014 //   // number of components = 1 and its value type has to be set to MED_REEL64
4015 //   // (ie a FIELD<double>)
4016
4017 //   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
4018 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to 1 and a value type set to MED_REEL64"));
4019
4020 //   END_OF_MED(LOC);
4021 // }
4022
4023 // /*!
4024 //   Fill values array with normal values.
4025 // */
4026 // template <class T, class INTERLACING_TAG>
4027 // void FIELD<T, INTERLACING_TAG>::getNormal() const throw (MEDEXCEPTION)
4028 // {
4029 //   const char * LOC = "FIELD<double>::getNormal() const : ";
4030 //   BEGIN_OF_MED(LOC);
4031
4032 //   // The field has to be initilised by a non empty support and a
4033 //   // number of components = 1 and its value type has to be set to MED_REEL64
4034 //   // (ie a FIELD<double>)
4035
4036 //   if (_support == (SUPPORT *) NULL)
4037 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
4038
4039 //   int dim_space = _support->getMesh()->getSpaceDimension();
4040
4041 //   if ((_numberOfComponents != dim_space) || (_valueType != MED_EN::MED_REEL64))
4042 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
4043
4044 //   END_OF_MED(LOC);
4045 // }
4046
4047 // /*!
4048 //   Fill values array with barycenter values.
4049 // */
4050 // template <class T, class INTERLACING_TAG>
4051 // void FIELD<T, INTERLACING_TAG>::getBarycenter() const throw (MEDEXCEPTION)
4052 // {
4053 //   const char * LOC = "FIELD<double>::getBarycenter() const : ";
4054 //   BEGIN_OF_MED(LOC);
4055
4056 //   // The field has to be initilised by a non empty support and a number of
4057 //   //components = space dimension and its value type has to be set to MED_REEL64
4058 //   // (ie a FIELD<double>)
4059
4060 //   if (_support == (SUPPORT *) NULL)
4061 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
4062
4063 //   int dim_space = _support->getMesh()->getSpaceDimension();
4064
4065 //   if ((_numberOfComponents != dim_space) || (_valueType != MED_EN::MED_REEL64))
4066 //       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
4067
4068 //   END_OF_MED(LOC);
4069 // }
4070
4071 /*!
4072   Fill array by using T_Analytic.
4073   WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
4074   Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
4075  */
4076 template <class T, class INTERLACING_TAG>
4077 void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION)
4078 {
4079   const char * LOC = "void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) : ";
4080   int i,j;
4081   if (_support == (SUPPORT *) NULL)
4082       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No Support defined."));
4083
4084   MESH * mesh = _support->getMesh();
4085   int spaceDim = mesh->getSpaceDimension();
4086   const double * coord;
4087
4088   const double * bary;
4089   FIELD<double,FullInterlace> * barycenterField=0;
4090
4091   double ** xyz=new double* [spaceDim];
4092   bool deallocateXyz=false;
4093   if(_support->getEntity()==MED_EN::MED_NODE)
4094     {
4095       if (_support->isOnAllElements())
4096         {
4097           coord=mesh->getCoordinates(MED_EN::MED_NO_INTERLACE);
4098           for(i=0; i<spaceDim; i++)
4099             xyz[i]=(double *)coord+i*_numberOfValues;
4100         }
4101       else
4102         {
4103           coord = mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
4104           const int * nodesNumber=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
4105           for(i=0; i<spaceDim; i++)
4106             xyz[i]=new double[_numberOfValues];
4107           deallocateXyz=true;
4108           for(i=0;i<_numberOfValues;i++)
4109             {
4110               for(j=0;j<spaceDim;j++)
4111                 xyz[j][i]=coord[(nodesNumber[i]-1)*spaceDim+j];
4112             }
4113         }
4114     }
4115   else
4116     {
4117       barycenterField = mesh->getBarycenter(_support);
4118       bary = barycenterField->getValue();
4119       //for(i=0; i<spaceDim; i++)
4120       //  xyz[i]=(double *)(bary+i*_numberOfValues);
4121       for(i=0; i<spaceDim; i++)
4122         xyz[i]=new double[_numberOfValues];
4123       deallocateXyz=true;
4124       for(i=0;i<_numberOfValues;i++) {
4125         for(j=0;j<spaceDim;j++)
4126           xyz[j][i]=bary[i*spaceDim+j];
4127       }
4128     }
4129   T* valsToSet=(T*)getValue();
4130   double *temp=new double[spaceDim];
4131   for(i=0;i<_numberOfValues;i++)
4132   {
4133     for(j=0;j<spaceDim;j++)
4134       temp[j]=xyz[j][i];
4135     f(temp,valsToSet+i*_numberOfComponents);
4136   }
4137   delete [] temp;
4138   if(barycenterField)
4139     delete barycenterField;
4140   if(deallocateXyz)
4141     for(j=0;j<spaceDim;j++)
4142       delete [] xyz[j];
4143   delete [] xyz;
4144 }
4145
4146 /*!
4147   Execute a function on _values on 'this' and put the result on a newly created field that has to be deallocated.
4148   WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
4149   Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
4150  */
4151 template <class T, class INTERLACING_TAG>
4152 FIELD<T,INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::execFunc(int nbOfComponents,
4153                                                               myFuncType2 f) throw (MEDEXCEPTION)
4154 {
4155   FIELD<T,INTERLACING_TAG> *ret=new FIELD<T,INTERLACING_TAG>(_support,nbOfComponents);
4156   const T* valsInput=getValue();
4157   T* valsOutPut=(T*)ret->getValue();
4158   int i;
4159   for(i=0;i<_numberOfValues;i++)
4160     f(valsInput+i*_numberOfComponents,valsOutPut+i*nbOfComponents);
4161   return ret;
4162 }
4163
4164 }//End namespace MEDMEM
4165
4166 #endif /* FIELD_HXX */