Back to index

salome-med  6.5.0
MEDMEM_Field.hxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 /*
00024  File Field.hxx
00025 */
00026 
00027 #ifndef FIELD_HXX
00028 #define FIELD_HXX
00029 
00030 #include "MEDMEM.hxx"
00031 
00032 #include "MEDMEM_Utilities.hxx"
00033 #include "MEDMEM_MedVersion.hxx"
00034 #include "MEDMEM_Exception.hxx"
00035 #include "MEDMEM_define.hxx"
00036 #include "MEDMEM_Support.hxx"
00037 #include "MEDMEM_Unit.hxx"
00038 #include "MEDMEM_nArray.hxx"
00039 #include "MEDMEM_GenDriver.hxx"
00040 #include "MEDMEM_RCBase.hxx"
00041 #include "MEDMEM_ArrayInterface.hxx"
00042 #include "MEDMEM_SetInterlacingType.hxx"
00043 #include "MEDMEM_FieldForward.hxx"
00044 #include "MEDMEM_GaussLocalization.hxx"
00045 #include "InterpKernelGaussCoords.hxx"
00046 #include "PointLocator.hxx"
00047 
00048 #include <vector>
00049 #include <map>
00050 #include <algorithm>
00051 #include <memory>
00052 #include <math.h>
00053 #include <cmath>
00054 #include <float.h>
00055 
00056 namespace MEDMEM {
00057 
00058   template<class T>
00059   struct MinMax {
00060   };
00061 
00062   template<>
00063   struct MinMax<double> {
00064     static double getMin() { return DBL_MIN; }
00065     static double getMax() { return DBL_MAX; }
00066   };
00067 
00068   template<>
00069   struct MinMax<int> {
00070     static int getMin() { return INT_MIN; }
00071     static int getMax() { return INT_MAX; }
00072   };  
00073 
00074   template < typename T > struct SET_VALUE_TYPE {
00075     static const MED_EN::med_type_champ _valueType = MED_EN::MED_UNDEFINED_TYPE;};
00076   template < > struct SET_VALUE_TYPE<double> {
00077     static const MED_EN::med_type_champ _valueType = MED_EN::MED_REEL64; };
00078   template < > struct SET_VALUE_TYPE<int> {
00079     static const MED_EN::med_type_champ _valueType = MED_EN::MED_INT32; };
00080 
00206   class MEDMEM_EXPORT FIELD_ : public RCBASE    // GENERIC POINTER TO a template <class T, class INTERLACING_TAG> class FIELD
00207 {
00208 protected:
00209 
00210   bool            _isRead ;
00211   bool            _isMinMax;
00212 
00218   string          _name ;
00224   string          _description ;
00230   const SUPPORT * _support ;
00231 
00237   int             _numberOfComponents ;
00245   int             _numberOfValues ;
00246 
00262   vector<int>     _componentsTypes ;
00269   vector<string>  _componentsNames;
00276   vector<string>  _componentsDescriptions;
00283   vector<UNIT>    _componentsUnits;
00290   vector<string>  _MEDComponentsUnits;
00296   int             _iterationNumber ;
00302   double          _time;
00308   int             _orderNumber ;
00316   MED_EN::med_type_champ _valueType ;
00324   MED_EN::medModeSwitch _interlacingType;
00325 
00326   vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
00327   static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
00328   static void _deepCheckFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
00329   void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL,
00330                                const bool           nodalAllowed = false) const  throw (MEDEXCEPTION);
00331   FIELD<double>* _getFieldSize(const SUPPORT *subSupport=NULL) const;
00332  public:
00336   FIELD_ ();
00341   FIELD_(const SUPPORT * Support, const int NumberOfComponents);
00346   FIELD_(const FIELD_ &m);
00347 
00351   virtual ~FIELD_();
00352 
00353 public:
00354 
00355   friend class MED_MED_RDONLY_DRIVER21;
00356   friend class MED_MED_WRONLY_DRIVER21;
00357   friend class MED_MED_RDWR_DRIVER21;
00358   friend class MED_MED_RDONLY_DRIVER22;
00359   friend class MED_MED_WRONLY_DRIVER22;
00360   friend class MED_MED_RDWR_DRIVER22;
00361   friend class VTK_MED_DRIVER;
00362 
00363  FIELD_& operator=(const FIELD_ &m);
00364 
00365   virtual  void     rmDriver(int index=0);
00366 
00380   virtual   int     addDriver(driverTypes driverType,
00381                               const string & fileName="Default File Name.med",
00382                               const string & driverFieldName="Default Field Nam",
00383                               MED_EN::med_mode_acces access=MED_EN::RDWR) ;
00384 
00385   virtual  int      addDriver( GENDRIVER & driver);
00386   virtual  void     read (driverTypes driverType, const std::string & fileName);
00387   virtual  void     read (const GENDRIVER &);
00388   virtual  void     read(int index=0);
00389   virtual  void     openAppend( void );
00390   virtual  void     write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
00391   virtual  void     write(driverTypes driverType,
00392                           const std::string & fileName,
00393                           MED_EN::med_mode_acces medMode=MED_EN::RDWR);
00394 
00397   virtual  void     write(int index=0);
00400   virtual  void     writeAppend(const GENDRIVER &);
00401   virtual  void     writeAppend(int index=0, const string & driverName="");
00402 
00403   inline void     setName(const string Name);
00404   inline string   getName() const;
00405   inline void     setDescription(const string Description);
00406   inline string   getDescription() const;
00407   inline const SUPPORT * getSupport() const;
00408   inline void     setSupport(const SUPPORT * support);
00409   inline void     setNumberOfComponents(const int NumberOfComponents);
00410   inline int      getNumberOfComponents() const;
00411   inline void     setNumberOfValues(const int NumberOfValues);
00412   inline int      getNumberOfValues() const;
00413   inline void     setComponentsNames(const string * ComponentsNames);
00414   inline void     setComponentName(int i, const string ComponentName);
00415   inline const string * getComponentsNames() const;
00416   inline string   getComponentName(int i) const;
00417   inline void     setComponentsDescriptions(const string * ComponentsDescriptions);
00418   inline void     setComponentDescription(int i, const string ComponentDescription);
00419   inline const string * getComponentsDescriptions() const;
00420   inline string   getComponentDescription(int i) const;
00421 
00422   // provisoire : en attendant de regler le probleme des unites !
00423   inline void     setComponentsUnits(const UNIT * ComponentsUnits);
00424   inline const UNIT *   getComponentsUnits() const;
00425   inline const UNIT *   getComponentUnit(int i) const;
00426   inline void     setMEDComponentsUnits(const string * MEDComponentsUnits);
00427   inline void     setMEDComponentUnit(int i, const string MEDComponentUnit);
00428   inline const string * getMEDComponentsUnits() const;
00429   inline string   getMEDComponentUnit(int i) const;
00430 
00431   inline void     setIterationNumber(int IterationNumber);
00432   inline int      getIterationNumber() const;
00433   inline void     setTime(double Time);
00434   inline double   getTime() const;
00435   inline void     setOrderNumber(int OrderNumber);
00436   inline int      getOrderNumber() const;
00437 
00438   inline MED_EN::med_type_champ getValueType () const;
00439   inline MED_EN::medModeSwitch  getInterlacingType() const;
00440   virtual inline bool getGaussPresence() const throw (MEDEXCEPTION);
00441 protected:
00442   void copyGlobalInfo(const FIELD_& m);
00443 };
00444 
00445 // -----------------
00446 // Methodes Inline
00447 // -----------------
00457 inline void FIELD_::setName(const string Name)
00458 {
00459   _name=Name;
00460 }
00464 inline string FIELD_::getName() const
00465 {
00466   return _name;
00467 }
00471 inline void FIELD_::setDescription(const string Description)
00472 {
00473   _description=Description;
00474 }
00478 inline string FIELD_::getDescription() const
00479 {
00480   return _description;
00481 }
00485 inline void FIELD_::setNumberOfComponents(const int NumberOfComponents)
00486 {
00487   _numberOfComponents=NumberOfComponents;
00488   _componentsTypes.resize(_numberOfComponents);
00489   _componentsNames.resize(_numberOfComponents);
00490   _componentsDescriptions.resize(_numberOfComponents);
00491   _componentsUnits.resize(_numberOfComponents);
00492   _MEDComponentsUnits.resize(_numberOfComponents);
00493 }
00497 inline int FIELD_::getNumberOfComponents() const
00498 {
00499   return _numberOfComponents ;
00500 }
00506 inline void FIELD_::setNumberOfValues(const int NumberOfValues)
00507 {
00508   _numberOfValues=NumberOfValues;
00509 }
00513 inline int FIELD_::getNumberOfValues() const
00514 {
00515   return _numberOfValues ;
00516 }
00517 
00524 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
00525 {
00526   _componentsNames.resize(_numberOfComponents);
00527   for (int i=0; i<_numberOfComponents; i++)
00528     _componentsNames[i]=ComponentsNames[i] ;
00529 }
00536 inline void FIELD_::setComponentName(int i, const string ComponentName)
00537 {
00538   const char * LOC = " FIELD_::setComponentName() : ";
00539   BEGIN_OF_MED(LOC);
00540   if( i<1 || i>_numberOfComponents )
00541     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
00542 
00543   _componentsNames[i-1]=ComponentName ;
00544 }
00550 inline const string * FIELD_::getComponentsNames() const
00551 {
00552   return &(_componentsNames[0]) ;
00553 }
00558 inline string FIELD_::getComponentName(int i) const
00559 {
00560   const char * LOC = " FIELD_::getComponentName() : ";
00561   BEGIN_OF_MED(LOC);
00562   if( i<1 || i>_numberOfComponents )
00563     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
00564 
00565   return _componentsNames[i-1] ;
00566 }
00574 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
00575 {
00576   _componentsDescriptions.resize(_numberOfComponents);
00577   for (int i=0; i<_numberOfComponents; i++)
00578     _componentsDescriptions[i]=ComponentsDescriptions[i] ;
00579 }
00586 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
00587 {
00588   const char * LOC = " FIELD_::setComponentDescription() : ";
00589   BEGIN_OF_MED(LOC);
00590   if( i<1 || i>_numberOfComponents )
00591     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
00592 
00593   _componentsDescriptions[i-1]=ComponentDescription ;
00594 }
00600 inline const string * FIELD_::getComponentsDescriptions() const
00601 {
00602   return &(_componentsDescriptions[0]);
00603 }
00608 inline string FIELD_::getComponentDescription(int i) const
00609 {
00610   const char * LOC = " FIELD_::setComponentDescription() : ";
00611   BEGIN_OF_MED(LOC);
00612   if( i<1 || i>_numberOfComponents )
00613     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
00614 
00615   return _componentsDescriptions[i-1];
00616 }
00617 
00627 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
00628 {
00629   _componentsUnits.resize(_numberOfComponents);
00630   for (int i=0; i<_numberOfComponents; i++)
00631     _componentsUnits[i]=ComponentsUnits[i] ;
00632 }
00639 inline const UNIT * FIELD_::getComponentsUnits() const
00640 {
00641   return &(_componentsUnits[0]);
00642 }
00647 inline const UNIT * FIELD_::getComponentUnit(int i) const
00648 {
00649   const char * LOC = " FIELD_::getComponentUnit() : ";
00650   BEGIN_OF_MED(LOC);
00651   if( i<1 || i>_numberOfComponents )
00652     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
00653 
00654   return &_componentsUnits[i-1] ;
00655 }
00664 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
00665 {
00666   _MEDComponentsUnits.resize(_numberOfComponents);
00667   for (int i=0; i<_numberOfComponents; i++)
00668     _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
00669 }
00676 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
00677 {
00678   const char * LOC = " FIELD_::setMEDComponentUnit() : ";
00679   BEGIN_OF_MED(LOC);
00680   if( i<1 || i>_numberOfComponents )
00681     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
00682 
00683   _MEDComponentsUnits[i-1]=MEDComponentUnit ;
00684 }
00690 inline const string * FIELD_::getMEDComponentsUnits() const
00691 {
00692   return &(_MEDComponentsUnits[0]);
00693 }
00698 inline string FIELD_::getMEDComponentUnit(int i) const
00699 {
00700   const char * LOC = " FIELD_::getMEDComponentUnit() : ";
00701   BEGIN_OF_MED(LOC);
00702   if( i<1 || i>_numberOfComponents )
00703     throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
00704 
00705   return _MEDComponentsUnits[i-1] ;
00706 }
00710 inline void FIELD_::setIterationNumber(int IterationNumber)
00711 {
00712   _iterationNumber=IterationNumber;
00713 }
00717 inline int FIELD_::getIterationNumber() const
00718 {
00719   return _iterationNumber ;
00720 }
00724 inline void FIELD_::setTime(double Time)
00725 {
00726   _time=Time ;
00727 }
00731 inline double FIELD_::getTime() const
00732 {
00733   return _time ;
00734 }
00740 inline void FIELD_::setOrderNumber(int OrderNumber)
00741 {
00742   _orderNumber=OrderNumber ;
00743 }
00747 inline int FIELD_::getOrderNumber() const
00748 {
00749   return _orderNumber ;
00750 }
00754 inline  const SUPPORT * FIELD_::getSupport() const
00755 {
00756   return _support ;
00757 }
00763 inline void FIELD_::setSupport(const SUPPORT * support)
00764 {
00765   //A.G. Addings for RC
00766   if(_support!=support)
00767     {
00768       if(_support)
00769         _support->removeReference();
00770       _support = support ;
00771       if(_support)
00772         _support->addReference();
00773     }
00774 }
00778 inline MED_EN::med_type_champ FIELD_::getValueType () const
00779 {
00780   return _valueType ;
00781 }
00782 
00786   inline MED_EN::medModeSwitch FIELD_::getInterlacingType () const
00787 {
00788   return _interlacingType ;
00789 }
00800   inline bool  FIELD_::getGaussPresence() const throw (MEDEXCEPTION)
00801 {
00802   const char * LOC = "FIELD_::getGaussPresence() : ";
00803   throw MEDEXCEPTION(STRING(LOC) << " This FIELD_ doesn't rely on a FIELD<T>" );
00804 }
00805 
00808 } //End namespace MEDMEM
00809 
00811 // END OF CLASS FIELD_ //
00813 
00823 namespace MEDMEM {
00824 
00825   template<class T2> class MED_FIELD_RDONLY_DRIVER;
00826   template<class T2> class MED_FIELD_WRONLY_DRIVER;
00827   template<class T2> class VTK_FIELD_DRIVER;
00828 
00829 
00830   template <class T,
00831             class INTERLACING_TAG
00832             > class FIELD : public FIELD_
00833 {
00834 protected:
00835 
00836   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array   ArrayNoGauss;
00837   typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array     ArrayGauss;
00838   typedef typename MEDMEM_ArrayInterface<T,NoInterlace,NoGauss>::Array       ArrayNo;
00839   typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array     ArrayFull;
00840   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,NoGauss>::Array ArrayNoByType;
00841   typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,Gauss>::Array   ArrayNoByTypeGauss;
00842   typedef MEDMEM_Array_ Array;
00843   typedef T ElementType;
00844   typedef INTERLACING_TAG InterlacingTag;
00845   typedef map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> locMap;
00846 
00847   // array of value of type T
00848   Array *_value ;
00849 
00850   // MESH, to be used for field reading from a file (if desired to link
00851   // to existing support instead of new support creation for the field)
00852   GMESH* _mesh;
00853 
00854   // extrema values
00855   T _vmin;
00856   T _vmax;
00857 
00858   map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> _gaussModel; //A changer quand les drivers seront template de l'entrelacement
00859 
00860   static T _scalarForPow;
00861   static T pow(T x);
00862 
00863 private:
00864   void _operation(const FIELD& m,const FIELD& n, const char* Op);
00865   void _operationInitialize(const FIELD& m,const FIELD& n, const char* Op);
00866   void _add_in_place(const FIELD& m,const FIELD& n);
00867   void _sub_in_place(const FIELD& m,const FIELD& n);
00868   void _mul_in_place(const FIELD& m,const FIELD& n);
00869   void _div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION);
00870 public:
00871   FIELD();
00872   FIELD(const FIELD &m);
00873   FIELD(const SUPPORT * Support, const int NumberOfComponents) throw (MEDEXCEPTION);
00874   FIELD(driverTypes driverType,
00875         const string & fileName, const string & fieldDriverName,
00876         const int iterationNumber=-1, const int orderNumber=-1,
00877         GMESH* mesh = 0)
00878     throw (MEDEXCEPTION);
00879   FIELD(const SUPPORT * Support, driverTypes driverType,
00880         const string & fileName="", const string & fieldName="",
00881         const int iterationNumber = -1, const int orderNumber = -1)
00882     throw (MEDEXCEPTION);
00883   ~FIELD();
00884 
00885 public:
00886   FIELD & operator=(const FIELD &m);
00887         FIELD & operator=(T value);
00888   FIELD *operator+(const FIELD& m) const;
00889   FIELD *operator-(const FIELD& m) const;
00890   FIELD *operator*(const FIELD& m) const;
00891   FIELD *operator/(const FIELD& m) const;
00892   FIELD *operator-() const;
00893   FIELD& operator+=(const FIELD& m);
00894   FIELD& operator-=(const FIELD& m);
00895   FIELD& operator*=(const FIELD& m);
00896   FIELD& operator/=(const FIELD& m);
00897 
00898   void          applyLin(T a, T b, int icomp);
00899   static FIELD* add(const FIELD& m, const FIELD& n);
00900   static FIELD* addDeep(const FIELD& m, const FIELD& n);
00901   static FIELD* sub(const FIELD& m, const FIELD& n);
00902   static FIELD* subDeep(const FIELD& m, const FIELD& n);
00903   static FIELD* mul(const FIELD& m, const FIELD& n);
00904   static FIELD* mulDeep(const FIELD& m, const FIELD& n);
00905   static FIELD* div(const FIELD& m, const FIELD& n);
00906   static FIELD* divDeep(const FIELD& m, const FIELD& n);
00907   double normMax() const throw (MEDEXCEPTION);
00908 
00909   //------- TDG and BS addings
00910 
00911   void getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION);
00912   vector<int> getHistogram(int &nbint) throw (MEDEXCEPTION);
00913   FIELD<double>* buildGradient() const throw (MEDEXCEPTION);
00914   FIELD<double>* buildNorm2Field() const throw (MEDEXCEPTION);
00915 
00916   //-------------------
00917 
00918   double norm2() const throw (MEDEXCEPTION);
00919   void   applyLin(T a, T b);
00920   template <T T_function(T)> void applyFunc();
00921   void applyPow(T scalar);
00922   static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
00923   double normL2(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
00924   double normL2(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
00925   double normL1(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
00926   double normL1(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
00927   double integral(const SUPPORT *subSupport=NULL) const throw (MEDEXCEPTION);
00928   FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
00929 
00930   friend class MED_FIELD_RDONLY_DRIVER<T>;
00931   friend class MED_FIELD_WRONLY_DRIVER<T>;
00932   friend class VTK_FIELD_DRIVER<T>;
00933 
00934   void init ();
00935   void rmDriver(int index=0);
00936   int  addDriver(driverTypes driverType,
00937                  const string & fileName="Default File Name.med",
00938                  const string & driverFieldName="Default Field Name",
00939                  MED_EN::med_mode_acces access=MED_EN::RDWR) ;
00940 
00941   int  addDriver(GENDRIVER & driver);
00942 
00943   void allocValue(const int NumberOfComponents);
00944   void allocValue(const int NumberOfComponents, const int LengthValue);
00945 
00946   void deallocValue();
00947 
00948   inline void read(int index=0);
00949   inline void read(const GENDRIVER & genDriver);
00950   inline void read(driverTypes driverType, const std::string& filename);
00951   inline void write(int index=0);
00952   inline void write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
00953   inline void write(driverTypes        driverType,
00954                     const std::string& filename,
00955                     MED_EN::med_mode_acces medMode=MED_EN::RDWR);
00956 
00957   inline void writeAppend(int index=0, const string & driverName = "");
00958   inline void writeAppend(const GENDRIVER &);
00959 
00960   inline MEDMEM_Array_  * getArray()        const throw (MEDEXCEPTION);
00961   inline ArrayGauss     * getArrayGauss()   const throw (MEDEXCEPTION);
00962   inline ArrayNoGauss   * getArrayNoGauss() const throw (MEDEXCEPTION);
00963   inline bool            getGaussPresence() const throw (MEDEXCEPTION);
00964 
00965   inline int          getValueLength() const throw (MEDEXCEPTION);
00966 
00972   inline const T*     getValue()       const throw (MEDEXCEPTION);
00973   inline const T*     getRow(int i)    const throw (MEDEXCEPTION);
00974   inline const T*     getColumn(int j) const throw (MEDEXCEPTION);
00984   inline T            getValueIJ(int i,int j) const throw (MEDEXCEPTION);
00985 
00993   inline T            getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION);
00994 
00995   inline int          getValueByTypeLength(int t)                const throw (MEDEXCEPTION);
00996   inline const T*     getValueByType(int t)                      const throw (MEDEXCEPTION);
00997   inline T            getValueIJByType(int i,int j,int t)        const throw (MEDEXCEPTION);
00998   inline T            getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION);
00999 
01007   bool                getValueOnElement(int eltIdInSup,T* retValues) const throw (MEDEXCEPTION);
01008   void                getValueOnPoint(const double* coords, double* output) const throw (MEDEXCEPTION);
01009   void                getValueOnPoints(int nb_points, const double* coords, double* output) const throw (MEDEXCEPTION);
01010 
01011   const int   getNumberOfGeometricTypes() const throw (MEDEXCEPTION);
01012   const GAUSS_LOCALIZATION<INTERLACING_TAG> & getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
01013   const GAUSS_LOCALIZATION<INTERLACING_TAG> * getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
01014   const GAUSS_LOCALIZATION_* getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
01015   void setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG> & gaussloc);
01016   void setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc);
01017   const int * getNumberOfGaussPoints() const throw (MEDEXCEPTION);
01018   const int   getNumberOfGaussPoints( MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
01019   const int   getNbGaussI(int i)          const throw (MEDEXCEPTION);
01020   const int * getNumberOfElements()       const throw (MEDEXCEPTION);
01021   const MED_EN::medGeometryElement  * getGeometricTypes()  const throw (MEDEXCEPTION);
01022   bool        isOnAllElements()           const throw (MEDEXCEPTION);
01023  
01024   inline void setArray(MEDMEM_Array_ *value) throw (MEDEXCEPTION);
01025 
01026   FIELD<double, FullInterlace>* getGaussPointsCoordinates() const throw (MEDEXCEPTION);
01027 
01040   inline void setValue( T* value) throw (MEDEXCEPTION);
01041   inline void setRow( int i, T* value) throw (MEDEXCEPTION);
01042   inline void setColumn( int i, T* value) throw (MEDEXCEPTION);
01046   inline void setValueIJ(int i, int j, T value) throw (MEDEXCEPTION);
01048   inline void setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION);
01049   inline void setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION);
01050   inline void setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION);
01051 
01052   typedef void (*myFuncType)(const double *,T*);
01053   void fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION);
01054   typedef void (*myFuncType2)(const T *,T*);
01055   FIELD<T,INTERLACING_TAG> *execFunc(int nbOfComponents, myFuncType2 f) throw (MEDEXCEPTION);
01056 };
01057 }
01058 
01059 #include "MEDMEM_DriverFactory.hxx"
01060 
01061 namespace MEDMEM {
01062 
01063 template <class T,class INTERLACING_TAG> T FIELD<T, INTERLACING_TAG>::_scalarForPow=1;
01064 
01065 // --------------------
01066 // Implemented Methods
01067 // --------------------
01068 
01072 template <class T, class INTERLACING_TAG>
01073 FIELD<T, INTERLACING_TAG>::FIELD():FIELD_()
01074 {
01075   MESSAGE_MED("Constructeur FIELD sans parametre");
01076 
01077   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
01078   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE);
01079   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
01080 
01081   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
01082   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE);
01083   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
01084 
01085   _value = ( ArrayNoGauss * ) NULL;
01086 
01087   _mesh  = ( MESH* ) NULL;
01088 }
01089 
01112 template <class T, class INTERLACING_TAG>
01113 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
01114                                  const int NumberOfComponents) throw (MEDEXCEPTION) :
01115   FIELD_(Support, NumberOfComponents),_value(NULL)
01116 {
01117   const char* LOC = "FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)";
01118   BEGIN_OF_MED(LOC);
01119   SCRUTE_MED(this);
01120 
01121   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
01122   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
01123     FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
01124 
01125   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
01126   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
01127     FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
01128 
01129   try
01130     {
01131       // becarefull about the numbre of gauss point
01132       _numberOfValues = Support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
01133     }
01134 #if defined(_DEBUG_) || defined(_DEBUG)
01135   catch (MEDEXCEPTION &ex)
01136 #else
01137   catch (MEDEXCEPTION )
01138 #endif
01139     {
01140       MESSAGE_MED("No value defined ! ("<<ex.what()<<")");
01141     }
01142   MESSAGE_MED("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
01143   if ( _numberOfValues > 0 )
01144     {
01145       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
01146         {
01147           const int * nbelgeo = Support->getNumberOfElements();
01148           vector<int> nbelgeoc( Support->getNumberOfTypes() + 1 );
01149           nbelgeoc[0] = 0;
01150           for ( int t = 1; t < (int)nbelgeoc.size(); ++t )
01151             nbelgeoc[t] = nbelgeoc[t-1] + nbelgeo[t-1];
01152           _value = new ArrayNoByType (_numberOfComponents,_numberOfValues,
01153                                       Support->getNumberOfTypes(), &nbelgeoc[0]);
01154         }
01155       else
01156         {
01157           _value = new ArrayNoGauss (_numberOfComponents,_numberOfValues);
01158         }
01159       _isRead = true ;
01160     }
01161   _mesh  = ( MESH* ) NULL;
01162 
01163   END_OF_MED(LOC);
01164 }
01172 template <class T, class INTERLACING_TAG> void FIELD<T, INTERLACING_TAG>::init ()
01173 {
01174 }
01175 
01179 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::FIELD(const FIELD & m):
01180   FIELD_(m)
01181 {
01182   MESSAGE_MED("Constructeur FIELD de recopie");
01183 
01184   // RECOPIE PROFONDE <> de l'operateur= Rmq from EF
01185   if (m._value != NULL)
01186     {
01187       if ( m.getGaussPresence() )
01188         _value = new ArrayGauss( *(static_cast< ArrayGauss * > (m._value) ) ,false);
01189       else
01190         _value = new ArrayNoGauss( *(static_cast< ArrayNoGauss * > (m._value)) ,false);
01191     }
01192   else
01193     _value = (ArrayNoGauss *) NULL;
01194   locMap::const_iterator it;
01195   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
01196     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
01197       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
01198                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
01199                                               );
01200 
01201   _valueType       = m._valueType;
01202   _interlacingType = m._interlacingType;
01203   //drivers = m._drivers;
01204   _mesh            = m._mesh;
01205   if(_mesh)
01206     _mesh->addReference();
01207 }
01208 
01214 template <class T, class INTERLACING_TAG>
01215 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(const FIELD &m)
01216 {
01217   MESSAGE_MED("Appel de FIELD<T>::operator=") ;
01218   if ( this == &m) return *this;
01219 
01220   // copy values array
01221   FIELD_::operator=(m);  // Driver are ignored & ?copie su pointeur de Support?
01222 
01223   _value           = m._value; //PROBLEME RECOPIE DES POINTEURS PAS COHERENT AVEC LE
01224                                //CONSTRUCTEUR PAR RECOPIE
01225                                //CF :Commentaire dans MEDMEM_Array
01226   locMap::const_iterator it;
01227   for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
01228     _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
01229       new GAUSS_LOCALIZATION<INTERLACING_TAG>(
01230                                               *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
01231                                               );
01232 
01233   _valueType       = m._valueType;
01234   _interlacingType = m._interlacingType;
01235   if(_mesh!=m._mesh)
01236     {
01237       if(_mesh)
01238         _mesh->removeReference();
01239       _mesh = m._mesh;
01240       if(_mesh)
01241         _mesh->addReference();
01242     }
01243   return *this;
01244 }
01245 
01249 template <class T, class INTERLACING_TAG>
01250 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(T value)
01251 {
01252   MESSAGE_MED("Appel de FIELD<T>::operator= T") ;
01253         int size=getNumberOfComponents()*getNumberOfValues();
01254         T* ptr= const_cast<T*>( getValue());
01255         for (int i=0; i< size; i++)
01256                 {*ptr++=value;}
01257 
01258   return *this;
01259 }
01260 
01285 template <class T, class INTERLACING_TAG>
01286 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator+(const FIELD & m) const
01287 {
01288   const char* LOC = "FIELD<T>::operator+(const FIELD & m)";
01289   BEGIN_OF_MED(LOC);
01290     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
01291 
01292     // Creation of the result - memory is allocated by FIELD constructor
01293     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01294     result->_operationInitialize(*this,m,"+"); // perform Atribute's initialization
01295     result->_add_in_place(*this,m); // perform addition
01296 
01297   END_OF_MED(LOC);
01298     return result;
01299 }
01300 
01305 template <class T, class INTERLACING_TAG>
01306 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator+=(const FIELD & m)
01307 {
01308   const char* LOC = "FIELD<T>::operator+=(const FIELD & m)";
01309   BEGIN_OF_MED(LOC);
01310     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
01311 
01312     const T* value1=m.getValue(); // get pointers to the values we are adding
01313 
01314     // get a non const pointer to the inside array of values and perform operation
01315     T * value=const_cast<T *> (getValue());
01316     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
01317     const T* endV=value+size; // pointer to the end of value
01318     for(;value!=endV; value1++,value++)
01319         *value += *value1;
01320   END_OF_MED(LOC);
01321     return *this;
01322 }
01323 
01324 
01330 template <class T, class INTERLACING_TAG>
01331 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::add(const FIELD& m, const FIELD& n)
01332 {
01333   const char* LOC = "FIELD<T>::add(const FIELD & m, const FIELD& n)";
01334   BEGIN_OF_MED(LOC);
01335     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
01336 
01337     // Creation of a new field
01338     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01339                                                                       m.getNumberOfComponents());
01340     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
01341     result->_add_in_place(m,n); // perform addition
01342 
01343   END_OF_MED(LOC);
01344     return result;
01345 }
01346 
01349 template <class T, class INTERLACING_TAG>
01350 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::addDeep(const FIELD& m, const FIELD& n)
01351 {
01352   const char* LOC = "FIELD<T>::addDeep(const FIELD & m, const FIELD& n)";
01353   BEGIN_OF_MED(LOC);
01354     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
01355 
01356     // Creation of a new field
01357     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01358                                                                       m.getNumberOfComponents());
01359     result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
01360     result->_add_in_place(m,n); // perform addition
01361 
01362   END_OF_MED(LOC);
01363     return result;
01364 }
01365 
01386 template <class T, class INTERLACING_TAG>
01387 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator-(const FIELD & m) const
01388 {
01389   const char* LOC = "FIELD<T>::operator-(const FIELD & m)";
01390   BEGIN_OF_MED(LOC);
01391     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
01392 
01393     // Creation of the result - memory is allocated by FIELD constructor
01394     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01395     result->_operationInitialize(*this,m,"-"); // perform Atribute's initialization
01396     result->_sub_in_place(*this,m); // perform substracion
01397 
01398   END_OF_MED(LOC);
01399     return result;
01400 }
01401 
01402 template <class T, class INTERLACING_TAG>
01403 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator-() const
01404 {
01405   const char* LOC = "FIELD<T>::operator-()";
01406   BEGIN_OF_MED(LOC);
01407 
01408     // Creation of the result - memory is allocated by FIELD constructor
01409   FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01410     // Atribute's initialization
01411     result->setName("- "+getName());
01412     result->setComponentsNames(getComponentsNames());
01413     // not yet implemented    setComponentType(getComponentType());
01414     result->setComponentsDescriptions(getComponentsDescriptions());
01415     result->setMEDComponentsUnits(getMEDComponentsUnits());
01416     result->setComponentsUnits(getComponentsUnits());
01417     result->setIterationNumber(getIterationNumber());
01418     result->setTime(getTime());
01419     result->setOrderNumber(getOrderNumber());
01420 
01421     const T* value1=getValue();
01422     // get a non const pointer to the inside array of values and perform operation
01423     T * value=const_cast<T *> (result->getValue());
01424     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
01425     const T* endV=value+size; // pointer to the end of value
01426 
01427     for(;value!=endV; value1++,value++)
01428         *value = -(*value1);
01429   END_OF_MED(LOC);
01430     return result;
01431 }
01432 
01437 template <class T, class INTERLACING_TAG>
01438 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator-=(const FIELD & m)
01439 {
01440   const char* LOC = "FIELD<T>::operator-=(const FIELD & m)";
01441   BEGIN_OF_MED(LOC);
01442     FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
01443 
01444     const T* value1=m.getValue();
01445 
01446     // get a non const pointer to the inside array of values and perform operation
01447     T * value=const_cast<T *> (getValue());
01448     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
01449     const T* endV=value+size; // pointer to the end of value
01450 
01451     for(;value!=endV; value1++,value++)
01452         *value -= *value1;
01453 
01454   END_OF_MED(LOC);
01455     return *this;
01456 }
01457 
01458 
01459 
01463 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b, int icomp)
01464 {
01465     // get a non const pointer to the inside array of values and perform operation in place
01466     T * value=const_cast<T *> (getValue());
01467          
01468     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
01469 
01470     if (size>0) // for a negative size, there is nothing to do
01471     {
01472                         value+=icomp-1;
01473                         const T* lastvalue=value+size; // pointer to the end of value
01474  
01475                         for(;value!=lastvalue; value+=getNumberOfComponents()) // apply linear transformation
01476                                 *value = a*(*value)+b;
01477     }
01478 }
01479 
01485 template <class T, class INTERLACING_TAG>
01486 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::sub(const FIELD& m, const FIELD& n)
01487 {
01488   const char* LOC = "FIELD<T>::sub(const FIELD & m, const FIELD& n)";
01489   BEGIN_OF_MED(LOC);
01490     FIELD_::_checkFieldCompatibility(m, n); // may throw exception
01491 
01492     // Creation of a new field
01493     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01494                                                                       m.getNumberOfComponents());
01495     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
01496     result->_sub_in_place(m,n); // perform substraction
01497 
01498   END_OF_MED(LOC);
01499     return result;
01500 }
01501 
01504 template <class T, class INTERLACING_TAG>
01505 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::subDeep(const FIELD& m, const FIELD& n)
01506 {
01507   const char* LOC = "FIELD<T>::subDeep(const FIELD & m, const FIELD& n)";
01508   BEGIN_OF_MED(LOC);
01509     FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
01510 
01511     // Creation of a new field
01512     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01513                                                                       m.getNumberOfComponents());
01514     result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
01515     result->_sub_in_place(m,n); // perform substraction
01516 
01517   END_OF_MED(LOC);
01518     return result;
01519 }
01520 
01541 template <class T, class INTERLACING_TAG>
01542 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator*(const FIELD & m) const
01543 {
01544   const char* LOC = "FIELD<T>::operator*(const FIELD & m)";
01545   BEGIN_OF_MED(LOC);
01546     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
01547 
01548     // Creation of the result - memory is allocated by FIELD constructor
01549     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01550     result->_operationInitialize(*this,m,"*"); // perform Atribute's initialization
01551     result->_mul_in_place(*this,m); // perform multiplication
01552 
01553   END_OF_MED(LOC);
01554     return result;
01555 }
01556 
01561 template <class T, class INTERLACING_TAG>
01562 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator*=(const FIELD & m)
01563 {
01564   const char* LOC = "FIELD<T>::operator*=(const FIELD & m)";
01565   BEGIN_OF_MED(LOC);
01566     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
01567 
01568     const T* value1=m.getValue();
01569 
01570     // get a non const pointer to the inside array of values and perform operation
01571     T * value=const_cast<T *> (getValue());
01572     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
01573     const T* endV=value+size; // pointer to the end of value
01574 
01575     for(;value!=endV; value1++,value++)
01576         *value *= *value1;
01577 
01578   END_OF_MED(LOC);
01579     return *this;
01580 }
01581 
01582 
01588 template <class T, class INTERLACING_TAG>
01589 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mul(const FIELD& m, const FIELD& n)
01590 {
01591   const char* LOC = "FIELD<T>::mul(const FIELD & m, const FIELD& n)";
01592   BEGIN_OF_MED(LOC);
01593     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
01594 
01595     // Creation of a new field
01596     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01597                                                                       m.getNumberOfComponents());
01598     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
01599     result->_mul_in_place(m,n); // perform multiplication
01600 
01601   END_OF_MED(LOC);
01602     return result;
01603 }
01604 
01607 template <class T, class INTERLACING_TAG>
01608 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mulDeep(const FIELD& m, const FIELD& n)
01609 {
01610   const char* LOC = "FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)";
01611   BEGIN_OF_MED(LOC);
01612     FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
01613 
01614     // Creation of a new field
01615     FIELD<T, INTERLACING_TAG>* result = new FIELD<T,INTERLACING_TAG>(m.getSupport(),
01616                                                                      m.getNumberOfComponents());
01617     result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
01618     result->_mul_in_place(m,n); // perform multiplication
01619 
01620   END_OF_MED(LOC);
01621     return result;
01622 }
01623 
01644 template <class T, class INTERLACING_TAG>
01645 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator/(const FIELD & m) const
01646 {
01647   const char* LOC = "FIELD<T>::operator/(const FIELD & m)";
01648   BEGIN_OF_MED(LOC);
01649     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
01650 
01651     // Creation of the result - memory is allocated by FIELD constructor
01652     FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01653     try
01654       {
01655         result->_operationInitialize(*this,m,"/"); // perform Atribute's initialization
01656         result->_div_in_place(*this,m); // perform division
01657       }
01658     catch(MEDEXCEPTION& e)
01659       {
01660         result->removeReference();
01661         throw e;
01662       }
01663 
01664   END_OF_MED(LOC);
01665     return result;
01666 }
01667 
01668 
01673 template <class T, class INTERLACING_TAG>
01674 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator/=(const FIELD & m)
01675 {
01676   const char* LOC = "FIELD<T>::operator/=(const FIELD & m)";
01677   BEGIN_OF_MED(LOC);
01678     FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
01679 
01680     const T* value1=m.getValue(); // get pointers to the values we are adding
01681 
01682     // get a non const pointer to the inside array of values and perform operation
01683     T * value=const_cast<T *> (getValue());
01684     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
01685     const T* endV=value+size; // pointer to the end of value
01686 
01687     for(;value!=endV; value1++,value++)
01688         *value /= *value1;
01689 
01690   END_OF_MED(LOC);
01691     return *this;
01692 }
01693 
01694 
01700 template <class T, class INTERLACING_TAG>
01701 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::div(const FIELD& m, const FIELD& n)
01702 {
01703   const char* LOC = "FIELD<T>::div(const FIELD & m, const FIELD& n)";
01704   BEGIN_OF_MED(LOC);
01705     FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
01706 
01707     // Creation of a new field
01708     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01709                                                                       m.getNumberOfComponents());
01710     try
01711       {
01712         result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
01713         result->_div_in_place(m,n); // perform division
01714       }
01715     catch(MEDEXCEPTION& e)
01716       {
01717         result->removeReference();
01718         throw e;
01719       }
01720   END_OF_MED(LOC);
01721     return result;
01722 }
01723 
01726 template <class T,class INTERLACING_TAG>
01727 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::divDeep(const FIELD& m, const FIELD& n)
01728 {
01729   const char* LOC = "FIELD<T>::divDeep(const FIELD & m, const FIELD& n)";
01730   BEGIN_OF_MED(LOC);
01731   FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
01732 
01733   // Creation of a new field
01734   FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01735                                                                     m.getNumberOfComponents());
01736   try
01737     {
01738       result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
01739       result->_div_in_place(m,n); // perform division
01740     }
01741   catch(MEDEXCEPTION& e)
01742       {
01743         result->removeReference();
01744         throw e;
01745       }
01746   END_OF_MED(LOC);
01747   return result;
01748 }
01749 
01750 
01762 template <class T, class INTERLACING_TAG>
01763 void FIELD<T, INTERLACING_TAG>::_operationInitialize(const FIELD& m,const FIELD& n, const char* Op)
01764 {
01765     MESSAGE_MED("Appel methode interne " << Op);
01766 
01767     // Atribute's initialization (copy of the first field's attributes)
01768     // Other data members (_support, _numberOfValues) are initialized in the field's constr.
01769     setName(m.getName()+" "+Op+" "+n.getName());
01770     setComponentsNames(m.getComponentsNames());
01771     setComponentsDescriptions(m.getComponentsDescriptions());
01772     setMEDComponentsUnits(m.getMEDComponentsUnits());
01773 
01774     // The following data member may differ from field m to n.
01775     // The initialization is done based on the first field.
01776 
01777     setComponentsUnits(m.getComponentsUnits());
01778 
01779     setIterationNumber(m.getIterationNumber());
01780     setTime(m.getTime());
01781     setOrderNumber(m.getOrderNumber());
01782 }
01783 
01784 
01793 template <class T, class INTERLACING_TAG>
01794 void FIELD<T, INTERLACING_TAG>::_add_in_place(const FIELD& m,const FIELD& n)
01795 {
01796     // get pointers to the values we are adding
01797     const T* value1=m.getValue();
01798     const T* value2=n.getValue();
01799     // get a non const pointer to the inside array of values and perform operation
01800     T * value=const_cast<T *> (getValue());
01801 
01802     const int size=getNumberOfValues()*getNumberOfComponents();
01803     SCRUTE_MED(size);
01804     const T* endV1=value1+size;
01805     for(;value1!=endV1; value1++,value2++,value++)
01806         *value=(*value1)+(*value2);
01807 }
01808 
01817 template <class T, class INTERLACING_TAG>
01818 void FIELD<T, INTERLACING_TAG>::_sub_in_place(const FIELD& m,const FIELD& n)
01819 {
01820     // get pointers to the values we are adding
01821     const T* value1=m.getValue();
01822     const T* value2=n.getValue();
01823     // get a non const pointer to the inside array of values and perform operation
01824     T * value=const_cast<T *> (getValue());
01825 
01826     const int size=getNumberOfValues()*getNumberOfComponents();
01827     SCRUTE_MED(size);
01828     const T* endV1=value1+size;
01829     for(;value1!=endV1; value1++,value2++,value++)
01830         *value=(*value1)-(*value2);
01831 }
01832 
01841 template <class T, class INTERLACING_TAG>
01842 void FIELD<T, INTERLACING_TAG>::_mul_in_place(const FIELD& m,const FIELD& n)
01843 {
01844     // get pointers to the values we are adding
01845     const T* value1=m.getValue();
01846     const T* value2=n.getValue();
01847     // get a non const pointer to the inside array of values and perform operation
01848     T * value=const_cast<T *> (getValue());
01849 
01850     const int size=getNumberOfValues()*getNumberOfComponents();
01851     SCRUTE_MED(size);
01852     const T* endV1=value1+size;
01853     for(;value1!=endV1; value1++,value2++,value++)
01854         *value=(*value1)*(*value2);
01855 }
01856 
01865 template <class T, class INTERLACING_TAG>
01866 void FIELD<T, INTERLACING_TAG>::_div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION)
01867 {
01868     // get pointers to the values we are adding
01869     const T* value1=m.getValue();
01870     const T* value2=n.getValue();
01871     // get a non const pointer to the inside array of values and perform operation
01872     T * value=const_cast<T *> (getValue());
01873 
01874     const int size=getNumberOfValues()*getNumberOfComponents();
01875     SCRUTE_MED(size);
01876     const T* endV1=value1+size;
01877     for(;value1!=endV1; value1++,value2++,value++){
01878       if ( *value2 == 0 ) { // FAIRE PLUTOT UN TRY CATCH Rmq from EF
01879           string diagnosis;
01880           diagnosis="FIELD<T,INTERLACING_TAG>::_div_in_place(...) : Divide by zero !";
01881           throw MEDEXCEPTION(diagnosis.c_str());
01882         }
01883         *value=(*value1)/(*value2);
01884     }
01885 }
01886 
01894 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::normMax() const throw (MEDEXCEPTION)
01895 {
01896     const T* value=getValue(); // get pointer to the values
01897     const int size=getNumberOfValues()*getNumberOfComponents();
01898     if (size <= 0) // Size of array has to be strictly positive
01899     {
01900         string diagnosis;
01901         diagnosis="FIELD<T,INTERLACIN_TAG>::normMax() : cannot compute the norm of "+getName()+
01902             " : it size is non positive!";
01903         throw MEDEXCEPTION(diagnosis.c_str());
01904     }
01905     const T* lastvalue=value+size; // get pointer just after last value
01906     const T* pMax=value; // pointer to the max value
01907     const T* pMin=value; // pointer to the min value
01908 
01909     // get pointers to the max & min value of array
01910     while ( ++value != lastvalue )
01911     {
01912         if ( *pMin > *value )
01913             pMin=value;
01914         if ( *pMax < *value )
01915             pMax=value;
01916     }
01917 
01918     T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
01919     T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
01920 
01921     return Max>Min ? double(Max) : double(Min);
01922 }
01923 
01926 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::norm2() const throw (MEDEXCEPTION)
01927 {
01928     const T* value=this->getValue(); // get const pointer to the values
01929     const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
01930     if (size <= 0) // Size of array has to be strictly positive
01931     {
01932         string diagnosis;
01933         diagnosis="FIELD<T,INTERLACIN_TAG>::norm2() : cannot compute the norm of "+getName()+
01934             " : it size is non positive!";
01935         throw MEDEXCEPTION(diagnosis.c_str());
01936     }
01937     const T* lastvalue=value+size; // point just after last value
01938 
01939     T result((T)0); // init
01940     for( ; value!=lastvalue ; ++value)
01941         result += (*value) * (*value);
01942 
01943     return std::sqrt(double(result));
01944 }
01945 
01946 
01947 //------------- TDG and BS addings 
01948 
01951  template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION)
01952 {
01953   const T* value=getValue(); // get pointer to the values
01954   const int size=getNumberOfValues()*getNumberOfComponents();
01955   const T* lastvalue=value+size; // point just after last value
01956     
01957   if (size <= 0){ // Size of array has to be strictly positive
01958       
01959     string diagnosis;
01960     diagnosis="FIELD<T,INTERLACIN_TAG>::getMinMax() : cannot compute the extremums of "+getName()+
01961       " : its size is non positive!";
01962     throw MEDEXCEPTION(diagnosis.c_str());
01963   }
01964     
01965   if (!_isMinMax){
01966     vmax=MinMax<T>::getMin(); // init a max value
01967     vmin=MinMax<T>::getMax(); // init a min value
01968       
01969     for( ; value!=lastvalue ; ++value){
01970       if ( vmin > *value )
01971         vmin=*value;
01972       if ( vmax < *value )
01973         vmax=*value;
01974     }
01975     _isMinMax=true;
01976     _vmin=vmin;
01977     _vmax=vmax;
01978   }
01979   else{
01980     vmin = _vmin;
01981     vmax = _vmax;
01982   }
01983 
01984 }
01985 
01988  template <class T, class INTERLACIN_TAG> vector<int> FIELD<T, INTERLACIN_TAG>::getHistogram(int &nbint) throw (MEDEXCEPTION)
01989 {
01990   const T* value=getValue(); // get pointer to the values
01991   const int size=getNumberOfValues()*getNumberOfComponents();
01992   const T* lastvalue=value+size; // point just after last value
01993 
01994   if (size <= 0){ // Size of array has to be strictly positive
01995 
01996     string diagnosis;
01997     diagnosis="FIELD<T,INTERLACIN_TAG>::getHistogram() : cannot compute the histogram of "+getName()+
01998       " : it size is non positive!";
01999     throw MEDEXCEPTION(diagnosis.c_str());
02000   }
02001 
02002   vector<int> Histogram(nbint) ;
02003   T vmin,vmax;
02004   int j;
02005 
02006   for( j=0 ; j!=nbint ; j++) Histogram[j]=0 ;
02007     
02008   getMinMax(vmin,vmax);
02009   for( ; value!=lastvalue ; ++value){
02010     if(*value==vmax) j = nbint-1;
02011     else j = (int)(((double)nbint * (*value-vmin))/(vmax-vmin));
02012     Histogram[j]+=1 ;
02013   }
02014 
02015   return Histogram ;
02016 
02017 }
02018 
02021 template <class T, class INTERLACIN_TAG> 
02022 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildGradient() const throw (MEDEXCEPTION)
02023 {
02024   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildGradient() : ";
02025   BEGIN_OF_MED(LOC);
02026 
02027   // space dimension of input mesh
02028   int spaceDim = getSupport()->getMesh()->getSpaceDimension();
02029   double *x = new double[spaceDim];
02030 
02031   FIELD<double, FullInterlace>* Gradient =
02032     new FIELD<double, FullInterlace>(getSupport(),spaceDim);
02033 
02034   string name("gradient of ");
02035   name += getName();
02036   Gradient->setName(name);
02037   string descr("gradient of ");
02038   descr += getDescription();
02039   Gradient->setDescription(descr);
02040 
02041   if( _numberOfComponents > 1 ){
02042     delete Gradient;
02043     delete [] x;
02044     throw MEDEXCEPTION("gradient calculation only on scalar field");
02045   }
02046 
02047   for(int i=1;i<=spaceDim;i++){
02048     string nameC("gradient of ");
02049     nameC += getName();
02050     Gradient->setComponentName(i,nameC);
02051     Gradient->setComponentDescription(i,"gradient");
02052     string MEDComponentUnit = getMEDComponentUnit(1)+getSupport()->getMesh()->getCoordinatesUnits()[i-1];
02053     Gradient->setMEDComponentUnit(i,MEDComponentUnit);
02054   }
02055 
02056   Gradient->setIterationNumber(getIterationNumber());
02057   Gradient->setOrderNumber(getOrderNumber());
02058   Gradient->setTime(getTime());
02059 
02060   // typ of entity on what is field
02061   MED_EN::medEntityMesh typ = getSupport()->getEntity();
02062 
02063   const int *C;
02064   const int *iC;  
02065   const int *revC;
02066   const int *indC;
02067   const double *coord;
02068   int NumberOf;
02069 
02070   switch (typ) {
02071   case MED_CELL:
02072   case MED_FACE:
02073   case MED_EDGE:
02074     {
02075       // read connectivity array to have the list of nodes contained by an element
02076       C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,typ,MED_ALL_ELEMENTS);
02077       iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,typ);
02078       // calculate reverse connectivity to have the list of elements which contains node i
02079       revC = getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
02080       indC = getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
02081       // coordinates of each node
02082       coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
02083       // number of elements
02084       NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
02085       // barycenter field of elements
02086       FIELD<double, FullInterlace>* barycenter = getSupport()->getMesh()->getBarycenter(getSupport());
02087 
02088       // calculate gradient vector for each element i
02089       for (int i = 1; i < NumberOf + 1; i++) {
02090 
02091         // listElements contains elements which contains a node of element i
02092         set <int> listElements;
02093         set <int>::iterator elemIt;
02094         listElements.clear();
02095 
02096         // for each node j of element i
02097         for (int ij = iC[i-1]; ij < iC[i]; ij++) {
02098           int j = C[ij-1];
02099           for (int k = indC[j-1]; k < indC[j]; k++) {
02100             // c element contains node j
02101             int c = revC[k-1];
02102             // we put the elements in set
02103             if (c != i)
02104               listElements.insert(c);
02105           }
02106         }
02107         // coordinates of barycentre of element i in space of dimension spaceDim
02108         for (int j = 0; j < spaceDim; j++)
02109           x[j] = barycenter->getValueIJ(i,j+1);
02110 
02111         for (int j = 0; j < spaceDim; j++) {
02112           // value of field of element i
02113           double val = getValueIJ(i,1);
02114           double grad = 0.;
02115           // calculate gradient for each neighbor element
02116           for (elemIt = listElements.begin(); elemIt != listElements.end(); elemIt++) {
02117             int elem = *elemIt;
02118             double d2 = 0.;
02119             for (int l = 0; l < spaceDim; l++) {
02120               // coordinate of barycenter of element elem
02121               double xx = barycenter->getValueIJ(elem, l+1);
02122               d2 += (x[l]-xx) * (x[l]-xx);
02123             }
02124             grad += (barycenter->getValueIJ(elem,j+1)-x[j])*(getValueIJ(elem,1)-val)/sqrt(d2);
02125           }
02126           if (listElements.size() != 0) grad /= listElements.size();
02127           Gradient->setValueIJ(i,j+1,grad);
02128         }
02129       }
02130       delete barycenter;
02131     }
02132     break;
02133   case MED_NODE:
02134     // read connectivity array to have the list of nodes contained by an element
02135     C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
02136     iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,MED_CELL);
02137     // calculate reverse connectivity to have the list of elements which contains node i
02138     revC=getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
02139     indC=getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
02140     // coordinates of each node
02141     coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
02142 
02143     // calculate gradient for each node
02144     NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
02145     for (int i=1; i<NumberOf+1; i++){
02146       // listNodes contains nodes neigbor of node i 
02147       set <int> listNodes;
02148       set <int>::iterator nodeIt ;
02149       listNodes.clear();
02150       for(int j=indC[i-1];j<indC[i];j++){
02151         // c element contains node i
02152         int c=revC[j-1];
02153         // we put the nodes of c element in set
02154         for(int k=iC[c-1];k<iC[c];k++)
02155           if(C[k-1] != i)
02156             listNodes.insert(C[k-1]);
02157       }
02158       // coordinates of node i in space of dimension spaceDim
02159       for(int j=0;j<spaceDim;j++)
02160         x[j] = coord[(i-1)*spaceDim+j];
02161       
02162       for(int j=0;j<spaceDim;j++){
02163         // value of field
02164         double val = getValueIJ(i,1);
02165         double grad = 0.;
02166         // calculate gradient for each neighbor node
02167         for(nodeIt=listNodes.begin();nodeIt!=listNodes.end();nodeIt++){
02168           int node = *nodeIt;
02169           double d2 = 0.;
02170           for(int l=0;l<spaceDim;l++){
02171             double xx = coord[(node-1)*spaceDim+l];
02172             d2 += (x[l]-xx) * (x[l]-xx);
02173           }
02174           grad += (coord[(node-1)*spaceDim+j]-x[j])*(getValueIJ(node,1)-val)/sqrt(d2);
02175         }
02176         if(listNodes.size() != 0) grad /= listNodes.size();
02177         Gradient->setValueIJ(i,j+1,grad);
02178       }
02179     }
02180     break;
02181   case MED_ALL_ENTITIES:
02182     delete [] x;
02183     delete Gradient;
02184     throw MEDEXCEPTION("gradient calculation not yet implemented on all elements");
02185     break;
02186   }
02187 
02188   delete [] x;
02189 
02190   END_OF_MED(LOC);
02191   return Gradient;
02192 }
02193 
02196 template <class T, class INTERLACIN_TAG> 
02197 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildNorm2Field() const throw (MEDEXCEPTION)
02198 {
02199   const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildNorm2Field() : ";
02200   BEGIN_OF_MED(LOC);
02201 
02202   FIELD<double, FullInterlace>* Norm2Field =
02203     new FIELD<double, FullInterlace>(getSupport(),1);
02204 
02205   string name("norm2 of ");
02206   name += getName();
02207   Norm2Field->setName(name);
02208   string descr("norm2 of ");
02209   descr += getDescription();
02210   Norm2Field->setDescription(descr);
02211 
02212   string nameC("norm2 of ");
02213   nameC += getName();
02214   Norm2Field->setComponentName(1,nameC);
02215   Norm2Field->setComponentDescription(1,"norm2");
02216   string MEDComponentUnit = getMEDComponentUnit(1);
02217   Norm2Field->setMEDComponentUnit(1,MEDComponentUnit);
02218 
02219   Norm2Field->setIterationNumber(getIterationNumber());
02220   Norm2Field->setOrderNumber(getOrderNumber());
02221   Norm2Field->setTime(getTime());
02222 
02223   // calculate nom2 for each element
02224   int NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
02225   for (int i=1; i<NumberOf+1; i++){
02226     double norm2 = 0.;
02227     for(int j=1;j<=getNumberOfComponents();j++)
02228       norm2 += getValueIJ(i,j)*getValueIJ(i,j);
02229     Norm2Field->setValueIJ(i,1,sqrt(norm2));
02230   }
02231 
02232   END_OF_MED(LOC);
02233   return Norm2Field;
02234 
02235 }
02236 
02246 template <class T, class INTERLACIN_TAG> template <T T_function(T)>
02247 void FIELD<T, INTERLACIN_TAG>::applyFunc()
02248 {
02249   // get a non const pointer to the inside array of values and perform operation
02250     T * value=const_cast<T *> (getValue());
02251     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
02252 
02253     if (size>0) // for a negative size, there is nothing to do
02254     {
02255         const T* lastvalue=value+size; // pointer to the end of value
02256         for(;value!=lastvalue; ++value) // apply linear transformation
02257             *value = T_function(*value);
02258     }
02259 }
02260 
02261 template <class T, class INTERLACIN_TAG> T FIELD<T, INTERLACIN_TAG>::pow(T x)
02262 {
02263   return (T)::pow((double)x,FIELD<T, INTERLACIN_TAG>::_scalarForPow);
02264 }
02265 
02273 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyPow(T scalar)
02274 {
02275   FIELD<T, INTERLACIN_TAG>::_scalarForPow=scalar;
02276   applyFunc<FIELD<T, INTERLACIN_TAG>::pow>();
02277 }
02278 
02282 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b)
02283 {
02284     // get a non const pointer to the inside array of values and perform operation in place
02285     T * value=const_cast<T *> (getValue());
02286     const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
02287 
02288     if (size>0) // for a negative size, there is nothing to do
02289     {
02290         const T* lastvalue=value+size; // pointer to the end of value
02291         for(;value!=lastvalue; ++value) // apply linear transformation
02292             *value = a*(*value)+b;
02293     }
02294 }
02295 
02296 
02311 template <class T, class INTERLACING_TAG>
02312 FIELD<T, INTERLACING_TAG>*
02313 FIELD<T, INTERLACING_TAG>::scalarProduct(const FIELD & m, const FIELD & n, bool deepCheck)
02314 {
02315     if(!deepCheck)
02316       FIELD_::_checkFieldCompatibility( m, n, false); // may throw exception
02317     else
02318       FIELD_::_deepCheckFieldCompatibility(m, n, false);
02319 
02320     // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
02321     // result type imply INTERLACING_TAG=FullInterlace for m & n
02322 
02323     const int numberOfElements=m.getNumberOfValues(); // strictly positive
02324     const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
02325 
02326     // Creation & init of a the result field on the same support, with one component
02327     // You have to be careful about the interlacing mode, because in the computation step,
02328     // it seems to assume the the interlacing mode is the FullInterlacing
02329 
02330     FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),1);
02331     result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
02332     result->setIterationNumber(m.getIterationNumber());
02333     result->setTime(m.getTime());
02334     result->setOrderNumber(m.getOrderNumber());
02335 
02336     const T* value1=m.getValue(); // get const pointer to the values
02337     const T* value2=n.getValue(); // get const pointer to the values
02338     // get a non const pointer to the inside array of values and perform operation
02339     T * value=const_cast<T *> (result->getValue());
02340 
02341     const T* lastvalue=value+numberOfElements; // pointing just after last value of result
02342     for ( ; value!=lastvalue ; ++value ) // loop on all elements
02343     {
02344         *value=(T)0; // initialize value
02345         const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
02346         for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
02347             *value += (*value1) * (*value2);
02348     }
02349     return result;
02350 }
02351 
02357 template <class T, class INTERLACING_TAG>
02358 double FIELD<T, INTERLACING_TAG>::normL2(int component,
02359                                          const FIELD<double, FullInterlace> * p_field_volume) const
02360 {
02361     _checkNormCompatibility(p_field_volume, /*nodalAllowed=*/true); // may throw exception
02362     if ( component<1 || component>getNumberOfComponents() )
02363         throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
02364 
02365     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
02366     if(!p_field_volume) // if the user don't supply the volume
02367       p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
02368     else
02369       p_field_size->addReference();
02370     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
02371     const double* vol=p_field_size->getValue();
02372     // Il n'est vraiment pas optimal de mixer des champs dans des modes d'entrelacement
02373     // different juste pour le calcul
02374 
02375 
02376     double integrale=0.0;
02377     double totVol=0.0;
02378 
02379     if ( getSupport()->getEntity() == MED_NODE ) // issue 20120: [CEA 206] normL2 on NODE field
02380     {
02381       //Most frequently the FIELD is on the whole mesh and
02382       // there is no need in optimizing iterations from supporting nodes-> back to cells,
02383       // so we iterate just on all cells
02384       const MESH * mesh = getSupport()->getMesh()->convertInMESH();
02385       const int nbCells = mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
02386       const int *C = mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
02387       const int *iC = mesh->getConnectivityIndex(MED_NODAL,MED_CELL);
02388       for (int i = 0; i < nbCells; ++i, ++vol) {
02389         // calculate integral on current element as average summ of values on all it's nodes
02390         double curCellValue = 0;
02391         try { // we expect exception with partial fields for nodes w/o values
02392           for (int ij = iC[i]; ij < iC[i+1]; ij++) {
02393             int node = C[ij-1];
02394             curCellValue += getValueIJ( node, component );
02395           }
02396         }
02397         catch ( MEDEXCEPTION ) {
02398           continue;
02399         }
02400         int nbNodes = iC[i+1]-iC[i];
02401         curCellValue /= nbNodes;
02402         integrale += (curCellValue * curCellValue) * std::abs(*vol);
02403         totVol+=std::abs(*vol);
02404       }
02405       mesh->removeReference();
02406     }
02407     else
02408     {
02409       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
02410         const T* value = getValue();
02411         value = value + (component-1) * getNumberOfValues();
02412         const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
02413         for (; value!=lastvalue ; ++value ,++vol) {
02414           integrale += double((*value) * (*value)) * std::abs(*vol);
02415           totVol+=std::abs(*vol);
02416         }
02417       }
02418       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
02419         ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
02420         for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
02421           integrale += anArray->getIJ(i,component) * anArray->getIJ(i,component) * std::abs(*vol);
02422           totVol+=std::abs(*vol);
02423         }
02424       }
02425       else { // FULL_INTERLACE
02426         ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
02427         for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
02428           integrale += anArray->getIJ(i,component) * anArray->getIJ(i,component) * std::abs(*vol);
02429           totVol+=std::abs(*vol);
02430         }
02431       }
02432     }
02433 
02434     if(p_field_size)
02435       p_field_size->removeReference(); // delete temporary volume field
02436 
02437     if( totVol <= 0)
02438         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
02439     return integrale/totVol;
02440 }
02441 
02446 template <class T, class INTERLACING_TAG>
02447 double FIELD<T, INTERLACING_TAG>::normL2(const FIELD<double, FullInterlace> * p_field_volume) const
02448 {
02449     _checkNormCompatibility(p_field_volume, /*nodalAllowed=*/true); // may throw exception
02450     const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
02451     if(!p_field_volume) // if the user don't supply the volume
02452       p_field_size=_getFieldSize(); // we calculate the volume
02453     else
02454       p_field_size->addReference();
02455     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
02456     const double* vol=p_field_size->getValue();
02457     const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
02458 
02459 
02460     double integrale=0.0;
02461     double totVol=0.0;
02462 
02463     if ( getSupport()->getEntity() == MED_NODE ) // issue 20120: [CEA 206] normL2 on NODE field
02464     {
02465       //Most frequently the FIELD is on the whole mesh and
02466       // there is no need in optimizing iterations from supporting nodes-> back to cells,
02467       // so we iterate just on all cells
02468       const MESH * mesh = getSupport()->getMesh()->convertInMESH();
02469       const int nbCells = mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
02470       const int *C = mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
02471       const int *iC = mesh->getConnectivityIndex(MED_NODAL,MED_CELL);
02472       int nbComp = getNumberOfComponents();
02473       for (int i = 0; i < nbCells; ++i, ++vol) {
02474         // calculate integral on current element as average summ of values on all it's nodes
02475         int nbNodes = iC[i+1]-iC[i];
02476         vector< double > curCellValue( nbComp, 0 );
02477         try { // we expect exception with partial fields for nodes w/o values
02478           for (int ij = iC[i]; ij < iC[i+1]; ij++) {
02479             int node = C[ij-1];
02480             for ( int j = 0; j < nbComp; ++j )
02481               curCellValue[ j ] += getValueIJ( node, j+1 ) / nbNodes;
02482           }
02483         }
02484         catch ( MEDEXCEPTION ) {
02485           continue;
02486         }
02487 
02488         for ( int j = 0; j < nbComp; ++j ) {
02489           integrale += (curCellValue[j] * curCellValue[j]) * std::abs(*vol);
02490         }
02491         totVol+=std::abs(*vol);
02492       }
02493       mesh->removeReference();
02494       if ( nbCells > 0 && totVol == 0.)
02495         throw MEDEXCEPTION("can't compute sobolev norm : "
02496                            "none of elements has values on all it's nodes");
02497     }
02498     else
02499     {
02500       const double* p_vol=vol;
02501       for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
02502         totVol+=std::abs(*p_vol);
02503 
02504       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
02505         const T* value = getValue();
02506         for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
02507           for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
02508             integrale += (*value) * (*value) * std::abs(*p_vol);
02509           }
02510         }
02511       }
02512       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
02513         ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
02514         for (int j=1; j<=anArray->getDim(); j++) {
02515           int i = 1;
02516           for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
02517             integrale += anArray->getIJ(i,j) * anArray->getIJ(i,j) * std::abs(*p_vol);
02518           }
02519         }
02520       }
02521       else { // FULL_INTERLACE
02522         ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
02523         for (int j=1; j<=anArray->getDim(); j++) {
02524           int i = 1;
02525           for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
02526             integrale += anArray->getIJ(i,j) * anArray->getIJ(i,j) * std::abs(*p_vol);
02527           }
02528         }
02529       }
02530     }
02531     if(p_field_size)
02532         p_field_size->removeReference(); // delete temporary volume field
02533 
02534     if( totVol <= 0)
02535       throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
02536     return integrale/totVol;
02537 }
02538 
02542 template <class T, class INTERLACING_TAG>
02543 double FIELD<T, INTERLACING_TAG>::normL1(int component,
02544                                          const FIELD<double, FullInterlace > * p_field_volume) const
02545 {
02546     _checkNormCompatibility(p_field_volume); // may throw exception
02547     if ( component<1 || component>getNumberOfComponents() )
02548         throw MEDEXCEPTION(STRING("FIELD<T,INTERLACING_TAG>::normL1() : The component argument should be between 1 and the number of components"));
02549 
02550     const FIELD<double,FullInterlace> * p_field_size=p_field_volume;
02551     if(!p_field_volume) // if the user don't supply the volume
02552       p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'impl√√҆mentation dans mesh]
02553     else
02554       p_field_size->addReference();
02555     // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
02556     const double* vol = p_field_size->getValue();
02557 
02558     double integrale=0.0;
02559     double totVol=0.0;
02560 
02561     if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
02562       const T* value = getValue();
02563       const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
02564       for (; value!=lastvalue ; ++value ,++vol) {
02565         integrale += std::abs( *value * *vol );
02566         totVol+=std::abs(*vol);
02567       }
02568     }
02569     else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
02570       ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
02571       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
02572         integrale += std::abs( anArray->getIJ(i,component) * (*vol));
02573         totVol+=std::abs(*vol);
02574       }
02575     }
02576     else { // FULL_INTERLACE
02577       ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
02578       for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
02579         integrale += std::abs( anArray->getIJ(i,component) * *vol);
02580         totVol+=std::abs(*vol);
02581       }
02582     }
02583     
02584     if(p_field_size)
02585       p_field_size->removeReference(); // delete temporary volume field
02586     if( totVol <= 0)
02587         throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
02588     return integrale/totVol;
02589 }
02590 
02595 template <class T, class INTERLACING_TAG>
02596 double FIELD<T, INTERLACING_TAG>::normL1(const FIELD<double, FullInterlace> * p_field_volume) const
02597 {
02598   _checkNormCompatibility(p_field_volume); // may throw exception
02599   const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
02600   if(!p_field_volume) // if the user don't supply the volume
02601     p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'impl√¬©mentation dans mesh]
02602   else
02603     p_field_size->addReference();
02604   // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
02605   const double* vol = p_field_size->getValue();
02606   const double* lastvol = vol+getNumberOfValues(); // pointing just after the end of vol
02607 
02608   double integrale=0.0;
02609   double totVol=0.0;
02610   const double* p_vol=vol;
02611   for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
02612     totVol+=std::abs(*p_vol);
02613 
02614   if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
02615     const T* value = getValue();
02616     for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
02617       for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
02618         integrale += std::abs( *value * *p_vol );
02619       }
02620     }
02621   }
02622   else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
02623     ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
02624     for (int j=1; j<=anArray->getDim(); j++) {
02625       int i = 1;
02626       for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
02627         integrale += std::abs( anArray->getIJ(i,j) * *p_vol );
02628       }
02629     }
02630   }
02631   else { // FULL_INTERLACE
02632     ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
02633     for (int j=1; j<=anArray->getDim(); j++) {
02634       int i = 1;
02635       for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
02636         integrale += std::abs( anArray->getIJ(i,j) * *p_vol );
02637       }
02638     }
02639   }
02640   if(p_field_size)
02641     p_field_size->removeReference();
02642   if( totVol <= 0)
02643     throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
02644   return integrale/totVol;
02645 }
02646 
02652 template <class T, class INTERLACING_TAG>
02653 double FIELD<T, INTERLACING_TAG>::integral(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
02654 {
02655   const char* LOC = "FIELD<>::integral(subSupport): ";
02656 
02657   double integrale = 0;
02658 
02659   if (!subSupport ) subSupport = _support;
02660 
02661   // check feasibility
02662   if ( getGaussPresence() )
02663     throw MEDEXCEPTION(STRING(LOC)<<"Gauss numbers greater than one are not yet implemented!");
02664   if ( subSupport->getEntity() != _support->getEntity())
02665     throw MEDEXCEPTION(STRING(LOC)<<"Different support entity of this field and subSupport");
02666   if ( subSupport->getEntity() == MED_EN::MED_NODE )
02667     throw MEDEXCEPTION(STRING(LOC)<<"Integral of nodal field not yet supported");
02668 
02669   // analyze support
02670   const int nbElems = subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
02671   const bool subOnAll = ( subSupport->isOnAllElements() );
02672   const bool  myOnAll = ( _support->isOnAllElements() );
02673   const int* subNums = !subOnAll ? subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS) : 0;
02674   const int*   myNums = !myOnAll ? _support->getNumber(MED_EN::MED_ALL_ELEMENTS) : 0;
02675   if ( !subOnAll && !subNums )
02676     throw MEDEXCEPTION(STRING(LOC)<<"Invalid support: no element numbers");
02677   if ( !myOnAll && !myNums )
02678     throw MEDEXCEPTION(STRING(LOC)<<"Invalid field support: no element numbers");
02679   if ( subOnAll && !myOnAll )
02680     return integral(NULL);
02681 
02682   // get size of elements
02683   const FIELD<double, FullInterlace> * cellSize=_getFieldSize(subSupport);
02684   const double* size = cellSize->getValue();
02685   const double* lastSize = size + nbElems; // pointing just after the end of size
02686 
02687   const T* value = getValue();
02688 
02689   // calculate integrale
02690   if ( (subOnAll && _support->isOnAllElements()) || subSupport == _support )
02691     {
02692       const double* p_vol;
02693       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
02694         {
02695           for (int j=1; j<=getNumberOfComponents(); ++j)
02696             for ( p_vol=size; p_vol != lastSize; ++value ,++p_vol)
02697               integrale += std::abs( *value * *p_vol );
02698         }
02699       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
02700         {
02701           typename ArrayNoByType::InterlacingPolicy* indexer =
02702             dynamic_cast< typename ArrayNoByType::InterlacingPolicy * > ( getArrayNoGauss() );
02703           for (int i, j=1; j<=getNumberOfComponents(); j++)
02704             for (i = 1, p_vol=size; p_vol!=lastSize; i++, ++p_vol )
02705               integrale += std::abs( value[indexer->getIndex(i,j)] * *p_vol );
02706         }
02707       else  // FULL_INTERLACE
02708         {
02709           for ( p_vol=size; p_vol != lastSize; ++p_vol)
02710             for (int j=0; j<getNumberOfComponents(); ++j, ++value)
02711               integrale += std::abs( *value * *p_vol );
02712         }
02713     }
02714   else
02715     {
02716       // find index for each element of subSupport
02717       PointerOf<int> index;
02718       if ( _support->isOnAllElements() )
02719         {
02720           index.set( subNums );
02721         }
02722       else // find common of two partial supports
02723         {
02724           // hope that numbers are in increasing order
02725           index.set( nbElems );
02726           for (int ii = 0; ii < nbElems; ii++)
02727             index[ii] = 0;
02728           bool allNumsFound = true;
02729           int i = 0, iSub = 0;
02730           for ( ; iSub < nbElems; ++iSub )
02731             {
02732               while ( i < getNumberOfValues() && subNums[iSub] > myNums[i] )
02733                 ++i;
02734               if (i == getNumberOfValues() /*subNums[iSub] > myNums[i]*/) // no more myNums
02735                 {
02736                   index[iSub] = 0; // no such number in myNums
02737                   break;
02738                 }
02739               else if ( subNums[iSub] == myNums[i] ) // elem number found
02740                 index[iSub] = ++i; // -- index counts from 1
02741               else // subNums[iSub] < myNums[i]
02742                 allNumsFound = (index[iSub] = 0); // no such number in myNums
02743             }
02744           if ( iSub != nbElems || !allNumsFound )
02745             {
02746               // check if numbers are in increasing order
02747               bool increasingOrder = true;
02748               for ( iSub = 1; iSub < nbElems && increasingOrder; ++iSub )
02749                 increasingOrder = ( subNums[iSub-1] < subNums[iSub] );
02750               for ( i = 1; i < getNumberOfValues() && increasingOrder; ++i )
02751                 increasingOrder = ( myNums[i-1] < myNums[i] );
02752 
02753               if ( !increasingOrder )
02754                 for ( iSub = 0; iSub < nbElems; ++iSub )
02755                   try
02756                     {
02757                       index[iSub] = _support->getValIndFromGlobalNumber( subNums[iSub] );
02758                     }
02759                   catch (MEDEXCEPTION)
02760                     {
02761                       index[iSub] = 0;
02762                     }
02763             }
02764         }
02765 
02766       // calculation
02767       if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
02768         {
02769           for (int j=0; j<getNumberOfComponents(); ++j)
02770             {
02771               value = getValue() + j * getNumberOfValues();
02772               for ( int i = 0; i < nbElems; ++i )
02773                 if ( index[i] )
02774                   integrale += std::abs( value[ index[i]-1 ] * size[i] );
02775             }
02776         }
02777       else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
02778         {
02779           typename ArrayNoByType::InterlacingPolicy* indexer =
02780             dynamic_cast< typename ArrayNoByType::InterlacingPolicy * > ( getArrayNoGauss() );
02781           for (int j=1; j<=getNumberOfComponents(); j++)
02782             for ( int i = 0; i < nbElems; ++i )
02783               if ( index[i] )
02784                 integrale += std::abs( value[indexer->getIndex(index[i],j)] * size[i] );
02785         }
02786       else  // FULL_INTERLACE
02787         {
02788           const int dim = getNumberOfComponents();
02789           for ( int i = 0; i < nbElems; ++i )
02790             if ( index[i] )
02791               for (int j=0; j<dim; ++j)
02792                 integrale += std::abs( value[ dim*(index[i]-1) + j] * size[i] );
02793         }
02794     }
02795   cellSize->removeReference();
02796   return integrale;
02797 }
02798 
02802 template <class T, class INTERLACING_TAG>
02803 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
02804 {
02805   if(!subSupport->belongsTo(*_support))
02806     throw MEDEXCEPTION("FIELD<T>::extract : subSupport not included in this->_support !");
02807   if(_support->isOnAllElements() && subSupport->isOnAllElements())
02808     return new FIELD<T, INTERLACING_TAG>(*this);
02809 
02810   FIELD<T, INTERLACING_TAG> *ret = new FIELD<T, INTERLACING_TAG>(subSupport,
02811                                                                  _numberOfComponents);
02812 
02813   if(!ret->_value)
02814     throw MEDEXCEPTION("FIELD<T>::extract : invalid support detected !");
02815 
02816   T* valuesToSet=(T*)ret->getValue();
02817 
02818   int nbOfEltsSub=subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
02819   const int *eltsSub=subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS);
02820   T* tempVals=new T[_numberOfComponents];
02821   for(int i=0;i<nbOfEltsSub;i++)
02822     {
02823       if(!getValueOnElement(eltsSub[i],tempVals))
02824         throw MEDEXCEPTION("Problem in belongsTo function !!!");
02825       for(int j=0;j<_numberOfComponents;j++)
02826         valuesToSet[i*_numberOfComponents+j]=tempVals[j];
02827     }
02828   delete [] tempVals;
02829 
02830   ret->copyGlobalInfo(*this);
02831   return ret;
02832 }
02849 template <class T, class INTERLACING_TAG>
02850 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
02851                                  driverTypes driverType,
02852                                  const string & fileName/*=""*/,
02853                                  const string & fieldDriverName/*=""*/,
02854                                  const int iterationNumber,
02855                                  const int orderNumber) throw (MEDEXCEPTION)
02856 {
02857   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) : ";
02858   BEGIN_OF_MED(LOC);
02859 
02860   int current;
02861 
02862   init();
02863 
02864   _mesh  = ( MESH* ) NULL;
02865 
02866   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
02867   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
02868   FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
02869 
02870   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
02871   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
02872   FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
02873 
02874   _support = Support;
02875   //A.G. Addings for RC
02876   if(_support)
02877     _support->addReference();
02878   // OCC 10/03/2006 -- According to the rules defined with help of 
02879   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
02880   // MEDMEM_Array<> template using INTERLACING_TAG parameter of 
02881   // FIELD template - MSVC++ 2003 compiler generated an error here.
02882   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
02883   _value = NULL;
02884 
02885   _iterationNumber = iterationNumber;
02886   _time = 0.0;
02887   _orderNumber = orderNumber;
02888 
02889   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
02890 
02891   _drivers[current]->open();
02892   _drivers[current]->read();
02893   _drivers[current]->close();
02894 
02895   END_OF_MED(LOC);
02896 }
02897 
02909 template <class T, class INTERLACING_TAG>
02910 FIELD<T,INTERLACING_TAG>::FIELD(driverTypes    driverType,
02911                                 const string & fileName,
02912                                 const string & fieldDriverName,
02913                                 const int      iterationNumber,
02914                                 const int      orderNumber,
02915                                 GMESH*         mesh)
02916   throw (MEDEXCEPTION) :FIELD_()
02917 {
02918   int current;
02919   const char* LOC = "FIELD<T,INTERLACING_TAG>::FIELD( driverTypes driverType, const string & fileName, string & fieldDriverName, int iterationNumber, int orderNumber) : ";
02920   BEGIN_OF_MED(LOC);
02921 
02922   init();
02923 
02924   _mesh = mesh;
02925   if(_mesh)
02926     _mesh->addReference();
02927 
02928   //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
02929   ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
02930   FIELD_::_valueType = SET_VALUE_TYPE<T>::_valueType;
02931 
02932   //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
02933   ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
02934   FIELD_::_interlacingType = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
02935 
02936   _support = (SUPPORT *) NULL;
02937   // OCC 10/03/2006 -- According to the rules defined with help of 
02938   // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
02939   // MEDMEM_Array<> template using INTERLACING_TAG parameter of 
02940   // FIELD template - MSVC++ 2003 compiler generated an error here.
02941   // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
02942   _value = NULL;
02943 
02944   _iterationNumber = iterationNumber;
02945   _time = 0.0;
02946   _orderNumber = orderNumber;
02947 
02948   current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
02949 
02950   _drivers[current]->open();
02951   _drivers[current]->read();
02952   _drivers[current]->close();
02953 
02954   END_OF_MED(LOC);
02955 }
02963 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::~FIELD()
02964 {
02965   const char* LOC = " Destructeur FIELD<T, INTERLACING_TAG>::~FIELD()";
02966   BEGIN_OF_MED(LOC);
02967   SCRUTE_MED(this);
02968   if (_value) delete _value; _value=0;
02969   locMap::const_iterator it;
02970   for ( it = _gaussModel.begin();it != _gaussModel.end(); it++ )
02971     delete (*it).second;
02972   _gaussModel.clear();
02973   if(_mesh)
02974     _mesh->removeReference();
02975   _mesh=0;
02976   END_OF_MED(LOC);
02977 }
02978 
02982 template <class T, class INTERLACING_TAG>
02983 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)
02984 {
02985   const char* LOC = "FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)";
02986   BEGIN_OF_MED(LOC);
02987 
02988   _numberOfComponents = NumberOfComponents ;
02989   _componentsTypes.resize(NumberOfComponents);
02990   _componentsNames.resize(NumberOfComponents);
02991   _componentsDescriptions.resize(NumberOfComponents);
02992   _componentsUnits.resize(NumberOfComponents);
02993   _MEDComponentsUnits.resize(NumberOfComponents);
02994   for (int i=0;i<NumberOfComponents;i++) {
02995     _componentsTypes[i] = 0 ;
02996   }
02997   delete _value;
02998   try {
02999     // becarefull about the number of gauss point
03000     _numberOfValues = _support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
03001     MESSAGE_MED(PREFIX_MED <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
03002 
03003     //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
03004     _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
03005 
03006     _isRead = true ;
03007   }
03008   catch (MEDEXCEPTION &ex) {
03009     MESSAGE_MED("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
03010     // OCC 10/03/2006 -- According to the rules defined with help of 
03011     // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
03012     // MEDMEM_Array<> template using INTERLACING_TAG parameter of 
03013     // FIELD template - MSVC++ 2003 compiler generated an error here.
03014     // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
03015     _value = NULL;
03016   }
03017 
03018   SCRUTE_MED(_value);
03019   END_OF_MED(LOC);
03020 }
03021 
03025 template <class T, class INTERLACING_TAG>
03026 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents,
03027                                            const int LengthValue)
03028 {
03029   const char* LOC = "void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)";
03030   BEGIN_OF_MED(LOC);
03031 
03032   _numberOfComponents = NumberOfComponents ;
03033   _componentsTypes.resize(NumberOfComponents);
03034   _componentsNames.resize(NumberOfComponents);
03035   _componentsDescriptions.resize(NumberOfComponents);
03036   _componentsUnits.resize(NumberOfComponents);
03037   _MEDComponentsUnits.resize(NumberOfComponents);
03038   for (int i=0;i<NumberOfComponents;i++) {
03039     _componentsTypes[i] = 0 ;
03040   }
03041 
03042   MESSAGE_MED("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
03043   _numberOfValues = LengthValue ;
03044   delete _value;
03045   //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
03046   _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
03047 
03048   _isRead = true ;
03049 
03050   SCRUTE_MED(_value);
03051   END_OF_MED(LOC);
03052 }
03053 
03057 template <class T, class INTERLACING_TAG>
03058 void FIELD<T, INTERLACING_TAG>::deallocValue()
03059 {
03060   const char* LOC = "void FIELD<T, INTERLACING_TAG>::deallocValue()";
03061   BEGIN_OF_MED(LOC);
03062   _numberOfValues = 0 ;
03063   _numberOfComponents = 0 ;
03064   if (_value != NULL) {
03065     delete _value;
03066     _value = NULL;
03067   }
03068 
03069   END_OF_MED(LOC);
03070 }
03071 
03072 
03073 
03079 // -----------------
03080 // Methodes Inline
03081 // -----------------
03082 
03094 template <class T, class INTERLACING_TAG>
03095 int FIELD<T, INTERLACING_TAG>::addDriver(driverTypes driverType,
03096                                          const string & fileName/*="Default File Name.med"*/,
03097                                          const string & driverName/*="Default Field Name"*/,
03098                                          MED_EN::med_mode_acces access)
03099 {
03100   //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) : ";
03101 
03102   GENDRIVER * driver;
03103 
03104   const char* LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName,const string & driverName,MED_EN::med_mode_acces access) :";
03105   BEGIN_OF_MED(LOC);
03106 
03107   SCRUTE_MED(driverType);
03108 
03109   driver = DRIVERFACTORY::buildDriverForField(driverType,fileName,this,access);
03110 
03111   _drivers.push_back(driver);
03112 
03113   int current = _drivers.size()-1;
03114 
03115   _drivers[current]->setFieldName(driverName);
03116 
03117   END_OF_MED(LOC);
03118 
03119   return current;
03120 }
03121 
03125 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(int index/*=0*/)
03126 {
03127   const char * LOC = "FIELD<T, INTERLACING_TAG>::read(int index=0) : ";
03128   BEGIN_OF_MED(LOC);
03129 
03130   if ( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
03131     _drivers[index]->open();
03132     try
03133     {
03134       _drivers[index]->read();
03135     }
03136     catch ( const MEDEXCEPTION& ex )
03137     {
03138       _drivers[index]->close();
03139       throw ex;
03140     }
03141     _drivers[index]->close();
03142   }
03143   else
03144     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
03145                                      << "The index given is invalid, index must be between  0 and |"
03146                                      << _drivers.size()
03147                                      )
03148                           );
03149   END_OF_MED(LOC);
03150 }
03151 
03156 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(const GENDRIVER & driver )
03157 {
03158   const char* LOC = " FIELD<T, INTERLACING_TAG>::read(const GENDRIVER &) : ";
03159   BEGIN_OF_MED(LOC);
03160 
03161   // For the case where driver does not know about me since it has been created through
03162   // constructor witout parameters: create newDriver knowing me and get missing data
03163   // from driver using merge()
03164   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
03165                                                                          driver.getFileName(),
03166                                                                          this, RDONLY));
03167   newDriver->merge( driver );
03168 
03169   newDriver->open();
03170   try
03171   {
03172     newDriver->read();
03173   }
03174   catch ( const MEDEXCEPTION& ex )
03175   {
03176     newDriver->close();
03177     throw ex;
03178   }
03179   newDriver->close();
03180 
03181   END_OF_MED(LOC);
03182 }
03183 
03188 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(driverTypes driverType, const std::string& filename)
03189 {
03190   const char* LOC = " FIELD<T, INTERLACING_TAG>::read(driverTypes driverType, const std::string& filename) : ";
03191   BEGIN_OF_MED(LOC);
03192 
03193   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driverType, filename,
03194                                                                          this, RDONLY));
03195   newDriver->open();
03196   try
03197   {
03198     newDriver->read();
03199   }
03200   catch ( const MEDEXCEPTION& ex )
03201   {
03202     newDriver->close();
03203     throw ex;
03204   }
03205   newDriver->close();
03206 
03207   END_OF_MED(LOC);
03208 }
03209 
03216 template <class T, class INTERLACING_TAG>
03217 inline int FIELD<T, INTERLACING_TAG>::addDriver (GENDRIVER & driver )
03218 {
03219   int current;
03220 
03221   const char* LOC = "FIELD<T, INTERLACING_TAG>::addDriver(GENDRIVER &) : ";
03222   BEGIN_OF_MED(LOC);
03223 
03224   GENDRIVER * newDriver = 
03225     DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
03226                                        driver.getFileName(), this,
03227                                        driver.getAccessMode());
03228   _drivers.push_back(newDriver);
03229 
03230   current = _drivers.size()-1;
03231   SCRUTE_MED(current);
03232   driver.setId(current);
03233 
03234   newDriver->merge( driver );
03235   newDriver->setId(current);
03236 
03237   return current ;
03238 }
03239 
03243 template <class T, class INTERLACING_TAG>
03244 void FIELD<T, INTERLACING_TAG>::rmDriver (int index/*=0*/)
03245 {
03246   const char * LOC = "FIELD<T, INTERLACING_TAG>::rmDriver (int index=0): ";
03247   BEGIN_OF_MED(LOC);
03248 
03249   if ( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
03250     MESSAGE_MED ("detruire");
03251   }
03252   else
03253     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
03254                                      << "The <index given is invalid, index must be between  0 and  |"
03255                                      << _drivers.size()
03256                                      )
03257                           );
03258 
03259   END_OF_MED(LOC);
03260 }
03261 
03279 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(int index/*=0*/)
03280 {
03281   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
03282   BEGIN_OF_MED(LOC);
03283 
03284   if( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
03285     _drivers[index]->open();
03286     try
03287     {
03288       _drivers[index]->write();
03289     }
03290     catch ( const MEDEXCEPTION& ex )
03291     {
03292       _drivers[index]->close();
03293       throw ex;
03294     }
03295     _drivers[index]->close();
03296   }
03297   else
03298     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
03299                                      << "The index given is invalid, index must be between  0 and |"
03300                                      << _drivers.size()
03301                                      )
03302                           );
03303   END_OF_MED(LOC);
03304 }
03308 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*/)
03309 {
03310   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
03311   BEGIN_OF_MED(LOC);
03312 
03313   // For the case where driver does not know about me since it has been created through
03314   // constructor witout parameters: create newDriver knowing me and get missing data
03315   // from driver using merge()
03316   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
03317                                                                          driver.getFileName(),
03318                                                                          this, WRONLY));
03319   newDriver->merge( driver );
03320   if ( newDriver->getDriverType() == MED_DRIVER )
03321     newDriver->setAccessMode( MED_EN::med_mode_acces( getMedAccessMode( medMode ) ));
03322 
03323   newDriver->open();
03324   try
03325   {
03326     newDriver->write();
03327   }
03328   catch ( const MEDEXCEPTION& ex )
03329   {
03330     newDriver->close();
03331     throw ex;
03332   }
03333   newDriver->close();
03334 
03335   END_OF_MED(LOC);
03336 }
03337 
03341 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*/)
03342 {
03343   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(driverTypes driverType, const std::string& filename) : ";
03344   BEGIN_OF_MED(LOC);
03345 
03346   std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driverType, filename,
03347                                                                          this, WRONLY));
03348   if ( newDriver->getDriverType() == MED_DRIVER )
03349     newDriver->setAccessMode( MED_EN::med_mode_acces( getMedAccessMode( medMode ) ));
03350 
03351   newDriver->open();
03352   try
03353   {
03354     newDriver->write();
03355   }
03356   catch ( const MEDEXCEPTION& ex )
03357   {
03358     newDriver->close();
03359     throw ex;
03360   }
03361   newDriver->close();
03362 
03363   END_OF_MED(LOC);
03364 }
03365 
03371 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
03372 {
03373   const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
03374   BEGIN_OF_MED(LOC);
03375 
03376   if( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
03377     _drivers[index]->openAppend();
03378     if (driverName != "") _drivers[index]->setFieldName(driverName);
03379     try
03380     {
03381       _drivers[index]->writeAppend();
03382     }
03383     catch ( const MEDEXCEPTION& ex )
03384     {
03385       _drivers[index]->close();
03386       throw ex;
03387     }
03388     _drivers[index]->close();
03389   }
03390   else
03391     throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
03392                                      << "The index given is invalid, index must be between  0 and |"
03393                                      << _drivers.size()
03394                                      )
03395                           );
03396   END_OF_MED(LOC);
03397 }
03398 
03405 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(const GENDRIVER & genDriver)
03406 {
03407   const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
03408   BEGIN_OF_MED(LOC);
03409 
03410   for (unsigned int index=0; index < _drivers.size(); index++ )
03411     if ( *_drivers[index] == genDriver ) {
03412       _drivers[index]->openAppend();
03413       try
03414       {
03415         _drivers[index]->writeAppend();
03416       }
03417       catch ( const MEDEXCEPTION& ex )
03418       {
03419         _drivers[index]->close();
03420         throw ex;
03421       }
03422       _drivers[index]->close();
03423     }
03424 
03425   END_OF_MED(LOC);
03426 
03427 }
03428 
03433 template <class T, class INTERLACING_TAG>
03434 bool FIELD<T, INTERLACING_TAG>::getValueOnElement(int eltIdInSup,T* retValues)
03435   const throw (MEDEXCEPTION)
03436 {
03437 
03438   if(eltIdInSup<1)
03439     return false;
03440   if(_support->isOnAllElements())
03441     {
03442       int nbOfEltsThis=_support->getMesh()->getNumberOfElements(_support->getEntity(),MED_EN::MED_ALL_ELEMENTS);
03443       if(eltIdInSup>nbOfEltsThis)
03444         return false;
03445       const T* valsThis=getValue();
03446       for(int j=0;j<_numberOfComponents;j++)
03447         retValues[j]=valsThis[(eltIdInSup-1)*_numberOfComponents+j];
03448       return true;
03449     }
03450   else
03451     {
03452       int nbOfEltsThis=_support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
03453       const int *eltsThis=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
03454       int iThis;
03455       bool found=false;
03456       for(iThis=0;iThis<nbOfEltsThis && !found;)
03457         if(eltsThis[iThis]==eltIdInSup)
03458           found=true;
03459         else
03460           iThis++;
03461       if(!found)
03462         return false;
03463       const T* valsThis=getValue();
03464       for(int j=0;j<_numberOfComponents;j++)
03465         retValues[j]=valsThis[iThis*_numberOfComponents+j];
03466       return true;
03467     }
03468 }
03469 
03475   template <class T, class INTERLACING_TAG>
03476   void FIELD<T, INTERLACING_TAG>::getValueOnPoint(const double* coords, double* output) const throw (MEDEXCEPTION)
03477   {
03478     getValueOnPoints(1, coords, output);
03479   }
03480 
03487   template <class T, class INTERLACING_TAG>
03488   void FIELD<T, INTERLACING_TAG>::getValueOnPoints(int nb_points, const double* coords, double* output) const throw (MEDEXCEPTION)
03489   {
03490     const char* LOC = " FIELD<T, INTERLACING_TAG>::getValueOnPoints(int nb_points, const double* coords, double* output) : ";
03491     // check operation feasibility
03492     if ( !getSupport() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL Support"));
03493     if ( !getSupport()->getMesh() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL Mesh"));
03494     if ( !_value ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL _value"));
03495     if ( getGaussPresence() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not implemeneted for Gauss points"));
03496 
03497     MED_EN::medEntityMesh entity = getSupport()->getEntity();
03498     if ( entity != MED_EN::MED_CELL &&
03499          entity != MED_EN::MED_NODE )
03500       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on CELLs or NODEs"));
03501 
03502     // initialize output value
03503     for ( int j = 0; j < nb_points*getNumberOfComponents(); ++j )
03504       output[j] = 0.0;
03505 
03506     const MEDMEM::MESH* mesh = getSupport()->getMesh()->convertInMESH();
03507     MEDMEM::AutoDeref derefMesh( mesh );
03508 
03509     const double* point = coords;
03510     double* value = output;
03511 
03512     if ( entity == MED_EN::MED_CELL )
03513       {
03514         MEDMEM::PointLocator pLocator (*mesh);
03515         for ( int i = 0; i < nb_points; ++i)
03516           {
03517             // find the cell enclosing the point
03518             std::list<int> cellIds = pLocator.locate( point );
03519             int nbCells = cellIds.size();
03520             if ( nbCells < 1 )
03521               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Point is out of mesh"));
03522 
03523             // retrieve value
03524             std::list<int>::iterator iCell = cellIds.begin();
03525             for ( ; iCell != cellIds.end(); ++iCell )
03526               for ( int j = 1; j <= getNumberOfComponents(); ++j )
03527                 value[j-1] += getValueIJ( *iCell, j ) / nbCells;
03528 
03529             // next point
03530             point += mesh->getSpaceDimension();
03531             value += getNumberOfComponents();
03532           }
03533       }
03534     else // MED_EN::MED_NODE
03535       {
03536         const double * allCoords = mesh->getCoordinates( MED_EN::MED_FULL_INTERLACE );
03537 
03538         MEDMEM::PointLocatorInSimplex pLocator (*mesh);
03539         for ( int i = 0; i < nb_points; ++i)
03540           {
03541             // find nodes of the simplex enclosing the point
03542             std::list<int> nodeIds = pLocator.locate( point );
03543             int nbNodes = nodeIds.size();
03544             if ( nbNodes < 1 )
03545               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Point is out of mesh"));
03546             if ( nbNodes != mesh->getMeshDimension() + 1 )
03547               throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid nb of points of simplex: "<<nbNodes));
03548 
03549             // get coordinates of simplex nodes
03550             std::vector<const double*> nodeCoords( nbNodes );
03551             std::list<int>::iterator iNode = nodeIds.begin();
03552             int n = 0;
03553             for ( ; n < nbNodes; ++iNode, ++n )
03554               nodeCoords[n] = allCoords + (*iNode-1) * mesh->getSpaceDimension();
03555 
03556             // compute wegths of simplex nodes
03557             double nodeWgt[4];
03558             pLocator.getNodeWightsInSimplex( nodeCoords, coords, nodeWgt );
03559 
03560             // retrieve value
03561             for ( n = 0, iNode = nodeIds.begin(); iNode != nodeIds.end(); ++iNode, ++n )
03562               for ( int j = 1; j <= getNumberOfComponents(); ++j )
03563                 value[j-1] += getValueIJ( *iNode, j ) * nodeWgt[ n ];
03564 
03565             // next point
03566             point += mesh->getSpaceDimension();
03567             value += getNumberOfComponents();
03568           }
03569       }
03570   }
03571 
03572 
03577 template <class T, class INTERLACING_TAG>
03578 FIELD<double, FullInterlace>* FIELD<T, INTERLACING_TAG>::getGaussPointsCoordinates() const
03579   throw (MEDEXCEPTION)
03580 {
03581   const char * LOC = "FIELD::getGaussPointsCoordinates() : ";
03582   BEGIN_OF_MED(LOC);
03583 
03584   if (!getSupport())
03585     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
03586 
03587   const MEDMEM::MESH* mesh = getSupport()->getMesh()->convertInMESH();
03588   MEDMEM::AutoDeref derefMesh( mesh );
03589 
03590   const double * coord = mesh->getCoordinates(MED_FULL_INTERLACE);
03591   int spaceDim         = mesh->getSpaceDimension();
03592 
03593   //Init calculator of the gauss point coordinates
03594   INTERP_KERNEL::GaussCoords calculator;
03595   locMap::const_iterator it;
03596 
03597   int nb_type                     = getSupport()->getNumberOfTypes();
03598   int length_values               = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
03599   const medGeometryElement* types = getSupport()->getTypes();
03600   medEntityMesh support_entity    = getSupport()->getEntity();
03601   bool isOnAll                    = getSupport()->isOnAllElements();
03602 
03603   const int* global_connectivity  = 0;
03604   const GAUSS_LOCALIZATION<INTERLACING_TAG>* gaussLock = NULL;
03605 
03606   typedef typename MEDMEM_ArrayInterface<double,INTERLACING_TAG,NoGauss>::Array ArrayCoord;
03607   typedef typename MEDMEM_ArrayInterface<double,INTERLACING_TAG,Gauss>::Array TArrayGauss;
03608 
03609   vector<int>  nbelgeoc, nbgaussgeo;
03610 
03611   nbelgeoc.resize(nb_type+1, 0);
03612   nbgaussgeo.resize(nb_type+1, 0);
03613 
03614   for ( int iType = 0 ; iType < nb_type ; iType++ ) {
03615 
03616     medGeometryElement elem_type = types[iType] ;
03617     if(elem_type == MED_EN::MED_POLYGON && elem_type == MED_EN::MED_POLYHEDRA ) 
03618       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad cell type : "<<MED_EN::geoNames[elem_type]<<" !!! "));
03619 
03620     it = _gaussModel.find(elem_type);
03621 
03622     if(it == _gaussModel.end())
03623       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Gauss localization not defined for "<<MED_EN::geoNames[elem_type]<<" type!!! "));
03624     gaussLock = static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second);
03625 
03626     ArrayCoord coord = gaussLock->getGsCoo();
03627     double* gaussCoord = new double[coord.getNbElem()*coord.getDim()];
03628     int idx = 0;
03629     for( int i = 1 ; i <= coord.getNbElem() ; i++ ) {
03630       for( int j = 1 ; j <= coord.getDim() ; j++ ) {
03631         gaussCoord[idx++] = coord.getIJ(i,j);
03632       }
03633     }
03634 
03635     idx = 0;
03636     ArrayCoord ref = gaussLock->getRefCoo();
03637     double* refCoord = new double[ref.getNbElem()*ref.getDim()];
03638     for( int i = 1 ; i <= ref.getNbElem() ; i++ ) {
03639       for( int j = 1 ; j <= ref.getDim() ; j++ ) {
03640         refCoord[idx++] = ref.getIJ(i,j);
03641       }
03642     }
03643       
03644     INTERP_KERNEL::NormalizedCellType normType;
03645     switch(elem_type) {
03646     case MED_EN::MED_SEG2 : normType = INTERP_KERNEL::NORM_SEG2;break;
03647     case MED_EN::MED_SEG3 : normType = INTERP_KERNEL::NORM_SEG3;break;
03648     default : normType = (INTERP_KERNEL::NormalizedCellType) ((((unsigned long)elem_type/100-2)*10) + ((unsigned long)elem_type%100));
03649       break;
03650     }
03651       
03652     calculator.addGaussInfo(normType,
03653                             elem_type/100,
03654                             gaussCoord,
03655                             gaussLock->getNbGauss(),
03656                             refCoord,
03657                             elem_type%100
03658                             );
03659     // Preapre Info for the gauss array
03660     nbelgeoc   [ iType+1 ] = nbelgeoc[ iType ] + getSupport()->getNumberOfElements(elem_type);
03661     nbgaussgeo [ iType+1 ] = gaussLock->getNbGauss();
03662     
03663     delete [] gaussCoord;
03664     delete [] refCoord;
03665   }
03666 
03667   FIELD<double, FullInterlace>* gpCoord =
03668     new FIELD<double, FullInterlace>(getSupport(),spaceDim);
03669   gpCoord->setName("Gauss Points Coordinates");
03670   gpCoord->setDescription("Gauss Points Coordinates");
03671   
03672   for(int dimId = 1 ; dimId <= spaceDim; dimId++) {
03673     switch(dimId) {
03674     case 1:
03675       gpCoord->setComponentName(dimId,"X");
03676       gpCoord->setComponentDescription(dimId,"X coordinate of the gauss point");
03677       break;
03678     case 2:
03679       gpCoord->setComponentName(dimId,"Y");
03680       gpCoord->setComponentDescription(dimId,"Y coordinate of the gauss point");
03681       break;
03682     case 3:
03683       gpCoord->setComponentName(dimId,"Z");
03684       gpCoord->setComponentDescription(dimId,"Z coordinate of the gauss point");
03685       break;
03686     }
03687     
03688     gpCoord->setMEDComponentUnit(dimId, mesh->getCoordinatesUnits()[dimId-1]);    
03689   }
03690   
03691   gpCoord->setIterationNumber(getIterationNumber());
03692   gpCoord->setOrderNumber(getOrderNumber());
03693   gpCoord->setTime(getTime());
03694 
03695   TArrayGauss *arrayGauss = new TArrayGauss(spaceDim, length_values,
03696                                             nb_type, & nbelgeoc[0], & nbgaussgeo[0]);
03697   gpCoord->setArray(arrayGauss);
03698 
03699 
03700   
03701     
03702   //Calculation of the coordinates  
03703   int index = 1;
03704   for ( int i = 0 ; i < nb_type ; i++ ) {
03705     
03706     medGeometryElement type = types[i] ;
03707     INTERP_KERNEL::NormalizedCellType normType;
03708     switch(type) {
03709     case MED_EN::MED_SEG2 : normType = INTERP_KERNEL::NORM_SEG2;break;
03710     case MED_EN::MED_SEG3 : normType = INTERP_KERNEL::NORM_SEG3;break;
03711     default : normType = (INTERP_KERNEL::NormalizedCellType) ((((unsigned long)type/100-2)*10) + ((unsigned long)type%100));
03712       break;
03713     }
03714     
03715     it = _gaussModel.find(type);
03716     
03717     gaussLock = static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second);
03718     int nb_entity_type = getSupport()->getNumberOfElements(type);
03719     
03720     
03721     if (isOnAll) {
03722       global_connectivity = mesh->getConnectivity(MED_NODAL,support_entity,type);
03723     }
03724     else {
03725       const int * supp_number = getSupport()->getNumber(type);
03726       const int * connectivity = mesh->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
03727       const int * connectivityIndex = mesh->getConnectivityIndex(MED_NODAL,support_entity);
03728       int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
03729       
03730       for (int k_type = 0; k_type<nb_entity_type; k_type++) {
03731         for (int j_ent = 0; j_ent<(type%100); j_ent++) {
03732           global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
03733         }
03734       }
03735       global_connectivity = global_connectivity_tmp;
03736     }
03737 
03738     int nbNodes = (type%100);
03739     double* gCoord = NULL;
03740     int* Ni = NULL; 
03741     
03742     for ( int elem = 0; elem < nb_entity_type; elem++ ) {
03743       int elem_index = nbNodes*elem;
03744       Ni = new int[nbNodes];
03745       for( int idx = 0 ; idx < nbNodes; idx++ ) {
03746         Ni[idx] = global_connectivity[ elem_index+idx ] - 1;
03747       }
03748       
03749       gCoord = calculator.calculateCoords(normType,
03750                                           coord,
03751                                           spaceDim,
03752                                           Ni);
03753       int resultIndex = 0;
03754       for( int k = 0; k < gaussLock->getNbGauss(); k++ ) {
03755         for( int dimId = 1; dimId <= spaceDim; dimId++ ) {
03756           gpCoord->setValueIJK(index,dimId,(k+1),gCoord[resultIndex]);
03757           resultIndex++;
03758         }
03759       }
03760       delete [] gCoord;
03761       delete [] Ni;
03762       index++;
03763     }
03764     if (!isOnAll && type != MED_EN::MED_POLYHEDRA && type != MED_EN::MED_POLYGON) {
03765       delete [] global_connectivity ;
03766     }
03767   }
03768   END_OF_MED(LOC);
03769   return gpCoord;
03770 }
03771 
03777 template <class T, class INTERLACING_TAG>
03778 inline void FIELD<T, INTERLACING_TAG>::setArray(MEDMEM_Array_ * Value)
03779   throw (MEDEXCEPTION)
03780 {
03781   if (NULL != _value) delete _value ;
03782   _value=Value ;
03783 }
03784 
03790 template <class T, class INTERLACING_TAG>
03791 inline MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() const throw (MEDEXCEPTION)
03792 {
03793   const char* LOC = "MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() : ";
03794   BEGIN_OF_MED(LOC);
03795   END_OF_MED(LOC);
03796   return _value ;
03797 }
03798 template <class T,class INTERLACING_TAG>  inline
03799 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array *
03800 FIELD<T, INTERLACING_TAG>::getArrayGauss() const throw (MEDEXCEPTION)
03801 {
03802   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayGauss() : ";
03803   BEGIN_OF_MED(LOC);
03804 
03805   if ( getGaussPresence() )
03806     return static_cast<ArrayGauss *> (_value);
03807   else
03808     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
03809                                  "The field has no Gauss Point"));
03810 
03811   END_OF_MED(LOC);
03812 
03813 }
03814 
03815 template <class T,class INTERLACING_TAG>  inline
03816 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array *
03817 FIELD<T, INTERLACING_TAG>::getArrayNoGauss() const throw (MEDEXCEPTION)
03818 {
03819   const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayNoGauss() : ";
03820   BEGIN_OF_MED(LOC);
03821 
03822   if ( ! getGaussPresence() )
03823     return static_cast < ArrayNoGauss * > (_value);
03824   else
03825     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
03826                                  "The field has Gauss Point"));
03827 
03828   END_OF_MED(LOC);
03829 }
03830 
03831 
03832 template <class T,class INTERLACING_TAG> inline bool
03833 FIELD<T, INTERLACING_TAG>::getGaussPresence() const throw (MEDEXCEPTION)
03834 {
03835   if (_value != NULL)
03836     return _value->getGaussPresence();
03837   else
03838     throw MEDEXCEPTION("FIELD<T, INTERLACING_TAG>::getGaussPresence() const : Can't call getGaussPresence on a null _value");
03839 }
03840 
03845 template <class T, class INTERLACING_TAG>
03846 inline int FIELD<T, INTERLACING_TAG>::getValueLength() const
03847   throw (MEDEXCEPTION)
03848 {
03849   if ( getGaussPresence() )
03850     return static_cast<ArrayGauss *>(_value)->getArraySize() ;
03851   else
03852     return static_cast<ArrayNoGauss *>(_value)->getArraySize() ;
03853 }
03854 
03855 
03883 template <class T, class INTERLACIN_TAG>
03884 inline const T* FIELD<T, INTERLACIN_TAG>::getValue() const throw (MEDEXCEPTION)
03885 {
03886   const char* LOC = "FIELD<T, INTERLACING_TAG>::getValue() : ";
03887   BEGIN_OF_MED(LOC);
03888   if ( getGaussPresence() )
03889     return static_cast<ArrayGauss *>(_value)->getPtr() ;
03890   else
03891     return static_cast<ArrayNoGauss *>(_value)->getPtr() ;
03892 }
03901 template <class T,class INTERLACING_TAG> inline
03902 const T*
03903 FIELD<T,INTERLACING_TAG>::getRow(int i) const throw (MEDEXCEPTION)
03904 {
03905   const char* LOC; LOC = "FIELD<T,INTERLACING_TAG>::getRow(int i) : ";
03906 
03907   int valIndex=-1;
03908   if (_support)
03909     valIndex = _support->getValIndFromGlobalNumber(i);
03910   else
03911     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
03912 
03913   if ( getGaussPresence() )
03914     return static_cast<ArrayGauss *>(_value)->getRow(valIndex) ;
03915   else
03916     return static_cast<ArrayNoGauss *>(_value)->getRow(valIndex) ;
03917 }
03918 
03923 template <class T,class INTERLACING_TAG> inline const T*
03924 FIELD<T,INTERLACING_TAG>::getColumn(int j) const throw (MEDEXCEPTION)
03925 {
03926   if ( getGaussPresence() )
03927     return static_cast<ArrayGauss *>(_value)->getColumn(j) ;
03928   else
03929     return static_cast<ArrayNoGauss *>(_value)->getColumn(j) ;
03930 }
03931 
03940 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJ(int i,int j) const throw (MEDEXCEPTION)
03941 {
03942   const char * LOC = "getValueIJ(..)";
03943   int valIndex=-1;
03944   if (_support)
03945     valIndex = _support->getValIndFromGlobalNumber(i);
03946   else
03947     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
03948 
03949   if ( getGaussPresence() )
03950     return static_cast<ArrayGauss *>(_value)->getIJ(valIndex,j) ;
03951   else
03952     return static_cast<ArrayNoGauss *>(_value)->getIJ(valIndex,j) ;
03953 }
03954 
03962 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION)
03963 {
03964   const char * LOC = "getValueIJK(..)";
03965   int valIndex=-1;
03966   if (_support)
03967     valIndex = _support->getValIndFromGlobalNumber(i);
03968   else
03969     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
03970 
03971   if ( getGaussPresence() )
03972     return static_cast<ArrayGauss *>(_value)->getIJK(valIndex,j,k) ;
03973   else
03974     return static_cast<ArrayNoGauss *>(_value)->getIJK(valIndex,j,k) ;
03975 }
03981 template <class T, class INTERLACIN_TAG>
03982 inline int FIELD<T, INTERLACIN_TAG>::getValueByTypeLength(int t) const throw (MEDEXCEPTION)
03983 {
03984   const char * LOC ="getValueByTypeLength() : ";
03985   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
03986     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
03987 
03988   if ( getGaussPresence() ) {
03989     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
03990     if (  t < 1 || t > array->getNbGeoType() )
03991       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
03992     return array->getLengthOfType( t );
03993   }
03994   else {
03995     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
03996     if (  t < 1 || t > array->getNbGeoType() )
03997       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
03998     return array->getLengthOfType( t );
03999   }
04000 }
04001 
04005 template <class T, class INTERLACIN_TAG>
04006 inline const T* FIELD<T, INTERLACIN_TAG>::getValueByType(int t) const throw (MEDEXCEPTION)
04007 {
04008   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
04009     throw MEDEXCEPTION(LOCALIZED("getValueByType() : not MED_NO_INTERLACE_BY_TYPE field" ));
04010 
04011   if ( getGaussPresence() ) {
04012     ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
04013     return array->getPtr() + array->getIndex( t );
04014   }
04015   else {
04016     ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
04017     return array->getPtr() + array->getIndex( t );
04018   }
04019 }
04020 
04024 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJByType(int i,int j, int t) const throw (MEDEXCEPTION)
04025 {
04026   const char * LOC = "getValueIJByType(..)";
04027   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
04028     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
04029     
04030   if ( getGaussPresence() )
04031     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJByType(i,j,t) ;
04032   else
04033     return static_cast<ArrayNoByType *>(_value)->getIJByType(i,j,t) ;
04034 }
04035 
04039 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION)
04040 {
04041   const char * LOC = "getValueIJKByType(..)";
04042   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
04043     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
04044 
04045   if ( getGaussPresence() )
04046     return static_cast<ArrayNoByTypeGauss *>(_value)->getIJKByType(i,j,k,t) ;
04047   else
04048     return static_cast<ArrayNoByType      *>(_value)->getIJKByType(i,j,k,t) ;
04049 }
04050 
04051 
04052 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGeometricTypes() const throw (MEDEXCEPTION)
04053 {
04054   const char * LOC = "getNumberOfGeometricTypes(..)";
04055   BEGIN_OF_MED(LOC);
04056   if (_support)
04057     return _support->getNumberOfTypes();
04058   else
04059     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
04060   END_OF_MED(LOC);
04061 }
04062 
04069 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> &
04070 FIELD<T,INTERLACING_TAG>::getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
04071 {
04072   const char * LOC ="getGaussLocalization(MED_EN::medGeometryElement geomElement) : ";
04073   const GAUSS_LOCALIZATION_ * locPtr=0;
04074 
04075   locMap::const_iterator it;
04076   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
04077         locPtr = (*it).second;
04078         return  *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
04079   }
04080   else
04081     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
04082 
04083 }
04084 
04085 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> *
04086 FIELD<T,INTERLACING_TAG>::getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
04087 {
04088   const char * LOC ="getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) : ";
04089   const GAUSS_LOCALIZATION_ * locPtr=0;
04090 
04091   locMap::const_iterator it;
04092   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
04093         locPtr = (*it).second;
04094         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
04095   }
04096   else
04097     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
04098 
04099 }
04100 
04104 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION_ *
04105 FIELD<T,INTERLACING_TAG>::getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
04106 {
04107   const char * LOC ="getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) : ";
04108 
04109   locMap::const_iterator it;
04110   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
04111         return (*it).second;
04112   }
04113   else
04114     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type: "<< geomElement ));
04115 
04116 }
04117 
04121 template <class T,class INTERLACING_TAG> void
04122 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc)
04123 {
04124   locMap::iterator it = _gaussModel.find(geomElement);
04125   if ( it != _gaussModel.end() ) {
04126     delete it->second;
04127     it->second = gaussloc;
04128   }
04129   else {
04130     _gaussModel[ geomElement ] = gaussloc;
04131   }
04132 }
04133 
04134 
04135 template <class T,class INTERLACING_TAG> void
04136 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG>& gaussloc)
04137 {
04138   locMap::iterator it = _gaussModel.find(geomElement);
04139   if ( it != _gaussModel.end() ) {
04140     delete it->second;
04141     it->second = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
04142   }
04143   else {
04144     _gaussModel[ geomElement ] = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
04145   }
04146 }
04147 
04156 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) const
04157   throw (MEDEXCEPTION)
04158 {
04159   const char * LOC ="getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
04160   const GAUSS_LOCALIZATION_ * locPtr=0;
04161 
04162   locMap::const_iterator it;
04163   if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
04164         locPtr = (*it).second;
04165         return  static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr)->getNbGauss();
04166   }
04167   else
04168     if (_support)
04169       try {
04170         if ( _support->getNumberOfElements(geomElement) ) return 1;
04171       } catch ( MEDEXCEPTION & ex) {
04172         throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<< "GeometricType not found !" )) ;
04173       }
04174     else
04175       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
04176 
04177   throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Should never execute this!" ));
04178 
04179 }
04180 
04188 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints() const
04189   throw (MEDEXCEPTION)
04190 {
04191   const char * LOC ="const int * getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
04192 
04193   if (_value)
04194     if ( getGaussPresence() ) {
04195       return static_cast<ArrayGauss *>(_value)->getNbGaussGeo()+1;
04196     } else
04197       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"value hasn't Gauss points " ));
04198 
04199     else
04200       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Value not defined" ));
04201 }
04202 
04208 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNbGaussI(int i) const throw (MEDEXCEPTION)
04209 {
04210   const char * LOC = "getNbGaussI(..)";
04211 
04212   int valIndex=-1;
04213   if (_support)
04214     valIndex = _support->getValIndFromGlobalNumber(i);
04215   else
04216     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
04217 
04218   if (_value)
04219    if ( getGaussPresence() )
04220      return static_cast<ArrayGauss *>(_value)->getNbGauss(valIndex) ;
04221    else
04222      return static_cast<ArrayNoGauss *>(_value)->getNbGauss(valIndex) ;
04223  else
04224    throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_value not defined" ));
04225 }
04231 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfElements() const throw (MEDEXCEPTION)
04232 {
04233   const char * LOC = "getNumberOfElements(..)";
04234   BEGIN_OF_MED(LOC);
04235   if (_support)
04236     return _support->getNumberOfElements();
04237   else
04238     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
04239   END_OF_MED(LOC);
04240 }
04241 
04242 template <class T,class INTERLACING_TAG> const MED_EN::medGeometryElement  * FIELD<T,INTERLACING_TAG>::getGeometricTypes()  const throw (MEDEXCEPTION)
04243 {
04244   const char * LOC = "getGeometricTypes(..)";
04245   BEGIN_OF_MED(LOC);
04246   if (_support)
04247     return _support->getTypes();
04248   else
04249     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
04250   END_OF_MED(LOC);
04251 }
04252 template <class T,class INTERLACING_TAG> bool  FIELD<T,INTERLACING_TAG>::isOnAllElements() const throw (MEDEXCEPTION)
04253 {
04254   const char * LOC = "isOnAllElements(..)";
04255   BEGIN_OF_MED(LOC);
04256   if (_support)
04257     return _support->isOnAllElements();
04258   else
04259     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
04260   END_OF_MED(LOC);
04261 }
04262 
04263 
04276 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValue( T* value) throw (MEDEXCEPTION) 
04277 {
04278   if ( getGaussPresence() )
04279     static_cast<ArrayGauss *>(_value)->setPtr(value) ;
04280   else
04281     static_cast<ArrayNoGauss *>(_value)->setPtr(value) ;
04282 }
04283 
04288 template <class T,class INTERLACING_TAG>
04289 inline void FIELD<T,INTERLACING_TAG>::setRow( int i, T* value) throw (MEDEXCEPTION)
04290 {
04291   const char * LOC = "FIELD<T,INTERLACING_TAG>::setRow(int i, T* value) : ";
04292   int valIndex=i;
04293   if (_support)
04294     valIndex = _support->getValIndFromGlobalNumber(i);
04295   else
04296     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
04297 
04298   if ( getGaussPresence() )
04299     static_cast<ArrayGauss *>(_value)->setRow(valIndex, value) ;
04300   else
04301     static_cast<ArrayNoGauss *>(_value)->setRow(valIndex, value) ;
04302 }
04303 
04308 template <class T,class INTERLACING_TAG>
04309 inline void FIELD<T,INTERLACING_TAG>::setColumn( int j, T* value) throw (MEDEXCEPTION)
04310 {
04311   if ( getGaussPresence() )
04312     static_cast<ArrayGauss *>(_value)->setColumn(j, value) ;
04313   else
04314     static_cast<ArrayNoGauss *>(_value)->setColumn(j, value) ;
04315 }
04316 
04320 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) throw (MEDEXCEPTION) 
04321 {
04322   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
04323   int valIndex=-1;
04324   if (_support)
04325     valIndex = _support->getValIndFromGlobalNumber(i);
04326   else
04327     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
04328 
04329   if ( getGaussPresence() )
04330     static_cast<ArrayGauss *>(_value)->setIJ(valIndex,j,value) ;
04331   else
04332     static_cast<ArrayNoGauss *>(_value)->setIJ(valIndex,j,value) ;
04333 }
04334 
04338 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION) 
04339 {
04340   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
04341   int valIndex=-1;
04342   if (_support)
04343     valIndex = _support->getValIndFromGlobalNumber(i);
04344   else
04345     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
04346 
04347   if ( getGaussPresence() )
04348     static_cast<ArrayGauss *>(_value)->setIJK(valIndex,j,k,value) ;
04349   else
04350     static_cast<ArrayNoGauss *>(_value)->setIJK(valIndex,j,k,value) ;
04351 }
04352 
04356 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION) 
04357 {
04358   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) : ";
04359   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
04360     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
04361 
04362   if ( getGaussPresence() )
04363     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJByType(i,j,t,value) ;
04364   else
04365     return static_cast<ArrayNoByType *>(_value)->setIJByType(i,j,t,value) ;
04366 }
04367 
04371 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) 
04372 {
04373   const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int t, int k, T value) : ";
04374   if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
04375     throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
04376 
04377   if ( getGaussPresence() )
04378     return static_cast<ArrayNoByTypeGauss *>(_value)->setIJKByType(i,j,k,t,value) ;
04379   else
04380     return static_cast<ArrayNoByType *>(_value)->setIJKByType(i,j,k,t,value) ;
04381 }
04384 /*
04385   METHODS
04386 */
04387 
04393 template <class T, class INTERLACING_TAG>
04394 void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION)
04395 {
04396   const char * LOC = "void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) : ";
04397   int i,j;
04398   if (_support == (SUPPORT *) NULL)
04399       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No Support defined."));
04400 
04401   const GMESH * mesh = _support->getMesh();
04402   int spaceDim = mesh->getSpaceDimension();
04403   const double * coord;
04404 
04405   const double * bary;
04406   FIELD<double,FullInterlace> * barycenterField=0;
04407 
04408   double ** xyz=new double* [spaceDim];
04409   bool deallocateXyz=false;
04410   if(_support->getEntity()==MED_EN::MED_NODE)
04411     {
04412       const MESH * unstructured = _support->getMesh()->convertInMESH();
04413       if (_support->isOnAllElements())
04414         {
04415           coord=unstructured->getCoordinates(MED_EN::MED_NO_INTERLACE);
04416           for(i=0; i<spaceDim; i++)
04417             xyz[i]=(double *)coord+i*_numberOfValues;
04418         }
04419       else
04420         {
04421           coord = unstructured->getCoordinates(MED_EN::MED_FULL_INTERLACE);
04422           const int * nodesNumber=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
04423           for(i=0; i<spaceDim; i++)
04424             xyz[i]=new double[_numberOfValues];
04425           deallocateXyz=true;
04426           for(i=0;i<_numberOfValues;i++)
04427             {
04428               for(j=0;j<spaceDim;j++)
04429                 xyz[j][i]=coord[(nodesNumber[i]-1)*spaceDim+j];
04430             }
04431         }
04432       unstructured->removeReference();
04433     }
04434   else
04435     {
04436       barycenterField = mesh->getBarycenter(_support);
04437       bary = barycenterField->getValue();
04438       for(i=0; i<spaceDim; i++)
04439         xyz[i]=new double[_numberOfValues];
04440       deallocateXyz=true;
04441       for(i=0;i<_numberOfValues;i++) {
04442         for(j=0;j<spaceDim;j++)
04443           xyz[j][i]=bary[i*spaceDim+j];
04444       }
04445     }
04446   T* valsToSet=(T*)getValue();
04447   double *temp=new double[spaceDim];
04448   for(i=0;i<_numberOfValues;i++)
04449   {
04450     for(j=0;j<spaceDim;j++)
04451       temp[j]=xyz[j][i];
04452     f(temp,valsToSet+i*_numberOfComponents);
04453   }
04454   delete [] temp;
04455   if(barycenterField)
04456     delete barycenterField;
04457   if(deallocateXyz)
04458     for(j=0;j<spaceDim;j++)
04459       delete [] xyz[j];
04460   delete [] xyz;
04461 }
04462 
04468 template <class T, class INTERLACING_TAG>
04469 FIELD<T,INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::execFunc(int nbOfComponents,
04470                                                               myFuncType2 f) throw (MEDEXCEPTION)
04471 {
04472   FIELD<T,INTERLACING_TAG> *ret=new FIELD<T,INTERLACING_TAG>(_support,nbOfComponents);
04473   const T* valsInput=getValue();
04474   T* valsOutPut=(T*)ret->getValue();
04475   int i;
04476   for(i=0;i<_numberOfValues;i++)
04477     f(valsInput+i*_numberOfComponents,valsOutPut+i*nbOfComponents);
04478   return ret;
04479 }
04480 
04481 }//End namespace MEDMEM
04482 
04483 #endif /* FIELD_HXX */