src/MEDMEM/MEDMEM_Field.hxx

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