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