00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
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
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;
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
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
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
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 }
00809
00811
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
00848 Array *_value ;
00849
00850
00851
00852 GMESH* _mesh;
00853
00854
00855 T _vmin;
00856 T _vmax;
00857
00858 map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> _gaussModel;
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
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
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
01078 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE);
01079 FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
01080
01081
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
01122 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
01123 FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
01124
01125
01126 ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
01127 FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
01128
01129 try
01130 {
01131
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
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
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
01221 FIELD_::operator=(m);
01222
01223 _value = m._value;
01224
01225
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);
01291
01292
01293 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01294 result->_operationInitialize(*this,m,"+");
01295 result->_add_in_place(*this,m);
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);
01311
01312 const T* value1=m.getValue();
01313
01314
01315 T * value=const_cast<T *> (getValue());
01316 const int size=getNumberOfValues()*getNumberOfComponents();
01317 const T* endV=value+size;
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);
01336
01337
01338 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01339 m.getNumberOfComponents());
01340 result->_operationInitialize(m,n,"+");
01341 result->_add_in_place(m,n);
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);
01355
01356
01357 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01358 m.getNumberOfComponents());
01359 result->_operationInitialize(m,n,"+");
01360 result->_add_in_place(m,n);
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);
01392
01393
01394 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01395 result->_operationInitialize(*this,m,"-");
01396 result->_sub_in_place(*this,m);
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
01409 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01410
01411 result->setName("- "+getName());
01412 result->setComponentsNames(getComponentsNames());
01413
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
01423 T * value=const_cast<T *> (result->getValue());
01424 const int size=getNumberOfValues()*getNumberOfComponents();
01425 const T* endV=value+size;
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);
01443
01444 const T* value1=m.getValue();
01445
01446
01447 T * value=const_cast<T *> (getValue());
01448 const int size=getNumberOfValues()*getNumberOfComponents();
01449 const T* endV=value+size;
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
01466 T * value=const_cast<T *> (getValue());
01467
01468 const int size=getNumberOfValues()*getNumberOfComponents();
01469
01470 if (size>0)
01471 {
01472 value+=icomp-1;
01473 const T* lastvalue=value+size;
01474
01475 for(;value!=lastvalue; value+=getNumberOfComponents())
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);
01491
01492
01493 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01494 m.getNumberOfComponents());
01495 result->_operationInitialize(m,n,"-");
01496 result->_sub_in_place(m,n);
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);
01510
01511
01512 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01513 m.getNumberOfComponents());
01514 result->_operationInitialize(m,n,"-");
01515 result->_sub_in_place(m,n);
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);
01547
01548
01549 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01550 result->_operationInitialize(*this,m,"*");
01551 result->_mul_in_place(*this,m);
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);
01567
01568 const T* value1=m.getValue();
01569
01570
01571 T * value=const_cast<T *> (getValue());
01572 const int size=getNumberOfValues()*getNumberOfComponents();
01573 const T* endV=value+size;
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);
01594
01595
01596 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01597 m.getNumberOfComponents());
01598 result->_operationInitialize(m,n,"*");
01599 result->_mul_in_place(m,n);
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);
01613
01614
01615 FIELD<T, INTERLACING_TAG>* result = new FIELD<T,INTERLACING_TAG>(m.getSupport(),
01616 m.getNumberOfComponents());
01617 result->_operationInitialize(m,n,"*");
01618 result->_mul_in_place(m,n);
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);
01650
01651
01652 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
01653 try
01654 {
01655 result->_operationInitialize(*this,m,"/");
01656 result->_div_in_place(*this,m);
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);
01679
01680 const T* value1=m.getValue();
01681
01682
01683 T * value=const_cast<T *> (getValue());
01684 const int size=getNumberOfValues()*getNumberOfComponents();
01685 const T* endV=value+size;
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);
01706
01707
01708 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01709 m.getNumberOfComponents());
01710 try
01711 {
01712 result->_operationInitialize(m,n,"/");
01713 result->_div_in_place(m,n);
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);
01732
01733
01734 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
01735 m.getNumberOfComponents());
01736 try
01737 {
01738 result->_operationInitialize(m,n,"/");
01739 result->_div_in_place(m,n);
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
01768
01769 setName(m.getName()+" "+Op+" "+n.getName());
01770 setComponentsNames(m.getComponentsNames());
01771 setComponentsDescriptions(m.getComponentsDescriptions());
01772 setMEDComponentsUnits(m.getMEDComponentsUnits());
01773
01774
01775
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
01797 const T* value1=m.getValue();
01798 const T* value2=n.getValue();
01799
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
01821 const T* value1=m.getValue();
01822 const T* value2=n.getValue();
01823
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
01845 const T* value1=m.getValue();
01846 const T* value2=n.getValue();
01847
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
01869 const T* value1=m.getValue();
01870 const T* value2=n.getValue();
01871
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 ) {
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();
01897 const int size=getNumberOfValues()*getNumberOfComponents();
01898 if (size <= 0)
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;
01906 const T* pMax=value;
01907 const T* pMin=value;
01908
01909
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;
01919 T Min= *pMin>(T) 0 ? *pMin : -*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();
01929 const int size=getNumberOfValues()*getNumberOfComponents();
01930 if (size <= 0)
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;
01938
01939 T result((T)0);
01940 for( ; value!=lastvalue ; ++value)
01941 result += (*value) * (*value);
01942
01943 return std::sqrt(double(result));
01944 }
01945
01946
01947
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();
01954 const int size=getNumberOfValues()*getNumberOfComponents();
01955 const T* lastvalue=value+size;
01956
01957 if (size <= 0){
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();
01967 vmin=MinMax<T>::getMax();
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();
01991 const int size=getNumberOfValues()*getNumberOfComponents();
01992 const T* lastvalue=value+size;
01993
01994 if (size <= 0){
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
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
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
02076 C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,typ,MED_ALL_ELEMENTS);
02077 iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,typ);
02078
02079 revC = getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
02080 indC = getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
02081
02082 coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
02083
02084 NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
02085
02086 FIELD<double, FullInterlace>* barycenter = getSupport()->getMesh()->getBarycenter(getSupport());
02087
02088
02089 for (int i = 1; i < NumberOf + 1; i++) {
02090
02091
02092 set <int> listElements;
02093 set <int>::iterator elemIt;
02094 listElements.clear();
02095
02096
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
02101 int c = revC[k-1];
02102
02103 if (c != i)
02104 listElements.insert(c);
02105 }
02106 }
02107
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
02113 double val = getValueIJ(i,1);
02114 double grad = 0.;
02115
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
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
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
02138 revC=getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
02139 indC=getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
02140
02141 coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
02142
02143
02144 NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
02145 for (int i=1; i<NumberOf+1; i++){
02146
02147 set <int> listNodes;
02148 set <int>::iterator nodeIt ;
02149 listNodes.clear();
02150 for(int j=indC[i-1];j<indC[i];j++){
02151
02152 int c=revC[j-1];
02153
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
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
02164 double val = getValueIJ(i,1);
02165 double grad = 0.;
02166
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
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
02250 T * value=const_cast<T *> (getValue());
02251 const int size=getNumberOfValues()*getNumberOfComponents();
02252
02253 if (size>0)
02254 {
02255 const T* lastvalue=value+size;
02256 for(;value!=lastvalue; ++value)
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
02285 T * value=const_cast<T *> (getValue());
02286 const int size=getNumberOfValues()*getNumberOfComponents();
02287
02288 if (size>0)
02289 {
02290 const T* lastvalue=value+size;
02291 for(;value!=lastvalue; ++value)
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);
02317 else
02318 FIELD_::_deepCheckFieldCompatibility(m, n, false);
02319
02320
02321
02322
02323 const int numberOfElements=m.getNumberOfValues();
02324 const int NumberOfComponents=m.getNumberOfComponents();
02325
02326
02327
02328
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();
02337 const T* value2=n.getValue();
02338
02339 T * value=const_cast<T *> (result->getValue());
02340
02341 const T* lastvalue=value+numberOfElements;
02342 for ( ; value!=lastvalue ; ++value )
02343 {
02344 *value=(T)0;
02345 const T* endofRow=value1+NumberOfComponents;
02346 for ( ; value1 != endofRow; ++value1, ++value2)
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, true);
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)
02367 p_field_size=_getFieldSize();
02368 else
02369 p_field_size->addReference();
02370
02371 const double* vol=p_field_size->getValue();
02372
02373
02374
02375
02376 double integrale=0.0;
02377 double totVol=0.0;
02378
02379 if ( getSupport()->getEntity() == MED_NODE )
02380 {
02381
02382
02383
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
02390 double curCellValue = 0;
02391 try {
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();
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 {
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();
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, true);
02450 const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
02451 if(!p_field_volume)
02452 p_field_size=_getFieldSize();
02453 else
02454 p_field_size->addReference();
02455
02456 const double* vol=p_field_size->getValue();
02457 const double* lastvol=vol+getNumberOfValues();
02458
02459
02460 double integrale=0.0;
02461 double totVol=0.0;
02462
02463 if ( getSupport()->getEntity() == MED_NODE )
02464 {
02465
02466
02467
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
02475 int nbNodes = iC[i+1]-iC[i];
02476 vector< double > curCellValue( nbComp, 0 );
02477 try {
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)
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) {
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 {
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();
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);
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)
02552 p_field_size=_getFieldSize();
02553 else
02554 p_field_size->addReference();
02555
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();
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 {
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();
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);
02599 const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
02600 if(!p_field_volume)
02601 p_field_size=_getFieldSize();
02602 else
02603 p_field_size->addReference();
02604
02605 const double* vol = p_field_size->getValue();
02606 const double* lastvol = vol+getNumberOfValues();
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)
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) {
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 {
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
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
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
02683 const FIELD<double, FullInterlace> * cellSize=_getFieldSize(subSupport);
02684 const double* size = cellSize->getValue();
02685 const double* lastSize = size + nbElems;
02686
02687 const T* value = getValue();
02688
02689
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
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
02717 PointerOf<int> index;
02718 if ( _support->isOnAllElements() )
02719 {
02720 index.set( subNums );
02721 }
02722 else
02723 {
02724
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() )
02735 {
02736 index[iSub] = 0;
02737 break;
02738 }
02739 else if ( subNums[iSub] == myNums[i] )
02740 index[iSub] = ++i;
02741 else
02742 allNumsFound = (index[iSub] = 0);
02743 }
02744 if ( iSub != nbElems || !allNumsFound )
02745 {
02746
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
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
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
02867 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
02868 FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
02869
02870
02871 ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
02872 FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
02873
02874 _support = Support;
02875
02876 if(_support)
02877 _support->addReference();
02878
02879
02880
02881
02882
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
02929 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
02930 FIELD_::_valueType = SET_VALUE_TYPE<T>::_valueType;
02931
02932
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
02938
02939
02940
02941
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
03000 _numberOfValues = _support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
03001 MESSAGE_MED(PREFIX_MED <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
03002
03003
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
03011
03012
03013
03014
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
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
03081
03082
03094 template <class T, class INTERLACING_TAG>
03095 int FIELD<T, INTERLACING_TAG>::addDriver(driverTypes driverType,
03096 const string & fileName,
03097 const string & driverName,
03098 MED_EN::med_mode_acces access)
03099 {
03100
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)
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
03162
03163
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)
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)
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)
03309 {
03310 const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
03311 BEGIN_OF_MED(LOC);
03312
03313
03314
03315
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)
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, 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
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
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
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
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
03530 point += mesh->getSpaceDimension();
03531 value += getNumberOfComponents();
03532 }
03533 }
03534 else
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
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
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
03557 double nodeWgt[4];
03558 pLocator.getNodeWightsInSimplex( nodeCoords, coords, nodeWgt );
03559
03560
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
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
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
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
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
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 }
04482
04483 #endif