13 #include "utilities.h"
14 #include "MEDMEM_Exception.hxx"
15 #include "MEDMEM_define.hxx"
17 #include "MEDMEM_Support.hxx"
18 #include "MEDMEM_Unit.hxx"
19 #include "MEDMEM_Array.hxx"
20 #include "MEDMEM_GenDriver.hxx"
22 #include "MEDMEM_MedFieldDriver.hxx"
23 #include "MEDMEM_MedMedDriver.hxx"
25 #include "MEDMEM_VtkFieldDriver.hxx"
26 #include "MEDMEM_VtkMedDriver.hxx"
28 using namespace MED_EN;
32 This class contains all the informations related with a template class FIELD :
33 - Components descriptions
34 - Time step description
35 - Location of the values (a SUPPORT class)
39 class FIELD_ // GENERIC POINTER TO a template <class T> class FIELD
59 Pointer to the support the field deals with.
62 const SUPPORT * _support ;
66 Number of field's components.
69 int _numberOfComponents ;
72 Number of field's values.
79 Array of size _numberOfComponents. /n
80 (constant, scalar, vector, tensor)/n
81 We could use an array of integer to store
84 - space dimension for vector,/n
85 - space dimension square for tensor./n
86 So numbers of values per entities would be
87 sum of _componentsTypes array.
89 Not implemented yet! All type are scalar !
92 int * _componentsTypes ;
95 Array of size _numberOfComponents
96 storing components names if any.
99 string * _componentsNames;
102 Array of size _numberOfComponents
103 storing components descriptions if any.
106 string * _componentsDescriptions;
109 Array of size _numberOfComponents
110 storing components units if any.
113 UNIT * _componentsUnits;
116 Array of size _numberOfComponents
117 storing components units if any.
120 string * _MEDComponentsUnits;
121 int _iterationNumber ;
125 // _valueType should be a static const. Here is an initialization exemple
126 // template < classType T > struct SET_VALUE_TYPE { static const med_type_champ _valueType = 0; }
127 // template < > struct SET_VALUE_TYPE<double> { static const med_type_champ _valueType = MED_EN::MED_REEL64; }
128 // template < > struct SET_VALUE_TYPE<int> { static const med_type_champ _valueType = MED_EN::MED_INT32; }
129 // static const med_type_champ _valueType = SET_VALUE_TYPE <T>::_valueType;
130 med_type_champ _valueType ;
132 vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
133 static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n ) throw (MEDEXCEPTION);
134 void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL) const throw (MEDEXCEPTION);
135 FIELD<double>* _getFieldSize() const;
139 friend class MED_MED_RDONLY_DRIVER;
140 friend class MED_MED_WRONLY_DRIVER;
141 friend class MED_MED_RDWR_DRIVER;
143 friend class VTK_MED_DRIVER;
152 FIELD_(const SUPPORT * Support, const int NumberOfComponents);
156 FIELD_(const FIELD_ &m);
163 // virtual void setIterationNumber (int IterationNumber);
164 // virtual void setOrderNumber (int OrderNumber);
165 // virtual void setFieldName (string& fieldName);
167 virtual void rmDriver(int index);
168 virtual int addDriver(driverTypes driverType,
169 const string & fileName="Default File Name.med",
170 const string & driverFieldName="Default Field Nam") ;
171 virtual int addDriver( GENDRIVER & driver);
172 virtual void read (const GENDRIVER &);
173 virtual void read(int index=0);
174 virtual void openAppend( void );
175 virtual void write(const GENDRIVER &);
176 virtual void write(int index=0, const string & driverName="");
178 virtual void writeAppend(const GENDRIVER &);
179 virtual void writeAppend(int index=0, const string & driverName="");
181 inline void setName(const string Name);
182 inline string getName() const;
183 inline void setDescription(const string Description);
184 inline string getDescription() const;
185 inline const SUPPORT * getSupport() const;
186 inline void setSupport(const SUPPORT * support);
187 inline void setNumberOfComponents(const int NumberOfComponents);
188 inline int getNumberOfComponents() const;
189 inline void setNumberOfValues(const int NumberOfValues);
190 inline int getNumberOfValues() const;
191 // inline void setComponentType(int *ComponentType);
192 // inline int * getComponentType() const;
193 // inline int getComponentTypeI(int i) const;
194 inline void setComponentsNames(const string * ComponentsNames);
195 inline void setComponentName(int i, const string ComponentName);
196 inline const string * getComponentsNames() const;
197 inline string getComponentName(int i) const;
198 inline void setComponentsDescriptions(const string * ComponentsDescriptions);
199 inline void setComponentDescription(int i, const string ComponentDescription);
200 inline const string * getComponentsDescriptions() const;
201 inline string getComponentDescription(int i) const;
203 // provisoire : en attendant de regler le probleme des unites !
204 inline void setComponentsUnits(const UNIT * ComponentsUnits);
205 inline const UNIT * getComponentsUnits() const;
206 inline const UNIT * getComponentUnit(int i) const;
207 inline void setMEDComponentsUnits(const string * MEDComponentsUnits);
208 inline void setMEDComponentUnit(int i, const string MEDComponentUnit);
209 inline const string * getMEDComponentsUnits() const;
210 inline string getMEDComponentUnit(int i) const;
212 inline void setIterationNumber(int IterationNumber);
213 inline int getIterationNumber() const;
214 inline void setTime(double Time);
215 inline double getTime() const;
216 inline void setOrderNumber(int OrderNumber);
217 inline int getOrderNumber() const;
219 inline void setValueType (const med_type_champ ValueType) ;
220 inline med_type_champ getValueType () const;
224 // ---------------------------------
225 // Implemented Methods : constructor
226 // ---------------------------------
234 inline void FIELD_::setName(const string Name)
241 inline string FIELD_::getName() const
246 Set FIELD description.
248 inline void FIELD_::setDescription(const string Description)
250 _description=Description;
253 Get FIELD description.
255 inline string FIELD_::getDescription() const
260 Set FIELD number of components.
262 inline void FIELD_::setNumberOfComponents(int NumberOfComponents)
264 _numberOfComponents=NumberOfComponents;
267 Get FIELD number of components.
269 inline int FIELD_::getNumberOfComponents() const
271 return _numberOfComponents ;
274 Set FIELD number of values.
276 It must be the same than in the associated SUPPORT object.
278 inline void FIELD_::setNumberOfValues(int NumberOfValues)
280 _numberOfValues=NumberOfValues;
283 Get FIELD number of value.
285 inline int FIELD_::getNumberOfValues() const
287 return _numberOfValues ;
290 // inline void FIELD_::setComponentType(int *ComponentType)
292 // _componentsTypes=ComponentType ;
294 // inline int * FIELD_::getComponentType() const
296 // return _componentsTypes ;
298 // inline int FIELD_::getComponentTypeI(int i) const
300 // return _componentsTypes[i-1] ;
304 Set FIELD components names.
306 Duplicate the ComponentsNames string array to put components names in
307 FIELD. ComponentsNames size must be equal to number of components.
309 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
311 if (NULL == _componentsNames)
312 _componentsNames = new string[_numberOfComponents] ;
313 for (int i=0; i<_numberOfComponents; i++)
314 _componentsNames[i]=ComponentsNames[i] ;
317 Set FIELD i^th component name.
319 i must be >=1 and <= number of components.
321 inline void FIELD_::setComponentName(int i, const string ComponentName)
323 _componentsNames[i-1]=ComponentName ;
326 Get a reference to the string array which contain the components names.
328 This Array size is equal to number of components
330 inline const string * FIELD_::getComponentsNames() const
332 return _componentsNames ;
335 Get the name of the i^th component.
337 inline string FIELD_::getComponentName(int i) const
339 return _componentsNames[i-1] ;
342 Set FIELD components descriptions.
344 Duplicate the ComponentsDescriptions string array to put components
345 descriptions in FIELD.
346 ComponentsDescriptions size must be equal to number of components.
348 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
350 if (NULL == _componentsDescriptions)
351 _componentsDescriptions = new string[_numberOfComponents] ;
352 for (int i=0; i<_numberOfComponents; i++)
353 _componentsDescriptions[i]=ComponentsDescriptions[i] ;
356 Set FIELD i^th component description.
358 i must be >=1 and <= number of components.
360 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
362 _componentsDescriptions[i-1]=ComponentDescription ;
365 Get a reference to the string array which contain the components descriptions.
367 This Array size is equal to number of components
369 inline const string * FIELD_::getComponentsDescriptions() const
371 return _componentsDescriptions ;
374 Get the description of the i^th component.
376 inline string FIELD_::getComponentDescription(int i) const
378 return _componentsDescriptions[i-1];
383 Set FIELD components UNIT.
385 Duplicate the ComponentsUnits UNIT array to put components
387 ComponentsUnits size must be equal to number of components.
389 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
391 if (NULL == _componentsUnits)
392 _componentsUnits = new UNIT[_numberOfComponents] ;
393 for (int i=0; i<_numberOfComponents; i++)
394 _componentsUnits[i]=ComponentsUnits[i] ;
397 Get a reference to the UNIT array which contain the components units.
399 This Array size is equal to number of components
401 inline const UNIT * FIELD_::getComponentsUnits() const
403 return _componentsUnits ;
406 Get the UNIT of the i^th component.
408 inline const UNIT * FIELD_::getComponentUnit(int i) const
410 return &_componentsUnits[i-1] ;
413 Set FIELD components unit.
415 Duplicate the MEDComponentsUnits string array to put components
417 MEDComponentsUnits size must be equal to number of components.
420 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
422 if (NULL == _MEDComponentsUnits)
423 _MEDComponentsUnits = new string[_numberOfComponents] ;
424 for (int i=0; i<_numberOfComponents; i++)
425 _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
428 Set FIELD i^th component unit.
430 i must be >=1 and <= number of components.
432 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
434 _MEDComponentsUnits[i-1]=MEDComponentUnit ;
437 Get a reference to the string array which contain the components units.
439 This Array size is equal to number of components
441 inline const string * FIELD_::getMEDComponentsUnits() const
443 return _MEDComponentsUnits ;
446 Get the string for unit of the i^th component.
448 inline string FIELD_::getMEDComponentUnit(int i) const
450 return _MEDComponentsUnits[i-1] ;
453 Set the iteration number where FIELD has been calculated.
455 inline void FIELD_::setIterationNumber(int IterationNumber)
457 _iterationNumber=IterationNumber;
460 Get the iteration number where FIELD has been calculated.
462 inline int FIELD_::getIterationNumber() const
464 return _iterationNumber ;
467 Set the time (in second) where FIELD has been calculated.
469 inline void FIELD_::setTime(double Time)
474 Get the time (in second) where FIELD has been calculated.
476 inline double FIELD_::getTime() const
481 Set the order number where FIELD has been calculated.
483 It corresponds to internal iteration during one time step.
485 inline void FIELD_::setOrderNumber(int OrderNumber)
487 _orderNumber=OrderNumber ;
490 Get the order number where FIELD has been calculated.
492 inline int FIELD_::getOrderNumber() const
494 return _orderNumber ;
497 Get a reference to the SUPPORT object associated to FIELD.
499 inline const SUPPORT * FIELD_::getSupport() const
504 Set the reference to the SUPPORT object associated to FIELD.
506 Reference is not duplicate, so it must not be deleted.
508 inline void FIELD_::setSupport(const SUPPORT * support)
513 Get the FIELD med value type (MED_INT32 or MED_REEL64).
515 inline med_type_champ FIELD_::getValueType () const
520 Set the FIELD med value type (MED_INT32 or MED_REEL64).
522 inline void FIELD_::setValueType (const med_type_champ ValueType)
524 _valueType = ValueType ;
527 /////////////////////////
528 // END OF CLASS FIELD_ //
529 /////////////////////////
533 This template class contains informations related with a FIELD :
534 - Values of the field
538 template <class T> class FIELD : public FIELD_
541 // ------- Drivers Management Part
544 //-----------------------//
546 //-----------------------//
549 virtual GENDRIVER * run(const string & fileName, FIELD<T> * ptrFIELD) const = 0 ;
552 //-------------------------------------------------------//
553 template <class T2> class INSTANCE_DE : public INSTANCE
554 //-------------------------------------------------------//
557 GENDRIVER * run(const string & fileName, FIELD<T> * ptrFIELD) const
559 return new T2(fileName,ptrFIELD);
563 // static INSTANCE_DE<MED_FIELD_RDONLY_DRIVER> inst_med_rdonly ;
564 static INSTANCE_DE<MED_FIELD_RDWR_DRIVER<T> > inst_med ;
565 static INSTANCE_DE<VTK_FIELD_DRIVER<T> > inst_vtk ;
567 static const INSTANCE * const instances[] ;
569 // ------ End of Drivers Management Part
571 // array of value of type T
572 MEDARRAY<T> *_value ;
575 void _operation(const FIELD& m,const FIELD& n, const medModeSwitch mode, char* Op);
576 void _operationInitialize(const FIELD& m,const FIELD& n, char* Op);
577 void _add_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
578 void _sub_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
579 void _mul_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
580 void _div_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
585 FIELD(const FIELD &m);
586 FIELD & operator=(const FIELD &m); // A FAIRE
587 FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode=MED_FULL_INTERLACE) throw (MEDEXCEPTION) ; // Ajout NB Constructeur FIELD avec allocation de memoire de tous ses attribut
588 FIELD(const SUPPORT * Support, driverTypes driverType,
589 const string & fileName="", const string & fieldName="",
590 const int iterationNumber = -1, const int orderNumber = -1)
591 throw (MEDEXCEPTION);
594 const FIELD operator+(const FIELD& m) const;
595 const FIELD operator-(const FIELD& m) const;
596 const FIELD operator*(const FIELD& m) const;
597 const FIELD operator/(const FIELD& m) const;
598 const FIELD operator-() const;
599 FIELD& operator+=(const FIELD& m);
600 FIELD& operator-=(const FIELD& m);
601 FIELD& operator*=(const FIELD& m);
602 FIELD& operator/=(const FIELD& m);
603 static FIELD* add(const FIELD& m, const FIELD& n);
604 static FIELD* sub(const FIELD& m, const FIELD& n);
605 static FIELD* mul(const FIELD& m, const FIELD& n);
606 static FIELD* div(const FIELD& m, const FIELD& n);
607 double normMax() const throw (MEDEXCEPTION);
608 double norm2() const throw (MEDEXCEPTION);
609 void applyLin(T a, T b);
610 template <T T_function(T)> void applyFunc();
611 static FIELD* scalarProduct(const FIELD& m, const FIELD& n);
612 double normL2(int component, const FIELD<double> * p_field_volume=NULL) const;
613 double normL2(const FIELD<double> * p_field_volume=NULL) const;
614 double normL1(int component, const FIELD<double> * p_field_volume=NULL) const;
615 double normL1(const FIELD<double> * p_field_volume=NULL) const;
617 friend class MED_FIELD_RDONLY_DRIVER<T>;
618 friend class MED_FIELD_WRONLY_DRIVER<T>;
620 friend class VTK_FIELD_DRIVER<T>;
621 //friend class MED_FIELD_RDWR_DRIVER <T>;
624 void rmDriver(int index=0);
625 int addDriver(driverTypes driverType,
626 const string & fileName="Default File Name.med",
627 const string & driverFieldName="Default Field Name") ;
628 int addDriver(GENDRIVER & driver);
630 void allocValue(const int NumberOfComponents);
631 void allocValue(const int NumberOfComponents, const int LengthValue);
635 inline void read(int index=0);
636 inline void read(const GENDRIVER & genDriver);
637 inline void write(int index=0, const string & driverName = "");
638 inline void write(const GENDRIVER &);
640 inline void writeAppend(int index=0, const string & driverName = "");
641 inline void writeAppend(const GENDRIVER &);
643 inline void setValue(MEDARRAY<T> *Value);
645 inline MEDARRAY<T>* getvalue() const;
646 inline const T* getValue(medModeSwitch Mode) const;
647 inline const T* getValueI(medModeSwitch Mode,int i) const;
648 inline T getValueIJ(int i,int j) const;
650 inline void setValue(medModeSwitch mode, T* value);
651 inline void setValueI(medModeSwitch mode, int i, T* value);
652 inline void setValueIJ(int i, int j, T value);
655 This fonction feeds the FIELD<double> private attributs _value with the
656 volume of each cells belonging to the argument Support. The field has to be
657 initialised via the constructor FIELD<double>(const SUPPORT * , const int )
658 with Support as SUPPORT argument, 1 has the number of components, and Support
659 has be a SUPPORT on 3D cells. This initialisation could be done by the empty
660 constructor followed by a setSupport and setNumberOfComponents call but it has
661 be followed by a setValueType(MED_REEL64) call.
663 void getVolume() const throw (MEDEXCEPTION) ;
665 This fonction feeds the FIELD<double> private attributs _value with the
666 area of each cells (or faces) belonging to the attribut _support. The field
667 has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
668 const int ) with 1 has the number of components, and _support has be a
669 SUPPORT on 2D cells or 3D faces. This initialisation could be done by the
670 empty constructor followed by a setSupport and setNumberOfComponents call but
671 it has be followed by a setValueType(MED_REEL64) call.
673 void getArea() const throw (MEDEXCEPTION) ;
675 This fonction feeds the FIELD<double> private attributs _value with the
676 length of each segments belonging to the attribut _support. The field has
677 to be initialised via the constructor FIELD<double>(const SUPPORT * ,
678 const int ) with 1 has the number of components, and _support has be a
679 SUPPORT on 3D edges or 2D faces. This initialisation could be done by the
680 empty constructor followed by a setSupport and setNumberOfComponents call but
681 it has be followed by a setValueType(MED_REEL64) call.
683 void getLength() const throw (MEDEXCEPTION) ;
685 This fonction feeds the FIELD<double> private attributs _value with the
686 normal vector of each faces belonging to the attribut _support. The field
687 has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
688 const int ) with the space dimension has the number of components, and
689 _support has be a SUPPORT on 3D or 2D faces. This initialisation could be done
690 by the empty constructor followed by a setSupport and setNumberOfComponents
691 call but it has be followed by a setValueType(MED_REEL64) call.
693 void getNormal() const throw (MEDEXCEPTION) ;
695 This fonction feeds the FIELD<double> private attributs _value with the
696 barycenter of each faces or cells or edges belonging to the attribut _support.
697 The field has to be initialised via the constructor
698 FIELD<double>(const SUPPORT * ,const int ) with the space dimension has the
699 number of components, and _support has be a SUPPORT on 3D cells or 2D faces.
700 This initialisation could be done by the empty constructor followed by a
701 setSupport and setNumberOfComponents call but it has be followed by a
702 setValueType(MED_REEL64) call.
704 void getBarycenter() const throw (MEDEXCEPTION) ;
707 // --------------------
708 // Implemented Methods
709 // --------------------
712 Constructor with no parameter, most of the attribut members are set to NULL.
714 template <class T> FIELD<T>::FIELD():
715 _value((MEDARRAY<T>*)NULL)
717 MESSAGE("Constructeur FIELD sans parametre");
721 Constructor with parameters such that all attrribut are set but the _value
722 attrribut is allocated but not set.
724 template <class T> FIELD<T>::FIELD(const SUPPORT * Support,
725 const int NumberOfComponents, const medModeSwitch Mode) throw (MEDEXCEPTION) :
726 FIELD_(Support, NumberOfComponents)
728 BEGIN_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
732 _numberOfValues = Support->getNumberOfElements(MED_ALL_ELEMENTS);
734 catch (MEDEXCEPTION &ex) {
735 MESSAGE("No value defined ! ("<<ex.what()<<")");
736 _value = (MEDARRAY<T>*)NULL ;
738 MESSAGE("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
739 if (0<_numberOfValues) {
740 _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues,Mode);
743 _value = (MEDARRAY<T>*)NULL ;
745 END_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
752 template <class T> void FIELD<T>::init ()
759 template <class T> FIELD<T>::FIELD(const FIELD & m):
762 MESSAGE("Constructeur FIELD de recopie");
763 if (m._value != NULL)
765 // copie only default !
766 _value = new MEDARRAY<T>(* m._value,false);
769 _value = (MEDARRAY<T> *) NULL;
770 //_drivers = m._drivers;
776 template <class T> FIELD<T> & FIELD<T>::operator=(const FIELD &m)
778 MESSAGE("Appel de FIELD<T>::operator=");
782 Overload addition operator.
783 This operation is authorized only for compatible fields that have the same support.
784 The compatibility checking includes equality tests of the folowing data members:/n
786 - _numberOfComponents
789 - _MEDComponentsUnits.
791 The data members of the returned field are initialized, based on the first field, except for the name,
792 which is the combination of the two field's names and the operator.
793 Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> /n
794 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
795 When using python, this operator calls the copy constructor in any case.
796 The user has to be aware that when using operator + in associatives expressions like
797 <tt> a = b + c + d +e; </tt> /n
798 no optimisation is performed : the evaluation of last expression requires the construction of
802 const FIELD<T> FIELD<T>::operator+(const FIELD & m) const
804 BEGIN_OF("FIELD<T>::operator+(const FIELD & m)");
805 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
807 // Select mode : avoid if possible any calculation of other mode for fields m or *this
809 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
810 mode=m.getvalue()->getMode();
812 mode=this->getvalue()->getMode();
814 // Creation of the result - memory is allocated by FIELD constructor
815 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
816 //result._operation(*this,m,mode,"+"); // perform Atribute's initialization & addition
817 result._operationInitialize(*this,m,"+"); // perform Atribute's initialization
818 result._add_in_place(*this,m,mode); // perform addition
820 END_OF("FIELD<T>::operator+(const FIELD & m)");
824 /*! Overloaded Operator +=
825 * Operations are directly performed in the first field's array.
826 * This operation is authorized only for compatible fields that have the same support.
829 FIELD<T>& FIELD<T>::operator+=(const FIELD & m)
831 BEGIN_OF("FIELD<T>::operator+=(const FIELD & m)");
832 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
834 // We choose to keep *this mode, even if it may cost a re-calculation for m
835 medModeSwitch mode=this->getvalue()->getMode();
836 const T* value1=m.getValue(mode); // get pointers to the values we are adding
838 // get a non const pointer to the inside array of values and perform operation
839 T * value=const_cast<T *> (getValue(mode));
840 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
841 const T* endV=value+size; // pointer to the end of value
842 for(;value!=endV; value1++,value++)
844 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
845 this->getvalue()->clearOtherMode();
846 END_OF("FIELD<T>::operator+=(const FIELD & m)");
851 /*! Addition of fields. Static member function.
852 * The function return a pointer to a new created field that holds the addition.
853 * Data members are checked for compatibility and initialized.
854 * The user is in charge of memory deallocation.
857 FIELD<T>* FIELD<T>::add(const FIELD& m, const FIELD& n)
859 BEGIN_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
860 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
862 // Select mode : avoid if possible any calculation of other mode for fields m or *this
864 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
865 mode=m.getvalue()->getMode();
867 mode=n.getvalue()->getMode();
869 // Creation of a new field
870 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
871 result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
872 result->_add_in_place(m,n,mode); // perform addition
874 END_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
879 Overload substraction operator.
880 This operation is authorized only for compatible fields that have the same support.
881 The compatibility checking includes equality tests of the folowing data members:/n
883 - _numberOfComponents
886 - _MEDComponentsUnits.
888 The data members of the returned field are initialized, based on the first field, except for the name,
889 which is the combination of the two field's names and the operator.
890 Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> /n
891 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
892 When using python, this operator calls the copy constructor in any case.
893 The user has to be aware that when using operator - in associatives expressions like
894 <tt> a = b - c - d -e; </tt> /n
895 no optimisation is performed : the evaluation of last expression requires the construction of
899 const FIELD<T> FIELD<T>::operator-(const FIELD & m) const
901 BEGIN_OF("FIELD<T>::operator-(const FIELD & m)");
902 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
904 // Select mode : avoid if possible any calculation of other mode for fields m or *this
906 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
907 mode=m.getvalue()->getMode();
909 mode=this->getvalue()->getMode();
911 // Creation of the result - memory is allocated by FIELD constructor
912 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
913 //result._operation(*this,m,mode,"-"); // perform Atribute's initialization & substraction
914 result._operationInitialize(*this,m,"-"); // perform Atribute's initialization
915 result._sub_in_place(*this,m,mode); // perform substracion
917 END_OF("FIELD<T>::operator-(const FIELD & m)");
922 const FIELD<T> FIELD<T>::operator-() const
924 BEGIN_OF("FIELD<T>::operator-()");
926 medModeSwitch mode=this->getvalue()->getMode();
927 // Creation of the result - memory is allocated by FIELD constructor
928 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
929 // Atribute's initialization
930 result.setName("- "+getName());
931 result.setComponentsNames(getComponentsNames());
932 // not yet implemented setComponentType(getComponentType());
933 result.setComponentsDescriptions(getComponentsDescriptions());
934 result.setMEDComponentsUnits(getMEDComponentsUnits());
935 result.setComponentsUnits(getComponentsUnits());
936 result.setIterationNumber(getIterationNumber());
937 result.setTime(getTime());
938 result.setOrderNumber(getOrderNumber());
939 result.setValueType(getValueType());
941 const T* value1=getValue(mode);
942 // get a non const pointer to the inside array of values and perform operation
943 T * value=const_cast<T *> (result.getValue(mode));
944 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
945 const T* endV=value+size; // pointer to the end of value
947 for(;value!=endV; value1++,value++)
949 END_OF("FIELD<T>::operator-=(const FIELD & m)");
953 /*! Overloaded Operator -=
954 * Operations are directly performed in the first field's array.
955 * This operation is authorized only for compatible fields that have the same support.
958 FIELD<T>& FIELD<T>::operator-=(const FIELD & m)
960 BEGIN_OF("FIELD<T>::operator-=(const FIELD & m)");
961 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
963 // We choose to keep *this mode, even if it may cost a re-calculation for m
964 medModeSwitch mode=this->getvalue()->getMode();
965 const T* value1=m.getValue(mode); // get pointers to the values we are adding
967 // get a non const pointer to the inside array of values and perform operation
968 T * value=const_cast<T *> (getValue(mode));
969 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
970 const T* endV=value+size; // pointer to the end of value
972 for(;value!=endV; value1++,value++)
974 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
975 this->getvalue()->clearOtherMode();
976 END_OF("FIELD<T>::operator-=(const FIELD & m)");
981 /*! Substraction of fields. Static member function.
982 * The function return a pointer to a new created field that holds the substraction.
983 * Data members are checked for compatibility and initialized.
984 * The user is in charge of memory deallocation.
987 FIELD<T>* FIELD<T>::sub(const FIELD& m, const FIELD& n)
989 BEGIN_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
990 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
992 // Select mode : avoid if possible any calculation of other mode for fields m or *this
994 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
995 mode=m.getvalue()->getMode();
997 mode=n.getvalue()->getMode();
999 // Creation of a new field
1000 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1001 result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1002 result->_sub_in_place(m,n,mode); // perform substraction
1004 END_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
1009 Overload multiplication operator.
1010 This operation is authorized only for compatible fields that have the same support.
1011 The compatibility checking includes equality tests of the folowing data members:/n
1013 - _numberOfComponents
1016 - _MEDComponentsUnits.
1018 The data members of the returned field are initialized, based on the first field, except for the name,
1019 which is the combination of the two field's names and the operator.
1020 Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> /n
1021 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1022 When using python, this operator calls the copy constructor in any case.
1023 The user has to be aware that when using operator * in associatives expressions like
1024 <tt> a = b * c * d *e; </tt> /n
1025 no optimisation is performed : the evaluation of last expression requires the construction of
1029 const FIELD<T> FIELD<T>::operator*(const FIELD & m) const
1031 BEGIN_OF("FIELD<T>::operator*(const FIELD & m)");
1032 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1034 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1036 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
1037 mode=m.getvalue()->getMode();
1039 mode=this->getvalue()->getMode();
1041 // Creation of the result - memory is allocated by FIELD constructor
1042 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
1043 //result._operation(*this,m,mode,"*"); // perform Atribute's initialization & multiplication
1044 result._operationInitialize(*this,m,"*"); // perform Atribute's initialization
1045 result._mul_in_place(*this,m,mode); // perform multiplication
1047 END_OF("FIELD<T>::operator*(const FIELD & m)");
1051 /*! Overloaded Operator *=
1052 * Operations are directly performed in the first field's array.
1053 * This operation is authorized only for compatible fields that have the same support.
1056 FIELD<T>& FIELD<T>::operator*=(const FIELD & m)
1058 BEGIN_OF("FIELD<T>::operator*=(const FIELD & m)");
1059 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1061 // We choose to keep *this mode, even if it may cost a re-calculation for m
1062 medModeSwitch mode=this->getvalue()->getMode();
1063 const T* value1=m.getValue(mode); // get pointers to the values we are adding
1065 // get a non const pointer to the inside array of values and perform operation
1066 T * value=const_cast<T *> (getValue(mode));
1067 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1068 const T* endV=value+size; // pointer to the end of value
1070 for(;value!=endV; value1++,value++)
1072 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1073 this->getvalue()->clearOtherMode();
1074 END_OF("FIELD<T>::operator*=(const FIELD & m)");
1079 /*! Multiplication of fields. Static member function.
1080 * The function return a pointer to a new created field that holds the multiplication.
1081 * Data members are checked for compatibility and initialized.
1082 * The user is in charge of memory deallocation.
1085 FIELD<T>* FIELD<T>::mul(const FIELD& m, const FIELD& n)
1087 BEGIN_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1088 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1090 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1092 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1093 mode=m.getvalue()->getMode();
1095 mode=n.getvalue()->getMode();
1097 // Creation of a new field
1098 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1099 result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1100 result->_mul_in_place(m,n,mode); // perform multiplication
1102 END_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1108 Overload division operator.
1109 This operation is authorized only for compatible fields that have the same support.
1110 The compatibility checking includes equality tests of the folowing data members:/n
1112 - _numberOfComponents
1115 - _MEDComponentsUnits.
1117 The data members of the returned field are initialized, based on the first field, except for the name,
1118 which is the combination of the two field's names and the operator.
1119 Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> /n
1120 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1121 When using python, this operator calls the copy constructor in any case.
1122 The user has to be aware that when using operator / in associatives expressions like
1123 <tt> a = b / c / d /e; </tt> /n
1124 no optimisation is performed : the evaluation of last expression requires the construction of
1128 const FIELD<T> FIELD<T>::operator/(const FIELD & m) const
1130 BEGIN_OF("FIELD<T>::operator/(const FIELD & m)");
1131 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1133 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1135 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
1136 mode=m.getvalue()->getMode();
1138 mode=this->getvalue()->getMode();
1140 // Creation of the result - memory is allocated by FIELD constructor
1141 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
1142 //result._operation(*this,m,mode,"/"); // perform Atribute's initialization & division
1143 result._operationInitialize(*this,m,"/"); // perform Atribute's initialization
1144 result._div_in_place(*this,m,mode); // perform division
1146 END_OF("FIELD<T>::operator/(const FIELD & m)");
1151 /*! Overloaded Operator /=
1152 * Operations are directly performed in the first field's array.
1153 * This operation is authorized only for compatible fields that have the same support.
1156 FIELD<T>& FIELD<T>::operator/=(const FIELD & m)
1158 BEGIN_OF("FIELD<T>::operator/=(const FIELD & m)");
1159 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1161 // We choose to keep *this mode, even if it may cost a re-calculation for m
1162 medModeSwitch mode=this->getvalue()->getMode();
1163 const T* value1=m.getValue(mode); // get pointers to the values we are adding
1165 // get a non const pointer to the inside array of values and perform operation
1166 T * value=const_cast<T *> (getValue(mode));
1167 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1168 const T* endV=value+size; // pointer to the end of value
1170 for(;value!=endV; value1++,value++)
1172 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1173 this->getvalue()->clearOtherMode();
1174 END_OF("FIELD<T>::operator/=(const FIELD & m)");
1179 /*! Division of fields. Static member function.
1180 * The function return a pointer to a new created field that holds the division.
1181 * Data members are checked for compatibility and initialized.
1182 * The user is in charge of memory deallocation.
1185 FIELD<T>* FIELD<T>::div(const FIELD& m, const FIELD& n)
1187 BEGIN_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1188 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1190 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1192 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1193 mode=m.getvalue()->getMode();
1195 mode=n.getvalue()->getMode();
1197 // Creation of a new field
1198 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1199 result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1200 result->_div_in_place(m,n,mode); // perform division
1202 END_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1209 This internal method initialize the members of a new field created to hold the result of the operation Op .
1210 Initialization is based on the first field, except for the name, which is the combination of the two field's names
1215 void FIELD<T>::_operationInitialize(const FIELD& m,const FIELD& n, char* Op)
1217 MESSAGE("Appel methode interne _add" << Op);
1219 // Atribute's initialization (copy of the first field's attributes)
1220 // Other data members (_support, _numberOfValues) are initialized in the field's constr.
1221 setName(m.getName()+" "+Op+" "+n.getName());
1222 setComponentsNames(m.getComponentsNames());
1223 // not yet implemented setComponentType(m.getComponentType());
1224 setComponentsDescriptions(m.getComponentsDescriptions());
1225 setMEDComponentsUnits(m.getMEDComponentsUnits());
1227 // The following data member may differ from field m to n.
1228 // The initialization is done based on the first field.
1229 setComponentsUnits(m.getComponentsUnits());
1230 setIterationNumber(m.getIterationNumber());
1231 setTime(m.getTime());
1232 setOrderNumber(m.getOrderNumber());
1233 setValueType(m.getValueType());
1239 Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
1240 This method is applied to a just created field with medModeSwitch mode.
1241 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1246 void FIELD<T>::_add_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1248 // get pointers to the values we are adding
1249 const T* value1=m.getValue(mode);
1250 const T* value2=n.getValue(mode);
1251 // get a non const pointer to the inside array of values and perform operation
1252 T * value=const_cast<T *> (getValue(mode));
1254 const int size=getNumberOfValues()*getNumberOfComponents();
1256 const T* endV1=value1+size;
1257 for(;value1!=endV1; value1++,value2++,value++)
1258 *value=(*value1)+(*value2);
1263 Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
1264 This method is applied to a just created field with medModeSwitch mode.
1265 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1270 void FIELD<T>::_sub_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1272 // get pointers to the values we are adding
1273 const T* value1=m.getValue(mode);
1274 const T* value2=n.getValue(mode);
1275 // get a non const pointer to the inside array of values and perform operation
1276 T * value=const_cast<T *> (getValue(mode));
1278 const int size=getNumberOfValues()*getNumberOfComponents();
1280 const T* endV1=value1+size;
1281 for(;value1!=endV1; value1++,value2++,value++)
1282 *value=(*value1)-(*value2);
1287 Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
1288 This method is applied to a just created field with medModeSwitch mode.
1289 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1294 void FIELD<T>::_mul_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1296 // get pointers to the values we are adding
1297 const T* value1=m.getValue(mode);
1298 const T* value2=n.getValue(mode);
1299 // get a non const pointer to the inside array of values and perform operation
1300 T * value=const_cast<T *> (getValue(mode));
1302 const int size=getNumberOfValues()*getNumberOfComponents();
1304 const T* endV1=value1+size;
1305 for(;value1!=endV1; value1++,value2++,value++)
1306 *value=(*value1)*(*value2);
1311 Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
1312 This method is applied to a just created field with medModeSwitch mode.
1313 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1318 void FIELD<T>::_div_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1320 // get pointers to the values we are adding
1321 const T* value1=m.getValue(mode);
1322 const T* value2=n.getValue(mode);
1323 // get a non const pointer to the inside array of values and perform operation
1324 T * value=const_cast<T *> (getValue(mode));
1326 const int size=getNumberOfValues()*getNumberOfComponents();
1328 const T* endV1=value1+size;
1329 for(;value1!=endV1; value1++,value2++,value++)
1330 *value=(*value1)/(*value2);
1335 template <class T> double FIELD<T>::normMax() const throw (MEDEXCEPTION)
1337 const T* value=getValue(getvalue()->getMode()); // get pointer to the values
1338 const int size=getNumberOfValues()*getNumberOfComponents();
1339 if (size <= 0) // Size of array has to be strictly positive
1342 diagnosis="FIELD<T>::normMax() : cannot compute the norm of "+getName()+
1343 " : it size is non positive!";
1344 throw MEDEXCEPTION(diagnosis.c_str());
1346 const T* lastvalue=value+size; // get pointer just after last value
1347 const T* pMax=value; // pointer to the max value
1348 const T* pMin=value; // pointer to the min value
1350 // get pointers to the max & min value of array
1351 while ( ++value != lastvalue )
1353 if ( *pMin > *value )
1355 if ( *pMax < *value )
1359 T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
1360 T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
1362 return Max>Min ? static_cast<double>(Max) : static_cast<double>(Min);
1365 /*! Return Euclidien norm
1367 template <class T> double FIELD<T>::norm2() const throw (MEDEXCEPTION)
1369 const T* value=this->getValue(this->getvalue()->getMode()); // get const pointer to the values
1370 const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
1371 if (size <= 0) // Size of array has to be strictly positive
1374 diagnosis="FIELD<T>::norm2() : cannot compute the norm of "+getName()+
1375 " : it size is non positive!";
1376 throw MEDEXCEPTION(diagnosis.c_str());
1378 const T* lastvalue=value+size; // point just after last value
1380 T result((T)0); // init
1381 for( ; value!=lastvalue ; ++value)
1382 result += (*value) * (*value);
1384 return std::sqrt(static_cast<double> (result));
1388 /*! Apply to each (scalar) field component the template parameter T_function,
1389 * which is a pointer to function.
1390 * Since the pointer is known at compile time, the function is inlined into the inner loop!
1391 * calculation is done "in place".
1394 * \code myField.applyFunc<std::sqrt>(); // apply sqare root function \endcode
1395 * \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
1397 template <class T> template <T T_function(T)>
1398 void FIELD<T>::applyFunc()
1400 medModeSwitch mode=getvalue()->getMode();
1402 // get a non const pointer to the inside array of values and perform operation
1403 T * value=const_cast<T *> (getValue(mode));
1404 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1406 if (size>0) // for a negative size, there is nothing to do
1408 const T* lastvalue=value+size; // pointer to the end of value
1409 for(;value!=lastvalue; ++value) // apply linear transformation
1410 *value = T_function(*value);
1411 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1412 getvalue()->clearOtherMode();
1417 /*! Apply to each (scalar) field component the linear function x -> ax+b.
1418 * calculation is done "in place".
1420 template <class T> void FIELD<T>::applyLin(T a, T b)
1422 medModeSwitch mode=getvalue()->getMode();
1424 // get a non const pointer to the inside array of values and perform operation in place
1425 T * value=const_cast<T *> (getValue(mode));
1426 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1428 if (size>0) // for a negative size, there is nothing to do
1430 const T* lastvalue=value+size; // pointer to the end of value
1431 for(;value!=lastvalue; ++value) // apply linear transformation
1432 *value = a*(*value)+b;
1433 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1434 getvalue()->clearOtherMode();
1440 * Return a pointer to a new field that holds the scalar product. Static member function.
1441 * This operation is authorized only for compatible fields that have the same support.
1442 * The compatibility checking includes equality tests of the folowing data members:/n
1444 * - _numberOfComponents
1446 * - _componentsTypes
1447 * - _MEDComponentsUnits.
1448 * Data members are initialized.
1449 * The new field point to the same support and has one component.
1450 * Each value of it is the scalar product of the two argument's fields.
1451 * The user is in charge of memory deallocation.
1453 template <class T> FIELD<T>* FIELD<T>::scalarProduct(const FIELD & m, const FIELD & n)
1455 FIELD_::_checkFieldCompatibility( m, n); // may throw exception
1456 // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
1457 const medModeSwitch mode=MED_FULL_INTERLACE;
1459 const int numberOfElements=m.getNumberOfValues(); // strictly positive
1460 const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
1462 // Creation & init of a the result field on the same suppot, with one component
1463 FIELD<T>* result = new FIELD<T>(m.getSupport(),1,mode);
1464 result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
1465 result->setIterationNumber(m.getIterationNumber());
1466 result->setTime(m.getTime());
1467 result->setOrderNumber(m.getOrderNumber());
1468 result->setValueType(m.getValueType());
1470 const T* value1=m.getValue(mode); // get const pointer to the values
1471 const T* value2=n.getValue(mode); // get const pointer to the values
1472 // get a non const pointer to the inside array of values and perform operation
1473 T * value=const_cast<T *> (result->getValue(mode));
1475 const T* lastvalue=value+numberOfElements; // pointing just after last value of result
1476 for ( ; value!=lastvalue ; ++value ) // loop on all elements
1478 *value=(T)0; // initialize value
1479 const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
1480 for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
1481 *value += (*value1) * (*value2);
1486 /*! Return L2 Norm of the field's component.
1487 * Cannot be applied to a field with a support on nodes.
1488 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1490 template <class T> double FIELD<T>::normL2(int component, const FIELD<double> * p_field_volume) const
1492 _checkNormCompatibility(p_field_volume); // may throw exception
1493 if ( component<1 || component>getNumberOfComponents() )
1494 throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
1496 const FIELD<double> * p_field_size=p_field_volume;
1497 if(!p_field_volume) // if the user don't supply the volume
1498 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1500 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1501 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1502 const T* value=getValueI( MED_NO_INTERLACE, component); // get pointer to the component's values
1503 const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
1505 double integrale=0.0;
1507 for (; value!=lastvalue ; ++value ,++vol)
1509 integrale += static_cast<double>((*value) * (*value)) * (*vol);
1513 if(!p_field_volume) // if the user didn't supply the volume
1514 delete p_field_size; // delete temporary volume field
1516 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1518 return integrale/totVol;
1521 /*! Return L2 Norm of the field.
1522 * Cannot be applied to a field with a support on nodes.
1523 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1525 template <class T> double FIELD<T>::normL2(const FIELD<double> * p_field_volume) const
1527 _checkNormCompatibility(p_field_volume); // may throw exception
1528 const FIELD<double> * p_field_size=p_field_volume;
1529 if(!p_field_volume) // if the user don't supply the volume
1530 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1532 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1533 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1534 const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
1535 const T* value=getValue( MED_NO_INTERLACE); // get pointer to the field's values
1538 const double* p_vol=vol;
1539 for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
1542 double integrale=0.0;
1543 for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
1544 for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
1545 integrale += static_cast<double>((*value) * (*value)) * (*p_vol);
1547 if(!p_field_volume) // if the user didn't supply the volume
1548 delete p_field_size; // delete temporary volume field
1550 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1552 return integrale/totVol;
1555 /*! Return L1 Norm of the field's component.
1556 * Cannot be applied to a field with a support on nodes.
1557 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1559 template <class T> double FIELD<T>::normL1(int component, const FIELD<double> * p_field_volume) const
1561 _checkNormCompatibility(p_field_volume); // may throw exception
1562 if ( component<1 || component>getNumberOfComponents() )
1563 throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
1565 const FIELD<double> * p_field_size=p_field_volume;
1566 if(!p_field_volume) // if the user don't supply the volume
1567 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1569 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1570 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1571 const T* value=getValueI( MED_NO_INTERLACE, component); // get pointer to the component's values
1572 const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
1574 double integrale=0.0;
1576 for (; value!=lastvalue ; ++value ,++vol)
1578 integrale += std::abs( static_cast<double>(*value) ) * (*vol);
1582 if(!p_field_volume) // if the user didn't supply the volume
1583 delete p_field_size; // delete temporary volume field
1585 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1587 return integrale/totVol;
1590 /*! Return L1 Norm of the field.
1591 * Cannot be applied to a field with a support on nodes.
1592 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1594 template <class T> double FIELD<T>::normL1(const FIELD<double> * p_field_volume) const
1596 _checkNormCompatibility(p_field_volume); // may throw exception
1597 const FIELD<double> * p_field_size=p_field_volume;
1598 if(!p_field_volume) // if the user don't supply the volume
1599 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1601 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1602 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1603 const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
1604 const T* value=getValue( MED_NO_INTERLACE); // get pointer to the field's values
1607 const double* p_vol=vol;
1608 for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
1611 double integrale=0.0;
1612 for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
1613 for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
1614 integrale += std::abs( static_cast<double>(*value) ) * (*p_vol);
1616 if(!p_field_volume) // if the user didn't supply the volume
1617 delete p_field_size; // delete temporary volume field
1619 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1621 return integrale/totVol;
1628 Constructor with parameters; the object is set via a file and its associated
1629 driver. For the moment only the MED_DRIVER is considered and if the last two
1630 argument (iterationNumber and orderNumber) are not set; their default value
1631 is -1. If the field fieldDriverName with the iteration number
1632 iterationNumber and the order number orderNumber does not exist in the file
1633 fieldDriverName; the constructor raises an exception.
1635 template <class T> FIELD<T>::FIELD(const SUPPORT * Support,
1636 driverTypes driverType,
1637 const string & fileName/*=""*/,
1638 const string & fieldDriverName/*=""*/,
1639 const int iterationNumber,
1640 const int orderNumber)
1641 throw (MEDEXCEPTION)
1643 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) : ";
1652 _value = (MEDARRAY<T>*)NULL;
1654 _iterationNumber = iterationNumber;
1656 _orderNumber = orderNumber;
1662 MED_FIELD_RDONLY_DRIVER<T> myDriver(fileName,this);
1663 myDriver.setFieldName(fieldDriverName);
1664 current = addDriver(myDriver);
1667 // current = addDriver(driverType,fileName,fieldDriverName);
1668 // switch(_drivers[current]->getAccessMode() ) {
1669 // case MED_WRONLY : {
1670 // MESSAGE("FIELD<T>::FIELD(driverTypes driverType, .....) : driverType must have a MED_RDONLY or MED_RDWR accessMode");
1671 // rmDriver(current);
1678 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Driver type unknown : can't create driver!"));
1683 _drivers[current]->open();
1684 _drivers[current]->read();
1685 _drivers[current]->close();
1693 template <class T> FIELD<T>::~FIELD()
1695 BEGIN_OF(" Destructeur FIELD<T>::~FIELD()");
1697 if (_value) delete _value;
1698 END_OF(" Destructeur FIELD<T>::~FIELD()");
1704 template <class T> void FIELD<T>::allocValue(const int NumberOfComponents)
1706 const char* LOC = "FIELD<T>::allocValue(const int NumberOfComponents)" ;
1709 _numberOfComponents = NumberOfComponents ;
1710 if (_componentsTypes == NULL)
1711 _componentsTypes = new int[NumberOfComponents] ;
1712 if (_componentsNames == NULL)
1713 _componentsNames = new string[NumberOfComponents];
1714 if (_componentsDescriptions == NULL)
1715 _componentsDescriptions = new string[NumberOfComponents];
1716 if (_componentsUnits == NULL)
1717 _componentsUnits = new UNIT[NumberOfComponents];
1718 if (_MEDComponentsUnits == NULL)
1719 _MEDComponentsUnits = new string[NumberOfComponents];
1720 for (int i=0;i<NumberOfComponents;i++) {
1721 _componentsTypes[i] = 0 ;
1725 _numberOfValues = _support->getNumberOfElements(MED_ALL_ELEMENTS);
1726 MESSAGE(LOC <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
1728 _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
1732 catch (MEDEXCEPTION &ex) {
1733 MESSAGE("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
1734 _value = (MEDARRAY<T>*)NULL ;
1738 END_OF("void FIELD<T>::allocValue(const int NumberOfComponents)");
1744 template <class T> void FIELD<T>::allocValue(const int NumberOfComponents, const int LengthValue)
1746 BEGIN_OF("void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)");
1748 _numberOfComponents = NumberOfComponents ;
1749 if (_componentsTypes == NULL)
1750 _componentsTypes = new int[NumberOfComponents] ;
1751 if (_componentsNames == NULL)
1752 _componentsNames = new string[NumberOfComponents];
1753 if (_componentsDescriptions == NULL)
1754 _componentsDescriptions = new string[NumberOfComponents];
1755 if (_componentsUnits == NULL)
1756 _componentsUnits = new UNIT[NumberOfComponents];
1757 if (_MEDComponentsUnits == NULL)
1758 _MEDComponentsUnits = new string[NumberOfComponents];
1759 for (int i=0;i<NumberOfComponents;i++) {
1760 _componentsTypes[i] = 0 ;
1763 MESSAGE("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
1764 _numberOfValues = LengthValue ;
1765 _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
1769 END_OF("void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)");
1775 template <class T> void FIELD<T>::deallocValue()
1777 BEGIN_OF("void FIELD<T>::deallocValue()");
1778 _numberOfValues = 0 ;
1779 _numberOfComponents = 0 ;
1783 END_OF("void FIELD<T>::deallocValue()");
1786 // -----------------
1788 // -----------------
1791 template <class T> FIELD<T>::INSTANCE_DE<MED_FIELD_RDWR_DRIVER<T> > FIELD<T>::inst_med ;
1793 template <class T> FIELD<T>::INSTANCE_DE<VTK_FIELD_DRIVER<T> > FIELD<T>::inst_vtk ;
1795 template <class T> const typename FIELD<T>::INSTANCE * const FIELD<T>::instances[] = { &FIELD<T>::inst_med, &FIELD<T>::inst_vtk} ;
1799 Create the specified driver and return its index reference to path to
1800 read or write methods.
1803 template <class T> int FIELD<T>::addDriver(driverTypes driverType,
1804 const string & fileName/*="Default File Name.med"*/,
1805 const string & driverName/*="Default Field Name"*/)
1807 const char * LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Field Name\") : ";
1811 int itDriver = (int) NO_DRIVER;
1817 SCRUTE(instances[driverType]);
1822 itDriver = (int) driverType ;
1831 case GIBI_DRIVER : {
1832 throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "driverType other than MED_DRIVER and VTK_DRIVER has been specified to the method which is not allowed for the object FIELD"));
1837 throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "driverType other than MED_DRIVER and VTK_DRIVER has been specified to the method which is not allowed for the object FIELD"));
1842 if (itDriver == ((int) NO_DRIVER))
1843 throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "driverType other than MED_DRIVER and VTK_DRIVER has been specified to the method which is not allowed for the object FIELD"));
1845 driver = instances[itDriver]->run(fileName, this) ;
1847 _drivers.push_back(driver);
1849 current = _drivers.size()-1;
1851 _drivers[current]->setFieldName(driverName);
1860 Duplicate the given driver and return its index reference to path to
1861 read or write methods.
1863 template <class T> inline int FIELD<T>::addDriver (GENDRIVER & driver )
1865 const char * LOC = "FIELD<T>::addDriver(GENDRIVER &) : ";
1868 // duplicate driver to delete it with destructor !
1869 GENDRIVER * newDriver = driver.copy() ;
1871 _drivers.push_back(newDriver);
1872 return _drivers.size() -1 ;
1878 Remove the driver referenced by its index.
1880 template <class T> void FIELD<T>::rmDriver (int index/*=0*/)
1882 const char * LOC = "FIELD<T>::rmDriver (int index=0): ";
1885 if ( _drivers[index] ) {
1886 //_drivers.erase(&_drivers[index]);
1888 MESSAGE ("detruire");
1891 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1892 << "The <index given is invalid, index must be between 0 and |"
1901 Read FIELD in the file specified in the driver given by its index.
1903 template <class T> inline void FIELD<T>::read(int index/*=0*/)
1905 const char * LOC = "FIELD<T>::read(int index=0) : ";
1908 if ( _drivers[index] ) {
1909 _drivers[index]->open();
1910 _drivers[index]->read();
1911 _drivers[index]->close();
1914 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1915 << "The index given is invalid, index must be between 0 and |"
1923 Write FIELD in the file specified in the driver given by its index.
1925 template <class T> inline void FIELD<T>::write(int index/*=0*/, const string & driverName /*= ""*/)
1927 const char * LOC = "FIELD<T>::write(int index=0, const string & driverName = \"\") : ";
1930 if( _drivers[index] ) {
1931 _drivers[index]->open();
1932 if (driverName != "") _drivers[index]->setFieldName(driverName);
1933 _drivers[index]->write();
1934 _drivers[index]->close();
1937 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1938 << "The index given is invalid, index must be between 0 and |"
1946 Write FIELD in the file specified in the driver given by its index. Use this
1947 method for ASCII drivers (e.g. VTK_DRIVER)
1949 template <class T> inline void FIELD<T>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
1951 const char * LOC = "FIELD<T>::write(int index=0, const string & driverName = \"\") : ";
1954 if( _drivers[index] ) {
1955 _drivers[index]->openAppend();
1956 if (driverName != "") _drivers[index]->setFieldName(driverName);
1957 _drivers[index]->writeAppend();
1958 _drivers[index]->close();
1961 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1962 << "The index given is invalid, index must be between 0 and |"
1971 Write FIELD with the driver which is equal to the given driver.
1975 template <class T> inline void FIELD<T>::write(const GENDRIVER & genDriver)
1977 const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
1980 for (unsigned int index=0; index < _drivers.size(); index++ )
1981 if ( *_drivers[index] == genDriver ) {
1982 _drivers[index]->open();
1983 _drivers[index]->write();
1984 _drivers[index]->close();
1993 Write FIELD with the driver which is equal to the given driver.
1995 Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
1997 template <class T> inline void FIELD<T>::writeAppend(const GENDRIVER & genDriver)
1999 const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
2002 for (unsigned int index=0; index < _drivers.size(); index++ )
2003 if ( *_drivers[index] == genDriver ) {
2004 _drivers[index]->openAppend();
2005 _drivers[index]->writeAppend();
2006 _drivers[index]->close();
2015 Read FIELD with the driver which is equal to the given driver.
2019 template <class T> inline void FIELD<T>::read(const GENDRIVER & genDriver)
2021 const char * LOC = " FIELD<T>::read(const GENDRIVER &) : ";
2024 for (unsigned int index=0; index < _drivers.size(); index++ )
2025 if ( *_drivers[index] == genDriver ) {
2026 _drivers[index]->open();
2027 _drivers[index]->read();
2028 _drivers[index]->close();
2037 Destroy the MEDARRAY<T> in FIELD and put the new one without copy.
2040 template <class T> inline void FIELD<T>::setValue(MEDARRAY<T> *Value)
2042 if (NULL != _value) delete _value ;
2048 Return a reference to the MEDARRAY<T> in FIELD.
2051 template <class T> inline MEDARRAY<T>* FIELD<T>::getvalue() const
2057 Return a reference to values array to read them.
2059 template <class T> inline const T* FIELD<T>::getValue(medModeSwitch Mode) const
2061 return _value->get(Mode) ;
2065 Return a reference to i^{th} row or column - component - (depend on Mode value)
2066 of FIELD values array.
2068 template <class T> inline const T* FIELD<T>::getValueI(medModeSwitch Mode,int i) const
2070 if ( Mode == MED_FULL_INTERLACE )
2072 return _value->getRow(i) ;
2074 ASSERT ( Mode == MED_NO_INTERLACE);
2075 return _value->getColumn(i);
2079 Return the value of i^{th} element and j^{th} component.
2081 template <class T> inline T FIELD<T>::getValueIJ(int i,int j) const
2083 return _value->getIJ(i,j) ;
2087 Copy new values array in FIELD according to the given mode.
2089 Array must have right size. If not results are unpredicable.
2091 template <class T> inline void FIELD<T>::setValue(medModeSwitch mode, T* value)
2093 _value->set(mode,value);
2097 Update values array in FIELD with the given ones according to specified mode.
2099 template <class T> inline void FIELD<T>::setValueI(medModeSwitch mode, int i, T* value)
2102 if (MED_FULL_INTERLACE == mode)
2103 _value->setI(i,value);
2104 else if (MED_NO_INTERLACE == mode)
2105 _value->setJ(i,value);
2107 throw MEDEXCEPTION(LOCALIZED("FIELD<T>::setValueI : bad medModeSwitch")) ;
2111 Set the value of i^{th} element and j^{th} component with the given one.
2113 template <class T> inline void FIELD<T>::setValueIJ(int i, int j, T value)
2115 _value->setIJ(i,j,value);
2123 Fill values array with volume values.
2125 template <class T> void FIELD<T>::getVolume() const throw (MEDEXCEPTION)
2127 const char * LOC = "FIELD<double>::getVolume() const : ";
2130 // The field has to be initilised by a non empty support and a
2131 // number of components = 1 and its value type has to be set to MED_REEL64
2132 // (ie a FIELD<double>)
2134 if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_REEL64))
2135 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to 1 and a value type set to MED_REEL64"));
2141 Fill values array with area values.
2143 template <class T> void FIELD<T>::getArea() const throw (MEDEXCEPTION)
2145 const char * LOC = "FIELD<double>::getArea() const : ";
2148 // The field has to be initilised by a non empty support and a
2149 // number of components = 1 and its value type has to be set to MED_REEL64
2150 // (ie a FIELD<double>)
2152 if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_REEL64))
2153 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to 1 and a value type set to MED_REEL64"));
2159 Fill values array with length values.
2161 template <class T> void FIELD<T>::getLength() const throw (MEDEXCEPTION)
2163 const char * LOC = "FIELD<double>::getLength() const : ";
2166 // The field has to be initilised by a non empty support and a
2167 // number of components = 1 and its value type has to be set to MED_REEL64
2168 // (ie a FIELD<double>)
2170 if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_REEL64))
2171 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to 1 and a value type set to MED_REEL64"));
2177 Fill values array with normal values.
2179 template <class T> void FIELD<T>::getNormal() const throw (MEDEXCEPTION)
2181 const char * LOC = "FIELD<double>::getNormal() const : ";
2184 // The field has to be initilised by a non empty support and a
2185 // number of components = 1 and its value type has to be set to MED_REEL64
2186 // (ie a FIELD<double>)
2188 if (_support == (SUPPORT *) NULL)
2189 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
2191 int dim_space = _support->getMesh()->getSpaceDimension();
2193 if ((_numberOfComponents != dim_space) || (_valueType != MED_REEL64))
2194 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
2200 Fill values array with barycenter values.
2202 template <class T> void FIELD<T>::getBarycenter() const throw (MEDEXCEPTION)
2204 const char * LOC = "FIELD<double>::getBarycenter() const : ";
2207 // The field has to be initilised by a non empty support and a number of
2208 //components = space dimension and its value type has to be set to MED_REEL64
2209 // (ie a FIELD<double>)
2211 if (_support == (SUPPORT *) NULL)
2212 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
2214 int dim_space = _support->getMesh()->getSpaceDimension();
2216 if ((_numberOfComponents != dim_space) || (_valueType != MED_REEL64))
2217 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The field has to be initialised with a non empty support, a number of components set to the space dimension and a value type set to MED_REEL64"));
2222 #endif /* FIELD_HXX */