src/MEDMEM/MEDMEM_Field.hxx

   1 /*
   2  File Field.hxx
   3  $Header$
   4 */
   5
   6 #ifndef FIELD_HXX
   7 #define FIELD_HXX
   8
   9 #include <vector>
  10 #include <algorithm>
  11 #include <cmath>
  12
  13 #include "utilities.h"
  14 #include "MEDMEM_Exception.hxx"
  15 #include "MEDMEM_define.hxx"
  16 #include "MEDMEM_Support.hxx"
  17 #include "MEDMEM_Unit.hxx"
  18 #include "MEDMEM_Array.hxx"
  19 #include "MEDMEM_GenDriver.hxx"
  20
  21 /*!
  22
  23   This class contains all the informations related with a template class FIELD :
  24   - Components descriptions
  25   - Time step description
  26   - Location of the values (a SUPPORT class)
  27
  28 */
  29
  30 namespace MEDMEM {
  31
  32   template<class T> class FIELD;
  33
  34 class FIELD_    // GENERIC POINTER TO a template <class T> class FIELD
  35 {
  36 protected:
  37
  38   bool            _isRead ;
  39
  40   /*!
  41     \if developper
  42     Field name.
  43     \endif
  44   */
  45   string          _name ;
  46   /*!
  47     \if developper
  48     Field description.
  49     \endif
  50   */
  51   string          _description ;
  52   /*!
  53     \if developper
  54     Pointer to the support the field deals with.
  55     \endif
  56   */
  57   const SUPPORT * _support ;
  58
  59   /*!
  60     \if developper
  61     Number of field's components.
  62     \endif
  63   */
  64   int             _numberOfComponents ;
  65   /*!
  66     \if developper
  67     Number of field's values.
  68     \endif
  69   */
  70   int             _numberOfValues ;
  71
  72   /*!
  73     \if developper
  74     Array of size _numberOfComponents. \n
  75     (constant, scalar, vector, tensor)\n
  76     We could use an array of integer to store
  77     numbers of values: \n
  78     - 1 for scalar,\n
  79     - space dimension for vector,\n
  80     - space dimension square for tensor.\n
  81     So numbers of values per entities would be
  82     sum of _componentsTypes array.
  83
  84     Not implemented yet! All type are scalar !
  85     \endif
  86   */
  87   int *           _componentsTypes ;
  88   /*!
  89     \if developper
  90     Array of size _numberOfComponents
  91     storing components names if any.
  92     \endif
  93   */
  94   string *        _componentsNames;
  95   /*!
  96     \if developper
  97     Array of size _numberOfComponents
  98     storing components descriptions if any.
  99     \endif
 100  */
 101   string *        _componentsDescriptions;
 102   /*!
 103     \if developper
 104     Array of size _numberOfComponents
 105     storing components units if any.
 106     \endif
 107  */
 108   UNIT *          _componentsUnits;
 109   /*!
 110     \if developper
 111     Array of size _numberOfComponents
 112     storing components units if any.
 113     \endif
 114   */
 115   string *        _MEDComponentsUnits;
 116   int             _iterationNumber ;
 117   double          _time;
 118   int             _orderNumber ;
 119
 120   // _valueType should be a static const. Here is an initialization exemple
 121   // template < classType T > struct SET_VALUE_TYPE { static const med_type_champ _valueType = 0; }
 122   // template < > struct SET_VALUE_TYPE<double> { static const med_type_champ _valueType = MED_EN::MED_REEL64; }
 123   // template < > struct SET_VALUE_TYPE<int> { static const med_type_champ _valueType = MED_EN::MED_INT32; }
 124   // static const med_type_champ _valueType = SET_VALUE_TYPE <T>::_valueType;
 125   MED_EN::med_type_champ _valueType ;
 126
 127   vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
 128   static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
 129   static void _deepCheckFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true ) throw (MEDEXCEPTION);
 130   void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL) const  throw (MEDEXCEPTION);
 131   FIELD<double>* _getFieldSize() const;
 132
 133 public:
 134
 135   friend class MED_MED_RDONLY_DRIVER21;
 136   friend class MED_MED_WRONLY_DRIVER21;
 137   friend class MED_MED_RDWR_DRIVER21;
 138   friend class MED_MED_RDONLY_DRIVER22;
 139   friend class MED_MED_WRONLY_DRIVER22;
 140   friend class MED_MED_RDWR_DRIVER22;
 141
 142   friend class VTK_MED_DRIVER;
 143
 144   /*!
 145     Constructor.
 146   */
 147   FIELD_ ();
 148   /*!
 149     Constructor.
 150   */
 151   FIELD_(const SUPPORT * Support, const int NumberOfComponents);
 152   /*!
 153     Copy constructor.
 154   */
 155   FIELD_(const FIELD_ &m);
 156
 157   /*!
 158     Destructor.
 159   */
 160   virtual ~FIELD_();
 161
 162 //    virtual  void     setIterationNumber (int IterationNumber);
 163 //    virtual  void     setOrderNumber     (int OrderNumber);
 164 //    virtual  void     setFieldName       (string& fieldName);
 165
 166   virtual  void     rmDriver(int index=0);
 167   virtual   int     addDriver(driverTypes driverType,
 168                               const string & fileName="Default File Name.med",
 169                               const string & driverFieldName="Default Field Nam",
 170                               MED_EN::med_mode_acces access=MED_EN::MED_REMP) ;
 171
 172   virtual  int      addDriver( GENDRIVER & driver);
 173   virtual  void     read (const GENDRIVER &);
 174   virtual  void     read(int index=0);
 175   virtual void openAppend( void );
 176   virtual  void     write(const GENDRIVER &);
 177   virtual  void     write(int index=0, const string & driverName="");
 178
 179   virtual  void     writeAppend(const GENDRIVER &);
 180   virtual  void     writeAppend(int index=0, const string & driverName="");
 181
 182   inline void     setName(const string Name);
 183   inline string   getName() const;
 184   inline void     setDescription(const string Description);
 185   inline string   getDescription() const;
 186   inline const SUPPORT * getSupport() const;
 187   inline void     setSupport(const SUPPORT * support);
 188   inline void     setNumberOfComponents(const int NumberOfComponents);
 189   inline int      getNumberOfComponents() const;
 190   inline void     setNumberOfValues(const int NumberOfValues);
 191   inline int      getNumberOfValues() const;
 192   //    inline void     setComponentType(int *ComponentType);
 193   //    inline int *    getComponentType() const;
 194   //    inline int      getComponentTypeI(int i) const;
 195   inline void     setComponentsNames(const string * ComponentsNames);
 196   inline void     setComponentName(int i, const string ComponentName);
 197   inline const string * getComponentsNames() const;
 198   inline string   getComponentName(int i) const;
 199   inline void     setComponentsDescriptions(const string * ComponentsDescriptions);
 200   inline void     setComponentDescription(int i, const string ComponentDescription);
 201   inline const string * getComponentsDescriptions() const;
 202   inline string   getComponentDescription(int i) const;
 203
 204   // provisoire : en attendant de regler le probleme des unites !
 205   inline void     setComponentsUnits(const UNIT * ComponentsUnits);
 206   inline const UNIT *   getComponentsUnits() const;
 207   inline const UNIT *   getComponentUnit(int i) const;
 208   inline void     setMEDComponentsUnits(const string * MEDComponentsUnits);
 209   inline void     setMEDComponentUnit(int i, const string MEDComponentUnit);
 210   inline const string * getMEDComponentsUnits() const;
 211   inline string   getMEDComponentUnit(int i) const;
 212
 213   inline void     setIterationNumber(int IterationNumber);
 214   inline int      getIterationNumber() const;
 215   inline void     setTime(double Time);
 216   inline double   getTime() const;
 217   inline void     setOrderNumber(int OrderNumber);
 218   inline int      getOrderNumber() const;
 219
 220   inline void     setValueType (const MED_EN::med_type_champ ValueType) ;
 221   inline MED_EN::med_type_champ getValueType () const;
 222 protected:
 223   void copyGlobalInfo(const FIELD_& m);
 224 };
 225
 226 // ---------------------------------
 227 // Implemented Methods : constructor
 228 // ---------------------------------
 229
 230 // -----------------
 231 // Methodes Inline
 232 // -----------------
 233 /*!
 234   Set FIELD name.
 235 */
 236 inline void FIELD_::setName(const string Name)
 237 {
 238   _name=Name;
 239 }
 240 /*!
 241   Get FIELD name.
 242 */
 243 inline string FIELD_::getName() const
 244 {
 245   return _name;
 246 }
 247 /*!
 248   Set FIELD description.
 249 */
 250 inline void FIELD_::setDescription(const string Description)
 251 {
 252   _description=Description;
 253 }
 254 /*!
 255   Get FIELD description.
 256 */
 257 inline string FIELD_::getDescription() const
 258 {
 259   return _description;
 260 }
 261 /*!
 262   Set FIELD number of components.
 263 */
 264 inline void FIELD_::setNumberOfComponents(const int NumberOfComponents)
 265 {
 266   _numberOfComponents=NumberOfComponents;
 267 }
 268 /*!
 269   Get FIELD number of components.
 270 */
 271 inline int FIELD_::getNumberOfComponents() const
 272 {
 273   return _numberOfComponents ;
 274 }
 275 /*!
 276   Set FIELD number of values.
 277
 278   It must be the same than in the associated SUPPORT object.
 279 */
 280 inline void FIELD_::setNumberOfValues(const int NumberOfValues)
 281 {
 282   _numberOfValues=NumberOfValues;
 283 }
 284 /*!
 285   Get FIELD number of value.
 286 */
 287 inline int FIELD_::getNumberOfValues() const
 288 {
 289   return _numberOfValues ;
 290 }
 291
 292 //  inline void FIELD_::setComponentType(int *ComponentType)
 293 //  {
 294 //    _componentsTypes=ComponentType ;
 295 //  }
 296 //  inline int * FIELD_::getComponentType() const
 297 //  {
 298 //    return _componentsTypes ;
 299 //  }
 300 //  inline int FIELD_::getComponentTypeI(int i) const
 301 //  {
 302 //    return _componentsTypes[i-1] ;
 303 //  }
 304
 305 /*!
 306   Set FIELD components names.
 307
 308   Duplicate the ComponentsNames string array to put components names in
 309   FIELD. ComponentsNames size must be equal to number of components.
 310 */
 311 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
 312 {
 313   if (NULL == _componentsNames)
 314     _componentsNames = new string[_numberOfComponents] ;
 315   for (int i=0; i<_numberOfComponents; i++)
 316     _componentsNames[i]=ComponentsNames[i] ;
 317 }
 318 /*!
 319   Set FIELD i^th component name.
 320
 321   i must be >=1 and <= number of components.
 322 */
 323 inline void FIELD_::setComponentName(int i, const string ComponentName)
 324 {
 325   _componentsNames[i-1]=ComponentName ;
 326 }
 327 /*!
 328   Get a reference to the string array which contain the components names.
 329
 330   This Array size is equal to number of components
 331 */
 332 inline const string * FIELD_::getComponentsNames() const
 333 {
 334   return _componentsNames ;
 335 }
 336 /*!
 337   Get the name of the i^th component.
 338 */
 339 inline string FIELD_::getComponentName(int i) const
 340 {
 341   return _componentsNames[i-1] ;
 342 }
 343 /*!
 344   Set FIELD components descriptions.
 345
 346   Duplicate the ComponentsDescriptions string array to put components
 347   descriptions in FIELD.
 348   ComponentsDescriptions size must be equal to number of components.
 349 */
 350 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
 351 {
 352   if (NULL == _componentsDescriptions)
 353     _componentsDescriptions = new string[_numberOfComponents] ;
 354   for (int i=0; i<_numberOfComponents; i++)
 355     _componentsDescriptions[i]=ComponentsDescriptions[i] ;
 356 }
 357 /*!
 358   Set FIELD i^th component description.
 359
 360   i must be >=1 and <= number of components.
 361 */
 362 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
 363 {
 364   _componentsDescriptions[i-1]=ComponentDescription ;
 365 }
 366 /*!
 367   Get a reference to the string array which contain the components descriptions.
 368
 369   This Array size is equal to number of components
 370 */
 371 inline const string * FIELD_::getComponentsDescriptions() const
 372 {
 373   return _componentsDescriptions ;
 374 }
 375 /*!
 376   Get the description of the i^th component.
 377 */
 378 inline string FIELD_::getComponentDescription(int i) const
 379 {
 380   return _componentsDescriptions[i-1];
 381 }
 382
 383 /*!
 384   \todo
 385   Set FIELD components UNIT.
 386
 387   Duplicate the ComponentsUnits UNIT array to put components
 388   units in FIELD.
 389   ComponentsUnits size must be equal to number of components.
 390 */
 391 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
 392 {
 393   if (NULL == _componentsUnits)
 394     _componentsUnits = new UNIT[_numberOfComponents] ;
 395   for (int i=0; i<_numberOfComponents; i++)
 396     _componentsUnits[i]=ComponentsUnits[i] ;
 397 }
 398 /*!
 399   Get a reference to the UNIT array which contain the components units.
 400
 401   This Array size is equal to number of components
 402 */
 403 inline const UNIT * FIELD_::getComponentsUnits() const
 404 {
 405   return _componentsUnits ;
 406 }
 407 /*!
 408   Get the UNIT of the i^th component.
 409 */
 410 inline const UNIT * FIELD_::getComponentUnit(int i) const
 411 {
 412   return &_componentsUnits[i-1] ;
 413 }
 414 /*!
 415   Set FIELD components unit.
 416
 417   Duplicate the MEDComponentsUnits string array to put components
 418   units in FIELD.
 419   MEDComponentsUnits size must be equal to number of components.
 420
 421 */
 422 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
 423 {
 424   if (NULL == _MEDComponentsUnits)
 425     _MEDComponentsUnits = new string[_numberOfComponents] ;
 426   for (int i=0; i<_numberOfComponents; i++)
 427     _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
 428 }
 429 /*!
 430   Set FIELD i^th component unit.
 431
 432   i must be >=1 and <= number of components.
 433 */
 434 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
 435 {
 436   _MEDComponentsUnits[i-1]=MEDComponentUnit ;
 437 }
 438 /*!
 439   Get a reference to the string array which contain the components units.
 440
 441   This Array size is equal to number of components
 442 */
 443 inline const string * FIELD_::getMEDComponentsUnits() const
 444 {
 445   return _MEDComponentsUnits ;
 446 }
 447 /*!
 448   Get the string for unit of the i^th component.
 449 */
 450 inline string FIELD_::getMEDComponentUnit(int i) const
 451 {
 452   return _MEDComponentsUnits[i-1] ;
 453 }
 454 /*!
 455   Set the iteration number where FIELD has been calculated.
 456 */
 457 inline void FIELD_::setIterationNumber(int IterationNumber)
 458 {
 459   _iterationNumber=IterationNumber;
 460 }
 461 /*!
 462   Get the iteration number where FIELD has been calculated.
 463 */
 464 inline int FIELD_::getIterationNumber() const
 465 {
 466   return _iterationNumber ;
 467 }
 468 /*!
 469   Set the time (in second) where FIELD has been calculated.
 470 */
 471 inline void FIELD_::setTime(double Time)
 472 {
 473   _time=Time ;
 474 }
 475 /*!
 476   Get the time (in second) where FIELD has been calculated.
 477 */
 478 inline double FIELD_::getTime() const
 479 {
 480   return _time ;
 481 }
 482 /*!
 483   Set the order number where FIELD has been calculated.
 484
 485   It corresponds to internal iteration during one time step.
 486 */
 487 inline void FIELD_::setOrderNumber(int OrderNumber)
 488 {
 489   _orderNumber=OrderNumber ;
 490 }
 491 /*!
 492   Get the order number where FIELD has been calculated.
 493 */
 494 inline int FIELD_::getOrderNumber() const
 495 {
 496   return _orderNumber ;
 497 }
 498 /*!
 499   Get a reference to the SUPPORT object associated to FIELD.
 500 */
 501 inline  const SUPPORT * FIELD_::getSupport() const
 502 {
 503   return _support ;
 504 }
 505 /*!
 506   Set the reference to the SUPPORT object associated to FIELD.
 507
 508   Reference is not duplicate, so it must not be deleted.
 509 */
 510 inline void FIELD_::setSupport(const SUPPORT * support)
 511 {
 512   _support = support ;
 513 }
 514 /*!
 515   Get the FIELD med value type (MED_INT32 or MED_REEL64).
 516 */
 517 inline MED_EN::med_type_champ FIELD_::getValueType () const
 518 {
 519   return _valueType ;
 520 }
 521 /*!
 522   Set the FIELD med value type (MED_INT32 or MED_REEL64).
 523 */
 524 inline void FIELD_::setValueType (const MED_EN::med_type_champ ValueType)
 525 {
 526   _valueType = ValueType ;
 527 }
 528
 529 }//End namespace MEDMEM
 530
 531 /////////////////////////
 532 // END OF CLASS FIELD_ //
 533 /////////////////////////
 534
 535 /*!
 536
 537   This template class contains informations related with a FIELD :
 538   - Values of the field
 539
 540 */
 541
 542 namespace MEDMEM {
 543
 544   template<class T2> class MED_FIELD_RDONLY_DRIVER21;
 545   template<class T2> class MED_FIELD_WRONLY_DRIVER21;
 546   template<class T2> class MED_FIELD_RDONLY_DRIVER22;
 547   template<class T2> class MED_FIELD_WRONLY_DRIVER22;
 548   template<class T2> class VTK_FIELD_DRIVER;
 549
 550 template <class T> class FIELD : public FIELD_
 551 {
 552 protected:
 553   // ------ End of Drivers Management Part
 554
 555   // array of value of type T
 556   MEDARRAY<T> *_value ;
 557   static T _scalarForPow;
 558   static T pow(T x);
 559 private:
 560   void _operation(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode, char* Op);
 561   void _operationInitialize(const FIELD& m,const FIELD& n, char* Op);
 562   void _add_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode);
 563   void _sub_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode);
 564   void _mul_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode);
 565   void _div_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode) throw (MEDEXCEPTION);
 566   //setValueType() ;
 567
 568   FIELD & operator=(const FIELD &m);    // A FAIRE
 569 public:
 570   FIELD();
 571   FIELD(const FIELD &m);
 572   FIELD(const SUPPORT * Support, const int NumberOfComponents, const MED_EN::medModeSwitch Mode=MED_EN::MED_FULL_INTERLACE)  throw (MEDEXCEPTION) ; // Ajout NB Constructeur FIELD avec allocation de memoire de tous ses attribut
 573   FIELD(const SUPPORT * Support, driverTypes driverType,
 574         const string & fileName="", const string & fieldName="",
 575         const int iterationNumber = -1, const int orderNumber = -1)
 576     throw (MEDEXCEPTION);
 577   ~FIELD();
 578
 579   const FIELD operator+(const FIELD& m) const;
 580   const FIELD operator-(const FIELD& m) const;
 581   const FIELD operator*(const FIELD& m) const;
 582   const FIELD operator/(const FIELD& m) const;
 583   const FIELD operator-() const;
 584   FIELD& operator+=(const FIELD& m);
 585   FIELD& operator-=(const FIELD& m);
 586   FIELD& operator*=(const FIELD& m);
 587   FIELD& operator/=(const FIELD& m);
 588   static FIELD* add(const FIELD& m, const FIELD& n);
 589   static FIELD* addDeep(const FIELD& m, const FIELD& n);
 590   static FIELD* sub(const FIELD& m, const FIELD& n);
 591   static FIELD* subDeep(const FIELD& m, const FIELD& n);
 592   static FIELD* mul(const FIELD& m, const FIELD& n);
 593   static FIELD* mulDeep(const FIELD& m, const FIELD& n);
 594   static FIELD* div(const FIELD& m, const FIELD& n);
 595   static FIELD* divDeep(const FIELD& m, const FIELD& n);
 596   double normMax() const throw (MEDEXCEPTION);
 597   double norm2() const throw (MEDEXCEPTION);
 598   void   applyLin(T a, T b);
 599   template <T T_function(T)> void applyFunc();
 600   void applyPow(T scalar);
 601   static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
 602   double normL2(int component, const FIELD<double> * p_field_volume=NULL) const;
 603   double normL2(const FIELD<double> * p_field_volume=NULL) const;
 604   double normL1(int component, const FIELD<double> * p_field_volume=NULL) const;
 605   double normL1(const FIELD<double> * p_field_volume=NULL) const;
 606   FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
 607
 608   friend class MED_FIELD_RDONLY_DRIVER21<T>;
 609   friend class MED_FIELD_WRONLY_DRIVER21<T>;
 610   friend class MED_FIELD_RDONLY_DRIVER22<T>;
 611   friend class MED_FIELD_WRONLY_DRIVER22<T>;
 612   friend class VTK_FIELD_DRIVER<T>;
 613   //friend class MED_FIELD_RDWR_DRIVER  <T>;
 614
 615   void init ();
 616   void rmDriver(int index=0);
 617   int  addDriver(driverTypes driverType,
 618                  const string & fileName="Default File Name.med",
 619                  const string & driverFieldName="Default Field Name",
 620                  MED_EN::med_mode_acces access=MED_EN::MED_REMP) ;
 621
 622   int  addDriver(GENDRIVER & driver);
 623
 624   void allocValue(const int NumberOfComponents);
 625   void allocValue(const int NumberOfComponents, const int LengthValue);
 626
 627   void deallocValue();
 628
 629   inline void read(int index=0);
 630   inline void read(const GENDRIVER & genDriver);
 631   inline void write(int index=0, const string & driverName = "");
 632   inline void write(const GENDRIVER &);
 633
 634   inline void writeAppend(int index=0, const string & driverName = "");
 635   inline void writeAppend(const GENDRIVER &);
 636
 637   inline void     setValue(MEDARRAY<T> *Value);
 638
 639   inline MEDARRAY<T>* getvalue() const;
 640   inline int getValueLength(MED_EN::medModeSwitch Mode) const;
 641   inline const T*       getValue(MED_EN::medModeSwitch Mode) const;
 642   inline const T*       getValueI(MED_EN::medModeSwitch Mode,int i) const;
 643   inline T        getValueIJ(int i,int j) const;
 644   bool getValueOnElement(int eltIdInSup,T* retValues) const;
 645
 646   inline void setValue(MED_EN::medModeSwitch mode, T* value);
 647   inline void setValueI(MED_EN::medModeSwitch mode, int i, T* value);
 648   inline void setValueIJ(int i, int j, T value);
 649
 650   /*!
 651     This fonction feeds the FIELD<double> private attributs _value with the
 652     volume of each cells belonging to the argument Support. The field has to be
 653     initialised via the constructor FIELD<double>(const SUPPORT * , const int )
 654     with Support as SUPPORT argument, 1 has the number of components, and Support
 655     has be a SUPPORT on 3D cells. This initialisation could be done by the empty
 656     constructor followed by a setSupport and setNumberOfComponents call but it has
 657     be followed by a setValueType(MED_REEL64) call.
 658    */
 659   void getVolume() const throw (MEDEXCEPTION) ;
 660   /*!
 661     This fonction feeds the FIELD<double> private attributs _value with the
 662     area of each cells (or faces) belonging to the attribut _support. The field
 663     has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
 664     const int ) with 1 has the number of components, and _support has be a
 665     SUPPORT on 2D cells or 3D faces. This initialisation could be done by the
 666     empty constructor followed by a setSupport and setNumberOfComponents call but
 667     it has be followed by a setValueType(MED_REEL64) call.
 668    */
 669   void getArea() const throw (MEDEXCEPTION) ;
 670   /*!
 671     This fonction feeds the FIELD<double> private attributs _value with the
 672     length of each segments belonging to the attribut _support. The field has
 673     to be initialised via the constructor FIELD<double>(const SUPPORT * ,
 674     const int ) with 1 has the number of components, and _support has be a
 675     SUPPORT on 3D edges or 2D faces. This initialisation could be done by the
 676     empty constructor followed by a setSupport and setNumberOfComponents call but
 677     it has be followed by a setValueType(MED_REEL64) call.
 678    */
 679   void getLength() const throw (MEDEXCEPTION) ;
 680   /*!
 681     This fonction feeds the FIELD<double> private attributs _value with the
 682     normal vector of each faces belonging to the attribut _support. The field
 683     has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
 684     const int ) with the space dimension has the number of components, and
 685     _support has be a SUPPORT on 3D or 2D faces. This initialisation could be done
 686     by the empty constructor followed by a setSupport and setNumberOfComponents
 687     call but it has be followed by a setValueType(MED_REEL64) call.
 688    */
 689   void getNormal() const throw (MEDEXCEPTION) ;
 690   /*!
 691     This fonction feeds the FIELD<double> private attributs _value with the
 692     barycenter of each faces or cells or edges belonging to the attribut _support.
 693     The field has to be initialised via the constructor
 694     FIELD<double>(const SUPPORT * ,const int ) with the space dimension has the
 695     number of components, and _support has be a SUPPORT on 3D cells or 2D faces.
 696     This initialisation could be done by the empty constructor followed by a
 697     setSupport and setNumberOfComponents call but it has be followed by a
 698     setValueType(MED_REEL64) call.
 699    */
 700   void getBarycenter() const throw (MEDEXCEPTION) ;
 701   template<void T_Analytic(const double *,T*)>
 702   void fillFromAnalytic();
 703 };
 704 }
 705
 706 #include "MEDMEM_DriverFactory.hxx"
 707
 708 namespace MEDMEM {
 709
 710 template <class T> T FIELD<T>::_scalarForPow=1;
 711
 712 // --------------------
 713 // Implemented Methods
 714 // --------------------
 715
 716 /*!
 717   Constructor with no parameter, most of the attribut members are set to NULL.
 718 */
 719 template <class T>  FIELD<T>::FIELD():
 720  _value((MEDARRAY<T>*)NULL)
 721 {
 722   MESSAGE("Constructeur FIELD sans parametre");
 723 }
 724
 725 /*!
 726   Constructor with parameters such that all attrribut are set but the _value
 727   attribut is allocated but not set.
 728 */
 729 template <class T>  FIELD<T>::FIELD(const SUPPORT * Support,
 730                                     const int NumberOfComponents, const MED_EN::medModeSwitch Mode) throw (MEDEXCEPTION) :
 731   FIELD_(Support, NumberOfComponents)
 732 {
 733   BEGIN_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
 734   SCRUTE(this);
 735
 736   try {
 737     _numberOfValues = Support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 738   }
 739   catch (MEDEXCEPTION &ex) {
 740     MESSAGE("No value defined ! ("<<ex.what()<<")");
 741     _value = (MEDARRAY<T>*)NULL ;
 742   }
 743   MESSAGE("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
 744   if (0<_numberOfValues) {
 745     _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues,Mode);
 746     _isRead = true ;
 747   } else
 748     _value = (MEDARRAY<T>*)NULL ;
 749
 750   END_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
 751 }
 752
 753 /*!
 754   \if developper
 755   \endif
 756 */
 757 template <class T> void FIELD<T>::init ()
 758 {
 759 }
 760
 761 /*!
 762   Copy constructor.
 763 */
 764 template <class T> FIELD<T>::FIELD(const FIELD & m):
 765   FIELD_((FIELD_) m)
 766 {
 767   MESSAGE("Constructeur FIELD de recopie");
 768   if (m._value != NULL)
 769     {
 770       // copie only default !
 771       _value = new MEDARRAY<T>(* m._value,false);
 772     }
 773   else
 774     _value = (MEDARRAY<T> *) NULL;
 775   //_drivers = m._drivers;
 776 }
 777
 778 /*!
 779   \if developper
 780   Not implemented.
 781   \endif
 782 */
 783 template <class T> FIELD<T> & FIELD<T>::operator=(const FIELD &m)
 784 {
 785   MESSAGE("Appel de FIELD<T>::operator=") ;
 786   // operator= on FIELD_
 787   // ignore driver
 788   // copy values array
 789 }
 790
 791 /*!
 792      Overload addition operator.
 793      This operation is authorized only for compatible fields that have the same support.
 794      The compatibility checking includes equality tests of the folowing data members:\n
 795      - _support
 796      - _numberOfComponents
 797      - _numberOfValues
 798      - _componentsTypes
 799      - _MEDComponentsUnits.
 800
 801      The data members of the returned field are initialized, based on the first field, except for the name,
 802      which is the combination of the two field's names and the operator.
 803      Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> \n
 804      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
 805      When using python, this operator calls the copy constructor in any case.
 806      The user has to be aware that when using operator + in associatives expressions like
 807      <tt> a = b + c + d +e; </tt> \n
 808      no optimisation is performed : the evaluation of last expression requires the construction of
 809      3 temporary fields.
 810 */
 811 template <class T>
 812 const FIELD<T> FIELD<T>::operator+(const FIELD & m) const
 813 {
 814     BEGIN_OF("FIELD<T>::operator+(const FIELD & m)");
 815     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
 816
 817     // Select mode : avoid if possible any calculation of other mode for fields m or *this
 818     MED_EN::medModeSwitch mode;
 819     if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
 820         mode=m.getvalue()->getMode();
 821     else
 822         mode=this->getvalue()->getMode();
 823
 824     // Creation of the result - memory is allocated by FIELD constructor
 825     FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
 826     //result._operation(*this,m,mode,"+"); // perform Atribute's initialization & addition
 827     result._operationInitialize(*this,m,"+"); // perform Atribute's initialization
 828     result._add_in_place(*this,m,mode); // perform addition
 829
 830     END_OF("FIELD<T>::operator+(const FIELD & m)");
 831     return result;
 832 }
 833
 834 /*!  Overloaded Operator +=
 835  *   Operations are directly performed in the first field's array.
 836  *   This operation is authorized only for compatible fields that have the same support.
 837  */
 838 template <class T>
 839 FIELD<T>& FIELD<T>::operator+=(const FIELD & m)
 840 {
 841     BEGIN_OF("FIELD<T>::operator+=(const FIELD & m)");
 842     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
 843
 844     // We choose to keep *this mode, even if it may cost a re-calculation for m
 845     MED_EN::medModeSwitch mode=this->getvalue()->getMode();
 846     const T* value1=m.getValue(mode); // get pointers to the values we are adding
 847
 848     // get a non const pointer to the inside array of values and perform operation
 849     T * value=const_cast<T *> (getValue(mode));
 850     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
 851     const T* endV=value+size; // pointer to the end of value
 852     for(;value!=endV; value1++,value++)
 853         *value += *value1;
 854     // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
 855     this->getvalue()->clearOtherMode();
 856     END_OF("FIELD<T>::operator+=(const FIELD & m)");
 857     return *this;
 858 }
 859
 860
 861 /*! Addition of fields. Static member function.
 862  *  The function return a pointer to a new created field that holds the addition.
 863  *  Data members are checked for compatibility and initialized.
 864  *  The user is in charge of memory deallocation.
 865  */
 866 template <class T>
 867 FIELD<T>* FIELD<T>::add(const FIELD& m, const FIELD& n)
 868 {
 869     BEGIN_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
 870     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
 871
 872     // Select mode : avoid if possible any calculation of other mode for fields m or *this
 873     MED_EN::medModeSwitch mode;
 874     if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
 875         mode=m.getvalue()->getMode();
 876     else
 877         mode=n.getvalue()->getMode();
 878
 879     // Creation of a new field
 880     FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
 881     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
 882     result->_add_in_place(m,n,mode); // perform addition
 883
 884     END_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
 885     return result;
 886 }
 887
 888 /*! Same as add method except that field check is deeper.
 889  */
 890 template <class T>
 891 FIELD<T>* FIELD<T>::addDeep(const FIELD& m, const FIELD& n)
 892 {
 893     BEGIN_OF("FIELD<T>::addDeep(const FIELD & m, const FIELD& n)");
 894     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
 895
 896     // Select mode : avoid if possible any calculation of other mode for fields m or *this
 897     MED_EN::medModeSwitch mode;
 898     if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
 899         mode=m.getvalue()->getMode();
 900     else
 901         mode=n.getvalue()->getMode();
 902
 903     // Creation of a new field
 904     FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
 905     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
 906     result->_add_in_place(m,n,mode); // perform addition
 907
 908     END_OF("FIELD<T>::addDeep(const FIELD & m, const FIELD& n)");
 909     return result;
 910 }
 911
 912 /*!
 913      Overload substraction operator.
 914      This operation is authorized only for compatible fields that have the same support.
 915      The compatibility checking includes equality tests of the folowing data members:\n
 916      - _support
 917      - _numberOfComponents
 918      - _numberOfValues
 919      - _componentsTypes
 920      - _MEDComponentsUnits.
 921
 922      The data members of the returned field are initialized, based on the first field, except for the name,
 923      which is the combination of the two field's names and the operator.
 924      Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> \n
 925      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
 926      When using python, this operator calls the copy constructor in any case.
 927      The user has to be aware that when using operator - in associatives expressions like
 928      <tt> a = b - c - d -e; </tt> \n
 929      no optimisation is performed : the evaluation of last expression requires the construction of
 930      3 temporary fields.
 931 */
 932 template <class T>
 933 const FIELD<T> FIELD<T>::operator-(const FIELD & m) const
 934 {
 935     BEGIN_OF("FIELD<T>::operator-(const FIELD & m)");
 936     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
 937
 938     MED_EN::medModeSwitch mode;
 939     if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
 940         mode=m.getvalue()->getMode();
 941     else
 942         mode=this->getvalue()->getMode();
 943
 944     // Creation of the result - memory is allocated by FIELD constructor
 945     FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
 946     //result._operation(*this,m,mode,"-"); // perform Atribute's initialization & substraction
 947     result._operationInitialize(*this,m,"-"); // perform Atribute's initialization
 948     result._sub_in_place(*this,m,mode); // perform substracion
 949
 950     END_OF("FIELD<T>::operator-(const FIELD & m)");
 951     return result;
 952 }
 953
 954 template <class T>
 955 const FIELD<T> FIELD<T>::operator-() const
 956 {
 957     BEGIN_OF("FIELD<T>::operator-()");
 958
 959     MED_EN::medModeSwitch mode=this->getvalue()->getMode();
 960     // Creation of the result - memory is allocated by FIELD constructor
 961     FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
 962     // Atribute's initialization
 963     result.setName("- "+getName());
 964     result.setComponentsNames(getComponentsNames());
 965     // not yet implemented    setComponentType(getComponentType());
 966     result.setComponentsDescriptions(getComponentsDescriptions());
 967     result.setMEDComponentsUnits(getMEDComponentsUnits());
 968     result.setComponentsUnits(getComponentsUnits());
 969     result.setIterationNumber(getIterationNumber());
 970     result.setTime(getTime());
 971     result.setOrderNumber(getOrderNumber());
 972     result.setValueType(getValueType());
 973
 974     const T* value1=getValue(mode);
 975     // get a non const pointer to the inside array of values and perform operation
 976     T * value=const_cast<T *> (result.getValue(mode));
 977     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
 978     const T* endV=value+size; // pointer to the end of value
 979
 980     for(;value!=endV; value1++,value++)
 981         *value = -(*value1);
 982     END_OF("FIELD<T>::operator-=(const FIELD & m)");
 983     return result;
 984 }
 985
 986 /*!  Overloaded Operator -=
 987  *   Operations are directly performed in the first field's array.
 988  *   This operation is authorized only for compatible fields that have the same support.
 989  */
 990 template <class T>
 991 FIELD<T>& FIELD<T>::operator-=(const FIELD & m)
 992 {
 993     BEGIN_OF("FIELD<T>::operator-=(const FIELD & m)");
 994     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
 995
 996     // We choose to keep *this mode, even if it may cost a re-calculation for m
 997     MED_EN::medModeSwitch mode=this->getvalue()->getMode();
 998     const T* value1=m.getValue(mode); // get pointers to the values we are adding
 999
1000     // get a non const pointer to the inside array of values and perform operation
1001     T * value=const_cast<T *> (getValue(mode));
1002     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1003     const T* endV=value+size; // pointer to the end of value
1004
1005     for(;value!=endV; value1++,value++)
1006         *value -= *value1;
1007     // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1008     this->getvalue()->clearOtherMode();
1009     END_OF("FIELD<T>::operator-=(const FIELD & m)");
1010     return *this;
1011 }
1012
1013
1014 /*! Substraction of fields. Static member function.
1015  *  The function return a pointer to a new created field that holds the substraction.
1016  *  Data members are checked for compatibility and initialized.
1017  *  The user is in charge of memory deallocation.
1018  */
1019 template <class T>
1020 FIELD<T>* FIELD<T>::sub(const FIELD& m, const FIELD& n)
1021 {
1022     BEGIN_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
1023     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1024
1025     // Select mode : avoid if possible any calculation of other mode for fields m or *this
1026     MED_EN::medModeSwitch mode;
1027     if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1028         mode=m.getvalue()->getMode();
1029     else
1030         mode=n.getvalue()->getMode();
1031
1032     // Creation of a new field
1033     FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1034     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1035     result->_sub_in_place(m,n,mode); // perform substraction
1036
1037     END_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
1038     return result;
1039 }
1040
1041 /*! Same as sub method except that field check is deeper.
1042  */
1043 template <class T>
1044 FIELD<T>* FIELD<T>::subDeep(const FIELD& m, const FIELD& n)
1045 {
1046     BEGIN_OF("FIELD<T>::subDeep(const FIELD & m, const FIELD& n)");
1047     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1048
1049     // Select mode : avoid if possible any calculation of other mode for fields m or *this
1050     MED_EN::medModeSwitch mode;
1051     if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1052         mode=m.getvalue()->getMode();
1053     else
1054         mode=n.getvalue()->getMode();
1055
1056     // Creation of a new field
1057     FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1058     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1059     result->_sub_in_place(m,n,mode); // perform substraction
1060
1061     END_OF("FIELD<T>::subDeep(const FIELD & m, const FIELD& n)");
1062     return result;
1063 }
1064
1065 /*!
1066      Overload multiplication operator.
1067      This operation is authorized only for compatible fields that have the same support.
1068      The compatibility checking includes equality tests of the folowing data members:\n
1069      - _support
1070      - _numberOfComponents
1071      - _numberOfValues
1072      - _componentsTypes
1073      - _MEDComponentsUnits.
1074
1075      The data members of the returned field are initialized, based on the first field, except for the name,
1076      which is the combination of the two field's names and the operator.
1077      Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> \n
1078      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1079      When using python, this operator calls the copy constructor in any case.
1080      The user has to be aware that when using operator * in associatives expressions like
1081      <tt> a = b * c * d *e; </tt> \n
1082      no optimisation is performed : the evaluation of last expression requires the construction of
1083      3 temporary fields.
1084 */
1085 template <class T>
1086 const FIELD<T> FIELD<T>::operator*(const FIELD & m) const
1087 {
1088     BEGIN_OF("FIELD<T>::operator*(const FIELD & m)");
1089     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1090
1091     // Select mode : avoid if possible any calculation of other mode for fields m or *this
1092     MED_EN::medModeSwitch mode;
1093     if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
1094         mode=m.getvalue()->getMode();
1095     else
1096         mode=this->getvalue()->getMode();
1097
1098     // Creation of the result - memory is allocated by FIELD constructor
1099     FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
1100     //result._operation(*this,m,mode,"*"); // perform Atribute's initialization & multiplication
1101     result._operationInitialize(*this,m,"*"); // perform Atribute's initialization
1102     result._mul_in_place(*this,m,mode); // perform multiplication
1103
1104     END_OF("FIELD<T>::operator*(const FIELD & m)");
1105     return result;
1106 }
1107
1108 /*!  Overloaded Operator *=
1109  *   Operations are directly performed in the first field's array.
1110  *   This operation is authorized only for compatible fields that have the same support.
1111  */
1112 template <class T>
1113 FIELD<T>& FIELD<T>::operator*=(const FIELD & m)
1114 {
1115     BEGIN_OF("FIELD<T>::operator*=(const FIELD & m)");
1116     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1117
1118     // We choose to keep *this mode, even if it may cost a re-calculation for m
1119     MED_EN::medModeSwitch mode=this->getvalue()->getMode();
1120     const T* value1=m.getValue(mode); // get pointers to the values we are adding
1121
1122     // get a non const pointer to the inside array of values and perform operation
1123     T * value=const_cast<T *> (getValue(mode));
1124     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1125     const T* endV=value+size; // pointer to the end of value
1126
1127     for(;value!=endV; value1++,value++)
1128         *value *= *value1;
1129     // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1130     this->getvalue()->clearOtherMode();
1131     END_OF("FIELD<T>::operator*=(const FIELD & m)");
1132     return *this;
1133 }
1134
1135
1136 /*! Multiplication of fields. Static member function.
1137  *  The function return a pointer to a new created field that holds the multiplication.
1138  *  Data members are checked for compatibility and initialized.
1139  *  The user is in charge of memory deallocation.
1140  */
1141 template <class T>
1142 FIELD<T>* FIELD<T>::mul(const FIELD& m, const FIELD& n)
1143 {
1144     BEGIN_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1145     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1146
1147     // Select mode : avoid if possible any calculation of other mode for fields m or *this
1148     MED_EN::medModeSwitch mode;
1149     if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1150         mode=m.getvalue()->getMode();
1151     else
1152         mode=n.getvalue()->getMode();
1153
1154     // Creation of a new field
1155     FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1156     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1157     result->_mul_in_place(m,n,mode); // perform multiplication
1158
1159     END_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1160     return result;
1161 }
1162
1163 /*! Same as mul method except that field check is deeper.
1164  */
1165 template <class T>
1166 FIELD<T>* FIELD<T>::mulDeep(const FIELD& m, const FIELD& n)
1167 {
1168     BEGIN_OF("FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)");
1169     FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1170
1171     // Select mode : avoid if possible any calculation of other mode for fields m or *this
1172     MED_EN::medModeSwitch mode;
1173     if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1174         mode=m.getvalue()->getMode();
1175     else
1176         mode=n.getvalue()->getMode();
1177
1178     // Creation of a new field
1179     FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1180     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1181     result->_mul_in_place(m,n,mode); // perform multiplication
1182
1183     END_OF("FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)");
1184     return result;
1185 }
1186
1187 /*!
1188      Overload division operator.
1189      This operation is authorized only for compatible fields that have the same support.
1190      The compatibility checking includes equality tests of the folowing data members:\n
1191      - _support
1192      - _numberOfComponents
1193      - _numberOfValues
1194      - _componentsTypes
1195      - _MEDComponentsUnits.
1196
1197      The data members of the returned field are initialized, based on the first field, except for the name,
1198      which is the combination of the two field's names and the operator.
1199      Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> \n
1200      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1201      When using python, this operator calls the copy constructor in any case.
1202      The user has to be aware that when using operator / in associatives expressions like
1203      <tt> a = b / c / d /e; </tt> \n
1204      no optimisation is performed : the evaluation of last expression requires the construction of
1205      3 temporary fields.
1206 */
1207 template <class T>
1208 const FIELD<T> FIELD<T>::operator/(const FIELD & m) const
1209 {
1210     BEGIN_OF("FIELD<T>::operator/(const FIELD & m)");
1211     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1212
1213     // Select mode : avoid if possible any calculation of other mode for fields m or *this
1214     MED_EN::medModeSwitch mode;
1215     if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
1216         mode=m.getvalue()->getMode();
1217     else
1218         mode=this->getvalue()->getMode();
1219
1220     // Creation of the result - memory is allocated by FIELD constructor
1221     FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
1222     //result._operation(*this,m,mode,"/"); // perform Atribute's initialization & division
1223     result._operationInitialize(*this,m,"/"); // perform Atribute's initialization
1224     result._div_in_place(*this,m,mode); // perform division
1225
1226     END_OF("FIELD<T>::operator/(const FIELD & m)");
1227     return result;
1228 }
1229
1230
1231 /*!  Overloaded Operator /=
1232  *   Operations are directly performed in the first field's array.
1233  *   This operation is authorized only for compatible fields that have the same support.
1234  */
1235 template <class T>
1236 FIELD<T>& FIELD<T>::operator/=(const FIELD & m)
1237 {
1238     BEGIN_OF("FIELD<T>::operator/=(const FIELD & m)");
1239     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1240
1241     // We choose to keep *this mode, even if it may cost a re-calculation for m
1242     MED_EN::medModeSwitch mode=this->getvalue()->getMode();
1243     const T* value1=m.getValue(mode); // get pointers to the values we are adding
1244
1245     // get a non const pointer to the inside array of values and perform operation
1246     T * value=const_cast<T *> (getValue(mode));
1247     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1248     const T* endV=value+size; // pointer to the end of value
1249
1250     for(;value!=endV; value1++,value++)
1251         *value /= *value1;
1252     // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1253     this->getvalue()->clearOtherMode();
1254     END_OF("FIELD<T>::operator/=(const FIELD & m)");
1255     return *this;
1256 }
1257
1258
1259 /*! Division of fields. Static member function.
1260  *  The function return a pointer to a new created field that holds the division.
1261  *  Data members are checked for compatibility and initialized.
1262  *  The user is in charge of memory deallocation.
1263  */
1264 template <class T>
1265 FIELD<T>* FIELD<T>::div(const FIELD& m, const FIELD& n)
1266 {
1267     BEGIN_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1268     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1269
1270     // Select mode : avoid if possible any calculation of other mode for fields m or *this
1271     MED_EN::medModeSwitch mode;
1272     if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1273         mode=m.getvalue()->getMode();
1274     else
1275         mode=n.getvalue()->getMode();
1276
1277     // Creation of a new field
1278     FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1279     result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1280     result->_div_in_place(m,n,mode); // perform division
1281
1282     END_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1283     return result;
1284 }
1285
1286 /*! Same as div method except that field check is deeper.
1287  */
1288 template <class T>
1289 FIELD<T>* FIELD<T>::divDeep(const FIELD& m, const FIELD& n)
1290 {
1291   BEGIN_OF("FIELD<T>::divDeep(const FIELD & m, const FIELD& n)");
1292   FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1293
1294   // Select mode : avoid if possible any calculation of other mode for fields m or *this
1295   MED_EN::medModeSwitch mode;
1296   if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1297     mode=m.getvalue()->getMode();
1298   else
1299     mode=n.getvalue()->getMode();
1300
1301   // Creation of a new field
1302   FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1303   result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1304   result->_div_in_place(m,n,mode); // perform division
1305
1306   END_OF("FIELD<T>::divDeep(const FIELD & m, const FIELD& n)");
1307   return result;
1308 }
1309
1310 /*!
1311   \if developper
1312   This internal method initialize the members of a new field created to hold the result of the operation Op .
1313   Initialization is based on the first field, except for the name, which is the combination of the two field's names
1314   and the operator.
1315   \endif
1316 */
1317 template <class T>
1318 void FIELD<T>::_operationInitialize(const FIELD& m,const FIELD& n, char* Op)
1319 {
1320     MESSAGE("Appel methode interne _add" << Op);
1321
1322     // Atribute's initialization (copy of the first field's attributes)
1323     // Other data members (_support, _numberOfValues) are initialized in the field's constr.
1324     setName(m.getName()+" "+Op+" "+n.getName());
1325     setComponentsNames(m.getComponentsNames());
1326     // not yet implemented    setComponentType(m.getComponentType());
1327     setComponentsDescriptions(m.getComponentsDescriptions());
1328     setMEDComponentsUnits(m.getMEDComponentsUnits());
1329
1330     // The following data member may differ from field m to n.
1331     // The initialization is done based on the first field.
1332     setComponentsUnits(m.getComponentsUnits());
1333     setIterationNumber(m.getIterationNumber());
1334     setTime(m.getTime());
1335     setOrderNumber(m.getOrderNumber());
1336     setValueType(m.getValueType());
1337 }
1338
1339
1340 /*!
1341   \if developper
1342   Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
1343   This method is applied to a just created field with medModeSwitch mode.
1344   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1345   it doesn't exist!
1346   \endif
1347 */
1348 template <class T>
1349 void FIELD<T>::_add_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode)
1350 {
1351     // get pointers to the values we are adding
1352     const T* value1=m.getValue(mode);
1353     const T* value2=n.getValue(mode);
1354     // get a non const pointer to the inside array of values and perform operation
1355     T * value=const_cast<T *> (getValue(mode));
1356
1357     const int size=getNumberOfValues()*getNumberOfComponents();
1358     SCRUTE(size);
1359     const T* endV1=value1+size;
1360     for(;value1!=endV1; value1++,value2++,value++)
1361         *value=(*value1)+(*value2);
1362 }
1363
1364 /*!
1365   \if developper
1366   Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
1367   This method is applied to a just created field with medModeSwitch mode.
1368   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1369   it doesn't exist!
1370   \endif
1371 */
1372 template <class T>
1373 void FIELD<T>::_sub_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode)
1374 {
1375     // get pointers to the values we are adding
1376     const T* value1=m.getValue(mode);
1377     const T* value2=n.getValue(mode);
1378     // get a non const pointer to the inside array of values and perform operation
1379     T * value=const_cast<T *> (getValue(mode));
1380
1381     const int size=getNumberOfValues()*getNumberOfComponents();
1382     SCRUTE(size);
1383     const T* endV1=value1+size;
1384     for(;value1!=endV1; value1++,value2++,value++)
1385         *value=(*value1)-(*value2);
1386 }
1387
1388 /*!
1389   \if developper
1390   Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
1391   This method is applied to a just created field with medModeSwitch mode.
1392   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1393   it doesn't exist!
1394   \endif
1395 */
1396 template <class T>
1397 void FIELD<T>::_mul_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode)
1398 {
1399     // get pointers to the values we are adding
1400     const T* value1=m.getValue(mode);
1401     const T* value2=n.getValue(mode);
1402     // get a non const pointer to the inside array of values and perform operation
1403     T * value=const_cast<T *> (getValue(mode));
1404
1405     const int size=getNumberOfValues()*getNumberOfComponents();
1406     SCRUTE(size);
1407     const T* endV1=value1+size;
1408     for(;value1!=endV1; value1++,value2++,value++)
1409         *value=(*value1)*(*value2);
1410 }
1411
1412 /*!
1413   \if developper
1414   Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
1415   This method is applied to a just created field with medModeSwitch mode.
1416   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1417   it doesn't exist!
1418   \endif
1419 */
1420 template <class T>
1421 void FIELD<T>::_div_in_place(const FIELD& m,const FIELD& n, const MED_EN::medModeSwitch mode) throw (MEDEXCEPTION)
1422 {
1423     // get pointers to the values we are adding
1424     const T* value1=m.getValue(mode);
1425     const T* value2=n.getValue(mode);
1426     // get a non const pointer to the inside array of values and perform operation
1427     T * value=const_cast<T *> (getValue(mode));
1428
1429     const int size=getNumberOfValues()*getNumberOfComponents();
1430     SCRUTE(size);
1431     const T* endV1=value1+size;
1432     for(;value1!=endV1; value1++,value2++,value++){
1433       if ( *value2 == 0 ) {
1434           string diagnosis;
1435           diagnosis="FIELD<T>::_div_in_place(...) : Divide by zero !";
1436           throw MEDEXCEPTION(diagnosis.c_str());
1437         }
1438         *value=(*value1)/(*value2);
1439     }
1440 }
1441
1442 /*!  Return Max Norm
1443  */
1444 template <class T> double FIELD<T>::normMax() const throw (MEDEXCEPTION)
1445 {
1446     const T* value=getValue(getvalue()->getMode()); // get pointer to the values
1447     const int size=getNumberOfValues()*getNumberOfComponents();
1448     if (size <= 0) // Size of array has to be strictly positive
1449     {
1450         string diagnosis;
1451         diagnosis="FIELD<T>::normMax() : cannot compute the norm of "+getName()+
1452             " : it size is non positive!";
1453         throw MEDEXCEPTION(diagnosis.c_str());
1454     }
1455     const T* lastvalue=value+size; // get pointer just after last value
1456     const T* pMax=value; // pointer to the max value
1457     const T* pMin=value; // pointer to the min value
1458
1459     // get pointers to the max & min value of array
1460     while ( ++value != lastvalue )
1461     {
1462         if ( *pMin > *value )
1463             pMin=value;
1464         if ( *pMax < *value )
1465             pMax=value;
1466     }
1467
1468     T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
1469     T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
1470
1471     return Max>Min ? static_cast<double>(Max) : static_cast<double>(Min);
1472 }
1473
1474 /*!  Return Euclidien norm
1475  */
1476 template <class T> double FIELD<T>::norm2() const throw (MEDEXCEPTION)
1477 {
1478     const T* value=this->getValue(this->getvalue()->getMode()); // get const pointer to the values
1479     const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
1480     if (size <= 0) // Size of array has to be strictly positive
1481     {
1482         string diagnosis;
1483         diagnosis="FIELD<T>::norm2() : cannot compute the norm of "+getName()+
1484             " : it size is non positive!";
1485         throw MEDEXCEPTION(diagnosis.c_str());
1486     }
1487     const T* lastvalue=value+size; // point just after last value
1488
1489     T result((T)0); // init
1490     for( ; value!=lastvalue ; ++value)
1491         result += (*value) * (*value);
1492
1493     return std::sqrt(static_cast<double> (result));
1494 }
1495
1496
1497 /*!  Apply to each (scalar) field component the template parameter T_function,
1498  *   which is a pointer to function.
1499  *   Since the pointer is known at compile time, the function is inlined into the inner loop!
1500  *   calculation is done "in place".
1501  *   Use examples :
1502  *
1503  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
1504  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
1505  */
1506 template <class T> template <T T_function(T)>
1507 void FIELD<T>::applyFunc()
1508 {
1509     MED_EN::medModeSwitch mode=getvalue()->getMode();
1510
1511     // get a non const pointer to the inside array of values and perform operation
1512     T * value=const_cast<T *> (getValue(mode));
1513     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1514
1515     if (size>0) // for a negative size, there is nothing to do
1516     {
1517         const T* lastvalue=value+size; // pointer to the end of value
1518         for(;value!=lastvalue; ++value) // apply linear transformation
1519             *value = T_function(*value);
1520         // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1521         getvalue()->clearOtherMode();
1522     }
1523 }
1524
1525 template <class T> T FIELD<T>::pow(T x)
1526 {
1527   return (T)::pow(x,FIELD<T>::_scalarForPow);
1528 }
1529
1530 /*!  Apply to each (scalar) field component the math function pow.
1531  *   calculation is done "in place".
1532  *   Use examples :
1533  *
1534  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
1535  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
1536  */
1537 template <class T> void FIELD<T>::applyPow(T scalar)
1538 {
1539   FIELD<T>::_scalarForPow=scalar;
1540   applyFunc<FIELD<T>::pow>();
1541 }
1542
1543 /*!  Apply to each (scalar) field component the linear function x -> ax+b.
1544  *   calculation is done "in place".
1545  */
1546 template <class T> void FIELD<T>::applyLin(T a, T b)
1547 {
1548     MED_EN::medModeSwitch mode=getvalue()->getMode();
1549
1550     // get a non const pointer to the inside array of values and perform operation in place
1551     T * value=const_cast<T *> (getValue(mode));
1552     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1553
1554     if (size>0) // for a negative size, there is nothing to do
1555     {
1556         const T* lastvalue=value+size; // pointer to the end of value
1557         for(;value!=lastvalue; ++value) // apply linear transformation
1558             *value = a*(*value)+b;
1559         // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1560         getvalue()->clearOtherMode();
1561     }
1562 }
1563
1564
1565 /*!
1566  *   Return a pointer to a new field that holds the scalar product. Static member function.
1567  *   This operation is authorized only for compatible fields that have the same support.
1568  *   The compatibility checking includes equality tests of the folowing data members:\n
1569  *   - _support
1570  *   - _numberOfComponents
1571  *   - _numberOfValues
1572  *   - _componentsTypes
1573  *   - _MEDComponentsUnits.
1574  *   Data members are initialized.
1575  *   The new field point to the same support and has one component.
1576  *   Each value of it is the scalar product of the two argument's fields.
1577  *   The user is in charge of memory deallocation.
1578  */
1579 template <class T> FIELD<T>* FIELD<T>::scalarProduct(const FIELD & m, const FIELD & n, bool deepCheck)
1580 {
1581   if(!deepCheck)
1582     FIELD_::_checkFieldCompatibility( m, n, false); // may throw exception
1583   else
1584     FIELD_::_deepCheckFieldCompatibility(m, n, false);
1585   // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
1586     const MED_EN::medModeSwitch mode=MED_EN::MED_FULL_INTERLACE;
1587
1588     const int numberOfElements=m.getNumberOfValues(); // strictly positive
1589     const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
1590
1591     // Creation & init of a the result field on the same suppot, with one component
1592     FIELD<T>* result = new FIELD<T>(m.getSupport(),1,mode);
1593     result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
1594     result->setIterationNumber(m.getIterationNumber());
1595     result->setTime(m.getTime());
1596     result->setOrderNumber(m.getOrderNumber());
1597     result->setValueType(m.getValueType());
1598
1599     const T* value1=m.getValue(mode); // get const pointer to the values
1600     const T* value2=n.getValue(mode); // get const pointer to the values
1601     // get a non const pointer to the inside array of values and perform operation
1602     T * value=const_cast<T *> (result->getValue(mode));
1603
1604     const T* lastvalue=value+numberOfElements; // pointing just after last value of result
1605     for ( ; value!=lastvalue ; ++value ) // loop on all elements
1606     {
1607         *value=(T)0; // initialize value
1608         const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
1609         for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
1610             *value += (*value1) * (*value2);
1611     }
1612     return result;
1613 }
1614
1615 /*!  Return L2 Norm  of the field's component.
1616  *   Cannot be applied to a field with a support on nodes.
1617  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1618  */
1619 template <class T> double FIELD<T>::normL2(int component, const FIELD<double> * p_field_volume) const
1620 {
1621     _checkNormCompatibility(p_field_volume); // may throw exception
1622     if ( component<1 || component>getNumberOfComponents() )
1623         throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
1624
1625     const FIELD<double> * p_field_size=p_field_volume;
1626     if(!p_field_volume) // if the user don't supply the volume
1627         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1628
1629     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1630     const double* vol=p_field_size->getValue(MED_EN::MED_FULL_INTERLACE);
1631     const T* value=getValueI( MED_EN::MED_NO_INTERLACE, component); // get pointer to the component's values
1632     const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
1633
1634     double integrale=0.0;
1635     double totVol=0.0;
1636     for (; value!=lastvalue ; ++value ,++vol)
1637     {
1638         integrale += static_cast<double>((*value) * (*value)) * (*vol);
1639         totVol+=*vol;
1640     }
1641
1642     if(!p_field_volume) // if the user didn't supply the volume
1643         delete p_field_size; // delete temporary volume field
1644     if( totVol <= 0)
1645         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1646
1647     return integrale/totVol;
1648 }
1649
1650 /*!  Return L2 Norm  of the field.
1651  *   Cannot be applied to a field with a support on nodes.
1652  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1653  */
1654 template <class T> double FIELD<T>::normL2(const FIELD<double> * p_field_volume) const
1655 {
1656     _checkNormCompatibility(p_field_volume); // may throw exception
1657     const FIELD<double> * p_field_size=p_field_volume;
1658     if(!p_field_volume) // if the user don't supply the volume
1659         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1660
1661     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1662     const double* vol=p_field_size->getValue(MED_EN::MED_FULL_INTERLACE);
1663     const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
1664     const T* value=getValue( MED_EN::MED_NO_INTERLACE); // get pointer to the field's values
1665
1666     double totVol=0.0;
1667     const double* p_vol=vol;
1668     for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
1669         totVol+=*p_vol;
1670
1671     double integrale=0.0;
1672     for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
1673         for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
1674             integrale += static_cast<double>((*value) * (*value)) * (*p_vol);
1675
1676     if(!p_field_volume) // if the user didn't supply the volume
1677         delete p_field_size; // delete temporary volume field
1678     if( totVol <= 0)
1679         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1680
1681     return integrale/totVol;
1682 }
1683
1684 /*!  Return L1 Norm  of the field's component.
1685  *   Cannot be applied to a field with a support on nodes.
1686  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1687  */
1688 template <class T> double FIELD<T>::normL1(int component, const FIELD<double> * p_field_volume) const
1689 {
1690     _checkNormCompatibility(p_field_volume); // may throw exception
1691     if ( component<1 || component>getNumberOfComponents() )
1692         throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
1693
1694     const FIELD<double> * p_field_size=p_field_volume;
1695     if(!p_field_volume) // if the user don't supply the volume
1696         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1697
1698     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1699     const double* vol=p_field_size->getValue(MED_EN::MED_FULL_INTERLACE);
1700     const T* value=getValueI( MED_EN::MED_NO_INTERLACE, component); // get pointer to the component's values
1701     const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
1702
1703     double integrale=0.0;
1704     double totVol=0.0;
1705     for (; value!=lastvalue ; ++value ,++vol)
1706     {
1707         integrale += std::abs( static_cast<double>(*value) ) * (*vol);
1708         totVol+=*vol;
1709     }
1710
1711     if(!p_field_volume) // if the user didn't supply the volume
1712         delete p_field_size; // delete temporary volume field
1713     if( totVol <= 0)
1714         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1715
1716     return integrale/totVol;
1717 }
1718
1719 /*!  Return L1 Norm  of the field.
1720  *   Cannot be applied to a field with a support on nodes.
1721  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1722  */
1723 template <class T> double FIELD<T>::normL1(const FIELD<double> * p_field_volume) const
1724 {
1725     _checkNormCompatibility(p_field_volume); // may throw exception
1726     const FIELD<double> * p_field_size=p_field_volume;
1727     if(!p_field_volume) // if the user don't supply the volume
1728         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1729
1730     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1731     const double* vol=p_field_size->getValue(MED_EN::MED_FULL_INTERLACE);
1732     const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
1733     const T* value=getValue( MED_EN::MED_NO_INTERLACE); // get pointer to the field's values
1734
1735     double totVol=0.0;
1736     const double* p_vol=vol;
1737     for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
1738         totVol+=*p_vol;
1739
1740     double integrale=0.0;
1741     for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
1742         for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
1743             integrale += std::abs( static_cast<double>(*value) ) * (*p_vol);
1744
1745     if(!p_field_volume) // if the user didn't supply the volume
1746         delete p_field_size; // delete temporary volume field
1747     if( totVol <= 0)
1748         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1749
1750     return integrale/totVol;
1751 }
1752
1753 /*! Return a new field (to deallocate with delete) lying on subSupport that is included by
1754  *   this->_support with corresponding values extracting from this->_value.
1755  */
1756 template <class T> FIELD<T>* FIELD<T>::extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
1757 {
1758   if(!subSupport->belongsTo(*_support))
1759     throw MEDEXCEPTION("FIELD<T>::extract : subSupport not included in this->_support !");
1760   if(_support->isOnAllElements() && subSupport->isOnAllElements())
1761     return new FIELD<T>(*this);
1762   FIELD<T> *ret=new FIELD<T>(subSupport,_numberOfComponents,MED_EN::MED_FULL_INTERLACE);
1763   if(!ret->_value)
1764     throw MEDEXCEPTION("FIELD<T>::extract : unvalid support detected !");
1765   T* valuesToSet=(T*)ret->_value->get(MED_EN::MED_FULL_INTERLACE);
1766   int nbOfEltsSub=subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
1767   const int *eltsSub=subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS);
1768   T* tempVals=new T[_numberOfComponents];
1769   for(int i=0;i<nbOfEltsSub;i++)
1770     {
1771       if(!getValueOnElement(eltsSub[i],tempVals))
1772         throw MEDEXCEPTION("Problem in belongsTo function !!!");
1773       for(int j=0;j<_numberOfComponents;j++)
1774         valuesToSet[i*_numberOfComponents+j]=tempVals[j];
1775     }
1776   delete [] tempVals;
1777   ret->setValueType(_valueType);
1778   ret->copyGlobalInfo(*this);
1779   return ret;
1780 }
1781
1782 /*!
1783   Constructor with parameters; the object is set via a file and its associated
1784   driver. For the moment only the MED_DRIVER is considered and if the last two
1785   argument (iterationNumber and orderNumber) are not set; their default value
1786   is -1. If the field fieldDriverName with the iteration number
1787   iterationNumber and the order number orderNumber does not exist in the file
1788   fieldDriverName; the constructor raises an exception.
1789 */
1790 template <class T> FIELD<T>::FIELD(const SUPPORT * Support,
1791                                    driverTypes driverType,
1792                                    const string & fileName/*=""*/,
1793                                    const string & fieldDriverName/*=""*/,
1794                                    const int iterationNumber,
1795                                    const int orderNumber)
1796   throw (MEDEXCEPTION)
1797 {
1798   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) : ";
1799
1800   int current;
1801
1802   BEGIN_OF(LOC);
1803
1804   init();
1805
1806   _support = Support;
1807   _value = (MEDARRAY<T>*)NULL;
1808
1809   _iterationNumber = iterationNumber;
1810   _time = 0.0;
1811   _orderNumber = orderNumber;
1812
1813   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::MED_LECT);
1814
1815 //   switch(driverType)
1816 //     {
1817 //     case MED_DRIVER :
1818 //       {
1819 //      MED_FIELD_RDONLY_DRIVER<T> myDriver(fileName,this);
1820 //      myDriver.setFieldName(fieldDriverName);
1821 //      current = addDriver(myDriver);
1822 //      break;
1823 //       }
1824 //   current = addDriver(driverType,fileName,fieldDriverName);
1825 //   switch(_drivers[current]->getAccessMode() ) {
1826 //   case MED_WRONLY : {
1827 //     MESSAGE("FIELD<T>::FIELD(driverTypes driverType, .....) : driverType must have a MED_RDONLY or MED_RDWR accessMode");
1828 //     rmDriver(current);
1829 //     break;}
1830 //   default : {
1831 //   }
1832 //   }
1833 //     default :
1834 //       {
1835 //      throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Driver type unknown : can't create driver!"));
1836 //      break;
1837 //       }
1838 //     }
1839
1840   _drivers[current]->open();
1841   _drivers[current]->read();
1842   _drivers[current]->close();
1843
1844   END_OF(LOC);
1845 }
1846
1847 /*!
1848   Destructor.
1849 */
1850 template <class T> FIELD<T>::~FIELD()
1851 {
1852   BEGIN_OF(" Destructeur FIELD<T>::~FIELD()");
1853   SCRUTE(this);
1854   if (_value) delete _value;
1855   END_OF(" Destructeur FIELD<T>::~FIELD()");
1856 }
1857
1858 /*!
1859
1860 */
1861 template <class T> void FIELD<T>::allocValue(const int NumberOfComponents)
1862 {
1863   const char* LOC = "FIELD<T>::allocValue(const int NumberOfComponents)" ;
1864   BEGIN_OF(LOC);
1865
1866   _numberOfComponents = NumberOfComponents ;
1867   if (_componentsTypes == NULL)
1868     _componentsTypes = new int[NumberOfComponents] ;
1869   if (_componentsNames == NULL)
1870     _componentsNames = new string[NumberOfComponents];
1871   if (_componentsDescriptions == NULL)
1872     _componentsDescriptions = new string[NumberOfComponents];
1873   if (_componentsUnits == NULL)
1874     _componentsUnits = new UNIT[NumberOfComponents];
1875   if (_MEDComponentsUnits == NULL)
1876     _MEDComponentsUnits = new string[NumberOfComponents];
1877   for (int i=0;i<NumberOfComponents;i++) {
1878     _componentsTypes[i] = 0 ;
1879   }
1880
1881   try {
1882     _numberOfValues = _support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
1883     MESSAGE(LOC <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
1884
1885     _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
1886
1887     _isRead = true ;
1888   }
1889   catch (MEDEXCEPTION &ex) {
1890     MESSAGE("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
1891     _value = (MEDARRAY<T>*)NULL ;
1892   }
1893
1894   SCRUTE(_value);
1895   END_OF("void FIELD<T>::allocValue(const int NumberOfComponents)");
1896 }
1897
1898 /*!
1899
1900 */
1901 template <class T> void FIELD<T>::allocValue(const int NumberOfComponents, const int LengthValue)
1902 {
1903   BEGIN_OF("void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)");
1904
1905   _numberOfComponents = NumberOfComponents ;
1906   if (_componentsTypes == NULL)
1907     _componentsTypes = new int[NumberOfComponents] ;
1908   if (_componentsNames == NULL)
1909     _componentsNames = new string[NumberOfComponents];
1910   if (_componentsDescriptions == NULL)
1911     _componentsDescriptions = new string[NumberOfComponents];
1912   if (_componentsUnits == NULL)
1913     _componentsUnits = new UNIT[NumberOfComponents];
1914   if (_MEDComponentsUnits == NULL)
1915     _MEDComponentsUnits = new string[NumberOfComponents];
1916   for (int i=0;i<NumberOfComponents;i++) {
1917     _componentsTypes[i] = 0 ;
1918   }
1919
1920   MESSAGE("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
1921   _numberOfValues = LengthValue ;
1922   _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
1923   _isRead = true ;
1924
1925   SCRUTE(_value);
1926   END_OF("void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)");
1927 }
1928
1929 /*!
1930
1931 */
1932 template <class T> void FIELD<T>::deallocValue()
1933 {
1934   BEGIN_OF("void FIELD<T>::deallocValue()");
1935   _numberOfValues = 0 ;
1936   _numberOfComponents = 0 ;
1937   if (_value != NULL)
1938     delete _value;
1939
1940   END_OF("void FIELD<T>::deallocValue()");
1941 }
1942
1943 // -----------------
1944 // Methodes Inline
1945 // -----------------
1946
1947 /*!
1948   Create the specified driver and return its index reference to path to
1949   read or write methods.
1950 */
1951
1952 template <class T> int FIELD<T>::addDriver(driverTypes driverType,
1953                                            const string & fileName/*="Default File Name.med"*/,
1954                                            const string & driverName/*="Default Field Name"*/,
1955                                            MED_EN::med_mode_acces access)
1956 {
1957   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) : ";
1958
1959   GENDRIVER * driver;
1960
1961   BEGIN_OF(LOC);
1962
1963   SCRUTE(driverType);
1964
1965   driver = DRIVERFACTORY::buildDriverForField(driverType,fileName,this,access);
1966
1967   _drivers.push_back(driver);
1968
1969   int current = _drivers.size()-1;
1970
1971   _drivers[current]->setFieldName(driverName);
1972
1973   END_OF(LOC);
1974
1975   return current;
1976 }
1977
1978
1979 /*!
1980   Duplicate the given driver and return its index reference to path to
1981   read or write methods.
1982 */
1983 template <class T> inline int FIELD<T>::addDriver (GENDRIVER & driver )
1984 {
1985   const char * LOC = "FIELD<T>::addDriver(GENDRIVER &) : ";
1986   int current;
1987
1988   BEGIN_OF(LOC);
1989
1990   // duplicate driver to delete it with destructor !
1991   GENDRIVER * newDriver = driver.copy() ;
1992
1993   _drivers.push_back(newDriver);
1994
1995   current = _drivers.size()-1;
1996   SCRUTE(current);
1997   driver.setId(current);
1998
1999   MESSAGE(LOC << " je suis la 1");
2000   END_OF(LOC);
2001   MESSAGE(LOC << " je suis la 2");
2002
2003   return current ;
2004 };
2005
2006 /*!
2007   Remove the driver referenced by its index.
2008 */
2009 template <class T> void FIELD<T>::rmDriver (int index/*=0*/)
2010 {
2011   const char * LOC = "FIELD<T>::rmDriver (int index=0): ";
2012   BEGIN_OF(LOC);
2013
2014   if ( _drivers[index] ) {
2015     //_drivers.erase(&_drivers[index]);
2016     // why not ????
2017     MESSAGE ("detruire");
2018   }
2019   else
2020     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2021                                      << "The <index given is invalid, index must be between  0 and  |"
2022                                      << _drivers.size()
2023                                      )
2024                           );
2025
2026   END_OF(LOC);
2027 }
2028
2029 /*!
2030   Read FIELD in the file specified in the driver given by its index.
2031 */
2032 template <class T> inline  void FIELD<T>::read(int index/*=0*/)
2033 {
2034   const char * LOC = "FIELD<T>::read(int index=0) : ";
2035   BEGIN_OF(LOC);
2036
2037   if ( _drivers[index] ) {
2038     _drivers[index]->open();
2039     _drivers[index]->read();
2040     _drivers[index]->close();
2041   }
2042   else
2043     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2044                                      << "The index given is invalid, index must be between  0 and |"
2045                                      << _drivers.size()
2046                                      )
2047                           );
2048   END_OF(LOC);
2049 }
2050
2051 /*!
2052   Write FIELD in the file specified in the driver given by its index.
2053 */
2054 template <class T> inline void FIELD<T>::write(int index/*=0*/, const string & driverName /*= ""*/)
2055 {
2056   const char * LOC = "FIELD<T>::write(int index=0, const string & driverName = \"\") : ";
2057   BEGIN_OF(LOC);
2058
2059   if( _drivers[index] ) {
2060     _drivers[index]->open();
2061     if (driverName != "") _drivers[index]->setFieldName(driverName);
2062     _drivers[index]->write();
2063     _drivers[index]->close();
2064   }
2065   else
2066     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2067                                      << "The index given is invalid, index must be between  0 and |"
2068                                      << _drivers.size()
2069                                      )
2070                           );
2071   END_OF(LOC);
2072 }
2073
2074 /*!
2075   Write FIELD in the file specified in the driver given by its index. Use this
2076   method for ASCII drivers (e.g. VTK_DRIVER)
2077 */
2078 template <class T> inline void FIELD<T>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
2079 {
2080   const char * LOC = "FIELD<T>::write(int index=0, const string & driverName = \"\") : ";
2081   BEGIN_OF(LOC);
2082
2083   if( _drivers[index] ) {
2084     _drivers[index]->openAppend();
2085     if (driverName != "") _drivers[index]->setFieldName(driverName);
2086     _drivers[index]->writeAppend();
2087     _drivers[index]->close();
2088   }
2089   else
2090     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2091                                      << "The index given is invalid, index must be between  0 and |"
2092                                      << _drivers.size()
2093                                      )
2094                           );
2095   END_OF(LOC);
2096 }
2097
2098 /*!
2099   \internal
2100   Write FIELD with the driver which is equal to the given driver.
2101
2102   Use by MED object.
2103 */
2104 template <class T> inline void FIELD<T>::write(const GENDRIVER & genDriver)
2105 {
2106   const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
2107   BEGIN_OF(LOC);
2108
2109   for (unsigned int index=0; index < _drivers.size(); index++ )
2110     if ( *_drivers[index] == genDriver ) {
2111       _drivers[index]->open();
2112       _drivers[index]->write();
2113       _drivers[index]->close();
2114     }
2115
2116   END_OF(LOC);
2117
2118 }
2119
2120 /*!
2121   \internal
2122   Write FIELD with the driver which is equal to the given driver.
2123
2124   Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
2125 */
2126 template <class T> inline void FIELD<T>::writeAppend(const GENDRIVER & genDriver)
2127 {
2128   const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
2129   BEGIN_OF(LOC);
2130
2131   for (unsigned int index=0; index < _drivers.size(); index++ )
2132     if ( *_drivers[index] == genDriver ) {
2133       _drivers[index]->openAppend();
2134       _drivers[index]->writeAppend();
2135       _drivers[index]->close();
2136     }
2137
2138   END_OF(LOC);
2139
2140 }
2141
2142 /*!
2143   \internal
2144   Read FIELD with the driver which is equal to the given driver.
2145
2146   Use by MED object.
2147 */
2148 template <class T> inline void FIELD<T>::read(const GENDRIVER & genDriver)
2149 {
2150   const char * LOC = " FIELD<T>::read(const GENDRIVER &) : ";
2151   BEGIN_OF(LOC);
2152
2153   for (unsigned int index=0; index < _drivers.size(); index++ )
2154     if ( *_drivers[index] == genDriver ) {
2155       _drivers[index]->open();
2156       _drivers[index]->read();
2157       _drivers[index]->close();
2158     }
2159
2160   END_OF(LOC);
2161
2162 }
2163
2164 /*!
2165   Fills in already allocated retValues array the values related to eltIdInSup.
2166   If the element does not exist in this->_support false is returned, true otherwise.
2167  */
2168 template <class T> bool FIELD<T>::getValueOnElement(int eltIdInSup,T* retValues) const
2169 {
2170   if(eltIdInSup<1)
2171     return false;
2172   if(_support->isOnAllElements())
2173     {
2174       int nbOfEltsThis=_support->getMesh()->getNumberOfElements(_support->getEntity(),MED_EN::MED_ALL_ELEMENTS);
2175       if(eltIdInSup>nbOfEltsThis)
2176         return false;
2177       const T* valsThis=getValue(MED_EN::MED_FULL_INTERLACE);
2178       for(int j=0;j<_numberOfComponents;j++)
2179         retValues[j]=valsThis[(eltIdInSup-1)*_numberOfComponents+j];
2180       return true;
2181     }
2182   else
2183     {
2184       int nbOfEltsThis=_support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2185       const int *eltsThis=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
2186       int iThis;
2187       bool found=false;
2188       for(iThis=0;iThis<nbOfEltsThis && !found;)
2189         if(eltsThis[iThis]==eltIdInSup)
2190           found=true;
2191         else
2192           iThis++;
2193       if(!found)
2194         return false;
2195       const T* valsThis=getValue(MED_EN::MED_FULL_INTERLACE);
2196       for(int j=0;j<_numberOfComponents;j++)
2197         retValues[j]=valsThis[iThis*_numberOfComponents+j];
2198       return true;
2199     }
2200 }
2201
2202 /*!
2203   \if developper
2204   Destroy the MEDARRAY<T> in FIELD and put the new one without copy.
2205   \endif
2206 */
2207 template <class T> inline void FIELD<T>::setValue(MEDARRAY<T> *Value)
2208 {
2209   if (NULL != _value) delete _value ;
2210   _value=Value ;
2211 }
2212
2213 /*!
2214   \if developper
2215   Return a reference to  the MEDARRAY<T> in FIELD.
2216   \endif
2217 */
2218 template <class T> inline MEDARRAY<T>* FIELD<T>::getvalue() const
2219 {
2220   return _value ;
2221 }
2222
2223 /*!
2224   Return the actual length of the reference to values array returned by getValue.
2225 */
2226 template <class T> inline int FIELD<T>::getValueLength(MED_EN::medModeSwitch Mode) const{
2227   return _numberOfComponents*_numberOfValues;
2228 }
2229
2230 /*!
2231   Return a reference to values array to read them.
2232 */
2233 template <class T> inline const T* FIELD<T>::getValue(MED_EN::medModeSwitch Mode) const
2234 {
2235   return _value->get(Mode) ;
2236 }
2237
2238 /*!
2239   Return a reference to i^{th} row or column - component - (depend on Mode value)
2240   of FIELD values array.
2241 */
2242 template <class T> inline const T* FIELD<T>::getValueI(MED_EN::medModeSwitch Mode,int i) const
2243 {
2244  if ( Mode == MED_EN::MED_FULL_INTERLACE )
2245  {
2246          return _value->getRow(i) ;
2247  }
2248  ASSERT (  Mode == MED_EN::MED_NO_INTERLACE);
2249  return _value->getColumn(i);
2250 }
2251
2252 /*!
2253   Return the value of i^{th} element and j^{th} component.
2254 */
2255 template <class T> inline T FIELD<T>::getValueIJ(int i,int j) const
2256 {
2257   return _value->getIJ(i,j) ;
2258 }
2259
2260 /*!
2261   Copy new values array in FIELD according to the given mode.
2262
2263   Array must have right size. If not results are unpredicable.
2264 */
2265 template <class T> inline void FIELD<T>::setValue(MED_EN::medModeSwitch mode, T* value)
2266 {
2267   _value->set(mode,value);
2268 }
2269
2270 /*!
2271   Update values array in FIELD with the given ones according to specified mode.
2272 */
2273 template <class T> inline void FIELD<T>::setValueI(MED_EN::medModeSwitch mode, int i, T* value)
2274 {
2275   // PROVISOIRE :
2276   if (MED_EN::MED_FULL_INTERLACE == mode)
2277     _value->setI(i,value);
2278   else if (MED_EN::MED_NO_INTERLACE == mode)
2279     _value->setJ(i,value);
2280   else
2281     throw MEDEXCEPTION(LOCALIZED("FIELD<T>::setValueI : bad medModeSwitch")) ;
2282 }
2283
2284 /*!
2285   Set the value of i^{th} element and j^{th} component with the given one.
2286 */
2287 template <class T> inline void FIELD<T>::setValueIJ(int i, int j, T value)
2288 {
2289   _value->setIJ(i,j,value);
2290 }
2291
2292 /*
2293   METHODS
2294 */
2295
2296 /*!
2297   Fill values array with volume values.
2298 */
2299 template <class T> void FIELD<T>::getVolume() const throw (MEDEXCEPTION)
2300 {
2301   const char * LOC = "FIELD<double>::getVolume() const : ";
2302   BEGIN_OF(LOC);
2303
2304   // The field has to be initilised by a non empty support and a
2305   // number of components = 1 and its value type has to be set to MED_REEL64
2306   // (ie a FIELD<double>)
2307
2308   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
2309       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"));
2310
2311   END_OF(LOC);
2312 }
2313
2314 /*!
2315   Fill values array with area values.
2316 */
2317 template <class T> void FIELD<T>::getArea() const throw (MEDEXCEPTION)
2318 {
2319   const char * LOC = "FIELD<double>::getArea() const : ";
2320   BEGIN_OF(LOC);
2321
2322   // The field has to be initilised by a non empty support and a
2323   // number of components = 1 and its value type has to be set to MED_REEL64
2324   // (ie a FIELD<double>)
2325
2326   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
2327       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"));
2328
2329   END_OF(LOC);
2330 }
2331
2332 /*!
2333   Fill values array with length values.
2334 */
2335 template <class T> void FIELD<T>::getLength() const throw (MEDEXCEPTION)
2336 {
2337   const char * LOC = "FIELD<double>::getLength() const : ";
2338   BEGIN_OF(LOC);
2339
2340   // The field has to be initilised by a non empty support and a
2341   // number of components = 1 and its value type has to be set to MED_REEL64
2342   // (ie a FIELD<double>)
2343
2344   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
2345       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"));
2346
2347   END_OF(LOC);
2348 }
2349
2350 /*!
2351   Fill values array with normal values.
2352 */
2353 template <class T> void FIELD<T>::getNormal() const throw (MEDEXCEPTION)
2354 {
2355   const char * LOC = "FIELD<double>::getNormal() const : ";
2356   BEGIN_OF(LOC);
2357
2358   // The field has to be initilised by a non empty support and a
2359   // number of components = 1 and its value type has to be set to MED_REEL64
2360   // (ie a FIELD<double>)
2361
2362   if (_support == (SUPPORT *) NULL)
2363       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"));
2364
2365   int dim_space = _support->getMesh()->getSpaceDimension();
2366
2367   if ((_numberOfComponents != dim_space) || (_valueType != MED_EN::MED_REEL64))
2368       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"));
2369
2370   END_OF(LOC);
2371 }
2372
2373 /*!
2374   Fill values array with barycenter values.
2375 */
2376 template <class T> void FIELD<T>::getBarycenter() const throw (MEDEXCEPTION)
2377 {
2378   const char * LOC = "FIELD<double>::getBarycenter() const : ";
2379   BEGIN_OF(LOC);
2380
2381   // The field has to be initilised by a non empty support and a number of
2382   //components = space dimension and its value type has to be set to MED_REEL64
2383   // (ie a FIELD<double>)
2384
2385   if (_support == (SUPPORT *) NULL)
2386       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"));
2387
2388   int dim_space = _support->getMesh()->getSpaceDimension();
2389
2390   if ((_numberOfComponents != dim_space) || (_valueType != MED_EN::MED_REEL64))
2391       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"));
2392
2393   END_OF(LOC);
2394 }
2395
2396 /*!
2397   Fill array by using T_Analytic.
2398   WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
2399   Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
2400  */
2401 template <class T>
2402 template<void T_Analytic(const double *,T*)>
2403 void FIELD<T>::fillFromAnalytic()
2404 {
2405   int i,j;
2406   MESH * mesh = _support->getMesh();
2407   int spaceDim = mesh->getSpaceDimension();
2408   const double * coord;
2409   FIELD<double> * barycenterField=0;
2410   double ** xyz=new double* [spaceDim];
2411   bool deallocateXyz=false;
2412   if(_support->getEntity()==MED_EN::MED_NODE)
2413     {
2414       if (_support->isOnAllElements())
2415         {
2416           coord=mesh->getCoordinates(MED_EN::MED_NO_INTERLACE);
2417           for(i=0; i<spaceDim; i++)
2418             xyz[i]=(double *)coord+i*_numberOfValues;
2419         }
2420       else
2421         {
2422           coord = mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
2423           const int * nodesNumber=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
2424           for(i=0; i<spaceDim; i++)
2425             xyz[i]=new double[_numberOfValues];
2426           deallocateXyz=true;
2427           for(i=0;i<_numberOfValues;i++)
2428             {
2429               for(j=0;j<spaceDim;j++)
2430                 xyz[j][i]=coord[(nodesNumber[i]-1)*spaceDim+j];
2431             }
2432         }
2433     }
2434   else
2435     {
2436       barycenterField = mesh->getBarycenter(_support);
2437       coord=barycenterField->getValue(MED_EN::MED_NO_INTERLACE);
2438       for(i=0; i<spaceDim; i++)
2439         xyz[i]=(double *)(coord+i*_numberOfValues);
2440     }
2441   T* valsToSet=(T*)getValue(MED_EN::MED_FULL_INTERLACE);
2442   double *temp=new double[spaceDim];
2443   for(i=0;i<_numberOfValues;i++)
2444   {
2445     for(j=0;j<spaceDim;j++)
2446       temp[j]=xyz[j][i];
2447     T_Analytic(temp,valsToSet+i*_numberOfComponents);
2448   }
2449   delete [] temp;
2450   if(barycenterField)
2451     delete barycenterField;
2452   if(deallocateXyz)
2453     for(j=0;j<spaceDim;j++)
2454       delete [] xyz[j];
2455   delete [] xyz;
2456 }
2457
2458 }//End namespace MEDMEM
2459
2460 #endif /* FIELD_HXX */