src/MEDMEM/MEDMEM_Field.hxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21
  22 /*
  23  File Field.hxx
  24 */
  25
  26 #ifndef FIELD_HXX
  27 #define FIELD_HXX
  28
  29 #include "MEDMEM.hxx"
  30
  31 #include "MEDMEM_Utilities.hxx"
  32 #include "MEDMEM_MedVersion.hxx"
  33 #include "MEDMEM_Exception.hxx"
  34 #include "MEDMEM_define.hxx"
  35 #include "MEDMEM_Support.hxx"
  36 #include "MEDMEM_Unit.hxx"
  37 #include "MEDMEM_nArray.hxx"
  38 #include "MEDMEM_GenDriver.hxx"
  39 #include "MEDMEM_RCBase.hxx"
  40 #include "MEDMEM_ArrayInterface.hxx"
  41 #include "MEDMEM_SetInterlacingType.hxx"
  42 #include "MEDMEM_FieldForward.hxx"
  43 #include "MEDMEM_GaussLocalization.hxx"
  44 #include "InterpKernelGaussCoords.hxx"
  45 #include "PointLocator.hxx"
  46
  47 #include <vector>
  48 #include <map>
  49 #include <algorithm>
  50 #include <memory>
  51 #include <math.h>
  52 #include <cmath>
  53 #include <float.h>
  54
  55 namespace MEDMEM {
  56
  57   template<class T>
  58   struct MinMax {
  59   };
  60
  61   template<>
  62   struct MinMax<double> {
  63     static double getMin() { return DBL_MIN; }
  64     static double getMax() { return DBL_MAX; }
  65   };
  66
  67   template<>
  68   struct MinMax<int> {
  69     static int getMin() { return INT_MIN; }
  70     static int getMax() { return INT_MAX; }
  71   };
  72
  73   template < typename T > struct SET_VALUE_TYPE {
  74     static const MED_EN::med_type_champ _valueType = MED_EN::MED_UNDEFINED_TYPE;};
  75   template < > struct SET_VALUE_TYPE<double> {
  76     static const MED_EN::med_type_champ _valueType = MED_EN::MED_REEL64; };
  77   template < > struct SET_VALUE_TYPE<int> {
  78     static const MED_EN::med_type_champ _valueType = MED_EN::MED_INT32; };
  79
  80         /*!\defgroup FIELD_io Reading and writing files
  81
  82 Fields can be read or written to/from MED files.
  83
  84 \par Reading fields
  85
  86  For reading a field a typical use consists in :
  87 - reading the mesh associated on which the field lies
  88 - read the field, specifying its time step and order number
  89
  90 As an example :
  91 \verbatim
  92 //reading mesh from file
  93 MESH mesh(MED_DRIVER, "file.med", "my_Mesh");
  94 //reading the field from the file
  95 FIELD<double> field(group,MED_DRIVER,"file.med","my_Field",1,1,&mesh);
  96 \endverbatim
  97
  98 It is also possible to read a field without specifying its support. In this case, the field constructor
  99 creates a support with no link to the initial mesh:
 100 \verbatim
 101 FIELD<double> field(MED_DRIVER, "file.med", "myField",1,1);
 102 SUPPORT* support= field->getSupport();
 103 \endverbatim
 104
 105 See also \ref FIELD_constructors
 106
 107 \par Writing fields
 108
 109 When it comes to write fields, it is enough to call write() method.
 110 A typical use will be :
 111
 112 \verbatim
 113 mesh.write(MED_DRIVER, "myResultFile.med");
 114 field.write(MED_DRIVER, "myResultFile.med");
 115 \endverbatim
 116
 117 \defgroup FIELD_constructors
 118
 119 The different field constructors correspond to the two main
 120 ways a field is used :
 121 - either it is read from a file to be consulted,
 122 - or it can be created from scratch with a link to a support on which the values will be built.
 123
 124 \defgroup FIELD_algo Numerical operations on fields
 125 This section groups together the different operators that enable the user to
 126 treat the FIELD objects as high-level numerical arrays, giving operators for
 127 numerical treatment (overloading of basic operators, algorithms, etc...)
 128
 129 \defgroup FIELD_getset Basic Get/Set operations
 130
 131 This sections groups together the basic operations
 132 that describe access to all the elements constitutive of the description of the field :
 133 - name (compulsory),
 134 - time iteration number(compulsory),
 135 - inner loop iteration number(compulsory),
 136 - time(compulsory),
 137 - description(optional),
 138 - number of components(compulsory),
 139 - components names(optional),
 140 - components description(optional).
 141
 142 Some of these items are compulsory because they are essential to the field in order to define
 143 its structure or to be identified inside a MED file during the write process. The other ones
 144 are there for additional information and can be overlooked if not necessary.
 145
 146 When creating a field by reading a file, all the parameters are set according to the file
 147 data and can be consulted via the get methods. When creating a file from scratch, the
 148  name and number of components are set by the constructor, but the other items have to be
 149  set via the setXXX methods.
 150
 151 \defgroup FIELD_gauss Gauss points
 152
 153 In MED, it is possible to declare a Gauss model that
 154 describes the location of Gauss points in a reference cell.
 155 This Gauss model definition is contained in the
 156 \a GAUSS_LOCALIZATION class. A \a GAUSS_LOCALIZATION object
 157 is associated to a field and to a type.
 158
 159 It is not permitted to define a Gauss model in a polygonal
 160 or polyhedric element.
 161
 162 The Gauss model can be :
 163 - loaded from a MED file,
 164 - written to a MED file,
 165 - used to define a FIELD with multiple localizations per element.
 166
 167 \section gauss_constructors Constructing a Gauss Model
 168
 169 A Gauss model can be constructed with the following constructor :
 170 \param locName defines a name associated with the gauss model
 171 \param typeGeo names the type to which the Gauss model is assocaited
 172 \param nGauss defines the number of Gauss points
 173 \param cooRef defines an array giving the coordinates of the nodes of the reference element (dimension : spaceDimension * number of nodes for type \a typeGeo)
 174 \param cooGauss defines an array giving the coordinates of the nodes of the Gauss points (dimension : spaceDimension * \a nGauss
 175 )
 176 \param wg weights associated with each Gauss point (dimension : \a nGauss)
 177
 178 Example : in 2D, a Gauss model definition for a triangle
 179 would be written as :
 180
 181 \code
 182 string locname("gauss model");
 183 double cooRef[6] ={0.0, 0.0, 1.0, 0.0, 0.0, 1.0};
 184 double cooGauss[6]={0.2, 0.2, 0.8, 0.1, 0.1, 0.8};
 185 double wg[3]={0.3334, 0.3334, 0.3334};
 186 GAUSS_LOCALIZATION model(locname,
 187                          MED_EN::MED_TRIA3,
 188                          3,
 189                          cooRef,
 190                          cooGauss,
 191                          wg);
 192 \endcode
 193
 194 */
 195
 196
 197 /*!
 198
 199   This class contains all the informations related with a template class FIELD :
 200   - Components descriptions
 201   - Time step description
 202   - Location of the values (a SUPPORT class)
 203
 204 */
 205   class MEDMEM_EXPORT FIELD_ : public RCBASE    // GENERIC POINTER TO a template <class T, class INTERLACING_TAG> class FIELD
 206 {
 207 protected:
 208
 209   bool            _isRead ;
 210   bool            _isMinMax;
 211
 212   /*!
 213     \if developper
 214     Field name.
 215     \endif
 216   */
 217   string          _name ;
 218   /*!
 219     \if developper
 220     Field description.
 221     \endif
 222   */
 223   string          _description ;
 224   /*!
 225     \if developper
 226     Pointer to the support the field deals with.
 227     \endif
 228   */
 229   const SUPPORT * _support ;
 230
 231   /*!
 232     \if developper
 233     Number of field's components.
 234     \endif
 235   */
 236   int             _numberOfComponents ;
 237   /*!
 238     \if developper
 239     Number of field's values.
 240     doesn't take care of _numberOfComponents
 241     and number of Gauss points.
 242     \endif
 243   */
 244   int             _numberOfValues ;
 245
 246   /*!
 247     \if developper
 248     Array of size _numberOfComponents. \n
 249     (constant, scalar, vector, tensor)\n
 250     We could use an array of integer to store
 251     numbers of values: \n
 252     - 1 for scalar,\n
 253     - space dimension for vector,\n
 254     - space dimension square for tensor.\n
 255     So numbers of values per entities would be
 256     sum of _componentsTypes array.
 257
 258     Not implemented yet! All type are scalar !
 259     \endif
 260   */
 261   vector<int>     _componentsTypes ;
 262   /*!
 263     \if developper
 264     Array of size _numberOfComponents
 265     storing components names if any.
 266     \endif
 267   */
 268   vector<string>  _componentsNames;
 269   /*!
 270     \if developper
 271     Array of size _numberOfComponents
 272     storing components descriptions if any.
 273     \endif
 274   */
 275   vector<string>  _componentsDescriptions;
 276   /*!
 277     \if developper
 278     Array of size _numberOfComponents
 279     storing components units if any.
 280     \endif
 281   */
 282   vector<UNIT>    _componentsUnits;
 283   /*!
 284     \if developper
 285     Array of size _numberOfComponents
 286     storing components units if any.
 287     \endif
 288   */
 289   vector<string>  _MEDComponentsUnits;
 290   /*!
 291     \if developper
 292     Iteration number of the field.
 293     \endif
 294   */
 295   int             _iterationNumber ;
 296   /*!
 297     \if developper
 298     Time of the field.
 299     \endif
 300   */
 301   double          _time;
 302   /*!
 303     \if developper
 304     Order number of the field.
 305     \endif
 306   */
 307   int             _orderNumber ;
 308   /*!
 309     \if developper
 310     At the initialization step of the field using the constructors; this attribute,
 311     the value type (integer or real) , is set automatically. There is a get method
 312     but not a set method for this attribute.
 313     \endif
 314   */
 315   MED_EN::med_type_champ _valueType ;
 316   /*!
 317     \if developper
 318     At the initialization step of the field using the constructors; this attribute,
 319     the interlacing type (full interlace or no interlace field value storage), is set
 320     automatically. There is a get method but not a set method for this attribute.
 321     \endif
 322   */
 323   MED_EN::medModeSwitch _interlacingType;
 324
 325   vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
 326   static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
 327   static void _deepCheckFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
 328   void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL,
 329                                const bool           nodalAllowed = false) const  throw (MEDEXCEPTION);
 330   FIELD<double>* _getFieldSize(const SUPPORT *subSupport=NULL) const;
 331  public:
 332   /*!
 333     Constructor.
 334   */
 335   FIELD_ ();
 336   /*! \ifnot MEDMEM_ug
 337     Constructor.
 338 \endif
 339   */
 340   FIELD_(const SUPPORT * Support, const int NumberOfComponents);
 341   /*!  \ifnot MEDMEM_ug
 342     Copy constructor.
 343 \endif
 344   */
 345   FIELD_(const FIELD_ &m);
 346
 347   /*!
 348     Destructor.
 349   */
 350   virtual ~FIELD_();
 351
 352 public:
 353
 354   friend class MED_MED_RDONLY_DRIVER22;
 355   friend class MED_MED_WRONLY_DRIVER22;
 356   friend class MED_MED_RDWR_DRIVER22;
 357   friend class VTK_MED_DRIVER;
 358
 359  FIELD_& operator=(const FIELD_ &m);
 360
 361   virtual  void     rmDriver(int index=0);
 362
 363   /*! \if MEDMEM_ug
 364     \addtogroup FIELD_io
 365     @{
 366     \endif
 367   */
 368
 369   /*! Creates a driver for reading/writing fields in a file.
 370     \param driverType specifies the file type (MED_DRIVER, VTK_DRIVER)
 371     \param fileName name of the output file
 372     \param driverFieldName name of the field
 373     \param access specifies whether the file is opened for read, write or both.
 374   */
 375
 376   virtual   int     addDriver(driverTypes driverType,
 377                               const string & fileName="Default File Name.med",
 378                               const string & driverFieldName="Default Field Nam",
 379                               MED_EN::med_mode_acces access=MED_EN::RDWR) ;
 380
 381   virtual  int      addDriver( GENDRIVER & driver);
 382   virtual  void     read (driverTypes driverType, const std::string & fileName);
 383   virtual  void     read (const GENDRIVER &);
 384   virtual  void     read(int index=0);
 385   virtual  void     openAppend( void );
 386   virtual  void     write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
 387   virtual  void     write(driverTypes driverType,
 388                           const std::string & fileName,
 389                           MED_EN::med_mode_acces medMode=MED_EN::RDWR);
 390
 391   /*! Triggers the writing of the field with respect to the driver handle
 392     \a index given by \a addDriver(...) method. */
 393   virtual  void     write(int index=0);
 394   /*!\if MEDMEM_ug @} \endif */
 395
 396   virtual  void     writeAppend(const GENDRIVER &);
 397   virtual  void     writeAppend(int index=0, const string & driverName="");
 398
 399   inline void     setName(const string Name);
 400   inline string   getName() const;
 401   inline void     setDescription(const string Description);
 402   inline string   getDescription() const;
 403   inline const SUPPORT * getSupport() const;
 404   inline void     setSupport(const SUPPORT * support);
 405   inline void     setNumberOfComponents(const int NumberOfComponents);
 406   inline int      getNumberOfComponents() const;
 407   inline void     setNumberOfValues(const int NumberOfValues);
 408   inline int      getNumberOfValues() const;
 409   inline void     setComponentsNames(const string * ComponentsNames);
 410   inline void     setComponentName(int i, const string ComponentName);
 411   inline const string * getComponentsNames() const;
 412   inline string   getComponentName(int i) const;
 413   inline void     setComponentsDescriptions(const string * ComponentsDescriptions);
 414   inline void     setComponentDescription(int i, const string ComponentDescription);
 415   inline const string * getComponentsDescriptions() const;
 416   inline string   getComponentDescription(int i) const;
 417
 418   // provisoire : en attendant de regler le probleme des unites !
 419   inline void     setComponentsUnits(const UNIT * ComponentsUnits);
 420   inline const UNIT *   getComponentsUnits() const;
 421   inline const UNIT *   getComponentUnit(int i) const;
 422   inline void     setMEDComponentsUnits(const string * MEDComponentsUnits);
 423   inline void     setMEDComponentUnit(int i, const string MEDComponentUnit);
 424   inline const string * getMEDComponentsUnits() const;
 425   inline string   getMEDComponentUnit(int i) const;
 426
 427   inline void     setIterationNumber(int IterationNumber);
 428   inline int      getIterationNumber() const;
 429   inline void     setTime(double Time);
 430   inline double   getTime() const;
 431   inline void     setOrderNumber(int OrderNumber);
 432   inline int      getOrderNumber() const;
 433
 434   inline MED_EN::med_type_champ getValueType () const;
 435   inline MED_EN::medModeSwitch  getInterlacingType() const;
 436   virtual inline bool getGaussPresence() const throw (MEDEXCEPTION);
 437 protected:
 438   void copyGlobalInfo(const FIELD_& m);
 439 };
 440
 441 // -----------------
 442 // Methodes Inline
 443 // -----------------
 444 /*! \if MEDMEM_ug
 445 \addtogroup FIELD_getset
 446 @{
 447 \endif
 448 */
 449 /*!
 450   Sets FIELD name. The length should not exceed MED_TAILLE_NOM
 451 as defined in Med (i.e. 32 characters).
 452 */
 453 inline void FIELD_::setName(const string Name)
 454 {
 455   _name=Name;
 456 }
 457 /*!
 458   Gets FIELD name.
 459 */
 460 inline string FIELD_::getName() const
 461 {
 462   return _name;
 463 }
 464 /*!
 465   Sets FIELD description. The length should not exceed MED_TAILLE_DESC as defined in Med (i.e. 200 characters).
 466 */
 467 inline void FIELD_::setDescription(const string Description)
 468 {
 469   _description=Description;
 470 }
 471 /*!
 472   Gets FIELD description.
 473 */
 474 inline string FIELD_::getDescription() const
 475 {
 476   return _description;
 477 }
 478 /*!
 479   Sets FIELD number of components.
 480 */
 481 inline void FIELD_::setNumberOfComponents(const int NumberOfComponents)
 482 {
 483   _numberOfComponents=NumberOfComponents;
 484   _componentsTypes.resize(_numberOfComponents);
 485   _componentsNames.resize(_numberOfComponents);
 486   _componentsDescriptions.resize(_numberOfComponents);
 487   _componentsUnits.resize(_numberOfComponents);
 488   _MEDComponentsUnits.resize(_numberOfComponents);
 489 }
 490 /*!
 491   Gets FIELD number of components.
 492 */
 493 inline int FIELD_::getNumberOfComponents() const
 494 {
 495   return _numberOfComponents ;
 496 }
 497 /*!
 498   Sets FIELD number of values.
 499
 500   It must be the same than in the associated SUPPORT object.
 501 */
 502 inline void FIELD_::setNumberOfValues(const int NumberOfValues)
 503 {
 504   _numberOfValues=NumberOfValues;
 505 }
 506 /*!
 507   Gets FIELD number of value.
 508 */
 509 inline int FIELD_::getNumberOfValues() const
 510 {
 511   return _numberOfValues ;
 512 }
 513
 514 /*!
 515   Sets FIELD components names.
 516
 517   Duplicates the ComponentsNames string array to put components names in
 518   FIELD. ComponentsNames size must be equal to number of components.
 519 */
 520 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
 521 {
 522   _componentsNames.resize(_numberOfComponents);
 523   for (int i=0; i<_numberOfComponents; i++)
 524     _componentsNames[i]=ComponentsNames[i] ;
 525 }
 526 /*! \ifnot MEDMEM_ug
 527   Sets FIELD i^th component name.
 528
 529   i must be >=1 and <= number of components.
 530 \endif
 531 */
 532 inline void FIELD_::setComponentName(int i, const string ComponentName)
 533 {
 534   const char * LOC = " FIELD_::setComponentName() : ";
 535   BEGIN_OF_MED(LOC);
 536   if( i<1 || i>_numberOfComponents )
 537     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 538
 539   _componentsNames[i-1]=ComponentName ;
 540 }
 541 /*!
 542   Gets a reference to the string array which contain the components names.
 543
 544   This Array size is equal to number of components
 545 */
 546 inline const string * FIELD_::getComponentsNames() const
 547 {
 548   return &(_componentsNames[0]) ;
 549 }
 550 /*!\ifnot MEDMEM_ug
 551   Gets the name of the i^th component.
 552 \endif
 553 */
 554 inline string FIELD_::getComponentName(int i) const
 555 {
 556   const char * LOC = " FIELD_::getComponentName() : ";
 557   BEGIN_OF_MED(LOC);
 558   if( i<1 || i>_numberOfComponents )
 559     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 560
 561   return _componentsNames[i-1] ;
 562 }
 563 /*!
 564   Sets FIELD components descriptions.
 565
 566   Duplicates the ComponentsDescriptions string array to put components
 567   descriptions in FIELD.
 568   ComponentsDescriptions size must be equal to number of components.
 569 */
 570 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
 571 {
 572   _componentsDescriptions.resize(_numberOfComponents);
 573   for (int i=0; i<_numberOfComponents; i++)
 574     _componentsDescriptions[i]=ComponentsDescriptions[i] ;
 575 }
 576 /*!\ifnot MEDMEM_ug
 577   Sets FIELD i^th component description.
 578
 579   i must be >=1 and <= number of components.
 580 \endif
 581 */
 582 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
 583 {
 584   const char * LOC = " FIELD_::setComponentDescription() : ";
 585   BEGIN_OF_MED(LOC);
 586   if( i<1 || i>_numberOfComponents )
 587     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 588
 589   _componentsDescriptions[i-1]=ComponentDescription ;
 590 }
 591 /*!
 592   Gets a reference to the string array which contain the components descriptions.
 593
 594   This Array size is equal to number of components
 595 */
 596 inline const string * FIELD_::getComponentsDescriptions() const
 597 {
 598   return &(_componentsDescriptions[0]);
 599 }
 600 /*!\ifnot MEDMEM_ug
 601   Gets the description of the i^th component.
 602 \endif
 603 */
 604 inline string FIELD_::getComponentDescription(int i) const
 605 {
 606   const char * LOC = " FIELD_::setComponentDescription() : ";
 607   BEGIN_OF_MED(LOC);
 608   if( i<1 || i>_numberOfComponents )
 609     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 610
 611   return _componentsDescriptions[i-1];
 612 }
 613
 614 /*!\ifnot MEDMEM_ug
 615   \todo
 616   Sets FIELD components UNIT.
 617
 618   Duplicates the ComponentsUnits UNIT array to put components
 619   units in FIELD.
 620   ComponentsUnits size must be equal to number of components.
 621 \endif
 622 */
 623 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
 624 {
 625   _componentsUnits.resize(_numberOfComponents);
 626   for (int i=0; i<_numberOfComponents; i++)
 627     _componentsUnits[i]=ComponentsUnits[i] ;
 628 }
 629 /*!\ifnot MEDMEM_ug
 630   Gets a reference to the UNIT array which contain the components units.
 631
 632   This array size is equal to number of components
 633 \endif
 634 */
 635 inline const UNIT * FIELD_::getComponentsUnits() const
 636 {
 637   return &(_componentsUnits[0]);
 638 }
 639 /*!\ifnot MEDMEM_ug
 640   Gets the UNIT of the i^th component.
 641 \endif
 642 */
 643 inline const UNIT * FIELD_::getComponentUnit(int i) const
 644 {
 645   const char * LOC = " FIELD_::getComponentUnit() : ";
 646   BEGIN_OF_MED(LOC);
 647   if( i<1 || i>_numberOfComponents )
 648     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 649
 650   return &_componentsUnits[i-1] ;
 651 }
 652 /*!
 653   Sets FIELD components unit.
 654
 655   Duplicates the MEDComponentsUnits string array to put components
 656   units in FIELD.
 657   MEDComponentsUnits size must be equal to number of components.
 658
 659 */
 660 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
 661 {
 662   _MEDComponentsUnits.resize(_numberOfComponents);
 663   for (int i=0; i<_numberOfComponents; i++)
 664     _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
 665 }
 666 /*!\ifnot MEDMEM_ug
 667   Sets FIELD i^th component unit.
 668
 669   i must be >=1 and <= number of components.
 670 \endif
 671 */
 672 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
 673 {
 674   const char * LOC = " FIELD_::setMEDComponentUnit() : ";
 675   BEGIN_OF_MED(LOC);
 676   if( i<1 || i>_numberOfComponents )
 677     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 678
 679   _MEDComponentsUnits[i-1]=MEDComponentUnit ;
 680 }
 681 /*!
 682   Gets a reference to the string array which contain the components units.
 683
 684   This array size is equal to number of components
 685 */
 686 inline const string * FIELD_::getMEDComponentsUnits() const
 687 {
 688   return &(_MEDComponentsUnits[0]);
 689 }
 690 /*! \ifnot MEDMEM_ug
 691   Gets the string for unit of the i^th component.
 692 \endif
 693 */
 694 inline string FIELD_::getMEDComponentUnit(int i) const
 695 {
 696   const char * LOC = " FIELD_::getMEDComponentUnit() : ";
 697   BEGIN_OF_MED(LOC);
 698   if( i<1 || i>_numberOfComponents )
 699     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
 700
 701   return _MEDComponentsUnits[i-1] ;
 702 }
 703 /*!
 704   Sets the iteration number where FIELD has been calculated.
 705 */
 706 inline void FIELD_::setIterationNumber(int IterationNumber)
 707 {
 708   _iterationNumber=IterationNumber;
 709 }
 710 /*!
 711   Gets the iteration number where FIELD has been calculated.
 712 */
 713 inline int FIELD_::getIterationNumber() const
 714 {
 715   return _iterationNumber ;
 716 }
 717 /*!
 718   Sets the time when FIELD has been calculated.
 719 */
 720 inline void FIELD_::setTime(double Time)
 721 {
 722   _time=Time ;
 723 }
 724 /*!
 725   Gets the time when FIELD has been calculated.
 726 */
 727 inline double FIELD_::getTime() const
 728 {
 729   return _time ;
 730 }
 731 /*!
 732   Sets the order number where FIELD has been calculated.
 733
 734   It corresponds to internal iteration during one time step.
 735 */
 736 inline void FIELD_::setOrderNumber(int OrderNumber)
 737 {
 738   _orderNumber=OrderNumber ;
 739 }
 740 /*!
 741   Gets the order number where FIELD has been calculated.
 742 */
 743 inline int FIELD_::getOrderNumber() const
 744 {
 745   return _orderNumber ;
 746 }
 747 /*!
 748   Gets a reference to the SUPPORT object associated to FIELD.
 749 */
 750 inline  const SUPPORT * FIELD_::getSupport() const
 751 {
 752   return _support ;
 753 }
 754 /*!
 755   Sets the reference to the SUPPORT object associated to FIELD.
 756
 757   Reference is not duplicate, so it must not be deleted.
 758 */
 759 inline void FIELD_::setSupport(const SUPPORT * support)
 760 {
 761   //A.G. Addings for RC
 762   if(_support!=support)
 763     {
 764       if(_support)
 765         _support->removeReference();
 766       _support = support ;
 767       if(_support)
 768         _support->addReference();
 769     }
 770 }
 771 /*!
 772   Gets the FIELD med value type (MED_INT32 or MED_REEL64).
 773 */
 774 inline MED_EN::med_type_champ FIELD_::getValueType () const
 775 {
 776   return _valueType ;
 777 }
 778
 779 /*!
 780   Gets the FIELD med interlacing type (MED_FULL_INTERLACE or MED_NO_INTERLACE).
 781 */
 782   inline MED_EN::medModeSwitch FIELD_::getInterlacingType () const
 783 {
 784   return _interlacingType ;
 785 }
 786  /*!\if MEDMEM_ug @} \endif*/
 787
 788 /*!\if MEDMEM_ug
 789 \addtogroup FIELD_gauss
 790 @{
 791 \endif */
 792
 793 /*!
 794  Determines whether the field stores several Gauss points per element.
 795 */
 796   inline bool  FIELD_::getGaussPresence() const throw (MEDEXCEPTION)
 797 {
 798   const char * LOC = "FIELD_::getGaussPresence() : ";
 799   throw MEDEXCEPTION(STRING(LOC) << " This FIELD_ doesn't rely on a FIELD<T>" );
 800 }
 801
 802   /*!\if MEDMEM_ug @} \endif*/
 803
 804 } //End namespace MEDMEM
 805
 806 /////////////////////////
 807 // END OF CLASS FIELD_ //
 808 /////////////////////////
 809
 810 /*!
 811
 812   This template class contains informations related with a FIELD :
 813   - Values of the field, their type (real or integer), the storage mode (full interlace,
 814     no interlace or no interlace by type).
 815
 816 */
 817
 818
 819 namespace MEDMEM {
 820
 821   template<class T2> class MED_FIELD_RDONLY_DRIVER;
 822   template<class T2> class MED_FIELD_WRONLY_DRIVER;
 823   template<class T2> class VTK_FIELD_DRIVER;
 824
 825
 826   template <class T,
 827             class INTERLACING_TAG
 828             > class FIELD : public FIELD_
 829 {
 830 protected:
 831
 832   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array   ArrayNoGauss;
 833   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array     ArrayGauss;
 834   typedef typename MEDMEM_ArrayInterface<T,NoInterlace,NoGauss>::Array       ArrayNo;
 835   typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array     ArrayFull;
 836   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,NoGauss>::Array ArrayNoByType;
 837   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,Gauss>::Array   ArrayNoByTypeGauss;
 838   typedef MEDMEM_Array_ Array;
 839   typedef T ElementType;
 840   typedef INTERLACING_TAG InterlacingTag;
 841   typedef map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> locMap;
 842
 843   // array of value of type T
 844   Array *_value ;
 845
 846   // MESH, to be used for field reading from a file (if desired to link
 847   // to existing support instead of new support creation for the field)
 848   GMESH* _mesh;
 849
 850   // extrema values
 851   T _vmin;
 852   T _vmax;
 853
 854   map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> _gaussModel; //A changer quand les drivers seront template de l'entrelacement
 855
 856   static T _scalarForPow;
 857   static T pow(T x);
 858
 859 private:
 860   void _operation(const FIELD& m,const FIELD& n, const char* Op);
 861   void _operationInitialize(const FIELD& m,const FIELD& n, const char* Op);
 862   void _add_in_place(const FIELD& m,const FIELD& n);
 863   void _sub_in_place(const FIELD& m,const FIELD& n);
 864   void _mul_in_place(const FIELD& m,const FIELD& n);
 865   void _div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION);
 866 public:
 867   FIELD();
 868   FIELD(const FIELD &m);
 869   FIELD(const SUPPORT * Support, const int NumberOfComponents) throw (MEDEXCEPTION);
 870   FIELD(driverTypes driverType,
 871         const string & fileName, const string & fieldDriverName,
 872         const int iterationNumber=-1, const int orderNumber=-1,
 873         GMESH* mesh = 0)
 874     throw (MEDEXCEPTION);
 875   FIELD(const SUPPORT * Support, driverTypes driverType,
 876         const string & fileName="", const string & fieldName="",
 877         const int iterationNumber = -1, const int orderNumber = -1)
 878     throw (MEDEXCEPTION);
 879   ~FIELD();
 880
 881 public:
 882   FIELD & operator=(const FIELD &m);
 883         FIELD & operator=(T value);
 884   FIELD *operator+(const FIELD& m) const;
 885   FIELD *operator-(const FIELD& m) const;
 886   FIELD *operator*(const FIELD& m) const;
 887   FIELD *operator/(const FIELD& m) const;
 888   FIELD *operator-() const;
 889   FIELD& operator+=(const FIELD& m);
 890   FIELD& operator-=(const FIELD& m);
 891   FIELD& operator*=(const FIELD& m);
 892   FIELD& operator/=(const FIELD& m);
 893
 894   void          applyLin(T a, T b, int icomp);
 895   static FIELD* add(const FIELD& m, const FIELD& n);
 896   static FIELD* addDeep(const FIELD& m, const FIELD& n);
 897   static FIELD* sub(const FIELD& m, const FIELD& n);
 898   static FIELD* subDeep(const FIELD& m, const FIELD& n);
 899   static FIELD* mul(const FIELD& m, const FIELD& n);
 900   static FIELD* mulDeep(const FIELD& m, const FIELD& n);
 901   static FIELD* div(const FIELD& m, const FIELD& n);
 902   static FIELD* divDeep(const FIELD& m, const FIELD& n);
 903   double normMax() const throw (MEDEXCEPTION);
 904
 905   //------- TDG and BS addings
 906
 907   void getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION);
 908   vector<int> getHistogram(int &nbint) throw (MEDEXCEPTION);
 909   FIELD<double>* buildGradient() const throw (MEDEXCEPTION);
 910   FIELD<double>* buildNorm2Field() const throw (MEDEXCEPTION);
 911
 912   //-------------------
 913
 914   double norm2() const throw (MEDEXCEPTION);
 915   void   applyLin(T a, T b);
 916   template <T T_function(T)> void applyFunc();
 917   void applyPow(T scalar);
 918   static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
 919   double normL2(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 920   double normL2(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 921   double normL1(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 922   double normL1(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 923   double integral(const SUPPORT *subSupport=NULL) const throw (MEDEXCEPTION);
 924   FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
 925
 926   friend class MED_FIELD_RDONLY_DRIVER<T>;
 927   friend class MED_FIELD_WRONLY_DRIVER<T>;
 928   friend class VTK_FIELD_DRIVER<T>;
 929
 930   void init ();
 931   void rmDriver(int index=0);
 932   int  addDriver(driverTypes driverType,
 933                  const string & fileName="Default File Name.med",
 934                  const string & driverFieldName="Default Field Name",
 935                  MED_EN::med_mode_acces access=MED_EN::RDWR) ;
 936
 937   int  addDriver(GENDRIVER & driver);
 938
 939   void allocValue(const int NumberOfComponents);
 940   void allocValue(const int NumberOfComponents, const int LengthValue);
 941
 942   void deallocValue();
 943
 944   inline void read(int index=0);
 945   inline void read(const GENDRIVER & genDriver);
 946   inline void read(driverTypes driverType, const std::string& filename);
 947   inline void write(int index=0);
 948   inline void write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
 949   inline void write(driverTypes        driverType,
 950                     const std::string& filename,
 951                     MED_EN::med_mode_acces medMode=MED_EN::RDWR);
 952
 953   inline void writeAppend(int index=0, const string & driverName = "");
 954   inline void writeAppend(const GENDRIVER &);
 955
 956   inline MEDMEM_Array_  * getArray()        const throw (MEDEXCEPTION);
 957   inline ArrayGauss     * getArrayGauss()   const throw (MEDEXCEPTION);
 958   inline ArrayNoGauss   * getArrayNoGauss() const throw (MEDEXCEPTION);
 959   inline bool            getGaussPresence() const throw (MEDEXCEPTION);
 960
 961   inline int          getValueLength() const throw (MEDEXCEPTION);
 962
 963   /*! \if MEDMEM_ug
 964 \addtogroup FIELD_value
 965 @{
 966 \endif  */
 967   /*! Returns a pointer to the value array.*/
 968   inline const T*     getValue()       const throw (MEDEXCEPTION);
 969   inline const T*     getRow(int i)    const throw (MEDEXCEPTION);
 970   inline const T*     getColumn(int j) const throw (MEDEXCEPTION);
 971 /*!
 972   Returns the value of \f$ i^{th} \f$ element and \f$ j^{th}\f$ component.
 973   This method only works with fields having no particular Gauss point
 974 definition (i.e., fields having one value per element).
 975  This method makes the retrieval of the value independent from the
 976   interlacing pattern, but it is slower than the complete retrieval
 977   obtained by the \b getValue() method.
 978 */
 979
 980   inline T            getValueIJ(int i,int j) const throw (MEDEXCEPTION);
 981
 982 /*!
 983   Returns the \f$ j^{th}\f$  component of \f$ k^{th}\f$  Gauss points of \f$ i^{th}\f$  value.
 984   This method is compatible with elements having more than one Gauss point.
 985   This method makes the retrieval of the value independent from the
 986   interlacing pattern, but it is slower than the complete retrieval
 987   obtained by the \b getValue() method.
 988 */
 989   inline T            getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION);
 990
 991   inline int          getValueByTypeLength(int t)                const throw (MEDEXCEPTION);
 992   inline const T*     getValueByType(int t)                      const throw (MEDEXCEPTION);
 993   inline T            getValueIJByType(int i,int j,int t)        const throw (MEDEXCEPTION);
 994   inline T            getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION);
 995
 996   /*!
 997                 The following example describes the creation of a FIELD.
 998
 999                 \example FIELDcreate.cxx
1000
1001  \if MEDMEM_ug @} \endif */
1002
1003   bool                getValueOnElement(int eltIdInSup,T* retValues) const throw (MEDEXCEPTION);
1004   void                getValueOnPoint(const double* coords, double* output) const throw (MEDEXCEPTION);
1005   void                getValueOnPoints(int nb_points, const double* coords, double* output) const throw (MEDEXCEPTION);
1006
1007   const int   getNumberOfGeometricTypes() const throw (MEDEXCEPTION);
1008   const GAUSS_LOCALIZATION<INTERLACING_TAG> & getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1009   const GAUSS_LOCALIZATION<INTERLACING_TAG> * getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1010   const GAUSS_LOCALIZATION_* getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1011   void setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG> & gaussloc);
1012   void setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc);
1013   const int * getNumberOfGaussPoints() const throw (MEDEXCEPTION);
1014   const int   getNumberOfGaussPoints( MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1015   const int   getNbGaussI(int i)          const throw (MEDEXCEPTION);
1016   const int * getNumberOfElements()       const throw (MEDEXCEPTION);
1017   const MED_EN::medGeometryElement  * getGeometricTypes()  const throw (MEDEXCEPTION);
1018   bool        isOnAllElements()           const throw (MEDEXCEPTION);
1019
1020   inline void setArray(MEDMEM_Array_ *value) throw (MEDEXCEPTION);
1021
1022   FIELD<double, FullInterlace>* getGaussPointsCoordinates() const throw (MEDEXCEPTION);
1023
1024   /*! \if MEDMEM_ug
1025  \addtogroup FIELD_value
1026 @{
1027 \endif
1028   */
1029 /*!
1030 This method makes it possible to have the field pointing to
1031 an existing value array. The ordering of the elements in the value array must
1032 conform to the MEDMEM ordering (I,K,J) : the outer loop is on the elements,
1033 the intermediate loop is on the Gauss points, the inner loop is on
1034 the components.
1035 */
1036   inline void setValue( T* value) throw (MEDEXCEPTION);
1037   inline void setRow( int i, T* value) throw (MEDEXCEPTION);
1038   inline void setColumn( int i, T* value) throw (MEDEXCEPTION);
1039 /*!
1040   Sets the value of \f$ i^{th} \f$ element and \f$ j^{th}\f$ component with \a value.
1041 */
1042   inline void setValueIJ(int i, int j, T value) throw (MEDEXCEPTION);
1043   /*! \if MEDMEM_ug @} \endif */
1044   inline void setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION);
1045   inline void setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION);
1046   inline void setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION);
1047
1048   typedef void (*myFuncType)(const double *,T*);
1049   void fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION);
1050   typedef void (*myFuncType2)(const T *,T*);
1051   FIELD<T,INTERLACING_TAG> *execFunc(int nbOfComponents, myFuncType2 f) throw (MEDEXCEPTION);
1052 };
1053 }
1054
1055 #include "MEDMEM_DriverFactory.hxx"
1056
1057 namespace MEDMEM {
1058
1059 template <class T,class INTERLACING_TAG> T FIELD<T, INTERLACING_TAG>::_scalarForPow=1;
1060
1061 // --------------------
1062 // Implemented Methods
1063 // --------------------
1064
1065 /*!
1066   Constructor with no parameter, most of the attribut members are set to NULL.
1067 */
1068 template <class T, class INTERLACING_TAG>
1069 FIELD<T, INTERLACING_TAG>::FIELD():FIELD_()
1070 {
1071   MESSAGE_MED("Constructeur FIELD sans parametre");
1072
1073   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
1074   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE);
1075   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
1076
1077   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
1078   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE);
1079   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
1080
1081   _value = ( ArrayNoGauss * ) NULL;
1082
1083   _mesh  = ( MESH* ) NULL;
1084 }
1085
1086         /*!
1087 \addtogroup FIELD_constructors FIELD<T> constructors
1088 @{
1089         */
1090
1091 /*!
1092   Constructor that allocates the value array with the dimensions provided by
1093 \a NumberOfComponents and the dimension of \a Support. The value array is
1094  allocated but not initialized.
1095 This constructor does not allow the creation of fields with Gauss points.
1096 \param Support support on which the field lies
1097 \param NumberOfComponents number of components of the variable stored. For instance,
1098 it will be 3 for a (vx,vy,vz) vector.
1099
1100 \code
1101 FIELD<double> field (support, 3);
1102 int nbelem = support->getNumberOfElements(MED_ALL_ELEMENTS);
1103 for (int i=1; i<=nbelem; i++)
1104    for (j=1; j<=3;j++)
1105        field->setValueIJ(i,j,0.0);
1106 \endcode
1107 */
1108 template <class T, class INTERLACING_TAG>
1109 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
1110                                  const int NumberOfComponents) throw (MEDEXCEPTION) :
1111   FIELD_(Support, NumberOfComponents),_value(NULL)
1112 {
1113   const char* LOC = "FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)";
1114   BEGIN_OF_MED(LOC);
1115   SCRUTE_MED(this);
1116
1117   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
1118   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
1119     FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
1120
1121   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
1122   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
1123     FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
1124
1125   try
1126     {
1127       // becarefull about the numbre of gauss point
1128       _numberOfValues = Support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
1129     }
1130 #if defined(_DEBUG_) || defined(_DEBUG)
1131   catch (MEDEXCEPTION &ex)
1132 #else
1133   catch (MEDEXCEPTION )
1134 #endif
1135     {
1136       MESSAGE_MED("No value defined ! ("<<ex.what()<<")");
1137     }
1138   MESSAGE_MED("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
1139   if ( _numberOfValues > 0 )
1140     {
1141       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
1142         {
1143           const int * nbelgeo = Support->getNumberOfElements();
1144           vector<int> nbelgeoc( Support->getNumberOfTypes() + 1 );
1145           nbelgeoc[0] = 0;
1146           for ( int t = 1; t < (int)nbelgeoc.size(); ++t )
1147             nbelgeoc[t] = nbelgeoc[t-1] + nbelgeo[t-1];
1148           _value = new ArrayNoByType (_numberOfComponents,_numberOfValues,
1149                                       Support->getNumberOfTypes(), &nbelgeoc[0]);
1150         }
1151       else
1152         {
1153           _value = new ArrayNoGauss (_numberOfComponents,_numberOfValues);
1154         }
1155       _isRead = true ;
1156     }
1157   _mesh  = ( MESH* ) NULL;
1158
1159   END_OF_MED(LOC);
1160 }
1161         /*!
1162 @}
1163         */
1164 /*!
1165   \if developper
1166   \endif
1167 */
1168 template <class T, class INTERLACING_TAG> void FIELD<T, INTERLACING_TAG>::init ()
1169 {
1170 }
1171
1172 /*!
1173   Copy constructor.
1174 */
1175 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::FIELD(const FIELD & m):
1176   FIELD_(m)
1177 {
1178   MESSAGE_MED("Constructeur FIELD de recopie");
1179
1180   // RECOPIE PROFONDE <> de l'operateur= Rmq from EF
1181   if (m._value != NULL)
1182     {
1183       if ( m.getGaussPresence() )
1184         _value = new ArrayGauss( *(static_cast< ArrayGauss * > (m._value) ) ,false);
1185       else
1186         _value = new ArrayNoGauss( *(static_cast< ArrayNoGauss * > (m._value)) ,false);
1187     }
1188   else
1189     _value = (ArrayNoGauss *) NULL;
1190   locMap::const_iterator it;
1191   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
1192     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
1193       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
1194                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
1195                                               );
1196
1197   _valueType       = m._valueType;
1198   _interlacingType = m._interlacingType;
1199   //drivers = m._drivers;
1200   _mesh            = m._mesh;
1201   if(_mesh)
1202     _mesh->addReference();
1203 }
1204
1205 /*!
1206   \if developper
1207   Not implemented.
1208   \endif
1209 */
1210 template <class T, class INTERLACING_TAG>
1211 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(const FIELD &m)
1212 {
1213   MESSAGE_MED("Appel de FIELD<T>::operator=") ;
1214   if ( this == &m) return *this;
1215
1216   // copy values array
1217   FIELD_::operator=(m);  // Driver are ignored & ?copie su pointeur de Support?
1218
1219   _value           = m._value; //PROBLEME RECOPIE DES POINTEURS PAS COHERENT AVEC LE
1220                                //CONSTRUCTEUR PAR RECOPIE
1221                                //CF :Commentaire dans MEDMEM_Array
1222   locMap::const_iterator it;
1223   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
1224     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
1225       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
1226                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
1227                                               );
1228
1229   _valueType       = m._valueType;
1230   _interlacingType = m._interlacingType;
1231   if(_mesh!=m._mesh)
1232     {
1233       if(_mesh)
1234         _mesh->removeReference();
1235       _mesh = m._mesh;
1236       if(_mesh)
1237         _mesh->addReference();
1238     }
1239   return *this;
1240 }
1241
1242 /*!
1243         Initializes all the field values to \a value
1244 */
1245 template <class T, class INTERLACING_TAG>
1246 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(T value)
1247 {
1248   MESSAGE_MED("Appel de FIELD<T>::operator= T") ;
1249         int size=getNumberOfComponents()*getNumberOfValues();
1250         T* ptr= const_cast<T*>( getValue());
1251         for (int i=0; i< size; i++)
1252                 {*ptr++=value;}
1253
1254   return *this;
1255 }
1256
1257 /*!\addtogroup FIELD_algo
1258 @{
1259 */
1260
1261 /*!
1262      Overload addition operator.
1263      This operation is authorized only for compatible fields that have the same support.
1264      The compatibility checking includes equality tests of the folowing data members:\n
1265      - _support
1266      - _numberOfComponents
1267      - _numberOfValues
1268      - _componentsTypes
1269      - _MEDComponentsUnits.
1270
1271      The data members of the returned field are initialized, based on the first field, except for the name,
1272      which is the combination of the two field's names and the operator.
1273      Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> \n
1274      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1275      When using python, this operator calls the copy constructor in any case.
1276      The user has to be aware that when using operator + in associatives expressions like
1277      <tt> a = b + c + d +e; </tt> \n
1278      no optimisation is performed : the evaluation of last expression requires the construction of
1279      3 temporary fields.
1280 */
1281 template <class T, class INTERLACING_TAG>
1282 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator+(const FIELD & m) const
1283 {
1284   const char* LOC = "FIELD<T>::operator+(const FIELD & m)";
1285   BEGIN_OF_MED(LOC);
1286     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1287
1288     // Creation of the result - memory is allocated by FIELD constructor
1289     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1290     result->_operationInitialize(*this,m,"+"); // perform Atribute's initialization
1291     result->_add_in_place(*this,m); // perform addition
1292
1293   END_OF_MED(LOC);
1294     return result;
1295 }
1296
1297 /*!  Overloaded Operator +=
1298  *   Operations are directly performed in the first field's array.
1299  *   This operation is authorized only for compatible fields that have the same support.
1300  */
1301 template <class T, class INTERLACING_TAG>
1302 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator+=(const FIELD & m)
1303 {
1304   const char* LOC = "FIELD<T>::operator+=(const FIELD & m)";
1305   BEGIN_OF_MED(LOC);
1306     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1307
1308     const T* value1=m.getValue(); // get pointers to the values we are adding
1309
1310     // get a non const pointer to the inside array of values and perform operation
1311     T * value=const_cast<T *> (getValue());
1312     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1313     const T* endV=value+size; // pointer to the end of value
1314     for(;value!=endV; value1++,value++)
1315         *value += *value1;
1316   END_OF_MED(LOC);
1317     return *this;
1318 }
1319
1320
1321 /*! Addition of fields. Static member function.
1322  *  The function return a pointer to a new created field that holds the addition.
1323  *  Data members are checked for compatibility and initialized.
1324  *  The user is in charge of memory deallocation.
1325  */
1326 template <class T, class INTERLACING_TAG>
1327 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::add(const FIELD& m, const FIELD& n)
1328 {
1329   const char* LOC = "FIELD<T>::add(const FIELD & m, const FIELD& n)";
1330   BEGIN_OF_MED(LOC);
1331     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1332
1333     // Creation of a new field
1334     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1335                                                                       m.getNumberOfComponents());
1336     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1337     result->_add_in_place(m,n); // perform addition
1338
1339   END_OF_MED(LOC);
1340     return result;
1341 }
1342
1343 /*! Same as add method except that field check is deeper.
1344  */
1345 template <class T, class INTERLACING_TAG>
1346 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::addDeep(const FIELD& m, const FIELD& n)
1347 {
1348   const char* LOC = "FIELD<T>::addDeep(const FIELD & m, const FIELD& n)";
1349   BEGIN_OF_MED(LOC);
1350     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1351
1352     // Creation of a new field
1353     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1354                                                                       m.getNumberOfComponents());
1355     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1356     result->_add_in_place(m,n); // perform addition
1357
1358   END_OF_MED(LOC);
1359     return result;
1360 }
1361
1362 /*!
1363      Overload substraction operator.
1364      This operation is authorized only for compatible fields that have the same support.
1365      The compatibility checking includes equality tests of the folowing data members:\n
1366      - _support
1367      - _numberOfComponents
1368      - _numberOfValues
1369      - _componentsTypes
1370      - _MEDComponentsUnits.
1371
1372      The data members of the returned field are initialized, based on the first field, except for the name,
1373      which is the combination of the two field's names and the operator.
1374      Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> \n
1375      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1376      When using python, this operator calls the copy constructor in any case.
1377      The user has to be aware that when using operator - in associatives expressions like
1378      <tt> a = b - c - d -e; </tt> \n
1379      no optimisation is performed : the evaluation of last expression requires the construction of
1380      3 temporary fields.
1381 */
1382 template <class T, class INTERLACING_TAG>
1383 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator-(const FIELD & m) const
1384 {
1385   const char* LOC = "FIELD<T>::operator-(const FIELD & m)";
1386   BEGIN_OF_MED(LOC);
1387     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1388
1389     // Creation of the result - memory is allocated by FIELD constructor
1390     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1391     result->_operationInitialize(*this,m,"-"); // perform Atribute's initialization
1392     result->_sub_in_place(*this,m); // perform substracion
1393
1394   END_OF_MED(LOC);
1395     return result;
1396 }
1397
1398 template <class T, class INTERLACING_TAG>
1399 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator-() const
1400 {
1401   const char* LOC = "FIELD<T>::operator-()";
1402   BEGIN_OF_MED(LOC);
1403
1404     // Creation of the result - memory is allocated by FIELD constructor
1405   FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1406     // Atribute's initialization
1407     result->setName("- "+getName());
1408     result->setComponentsNames(getComponentsNames());
1409     // not yet implemented    setComponentType(getComponentType());
1410     result->setComponentsDescriptions(getComponentsDescriptions());
1411     result->setMEDComponentsUnits(getMEDComponentsUnits());
1412     result->setComponentsUnits(getComponentsUnits());
1413     result->setIterationNumber(getIterationNumber());
1414     result->setTime(getTime());
1415     result->setOrderNumber(getOrderNumber());
1416
1417     const T* value1=getValue();
1418     // get a non const pointer to the inside array of values and perform operation
1419     T * value=const_cast<T *> (result->getValue());
1420     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1421     const T* endV=value+size; // pointer to the end of value
1422
1423     for(;value!=endV; value1++,value++)
1424         *value = -(*value1);
1425   END_OF_MED(LOC);
1426     return result;
1427 }
1428
1429 /*!  Overloaded Operator -=
1430  *   Operations are directly performed in the first field's array.
1431  *   This operation is authorized only for compatible fields that have the same support.
1432  */
1433 template <class T, class INTERLACING_TAG>
1434 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator-=(const FIELD & m)
1435 {
1436   const char* LOC = "FIELD<T>::operator-=(const FIELD & m)";
1437   BEGIN_OF_MED(LOC);
1438     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1439
1440     const T* value1=m.getValue();
1441
1442     // get a non const pointer to the inside array of values and perform operation
1443     T * value=const_cast<T *> (getValue());
1444     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1445     const T* endV=value+size; // pointer to the end of value
1446
1447     for(;value!=endV; value1++,value++)
1448         *value -= *value1;
1449
1450   END_OF_MED(LOC);
1451     return *this;
1452 }
1453
1454
1455
1456 /*!  Apply to a given field component the linear function x -> ax+b.
1457  *   calculation is done "in place".
1458  */
1459 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b, int icomp)
1460 {
1461     // get a non const pointer to the inside array of values and perform operation in place
1462     T * value=const_cast<T *> (getValue());
1463
1464     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1465
1466     if (size>0) // for a negative size, there is nothing to do
1467     {
1468                         value+=icomp-1;
1469                         const T* lastvalue=value+size; // pointer to the end of value
1470
1471                         for(;value!=lastvalue; value+=getNumberOfComponents()) // apply linear transformation
1472                                 *value = a*(*value)+b;
1473     }
1474 }
1475
1476 /*! Substraction of fields. Static member function.
1477  *  The function return a pointer to a new created field that holds the substraction.
1478  *  Data members are checked for compatibility and initialized.
1479  *  The user is in charge of memory deallocation.
1480  */
1481 template <class T, class INTERLACING_TAG>
1482 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::sub(const FIELD& m, const FIELD& n)
1483 {
1484   const char* LOC = "FIELD<T>::sub(const FIELD & m, const FIELD& n)";
1485   BEGIN_OF_MED(LOC);
1486     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1487
1488     // Creation of a new field
1489     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1490                                                                       m.getNumberOfComponents());
1491     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1492     result->_sub_in_place(m,n); // perform substraction
1493
1494   END_OF_MED(LOC);
1495     return result;
1496 }
1497
1498 /*! Same as sub method except that field check is deeper.
1499  */
1500 template <class T, class INTERLACING_TAG>
1501 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::subDeep(const FIELD& m, const FIELD& n)
1502 {
1503   const char* LOC = "FIELD<T>::subDeep(const FIELD & m, const FIELD& n)";
1504   BEGIN_OF_MED(LOC);
1505     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1506
1507     // Creation of a new field
1508     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1509                                                                       m.getNumberOfComponents());
1510     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1511     result->_sub_in_place(m,n); // perform substraction
1512
1513   END_OF_MED(LOC);
1514     return result;
1515 }
1516
1517 /*!
1518      Overload multiplication operator.
1519      This operation is authorized only for compatible fields that have the same support.
1520      The compatibility checking includes equality tests of the folowing data members:\n
1521      - _support
1522      - _numberOfComponents
1523      - _numberOfValues
1524      - _componentsTypes
1525      - _MEDComponentsUnits.
1526
1527      The data members of the returned field are initialized, based on the first field, except for the name,
1528      which is the combination of the two field's names and the operator.
1529      Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> \n
1530      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1531      When using python, this operator calls the copy constructor in any case.
1532      The user has to be aware that when using operator * in associatives expressions like
1533      <tt> a = b * c * d *e; </tt> \n
1534      no optimisation is performed : the evaluation of last expression requires the construction of
1535      3 temporary fields.
1536 */
1537 template <class T, class INTERLACING_TAG>
1538 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator*(const FIELD & m) const
1539 {
1540   const char* LOC = "FIELD<T>::operator*(const FIELD & m)";
1541   BEGIN_OF_MED(LOC);
1542     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1543
1544     // Creation of the result - memory is allocated by FIELD constructor
1545     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1546     result->_operationInitialize(*this,m,"*"); // perform Atribute's initialization
1547     result->_mul_in_place(*this,m); // perform multiplication
1548
1549   END_OF_MED(LOC);
1550     return result;
1551 }
1552
1553 /*!  Overloaded Operator *=
1554  *   Operations are directly performed in the first field's array.
1555  *   This operation is authorized only for compatible fields that have the same support.
1556  */
1557 template <class T, class INTERLACING_TAG>
1558 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator*=(const FIELD & m)
1559 {
1560   const char* LOC = "FIELD<T>::operator*=(const FIELD & m)";
1561   BEGIN_OF_MED(LOC);
1562     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1563
1564     const T* value1=m.getValue();
1565
1566     // get a non const pointer to the inside array of values and perform operation
1567     T * value=const_cast<T *> (getValue());
1568     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1569     const T* endV=value+size; // pointer to the end of value
1570
1571     for(;value!=endV; value1++,value++)
1572         *value *= *value1;
1573
1574   END_OF_MED(LOC);
1575     return *this;
1576 }
1577
1578
1579 /*! Multiplication of fields. Static member function.
1580  *  The function return a pointer to a new created field that holds the multiplication.
1581  *  Data members are checked for compatibility and initialized.
1582  *  The user is in charge of memory deallocation.
1583  */
1584 template <class T, class INTERLACING_TAG>
1585 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mul(const FIELD& m, const FIELD& n)
1586 {
1587   const char* LOC = "FIELD<T>::mul(const FIELD & m, const FIELD& n)";
1588   BEGIN_OF_MED(LOC);
1589     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1590
1591     // Creation of a new field
1592     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1593                                                                       m.getNumberOfComponents());
1594     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1595     result->_mul_in_place(m,n); // perform multiplication
1596
1597   END_OF_MED(LOC);
1598     return result;
1599 }
1600
1601 /*! Same as mul method except that field check is deeper.
1602  */
1603 template <class T, class INTERLACING_TAG>
1604 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mulDeep(const FIELD& m, const FIELD& n)
1605 {
1606   const char* LOC = "FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)";
1607   BEGIN_OF_MED(LOC);
1608     FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1609
1610     // Creation of a new field
1611     FIELD<T, INTERLACING_TAG>* result = new FIELD<T,INTERLACING_TAG>(m.getSupport(),
1612                                                                      m.getNumberOfComponents());
1613     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1614     result->_mul_in_place(m,n); // perform multiplication
1615
1616   END_OF_MED(LOC);
1617     return result;
1618 }
1619
1620 /*!
1621      Overload division operator.
1622      This operation is authorized only for compatible fields that have the same support.
1623      The compatibility checking includes equality tests of the folowing data members:\n
1624      - _support
1625      - _numberOfComponents
1626      - _numberOfValues
1627      - _componentsTypes
1628      - _MEDComponentsUnits.
1629
1630      The data members of the returned field are initialized, based on the first field, except for the name,
1631      which is the combination of the two field's names and the operator.
1632      Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> \n
1633      In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1634      When using python, this operator calls the copy constructor in any case.
1635      The user has to be aware that when using operator / in associatives expressions like
1636      <tt> a = b / c / d /e; </tt> \n
1637      no optimisation is performed : the evaluation of last expression requires the construction of
1638      3 temporary fields.
1639 */
1640 template <class T, class INTERLACING_TAG>
1641 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator/(const FIELD & m) const
1642 {
1643   const char* LOC = "FIELD<T>::operator/(const FIELD & m)";
1644   BEGIN_OF_MED(LOC);
1645     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1646
1647     // Creation of the result - memory is allocated by FIELD constructor
1648     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1649     try
1650       {
1651         result->_operationInitialize(*this,m,"/"); // perform Atribute's initialization
1652         result->_div_in_place(*this,m); // perform division
1653       }
1654     catch(MEDEXCEPTION& e)
1655       {
1656         result->removeReference();
1657         throw e;
1658       }
1659
1660   END_OF_MED(LOC);
1661     return result;
1662 }
1663
1664
1665 /*!  Overloaded Operator /=
1666  *   Operations are directly performed in the first field's array.
1667  *   This operation is authorized only for compatible fields that have the same support.
1668  */
1669 template <class T, class INTERLACING_TAG>
1670 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator/=(const FIELD & m)
1671 {
1672   const char* LOC = "FIELD<T>::operator/=(const FIELD & m)";
1673   BEGIN_OF_MED(LOC);
1674     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1675
1676     const T* value1=m.getValue(); // get pointers to the values we are adding
1677
1678     // get a non const pointer to the inside array of values and perform operation
1679     T * value=const_cast<T *> (getValue());
1680     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1681     const T* endV=value+size; // pointer to the end of value
1682
1683     for(;value!=endV; value1++,value++)
1684         *value /= *value1;
1685
1686   END_OF_MED(LOC);
1687     return *this;
1688 }
1689
1690
1691 /*! Division of fields. Static member function.
1692  *  The function return a pointer to a new created field that holds the division.
1693  *  Data members are checked for compatibility and initialized.
1694  *  The user is in charge of memory deallocation.
1695  */
1696 template <class T, class INTERLACING_TAG>
1697 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::div(const FIELD& m, const FIELD& n)
1698 {
1699   const char* LOC = "FIELD<T>::div(const FIELD & m, const FIELD& n)";
1700   BEGIN_OF_MED(LOC);
1701     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1702
1703     // Creation of a new field
1704     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1705                                                                       m.getNumberOfComponents());
1706     try
1707       {
1708         result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1709         result->_div_in_place(m,n); // perform division
1710       }
1711     catch(MEDEXCEPTION& e)
1712       {
1713         result->removeReference();
1714         throw e;
1715       }
1716   END_OF_MED(LOC);
1717     return result;
1718 }
1719
1720 /*! Same as div method except that field check is deeper.
1721  */
1722 template <class T,class INTERLACING_TAG>
1723 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::divDeep(const FIELD& m, const FIELD& n)
1724 {
1725   const char* LOC = "FIELD<T>::divDeep(const FIELD & m, const FIELD& n)";
1726   BEGIN_OF_MED(LOC);
1727   FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1728
1729   // Creation of a new field
1730   FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1731                                                                     m.getNumberOfComponents());
1732   try
1733     {
1734       result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1735       result->_div_in_place(m,n); // perform division
1736     }
1737   catch(MEDEXCEPTION& e)
1738       {
1739         result->removeReference();
1740         throw e;
1741       }
1742   END_OF_MED(LOC);
1743   return result;
1744 }
1745
1746
1747 /*!
1748 @}
1749 */
1750
1751 /*!
1752   \if developper
1753   This internal method initialize the members of a new field created to hold the result of the operation Op .
1754   Initialization is based on the first field, except for the name, which is the combination of the two field's names
1755   and the operator.
1756   \endif
1757 */
1758 template <class T, class INTERLACING_TAG>
1759 void FIELD<T, INTERLACING_TAG>::_operationInitialize(const FIELD& m,const FIELD& n, const char* Op)
1760 {
1761     MESSAGE_MED("Appel methode interne " << Op);
1762
1763     // Atribute's initialization (copy of the first field's attributes)
1764     // Other data members (_support, _numberOfValues) are initialized in the field's constr.
1765     setName(m.getName()+" "+Op+" "+n.getName());
1766     setComponentsNames(m.getComponentsNames());
1767     setComponentsDescriptions(m.getComponentsDescriptions());
1768     setMEDComponentsUnits(m.getMEDComponentsUnits());
1769
1770     // The following data member may differ from field m to n.
1771     // The initialization is done based on the first field.
1772
1773     setComponentsUnits(m.getComponentsUnits());
1774
1775     setIterationNumber(m.getIterationNumber());
1776     setTime(m.getTime());
1777     setOrderNumber(m.getOrderNumber());
1778 }
1779
1780
1781 /*!
1782   \if developper
1783   Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
1784   This method is applied to a just created field with medModeSwitch mode.
1785   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1786   it doesn't exist!
1787   \endif
1788 */
1789 template <class T, class INTERLACING_TAG>
1790 void FIELD<T, INTERLACING_TAG>::_add_in_place(const FIELD& m,const FIELD& n)
1791 {
1792     // get pointers to the values we are adding
1793     const T* value1=m.getValue();
1794     const T* value2=n.getValue();
1795     // get a non const pointer to the inside array of values and perform operation
1796     T * value=const_cast<T *> (getValue());
1797
1798     const int size=getNumberOfValues()*getNumberOfComponents();
1799     SCRUTE_MED(size);
1800     const T* endV1=value1+size;
1801     for(;value1!=endV1; value1++,value2++,value++)
1802         *value=(*value1)+(*value2);
1803 }
1804
1805 /*!
1806   \if developper
1807   Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
1808   This method is applied to a just created field with medModeSwitch mode.
1809   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1810   it doesn't exist!
1811   \endif
1812 */
1813 template <class T, class INTERLACING_TAG>
1814 void FIELD<T, INTERLACING_TAG>::_sub_in_place(const FIELD& m,const FIELD& n)
1815 {
1816     // get pointers to the values we are adding
1817     const T* value1=m.getValue();
1818     const T* value2=n.getValue();
1819     // get a non const pointer to the inside array of values and perform operation
1820     T * value=const_cast<T *> (getValue());
1821
1822     const int size=getNumberOfValues()*getNumberOfComponents();
1823     SCRUTE_MED(size);
1824     const T* endV1=value1+size;
1825     for(;value1!=endV1; value1++,value2++,value++)
1826         *value=(*value1)-(*value2);
1827 }
1828
1829 /*!
1830   \if developper
1831   Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
1832   This method is applied to a just created field with medModeSwitch mode.
1833   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1834   it doesn't exist!
1835   \endif
1836 */
1837 template <class T, class INTERLACING_TAG>
1838 void FIELD<T, INTERLACING_TAG>::_mul_in_place(const FIELD& m,const FIELD& n)
1839 {
1840     // get pointers to the values we are adding
1841     const T* value1=m.getValue();
1842     const T* value2=n.getValue();
1843     // get a non const pointer to the inside array of values and perform operation
1844     T * value=const_cast<T *> (getValue());
1845
1846     const int size=getNumberOfValues()*getNumberOfComponents();
1847     SCRUTE_MED(size);
1848     const T* endV1=value1+size;
1849     for(;value1!=endV1; value1++,value2++,value++)
1850         *value=(*value1)*(*value2);
1851 }
1852
1853 /*!
1854   \if developper
1855   Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
1856   This method is applied to a just created field with medModeSwitch mode.
1857   For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1858   it doesn't exist!
1859   \endif
1860 */
1861 template <class T, class INTERLACING_TAG>
1862 void FIELD<T, INTERLACING_TAG>::_div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION)
1863 {
1864     // get pointers to the values we are adding
1865     const T* value1=m.getValue();
1866     const T* value2=n.getValue();
1867     // get a non const pointer to the inside array of values and perform operation
1868     T * value=const_cast<T *> (getValue());
1869
1870     const int size=getNumberOfValues()*getNumberOfComponents();
1871     SCRUTE_MED(size);
1872     const T* endV1=value1+size;
1873     for(;value1!=endV1; value1++,value2++,value++){
1874       if ( *value2 == 0 ) { // FAIRE PLUTOT UN TRY CATCH Rmq from EF
1875           string diagnosis;
1876           diagnosis="FIELD<T,INTERLACING_TAG>::_div_in_place(...) : Divide by zero !";
1877           throw MEDEXCEPTION(diagnosis.c_str());
1878         }
1879         *value=(*value1)/(*value2);
1880     }
1881 }
1882
1883 /*!
1884 \addtogroup FIELD_algo
1885 @{
1886 */
1887
1888 /*!  Return maximum of all absolute values contained in the array (all elements and all components are browsed).
1889  */
1890 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::normMax() const throw (MEDEXCEPTION)
1891 {
1892     const T* value=getValue(); // get pointer to the values
1893     const int size=getNumberOfValues()*getNumberOfComponents();
1894     if (size <= 0) // Size of array has to be strictly positive
1895     {
1896         string diagnosis;
1897         diagnosis="FIELD<T,INTERLACIN_TAG>::normMax() : cannot compute the norm of "+getName()+
1898             " : it size is non positive!";
1899         throw MEDEXCEPTION(diagnosis.c_str());
1900     }
1901     const T* lastvalue=value+size; // get pointer just after last value
1902     const T* pMax=value; // pointer to the max value
1903     const T* pMin=value; // pointer to the min value
1904
1905     // get pointers to the max & min value of array
1906     while ( ++value != lastvalue )
1907     {
1908         if ( *pMin > *value )
1909             pMin=value;
1910         if ( *pMax < *value )
1911             pMax=value;
1912     }
1913
1914     T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
1915     T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
1916
1917     return Max>Min ? double(Max) : double(Min);
1918 }
1919
1920 /*!  Return Euclidian norm for all elements of the array.
1921  */
1922 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::norm2() const throw (MEDEXCEPTION)
1923 {
1924     const T* value=this->getValue(); // get const pointer to the values
1925     const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
1926     if (size <= 0) // Size of array has to be strictly positive
1927     {
1928         string diagnosis;
1929         diagnosis="FIELD<T,INTERLACIN_TAG>::norm2() : cannot compute the norm of "+getName()+
1930             " : it size is non positive!";
1931         throw MEDEXCEPTION(diagnosis.c_str());
1932     }
1933     const T* lastvalue=value+size; // point just after last value
1934
1935     T result((T)0); // init
1936     for( ; value!=lastvalue ; ++value)
1937         result += (*value) * (*value);
1938
1939     return std::sqrt(double(result));
1940 }
1941
1942
1943 //------------- TDG and BS addings
1944
1945 /*!  Return Extrema of field
1946  */
1947  template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION)
1948 {
1949   const T* value=getValue(); // get pointer to the values
1950   const int size=getNumberOfValues()*getNumberOfComponents();
1951   const T* lastvalue=value+size; // point just after last value
1952
1953   if (size <= 0){ // Size of array has to be strictly positive
1954
1955     string diagnosis;
1956     diagnosis="FIELD<T,INTERLACIN_TAG>::getMinMax() : cannot compute the extremums of "+getName()+
1957       " : its size is non positive!";
1958     throw MEDEXCEPTION(diagnosis.c_str());
1959   }
1960
1961   if (!_isMinMax){
1962     vmax=MinMax<T>::getMin(); // init a max value
1963     vmin=MinMax<T>::getMax(); // init a min value
1964
1965     for( ; value!=lastvalue ; ++value){
1966       if ( vmin > *value )
1967         vmin=*value;
1968       if ( vmax < *value )
1969         vmax=*value;
1970     }
1971     _isMinMax=true;
1972     _vmin=vmin;
1973     _vmax=vmax;
1974   }
1975   else{
1976     vmin = _vmin;
1977     vmax = _vmax;
1978   }
1979
1980 }
1981
1982 /*!  Return Histogram of field
1983  */
1984  template <class T, class INTERLACIN_TAG> vector<int> FIELD<T, INTERLACIN_TAG>::getHistogram(int &nbint) throw (MEDEXCEPTION)
1985 {
1986   const T* value=getValue(); // get pointer to the values
1987   const int size=getNumberOfValues()*getNumberOfComponents();
1988   const T* lastvalue=value+size; // point just after last value
1989
1990   if (size <= 0){ // Size of array has to be strictly positive
1991
1992     string diagnosis;
1993     diagnosis="FIELD<T,INTERLACIN_TAG>::getHistogram() : cannot compute the histogram of "+getName()+
1994       " : it size is non positive!";
1995     throw MEDEXCEPTION(diagnosis.c_str());
1996   }
1997
1998   vector<int> Histogram(nbint) ;
1999   T vmin,vmax;
2000   int j;
2001
2002   for( j=0 ; j!=nbint ; j++) Histogram[j]=0 ;
2003
2004   getMinMax(vmin,vmax);
2005   for( ; value!=lastvalue ; ++value){
2006     if(*value==vmax) j = nbint-1;
2007     else j = (int)(((double)nbint * (*value-vmin))/(vmax-vmin));
2008     Histogram[j]+=1 ;
2009   }
2010
2011   return Histogram ;
2012
2013 }
2014
2015 /*!  Return vectorial gradient field
2016  */
2017 template <class T, class INTERLACIN_TAG>
2018 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildGradient() const throw (MEDEXCEPTION)
2019 {
2020   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildGradient() : ";
2021   BEGIN_OF_MED(LOC);
2022
2023   // space dimension of input mesh
2024   int spaceDim = getSupport()->getMesh()->getSpaceDimension();
2025   double *x = new double[spaceDim];
2026
2027   FIELD<double, FullInterlace>* Gradient =
2028     new FIELD<double, FullInterlace>(getSupport(),spaceDim);
2029
2030   string name("gradient of ");
2031   name += getName();
2032   Gradient->setName(name);
2033   string descr("gradient of ");
2034   descr += getDescription();
2035   Gradient->setDescription(descr);
2036
2037   if( _numberOfComponents > 1 ){
2038     delete Gradient;
2039     delete [] x;
2040     throw MEDEXCEPTION("gradient calculation only on scalar field");
2041   }
2042
2043   for(int i=1;i<=spaceDim;i++){
2044     string nameC("gradient of ");
2045     nameC += getName();
2046     Gradient->setComponentName(i,nameC);
2047     Gradient->setComponentDescription(i,"gradient");
2048     string MEDComponentUnit = getMEDComponentUnit(1)+getSupport()->getMesh()->getCoordinatesUnits()[i-1];
2049     Gradient->setMEDComponentUnit(i,MEDComponentUnit);
2050   }
2051
2052   Gradient->setIterationNumber(getIterationNumber());
2053   Gradient->setOrderNumber(getOrderNumber());
2054   Gradient->setTime(getTime());
2055
2056   // typ of entity on what is field
2057   MED_EN::medEntityMesh typ = getSupport()->getEntity();
2058
2059   const int *C;
2060   const int *iC;
2061   const int *revC;
2062   const int *indC;
2063   const double *coord;
2064   int NumberOf;
2065
2066   switch (typ) {
2067   case MED_CELL:
2068   case MED_FACE:
2069   case MED_EDGE:
2070     {
2071       // read connectivity array to have the list of nodes contained by an element
2072       C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,typ,MED_ALL_ELEMENTS);
2073       iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,typ);
2074       // calculate reverse connectivity to have the list of elements which contains node i
2075       revC = getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
2076       indC = getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
2077       // coordinates of each node
2078       coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
2079       // number of elements
2080       NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2081       // barycenter field of elements
2082       FIELD<double, FullInterlace>* barycenter = getSupport()->getMesh()->getBarycenter(getSupport());
2083
2084       // calculate gradient vector for each element i
2085       for (int i = 1; i < NumberOf + 1; i++) {
2086
2087         // listElements contains elements which contains a node of element i
2088         set <int> listElements;
2089         set <int>::iterator elemIt;
2090         listElements.clear();
2091
2092         // for each node j of element i
2093         for (int ij = iC[i-1]; ij < iC[i]; ij++) {
2094           int j = C[ij-1];
2095           for (int k = indC[j-1]; k < indC[j]; k++) {
2096             // c element contains node j
2097             int c = revC[k-1];
2098             // we put the elements in set
2099             if (c != i)
2100               listElements.insert(c);
2101           }
2102         }
2103         // coordinates of barycentre of element i in space of dimension spaceDim
2104         for (int j = 0; j < spaceDim; j++)
2105           x[j] = barycenter->getValueIJ(i,j+1);
2106
2107         for (int j = 0; j < spaceDim; j++) {
2108           // value of field of element i
2109           double val = getValueIJ(i,1);
2110           double grad = 0.;
2111           // calculate gradient for each neighbor element
2112           for (elemIt = listElements.begin(); elemIt != listElements.end(); elemIt++) {
2113             int elem = *elemIt;
2114             double d2 = 0.;
2115             for (int l = 0; l < spaceDim; l++) {
2116               // coordinate of barycenter of element elem
2117               double xx = barycenter->getValueIJ(elem, l+1);
2118               d2 += (x[l]-xx) * (x[l]-xx);
2119             }
2120             grad += (barycenter->getValueIJ(elem,j+1)-x[j])*(getValueIJ(elem,1)-val)/sqrt(d2);
2121           }
2122           if (listElements.size() != 0) grad /= listElements.size();
2123           Gradient->setValueIJ(i,j+1,grad);
2124         }
2125       }
2126       delete barycenter;
2127     }
2128     break;
2129   case MED_NODE:
2130     // read connectivity array to have the list of nodes contained by an element
2131     C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
2132     iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,MED_CELL);
2133     // calculate reverse connectivity to have the list of elements which contains node i
2134     revC=getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
2135     indC=getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
2136     // coordinates of each node
2137     coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
2138
2139     // calculate gradient for each node
2140     NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2141     for (int i=1; i<NumberOf+1; i++){
2142       // listNodes contains nodes neigbor of node i
2143       set <int> listNodes;
2144       set <int>::iterator nodeIt ;
2145       listNodes.clear();
2146       for(int j=indC[i-1];j<indC[i];j++){
2147         // c element contains node i
2148         int c=revC[j-1];
2149         // we put the nodes of c element in set
2150         for(int k=iC[c-1];k<iC[c];k++)
2151           if(C[k-1] != i)
2152             listNodes.insert(C[k-1]);
2153       }
2154       // coordinates of node i in space of dimension spaceDim
2155       for(int j=0;j<spaceDim;j++)
2156         x[j] = coord[(i-1)*spaceDim+j];
2157
2158       for(int j=0;j<spaceDim;j++){
2159         // value of field
2160         double val = getValueIJ(i,1);
2161         double grad = 0.;
2162         // calculate gradient for each neighbor node
2163         for(nodeIt=listNodes.begin();nodeIt!=listNodes.end();nodeIt++){
2164           int node = *nodeIt;
2165           double d2 = 0.;
2166           for(int l=0;l<spaceDim;l++){
2167             double xx = coord[(node-1)*spaceDim+l];
2168             d2 += (x[l]-xx) * (x[l]-xx);
2169           }
2170           grad += (coord[(node-1)*spaceDim+j]-x[j])*(getValueIJ(node,1)-val)/sqrt(d2);
2171         }
2172         if(listNodes.size() != 0) grad /= listNodes.size();
2173         Gradient->setValueIJ(i,j+1,grad);
2174       }
2175     }
2176     break;
2177   case MED_ALL_ENTITIES:
2178     delete [] x;
2179     delete Gradient;
2180     throw MEDEXCEPTION("gradient calculation not yet implemented on all elements");
2181     break;
2182   }
2183
2184   delete [] x;
2185
2186   END_OF_MED(LOC);
2187   return Gradient;
2188 }
2189
2190 /*!  Return scalar norm2 field
2191  */
2192 template <class T, class INTERLACIN_TAG>
2193 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildNorm2Field() const throw (MEDEXCEPTION)
2194 {
2195   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildNorm2Field() : ";
2196   BEGIN_OF_MED(LOC);
2197
2198   FIELD<double, FullInterlace>* Norm2Field =
2199     new FIELD<double, FullInterlace>(getSupport(),1);
2200
2201   string name("norm2 of ");
2202   name += getName();
2203   Norm2Field->setName(name);
2204   string descr("norm2 of ");
2205   descr += getDescription();
2206   Norm2Field->setDescription(descr);
2207
2208   string nameC("norm2 of ");
2209   nameC += getName();
2210   Norm2Field->setComponentName(1,nameC);
2211   Norm2Field->setComponentDescription(1,"norm2");
2212   string MEDComponentUnit = getMEDComponentUnit(1);
2213   Norm2Field->setMEDComponentUnit(1,MEDComponentUnit);
2214
2215   Norm2Field->setIterationNumber(getIterationNumber());
2216   Norm2Field->setOrderNumber(getOrderNumber());
2217   Norm2Field->setTime(getTime());
2218
2219   // calculate nom2 for each element
2220   int NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2221   for (int i=1; i<NumberOf+1; i++){
2222     double norm2 = 0.;
2223     for(int j=1;j<=getNumberOfComponents();j++)
2224       norm2 += getValueIJ(i,j)*getValueIJ(i,j);
2225     Norm2Field->setValueIJ(i,1,sqrt(norm2));
2226   }
2227
2228   END_OF_MED(LOC);
2229   return Norm2Field;
2230
2231 }
2232
2233 /*!  Apply to each (scalar) field component the template parameter T_function,
2234  *   which is a pointer to function.
2235  *   Since the pointer is known at compile time, the function is inlined into the inner loop!
2236  *   calculation is done "in place".
2237  *   Use examples :
2238  *
2239  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
2240  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
2241  */
2242 template <class T, class INTERLACIN_TAG> template <T T_function(T)>
2243 void FIELD<T, INTERLACIN_TAG>::applyFunc()
2244 {
2245   // get a non const pointer to the inside array of values and perform operation
2246     T * value=const_cast<T *> (getValue());
2247     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
2248
2249     if (size>0) // for a negative size, there is nothing to do
2250     {
2251         const T* lastvalue=value+size; // pointer to the end of value
2252         for(;value!=lastvalue; ++value) // apply linear transformation
2253             *value = T_function(*value);
2254     }
2255 }
2256
2257 template <class T, class INTERLACIN_TAG> T FIELD<T, INTERLACIN_TAG>::pow(T x)
2258 {
2259   return (T)::pow((double)x,FIELD<T, INTERLACIN_TAG>::_scalarForPow);
2260 }
2261
2262 /*!  Apply to each (scalar) field component the math function pow.
2263  *   calculation is done "in place".
2264  *   Use examples :
2265  *
2266  *   \code  myField.applyFunc<std::sqrt>();  // apply sqare root function \endcode
2267  *     \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
2268  */
2269 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyPow(T scalar)
2270 {
2271   FIELD<T, INTERLACIN_TAG>::_scalarForPow=scalar;
2272   applyFunc<FIELD<T, INTERLACIN_TAG>::pow>();
2273 }
2274
2275 /*!  Apply to each (scalar) field component the linear function x -> ax+b.
2276  *   calculation is done "in place".
2277  */
2278 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b)
2279 {
2280     // get a non const pointer to the inside array of values and perform operation in place
2281     T * value=const_cast<T *> (getValue());
2282     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
2283
2284     if (size>0) // for a negative size, there is nothing to do
2285     {
2286         const T* lastvalue=value+size; // pointer to the end of value
2287         for(;value!=lastvalue; ++value) // apply linear transformation
2288             *value = a*(*value)+b;
2289     }
2290 }
2291
2292
2293 /*!
2294  *   Return a pointer to a new field that holds the scalar product. Static member function.
2295  *   This operation is authorized only for compatible fields that have the same support.
2296  *   The compatibility checking includes equality tests of the folowing data members:\n
2297  *   - _support
2298  *   - _numberOfComponents
2299  *   - _numberOfValues
2300  *   - _componentsTypes
2301  *   - _MEDComponentsUnits.
2302  *   Data members are initialized.
2303  *   The new field point to the same support and has one component.
2304  *   Each value of it is the scalar product of the two argument's fields.
2305  *   The user is in charge of memory deallocation.
2306  */
2307 template <class T, class INTERLACING_TAG>
2308 FIELD<T, INTERLACING_TAG>*
2309 FIELD<T, INTERLACING_TAG>::scalarProduct(const FIELD & m, const FIELD & n, bool deepCheck)
2310 {
2311     if(!deepCheck)
2312       FIELD_::_checkFieldCompatibility( m, n, false); // may throw exception
2313     else
2314       FIELD_::_deepCheckFieldCompatibility(m, n, false);
2315
2316     // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
2317     // result type imply INTERLACING_TAG=FullInterlace for m & n
2318
2319     const int numberOfElements=m.getNumberOfValues(); // strictly positive
2320     const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
2321
2322     // Creation & init of a the result field on the same support, with one component
2323     // You have to be careful about the interlacing mode, because in the computation step,
2324     // it seems to assume the the interlacing mode is the FullInterlacing
2325
2326     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),1);
2327     result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
2328     result->setIterationNumber(m.getIterationNumber());
2329     result->setTime(m.getTime());
2330     result->setOrderNumber(m.getOrderNumber());
2331
2332     const T* value1=m.getValue(); // get const pointer to the values
2333     const T* value2=n.getValue(); // get const pointer to the values
2334     // get a non const pointer to the inside array of values and perform operation
2335     T * value=const_cast<T *> (result->getValue());
2336
2337     const T* lastvalue=value+numberOfElements; // pointing just after last value of result
2338     for ( ; value!=lastvalue ; ++value ) // loop on all elements
2339     {
2340         *value=(T)0; // initialize value
2341         const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
2342         for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
2343             *value += (*value1) * (*value2);
2344     }
2345     return result;
2346 }
2347
2348 /*!  Return L2 Norm  of the field's component.
2349  *   Cannot be applied to a field with a support on nodes.
2350  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2351  *   For the nodal field, p_field_volume must be for all cells even if the field is partial.
2352  */
2353 template <class T, class INTERLACING_TAG>
2354 double FIELD<T, INTERLACING_TAG>::normL2(int component,
2355                                          const FIELD<double, FullInterlace> * p_field_volume) const
2356 {
2357     _checkNormCompatibility(p_field_volume, /*nodalAllowed=*/true); // may throw exception
2358     if ( component<1 || component>getNumberOfComponents() )
2359         throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
2360
2361     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2362     if(!p_field_volume) // if the user don't supply the volume
2363       p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
2364     else
2365       p_field_size->addReference();
2366     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2367     const double* vol=p_field_size->getValue();
2368     // Il n'est vraiment pas optimal de mixer des champs dans des modes d'entrelacement
2369     // different juste pour le calcul
2370
2371
2372     double integrale=0.0;
2373     double totVol=0.0;
2374
2375     if ( getSupport()->getEntity() == MED_NODE ) // issue 20120: [CEA 206] normL2 on NODE field
2376     {
2377       //Most frequently the FIELD is on the whole mesh and
2378       // there is no need in optimizing iterations from supporting nodes-> back to cells,
2379       // so we iterate just on all cells
2380       const MESH * mesh = getSupport()->getMesh()->convertInMESH();
2381       const int nbCells = mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
2382       const int *C = mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
2383       const int *iC = mesh->getConnectivityIndex(MED_NODAL,MED_CELL);
2384       for (int i = 0; i < nbCells; ++i, ++vol) {
2385         // calculate integral on current element as average summ of values on all it's nodes
2386         double curCellValue = 0;
2387         try { // we expect exception with partial fields for nodes w/o values
2388           for (int ij = iC[i]; ij < iC[i+1]; ij++) {
2389             int node = C[ij-1];
2390             curCellValue += getValueIJ( node, component );
2391           }
2392         }
2393         catch ( MEDEXCEPTION ) {
2394           continue;
2395         }
2396         int nbNodes = iC[i+1]-iC[i];
2397         curCellValue /= nbNodes;
2398         integrale += (curCellValue * curCellValue) * std::abs(*vol);
2399         totVol+=std::abs(*vol);
2400       }
2401       mesh->removeReference();
2402     }
2403     else
2404     {
2405       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2406         const T* value = getValue();
2407         value = value + (component-1) * getNumberOfValues();
2408         const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
2409         for (; value!=lastvalue ; ++value ,++vol) {
2410           integrale += double((*value) * (*value)) * std::abs(*vol);
2411           totVol+=std::abs(*vol);
2412         }
2413       }
2414       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2415         ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2416         for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2417           integrale += anArray->getIJ(i,component) * anArray->getIJ(i,component) * std::abs(*vol);
2418           totVol+=std::abs(*vol);
2419         }
2420       }
2421       else { // FULL_INTERLACE
2422         ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2423         for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2424           integrale += anArray->getIJ(i,component) * anArray->getIJ(i,component) * std::abs(*vol);
2425           totVol+=std::abs(*vol);
2426         }
2427       }
2428     }
2429
2430     if(p_field_size)
2431       p_field_size->removeReference(); // delete temporary volume field
2432
2433     if( totVol <= 0)
2434         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2435     return integrale/totVol;
2436 }
2437
2438 /*!  Return L2 Norm  of the field.
2439  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2440  *   For the nodal field, p_field_volume must be for all cells even if the field is partial.
2441  */
2442 template <class T, class INTERLACING_TAG>
2443 double FIELD<T, INTERLACING_TAG>::normL2(const FIELD<double, FullInterlace> * p_field_volume) const
2444 {
2445     _checkNormCompatibility(p_field_volume, /*nodalAllowed=*/true); // may throw exception
2446     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2447     if(!p_field_volume) // if the user don't supply the volume
2448       p_field_size=_getFieldSize(); // we calculate the volume
2449     else
2450       p_field_size->addReference();
2451     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2452     const double* vol=p_field_size->getValue();
2453     const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
2454
2455
2456     double integrale=0.0;
2457     double totVol=0.0;
2458
2459     if ( getSupport()->getEntity() == MED_NODE ) // issue 20120: [CEA 206] normL2 on NODE field
2460     {
2461       //Most frequently the FIELD is on the whole mesh and
2462       // there is no need in optimizing iterations from supporting nodes-> back to cells,
2463       // so we iterate just on all cells
2464       const MESH * mesh = getSupport()->getMesh()->convertInMESH();
2465       const int nbCells = mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
2466       const int *C = mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
2467       const int *iC = mesh->getConnectivityIndex(MED_NODAL,MED_CELL);
2468       int nbComp = getNumberOfComponents();
2469       for (int i = 0; i < nbCells; ++i, ++vol) {
2470         // calculate integral on current element as average summ of values on all it's nodes
2471         int nbNodes = iC[i+1]-iC[i];
2472         vector< double > curCellValue( nbComp, 0 );
2473         try { // we expect exception with partial fields for nodes w/o values
2474           for (int ij = iC[i]; ij < iC[i+1]; ij++) {
2475             int node = C[ij-1];
2476             for ( int j = 0; j < nbComp; ++j )
2477               curCellValue[ j ] += getValueIJ( node, j+1 ) / nbNodes;
2478           }
2479         }
2480         catch ( MEDEXCEPTION ) {
2481           continue;
2482         }
2483
2484         for ( int j = 0; j < nbComp; ++j ) {
2485           integrale += (curCellValue[j] * curCellValue[j]) * std::abs(*vol);
2486         }
2487         totVol+=std::abs(*vol);
2488       }
2489       mesh->removeReference();
2490       if ( nbCells > 0 && totVol == 0.)
2491         throw MEDEXCEPTION("can't compute sobolev norm : "
2492                            "none of elements has values on all it's nodes");
2493     }
2494     else
2495     {
2496       const double* p_vol=vol;
2497       for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2498         totVol+=std::abs(*p_vol);
2499
2500       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2501         const T* value = getValue();
2502         for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
2503           for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
2504             integrale += (*value) * (*value) * std::abs(*p_vol);
2505           }
2506         }
2507       }
2508       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2509         ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2510         for (int j=1; j<=anArray->getDim(); j++) {
2511           int i = 1;
2512           for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2513             integrale += anArray->getIJ(i,j) * anArray->getIJ(i,j) * std::abs(*p_vol);
2514           }
2515         }
2516       }
2517       else { // FULL_INTERLACE
2518         ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2519         for (int j=1; j<=anArray->getDim(); j++) {
2520           int i = 1;
2521           for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2522             integrale += anArray->getIJ(i,j) * anArray->getIJ(i,j) * std::abs(*p_vol);
2523           }
2524         }
2525       }
2526     }
2527     if(p_field_size)
2528         p_field_size->removeReference(); // delete temporary volume field
2529
2530     if( totVol <= 0)
2531       throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2532     return integrale/totVol;
2533 }
2534
2535 /*!  Return L1 Norm  of the field's component.
2536  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2537  */
2538 template <class T, class INTERLACING_TAG>
2539 double FIELD<T, INTERLACING_TAG>::normL1(int component,
2540                                          const FIELD<double, FullInterlace > * p_field_volume) const
2541 {
2542     _checkNormCompatibility(p_field_volume); // may throw exception
2543     if ( component<1 || component>getNumberOfComponents() )
2544         throw MEDEXCEPTION(STRING("FIELD<T,INTERLACING_TAG>::normL1() : The component argument should be between 1 and the number of components"));
2545
2546     const FIELD<double,FullInterlace> * p_field_size=p_field_volume;
2547     if(!p_field_volume) // if the user don't supply the volume
2548       p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃÃÅ mentation dans mesh]
2549     else
2550       p_field_size->addReference();
2551     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2552     const double* vol = p_field_size->getValue();
2553
2554     double integrale=0.0;
2555     double totVol=0.0;
2556
2557     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2558       const T* value = getValue();
2559       const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
2560       for (; value!=lastvalue ; ++value ,++vol) {
2561         integrale += std::abs( *value * *vol );
2562         totVol+=std::abs(*vol);
2563       }
2564     }
2565     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2566       ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2567       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2568         integrale += std::abs( anArray->getIJ(i,component) * (*vol));
2569         totVol+=std::abs(*vol);
2570       }
2571     }
2572     else { // FULL_INTERLACE
2573       ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2574       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2575         integrale += std::abs( anArray->getIJ(i,component) * *vol);
2576         totVol+=std::abs(*vol);
2577       }
2578     }
2579
2580     if(p_field_size)
2581       p_field_size->removeReference(); // delete temporary volume field
2582     if( totVol <= 0)
2583         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2584     return integrale/totVol;
2585 }
2586
2587 /*!  Return L1 Norm  of the field.
2588  *   Cannot be applied to a field with a support on nodes.
2589  *   If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2590  */
2591 template <class T, class INTERLACING_TAG>
2592 double FIELD<T, INTERLACING_TAG>::normL1(const FIELD<double, FullInterlace> * p_field_volume) const
2593 {
2594   _checkNormCompatibility(p_field_volume); // may throw exception
2595   const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2596   if(!p_field_volume) // if the user don't supply the volume
2597     p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃÂ©mentation dans mesh]
2598   else
2599     p_field_size->addReference();
2600   // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2601   const double* vol = p_field_size->getValue();
2602   const double* lastvol = vol+getNumberOfValues(); // pointing just after the end of vol
2603
2604   double integrale=0.0;
2605   double totVol=0.0;
2606   const double* p_vol=vol;
2607   for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2608     totVol+=std::abs(*p_vol);
2609
2610   if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2611     const T* value = getValue();
2612     for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
2613       for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
2614         integrale += std::abs( *value * *p_vol );
2615       }
2616     }
2617   }
2618   else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2619     ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2620     for (int j=1; j<=anArray->getDim(); j++) {
2621       int i = 1;
2622       for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2623         integrale += std::abs( anArray->getIJ(i,j) * *p_vol );
2624       }
2625     }
2626   }
2627   else { // FULL_INTERLACE
2628     ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2629     for (int j=1; j<=anArray->getDim(); j++) {
2630       int i = 1;
2631       for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2632         integrale += std::abs( anArray->getIJ(i,j) * *p_vol );
2633       }
2634     }
2635   }
2636   if(p_field_size)
2637     p_field_size->removeReference();
2638   if( totVol <= 0)
2639     throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2640   return integrale/totVol;
2641 }
2642
2643 /*!
2644  * \brief Return integral of the field.
2645  *  \param subSupport - optional part of a field to consider.
2646  *  \retval double - value of integral
2647  */
2648 template <class T, class INTERLACING_TAG>
2649 double FIELD<T, INTERLACING_TAG>::integral(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
2650 {
2651   const char* LOC = "FIELD<>::integral(subSupport): ";
2652
2653   double integrale = 0;
2654
2655   if (!subSupport ) subSupport = _support;
2656
2657   // check feasibility
2658   if ( getGaussPresence() )
2659     throw MEDEXCEPTION(STRING(LOC)<<"Gauss numbers greater than one are not yet implemented!");
2660   if ( subSupport->getEntity() != _support->getEntity())
2661     throw MEDEXCEPTION(STRING(LOC)<<"Different support entity of this field and subSupport");
2662   if ( subSupport->getEntity() == MED_EN::MED_NODE )
2663     throw MEDEXCEPTION(STRING(LOC)<<"Integral of nodal field not yet supported");
2664
2665   // analyze support
2666   const int nbElems = subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2667   const bool subOnAll = ( subSupport->isOnAllElements() );
2668   const bool  myOnAll = ( _support->isOnAllElements() );
2669   const int* subNums = !subOnAll ? subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS) : 0;
2670   const int*   myNums = !myOnAll ? _support->getNumber(MED_EN::MED_ALL_ELEMENTS) : 0;
2671   if ( !subOnAll && !subNums )
2672     throw MEDEXCEPTION(STRING(LOC)<<"Invalid support: no element numbers");
2673   if ( !myOnAll && !myNums )
2674     throw MEDEXCEPTION(STRING(LOC)<<"Invalid field support: no element numbers");
2675   if ( subOnAll && !myOnAll )
2676     return integral(NULL);
2677
2678   // get size of elements
2679   const FIELD<double, FullInterlace> * cellSize=_getFieldSize(subSupport);
2680   const double* size = cellSize->getValue();
2681   const double* lastSize = size + nbElems; // pointing just after the end of size
2682
2683   const T* value = getValue();
2684
2685   // calculate integrale
2686   if ( (subOnAll && _support->isOnAllElements()) || subSupport == _support )
2687     {
2688       const double* p_vol;
2689       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2690         {
2691           for (int j=1; j<=getNumberOfComponents(); ++j)
2692             for ( p_vol=size; p_vol != lastSize; ++value ,++p_vol)
2693               integrale += std::abs( *value * *p_vol );
2694         }
2695       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
2696         {
2697           typename ArrayNoByType::InterlacingPolicy* indexer =
2698             dynamic_cast< typename ArrayNoByType::InterlacingPolicy * > ( getArrayNoGauss() );
2699           for (int i, j=1; j<=getNumberOfComponents(); j++)
2700             for (i = 1, p_vol=size; p_vol!=lastSize; i++, ++p_vol )
2701               integrale += std::abs( value[indexer->getIndex(i,j)] * *p_vol );
2702         }
2703       else  // FULL_INTERLACE
2704         {
2705           for ( p_vol=size; p_vol != lastSize; ++p_vol)
2706             for (int j=0; j<getNumberOfComponents(); ++j, ++value)
2707               integrale += std::abs( *value * *p_vol );
2708         }
2709     }
2710   else
2711     {
2712       // find index for each element of subSupport
2713       PointerOf<int> index;
2714       if ( _support->isOnAllElements() )
2715         {
2716           index.set( subNums );
2717         }
2718       else // find common of two partial supports
2719         {
2720           // hope that numbers are in increasing order
2721           index.set( nbElems );
2722           for (int ii = 0; ii < nbElems; ii++)
2723             index[ii] = 0;
2724           bool allNumsFound = true;
2725           int i = 0, iSub = 0;
2726           for ( ; iSub < nbElems; ++iSub )
2727             {
2728               while ( i < getNumberOfValues() && subNums[iSub] > myNums[i] )
2729                 ++i;
2730               if (i == getNumberOfValues() /*subNums[iSub] > myNums[i]*/) // no more myNums
2731                 {
2732                   index[iSub] = 0; // no such number in myNums
2733                   break;
2734                 }
2735               else if ( subNums[iSub] == myNums[i] ) // elem number found
2736                 index[iSub] = ++i; // -- index counts from 1
2737               else // subNums[iSub] < myNums[i]
2738                 allNumsFound = (index[iSub] = 0); // no such number in myNums
2739             }
2740           if ( iSub != nbElems || !allNumsFound )
2741             {
2742               // check if numbers are in increasing order
2743               bool increasingOrder = true;
2744               for ( iSub = 1; iSub < nbElems && increasingOrder; ++iSub )
2745                 increasingOrder = ( subNums[iSub-1] < subNums[iSub] );
2746               for ( i = 1; i < getNumberOfValues() && increasingOrder; ++i )
2747                 increasingOrder = ( myNums[i-1] < myNums[i] );
2748
2749               if ( !increasingOrder )
2750                 for ( iSub = 0; iSub < nbElems; ++iSub )
2751                   try
2752                     {
2753                       index[iSub] = _support->getValIndFromGlobalNumber( subNums[iSub] );
2754                     }
2755                   catch (MEDEXCEPTION)
2756                     {
2757                       index[iSub] = 0;
2758                     }
2759             }
2760         }
2761
2762       // calculation
2763       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2764         {
2765           for (int j=0; j<getNumberOfComponents(); ++j)
2766             {
2767               value = getValue() + j * getNumberOfValues();
2768               for ( int i = 0; i < nbElems; ++i )
2769                 if ( index[i] )
2770                   integrale += std::abs( value[ index[i]-1 ] * size[i] );
2771             }
2772         }
2773       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
2774         {
2775           typename ArrayNoByType::InterlacingPolicy* indexer =
2776             dynamic_cast< typename ArrayNoByType::InterlacingPolicy * > ( getArrayNoGauss() );
2777           for (int j=1; j<=getNumberOfComponents(); j++)
2778             for ( int i = 0; i < nbElems; ++i )
2779               if ( index[i] )
2780                 integrale += std::abs( value[indexer->getIndex(index[i],j)] * size[i] );
2781         }
2782       else  // FULL_INTERLACE
2783         {
2784           const int dim = getNumberOfComponents();
2785           for ( int i = 0; i < nbElems; ++i )
2786             if ( index[i] )
2787               for (int j=0; j<dim; ++j)
2788                 integrale += std::abs( value[ dim*(index[i]-1) + j] * size[i] );
2789         }
2790     }
2791   cellSize->removeReference();
2792   return integrale;
2793 }
2794
2795 /*! Return a new field (to deallocate with delete) lying on subSupport that is included by
2796  *   this->_support with corresponding values extracting from this->_value.
2797  */
2798 template <class T, class INTERLACING_TAG>
2799 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
2800 {
2801   if(!subSupport->belongsTo(*_support))
2802     throw MEDEXCEPTION("FIELD<T>::extract : subSupport not included in this->_support !");
2803   if(_support->isOnAllElements() && subSupport->isOnAllElements())
2804     return new FIELD<T, INTERLACING_TAG>(*this);
2805
2806   FIELD<T, INTERLACING_TAG> *ret = new FIELD<T, INTERLACING_TAG>(subSupport,
2807                                                                  _numberOfComponents);
2808
2809   if(!ret->_value)
2810     throw MEDEXCEPTION("FIELD<T>::extract : invalid support detected !");
2811
2812   T* valuesToSet=(T*)ret->getValue();
2813
2814   int nbOfEltsSub=subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2815   const int *eltsSub=subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS);
2816   T* tempVals=new T[_numberOfComponents];
2817   for(int i=0;i<nbOfEltsSub;i++)
2818     {
2819       if(!getValueOnElement(eltsSub[i],tempVals))
2820         throw MEDEXCEPTION("Problem in belongsTo function !!!");
2821       for(int j=0;j<_numberOfComponents;j++)
2822         valuesToSet[i*_numberOfComponents+j]=tempVals[j];
2823     }
2824   delete [] tempVals;
2825
2826   ret->copyGlobalInfo(*this);
2827   return ret;
2828 }
2829 /*!
2830 @}
2831 */
2832
2833 /*!
2834         \addtogroup FIELD_io
2835         @{
2836 */
2837 /*!
2838   Constructor with parameters; the object is set via a file and its associated
2839   driver. For the moment only the MED_DRIVER is considered and if the last two
2840   argument (iterationNumber and orderNumber) are not set; their default value
2841   is -1. If the field fieldDriverName with the iteration number
2842   iterationNumber and the order number orderNumber does not exist in the file
2843   fieldDriverName; the constructor raises an exception.
2844 */
2845 template <class T, class INTERLACING_TAG>
2846 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
2847                                  driverTypes driverType,
2848                                  const string & fileName/*=""*/,
2849                                  const string & fieldDriverName/*=""*/,
2850                                  const int iterationNumber,
2851                                  const int orderNumber) throw (MEDEXCEPTION)
2852 {
2853   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) : ";
2854   BEGIN_OF_MED(LOC);
2855
2856   int current;
2857
2858   init();
2859
2860   _mesh  = ( MESH* ) NULL;
2861
2862   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2863   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2864   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
2865
2866   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2867   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2868   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2869
2870   _support = Support;
2871   //A.G. Addings for RC
2872   if(_support)
2873     _support->addReference();
2874   // OCC 10/03/2006 -- According to the rules defined with help of
2875   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2876   // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2877   // FIELD template - MSVC++ 2003 compiler generated an error here.
2878   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2879   _value = NULL;
2880
2881   _iterationNumber = iterationNumber;
2882   _time = 0.0;
2883   _orderNumber = orderNumber;
2884
2885   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
2886
2887   _drivers[current]->open();
2888   _drivers[current]->read();
2889   _drivers[current]->close();
2890
2891   END_OF_MED(LOC);
2892 }
2893
2894 /*!
2895   If the mesh argument is not initialized or passed NULL,
2896   this constructor, at least, allows to create a FIELD without creating any
2897   SUPPORT then without having to load a MESH object, a support is created. It
2898   provides the meshName related mesh but doesn't not set a mesh in the created
2899   support.
2900   If the passed mesh contains corresponding support, this support will be used
2901   for the field. This support will be found in mesh by name of one of profiles,
2902   on which the FIELD lays in MED-file. This has sense for the case, then MED-file
2903   was created by MEDMEM, and so name of profile contains name of corresponding support.
2904 */
2905 template <class T, class INTERLACING_TAG>
2906 FIELD<T,INTERLACING_TAG>::FIELD(driverTypes    driverType,
2907                                 const string & fileName,
2908                                 const string & fieldDriverName,
2909                                 const int      iterationNumber,
2910                                 const int      orderNumber,
2911                                 GMESH*         mesh)
2912   throw (MEDEXCEPTION) :FIELD_()
2913 {
2914   int current;
2915   const char* LOC = "FIELD<T,INTERLACING_TAG>::FIELD( driverTypes driverType, const string & fileName, string & fieldDriverName, int iterationNumber, int orderNumber) : ";
2916   BEGIN_OF_MED(LOC);
2917
2918   init();
2919
2920   _mesh = mesh;
2921   if(_mesh)
2922     _mesh->addReference();
2923
2924   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2925   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2926   FIELD_::_valueType = SET_VALUE_TYPE<T>::_valueType;
2927
2928   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2929   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2930   FIELD_::_interlacingType = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2931
2932   _support = (SUPPORT *) NULL;
2933   // OCC 10/03/2006 -- According to the rules defined with help of
2934   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2935   // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2936   // FIELD template - MSVC++ 2003 compiler generated an error here.
2937   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2938   _value = NULL;
2939
2940   _iterationNumber = iterationNumber;
2941   _time = 0.0;
2942   _orderNumber = orderNumber;
2943
2944   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
2945
2946   _drivers[current]->open();
2947   _drivers[current]->read();
2948   _drivers[current]->close();
2949
2950   END_OF_MED(LOC);
2951 }
2952 /*!
2953 @}
2954 */
2955
2956 /*!
2957   Destructor.
2958 */
2959 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::~FIELD()
2960 {
2961   const char* LOC = " Destructeur FIELD<T, INTERLACING_TAG>::~FIELD()";
2962   BEGIN_OF_MED(LOC);
2963   SCRUTE_MED(this);
2964   if (_value) delete _value; _value=0;
2965   locMap::const_iterator it;
2966   for ( it = _gaussModel.begin();it != _gaussModel.end(); it++ )
2967     delete (*it).second;
2968   _gaussModel.clear();
2969   if(_mesh)
2970     _mesh->removeReference();
2971   _mesh=0;
2972   END_OF_MED(LOC);
2973 }
2974
2975 /*!
2976
2977 */
2978 template <class T, class INTERLACING_TAG>
2979 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)
2980 {
2981   const char* LOC = "FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)";
2982   BEGIN_OF_MED(LOC);
2983
2984   _numberOfComponents = NumberOfComponents ;
2985   _componentsTypes.resize(NumberOfComponents);
2986   _componentsNames.resize(NumberOfComponents);
2987   _componentsDescriptions.resize(NumberOfComponents);
2988   _componentsUnits.resize(NumberOfComponents);
2989   _MEDComponentsUnits.resize(NumberOfComponents);
2990   for (int i=0;i<NumberOfComponents;i++) {
2991     _componentsTypes[i] = 0 ;
2992   }
2993   delete _value;
2994   try {
2995     // becarefull about the number of gauss point
2996     _numberOfValues = _support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2997     MESSAGE_MED(PREFIX_MED <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
2998
2999     //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
3000     _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
3001
3002     _isRead = true ;
3003   }
3004   catch (MEDEXCEPTION &ex) {
3005     MESSAGE_MED("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
3006     // OCC 10/03/2006 -- According to the rules defined with help of
3007     // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
3008     // MEDMEM_Array<> template using INTERLACING_TAG parameter of
3009     // FIELD template - MSVC++ 2003 compiler generated an error here.
3010     // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
3011     _value = NULL;
3012   }
3013
3014   SCRUTE_MED(_value);
3015   END_OF_MED(LOC);
3016 }
3017
3018 /*!
3019
3020 */
3021 template <class T, class INTERLACING_TAG>
3022 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents,
3023                                            const int LengthValue)
3024 {
3025   const char* LOC = "void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)";
3026   BEGIN_OF_MED(LOC);
3027
3028   _numberOfComponents = NumberOfComponents ;
3029   _componentsTypes.resize(NumberOfComponents);
3030   _componentsNames.resize(NumberOfComponents);
3031   _componentsDescriptions.resize(NumberOfComponents);
3032   _componentsUnits.resize(NumberOfComponents);
3033   _MEDComponentsUnits.resize(NumberOfComponents);
3034   for (int i=0;i<NumberOfComponents;i++) {
3035     _componentsTypes[i] = 0 ;
3036   }
3037
3038   MESSAGE_MED("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
3039   _numberOfValues = LengthValue ;
3040   delete _value;
3041   //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
3042   _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
3043
3044   _isRead = true ;
3045
3046   SCRUTE_MED(_value);
3047   END_OF_MED(LOC);
3048 }
3049
3050 /*!
3051
3052 */
3053 template <class T, class INTERLACING_TAG>
3054 void FIELD<T, INTERLACING_TAG>::deallocValue()
3055 {
3056   const char* LOC = "void FIELD<T, INTERLACING_TAG>::deallocValue()";
3057   BEGIN_OF_MED(LOC);
3058   _numberOfValues = 0 ;
3059   _numberOfComponents = 0 ;
3060   if (_value != NULL) {
3061     delete _value;
3062     _value = NULL;
3063   }
3064
3065   END_OF_MED(LOC);
3066 }
3067
3068
3069
3070 /*!\if MEDMEM_ug
3071
3072 \addtogroup FIELD_io
3073 @{
3074 \endif */
3075 // -----------------
3076 // Methodes Inline
3077 // -----------------
3078
3079 /*!
3080   Creates the specified driver and return its index reference to path to
3081   read or write methods.
3082
3083 \param driverType specifies the file type (MED_DRIVER or VTK_DRIVER)
3084 \param fileName name of the output file
3085 \param driverName name of the field
3086 \param access access type (read, write or both)
3087
3088 */
3089
3090 template <class T, class INTERLACING_TAG>
3091 int FIELD<T, INTERLACING_TAG>::addDriver(driverTypes driverType,
3092                                          const string & fileName/*="Default File Name.med"*/,
3093                                          const string & driverName/*="Default Field Name"*/,
3094                                          MED_EN::med_mode_acces access)
3095 {
3096   //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) : ";
3097
3098   GENDRIVER * driver;
3099
3100   const char* LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName,const string & driverName,MED_EN::med_mode_acces access) :";
3101   BEGIN_OF_MED(LOC);
3102
3103   SCRUTE_MED(driverType);
3104
3105   driver = DRIVERFACTORY::buildDriverForField(driverType,fileName,this,access);
3106
3107   _drivers.push_back(driver);
3108
3109   int current = _drivers.size()-1;
3110
3111   _drivers[current]->setFieldName(driverName);
3112
3113   END_OF_MED(LOC);
3114
3115   return current;
3116 }
3117
3118 /*!
3119   Read FIELD in the file specified in the driver given by its index.
3120 */
3121 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(int index/*=0*/)
3122 {
3123   const char * LOC = "FIELD<T, INTERLACING_TAG>::read(int index=0) : ";
3124   BEGIN_OF_MED(LOC);
3125
3126   if ( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3127     _drivers[index]->open();
3128     try
3129     {
3130       _drivers[index]->read();
3131     }
3132     catch ( const MEDEXCEPTION& ex )
3133     {
3134       _drivers[index]->close();
3135       throw ex;
3136     }
3137     _drivers[index]->close();
3138   }
3139   else
3140     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3141                                      << "The index given is invalid, index must be between  0 and |"
3142                                      << _drivers.size()
3143                                      )
3144                           );
3145   END_OF_MED(LOC);
3146 }
3147
3148 /*!
3149   Read FIELD with the driver. Additional information (name etc.) to select a field
3150   must be set to the field.
3151 */
3152 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(const GENDRIVER & driver )
3153 {
3154   const char* LOC = " FIELD<T, INTERLACING_TAG>::read(const GENDRIVER &) : ";
3155   BEGIN_OF_MED(LOC);
3156
3157   // For the case where driver does not know about me since it has been created through
3158   // constructor witout parameters: create newDriver knowing me and get missing data
3159   // from driver using merge()
3160   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
3161                                                                          driver.getFileName(),
3162                                                                          this, RDONLY));
3163   newDriver->merge( driver );
3164
3165   newDriver->open();
3166   try
3167   {
3168     newDriver->read();
3169   }
3170   catch ( const MEDEXCEPTION& ex )
3171   {
3172     newDriver->close();
3173     throw ex;
3174   }
3175   newDriver->close();
3176
3177   END_OF_MED(LOC);
3178 }
3179
3180 /*!
3181   Read FIELD with driver of the given type. Additional information (name etc.) to select a field
3182   must be set to the field.
3183 */
3184 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(driverTypes driverType, const std::string& filename)
3185 {
3186   const char* LOC = " FIELD<T, INTERLACING_TAG>::read(driverTypes driverType, const std::string& filename) : ";
3187   BEGIN_OF_MED(LOC);
3188
3189   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driverType, filename,
3190                                                                          this, RDONLY));
3191   newDriver->open();
3192   try
3193   {
3194     newDriver->read();
3195   }
3196   catch ( const MEDEXCEPTION& ex )
3197   {
3198     newDriver->close();
3199     throw ex;
3200   }
3201   newDriver->close();
3202
3203   END_OF_MED(LOC);
3204 }
3205
3206 /*! \if MEDMEM_ug @} \endif */
3207
3208 /*!
3209   Duplicates the given driver and return its index reference to path to
3210   read or write methods.
3211 */
3212 template <class T, class INTERLACING_TAG>
3213 inline int FIELD<T, INTERLACING_TAG>::addDriver (GENDRIVER & driver )
3214 {
3215   int current;
3216
3217   const char* LOC = "FIELD<T, INTERLACING_TAG>::addDriver(GENDRIVER &) : ";
3218   BEGIN_OF_MED(LOC);
3219
3220   GENDRIVER * newDriver =
3221     DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
3222                                        driver.getFileName(), this,
3223                                        driver.getAccessMode());
3224   _drivers.push_back(newDriver);
3225
3226   current = _drivers.size()-1;
3227   SCRUTE_MED(current);
3228   driver.setId(current);
3229
3230   newDriver->merge( driver );
3231   newDriver->setId(current);
3232
3233   return current ;
3234 }
3235
3236 /*!
3237   Remove the driver referenced by its index.
3238 */
3239 template <class T, class INTERLACING_TAG>
3240 void FIELD<T, INTERLACING_TAG>::rmDriver (int index/*=0*/)
3241 {
3242   const char * LOC = "FIELD<T, INTERLACING_TAG>::rmDriver (int index=0): ";
3243   BEGIN_OF_MED(LOC);
3244
3245   if ( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3246     MESSAGE_MED ("detruire");
3247   }
3248   else
3249     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3250                                      << "The <index given is invalid, index must be between  0 and  |"
3251                                      << _drivers.size()
3252                                      )
3253                           );
3254
3255   END_OF_MED(LOC);
3256 }
3257
3258 /*! \if MEDMEM_ug
3259 \addtogroup FIELD_io
3260 @{
3261 \endif */
3262
3263 /*!
3264   Writes FIELD in the file specified by the driver handle \a index.
3265
3266 Example :
3267 \verbatim
3268 //...
3269 // Attaching the friver to file "output.med", meshname "Mesh"
3270 int driver_handle = mesh.addDriver(MED_DRIVER, "output.med", "Mesh");
3271 // Writing the content of mesh to the file
3272 mesh.write(driver_handle);
3273 \endverbatim
3274 */
3275 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(int index/*=0*/)
3276 {
3277   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
3278   BEGIN_OF_MED(LOC);
3279
3280   if( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3281     _drivers[index]->open();
3282     try
3283     {
3284       _drivers[index]->write();
3285     }
3286     catch ( const MEDEXCEPTION& ex )
3287     {
3288       _drivers[index]->close();
3289       throw ex;
3290     }
3291     _drivers[index]->close();
3292   }
3293   else
3294     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3295                                      << "The index given is invalid, index must be between  0 and |"
3296                                      << _drivers.size()
3297                                      )
3298                           );
3299   END_OF_MED(LOC);
3300 }
3301 /*!
3302   Write FIELD with the given driver.
3303 */
3304 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(const GENDRIVER & driver, MED_EN::med_mode_acces medMode/*=MED_EN::RDWR*/)
3305 {
3306   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
3307   BEGIN_OF_MED(LOC);
3308
3309   // For the case where driver does not know about me since it has been created through
3310   // constructor witout parameters: create newDriver knowing me and get missing data
3311   // from driver using merge()
3312   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
3313                                                                          driver.getFileName(),
3314                                                                          this, WRONLY));
3315   newDriver->merge( driver );
3316   if ( newDriver->getDriverType() == MED_DRIVER )
3317     newDriver->setAccessMode( MED_EN::med_mode_acces( getMedAccessMode( medMode ) ));
3318
3319   newDriver->open();
3320   try
3321   {
3322     newDriver->write();
3323   }
3324   catch ( const MEDEXCEPTION& ex )
3325   {
3326     newDriver->close();
3327     throw ex;
3328   }
3329   newDriver->close();
3330
3331   END_OF_MED(LOC);
3332 }
3333
3334 /*!
3335   Write FIELD with driver of the given type.
3336 */
3337 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(driverTypes driverType, const std::string& filename, MED_EN::med_mode_acces medMode/*=MED_EN::RDWR*/)
3338 {
3339   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(driverTypes driverType, const std::string& filename) : ";
3340   BEGIN_OF_MED(LOC);
3341
3342   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driverType, filename,
3343                                                                          this, WRONLY));
3344   if ( newDriver->getDriverType() == MED_DRIVER )
3345     newDriver->setAccessMode( MED_EN::med_mode_acces( getMedAccessMode( medMode ) ));
3346
3347   newDriver->open();
3348   try
3349   {
3350     newDriver->write();
3351   }
3352   catch ( const MEDEXCEPTION& ex )
3353   {
3354     newDriver->close();
3355     throw ex;
3356   }
3357   newDriver->close();
3358
3359   END_OF_MED(LOC);
3360 }
3361
3362 /*! \if MEDMEM_ug @} \endif */
3363 /*!
3364   Write FIELD in the file specified in the driver given by its index. Use this
3365   method for ASCII drivers (e.g. VTK_DRIVER)
3366 */
3367 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
3368 {
3369   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
3370   BEGIN_OF_MED(LOC);
3371
3372   if( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3373     _drivers[index]->openAppend();
3374     if (driverName != "") _drivers[index]->setFieldName(driverName);
3375     try
3376     {
3377       _drivers[index]->writeAppend();
3378     }
3379     catch ( const MEDEXCEPTION& ex )
3380     {
3381       _drivers[index]->close();
3382       throw ex;
3383     }
3384     _drivers[index]->close();
3385   }
3386   else
3387     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3388                                      << "The index given is invalid, index must be between  0 and |"
3389                                      << _drivers.size()
3390                                      )
3391                           );
3392   END_OF_MED(LOC);
3393 }
3394
3395 /*!
3396   \internal
3397   Write FIELD with the driver which is equal to the given driver.
3398
3399   Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
3400 */
3401 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(const GENDRIVER & genDriver)
3402 {
3403   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
3404   BEGIN_OF_MED(LOC);
3405
3406   for (unsigned int index=0; index < _drivers.size(); index++ )
3407     if ( *_drivers[index] == genDriver ) {
3408       _drivers[index]->openAppend();
3409       try
3410       {
3411         _drivers[index]->writeAppend();
3412       }
3413       catch ( const MEDEXCEPTION& ex )
3414       {
3415         _drivers[index]->close();
3416         throw ex;
3417       }
3418       _drivers[index]->close();
3419     }
3420
3421   END_OF_MED(LOC);
3422
3423 }
3424
3425 /*!
3426   Fills in already allocated retValues array the values related to eltIdInSup.
3427   If the element does not exist in this->_support false is returned, true otherwise.
3428 */
3429 template <class T, class INTERLACING_TAG>
3430 bool FIELD<T, INTERLACING_TAG>::getValueOnElement(int eltIdInSup,T* retValues)
3431   const throw (MEDEXCEPTION)
3432 {
3433
3434   if(eltIdInSup<1)
3435     return false;
3436   if(_support->isOnAllElements())
3437     {
3438       int nbOfEltsThis=_support->getMesh()->getNumberOfElements(_support->getEntity(),MED_EN::MED_ALL_ELEMENTS);
3439       if(eltIdInSup>nbOfEltsThis)
3440         return false;
3441       const T* valsThis=getValue();
3442       for(int j=0;j<_numberOfComponents;j++)
3443         retValues[j]=valsThis[(eltIdInSup-1)*_numberOfComponents+j];
3444       return true;
3445     }
3446   else
3447     {
3448       int nbOfEltsThis=_support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
3449       const int *eltsThis=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
3450       int iThis;
3451       bool found=false;
3452       for(iThis=0;iThis<nbOfEltsThis && !found;)
3453         if(eltsThis[iThis]==eltIdInSup)
3454           found=true;
3455         else
3456           iThis++;
3457       if(!found)
3458         return false;
3459       const T* valsThis=getValue();
3460       for(int j=0;j<_numberOfComponents;j++)
3461         retValues[j]=valsThis[iThis*_numberOfComponents+j];
3462       return true;
3463     }
3464 }
3465
3466   /*!
3467    * \brief Retrieve value in a given point
3468    *  \param coords - point coordinates
3469    *  \param output - output buffer
3470    */
3471   template <class T, class INTERLACING_TAG>
3472   void FIELD<T, INTERLACING_TAG>::getValueOnPoint(const double* coords, double* output) const throw (MEDEXCEPTION)
3473   {
3474     getValueOnPoints(1, coords, output);
3475   }
3476
3477   /*!
3478    * \brief Retrieve values in given points
3479    *  \param nb_points - number of points
3480    *  \param coords - point coordinates
3481    *  \param output - output buffer
3482    */
3483   template <class T, class INTERLACING_TAG>
3484   void FIELD<T, INTERLACING_TAG>::getValueOnPoints(int nb_points, const double* coords, double* output) const throw (MEDEXCEPTION)
3485   {
3486     const char* LOC = " FIELD<T, INTERLACING_TAG>::getValueOnPoints(int nb_points, const double* coords, double* output) : ";
3487     // check operation feasibility
3488     if ( !getSupport() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL Support"));
3489     if ( !getSupport()->getMesh() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL Mesh"));
3490     if ( !_value ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL _value"));
3491     if ( getGaussPresence() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not implemeneted for Gauss points"));
3492
3493     MED_EN::medEntityMesh entity = getSupport()->getEntity();
3494     if ( entity != MED_EN::MED_CELL &&
3495          entity != MED_EN::MED_NODE )
3496       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on CELLs or NODEs"));
3497
3498     // initialize output value
3499     for ( int j = 0; j < nb_points*getNumberOfComponents(); ++j )
3500       output[j] = 0.0;
3501
3502     const MEDMEM::MESH* mesh = getSupport()->getMesh()->convertInMESH();
3503     MEDMEM::AutoDeref derefMesh( mesh );
3504
3505     const double* point = coords;
3506     double* value = output;
3507
3508     if ( entity == MED_EN::MED_CELL )
3509       {
3510         MEDMEM::PointLocator pLocator (*mesh);
3511         for ( int i = 0; i < nb_points; ++i)
3512           {
3513             // find the cell enclosing the point
3514             std::list<int> cellIds = pLocator.locate( point );
3515             int nbCells = cellIds.size();
3516             if ( nbCells < 1 )
3517               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Point is out of mesh"));
3518
3519             // retrieve value
3520             std::list<int>::iterator iCell = cellIds.begin();
3521             for ( ; iCell != cellIds.end(); ++iCell )
3522               for ( int j = 1; j <= getNumberOfComponents(); ++j )
3523                 value[j-1] += getValueIJ( *iCell, j ) / nbCells;
3524
3525             // next point
3526             point += mesh->getSpaceDimension();
3527             value += getNumberOfComponents();
3528           }
3529       }
3530     else // MED_EN::MED_NODE
3531       {
3532         const double * allCoords = mesh->getCoordinates( MED_EN::MED_FULL_INTERLACE );
3533
3534         MEDMEM::PointLocatorInSimplex pLocator (*mesh);
3535         for ( int i = 0; i < nb_points; ++i)
3536           {
3537             // find nodes of the simplex enclosing the point
3538             std::list<int> nodeIds = pLocator.locate( point );
3539             int nbNodes = nodeIds.size();
3540             if ( nbNodes < 1 )
3541               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Point is out of mesh"));
3542             if ( nbNodes != mesh->getMeshDimension() + 1 )
3543               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid nb of points of simplex: "<<nbNodes));
3544
3545             // get coordinates of simplex nodes
3546             std::vector<const double*> nodeCoords( nbNodes );
3547             std::list<int>::iterator iNode = nodeIds.begin();
3548             int n = 0;
3549             for ( ; n < nbNodes; ++iNode, ++n )
3550               nodeCoords[n] = allCoords + (*iNode-1) * mesh->getSpaceDimension();
3551
3552             // compute wegths of simplex nodes
3553             double nodeWgt[4];
3554             pLocator.getNodeWightsInSimplex( nodeCoords, coords, nodeWgt );
3555
3556             // retrieve value
3557             for ( n = 0, iNode = nodeIds.begin(); iNode != nodeIds.end(); ++iNode, ++n )
3558               for ( int j = 1; j <= getNumberOfComponents(); ++j )
3559                 value[j-1] += getValueIJ( *iNode, j ) * nodeWgt[ n ];
3560
3561             // next point
3562             point += mesh->getSpaceDimension();
3563             value += getNumberOfComponents();
3564           }
3565       }
3566   }
3567
3568
3569 /*!
3570   Return the coordinates of the gauss points
3571   The returned field has SPACEDIM components
3572 */
3573 template <class T, class INTERLACING_TAG>
3574 FIELD<double, FullInterlace>* FIELD<T, INTERLACING_TAG>::getGaussPointsCoordinates() const
3575   throw (MEDEXCEPTION)
3576 {
3577   const char * LOC = "FIELD::getGaussPointsCoordinates() : ";
3578   BEGIN_OF_MED(LOC);
3579
3580   if (!getSupport())
3581     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3582
3583   const MEDMEM::MESH* mesh = getSupport()->getMesh()->convertInMESH();
3584   MEDMEM::AutoDeref derefMesh( mesh );
3585
3586   const double * coord = mesh->getCoordinates(MED_FULL_INTERLACE);
3587   int spaceDim         = mesh->getSpaceDimension();
3588
3589   //Init calculator of the gauss point coordinates
3590   INTERP_KERNEL::GaussCoords calculator;
3591   locMap::const_iterator it;
3592
3593   int nb_type                     = getSupport()->getNumberOfTypes();
3594   int length_values               = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
3595   const medGeometryElement* types = getSupport()->getTypes();
3596   medEntityMesh support_entity    = getSupport()->getEntity();
3597   bool isOnAll                    = getSupport()->isOnAllElements();
3598
3599   const int* global_connectivity  = 0;
3600   const GAUSS_LOCALIZATION<INTERLACING_TAG>* gaussLock = NULL;
3601
3602   typedef typename MEDMEM_ArrayInterface<double,INTERLACING_TAG,NoGauss>::Array ArrayCoord;
3603   typedef typename MEDMEM_ArrayInterface<double,INTERLACING_TAG,Gauss>::Array TArrayGauss;
3604
3605   vector<int>  nbelgeoc, nbgaussgeo;
3606
3607   nbelgeoc.resize(nb_type+1, 0);
3608   nbgaussgeo.resize(nb_type+1, 0);
3609
3610   for ( int iType = 0 ; iType < nb_type ; iType++ ) {
3611
3612     medGeometryElement elem_type = types[iType] ;
3613     if(elem_type == MED_EN::MED_POLYGON && elem_type == MED_EN::MED_POLYHEDRA )
3614       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad cell type : "<<MED_EN::geoNames[elem_type]<<" !!! "));
3615
3616     it = _gaussModel.find(elem_type);
3617
3618     if(it == _gaussModel.end())
3619       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Gauss localization not defined for "<<MED_EN::geoNames[elem_type]<<" type!!! "));
3620     gaussLock = static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second);
3621
3622     ArrayCoord coord = gaussLock->getGsCoo();
3623     double* gaussCoord = new double[coord.getNbElem()*coord.getDim()];
3624     int idx = 0;
3625     for( int i = 1 ; i <= coord.getNbElem() ; i++ ) {
3626       for( int j = 1 ; j <= coord.getDim() ; j++ ) {
3627         gaussCoord[idx++] = coord.getIJ(i,j);
3628       }
3629     }
3630
3631     idx = 0;
3632     ArrayCoord ref = gaussLock->getRefCoo();
3633     double* refCoord = new double[ref.getNbElem()*ref.getDim()];
3634     for( int i = 1 ; i <= ref.getNbElem() ; i++ ) {
3635       for( int j = 1 ; j <= ref.getDim() ; j++ ) {
3636         refCoord[idx++] = ref.getIJ(i,j);
3637       }
3638     }
3639
3640     INTERP_KERNEL::NormalizedCellType normType;
3641     switch(elem_type) {
3642     case MED_EN::MED_SEG2 : normType = INTERP_KERNEL::NORM_SEG2;break;
3643     case MED_EN::MED_SEG3 : normType = INTERP_KERNEL::NORM_SEG3;break;
3644     default : normType = (INTERP_KERNEL::NormalizedCellType) ((((unsigned long)elem_type/100-2)*10) + ((unsigned long)elem_type%100));
3645       break;
3646     }
3647
3648     calculator.addGaussInfo(normType,
3649                             elem_type/100,
3650                             gaussCoord,
3651                             gaussLock->getNbGauss(),
3652                             refCoord,
3653                             elem_type%100
3654                             );
3655     // Preapre Info for the gauss array
3656     nbelgeoc   [ iType+1 ] = nbelgeoc[ iType ] + getSupport()->getNumberOfElements(elem_type);
3657     nbgaussgeo [ iType+1 ] = gaussLock->getNbGauss();
3658
3659     delete [] gaussCoord;
3660     delete [] refCoord;
3661   }
3662
3663   FIELD<double, FullInterlace>* gpCoord =
3664     new FIELD<double, FullInterlace>(getSupport(),spaceDim);
3665   gpCoord->setName("Gauss Points Coordinates");
3666   gpCoord->setDescription("Gauss Points Coordinates");
3667
3668   for(int dimId = 1 ; dimId <= spaceDim; dimId++) {
3669     switch(dimId) {
3670     case 1:
3671       gpCoord->setComponentName(dimId,"X");
3672       gpCoord->setComponentDescription(dimId,"X coordinate of the gauss point");
3673       break;
3674     case 2:
3675       gpCoord->setComponentName(dimId,"Y");
3676       gpCoord->setComponentDescription(dimId,"Y coordinate of the gauss point");
3677       break;
3678     case 3:
3679       gpCoord->setComponentName(dimId,"Z");
3680       gpCoord->setComponentDescription(dimId,"Z coordinate of the gauss point");
3681       break;
3682     }
3683
3684     gpCoord->setMEDComponentUnit(dimId, mesh->getCoordinatesUnits()[dimId-1]);
3685   }
3686
3687   gpCoord->setIterationNumber(getIterationNumber());
3688   gpCoord->setOrderNumber(getOrderNumber());
3689   gpCoord->setTime(getTime());
3690
3691   TArrayGauss *arrayGauss = new TArrayGauss(spaceDim, length_values,
3692                                             nb_type, & nbelgeoc[0], & nbgaussgeo[0]);
3693   gpCoord->setArray(arrayGauss);
3694
3695
3696
3697
3698   //Calculation of the coordinates
3699   int index = 1;
3700   for ( int i = 0 ; i < nb_type ; i++ ) {
3701
3702     medGeometryElement type = types[i] ;
3703     INTERP_KERNEL::NormalizedCellType normType;
3704     switch(type) {
3705     case MED_EN::MED_SEG2 : normType = INTERP_KERNEL::NORM_SEG2;break;
3706     case MED_EN::MED_SEG3 : normType = INTERP_KERNEL::NORM_SEG3;break;
3707     default : normType = (INTERP_KERNEL::NormalizedCellType) ((((unsigned long)type/100-2)*10) + ((unsigned long)type%100));
3708       break;
3709     }
3710
3711     it = _gaussModel.find(type);
3712
3713     gaussLock = static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second);
3714     int nb_entity_type = getSupport()->getNumberOfElements(type);
3715
3716
3717     if (isOnAll) {
3718       global_connectivity = mesh->getConnectivity(MED_NODAL,support_entity,type);
3719     }
3720     else {
3721       const int * supp_number = getSupport()->getNumber(type);
3722       const int * connectivity = mesh->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
3723       const int * connectivityIndex = mesh->getConnectivityIndex(MED_NODAL,support_entity);
3724       int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
3725
3726       for (int k_type = 0; k_type<nb_entity_type; k_type++) {
3727         for (int j_ent = 0; j_ent<(type%100); j_ent++) {
3728           global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
3729         }
3730       }
3731       global_connectivity = global_connectivity_tmp;
3732     }
3733
3734     int nbNodes = (type%100);
3735     double* gCoord = NULL;
3736     int* Ni = NULL;
3737
3738     for ( int elem = 0; elem < nb_entity_type; elem++ ) {
3739       int elem_index = nbNodes*elem;
3740       Ni = new int[nbNodes];
3741       for( int idx = 0 ; idx < nbNodes; idx++ ) {
3742         Ni[idx] = global_connectivity[ elem_index+idx ] - 1;
3743       }
3744
3745       gCoord = calculator.calculateCoords(normType,
3746                                           coord,
3747                                           spaceDim,
3748                                           Ni);
3749       int resultIndex = 0;
3750       for( int k = 0; k < gaussLock->getNbGauss(); k++ ) {
3751         for( int dimId = 1; dimId <= spaceDim; dimId++ ) {
3752           gpCoord->setValueIJK(index,dimId,(k+1),gCoord[resultIndex]);
3753           resultIndex++;
3754         }
3755       }
3756       delete [] gCoord;
3757       delete [] Ni;
3758       index++;
3759     }
3760     if (!isOnAll && type != MED_EN::MED_POLYHEDRA && type != MED_EN::MED_POLYGON) {
3761       delete [] global_connectivity ;
3762     }
3763   }
3764   END_OF_MED(LOC);
3765   return gpCoord;
3766 }
3767
3768 /*!
3769   \if developper
3770   Destroy the MEDARRAY<T> in FIELD and put the new one without copy.
3771   \endif
3772 */
3773 template <class T, class INTERLACING_TAG>
3774 inline void FIELD<T, INTERLACING_TAG>::setArray(MEDMEM_Array_ * Value)
3775   throw (MEDEXCEPTION)
3776 {
3777   if (NULL != _value) delete _value ;
3778   _value=Value ;
3779 }
3780
3781 /*!
3782   \if developper
3783   Return a reference to  the MEDARRAY<T, INTERLACING_TAG> in FIELD.
3784   \endif
3785 */
3786 template <class T, class INTERLACING_TAG>
3787 inline MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() const throw (MEDEXCEPTION)
3788 {
3789   const char* LOC = "MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() : ";
3790   BEGIN_OF_MED(LOC);
3791   END_OF_MED(LOC);
3792   return _value ;
3793 }
3794 template <class T,class INTERLACING_TAG>  inline
3795 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array *
3796 FIELD<T, INTERLACING_TAG>::getArrayGauss() const throw (MEDEXCEPTION)
3797 {
3798   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayGauss() : ";
3799   BEGIN_OF_MED(LOC);
3800
3801   if ( getGaussPresence() )
3802     return static_cast<ArrayGauss *> (_value);
3803   else
3804     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
3805                                  "The field has no Gauss Point"));
3806
3807   END_OF_MED(LOC);
3808
3809 }
3810
3811 template <class T,class INTERLACING_TAG>  inline
3812 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array *
3813 FIELD<T, INTERLACING_TAG>::getArrayNoGauss() const throw (MEDEXCEPTION)
3814 {
3815   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayNoGauss() : ";
3816   BEGIN_OF_MED(LOC);
3817
3818   if ( ! getGaussPresence() )
3819     return static_cast < ArrayNoGauss * > (_value);
3820   else
3821     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
3822                                  "The field has Gauss Point"));
3823
3824   END_OF_MED(LOC);
3825 }
3826
3827
3828 template <class T,class INTERLACING_TAG> inline bool
3829 FIELD<T, INTERLACING_TAG>::getGaussPresence() const throw (MEDEXCEPTION)
3830 {
3831   if (_value != NULL)
3832     return _value->getGaussPresence();
3833   else
3834     throw MEDEXCEPTION("FIELD<T, INTERLACING_TAG>::getGaussPresence() const : Can't call getGaussPresence on a null _value");
3835 }
3836
3837 /*!
3838   Return the actual length of the reference to values array returned by getValue.
3839   Take care of number of components and number of Gauss points by geometric type
3840 */
3841 template <class T, class INTERLACING_TAG>
3842 inline int FIELD<T, INTERLACING_TAG>::getValueLength() const
3843   throw (MEDEXCEPTION)
3844 {
3845   if ( getGaussPresence() )
3846     return static_cast<ArrayGauss *>(_value)->getArraySize() ;
3847   else
3848     return static_cast<ArrayNoGauss *>(_value)->getArraySize() ;
3849 }
3850
3851
3852 /*! \if MEDMEM_ug
3853 \defgroup FIELD_value Field values
3854
3855 These methods are provided for accessing the values of a field.
3856 There are two ways to do so : one consists in using accessors
3857 that retrieve elements or group of elements from the entire field. Typical use is
3858 \verbatim
3859 FIELD field(MED_DRIVER, "result.med","Pressure");
3860 double P0=field.getValueIJ(1,1);
3861 \endverbatim
3862
3863 Another way is to retrieve the pointer to the array that contains the
3864 variable values. In this case, the user should be aware of the interlacing mode
3865 so that no mistakes are made when retrieving the values.
3866
3867 \verbatim
3868 FIELD field(MED_DRIVER, "result.med","Pressure");
3869 double* ptrP=field.getValue();
3870 double P0=ptrP[0];
3871 \endverbatim
3872
3873 @{
3874 \endif
3875 */
3876 /*!
3877   Returns a reference to values array to read them.
3878 */
3879 template <class T, class INTERLACIN_TAG>
3880 inline const T* FIELD<T, INTERLACIN_TAG>::getValue() const throw (MEDEXCEPTION)
3881 {
3882   const char* LOC = "FIELD<T, INTERLACING_TAG>::getValue() : ";
3883   BEGIN_OF_MED(LOC);
3884   if ( getGaussPresence() )
3885     return static_cast<ArrayGauss *>(_value)->getPtr() ;
3886   else
3887     return static_cast<ArrayNoGauss *>(_value)->getPtr() ;
3888 }
3889 /*!
3890   Returns a reference to \f$ i^{th} \f$ row
3891   of FIELD values array.
3892   If a faster accessor is intended you may use getArray() once,
3893   then MEDMEM_Array accessors.
3894   Be careful if field support is not on all elements getRow
3895   use support->getValIndFromGlobalNumber(i).
3896 */
3897 template <class T,class INTERLACING_TAG> inline
3898 const T*
3899 FIELD<T,INTERLACING_TAG>::getRow(int i) const throw (MEDEXCEPTION)
3900 {
3901   const char* LOC; LOC = "FIELD<T,INTERLACING_TAG>::getRow(int i) : ";
3902
3903   int valIndex=-1;
3904   if (_support)
3905     valIndex = _support->getValIndFromGlobalNumber(i);
3906   else
3907     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3908
3909   if ( getGaussPresence() )
3910     return static_cast<ArrayGauss *>(_value)->getRow(valIndex) ;
3911   else
3912     return static_cast<ArrayNoGauss *>(_value)->getRow(valIndex) ;
3913 }
3914
3915 /*!
3916   Returns a reference to $j^{th}$ column
3917   of FIELD values array.
3918 */
3919 template <class T,class INTERLACING_TAG> inline const T*
3920 FIELD<T,INTERLACING_TAG>::getColumn(int j) const throw (MEDEXCEPTION)
3921 {
3922   if ( getGaussPresence() )
3923     return static_cast<ArrayGauss *>(_value)->getColumn(j) ;
3924   else
3925     return static_cast<ArrayNoGauss *>(_value)->getColumn(j) ;
3926 }
3927
3928 /*!
3929   Returns the value of $i^{th}$ element and $j^{th}$ component.
3930   This method only works with fields having no particular Gauss point
3931 definition (i.e., fields having one value per element).
3932  This method makes the retrieval of the value independent from the
3933   interlacing pattern, but it is slower than the complete retrieval
3934   obtained by the \b getValue() method.
3935 */
3936 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJ(int i,int j) const throw (MEDEXCEPTION)
3937 {
3938   const char * LOC = "getValueIJ(..)";
3939   int valIndex=-1;
3940   if (_support)
3941     valIndex = _support->getValIndFromGlobalNumber(i);
3942   else
3943     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3944
3945   if ( getGaussPresence() )
3946     return static_cast<ArrayGauss *>(_value)->getIJ(valIndex,j) ;
3947   else
3948     return static_cast<ArrayNoGauss *>(_value)->getIJ(valIndex,j) ;
3949 }
3950
3951 /*!
3952   Returns the $j^{th}$ component of $k^{th}$ Gauss points of $i^{th}$ value.
3953   This method is compatible with elements having more than one Gauss point.
3954   This method makes the retrieval of the value independent from the
3955   interlacing pattern, but it is slower than the complete retrieval
3956   obtained by the \b getValue() method.
3957 */
3958 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION)
3959 {
3960   const char * LOC = "getValueIJK(..)";
3961   int valIndex=-1;
3962   if (_support)
3963     valIndex = _support->getValIndFromGlobalNumber(i);
3964   else
3965     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3966
3967   if ( getGaussPresence() )
3968     return static_cast<ArrayGauss *>(_value)->getIJK(valIndex,j,k) ;
3969   else
3970     return static_cast<ArrayNoGauss *>(_value)->getIJK(valIndex,j,k) ;
3971 }
3972 /*! \if MEDMEM_ug @} \endif */
3973
3974 /*!
3975   Return number of values of a geomertic type in NoInterlaceByType mode
3976 */
3977 template <class T, class INTERLACIN_TAG>
3978 inline int FIELD<T, INTERLACIN_TAG>::getValueByTypeLength(int t) const throw (MEDEXCEPTION)
3979 {
3980   const char * LOC ="getValueByTypeLength() : ";
3981   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3982     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3983
3984   if ( getGaussPresence() ) {
3985     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
3986     if (  t < 1 || t > array->getNbGeoType() )
3987       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
3988     return array->getLengthOfType( t );
3989   }
3990   else {
3991     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
3992     if (  t < 1 || t > array->getNbGeoType() )
3993       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
3994     return array->getLengthOfType( t );
3995   }
3996 }
3997
3998 /*!
3999   Return a reference to values array to read them.
4000 */
4001 template <class T, class INTERLACIN_TAG>
4002 inline const T* FIELD<T, INTERLACIN_TAG>::getValueByType(int t) const throw (MEDEXCEPTION)
4003 {
4004   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4005     throw MEDEXCEPTION(LOCALIZED("getValueByType() : not MED_NO_INTERLACE_BY_TYPE field" ));
4006
4007   if ( getGaussPresence() ) {
4008     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
4009     return array->getPtr() + array->getIndex( t );
4010   }
4011   else {
4012     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
4013     return array->getPtr() + array->getIndex( t );
4014   }
4015 }
4016
4017 /*!
4018   Return the value of i^{th} element in indicated type t and j^{th} component.
4019 */
4020 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJByType(int i,int j, int t) const throw (MEDEXCEPTION)
4021 {
4022   const char * LOC = "getValueIJByType(..)";
4023   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4024     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4025
4026   if ( getGaussPresence() )
4027     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJByType(i,j,t) ;
4028   else
4029     return static_cast<ArrayNoByType *>(_value)->getIJByType(i,j,t) ;
4030 }
4031
4032 /*!
4033   Return the j^{th} component of k^{th} gauss points of i^{th} value with type t.
4034 */
4035 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION)
4036 {
4037   const char * LOC = "getValueIJKByType(..)";
4038   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4039     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4040
4041   if ( getGaussPresence() )
4042     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJKByType(i,j,k,t) ;
4043   else
4044     return static_cast<ArrayNoByType      *>(_value)->getIJKByType(i,j,k,t) ;
4045 }
4046
4047
4048 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGeometricTypes() const throw (MEDEXCEPTION)
4049 {
4050   const char * LOC = "getNumberOfGeometricTypes(..)";
4051   BEGIN_OF_MED(LOC);
4052   if (_support)
4053     return _support->getNumberOfTypes();
4054   else
4055     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4056   END_OF_MED(LOC);
4057 }
4058
4059 /*! \if MEDMEM_ug
4060 \addtogroup FIELD_gauss
4061 @{
4062 \endif
4063  */
4064
4065 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> &
4066 FIELD<T,INTERLACING_TAG>::getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
4067 {
4068   const char * LOC ="getGaussLocalization(MED_EN::medGeometryElement geomElement) : ";
4069   const GAUSS_LOCALIZATION_ * locPtr=0;
4070
4071   locMap::const_iterator it;
4072   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4073         locPtr = (*it).second;
4074         return  *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
4075   }
4076   else
4077     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
4078
4079 }
4080
4081 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> *
4082 FIELD<T,INTERLACING_TAG>::getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
4083 {
4084   const char * LOC ="getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) : ";
4085   const GAUSS_LOCALIZATION_ * locPtr=0;
4086
4087   locMap::const_iterator it;
4088   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4089         locPtr = (*it).second;
4090         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
4091   }
4092   else
4093     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
4094
4095 }
4096
4097 /*!
4098  * \brief Return GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
4099  */
4100 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION_ *
4101 FIELD<T,INTERLACING_TAG>::getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
4102 {
4103   const char * LOC ="getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) : ";
4104
4105   locMap::const_iterator it;
4106   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4107         return (*it).second;
4108   }
4109   else
4110     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type: "<< geomElement ));
4111
4112 }
4113
4114 /*!
4115  * \brief Take onership of GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
4116  */
4117 template <class T,class INTERLACING_TAG> void
4118 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc)
4119 {
4120   locMap::iterator it = _gaussModel.find(geomElement);
4121   if ( it != _gaussModel.end() ) {
4122     delete it->second;
4123     it->second = gaussloc;
4124   }
4125   else {
4126     _gaussModel[ geomElement ] = gaussloc;
4127   }
4128 }
4129
4130
4131 template <class T,class INTERLACING_TAG> void
4132 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG>& gaussloc)
4133 {
4134   locMap::iterator it = _gaussModel.find(geomElement);
4135   if ( it != _gaussModel.end() ) {
4136     delete it->second;
4137     it->second = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
4138   }
4139   else {
4140     _gaussModel[ geomElement ] = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
4141   }
4142 }
4143
4144 /*!
4145   Returns number of Gauss points for this medGeometryElement.
4146
4147   Note :
4148   if there is no GAUSS_LOCALIZATION having this medGeometryElement but
4149   the medGeometryElement exist in the SUPPORT, getNumberOfGaussPoints
4150   return 1
4151 */
4152 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) const
4153   throw (MEDEXCEPTION)
4154 {
4155   const char * LOC ="getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
4156   const GAUSS_LOCALIZATION_ * locPtr=0;
4157
4158   locMap::const_iterator it;
4159   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4160         locPtr = (*it).second;
4161         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr)->getNbGauss();
4162   }
4163   else
4164     if (_support)
4165       try {
4166         if ( _support->getNumberOfElements(geomElement) ) return 1;
4167       } catch ( MEDEXCEPTION & ex) {
4168         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<< "GeometricType not found !" )) ;
4169       }
4170     else
4171       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4172
4173   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Should never execute this!" ));
4174
4175 }
4176
4177 /*!
4178   Returns number of Gauss points for each geometric type.
4179
4180   Note :
4181   if there is no gauss points whatever the geometric type is
4182   it returns an exception. (renvoyer un tableau de 1 ?)
4183 */
4184 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints() const
4185   throw (MEDEXCEPTION)
4186 {
4187   const char * LOC ="const int * getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
4188
4189   if (_value)
4190     if ( getGaussPresence() ) {
4191       return static_cast<ArrayGauss *>(_value)->getNbGaussGeo()+1;
4192     } else
4193       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"value hasn't Gauss points " ));
4194
4195     else
4196       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Value not defined" ));
4197 }
4198
4199 /*!
4200   Returns number of Gauss points for element n°i.
4201   The i index is a global index (take care of previous element
4202   on different geometric type).
4203 */
4204 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNbGaussI(int i) const throw (MEDEXCEPTION)
4205 {
4206   const char * LOC = "getNbGaussI(..)";
4207
4208   int valIndex=-1;
4209   if (_support)
4210     valIndex = _support->getValIndFromGlobalNumber(i);
4211   else
4212     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4213
4214   if (_value)
4215    if ( getGaussPresence() )
4216      return static_cast<ArrayGauss *>(_value)->getNbGauss(valIndex) ;
4217    else
4218      return static_cast<ArrayNoGauss *>(_value)->getNbGauss(valIndex) ;
4219  else
4220    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_value not defined" ));
4221 }
4222 /*!
4223 @}
4224  */
4225
4226
4227 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfElements() const throw (MEDEXCEPTION)
4228 {
4229   const char * LOC = "getNumberOfElements(..)";
4230   BEGIN_OF_MED(LOC);
4231   if (_support)
4232     return _support->getNumberOfElements();
4233   else
4234     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4235   END_OF_MED(LOC);
4236 }
4237
4238 template <class T,class INTERLACING_TAG> const MED_EN::medGeometryElement  * FIELD<T,INTERLACING_TAG>::getGeometricTypes()  const throw (MEDEXCEPTION)
4239 {
4240   const char * LOC = "getGeometricTypes(..)";
4241   BEGIN_OF_MED(LOC);
4242   if (_support)
4243     return _support->getTypes();
4244   else
4245     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4246   END_OF_MED(LOC);
4247 }
4248 template <class T,class INTERLACING_TAG> bool  FIELD<T,INTERLACING_TAG>::isOnAllElements() const throw (MEDEXCEPTION)
4249 {
4250   const char * LOC = "isOnAllElements(..)";
4251   BEGIN_OF_MED(LOC);
4252   if (_support)
4253     return _support->isOnAllElements();
4254   else
4255     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4256   END_OF_MED(LOC);
4257 }
4258
4259
4260 /*! \if MEDMEM_ug
4261 \addtogroup FIELD_value
4262 @{
4263 \endif */
4264
4265 /*!
4266   Copy new values array in FIELD according to the given mode.
4267
4268   Array must have right size. If not results are unpredicable.
4269   In MED_FULL_INTERLACE mode, values are stored elementwise in X1,Y1,Z1,X2,Y2,Z2.. order.
4270 In MED_NO_INTERLACE mode, values are stored componentwise in X1,X2,X3,...,Y1,Y2,Y3,... order.
4271 */
4272 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValue( T* value) throw (MEDEXCEPTION)
4273 {
4274   if ( getGaussPresence() )
4275     static_cast<ArrayGauss *>(_value)->setPtr(value) ;
4276   else
4277     static_cast<ArrayNoGauss *>(_value)->setPtr(value) ;
4278 }
4279
4280 /*!
4281   Update values array in the j^{th} row of FIELD values array with the given ones and
4282   according to specified mode.
4283 */
4284 template <class T,class INTERLACING_TAG>
4285 inline void FIELD<T,INTERLACING_TAG>::setRow( int i, T* value) throw (MEDEXCEPTION)
4286 {
4287   const char * LOC = "FIELD<T,INTERLACING_TAG>::setRow(int i, T* value) : ";
4288   int valIndex=i;
4289   if (_support)
4290     valIndex = _support->getValIndFromGlobalNumber(i);
4291   else
4292     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
4293
4294   if ( getGaussPresence() )
4295     static_cast<ArrayGauss *>(_value)->setRow(valIndex, value) ;
4296   else
4297     static_cast<ArrayNoGauss *>(_value)->setRow(valIndex, value) ;
4298 }
4299
4300 /*!
4301   Update values array in the $j^{th}$ column of FIELD values array with the given ones and
4302   according to specified mode.
4303 */
4304 template <class T,class INTERLACING_TAG>
4305 inline void FIELD<T,INTERLACING_TAG>::setColumn( int j, T* value) throw (MEDEXCEPTION)
4306 {
4307   if ( getGaussPresence() )
4308     static_cast<ArrayGauss *>(_value)->setColumn(j, value) ;
4309   else
4310     static_cast<ArrayNoGauss *>(_value)->setColumn(j, value) ;
4311 }
4312
4313 /*!
4314   Sets the value of i^{th} element and j^{th} component with the given one.
4315 */
4316 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) throw (MEDEXCEPTION)
4317 {
4318   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
4319   int valIndex=-1;
4320   if (_support)
4321     valIndex = _support->getValIndFromGlobalNumber(i);
4322   else
4323     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
4324
4325   if ( getGaussPresence() )
4326     static_cast<ArrayGauss *>(_value)->setIJ(valIndex,j,value) ;
4327   else
4328     static_cast<ArrayNoGauss *>(_value)->setIJ(valIndex,j,value) ;
4329 }
4330
4331 /*!
4332   Set the value of i^{th} element, j^{th} component and k^{th} gauss point with the given one.
4333 */
4334 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION)
4335 {
4336   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
4337   int valIndex=-1;
4338   if (_support)
4339     valIndex = _support->getValIndFromGlobalNumber(i);
4340   else
4341     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
4342
4343   if ( getGaussPresence() )
4344     static_cast<ArrayGauss *>(_value)->setIJK(valIndex,j,k,value) ;
4345   else
4346     static_cast<ArrayNoGauss *>(_value)->setIJK(valIndex,j,k,value) ;
4347 }
4348
4349 /*!
4350   Set the value of i^{th} element and j^{th} component with the given one.
4351 */
4352 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION)
4353 {
4354   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) : ";
4355   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4356     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4357
4358   if ( getGaussPresence() )
4359     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJByType(i,j,t,value) ;
4360   else
4361     return static_cast<ArrayNoByType *>(_value)->setIJByType(i,j,t,value) ;
4362 }
4363
4364 /*!
4365   Set the value of component of k^{th} gauss points of i^{th} element and j^{th} component with the given one.
4366 */
4367 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)
4368 {
4369   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int t, int k, T value) : ";
4370   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4371     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4372
4373   if ( getGaussPresence() )
4374     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJKByType(i,j,k,t,value) ;
4375   else
4376     return static_cast<ArrayNoByType *>(_value)->setIJKByType(i,j,k,t,value) ;
4377 }
4378 /*! \if MEDMEM_ug @} \endif */
4379
4380 /*
4381   METHODS
4382 */
4383
4384 /*!
4385   Fill array by using T_Analytic.
4386   WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
4387   Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
4388  */
4389 template <class T, class INTERLACING_TAG>
4390 void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION)
4391 {
4392   const char * LOC = "void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) : ";
4393   int i,j;
4394   if (_support == (SUPPORT *) NULL)
4395       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No Support defined."));
4396
4397   const GMESH * mesh = _support->getMesh();
4398   int spaceDim = mesh->getSpaceDimension();
4399   const double * coord;
4400
4401   const double * bary;
4402   FIELD<double,FullInterlace> * barycenterField=0;
4403
4404   double ** xyz=new double* [spaceDim];
4405   bool deallocateXyz=false;
4406   if(_support->getEntity()==MED_EN::MED_NODE)
4407     {
4408       const MESH * unstructured = _support->getMesh()->convertInMESH();
4409       if (_support->isOnAllElements())
4410         {
4411           coord=unstructured->getCoordinates(MED_EN::MED_NO_INTERLACE);
4412           for(i=0; i<spaceDim; i++)
4413             xyz[i]=(double *)coord+i*_numberOfValues;
4414         }
4415       else
4416         {
4417           coord = unstructured->getCoordinates(MED_EN::MED_FULL_INTERLACE);
4418           const int * nodesNumber=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
4419           for(i=0; i<spaceDim; i++)
4420             xyz[i]=new double[_numberOfValues];
4421           deallocateXyz=true;
4422           for(i=0;i<_numberOfValues;i++)
4423             {
4424               for(j=0;j<spaceDim;j++)
4425                 xyz[j][i]=coord[(nodesNumber[i]-1)*spaceDim+j];
4426             }
4427         }
4428       unstructured->removeReference();
4429     }
4430   else
4431     {
4432       barycenterField = mesh->getBarycenter(_support);
4433       bary = barycenterField->getValue();
4434       for(i=0; i<spaceDim; i++)
4435         xyz[i]=new double[_numberOfValues];
4436       deallocateXyz=true;
4437       for(i=0;i<_numberOfValues;i++) {
4438         for(j=0;j<spaceDim;j++)
4439           xyz[j][i]=bary[i*spaceDim+j];
4440       }
4441     }
4442   T* valsToSet=(T*)getValue();
4443   double *temp=new double[spaceDim];
4444   for(i=0;i<_numberOfValues;i++)
4445   {
4446     for(j=0;j<spaceDim;j++)
4447       temp[j]=xyz[j][i];
4448     f(temp,valsToSet+i*_numberOfComponents);
4449   }
4450   delete [] temp;
4451   if(barycenterField)
4452     delete barycenterField;
4453   if(deallocateXyz)
4454     for(j=0;j<spaceDim;j++)
4455       delete [] xyz[j];
4456   delete [] xyz;
4457 }
4458
4459 /*!
4460   Execute a function on _values on 'this' and put the result on a newly created field that has to be deallocated.
4461   WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
4462   Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
4463  */
4464 template <class T, class INTERLACING_TAG>
4465 FIELD<T,INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::execFunc(int nbOfComponents,
4466                                                               myFuncType2 f) throw (MEDEXCEPTION)
4467 {
4468   FIELD<T,INTERLACING_TAG> *ret=new FIELD<T,INTERLACING_TAG>(_support,nbOfComponents);
4469   const T* valsInput=getValue();
4470   T* valsOutPut=(T*)ret->getValue();
4471   int i;
4472   for(i=0;i<_numberOfValues;i++)
4473     f(valsInput+i*_numberOfComponents,valsOutPut+i*nbOfComponents);
4474   return ret;
4475 }
4476
4477 }//End namespace MEDMEM
4478
4479 #endif /* FIELD_HXX */