src/MEDMEM/MEDMEM_Field.hxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 /*
  21  File Field.hxx
  22  $Header$
  23 */
  24
  25 #ifndef FIELD_HXX
  26 #define FIELD_HXX
  27
  28 #include "MEDMEM.hxx"
  29
  30 #include <vector>
  31 #include <map>
  32 #include <algorithm>
  33 #include <math.h>
  34 #include <cmath>
  35
  36 #include "MEDMEM_Utilities.hxx"
  37 #include "MEDMEM_Exception.hxx"
  38 #include "MEDMEM_define.hxx"
  39 #include "MEDMEM_Support.hxx"
  40 #include "MEDMEM_Unit.hxx"
  41 #include "MEDMEM_nArray.hxx"
  42 #include "MEDMEM_GenDriver.hxx"
  43 #include "MEDMEM_ArrayInterface.hxx"
  44 #include "MEDMEM_SetInterlacingType.hxx"
  45 #include "MEDMEM_FieldForward.hxx"
  46 #include "MEDMEM_GaussLocalization.hxx"
  47
  48 #define DBL_MAX 1.0E+308
  49 #define DBL_MIN -1.0E+308
  50
  51 /*!
  52
  53   This class contains all the informations related with a template class FIELD :
  54   - Components descriptions
  55   - Time step description
  56   - Location of the values (a SUPPORT class)
  57
  58 */
  59
  60 namespace MEDMEM {
  61
  62   template<class T>
  63   struct MinMax {
  64   };
  65
  66   template<>
  67   struct MinMax<double> {
  68     static double getMin() { return DBL_MIN; }
  69     static double getMax() { return DBL_MAX; }
  70   };
  71
  72   template<>
  73   struct MinMax<int> {
  74     static int getMin() { return INT_MIN; }
  75     static int getMax() { return INT_MAX; }
  76   };
  77
  78   template < typename T > struct SET_VALUE_TYPE {
  79     static const MED_EN::med_type_champ _valueType = MED_EN::MED_UNDEFINED_TYPE;};
  80   template < > struct SET_VALUE_TYPE<double> {
  81     static const MED_EN::med_type_champ _valueType = MED_EN::MED_REEL64; };
  82   template < > struct SET_VALUE_TYPE<int> {
  83     static const MED_EN::med_type_champ _valueType = MED_EN::MED_INT32; };
  84
  85 class MEDMEM_EXPORT FIELD_    // GENERIC POINTER TO a template <class T, class INTERLACING_TAG> class FIELD
  86 {
  87 protected:
  88
  89   bool            _isRead ;
  90   bool            _isMinMax;
  91
  92   /*!
  93     \if developper
  94     Field name.
  95     \endif
  96   */
  97   string          _name ;
  98   /*!
  99     \if developper
 100     Field description.
 101     \endif
 102   */
 103   string          _description ;
 104   /*!
 105     \if developper
 106     Pointer to the support the field deals with.
 107     \endif
 108   */
 109   const SUPPORT * _support ;
 110
 111   /*!
 112     \if developper
 113     Number of field's components.
 114     \endif
 115   */
 116   int             _numberOfComponents ;
 117   /*!
 118     \if developper
 119     Number of field's values.
 120     doesn't take care of _numberOfComponents
 121     and number of Gauss points.
 122     \endif
 123   */
 124   int             _numberOfValues ;
 125
 126   /*!
 127     \if developper
 128     Array of size _numberOfComponents. \n
 129     (constant, scalar, vector, tensor)\n
 130     We could use an array of integer to store
 131     numbers of values: \n
 132     - 1 for scalar,\n
 133     - space dimension for vector,\n
 134     - space dimension square for tensor.\n
 135     So numbers of values per entities would be
 136     sum of _componentsTypes array.
 137
 138     Not implemented yet! All type are scalar !
 139     \endif
 140   */
 141   int *           _componentsTypes ;
 142   /*!
 143     \if developper
 144     Array of size _numberOfComponents
 145     storing components names if any.
 146     \endif
 147   */
 148   string *        _componentsNames;
 149   /*!
 150     \if developper
 151     Array of size _numberOfComponents
 152     storing components descriptions if any.
 153     \endif
 154  */
 155   string *        _componentsDescriptions;
 156   /*!
 157     \if developper
 158     Array of size _numberOfComponents
 159     storing components units if any.
 160     \endif
 161  */
 162   UNIT *          _componentsUnits;
 163   /*!
 164     \if developper
 165     Array of size _numberOfComponents
 166     storing components units if any.
 167     \endif
 168   */
 169   string *        _MEDComponentsUnits;
 170   /*!
 171     \if developper
 172     Iteration number of the field.
 173     \endif
 174   */
 175   int             _iterationNumber ;
 176   /*!
 177     \if developper
 178     Time of the field.
 179     \endif
 180   */
 181   double          _time;
 182   /*!
 183     \if developper
 184     Order number of the field.
 185     \endif
 186   */
 187   int             _orderNumber ;
 188   /*!
 189     \if developper
 190     At the initialization step of the field using the constructors; this attribute,
 191     the value type (integer or real) , is set automatically. There is a get method
 192     but not a set method for this attribute.
 193     \endif
 194   */
 195   MED_EN::med_type_champ _valueType ;
 196   /*!
 197     \if developper
 198     At the initialization step of the field using the constructors; this attribute,
 199     the interlacing type (full interlace or no interlace field value storage), is set
 200     automatically. There is a get method but not a set method for this attribute.
 201     \endif
 202   */
 203   MED_EN::medModeSwitch _interlacingType;
 204
 205   vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
 206   static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
 207   static void _deepCheckFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
 208   void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL) const  throw (MEDEXCEPTION);
 209   FIELD<double>* _getFieldSize() const;
 210
 211 public:
 212
 213   friend class MED_MED_RDONLY_DRIVER21;
 214   friend class MED_MED_WRONLY_DRIVER21;
 215   friend class MED_MED_RDWR_DRIVER21;
 216   friend class MED_MED_RDONLY_DRIVER22;
 217   friend class MED_MED_WRONLY_DRIVER22;
 218   friend class MED_MED_RDWR_DRIVER22;
 219   friend class VTK_MED_DRIVER;
 220
 221   /*!
 222     Constructor.
 223   */
 224   FIELD_ ();
 225   /*!
 226     Constructor.
 227   */
 228   FIELD_(const SUPPORT * Support, const int NumberOfComponents);
 229   /*!
 230     Copy constructor.
 231   */
 232   FIELD_(const FIELD_ &m);
 233
 234   /*!
 235     Destructor.
 236   */
 237   virtual ~FIELD_();
 238
 239  FIELD_& operator=(const FIELD_ &m);
 240
 241   virtual  void     rmDriver(int index=0);
 242   virtual   int     addDriver(driverTypes driverType,
 243                               const string & fileName="Default File Name.med",
 244                               const string & driverFieldName="Default Field Nam",
 245                               MED_EN::med_mode_acces access=MED_EN::MED_REMP) ;
 246
 247   virtual  int      addDriver( GENDRIVER & driver);
 248   virtual  void     read (const GENDRIVER &);
 249   virtual  void     read(int index=0);
 250   virtual  void     openAppend( void );
 251   virtual  void     write(const GENDRIVER &);
 252   virtual  void     write(int index=0, const string & driverName="");
 253
 254   virtual  void     writeAppend(const GENDRIVER &);
 255   virtual  void     writeAppend(int index=0, const string & driverName="");
 256
 257   inline void     setName(const string Name);
 258   inline string   getName() const;
 259   inline void     setDescription(const string Description);
 260   inline string   getDescription() const;
 261   inline const SUPPORT * getSupport() const;
 262   inline void     setSupport(const SUPPORT * support);
 263   inline void     setNumberOfComponents(const int NumberOfComponents);
 264   inline int      getNumberOfComponents() const;
 265   inline void     setNumberOfValues(const int NumberOfValues);
 266   inline int      getNumberOfValues() const;
 267   //    inline void     setComponentType(int *ComponentType);
 268   //    inline int *    getComponentType() const;
 269   //    inline int      getComponentTypeI(int i) const;
 270   inline void     setComponentsNames(const string * ComponentsNames);
 271   inline void     setComponentName(int i, const string ComponentName);
 272   inline const string * getComponentsNames() const;
 273   inline string   getComponentName(int i) const;
 274   inline void     setComponentsDescriptions(const string * ComponentsDescriptions);
 275   inline void     setComponentDescription(int i, const string ComponentDescription);
 276   inline const string * getComponentsDescriptions() const;
 277   inline string   getComponentDescription(int i) const;
 278
 279   // provisoire : en attendant de regler le probleme des unites !
 280   inline void     setComponentsUnits(const UNIT * ComponentsUnits);
 281   inline const UNIT *   getComponentsUnits() const;
 282   inline const UNIT *   getComponentUnit(int i) const;
 283   inline void     setMEDComponentsUnits(const string * MEDComponentsUnits);
 284   inline void     setMEDComponentUnit(int i, const string MEDComponentUnit);
 285   inline const string * getMEDComponentsUnits() const;
 286   inline string   getMEDComponentUnit(int i) const;
 287
 288   inline void     setIterationNumber(int IterationNumber);
 289   inline int      getIterationNumber() const;
 290   inline void     setTime(double Time);
 291   inline double   getTime() const;
 292   inline void     setOrderNumber(int OrderNumber);
 293   inline int      getOrderNumber() const;
 294
 295   inline MED_EN::med_type_champ getValueType () const;
 296   inline MED_EN::medModeSwitch  getInterlacingType() const;
 297   virtual inline bool getGaussPresence() const throw (MEDEXCEPTION);
 298 protected:
 299   void copyGlobalInfo(const FIELD_& m);
 300 };
 301
 302 // ---------------------------------
 303 // Implemented Methods : constructor
 304 // ---------------------------------
 305
 306 // -----------------
 307 // Methodes Inline
 308 // -----------------
 309 /*!
 310   Set FIELD name.
 311 */
 312 inline void FIELD_::setName(const string Name)
 313 {
 314   _name=Name;
 315 }
 316 /*!
 317   Get FIELD name.
 318 */
 319 inline string FIELD_::getName() const
 320 {
 321   return _name;
 322 }
 323 /*!
 324   Set FIELD description.
 325 */
 326 inline void FIELD_::setDescription(const string Description)
 327 {
 328   _description=Description;
 329 }
 330 /*!
 331   Get FIELD description.
 332 */
 333 inline string FIELD_::getDescription() const
 334 {
 335   return _description;
 336 }
 337 /*!
 338   Set FIELD number of components.
 339 */
 340 inline void FIELD_::setNumberOfComponents(const int NumberOfComponents)
 341 {
 342   _numberOfComponents=NumberOfComponents;
 343 }
 344 /*!
 345   Get FIELD number of components.
 346 */
 347 inline int FIELD_::getNumberOfComponents() const
 348 {
 349   return _numberOfComponents ;
 350 }
 351 /*!
 352   Set FIELD number of values.
 353
 354   It must be the same than in the associated SUPPORT object.
 355 */
 356 inline void FIELD_::setNumberOfValues(const int NumberOfValues)
 357 {
 358   _numberOfValues=NumberOfValues;
 359 }
 360 /*!
 361   Get FIELD number of value.
 362 */
 363 inline int FIELD_::getNumberOfValues() const
 364 {
 365   return _numberOfValues ;
 366 }
 367
 368 //  inline void FIELD_::setComponentType(int *ComponentType)
 369 //  {
 370 //    _componentsTypes=ComponentType ;
 371 //  }
 372 //  inline int * FIELD_::getComponentType() const
 373 //  {
 374 //    return _componentsTypes ;
 375 //  }
 376 //  inline int FIELD_::getComponentTypeI(int i) const
 377 //  {
 378 //    return _componentsTypes[i-1] ;
 379 //  }
 380
 381 /*!
 382   Set FIELD components names.
 383
 384   Duplicate the ComponentsNames string array to put components names in
 385   FIELD. ComponentsNames size must be equal to number of components.
 386 */
 387 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
 388 {
 389   if (NULL == _componentsNames)
 390     _componentsNames = new string[_numberOfComponents] ;
 391   for (int i=0; i<_numberOfComponents; i++)
 392     _componentsNames[i]=ComponentsNames[i] ;
 393 }
 394 /*!
 395   Set FIELD i^th component name.
 396
 397   i must be >=1 and <= number of components.
 398 */
 399 inline void FIELD_::setComponentName(int i, const string ComponentName)
 400 {
 401   _componentsNames[i-1]=ComponentName ;
 402 }
 403 /*!
 404   Get a reference to the string array which contain the components names.
 405
 406   This Array size is equal to number of components
 407 */
 408 inline const string * FIELD_::getComponentsNames() const
 409 {
 410   return _componentsNames ;
 411 }
 412 /*!
 413   Get the name of the i^th component.
 414 */
 415 inline string FIELD_::getComponentName(int i) const
 416 {
 417   return _componentsNames[i-1] ;
 418 }
 419 /*!
 420   Set FIELD components descriptions.
 421
 422   Duplicate the ComponentsDescriptions string array to put components
 423   descriptions in FIELD.
 424   ComponentsDescriptions size must be equal to number of components.
 425 */
 426 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
 427 {
 428   if (NULL == _componentsDescriptions)
 429     _componentsDescriptions = new string[_numberOfComponents] ;
 430   for (int i=0; i<_numberOfComponents; i++)
 431     _componentsDescriptions[i]=ComponentsDescriptions[i] ;
 432 }
 433 /*!
 434   Set FIELD i^th component description.
 435
 436   i must be >=1 and <= number of components.
 437 */
 438 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
 439 {
 440   _componentsDescriptions[i-1]=ComponentDescription ;
 441 }
 442 /*!
 443   Get a reference to the string array which contain the components descriptions.
 444
 445   This Array size is equal to number of components
 446 */
 447 inline const string * FIELD_::getComponentsDescriptions() const
 448 {
 449   return _componentsDescriptions ;
 450 }
 451 /*!
 452   Get the description of the i^th component.
 453 */
 454 inline string FIELD_::getComponentDescription(int i) const
 455 {
 456   return _componentsDescriptions[i-1];
 457 }
 458
 459 /*!
 460   \todo
 461   Set FIELD components UNIT.
 462
 463   Duplicate the ComponentsUnits UNIT array to put components
 464   units in FIELD.
 465   ComponentsUnits size must be equal to number of components.
 466 */
 467 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
 468 {
 469   if (NULL == _componentsUnits)
 470     _componentsUnits = new UNIT[_numberOfComponents] ;
 471   for (int i=0; i<_numberOfComponents; i++)
 472     _componentsUnits[i]=ComponentsUnits[i] ;
 473 }
 474 /*!
 475   Get a reference to the UNIT array which contain the components units.
 476
 477   This Array size is equal to number of components
 478 */
 479 inline const UNIT * FIELD_::getComponentsUnits() const
 480 {
 481   return _componentsUnits ;
 482 }
 483 /*!
 484   Get the UNIT of the i^th component.
 485 */
 486 inline const UNIT * FIELD_::getComponentUnit(int i) const
 487 {
 488   return &_componentsUnits[i-1] ;
 489 }
 490 /*!
 491   Set FIELD components unit.
 492
 493   Duplicate the MEDComponentsUnits string array to put components
 494   units in FIELD.
 495   MEDComponentsUnits size must be equal to number of components.
 496
 497 */
 498 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
 499 {
 500   if (NULL == _MEDComponentsUnits)
 501     _MEDComponentsUnits = new string[_numberOfComponents] ;
 502   for (int i=0; i<_numberOfComponents; i++)
 503     _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
 504 }
 505 /*!
 506   Set FIELD i^th component unit.
 507
 508   i must be >=1 and <= number of components.
 509 */
 510 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
 511 {
 512   _MEDComponentsUnits[i-1]=MEDComponentUnit ;
 513 }
 514 /*!
 515   Get a reference to the string array which contain the components units.
 516
 517   This Array size is equal to number of components
 518 */
 519 inline const string * FIELD_::getMEDComponentsUnits() const
 520 {
 521   return _MEDComponentsUnits ;
 522 }
 523 /*!
 524   Get the string for unit of the i^th component.
 525 */
 526 inline string FIELD_::getMEDComponentUnit(int i) const
 527 {
 528   return _MEDComponentsUnits[i-1] ;
 529 }
 530 /*!
 531   Set the iteration number where FIELD has been calculated.
 532 */
 533 inline void FIELD_::setIterationNumber(int IterationNumber)
 534 {
 535   _iterationNumber=IterationNumber;
 536 }
 537 /*!
 538   Get the iteration number where FIELD has been calculated.
 539 */
 540 inline int FIELD_::getIterationNumber() const
 541 {
 542   return _iterationNumber ;
 543 }
 544 /*!
 545   Set the time (in second) where FIELD has been calculated.
 546 */
 547 inline void FIELD_::setTime(double Time)
 548 {
 549   _time=Time ;
 550 }
 551 /*!
 552   Get the time (in second) where FIELD has been calculated.
 553 */
 554 inline double FIELD_::getTime() const
 555 {
 556   return _time ;
 557 }
 558 /*!
 559   Set the order number where FIELD has been calculated.
 560
 561   It corresponds to internal iteration during one time step.
 562 */
 563 inline void FIELD_::setOrderNumber(int OrderNumber)
 564 {
 565   _orderNumber=OrderNumber ;
 566 }
 567 /*!
 568   Get the order number where FIELD has been calculated.
 569 */
 570 inline int FIELD_::getOrderNumber() const
 571 {
 572   return _orderNumber ;
 573 }
 574 /*!
 575   Get a reference to the SUPPORT object associated to FIELD.
 576 */
 577 inline  const SUPPORT * FIELD_::getSupport() const
 578 {
 579   return _support ;
 580 }
 581 /*!
 582   Set the reference to the SUPPORT object associated to FIELD.
 583
 584   Reference is not duplicate, so it must not be deleted.
 585 */
 586 inline void FIELD_::setSupport(const SUPPORT * support)
 587 {
 588   //A.G. Addings for RC
 589   if(_support)
 590     _support->removeReference();
 591   _support = support ;
 592   if(_support)
 593     _support->addReference();
 594 }
 595 /*!
 596   Get the FIELD med value type (MED_INT32 or MED_REEL64).
 597 */
 598 inline MED_EN::med_type_champ FIELD_::getValueType () const
 599 {
 600   return _valueType ;
 601 }
 602
 603 /*!
 604   Get the FIELD med interlacing type (MED_FULL_INTERLACE, MED_NO_INTERLACE or MED_NO_INTERLACE_BY_TYPE).
 605 */
 606   inline MED_EN::medModeSwitch FIELD_::getInterlacingType () const
 607 {
 608   return _interlacingType ;
 609 }
 610
 611 /*!
 612   Get the FIELD gauss presence.
 613 */
 614   inline bool  FIELD_::getGaussPresence() const throw (MEDEXCEPTION)
 615 {
 616   const char * LOC = "FIELD_::getGaussPresence() : ";
 617   throw MEDEXCEPTION(STRING(LOC) << " This FIELD_ doesn't rely on a FIELD<T>" );
 618 }
 619
 620 } //End namespace MEDMEM
 621
 622 /////////////////////////
 623 // END OF CLASS FIELD_ //
 624 /////////////////////////
 625
 626 /*!
 627
 628   This template class contains informations related with a FIELD :
 629   - Values of the field, their type (real or integer), the storage mode (full interlace,
 630     no interlace or no interlace by type).
 631
 632 */
 633
 634 namespace MEDMEM {
 635
 636   template<class T2> class MED_FIELD_RDONLY_DRIVER21;
 637   template<class T2> class MED_FIELD_WRONLY_DRIVER21;
 638   template<class T2> class MED_FIELD_RDONLY_DRIVER22;
 639   template<class T2> class MED_FIELD_WRONLY_DRIVER22;
 640   template<class T2> class VTK_FIELD_DRIVER;
 641
 642   template <class T,
 643             class INTERLACING_TAG
 644             > class FIELD : public FIELD_
 645 {
 646 protected:
 647
 648   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array   ArrayNoGauss;
 649   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array     ArrayGauss;
 650   typedef typename MEDMEM_ArrayInterface<T,NoInterlace,NoGauss>::Array       ArrayNo;
 651   typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array     ArrayFull;
 652   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,NoGauss>::Array ArrayNoByType;
 653   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,Gauss>::Array   ArrayNoByTypeGauss;
 654   typedef MEDMEM_Array_ Array;
 655   typedef T ElementType;
 656   typedef INTERLACING_TAG InterlacingTag;
 657   typedef map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> locMap;
 658
 659   // array of value of type T
 660   Array *_value ;
 661
 662   // extrema values
 663   T _vmin;
 664   T _vmax;
 665
 666   map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> _gaussModel; //A changer quand les drivers seront template de l'entrelacement
 667
 668   static T _scalarForPow;
 669   static T pow(T x);
 670
 671 private:
 672   void _operation(const FIELD& m,const FIELD& n, char* Op);
 673   void _operationInitialize(const FIELD& m,const FIELD& n, char* Op);
 674   void _add_in_place(const FIELD& m,const FIELD& n);
 675   void _sub_in_place(const FIELD& m,const FIELD& n);
 676   void _mul_in_place(const FIELD& m,const FIELD& n);
 677   void _div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION);
 678   //setValueType() ;
 679
 680 public:
 681   FIELD();
 682   FIELD(const FIELD &m);
 683   FIELD(const SUPPORT * Support, const int NumberOfComponents)  throw (MEDEXCEPTION) ;
 684   FIELD( driverTypes driverType,
 685          const string & fileName, const string & fieldDriverName,
 686          const int iterationNumber=-1, const int orderNumber=-1)
 687     throw (MEDEXCEPTION);
 688   FIELD(const SUPPORT * Support, driverTypes driverType,
 689         const string & fileName="", const string & fieldName="",
 690         const int iterationNumber = -1, const int orderNumber = -1)
 691     throw (MEDEXCEPTION);
 692   ~FIELD();
 693
 694   FIELD & operator=(const FIELD &m);
 695
 696   const FIELD operator+(const FIELD& m) const;
 697   const FIELD operator-(const FIELD& m) const;
 698   const FIELD operator*(const FIELD& m) const;
 699   const FIELD operator/(const FIELD& m) const;
 700   const FIELD operator-() const;
 701   FIELD& operator+=(const FIELD& m);
 702   FIELD& operator-=(const FIELD& m);
 703   FIELD& operator*=(const FIELD& m);
 704   FIELD& operator/=(const FIELD& m);
 705   static FIELD* add(const FIELD& m, const FIELD& n);
 706   static FIELD* addDeep(const FIELD& m, const FIELD& n);
 707   static FIELD* sub(const FIELD& m, const FIELD& n);
 708   static FIELD* subDeep(const FIELD& m, const FIELD& n);
 709   static FIELD* mul(const FIELD& m, const FIELD& n);
 710   static FIELD* mulDeep(const FIELD& m, const FIELD& n);
 711   static FIELD* div(const FIELD& m, const FIELD& n);
 712   static FIELD* divDeep(const FIELD& m, const FIELD& n);
 713   double normMax() const throw (MEDEXCEPTION);
 714
 715   //------- TDG and BS addings
 716
 717   void getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION);
 718   vector<int> getHistogram(int &nbint) throw (MEDEXCEPTION);
 719   FIELD<double>* buildGradient() const throw (MEDEXCEPTION);
 720   FIELD<double>* buildNorm2Field() const throw (MEDEXCEPTION);
 721
 722   //-------------------
 723
 724   double norm2() const throw (MEDEXCEPTION);
 725   void   applyLin(T a, T b);
 726   template <T T_function(T)> void applyFunc();
 727   void applyPow(T scalar);
 728   static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
 729   double normL2(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 730   double normL2(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 731   double normL1(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 732   double normL1(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 733   FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
 734
 735   friend class MED_FIELD_RDONLY_DRIVER21<T>;
 736   friend class MED_FIELD_WRONLY_DRIVER21<T>;
 737   friend class MED_FIELD_RDONLY_DRIVER22<T>;
 738   friend class MED_FIELD_WRONLY_DRIVER22<T>;
 739   friend class VTK_FIELD_DRIVER<T>;
 740   //friend class MED_FIELD_RDWR_DRIVER  <T>;
 741
 742   void init ();
 743   void rmDriver(int index=0);
 744   int  addDriver(driverTypes driverType,
 745                  const string & fileName="Default File Name.med",
 746                  const string & driverFieldName="Default Field Name",
 747                  MED_EN::med_mode_acces access=MED_EN::MED_REMP) ;
 748
 749   int  addDriver(GENDRIVER & driver);
 750
 751   void allocValue(const int NumberOfComponents);
 752   void allocValue(const int NumberOfComponents, const int LengthValue);
 753
 754   void deallocValue();
 755
 756   inline void read(int index=0);
 757   inline void read(const GENDRIVER & genDriver);
 758   inline void write(int index=0, const string & driverName = "");
 759   inline void write(const GENDRIVER &);
 760
 761   inline void writeAppend(int index=0, const string & driverName = "");
 762   inline void writeAppend(const GENDRIVER &);
 763
 764   inline MEDMEM_Array_  * getArray()        const throw (MEDEXCEPTION);
 765   inline ArrayGauss     * getArrayGauss()   const throw (MEDEXCEPTION);
 766   inline ArrayNoGauss   * getArrayNoGauss() const throw (MEDEXCEPTION);
 767   inline bool            getGaussPresence() const throw (MEDEXCEPTION);
 768
 769   inline int          getValueLength() const throw (MEDEXCEPTION);
 770   inline const T*     getValue()       const throw (MEDEXCEPTION);
 771   inline const T*     getRow(int i)    const throw (MEDEXCEPTION);
 772   inline const T*     getColumn(int j) const throw (MEDEXCEPTION);
 773   inline T            getValueIJ(int i,int j) const throw (MEDEXCEPTION);
 774   inline T            getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION);
 775
 776   inline int          getValueByTypeLength(int t)                const throw (MEDEXCEPTION);
 777   inline const T*     getValueByType(int t)                      const throw (MEDEXCEPTION);
 778   inline T            getValueIJByType(int i,int j,int t)        const throw (MEDEXCEPTION);
 779   inline T            getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION);
 780
 781   bool                getValueOnElement(int eltIdInSup,T* retValues) const throw (MEDEXCEPTION);
 782
 783   const int   getNumberOfGeometricTypes() const throw (MEDEXCEPTION);
 784   const GAUSS_LOCALIZATION<INTERLACING_TAG> & getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 785   const GAUSS_LOCALIZATION<INTERLACING_TAG> * getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 786   const GAUSS_LOCALIZATION_* getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 787   void setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG> & gaussloc);
 788   void setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc);
 789   const int * getNumberOfGaussPoints() const throw (MEDEXCEPTION);
 790   const int   getNumberOfGaussPoints( MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 791   const int   getNbGaussI(int i)          const throw (MEDEXCEPTION);
 792   const int * getNumberOfElements()       const throw (MEDEXCEPTION);
 793   const MED_EN::medGeometryElement  * getGeometricTypes()  const throw (MEDEXCEPTION);
 794   bool        isOnAllElements()           const throw (MEDEXCEPTION);
 795
 796   inline void setArray(MEDMEM_Array_ *value) throw (MEDEXCEPTION);
 797   inline void setValue( T* value) throw (MEDEXCEPTION);
 798   inline void setRow( int i, T* value) throw (MEDEXCEPTION);
 799   inline void setColumn( int i, T* value) throw (MEDEXCEPTION);
 800   inline void setValueIJ(int i, int j, T value) throw (MEDEXCEPTION);
 801   inline void setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION);
 802   inline void setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION);
 803   inline void setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION);
 804
 805   /*!
 806     This fonction feeds the FIELD<double> private attributs _value with the
 807     volume of each cells belonging to the argument Support. The field has to be
 808     initialised via the constructor FIELD<double>(const SUPPORT * , const int )
 809     with Support as SUPPORT argument, 1 has the number of components, and Support
 810     has be a SUPPORT on 3D cells. This initialisation could be done by the empty
 811     constructor followed by a setSupport and setNumberOfComponents call.
 812    */
 813   void getVolume() const throw (MEDEXCEPTION) ;
 814   /*!
 815     This fonction feeds the FIELD<double> private attributs _value with the
 816     area of each cells (or faces) belonging to the attribut _support. The field
 817     has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
 818     const int ) with 1 has the number of components, and _support has be a
 819     SUPPORT on 2D cells or 3D faces. This initialisation could be done by the
 820     empty constructor followed by a setSupport and setNumberOfComponents call.
 821    */
 822   void getArea() const throw (MEDEXCEPTION) ;
 823   /*!
 824     This fonction feeds the FIELD<double> private attributs _value with the
 825     length of each segments belonging to the attribut _support. The field has
 826     to be initialised via the constructor FIELD<double>(const SUPPORT * ,
 827     const int ) with 1 has the number of components, and _support has be a
 828     SUPPORT on 3D edges or 2D faces. This initialisation could be done by the
 829     empty constructor followed by a setSupport and setNumberOfComponents call.
 830    */
 831   void getLength() const throw (MEDEXCEPTION) ;
 832   /*!
 833     This fonction feeds the FIELD<double> private attributs _value with the
 834     normal vector of each faces belonging to the attribut _support. The field
 835     has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
 836     const int ) with the space dimension has the number of components, and
 837     _support has be a SUPPORT on 3D or 2D faces. This initialisation could be done
 838     by the empty constructor followed by a setSupport and setNumberOfComponents
 839     call.
 840    */
 841   void getNormal() const throw (MEDEXCEPTION) ;
 842   /*!
 843     This fonction feeds the FIELD<double> private attributs _value with the
 844     barycenter of each faces or cells or edges belonging to the attribut _support.
 845     The field has to be initialised via the constructor
 846     FIELD<double>(const SUPPORT * ,const int ) with the space dimension has the
 847     number of components, and _support has be a SUPPORT on 3D cells or 2D faces.
 848     This initialisation could be done by the empty constructor followed by a
 849     setSupport and setNumberOfComponents call.
 850    */
 851   void getBarycenter() const throw (MEDEXCEPTION) ;
 852
 853   typedef void (*myFuncType)(const double *,T*);
 854   void fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION);
 855   typedef void (*myFuncType2)(const T *,T*);
 856   FIELD<T,INTERLACING_TAG> *execFunc(int nbOfComponents, myFuncType2 f) throw (MEDEXCEPTION);
 857 };
 858 }
 859
 860 #include "MEDMEM_DriverFactory.hxx"
 861
 862 namespace MEDMEM {
 863
 864 template <class T,class INTERLACING_TAG> T FIELD<T, INTERLACING_TAG>::_scalarForPow=1;
 865
 866 // --------------------
 867 // Implemented Methods
 868 // --------------------
 869
 870 /*!
 871   Constructor with no parameter, most of the attribut members are set to NULL.
 872 */
 873 template <class T, class INTERLACING_TAG>
 874 FIELD<T, INTERLACING_TAG>::FIELD():FIELD_()
 875 {
 876   MESSAGE("Constructeur FIELD sans parametre");
 877
 878   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
 879   ASSERT(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE);
 880   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
 881
 882   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
 883   ASSERT(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE);
 884   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
 885
 886   _value = ( ArrayNoGauss * ) NULL;
 887 }
 888
 889 /*!
 890   Constructor with parameters such that all attrribut are set but the _value
 891   attribut is allocated but not set.
 892 */
 893 template <class T, class INTERLACING_TAG>
 894 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
 895                                     const int NumberOfComponents) throw (MEDEXCEPTION) :
 896   FIELD_(Support, NumberOfComponents),_value(NULL)
 897 {
 898   BEGIN_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)");
 899   SCRUTE(this);
 900
 901   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
 902   ASSERT(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
 903   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
 904
 905   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
 906   ASSERT(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
 907   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
 908
 909   try {
 910     // becarefull about the numbre of gauss point
 911     _numberOfValues = Support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 912   }
 913 #if defined(_DEBUG_) || defined(_DEBUG)
 914   catch (MEDEXCEPTION &ex) {
 915 #else
 916   catch (MEDEXCEPTION ) {
 917 #endif
 918     MESSAGE("No value defined ! ("<<ex.what()<<")");
 919   }
 920   MESSAGE("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
 921   if (0<_numberOfValues) {
 922     _value = new ArrayNoGauss (_numberOfComponents,_numberOfValues);
 923     _isRead = true ;
 924   }
 925
 926   END_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)");
 927 }
 928
 929 /*!
 930   \if developper
 931   \endif
 932 */
 933 template <class T, class INTERLACING_TAG> void FIELD<T, INTERLACING_TAG>::init ()
 934 {
 935 }
 936
 937 /*!
 938   Copy constructor.
 939 */
 940 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::FIELD(const FIELD & m):
 941   FIELD_((FIELD_) m)
 942 {
 943   MESSAGE("Constructeur FIELD de recopie");
 944
 945   // RECOPIE PROFONDE <> de l'operateur= Rmq from EF
 946   if (m._value != NULL)
 947     {
 948       if ( m.getGaussPresence() )
 949         _value = new ArrayGauss( *(static_cast< ArrayGauss * > (m._value) ) ,false);
 950       else
 951         _value = new ArrayNoGauss( *(static_cast< ArrayNoGauss * > (m._value)) ,false);
 952     }
 953   else
 954     _value = (ArrayNoGauss *) NULL;
 955   locMap::const_iterator it;
 956   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
 957     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
 958       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
 959                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
 960                                               );
 961
 962   _valueType       = m._valueType;
 963   _interlacingType = m._interlacingType;
 964   //drivers = m._drivers;
 965 }
 966
 967 /*!
 968   \if developper
 969   Not implemented.
 970   \endif
 971 */
 972 template <class T, class INTERLACING_TAG>
 973 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(const FIELD &m)
 974 {
 975   MESSAGE("Appel de FIELD<T>::operator=") ;
 976   if ( this == &m) return *this;
 977
 978   // copy values array
 979   FIELD_::operator=(m);  // Driver are ignored & ?copie su pointeur de Support?
 980
 981   _value           = m._value; //PROBLEME RECOPIE DES POINTEURS PAS COHERENT AVEC LE
 982                                //CONSTRUCTEUR PAR RECOPIE
 983                                //CF :Commentaire dans MEDMEM_Array
 984   locMap::const_iterator it;
 985   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
 986     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
 987       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
 988                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
 989                                               );
 990
 991   _valueType       = m._valueType;
 992   _interlacingType = m._interlacingType;
 993
 994   return *this;
 995 }
 996
 997 /*!
 998      Overload addition operator.
 999      This operation is authorized only for compatible fields that have the same support.
1000      The compatibility checking includes equality tests of the folowing data members:\n
1001      - _support
1002      - _numberOfComponents
1003      - _numberOfValues
1004      - _componentsTypes
1005      - _MEDComponentsUnits.
1006
1007      The data members of the returned field are initialized, based on the first field, except for the name,
1008      which is the combination of the two field's names and the operator.
1009      Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> \n
1010      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1011      When using python, this operator calls the copy constructor in any case.
1012      The user has to be aware that when using operator + in associatives expressions like
1013      <tt> a = b + c + d +e; </tt> \n
1014      no optimisation is performed : the evaluation of last expression requires the construction of
1015      3 temporary fields.
1016 */
1017 template <class T, class INTERLACING_TAG>
1018 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator+(const FIELD & m) const
1019 {
1020     BEGIN_OF("FIELD<T>::operator+(const FIELD & m)");
1021     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1022
1023     // Creation of the result - memory is allocated by FIELD constructor
1024     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1025     //result._operation(*this,m,mode,"+"); // perform Atribute's initialization & addition
1026     result._operationInitialize(*this,m,"+"); // perform Atribute's initialization
1027     result._add_in_place(*this,m); // perform addition
1028
1029     END_OF("FIELD<T>::operator+(const FIELD & m)");
1030     return result;
1031 }
1032
1033 /*!  Overloaded Operator +=
1034  *   Operations are directly performed in the first field's array.
1035  *   This operation is authorized only for compatible fields that have the same support.
1036  */
1037 template <class T, class INTERLACING_TAG>
1038 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator+=(const FIELD & m)
1039 {
1040     BEGIN_OF("FIELD<T>::operator+=(const FIELD & m)");
1041     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1042
1043     const T* value1=m.getValue(); // get pointers to the values we are adding
1044
1045     // get a non const pointer to the inside array of values and perform operation
1046     T * value=const_cast<T *> (getValue());
1047     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1048     const T* endV=value+size; // pointer to the end of value
1049     for(;value!=endV; value1++,value++)
1050         *value += *value1;
1051     END_OF("FIELD<T>::operator+=(const FIELD & m)");
1052     return *this;
1053 }
1054
1055
1056 /*! Addition of fields. Static member function.
1057  *  The function return a pointer to a new created field that holds the addition.
1058  *  Data members are checked for compatibility and initialized.
1059  *  The user is in charge of memory deallocation.
1060  */
1061 template <class T, class INTERLACING_TAG>
1062 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::add(const FIELD& m, const FIELD& n)
1063 {
1064     BEGIN_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
1065     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1066
1067     // Creation of a new field
1068     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1069                                                                       m.getNumberOfComponents());
1070     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1071     result->_add_in_place(m,n); // perform addition
1072
1073     END_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
1074     return result;
1075 }
1076
1077 /*! Same as add method except that field check is deeper.
1078  */
1079 template <class T, class INTERLACING_TAG>
1080 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::addDeep(const FIELD& m, const FIELD& n)
1081 {
1082     BEGIN_OF("FIELD<T>::addDeep(const FIELD & m, const FIELD& n)");
1083     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1084
1085     // Creation of a new field
1086     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1087                                                                       m.getNumberOfComponents());
1088     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1089     result->_add_in_place(m,n); // perform addition
1090
1091     END_OF("FIELD<T>::addDeep(const FIELD & m, const FIELD& n)");
1092     return result;
1093 }
1094
1095 /*!
1096      Overload substraction operator.
1097      This operation is authorized only for compatible fields that have the same support.
1098      The compatibility checking includes equality tests of the folowing data members:\n
1099      - _support
1100      - _numberOfComponents
1101      - _numberOfValues
1102      - _componentsTypes
1103      - _MEDComponentsUnits.
1104
1105      The data members of the returned field are initialized, based on the first field, except for the name,
1106      which is the combination of the two field's names and the operator.
1107      Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> \n
1108      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1109      When using python, this operator calls the copy constructor in any case.
1110      The user has to be aware that when using operator - in associatives expressions like
1111      <tt> a = b - c - d -e; </tt> \n
1112      no optimisation is performed : the evaluation of last expression requires the construction of
1113      3 temporary fields.
1114 */
1115 template <class T, class INTERLACING_TAG>
1116 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator-(const FIELD & m) const
1117 {
1118     BEGIN_OF("FIELD<T>::operator-(const FIELD & m)");
1119     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1120
1121     // Creation of the result - memory is allocated by FIELD constructor
1122     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1123     //result._operation(*this,m,mode,"-"); // perform Atribute's initialization & substraction
1124     result._operationInitialize(*this,m,"-"); // perform Atribute's initialization
1125     result._sub_in_place(*this,m); // perform substracion
1126
1127     END_OF("FIELD<T>::operator-(const FIELD & m)");
1128     return result;
1129 }
1130
1131 template <class T, class INTERLACING_TAG>
1132 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator-() const
1133 {
1134     BEGIN_OF("FIELD<T>::operator-()");
1135
1136     // Creation of the result - memory is allocated by FIELD constructor
1137     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1138     // Atribute's initialization
1139     result.setName("- "+getName());
1140     result.setComponentsNames(getComponentsNames());
1141     // not yet implemented    setComponentType(getComponentType());
1142     result.setComponentsDescriptions(getComponentsDescriptions());
1143     result.setMEDComponentsUnits(getMEDComponentsUnits());
1144     result.setComponentsUnits(getComponentsUnits());
1145     result.setIterationNumber(getIterationNumber());
1146     result.setTime(getTime());
1147     result.setOrderNumber(getOrderNumber());
1148
1149     const T* value1=getValue();
1150     // get a non const pointer to the inside array of values and perform operation
1151     T * value=const_cast<T *> (result.getValue());
1152     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1153     const T* endV=value+size; // pointer to the end of value
1154
1155     for(;value!=endV; value1++,value++)
1156         *value = -(*value1);
1157     END_OF("FIELD<T>::operator-=(const FIELD & m)");
1158     return result;
1159 }
1160
1161 /*!  Overloaded Operator -=
1162  *   Operations are directly performed in the first field's array.
1163  *   This operation is authorized only for compatible fields that have the same support.
1164  */
1165 template <class T, class INTERLACING_TAG>
1166 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator-=(const FIELD & m)
1167 {
1168     BEGIN_OF("FIELD<T>::operator-=(const FIELD & m)");
1169     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1170
1171     const T* value1=m.getValue();
1172
1173     // get a non const pointer to the inside array of values and perform operation
1174     T * value=const_cast<T *> (getValue());
1175     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1176     const T* endV=value+size; // pointer to the end of value
1177
1178     for(;value!=endV; value1++,value++)
1179         *value -= *value1;
1180
1181     END_OF("FIELD<T>::operator-=(const FIELD & m)");
1182     return *this;
1183 }
1184
1185
1186 /*! Substraction of fields. Static member function.
1187  *  The function return a pointer to a new created field that holds the substraction.
1188  *  Data members are checked for compatibility and initialized.
1189  *  The user is in charge of memory deallocation.
1190  */
1191 template <class T, class INTERLACING_TAG>
1192 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::sub(const FIELD& m, const FIELD& n)
1193 {
1194     BEGIN_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
1195     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1196
1197     // Creation of a new field
1198     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1199                                                                       m.getNumberOfComponents());
1200     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1201     result->_sub_in_place(m,n); // perform substraction
1202
1203     END_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
1204     return result;
1205 }
1206
1207 /*! Same as sub method except that field check is deeper.
1208  */
1209 template <class T, class INTERLACING_TAG>
1210 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::subDeep(const FIELD& m, const FIELD& n)
1211 {
1212     BEGIN_OF("FIELD<T>::subDeep(const FIELD & m, const FIELD& n)");
1213     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1214
1215     // Creation of a new field
1216     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1217                                                                       m.getNumberOfComponents());
1218     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1219     result->_sub_in_place(m,n); // perform substraction
1220
1221     END_OF("FIELD<T>::subDeep(const FIELD & m, const FIELD& n)");
1222     return result;
1223 }
1224
1225 /*!
1226      Overload multiplication operator.
1227      This operation is authorized only for compatible fields that have the same support.
1228      The compatibility checking includes equality tests of the folowing data members:\n
1229      - _support
1230      - _numberOfComponents
1231      - _numberOfValues
1232      - _componentsTypes
1233      - _MEDComponentsUnits.
1234
1235      The data members of the returned field are initialized, based on the first field, except for the name,
1236      which is the combination of the two field's names and the operator.
1237      Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> \n
1238      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1239      When using python, this operator calls the copy constructor in any case.
1240      The user has to be aware that when using operator * in associatives expressions like
1241      <tt> a = b * c * d *e; </tt> \n
1242      no optimisation is performed : the evaluation of last expression requires the construction of
1243      3 temporary fields.
1244 */
1245 template <class T, class INTERLACING_TAG>
1246 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator*(const FIELD & m) const
1247 {
1248     BEGIN_OF("FIELD<T>::operator*(const FIELD & m)");
1249     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1250
1251     // Creation of the result - memory is allocated by FIELD constructor
1252     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1253     //result._operation(*this,m,mode,"*"); // perform Atribute's initialization & multiplication
1254     result._operationInitialize(*this,m,"*"); // perform Atribute's initialization
1255     result._mul_in_place(*this,m); // perform multiplication
1256
1257     END_OF("FIELD<T>::operator*(const FIELD & m)");
1258     return result;
1259 }
1260
1261 /*!  Overloaded Operator *=
1262  *   Operations are directly performed in the first field's array.
1263  *   This operation is authorized only for compatible fields that have the same support.
1264  */
1265 template <class T, class INTERLACING_TAG>
1266 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator*=(const FIELD & m)
1267 {
1268     BEGIN_OF("FIELD<T>::operator*=(const FIELD & m)");
1269     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1270
1271     const T* value1=m.getValue();
1272
1273     // get a non const pointer to the inside array of values and perform operation
1274     T * value=const_cast<T *> (getValue());
1275     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1276     const T* endV=value+size; // pointer to the end of value
1277
1278     for(;value!=endV; value1++,value++)
1279         *value *= *value1;
1280
1281     END_OF("FIELD<T>::operator*=(const FIELD & m)");
1282     return *this;
1283 }
1284
1285
1286 /*! Multiplication of fields. Static member function.
1287  *  The function return a pointer to a new created field that holds the multiplication.
1288  *  Data members are checked for compatibility and initialized.
1289  *  The user is in charge of memory deallocation.
1290  */
1291 template <class T, class INTERLACING_TAG>
1292 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mul(const FIELD& m, const FIELD& n)
1293 {
1294     BEGIN_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1295     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1296
1297     // Creation of a new field
1298     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1299                                                                       m.getNumberOfComponents());
1300     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1301     result->_mul_in_place(m,n); // perform multiplication
1302
1303     END_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1304     return result;
1305 }
1306
1307 /*! Same as mul method except that field check is deeper.
1308  */
1309 template <class T, class INTERLACING_TAG>
1310 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mulDeep(const FIELD& m, const FIELD& n)
1311 {
1312     BEGIN_OF("FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)");
1313     FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1314
1315     // Creation of a new field
1316     FIELD<T, INTERLACING_TAG>* result = new FIELD<T,INTERLACING_TAG>(m.getSupport(),
1317                                                                      m.getNumberOfComponents());
1318     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1319     result->_mul_in_place(m,n); // perform multiplication
1320
1321     END_OF("FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)");
1322     return result;
1323 }
1324
1325 /*!
1326      Overload division operator.
1327      This operation is authorized only for compatible fields that have the same support.
1328      The compatibility checking includes equality tests of the folowing data members:\n
1329      - _support
1330      - _numberOfComponents
1331      - _numberOfValues
1332      - _componentsTypes
1333      - _MEDComponentsUnits.
1334
1335      The data members of the returned field are initialized, based on the first field, except for the name,
1336      which is the combination of the two field's names and the operator.
1337      Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> \n
1338      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1339      When using python, this operator calls the copy constructor in any case.
1340      The user has to be aware that when using operator / in associatives expressions like
1341      <tt> a = b / c / d /e; </tt> \n
1342      no optimisation is performed : the evaluation of last expression requires the construction of
1343      3 temporary fields.
1344 */
1345 template <class T, class INTERLACING_TAG>
1346 const FIELD<T, INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::operator/(const FIELD & m) const
1347 {
1348     BEGIN_OF("FIELD<T>::operator/(const FIELD & m)");
1349     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1350
1351     // Creation of the result - memory is allocated by FIELD constructor
1352     FIELD<T, INTERLACING_TAG> result(this->getSupport(),this->getNumberOfComponents());
1353     //result._operation(*this,m,mode,"/"); // perform Atribute's initialization & division
1354     result._operationInitialize(*this,m,"/"); // perform Atribute's initialization
1355     result._div_in_place(*this,m); // perform division
1356
1357     END_OF("FIELD<T>::operator/(const FIELD & m)");
1358     return result;
1359 }
1360
1361
1362 /*!  Overloaded Operator /=
1363  *   Operations are directly performed in the first field's array.
1364  *   This operation is authorized only for compatible fields that have the same support.
1365  */
1366 template <class T, class INTERLACING_TAG>
1367 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator/=(const FIELD & m)
1368 {
1369     BEGIN_OF("FIELD<T>::operator/=(const FIELD & m)");
1370     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1371
1372     const T* value1=m.getValue(); // get pointers to the values we are adding
1373
1374     // get a non const pointer to the inside array of values and perform operation
1375     T * value=const_cast<T *> (getValue());
1376     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1377     const T* endV=value+size; // pointer to the end of value
1378
1379     for(;value!=endV; value1++,value++)
1380         *value /= *value1;
1381
1382     END_OF("FIELD<T>::operator/=(const FIELD & m)");
1383     return *this;
1384 }
1385
1386
1387 /*! Division of fields. Static member function.
1388  *  The function return a pointer to a new created field that holds the division.
1389  *  Data members are checked for compatibility and initialized.
1390  *  The user is in charge of memory deallocation.
1391  */
1392 template <class T, class INTERLACING_TAG>
1393 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::div(const FIELD& m, const FIELD& n)
1394 {
1395     BEGIN_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1396     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1397
1398     // Creation of a new field
1399     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1400                                                                       m.getNumberOfComponents());
1401     result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1402     result->_div_in_place(m,n); // perform division
1403
1404     END_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1405     return result;
1406 }
1407
1408 /*! Same as div method except that field check is deeper.
1409  */
1410 template <class T,class INTERLACING_TAG>
1411 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::divDeep(const FIELD& m, const FIELD& n)
1412 {
1413   BEGIN_OF("FIELD<T>::divDeep(const FIELD & m, const FIELD& n)");
1414   FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1415
1416   // Creation of a new field
1417   FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1418                                                                     m.getNumberOfComponents());
1419   result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1420   result->_div_in_place(m,n); // perform division
1421
1422   END_OF("FIELD<T>::divDeep(const FIELD & m, const FIELD& n)");
1423   return result;
1424 }
1425
1426 /*!
1427   \if developper
1428   This internal method initialize the members of a new field created to hold the result of the operation Op .
1429   Initialization is based on the first field, except for the name, which is the combination of the two field's names
1430   and the operator.
1431   \endif
1432 */
1433 template <class T, class INTERLACING_TAG>
1434 void FIELD<T, INTERLACING_TAG>::_operationInitialize(const FIELD& m,const FIELD& n, char* Op)
1435 {
1436     MESSAGE("Appel methode interne " << Op);
1437
1438     // Atribute's initialization (copy of the first field's attributes)
1439     // Other data members (_support, _numberOfValues) are initialized in the field's constr.
1440     setName(m.getName()+" "+Op+" "+n.getName());
1441     setComponentsNames(m.getComponentsNames());
1442     // not yet implemented    setComponentType(m.getComponentType());
1443     setComponentsDescriptions(m.getComponentsDescriptions());
1444     setMEDComponentsUnits(m.getMEDComponentsUnits());
1445
1446     // The following data member may differ from field m to n.
1447     // The initialization is done based on the first field.
1448
1449     if(m.getComponentsUnits() != NULL)
1450       setComponentsUnits(m.getComponentsUnits());
1451     else
1452       _componentsUnits = (UNIT *) NULL;
1453
1454     setIterationNumber(m.getIterationNumber());
1455     setTime(m.getTime());
1456     setOrderNumber(m.getOrderNumber());
1457 }
1458
1459
1460 /*!
1461   \if developper
1462   Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
1463   This method is applied to a just created field with medModeSwitch mode.
1464   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1465   it doesn't exist!
1466   \endif
1467 */
1468 template <class T, class INTERLACING_TAG>
1469 void FIELD<T, INTERLACING_TAG>::_add_in_place(const FIELD& m,const FIELD& n)
1470 {
1471     // get pointers to the values we are adding
1472     const T* value1=m.getValue();
1473     const T* value2=n.getValue();
1474     // get a non const pointer to the inside array of values and perform operation
1475     T * value=const_cast<T *> (getValue());
1476
1477     const int size=getNumberOfValues()*getNumberOfComponents();
1478     SCRUTE(size);
1479     const T* endV1=value1+size;
1480     for(;value1!=endV1; value1++,value2++,value++)
1481         *value=(*value1)+(*value2);
1482 }
1483
1484 /*!
1485   \if developper
1486   Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
1487   This method is applied to a just created field with medModeSwitch mode.
1488   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1489   it doesn't exist!
1490   \endif
1491 */
1492 template <class T, class INTERLACING_TAG>
1493 void FIELD<T, INTERLACING_TAG>::_sub_in_place(const FIELD& m,const FIELD& n)
1494 {
1495     // get pointers to the values we are adding
1496     const T* value1=m.getValue();
1497     const T* value2=n.getValue();
1498     // get a non const pointer to the inside array of values and perform operation
1499     T * value=const_cast<T *> (getValue());
1500
1501     const int size=getNumberOfValues()*getNumberOfComponents();
1502     SCRUTE(size);
1503     const T* endV1=value1+size;
1504     for(;value1!=endV1; value1++,value2++,value++)
1505         *value=(*value1)-(*value2);
1506 }
1507
1508 /*!
1509   \if developper
1510   Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
1511   This method is applied to a just created field with medModeSwitch mode.
1512   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1513   it doesn't exist!
1514   \endif
1515 */
1516 template <class T, class INTERLACING_TAG>
1517 void FIELD<T, INTERLACING_TAG>::_mul_in_place(const FIELD& m,const FIELD& n)
1518 {
1519     // get pointers to the values we are adding
1520     const T* value1=m.getValue();
1521     const T* value2=n.getValue();
1522     // get a non const pointer to the inside array of values and perform operation
1523     T * value=const_cast<T *> (getValue());
1524
1525     const int size=getNumberOfValues()*getNumberOfComponents();
1526     SCRUTE(size);
1527     const T* endV1=value1+size;
1528     for(;value1!=endV1; value1++,value2++,value++)
1529         *value=(*value1)*(*value2);
1530 }
1531
1532 /*!
1533   \if developper
1534   Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
1535   This method is applied to a just created field with medModeSwitch mode.
1536   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1537   it doesn't exist!
1538   \endif
1539 */
1540 template <class T, class INTERLACING_TAG>
1541 void FIELD<T, INTERLACING_TAG>::_div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION)
1542 {
1543     // get pointers to the values we are adding
1544     const T* value1=m.getValue();
1545     const T* value2=n.getValue();
1546     // get a non const pointer to the inside array of values and perform operation
1547     T * value=const_cast<T *> (getValue());
1548
1549     const int size=getNumberOfValues()*getNumberOfComponents();
1550     SCRUTE(size);
1551     const T* endV1=value1+size;
1552     for(;value1!=endV1; value1++,value2++,value++){
1553       if ( *value2 == 0 ) { // FAIRE PLUTOT UN TRY CATCH Rmq from EF
1554           string diagnosis;
1555           diagnosis="FIELD<T,INTERLACING_TAG>::_div_in_place(...) : Divide by zero !";
1556           throw MEDEXCEPTION(diagnosis.c_str());
1557         }
1558         *value=(*value1)/(*value2);
1559     }
1560 }
1561
1562 /*!  Return Max Norm
1563  */
1564 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::normMax() const throw (MEDEXCEPTION)
1565 {
1566     const T* value=getValue(); // get pointer to the values
1567     const int size=getNumberOfValues()*getNumberOfComponents();
1568     if (size <= 0) // Size of array has to be strictly positive
1569     {
1570         string diagnosis;
1571         diagnosis="FIELD<T,INTERLACIN_TAG>::normMax() : cannot compute the norm of "+getName()+
1572             " : it size is non positive!";
1573         throw MEDEXCEPTION(diagnosis.c_str());
1574     }
1575     const T* lastvalue=value+size; // get pointer just after last value
1576     const T* pMax=value; // pointer to the max value
1577     const T* pMin=value; // pointer to the min value
1578
1579     // get pointers to the max & min value of array
1580     while ( ++value != lastvalue )
1581     {
1582         if ( *pMin > *value )
1583             pMin=value;
1584         if ( *pMax < *value )
1585             pMax=value;
1586     }
1587
1588     T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
1589     T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
1590
1591     return Max>Min ? static_cast<double>(Max) : static_cast<double>(Min);
1592 }
1593
1594 /*!  Return Euclidien norm
1595  */
1596 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::norm2() const throw (MEDEXCEPTION)
1597 {
1598     const T* value=this->getValue(); // get const pointer to the values
1599     const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
1600     if (size <= 0) // Size of array has to be strictly positive
1601     {
1602         string diagnosis;
1603         diagnosis="FIELD<T,INTERLACIN_TAG>::norm2() : cannot compute the norm of "+getName()+
1604             " : it size is non positive!";
1605         throw MEDEXCEPTION(diagnosis.c_str());
1606     }
1607     const T* lastvalue=value+size; // point just after last value
1608
1609     T result((T)0); // init
1610     for( ; value!=lastvalue ; ++value)
1611         result += (*value) * (*value);
1612
1613     return std::sqrt(static_cast<double> (result));
1614 }
1615
1616
1617 //------------- TDG and BS addings
1618
1619 /*!  Return Extremums of field
1620  */
1621  template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION)
1622 {
1623   const T* value=getValue(); // get pointer to the values
1624   const int size=getNumberOfValues()*getNumberOfComponents();
1625   const T* lastvalue=value+size; // point just after last value
1626
1627   if (size <= 0){ // Size of array has to be strictly positive
1628
1629     string diagnosis;
1630     diagnosis="FIELD<T,INTERLACIN_TAG>::getMinMax() : cannot compute the extremums of "+getName()+
1631       " : its size is non positive!";
1632     throw MEDEXCEPTION(diagnosis.c_str());
1633   }
1634
1635   if (!_isMinMax){
1636     vmax=MinMax<T>::getMin(); // init a max value
1637     vmin=MinMax<T>::getMax(); // init a min value
1638
1639     for( ; value!=lastvalue ; ++value){
1640       if ( vmin > *value )
1641         vmin=*value;
1642       if ( vmax < *value )
1643         vmax=*value;
1644     }
1645     _isMinMax=true;
1646     _vmin=vmin;
1647     _vmax=vmax;
1648   }
1649   else{
1650     vmin = _vmin;
1651     vmax = _vmax;
1652   }
1653
1654 }
1655
1656 /*!  Return Histogram of field
1657  */
1658  template <class T, class INTERLACIN_TAG> vector<int> FIELD<T, INTERLACIN_TAG>::getHistogram(int &nbint) throw (MEDEXCEPTION)
1659 {
1660   const T* value=getValue(); // get pointer to the values
1661   const int size=getNumberOfValues()*getNumberOfComponents();
1662   const T* lastvalue=value+size; // point just after last value
1663
1664   if (size <= 0){ // Size of array has to be strictly positive
1665
1666     string diagnosis;
1667     diagnosis="FIELD<T,INTERLACIN_TAG>::getHistogram() : cannot compute the histogram of "+getName()+
1668       " : it size is non positive!";
1669     throw MEDEXCEPTION(diagnosis.c_str());
1670   }
1671   //    return static_cast<ArrayGauss *>(_value)->getIJ(valIndex,j) ;
1672
1673   vector<int> Histogram(nbint) ;
1674   T vmin,vmax;
1675   int j;
1676
1677   for( j=0 ; j!=nbint ; j++) Histogram[j]=0 ;
1678
1679   getMinMax(vmin,vmax);
1680   for( ; value!=lastvalue ; ++value){
1681     if(*value==vmax) j = nbint-1;
1682     else j = (int)(((double)nbint * (*value-vmin))/(vmax-vmin));
1683     Histogram[j]+=1 ;
1684   }
1685
1686   return Histogram ;
1687
1688 }
1689
1690 /*!  Return vectorial gradient field
1691  */
1692 template <class T, class INTERLACIN_TAG>
1693 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildGradient() const throw (MEDEXCEPTION)
1694 {
1695   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildGradient() : ";
1696   BEGIN_OF(LOC);
1697
1698   // space dimension of input mesh
1699   int spaceDim = getSupport()->getMesh()->getSpaceDimension();
1700   double *x = new double[spaceDim];
1701
1702   FIELD<double, FullInterlace>* Gradient =
1703     new FIELD<double, FullInterlace>(getSupport(),spaceDim);
1704
1705   string name("gradient of ");
1706   name += getName();
1707   Gradient->setName(name);
1708   string descr("gradient of ");
1709   descr += getDescription();
1710   Gradient->setDescription(descr);
1711
1712   if( _numberOfComponents > 1 ){
1713     delete Gradient;
1714     delete [] x;
1715     throw MEDEXCEPTION("gradient calculation only on scalar field");
1716   }
1717
1718   for(int i=1;i<=spaceDim;i++){
1719     string nameC("gradient of ");
1720     nameC += getName();
1721     Gradient->setComponentName(i,nameC);
1722     Gradient->setComponentDescription(i,"gradient");
1723     string MEDComponentUnit = getMEDComponentUnit(1)+getSupport()->getMesh()->getCoordinatesUnits()[i-1];
1724     Gradient->setMEDComponentUnit(i,MEDComponentUnit);
1725   }
1726
1727   Gradient->setIterationNumber(getIterationNumber());
1728   Gradient->setOrderNumber(getOrderNumber());
1729   Gradient->setTime(getTime());
1730
1731   // typ of entity on what is field
1732   MED_EN::medEntityMesh typ = getSupport()->getEntity();
1733
1734   const int *C;
1735   const int *iC;
1736   const int *revC;
1737   const int *indC;
1738   const double *coord;
1739   int NumberOf;
1740
1741   switch (typ) {
1742   case MED_CELL:
1743   case MED_FACE:
1744   case MED_EDGE:
1745     {
1746       // read connectivity array to have the list of nodes contained by an element
1747       C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,typ,MED_ALL_ELEMENTS);
1748       iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,typ);
1749       // calculate reverse connectivity to have the list of elements which contains node i
1750       revC = getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
1751       indC = getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
1752       // coordinates of each node
1753       coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
1754       // number of elements
1755       NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
1756       // barycenter field of elements
1757       FIELD<double, FullInterlace>* barycenter = getSupport()->getMesh()->getBarycenter(getSupport());
1758
1759       // calculate gradient vector for each element i
1760       for (int i = 1; i < NumberOf + 1; i++) {
1761
1762         // listElements contains elements which contains a node of element i
1763         set <int> listElements;
1764         set <int>::iterator elemIt;
1765         listElements.clear();
1766
1767         // for each node j of element i
1768         for (int ij = iC[i-1]; ij < iC[i]; ij++) {
1769           int j = C[ij-1];
1770           for (int k = indC[j-1]; k < indC[j]; k++) {
1771             // c element contains node j
1772             int c = revC[k-1];
1773             // we put the elements in set
1774             if (c != i)
1775               listElements.insert(c);
1776           }
1777         }
1778         // coordinates of barycentre of element i in space of dimension spaceDim
1779         for (int j = 0; j < spaceDim; j++)
1780           x[j] = barycenter->getValueIJ(i,j+1);
1781
1782         for (int j = 0; j < spaceDim; j++) {
1783           // value of field of element i
1784           double val = getValueIJ(i,1);
1785           double grad = 0.;
1786           // calculate gradient for each neighbor element
1787           for (elemIt = listElements.begin(); elemIt != listElements.end(); elemIt++) {
1788             int elem = *elemIt;
1789             double d2 = 0.;
1790             for (int l = 0; l < spaceDim; l++) {
1791               // coordinate of barycenter of element elem
1792               double xx = barycenter->getValueIJ(elem, l+1);
1793               d2 += (x[l]-xx) * (x[l]-xx);
1794             }
1795             grad += (barycenter->getValueIJ(elem,j+1)-x[j])*(getValueIJ(elem,1)-val)/sqrt(d2);
1796           }
1797           if (listElements.size() != 0) grad /= listElements.size();
1798           Gradient->setValueIJ(i,j+1,grad);
1799         }
1800       }
1801       delete barycenter;
1802     }
1803     break;
1804   case MED_NODE:
1805     // read connectivity array to have the list of nodes contained by an element
1806     C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
1807     iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,MED_CELL);
1808     // calculate reverse connectivity to have the list of elements which contains node i
1809     revC=getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
1810     indC=getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
1811     // coordinates of each node
1812     coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
1813
1814     // calculate gradient for each node
1815     NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
1816     for (int i=1; i<NumberOf+1; i++){
1817       // listNodes contains nodes neigbor of node i
1818       set <int> listNodes;
1819       set <int>::iterator nodeIt ;
1820       listNodes.clear();
1821       for(int j=indC[i-1];j<indC[i];j++){
1822         // c element contains node i
1823         int c=revC[j-1];
1824         // we put the nodes of c element in set
1825         for(int k=iC[c-1];k<iC[c];k++)
1826           if(C[k-1] != i)
1827             listNodes.insert(C[k-1]);
1828       }
1829       // coordinates of node i in space of dimension spaceDim
1830       for(int j=0;j<spaceDim;j++)
1831         x[j] = coord[(i-1)*spaceDim+j];
1832
1833       for(int j=0;j<spaceDim;j++){
1834         // value of field
1835         double val = getValueIJ(i,1);
1836         double grad = 0.;
1837         // calculate gradient for each neighbor node
1838         for(nodeIt=listNodes.begin();nodeIt!=listNodes.end();nodeIt++){
1839           int node = *nodeIt;
1840           double d2 = 0.;
1841           for(int l=0;l<spaceDim;l++){
1842             double xx = coord[(node-1)*spaceDim+l];
1843             d2 += (x[l]-xx) * (x[l]-xx);
1844           }
1845           grad += (coord[(node-1)*spaceDim+j]-x[j])*(getValueIJ(node,1)-val)/sqrt(d2);
1846         }
1847         if(listNodes.size() != 0) grad /= listNodes.size();
1848         Gradient->setValueIJ(i,j+1,grad);
1849       }
1850     }
1851     break;
1852   case MED_ALL_ENTITIES:
1853     delete [] x;
1854     delete Gradient;
1855     throw MEDEXCEPTION("gradient calculation not yet implemented on all elements");
1856     break;
1857   }
1858
1859   delete [] x;
1860
1861   END_OF(LOC);
1862   return Gradient;
1863 }
1864
1865 /*!  Return scalar norm2 field
1866  */
1867 template <class T, class INTERLACIN_TAG>
1868 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildNorm2Field() const throw (MEDEXCEPTION)
1869 {
1870   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildNorm2Field() : ";
1871   BEGIN_OF(LOC);
1872
1873   FIELD<double, FullInterlace>* Norm2Field =
1874     new FIELD<double, FullInterlace>(getSupport(),1);
1875
1876   string name("norm2 of ");
1877   name += getName();
1878   Norm2Field->setName(name);
1879   string descr("norm2 of ");
1880   descr += getDescription();
1881   Norm2Field->setDescription(descr);
1882
1883   string nameC("norm2 of ");
1884   nameC += getName();
1885   Norm2Field->setComponentName(1,nameC);
1886   Norm2Field->setComponentDescription(1,"norm2");
1887   string MEDComponentUnit = getMEDComponentUnit(1);
1888   Norm2Field->setMEDComponentUnit(1,MEDComponentUnit);
1889
1890   Norm2Field->setIterationNumber(getIterationNumber());
1891   Norm2Field->setOrderNumber(getOrderNumber());
1892   Norm2Field->setTime(getTime());
1893
1894   // calculate nom2 for each element
1895   int NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
1896   for (int i=1; i<NumberOf+1; i++){
1897     double norm2 = 0.;
1898     for(int j=1;j<=getNumberOfComponents();j++)
1899       norm2 += getValueIJ(i,j)*getValueIJ(i,j);
1900     Norm2Field->setValueIJ(i,1,sqrt(norm2));
1901   }
1902
1903   END_OF(LOC);
1904   return Norm2Field;
1905
1906 }
1907
1908 /*!  Apply to each (scalar) field component the template parameter T_function,
1909  *   which is a pointer to function.
1910  *   Since the pointer is known at compile time, the function is inlined into the inner loop!
1911  *   calculation is done "in place".
1912  *   Use examples :
1913  *
1914  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
1915  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
1916  */
1917 template <class T, class INTERLACIN_TAG> template <T T_function(T)>
1918 void FIELD<T, INTERLACIN_TAG>::applyFunc()
1919 {
1920   // get a non const pointer to the inside array of values and perform operation
1921     T * value=const_cast<T *> (getValue());
1922     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1923
1924     if (size>0) // for a negative size, there is nothing to do
1925     {
1926         const T* lastvalue=value+size; // pointer to the end of value
1927         for(;value!=lastvalue; ++value) // apply linear transformation
1928             *value = T_function(*value);
1929     }
1930 }
1931
1932 template <class T, class INTERLACIN_TAG> T FIELD<T, INTERLACIN_TAG>::pow(T x)
1933 {
1934   return (T)::pow(x,FIELD<T, INTERLACIN_TAG>::_scalarForPow);
1935 }
1936
1937 /*!  Apply to each (scalar) field component the math function pow.
1938  *   calculation is done "in place".
1939  *   Use examples :
1940  *
1941  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
1942  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
1943  */
1944 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyPow(T scalar)
1945 {
1946   FIELD<T, INTERLACIN_TAG>::_scalarForPow=scalar;
1947   applyFunc<FIELD<T, INTERLACIN_TAG>::pow>();
1948 }
1949
1950 /*!  Apply to each (scalar) field component the linear function x -> ax+b.
1951  *   calculation is done "in place".
1952  */
1953 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b)
1954 {
1955     // get a non const pointer to the inside array of values and perform operation in place
1956     T * value=const_cast<T *> (getValue());
1957     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1958
1959     if (size>0) // for a negative size, there is nothing to do
1960     {
1961         const T* lastvalue=value+size; // pointer to the end of value
1962         for(;value!=lastvalue; ++value) // apply linear transformation
1963             *value = a*(*value)+b;
1964     }
1965 }
1966
1967
1968 /*!
1969  *   Return a pointer to a new field that holds the scalar product. Static member function.
1970  *   This operation is authorized only for compatible fields that have the same support.
1971  *   The compatibility checking includes equality tests of the folowing data members:\n
1972  *   - _support
1973  *   - _numberOfComponents
1974  *   - _numberOfValues
1975  *   - _componentsTypes
1976  *   - _MEDComponentsUnits.
1977  *   Data members are initialized.
1978  *   The new field point to the same support and has one component.
1979  *   Each value of it is the scalar product of the two argument's fields.
1980  *   The user is in charge of memory deallocation.
1981  */
1982 template <class T, class INTERLACING_TAG>
1983 FIELD<T, INTERLACING_TAG>*
1984 FIELD<T, INTERLACING_TAG>::scalarProduct(const FIELD & m, const FIELD & n, bool deepCheck)
1985 {
1986     if(!deepCheck)
1987       FIELD_::_checkFieldCompatibility( m, n, false); // may throw exception
1988     else
1989       FIELD_::_deepCheckFieldCompatibility(m, n, false);
1990
1991     // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
1992     // result type imply INTERLACING_TAG=FullInterlace for m & n
1993
1994     const int numberOfElements=m.getNumberOfValues(); // strictly positive
1995     const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
1996
1997     // Creation & init of a the result field on the same support, with one component
1998     // You have to be careful about the interlacing mode, because in the computation step,
1999     // it seems to assume the the interlacing mode is the FullInterlacing
2000
2001     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),1);
2002     result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
2003     result->setIterationNumber(m.getIterationNumber());
2004     result->setTime(m.getTime());
2005     result->setOrderNumber(m.getOrderNumber());
2006
2007     const T* value1=m.getValue(); // get const pointer to the values
2008     const T* value2=n.getValue(); // get const pointer to the values
2009     // get a non const pointer to the inside array of values and perform operation
2010     T * value=const_cast<T *> (result->getValue());
2011
2012     const T* lastvalue=value+numberOfElements; // pointing just after last value of result
2013     for ( ; value!=lastvalue ; ++value ) // loop on all elements
2014     {
2015         *value=(T)0; // initialize value
2016         const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
2017         for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
2018             *value += (*value1) * (*value2);
2019     }
2020     return result;
2021 }
2022
2023 /*!  Return L2 Norm  of the field's component.
2024  *   Cannot be applied to a field with a support on nodes.
2025  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2026  */
2027 template <class T, class INTERLACING_TAG>
2028 double FIELD<T, INTERLACING_TAG>::normL2(int component,
2029                                          const FIELD<double, FullInterlace> * p_field_volume) const
2030 {
2031     _checkNormCompatibility(p_field_volume); // may throw exception
2032     if ( component<1 || component>getNumberOfComponents() )
2033         throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
2034
2035     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2036     if(!p_field_volume) // if the user don't supply the volume
2037         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
2038
2039     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2040     const double* vol=p_field_size->getValue();
2041     // Il n'est vraiment pas optimal de mixer des champs dans des modes d'entrelacement
2042     // different juste pour le calcul
2043
2044     const T * value     = NULL;
2045     ArrayNo * myArray   = NULL;
2046     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2047       value = getValue();
2048     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2049       myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2050       value   = myArray->getPtr();
2051     }
2052     else {
2053       myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2054       value   = myArray->getPtr();
2055     }
2056
2057     value = value + (component-1) * getNumberOfValues();
2058     const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
2059
2060     double integrale=0.0;
2061     double totVol=0.0;
2062     for (; value!=lastvalue ; ++value ,++vol)
2063     {
2064         integrale += static_cast<double>((*value) * (*value)) * (*vol);
2065         totVol+=*vol;
2066     }
2067
2068     if(!p_field_volume) // if the user didn't supply the volume
2069         delete p_field_size; // delete temporary volume field
2070     if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2071     if( totVol <= 0)
2072         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2073
2074     return integrale/totVol;
2075 }
2076
2077 /*!  Return L2 Norm  of the field.
2078  *   Cannot be applied to a field with a support on nodes.
2079  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2080  */
2081 template <class T, class INTERLACING_TAG>
2082 double FIELD<T, INTERLACING_TAG>::normL2(const FIELD<double, FullInterlace> * p_field_volume) const
2083 {
2084     _checkNormCompatibility(p_field_volume); // may throw exception
2085     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2086     if(!p_field_volume) // if the user don't supply the volume
2087         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃÂ©mentation dans mesh]
2088
2089     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2090     const double* vol=p_field_size->getValue();
2091     const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
2092
2093     const T * value     = NULL;
2094     ArrayNo * myArray   = NULL;
2095     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2096       value = getValue();
2097     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ){
2098       myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2099       value   = myArray->getPtr();
2100     }
2101     else {
2102       myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2103       value   = myArray->getPtr();
2104     }
2105
2106     double totVol=0.0;
2107     const double* p_vol=vol;
2108     for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2109         totVol+=*p_vol;
2110
2111     double integrale=0.0;
2112     for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
2113         for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
2114             integrale += static_cast<double>((*value) * (*value)) * (*p_vol);
2115
2116     if(!p_field_volume) // if the user didn't supply the volume
2117         delete p_field_size; // delete temporary volume field
2118     if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2119     if( totVol <= 0)
2120         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2121
2122     return integrale/totVol;
2123 }
2124
2125 /*!  Return L1 Norm  of the field's component.
2126  *   Cannot be applied to a field with a support on nodes.
2127  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2128  */
2129 template <class T, class INTERLACING_TAG>
2130 double FIELD<T, INTERLACING_TAG>::normL1(int component,
2131                                          const FIELD<double, FullInterlace > * p_field_volume) const
2132 {
2133     _checkNormCompatibility(p_field_volume); // may throw exception
2134     if ( component<1 || component>getNumberOfComponents() )
2135         throw MEDEXCEPTION(STRING("FIELD<T,INTERLACING_TAG>::normL2() : The component argument should be between 1 and the number of components"));
2136
2137     const FIELD<double,FullInterlace> * p_field_size=p_field_volume;
2138     if(!p_field_volume) // if the user don't supply the volume
2139         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃÂ©mentation dans mesh]
2140
2141     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2142     const double* vol=p_field_size->getValue();
2143     const T * value     = NULL;
2144     ArrayNo * myArray   = NULL;
2145     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2146       value = getColumn(component);
2147     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2148       myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2149       value   = myArray->getColumn(component);
2150     }
2151     else {
2152       myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2153       value   = myArray->getColumn(component);
2154     }
2155
2156     const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
2157
2158     double integrale=0.0;
2159     double totVol=0.0;
2160     for (; value!=lastvalue ; ++value ,++vol)
2161     {
2162         integrale += std::abs( static_cast<double>(*value) ) * (*vol);
2163         totVol+=*vol;
2164     }
2165
2166     if(!p_field_volume) // if the user didn't supply the volume
2167         delete p_field_size; // delete temporary volume field
2168     if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2169     if( totVol <= 0)
2170         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2171
2172     return integrale/totVol;
2173 }
2174
2175 /*!  Return L1 Norm  of the field.
2176  *   Cannot be applied to a field with a support on nodes.
2177  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2178  */
2179 template <class T, class INTERLACING_TAG>
2180 double FIELD<T, INTERLACING_TAG>::normL1(const FIELD<double, FullInterlace> * p_field_volume) const
2181 {
2182     _checkNormCompatibility(p_field_volume); // may throw exception
2183     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2184     if(!p_field_volume) // if the user don't supply the volume
2185         p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃÂ©mentation dans mesh]
2186
2187     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2188     const double* vol=p_field_size->getValue();
2189     const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
2190
2191    const T * value     = NULL;
2192     ArrayNo * myArray   = NULL;
2193     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2194       value = getValue();
2195     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2196       myArray = ArrayConvert2No( *( dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() ) ));
2197       value   = myArray->getPtr();
2198     }
2199     else {
2200       myArray = ArrayConvert( *( dynamic_cast< ArrayFull * > ( getArrayNoGauss() ) ));
2201       value   = myArray->getPtr();
2202     }
2203
2204     double totVol=0.0;
2205     const double* p_vol=vol;
2206     for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2207         totVol+=*p_vol;
2208
2209     double integrale=0.0;
2210     for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
2211         for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
2212             integrale += std::abs( static_cast<double>(*value) ) * (*p_vol);
2213
2214     if(!p_field_volume) // if the user didn't supply the volume
2215         delete p_field_size; // delete temporary volume field
2216     if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE ) delete myArray;
2217     if( totVol <= 0)
2218         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2219
2220     return integrale/totVol;
2221 }
2222
2223 /*! Return a new field (to deallocate with delete) lying on subSupport that is included by
2224  *   this->_support with corresponding values extracting from this->_value.
2225  */
2226 template <class T, class INTERLACING_TAG>
2227 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
2228 {
2229   if(!subSupport->belongsTo(*_support))
2230     throw MEDEXCEPTION("FIELD<T>::extract : subSupport not included in this->_support !");
2231   if(_support->isOnAllElements() && subSupport->isOnAllElements())
2232     return new FIELD<T, INTERLACING_TAG>(*this);
2233
2234   FIELD<T, INTERLACING_TAG> *ret = new FIELD<T, INTERLACING_TAG>(subSupport,
2235                                                                  _numberOfComponents);
2236
2237   if(!ret->_value)
2238     throw MEDEXCEPTION("FIELD<T>::extract : unvalid support detected !");
2239
2240   T* valuesToSet=(T*)ret->getValue();
2241
2242   int nbOfEltsSub=subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2243   const int *eltsSub=subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS);
2244   T* tempVals=new T[_numberOfComponents];
2245   for(int i=0;i<nbOfEltsSub;i++)
2246     {
2247       if(!getValueOnElement(eltsSub[i],tempVals))
2248         throw MEDEXCEPTION("Problem in belongsTo function !!!");
2249       for(int j=0;j<_numberOfComponents;j++)
2250         valuesToSet[i*_numberOfComponents+j]=tempVals[j];
2251     }
2252   delete [] tempVals;
2253
2254   ret->copyGlobalInfo(*this);
2255   return ret;
2256 }
2257
2258 /*!
2259   Constructor with parameters; the object is set via a file and its associated
2260   driver. For the moment only the MED_DRIVER is considered and if the last two
2261   argument (iterationNumber and orderNumber) are not set; their default value
2262   is -1. If the field fieldDriverName with the iteration number
2263   iterationNumber and the order number orderNumber does not exist in the file
2264   fieldDriverName; the constructor raises an exception.
2265 */
2266 template <class T, class INTERLACING_TAG>
2267 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
2268                                  driverTypes driverType,
2269                                  const string & fileName/*=""*/,
2270                                  const string & fieldDriverName/*=""*/,
2271                                  const int iterationNumber,
2272                                  const int orderNumber) throw (MEDEXCEPTION)
2273 {
2274   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) : ";
2275
2276   int current;
2277
2278   BEGIN_OF(LOC);
2279
2280   init();
2281
2282   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2283   ASSERT(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2284   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
2285
2286   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2287   ASSERT(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2288   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2289
2290   _support = Support;
2291   //A.G. Addings for RC
2292   if(_support)
2293     _support->addReference();
2294   // OCC 10/03/2006 -- According to the rules defined with help of
2295   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2296   // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2297   // FIELD template - MSVC++ 2003 compiler generated an error here.
2298   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2299   _value = NULL;
2300
2301   _iterationNumber = iterationNumber;
2302   _time = 0.0;
2303   _orderNumber = orderNumber;
2304
2305   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::MED_LECT);
2306
2307   _drivers[current]->open();
2308   _drivers[current]->read();
2309   _drivers[current]->close();
2310
2311   END_OF(LOC);
2312 }
2313
2314 /*!
2315   This constructor, at least, allows to create a FIELD without creating any
2316   SUPPORT then without having to load a MESH object, a support is created. It
2317   provides the meshName related mesh but doesn't not set a mesh in the created
2318   support.
2319 */
2320 template <class T, class INTERLACING_TAG>
2321 FIELD<T,INTERLACING_TAG>::FIELD(driverTypes driverType,
2322                                 const string & fileName,
2323                                 const string & fieldDriverName,
2324                                 const int iterationNumber,
2325                                 const int orderNumber)
2326   throw (MEDEXCEPTION) :FIELD_()
2327 {
2328   int current;
2329   const char * LOC ="FIELD<T,INTERLACING_TAG>::FIELD( driverTypes driverType, const string & fileName, string & fieldDriverName, int iterationNumber, int orderNumber) : ";
2330   BEGIN_OF(LOC);
2331
2332   init();
2333
2334   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2335   ASSERT(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2336   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
2337
2338   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2339   ASSERT(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2340   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2341
2342   _support = (SUPPORT *) NULL;
2343   // OCC 10/03/2006 -- According to the rules defined with help of
2344   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2345   // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2346   // FIELD template - MSVC++ 2003 compiler generated an error here.
2347   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2348   _value = NULL;
2349
2350   _iterationNumber = iterationNumber;
2351   _time = 0.0;
2352   _orderNumber = orderNumber;
2353
2354   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::MED_LECT);
2355
2356   _drivers[current]->open();
2357   _drivers[current]->read();
2358   _drivers[current]->close();
2359
2360   END_OF(LOC);
2361 }
2362
2363 /*!
2364   Destructor.
2365 */
2366 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::~FIELD()
2367 {
2368   BEGIN_OF(" Destructeur FIELD<T, INTERLACING_TAG>::~FIELD()");
2369   SCRUTE(this);
2370   if (_value) delete _value;
2371   locMap::const_iterator it;
2372   for ( it = _gaussModel.begin();it != _gaussModel.end(); it++ )
2373     delete (*it).second;
2374
2375   END_OF(" Destructeur FIELD<T,INTERLACING_TAG>::~FIELD()");
2376 }
2377
2378 /*!
2379
2380 */
2381 template <class T, class INTERLACING_TAG>
2382 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)
2383 {
2384   const char* LOC = "FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)" ;
2385   BEGIN_OF(LOC);
2386
2387   _numberOfComponents = NumberOfComponents ;
2388   if (_componentsTypes == NULL)
2389     _componentsTypes = new int[NumberOfComponents] ;
2390   if (_componentsNames == NULL)
2391     _componentsNames = new string[NumberOfComponents];
2392   if (_componentsDescriptions == NULL)
2393     _componentsDescriptions = new string[NumberOfComponents];
2394   if (_componentsUnits == NULL)
2395     _componentsUnits = new UNIT[NumberOfComponents];
2396   if (_MEDComponentsUnits == NULL)
2397     _MEDComponentsUnits = new string[NumberOfComponents];
2398   for (int i=0;i<NumberOfComponents;i++) {
2399     _componentsTypes[i] = 0 ;
2400   }
2401
2402   try {
2403     // becarefull about the number of gauss point
2404     _numberOfValues = _support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2405     MESSAGE(LOC <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
2406
2407     //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
2408     _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
2409
2410     _isRead = true ;
2411   }
2412   catch (MEDEXCEPTION &ex) {
2413     MESSAGE("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
2414     // OCC 10/03/2006 -- According to the rules defined with help of
2415     // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2416     // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2417     // FIELD template - MSVC++ 2003 compiler generated an error here.
2418     // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2419     _value = NULL;
2420   }
2421
2422   SCRUTE(_value);
2423   END_OF("void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)");
2424 }
2425
2426 /*!
2427
2428 */
2429 template <class T, class INTERLACING_TAG>
2430 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents,
2431                                            const int LengthValue)
2432 {
2433   BEGIN_OF("void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)");
2434
2435   _numberOfComponents = NumberOfComponents ;
2436   if (_componentsTypes == NULL)
2437     _componentsTypes = new int[NumberOfComponents] ;
2438   if (_componentsNames == NULL)
2439     _componentsNames = new string[NumberOfComponents];
2440   if (_componentsDescriptions == NULL)
2441     _componentsDescriptions = new string[NumberOfComponents];
2442   if (_componentsUnits == NULL)
2443     _componentsUnits = new UNIT[NumberOfComponents];
2444   if (_MEDComponentsUnits == NULL)
2445     _MEDComponentsUnits = new string[NumberOfComponents];
2446   for (int i=0;i<NumberOfComponents;i++) {
2447     _componentsTypes[i] = 0 ;
2448   }
2449
2450   MESSAGE("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
2451   _numberOfValues = LengthValue ;
2452
2453   //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
2454   _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
2455
2456   _isRead = true ;
2457
2458   SCRUTE(_value);
2459   END_OF("void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents,const int LengthValue)");
2460 }
2461
2462 /*!
2463
2464 */
2465 template <class T, class INTERLACING_TAG>
2466 void FIELD<T, INTERLACING_TAG>::deallocValue()
2467 {
2468   BEGIN_OF("void FIELD<T, INTERLACING_TAG>::deallocValue()");
2469   _numberOfValues = 0 ;
2470   _numberOfComponents = 0 ;
2471   if (_value != NULL)
2472     delete _value;
2473
2474   END_OF("void FIELD<T, INTERLACING_TAG>::deallocValue()");
2475 }
2476
2477 // -----------------
2478 // Methodes Inline
2479 // -----------------
2480
2481 /*!
2482   Create the specified driver and return its index reference to path to
2483   read or write methods.
2484 */
2485
2486 template <class T, class INTERLACING_TAG>
2487 int FIELD<T, INTERLACING_TAG>::addDriver(driverTypes driverType,
2488                                          const string & fileName/*="Default File Name.med"*/,
2489                                          const string & driverName/*="Default Field Name"*/,
2490                                          MED_EN::med_mode_acces access)
2491 {
2492   //jfa tmp (as last argument has no default value):const char * LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Field Name\",MED_EN::med_mode_acces access) : ";
2493   const char * LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName,const string & driverName,MED_EN::med_mode_acces access) :";//jfa tmp
2494
2495   GENDRIVER * driver;
2496
2497   BEGIN_OF(LOC);
2498
2499   SCRUTE(driverType);
2500
2501   driver = DRIVERFACTORY::buildDriverForField(driverType,fileName,this,access);
2502
2503   _drivers.push_back(driver);
2504
2505   int current = _drivers.size()-1;
2506
2507   _drivers[current]->setFieldName(driverName);
2508
2509   END_OF(LOC);
2510
2511   return current;
2512 }
2513
2514
2515 /*!
2516   Duplicate the given driver and return its index reference to path to
2517   read or write methods.
2518 */
2519 template <class T, class INTERLACING_TAG>
2520 inline int FIELD<T, INTERLACING_TAG>::addDriver (GENDRIVER & driver )
2521 {
2522   const char * LOC = "FIELD<T, INTERLACING_TAG>::addDriver(GENDRIVER &) : ";
2523   int current;
2524
2525   BEGIN_OF(LOC);
2526
2527   // duplicate driver to delete it with destructor !
2528   GENDRIVER * newDriver = driver.copy() ;
2529
2530   _drivers.push_back(newDriver);
2531
2532   current = _drivers.size()-1;
2533   SCRUTE(current);
2534   driver.setId(current);
2535
2536   MESSAGE(LOC << " je suis la 1");
2537   END_OF(LOC);
2538   MESSAGE(LOC << " je suis la 2");
2539
2540   return current ;
2541 };
2542
2543 /*!
2544   Remove the driver referenced by its index.
2545 */
2546 template <class T, class INTERLACING_TAG>
2547 void FIELD<T, INTERLACING_TAG>::rmDriver (int index/*=0*/)
2548 {
2549   const char * LOC = "FIELD<T, INTERLACING_TAG>::rmDriver (int index=0): ";
2550   BEGIN_OF(LOC);
2551
2552   if ( _drivers[index] ) {
2553     //_drivers.erase(&_drivers[index]);
2554     // why not ????
2555     MESSAGE ("detruire");
2556   }
2557   else
2558     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2559                                      << "The <index given is invalid, index must be between  0 and  |"
2560                                      << _drivers.size()
2561                                      )
2562                           );
2563
2564   END_OF(LOC);
2565 }
2566
2567 /*!
2568   Read FIELD in the file specified in the driver given by its index.
2569 */
2570 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(int index/*=0*/)
2571 {
2572   const char * LOC = "FIELD<T, INTERLACING_TAG>::read(int index=0) : ";
2573   BEGIN_OF(LOC);
2574
2575   if ( _drivers[index] ) {
2576     _drivers[index]->open();
2577     _drivers[index]->read();
2578     _drivers[index]->close();
2579   }
2580   else
2581     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2582                                      << "The index given is invalid, index must be between  0 and |"
2583                                      << _drivers.size()
2584                                      )
2585                           );
2586   END_OF(LOC);
2587 }
2588
2589 /*!
2590   Write FIELD in the file specified in the driver given by its index.
2591 */
2592 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(int index/*=0*/, const string & driverName /*= ""*/)
2593 {
2594   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
2595   BEGIN_OF(LOC);
2596
2597   if( _drivers[index] ) {
2598     _drivers[index]->open();
2599     if (driverName != "") _drivers[index]->setFieldName(driverName);
2600     _drivers[index]->write();
2601     _drivers[index]->close();
2602   }
2603   else
2604     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2605                                      << "The index given is invalid, index must be between  0 and |"
2606                                      << _drivers.size()
2607                                      )
2608                           );
2609   END_OF(LOC);
2610 }
2611
2612 /*!
2613   Write FIELD in the file specified in the driver given by its index. Use this
2614   method for ASCII drivers (e.g. VTK_DRIVER)
2615 */
2616 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
2617 {
2618   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
2619   BEGIN_OF(LOC);
2620
2621   if( _drivers[index] ) {
2622     _drivers[index]->openAppend();
2623     if (driverName != "") _drivers[index]->setFieldName(driverName);
2624     _drivers[index]->writeAppend();
2625     _drivers[index]->close();
2626   }
2627   else
2628     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
2629                                      << "The index given is invalid, index must be between  0 and |"
2630                                      << _drivers.size()
2631                                      )
2632                           );
2633   END_OF(LOC);
2634 }
2635
2636 /*!
2637   \internal
2638   Write FIELD with the driver which is equal to the given driver.
2639
2640   Use by MED object.
2641 */
2642 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(const GENDRIVER & genDriver)
2643 {
2644   const char * LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
2645   BEGIN_OF(LOC);
2646
2647   for (unsigned int index=0; index < _drivers.size(); index++ )
2648     if ( *_drivers[index] == genDriver ) {
2649       _drivers[index]->open();
2650       _drivers[index]->write();
2651       _drivers[index]->close();
2652     }
2653
2654   END_OF(LOC);
2655
2656 }
2657
2658 /*!
2659   \internal
2660   Write FIELD with the driver which is equal to the given driver.
2661
2662   Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
2663 */
2664 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(const GENDRIVER & genDriver)
2665 {
2666   const char * LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
2667   BEGIN_OF(LOC);
2668
2669   for (unsigned int index=0; index < _drivers.size(); index++ )
2670     if ( *_drivers[index] == genDriver ) {
2671       _drivers[index]->openAppend();
2672       _drivers[index]->writeAppend();
2673       _drivers[index]->close();
2674     }
2675
2676   END_OF(LOC);
2677
2678 }
2679
2680 /*!
2681   \internal
2682   Read FIELD with the driver which is equal to the given driver.
2683
2684   Use by MED object.
2685 */
2686 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(const GENDRIVER & genDriver)
2687 {
2688   const char * LOC = " FIELD<T, INTERLACING_TAG>::read(const GENDRIVER &) : ";
2689   BEGIN_OF(LOC);
2690
2691   for (unsigned int index=0; index < _drivers.size(); index++ )
2692     if ( *_drivers[index] == genDriver ) {
2693       _drivers[index]->open();
2694       _drivers[index]->read();
2695       _drivers[index]->close();
2696     }
2697
2698   END_OF(LOC);
2699
2700 }
2701
2702 /*!
2703   Fills in already allocated retValues array the values related to eltIdInSup.
2704   If the element does not exist in this->_support false is returned, true otherwise.
2705  */
2706 template <class T, class INTERLACING_TAG>
2707 bool FIELD<T, INTERLACING_TAG>::getValueOnElement(int eltIdInSup,T* retValues)
2708   const throw (MEDEXCEPTION)
2709 {
2710
2711   if(eltIdInSup<1)
2712     return false;
2713   if(_support->isOnAllElements())
2714     {
2715       int nbOfEltsThis=_support->getMesh()->getNumberOfElements(_support->getEntity(),MED_EN::MED_ALL_ELEMENTS);
2716       if(eltIdInSup>nbOfEltsThis)
2717         return false;
2718       const T* valsThis=getValue();
2719       for(int j=0;j<_numberOfComponents;j++)
2720         retValues[j]=valsThis[(eltIdInSup-1)*_numberOfComponents+j];
2721       return true;
2722     }
2723   else
2724     {
2725       int nbOfEltsThis=_support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2726       const int *eltsThis=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
2727       int iThis;
2728       bool found=false;
2729       for(iThis=0;iThis<nbOfEltsThis && !found;)
2730         if(eltsThis[iThis]==eltIdInSup)
2731           found=true;
2732         else
2733           iThis++;
2734       if(!found)
2735         return false;
2736       const T* valsThis=getValue();
2737       for(int j=0;j<_numberOfComponents;j++)
2738         retValues[j]=valsThis[iThis*_numberOfComponents+j];
2739       return true;
2740     }
2741 }
2742
2743 /*!
2744   \if developper
2745   Destroy the MEDARRAY<T> in FIELD and put the new one without copy.
2746   \endif
2747 */
2748 template <class T, class INTERLACING_TAG>
2749 inline void FIELD<T, INTERLACING_TAG>::setArray(MEDMEM_Array_ * Value)
2750   throw (MEDEXCEPTION)
2751 {
2752   if (NULL != _value) delete _value ;
2753   _value=Value ;
2754 }
2755
2756 /*!
2757   \if developper
2758   Return a reference to  the MEDARRAY<T, INTERLACING_TAG> in FIELD.
2759   \endif
2760 */
2761 template <class T, class INTERLACING_TAG>
2762 inline MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() const throw (MEDEXCEPTION)
2763 {
2764   const char * LOC = "MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() : ";
2765   BEGIN_OF(LOC);
2766   END_OF(LOC);
2767   return _value ;
2768 }
2769 template <class T,class INTERLACING_TAG>  inline
2770 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array *
2771 FIELD<T, INTERLACING_TAG>::getArrayGauss() const throw (MEDEXCEPTION)
2772 {
2773   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayGauss() : ";
2774   BEGIN_OF(LOC);
2775
2776   if ( getGaussPresence() )
2777     return static_cast<ArrayGauss *> (_value);
2778   else
2779     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
2780                                  "The field has no Gauss Point"));
2781
2782   END_OF(LOC);
2783
2784 }
2785
2786 template <class T,class INTERLACING_TAG>  inline
2787 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array *
2788 FIELD<T, INTERLACING_TAG>::getArrayNoGauss() const throw (MEDEXCEPTION)
2789 {
2790   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayNoGauss() : ";
2791   BEGIN_OF(LOC);
2792
2793   if ( ! getGaussPresence() )
2794     return static_cast < ArrayNoGauss * > (_value);
2795   else
2796     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
2797                                  "The field has Gauss Point"));
2798
2799   END_OF(LOC);
2800 }
2801
2802
2803 template <class T,class INTERLACING_TAG> inline bool
2804 FIELD<T, INTERLACING_TAG>::getGaussPresence() const throw (MEDEXCEPTION)
2805 {
2806   const char * LOC = "FIELD<T, INTERLACING_TAG>::getGaussPresence() const :";
2807   //BEGIN_OF(LOC);
2808
2809   if (_value != NULL)
2810     return _value->getGaussPresence();
2811   else
2812     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't call getGaussPresence on a null _value"));
2813
2814   //END_OF(LOC);
2815 }
2816
2817 /*!
2818   Return the actual length of the reference to values array returned by getValue.
2819   Take care of number of components and number of Gauss points by geometric type
2820 */
2821 template <class T, class INTERLACING_TAG>
2822 inline int FIELD<T, INTERLACING_TAG>::getValueLength() const
2823   throw (MEDEXCEPTION)
2824 {
2825   if ( getGaussPresence() )
2826     return static_cast<ArrayGauss *>(_value)->getArraySize() ;
2827   else
2828     return static_cast<ArrayNoGauss *>(_value)->getArraySize() ;
2829 }
2830
2831 /*!
2832   Return a reference to values array to read them.
2833 */
2834 template <class T, class INTERLACIN_TAG>
2835 inline const T* FIELD<T, INTERLACIN_TAG>::getValue() const throw (MEDEXCEPTION)
2836 {
2837   const char * LOC ="FIELD<T, INTERLACING_TAG>::getValue() : ";
2838   BEGIN_OF(LOC);
2839   if ( getGaussPresence() )
2840     return static_cast<ArrayGauss *>(_value)->getPtr() ;
2841   else
2842     return static_cast<ArrayNoGauss *>(_value)->getPtr() ;
2843 }
2844 /*!
2845   Return a reference to i^{th} row
2846   of FIELD values array.
2847   If a faster accessor is intended you may use getArray() once,
2848   then MEDMEM_Array accessors.
2849   Be careful if field support is not on all elements getRow
2850   use support->getValIndFromGlobalNumber(i).
2851 */
2852 template <class T,class INTERLACING_TAG> inline
2853 const T*
2854 FIELD<T,INTERLACING_TAG>::getRow(int i) const throw (MEDEXCEPTION)
2855 {
2856   const char * LOC = "FIELD<T,INTERLACING_TAG>::getRow(int i) : ";
2857   //BEGIN_OF(LOC);
2858
2859   int valIndex=-1;
2860   if (_support)
2861     valIndex = _support->getValIndFromGlobalNumber(i);
2862   else
2863     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
2864
2865   //cout << endl << "getRow Valindex : " << valIndex << endl;
2866   if ( getGaussPresence() )
2867     return static_cast<ArrayGauss *>(_value)->getRow(valIndex) ;
2868   else
2869     return static_cast<ArrayNoGauss *>(_value)->getRow(valIndex) ;
2870   //END_OF(LOC);
2871 }
2872
2873 /*!
2874   Return a reference to j^{th} column
2875   of FIELD values array.
2876 */
2877 template <class T,class INTERLACING_TAG> inline const T*
2878 FIELD<T,INTERLACING_TAG>::getColumn(int j) const throw (MEDEXCEPTION)
2879 {
2880   //const char * LOC ="FIELD<T,INTERLACING_TAG>::getColumn(int j) : ";
2881   //BEGIN_OF(LOC);
2882   if ( getGaussPresence() )
2883     return static_cast<ArrayGauss *>(_value)->getColumn(j) ;
2884   else
2885     return static_cast<ArrayNoGauss *>(_value)->getColumn(j) ;
2886 }
2887
2888 /*!
2889   Return the value of i^{th} element and j^{th} component.
2890 */
2891 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJ(int i,int j) const throw (MEDEXCEPTION)
2892 {
2893   const char * LOC = "getValueIJ(..)";
2894   //BEGIN_OF(LOC);
2895   int valIndex=-1;
2896   if (_support)
2897     valIndex = _support->getValIndFromGlobalNumber(i);
2898   else
2899     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
2900
2901   if ( getGaussPresence() )
2902     return static_cast<ArrayGauss *>(_value)->getIJ(valIndex,j) ;
2903   else
2904     return static_cast<ArrayNoGauss *>(_value)->getIJ(valIndex,j) ;
2905 }
2906
2907 /*!
2908   Return the j^{th} component of k^{th} gauss points of i^{th} value.
2909 */
2910 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION)
2911 {
2912   const char * LOC = "getValueIJK(..)";
2913   //BEGIN_OF(LOC);
2914   int valIndex=-1;
2915   if (_support)
2916     valIndex = _support->getValIndFromGlobalNumber(i);
2917   else
2918     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
2919
2920   if ( getGaussPresence() )
2921     return static_cast<ArrayGauss *>(_value)->getIJK(valIndex,j,k) ;
2922   else
2923     return static_cast<ArrayNoGauss *>(_value)->getIJK(valIndex,j,k) ;
2924 }
2925
2926 /*!
2927   Return number of values of a geomertic type in NoInterlaceByType mode
2928 */
2929 template <class T, class INTERLACIN_TAG>
2930 inline int FIELD<T, INTERLACIN_TAG>::getValueByTypeLength(int t) const throw (MEDEXCEPTION)
2931 {
2932   const char * LOC ="getValueByTypeLength() : ";
2933   //BEGIN_OF(LOC);
2934   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
2935     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
2936
2937   if ( getGaussPresence() ) {
2938     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
2939     if (  t < 1 || t > array->getNbGeoType() )
2940       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
2941     return array->getLengthOfType( t );
2942   }
2943   else {
2944     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
2945     if (  t < 1 || t > array->getNbGeoType() )
2946       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
2947     return array->getLengthOfType( t );
2948   }
2949 }
2950
2951 /*!
2952   Return a reference to values array to read them.
2953 */
2954 template <class T, class INTERLACIN_TAG>
2955 inline const T* FIELD<T, INTERLACIN_TAG>::getValueByType(int t) const throw (MEDEXCEPTION)
2956 {
2957   const char * LOC ="getValueByType() : ";
2958   //BEGIN_OF(LOC);
2959   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
2960     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
2961
2962   if ( getGaussPresence() ) {
2963     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
2964     return array->getPtr() + array->getIndex( t );
2965   }
2966   else {
2967     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
2968     return array->getPtr() + array->getIndex( t );
2969   }
2970 }
2971
2972 /*!
2973   Return the value of i^{th} element in indicated type t and j^{th} component.
2974 */
2975 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJByType(int i,int j, int t) const throw (MEDEXCEPTION)
2976 {
2977   const char * LOC = "getValueIJByType(..)";
2978   //BEGIN_OF(LOC);
2979   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
2980     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
2981
2982   if ( getGaussPresence() )
2983     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJByType(i,j,t) ;
2984   else
2985     return static_cast<ArrayNoByType *>(_value)->getIJByType(i,j,t) ;
2986 }
2987
2988 /*!
2989   Return the j^{th} component of k^{th} gauss points of i^{th} value with type t.
2990 */
2991 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION)
2992 {
2993   const char * LOC = "getValueIJKByType(..)";
2994   //BEGIN_OF(LOC);
2995   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
2996     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
2997
2998   if ( getGaussPresence() )
2999     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJKByType(i,j,k,t) ;
3000   else
3001     return static_cast<ArrayNoByType      *>(_value)->getIJKByType(i,j,k,t) ;
3002 }
3003
3004
3005 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGeometricTypes() const throw (MEDEXCEPTION)
3006 {
3007   const char * LOC = "getNumberOfGeometricTypes(..)";
3008   BEGIN_OF(LOC);
3009   if (_support)
3010     return _support->getNumberOfTypes();
3011   else
3012     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3013   END_OF(LOC);
3014 };
3015
3016
3017 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> &
3018 FIELD<T,INTERLACING_TAG>::getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
3019 {
3020   const char * LOC ="getGaussLocalization(MED_EN::medGeometryElement geomElement) : ";
3021   const GAUSS_LOCALIZATION_ * locPtr=0;
3022
3023   locMap::const_iterator it;
3024   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3025         locPtr = (*it).second;
3026         return  *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
3027   }
3028   else
3029     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
3030
3031 };
3032
3033 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> *
3034 FIELD<T,INTERLACING_TAG>::getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
3035 {
3036   const char * LOC ="getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) : ";
3037   const GAUSS_LOCALIZATION_ * locPtr=0;
3038
3039   locMap::const_iterator it;
3040   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3041         locPtr = (*it).second;
3042         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
3043   }
3044   else
3045     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
3046
3047 };
3048
3049 /*!
3050  * \brief Return GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
3051  */
3052 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION_ *
3053 FIELD<T,INTERLACING_TAG>::getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
3054 {
3055   const char * LOC ="getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) : ";
3056
3057   locMap::const_iterator it;
3058   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3059         return (*it).second;
3060   }
3061   else
3062     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type: "<< geomElement ));
3063
3064 };
3065
3066 /*!
3067  * \brief Take onership of GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
3068  */
3069 template <class T,class INTERLACING_TAG> void
3070 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc)
3071 {
3072   locMap::iterator it = _gaussModel.find(geomElement);
3073   if ( it != _gaussModel.end() ) {
3074     delete it->second;
3075     it->second = gaussloc;
3076   }
3077   else {
3078     _gaussModel[ geomElement ] = gaussloc;
3079   }
3080 };
3081
3082
3083 template <class T,class INTERLACING_TAG> void
3084 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG>& gaussloc)
3085 {
3086   locMap::iterator it = _gaussModel.find(geomElement);
3087   if ( it != _gaussModel.end() ) {
3088     delete it->second;
3089     it->second = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
3090   }
3091   else {
3092     _gaussModel[ geomElement ] = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
3093   }
3094 };
3095
3096 /*!
3097   Returns number of Gauss points for this medGeometryElement.
3098
3099   Note :
3100   if there is no GAUSS_LOCALIZATION having this medGeometryElement but
3101   the medGeometryElement exist in the SUPPORT, getNumberOfGaussPoints
3102   return 1
3103 */
3104 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) const
3105   throw (MEDEXCEPTION)
3106 {
3107   const char * LOC ="getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
3108   const GAUSS_LOCALIZATION_ * locPtr=0;
3109
3110   locMap::const_iterator it;
3111   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
3112         locPtr = (*it).second;
3113         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr)->getNbGauss();
3114   }
3115   else
3116     if (_support)
3117       try {
3118         if ( _support->getNumberOfElements(geomElement) ) return 1;
3119       } catch ( MEDEXCEPTION & ex) {
3120         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<< "GeometricType not found !" )) ;
3121       }
3122     else
3123       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3124
3125   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Should never execute this!" ));
3126
3127 }
3128
3129 /*!
3130   Returns number of Gauss points for each geometric type.
3131
3132   Note :
3133   if there is no gauss points whatever the geometric type is
3134   it returns an exception. (renvoyer un tableau de 1 ?)
3135 */
3136 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints() const
3137   throw (MEDEXCEPTION)
3138 {
3139   const char * LOC ="const int * getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
3140
3141   if (_value)
3142     if ( getGaussPresence() ) {
3143       return static_cast<ArrayGauss *>(_value)->getNbGaussGeo()+1;
3144     } else
3145       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"value hasn't Gauss points " ));
3146
3147     else
3148       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Value not defined" ));
3149 }
3150
3151 /*!
3152   Returns number of Gauss points for element n°i.
3153   The i index is a global index (take care of previous element
3154   on different geometric type).
3155 */
3156 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNbGaussI(int i) const throw (MEDEXCEPTION)
3157 {
3158   const char * LOC = "getNbGaussI(..)";
3159 //   BEGIN_OF(LOC);
3160
3161   int valIndex=-1;
3162   if (_support)
3163     valIndex = _support->getValIndFromGlobalNumber(i);
3164   else
3165     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3166
3167   if (_value)
3168    if ( getGaussPresence() )
3169      return static_cast<ArrayGauss *>(_value)->getNbGauss(valIndex) ;
3170    else
3171      return static_cast<ArrayNoGauss *>(_value)->getNbGauss(valIndex) ;
3172  else
3173    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_value not defined" ));
3174 //   END_OF(LOC);
3175 };
3176
3177 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfElements() const throw (MEDEXCEPTION)
3178 {
3179   const char * LOC = "getNumberOfElements(..)";
3180   BEGIN_OF(LOC);
3181   if (_support)
3182     return _support->getNumberOfElements();
3183   else
3184     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3185    END_OF(LOC);
3186 };
3187
3188 template <class T,class INTERLACING_TAG> const MED_EN::medGeometryElement  * FIELD<T,INTERLACING_TAG>::getGeometricTypes()  const throw (MEDEXCEPTION)
3189 {
3190   const char * LOC = "getGeometricTypes(..)";
3191   BEGIN_OF(LOC);
3192   if (_support)
3193     return _support->getTypes();
3194   else
3195     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3196    END_OF(LOC);
3197 };
3198 template <class T,class INTERLACING_TAG> bool  FIELD<T,INTERLACING_TAG>::isOnAllElements() const throw (MEDEXCEPTION)
3199 {
3200   const char * LOC = "isOnAllElements(..)";
3201   BEGIN_OF(LOC);
3202   if (_support)
3203     return _support->isOnAllElements();
3204   else
3205     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3206   END_OF(LOC);
3207 };
3208
3209
3210 /*!
3211   Copy new values array in FIELD according to the given mode.
3212
3213   Array must have right size. If not results are unpredicable.
3214 */
3215 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValue( T* value) throw (MEDEXCEPTION)
3216 {
3217   if ( getGaussPresence() )
3218     static_cast<ArrayGauss *>(_value)->setPtr(value) ;
3219   else
3220     static_cast<ArrayNoGauss *>(_value)->setPtr(value) ;
3221 }
3222
3223 /*!
3224   Update values array in the j^{th} row of FIELD values array with the given ones and
3225   according to specified mode.
3226 */
3227 template <class T,class INTERLACING_TAG>
3228 inline void FIELD<T,INTERLACING_TAG>::setRow( int i, T* value) throw (MEDEXCEPTION)
3229 {
3230   const char * LOC = "FIELD<T,INTERLACING_TAG>::setRow(int i, T* value) : ";
3231   int valIndex=i;
3232   if (_support)
3233     valIndex = _support->getValIndFromGlobalNumber(i);
3234   else
3235     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
3236
3237   if ( getGaussPresence() )
3238     static_cast<ArrayGauss *>(_value)->setRow(valIndex, value) ;
3239   else
3240     static_cast<ArrayNoGauss *>(_value)->setRow(valIndex, value) ;
3241 }
3242
3243 /*!
3244   Update values array in the j^{th} column of FIELD values array with the given ones and
3245   according to specified mode.
3246 */
3247 template <class T,class INTERLACING_TAG>
3248 inline void FIELD<T,INTERLACING_TAG>::setColumn( int j, T* value) throw (MEDEXCEPTION)
3249 {
3250   if ( getGaussPresence() )
3251     static_cast<ArrayGauss *>(_value)->setColumn(j, value) ;
3252   else
3253     static_cast<ArrayNoGauss *>(_value)->setColumn(j, value) ;
3254 }
3255
3256 /*!
3257   Set the value of i^{th} element and j^{th} component with the given one.
3258 */
3259 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) throw (MEDEXCEPTION)
3260 {
3261   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
3262   int valIndex=-1;
3263   if (_support)
3264     valIndex = _support->getValIndFromGlobalNumber(i);
3265   else
3266     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
3267
3268   if ( getGaussPresence() )
3269     static_cast<ArrayGauss *>(_value)->setIJ(valIndex,j,value) ;
3270   else
3271     static_cast<ArrayNoGauss *>(_value)->setIJ(valIndex,j,value) ;
3272 }
3273
3274 /*!
3275   Set the value of i^{th} element, j^{th} component and k^{th} gauss point with the given one.
3276 */
3277 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION)
3278 {
3279   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
3280   int valIndex=-1;
3281   if (_support)
3282     valIndex = _support->getValIndFromGlobalNumber(i);
3283   else
3284     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
3285
3286   if ( getGaussPresence() )
3287     static_cast<ArrayGauss *>(_value)->setIJK(valIndex,j,k,value) ;
3288   else
3289     static_cast<ArrayNoGauss *>(_value)->setIJK(valIndex,j,k,value) ;
3290 }
3291
3292 /*!
3293   Set the value of i^{th} element and j^{th} component with the given one.
3294 */
3295 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION)
3296 {
3297   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) : ";
3298   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3299     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3300
3301   if ( getGaussPresence() )
3302     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJByType(i,j,t,value) ;
3303   else
3304     return static_cast<ArrayNoByType *>(_value)->setIJByType(i,j,t,value) ;
3305 }
3306
3307 /*!
3308   Set the value of component of k^{th} gauss points of i^{th} element and j^{th} component with the given one.
3309 */
3310 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION)
3311 {
3312   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int t, int k, T value) : ";
3313   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3314     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3315
3316   if ( getGaussPresence() )
3317     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJKByType(i,j,k,t,value) ;
3318   else
3319     return static_cast<ArrayNoByType *>(_value)->setIJKByType(i,j,k,t,value) ;
3320 }
3321
3322 /*
3323   METHODS
3324 */
3325
3326 /*!
3327   Fill values array with volume values.
3328 */
3329 template <class T, class INTERLACING_TAG>
3330 void FIELD<T, INTERLACING_TAG>::getVolume() const throw (MEDEXCEPTION)
3331 {
3332   const char * LOC = "FIELD<double>::getVolume() const : ";
3333   BEGIN_OF(LOC);
3334
3335   // The field has to be initilised by a non empty support and a
3336   // number of components = 1 and its value type has to be set to MED_REEL64
3337   // (ie a FIELD<double>)
3338
3339   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
3340       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"));
3341
3342   END_OF(LOC);
3343 }
3344
3345 /*!
3346   Fill values array with area values.
3347 */
3348 template <class T, class INTERLACING_TAG>
3349 void FIELD<T, INTERLACING_TAG>::getArea() const throw (MEDEXCEPTION)
3350 {
3351   const char * LOC = "FIELD<double>::getArea() const : ";
3352   BEGIN_OF(LOC);
3353
3354   // The field has to be initilised by a non empty support and a
3355   // number of components = 1 and its value type has to be set to MED_REEL64
3356   // (ie a FIELD<double>)
3357
3358   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
3359       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"));
3360
3361   END_OF(LOC);
3362 }
3363
3364 /*!
3365   Fill values array with length values.
3366 */
3367 template <class T, class INTERLACING_TAG>
3368 void FIELD<T, INTERLACING_TAG>::getLength() const throw (MEDEXCEPTION)
3369 {
3370   const char * LOC = "FIELD<double>::getLength() const : ";
3371   BEGIN_OF(LOC);
3372
3373   // The field has to be initilised by a non empty support and a
3374   // number of components = 1 and its value type has to be set to MED_REEL64
3375   // (ie a FIELD<double>)
3376
3377   if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_EN::MED_REEL64))
3378       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"));
3379
3380   END_OF(LOC);
3381 }
3382
3383 /*!
3384   Fill values array with normal values.
3385 */
3386 template <class T, class INTERLACING_TAG>
3387 void FIELD<T, INTERLACING_TAG>::getNormal() const throw (MEDEXCEPTION)
3388 {
3389   const char * LOC = "FIELD<double>::getNormal() const : ";
3390   BEGIN_OF(LOC);
3391
3392   // The field has to be initilised by a non empty support and a
3393   // number of components = 1 and its value type has to be set to MED_REEL64
3394   // (ie a FIELD<double>)
3395
3396   if (_support == (SUPPORT *) NULL)
3397       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"));
3398
3399   int dim_space = _support->getMesh()->getSpaceDimension();
3400
3401   if ((_numberOfComponents != dim_space) || (_valueType != MED_EN::MED_REEL64))
3402       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"));
3403
3404   END_OF(LOC);
3405 }
3406
3407 /*!
3408   Fill values array with barycenter values.
3409 */
3410 template <class T, class INTERLACING_TAG>
3411 void FIELD<T, INTERLACING_TAG>::getBarycenter() const throw (MEDEXCEPTION)
3412 {
3413   const char * LOC = "FIELD<double>::getBarycenter() const : ";
3414   BEGIN_OF(LOC);
3415
3416   // The field has to be initilised by a non empty support and a number of
3417   //components = space dimension and its value type has to be set to MED_REEL64
3418   // (ie a FIELD<double>)
3419
3420   if (_support == (SUPPORT *) NULL)
3421       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"));
3422
3423   int dim_space = _support->getMesh()->getSpaceDimension();
3424
3425   if ((_numberOfComponents != dim_space) || (_valueType != MED_EN::MED_REEL64))
3426       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"));
3427
3428   END_OF(LOC);
3429 }
3430
3431 /*!
3432   Fill array by using T_Analytic.
3433   WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
3434   Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
3435  */
3436 template <class T, class INTERLACING_TAG>
3437 void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION)
3438 {
3439   const char * LOC = "void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) : ";
3440   int i,j;
3441   if (_support == (SUPPORT *) NULL)
3442       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No Support defined."));
3443
3444   MESH * mesh = _support->getMesh();
3445   int spaceDim = mesh->getSpaceDimension();
3446   const double * coord;
3447
3448   const double * bary;
3449   FIELD<double,FullInterlace> * barycenterField=0;
3450
3451   double ** xyz=new double* [spaceDim];
3452   bool deallocateXyz=false;
3453   if(_support->getEntity()==MED_EN::MED_NODE)
3454     {
3455       if (_support->isOnAllElements())
3456         {
3457           coord=mesh->getCoordinates(MED_EN::MED_NO_INTERLACE);
3458           for(i=0; i<spaceDim; i++)
3459             xyz[i]=(double *)coord+i*_numberOfValues;
3460         }
3461       else
3462         {
3463           coord = mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
3464           const int * nodesNumber=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
3465           for(i=0; i<spaceDim; i++)
3466             xyz[i]=new double[_numberOfValues];
3467           deallocateXyz=true;
3468           for(i=0;i<_numberOfValues;i++)
3469             {
3470               for(j=0;j<spaceDim;j++)
3471                 xyz[j][i]=coord[(nodesNumber[i]-1)*spaceDim+j];
3472             }
3473         }
3474     }
3475   else
3476     {
3477       barycenterField = mesh->getBarycenter(_support);
3478       bary=barycenterField->getValue();
3479       for(i=0; i<spaceDim; i++)
3480         xyz[i]=(double *)(bary+i*_numberOfValues);
3481     }
3482   T* valsToSet=(T*)getValue();
3483   double *temp=new double[spaceDim];
3484   for(i=0;i<_numberOfValues;i++)
3485   {
3486     for(j=0;j<spaceDim;j++)
3487       temp[j]=xyz[j][i];
3488     f(temp,valsToSet+i*_numberOfComponents);
3489   }
3490   delete [] temp;
3491   if(barycenterField)
3492     delete barycenterField;
3493   if(deallocateXyz)
3494     for(j=0;j<spaceDim;j++)
3495       delete [] xyz[j];
3496   delete [] xyz;
3497 }
3498 /*!
3499   Execute a function on _values on 'this' and put the result on a newly created field that has to be deallocated.
3500   WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
3501   Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
3502  */
3503                                                                          template <class T, class INTERLACING_TAG>
3504                                                                          FIELD<T,INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::execFunc(int nbOfComponents, myFuncType2 f) throw (MEDEXCEPTION)
3505 {
3506   FIELD<T,INTERLACING_TAG> *ret=new FIELD<T,INTERLACING_TAG>(_support,nbOfComponents);
3507   const T* valsInput=getValue();
3508   T* valsOutPut=(T*)ret->getValue();
3509   int i;
3510   for(i=0;i<_numberOfValues;i++)
3511     f(valsInput+i*_numberOfComponents,valsOutPut+i*nbOfComponents);
3512   return ret;
3513 }
3514
3515 }//End namespace MEDMEM
3516
3517 #endif /* FIELD_HXX */