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