src/MEDMEM/MEDMEM_Field.hxx

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