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