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)
40 class FIELD_ // GENERIC POINTER TO a template <class T> class FIELD
60 Pointer to the support the field deals with.
63 const SUPPORT * _support ;
67 Number of field's components.
70 int _numberOfComponents ;
73 Number of field's values.
80 Array of size _numberOfComponents. /n
81 (constant, scalar, vector, tensor)/n
82 We could use an array of integer to store
85 - space dimension for vector,/n
86 - space dimension square for tensor./n
87 So numbers of values per entities would be
88 sum of _componentsTypes array.
90 Not implemented yet! All type are scalar !
93 int * _componentsTypes ;
96 Array of size _numberOfComponents
97 storing components names if any.
100 string * _componentsNames;
103 Array of size _numberOfComponents
104 storing components descriptions if any.
107 string * _componentsDescriptions;
110 Array of size _numberOfComponents
111 storing components units if any.
114 UNIT * _componentsUnits;
117 Array of size _numberOfComponents
118 storing components units if any.
121 string * _MEDComponentsUnits;
122 int _iterationNumber ;
126 // _valueType should be a static const. Here is an initialization exemple
127 // template < classType T > struct SET_VALUE_TYPE { static const med_type_champ _valueType = 0; }
128 // template < > struct SET_VALUE_TYPE<double> { static const med_type_champ _valueType = MED_EN::MED_REEL64; }
129 // template < > struct SET_VALUE_TYPE<int> { static const med_type_champ _valueType = MED_EN::MED_INT32; }
130 // static const med_type_champ _valueType = SET_VALUE_TYPE <T>::_valueType;
131 med_type_champ _valueType ;
133 vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
134 static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n ) throw (MEDEXCEPTION);
135 void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL) const throw (MEDEXCEPTION);
136 FIELD<double>* _getFieldSize() const;
140 friend class MED_MED_RDONLY_DRIVER;
141 friend class MED_MED_WRONLY_DRIVER;
142 friend class MED_MED_RDWR_DRIVER;
144 friend class VTK_MED_DRIVER;
153 FIELD_(const SUPPORT * Support, const int NumberOfComponents);
157 FIELD_(const FIELD_ &m);
164 // virtual void setIterationNumber (int IterationNumber);
165 // virtual void setOrderNumber (int OrderNumber);
166 // virtual void setFieldName (string& fieldName);
168 virtual void rmDriver(int index);
169 virtual int addDriver(driverTypes driverType,
170 const string & fileName="Default File Name.med",
171 const string & driverFieldName="Default Field Nam") ;
172 virtual int addDriver( GENDRIVER & driver);
173 virtual void read (const GENDRIVER &);
174 virtual void read(int index=0);
175 virtual void openAppend( void );
176 virtual void write(const GENDRIVER &);
177 virtual void write(int index=0, const string & driverName="");
179 virtual void writeAppend(const GENDRIVER &);
180 virtual void writeAppend(int index=0, const string & driverName="");
182 inline void setName(const string Name);
183 inline string getName() const;
184 inline void setDescription(const string Description);
185 inline string getDescription() const;
186 inline const SUPPORT * getSupport() const;
187 inline void setSupport(const SUPPORT * support);
188 inline void setNumberOfComponents(const int NumberOfComponents);
189 inline int getNumberOfComponents() const;
190 inline void setNumberOfValues(const int NumberOfValues);
191 inline int getNumberOfValues() const;
192 // inline void setComponentType(int *ComponentType);
193 // inline int * getComponentType() const;
194 // inline int getComponentTypeI(int i) const;
195 inline void setComponentsNames(const string * ComponentsNames);
196 inline void setComponentName(int i, const string ComponentName);
197 inline const string * getComponentsNames() const;
198 inline string getComponentName(int i) const;
199 inline void setComponentsDescriptions(const string * ComponentsDescriptions);
200 inline void setComponentDescription(int i, const string ComponentDescription);
201 inline const string * getComponentsDescriptions() const;
202 inline string getComponentDescription(int i) const;
204 // provisoire : en attendant de regler le probleme des unites !
205 inline void setComponentsUnits(const UNIT * ComponentsUnits);
206 inline const UNIT * getComponentsUnits() const;
207 inline const UNIT * getComponentUnit(int i) const;
208 inline void setMEDComponentsUnits(const string * MEDComponentsUnits);
209 inline void setMEDComponentUnit(int i, const string MEDComponentUnit);
210 inline const string * getMEDComponentsUnits() const;
211 inline string getMEDComponentUnit(int i) const;
213 inline void setIterationNumber(int IterationNumber);
214 inline int getIterationNumber() const;
215 inline void setTime(double Time);
216 inline double getTime() const;
217 inline void setOrderNumber(int OrderNumber);
218 inline int getOrderNumber() const;
220 inline void setValueType (const med_type_champ ValueType) ;
221 inline med_type_champ getValueType () const;
226 // ---------------------------------
227 // Implemented Methods : constructor
228 // ---------------------------------
229 using namespace MEDMEM;
237 inline void FIELD_::setName(const string Name)
244 inline string FIELD_::getName() const
249 Set FIELD description.
251 inline void FIELD_::setDescription(const string Description)
253 _description=Description;
256 Get FIELD description.
258 inline string FIELD_::getDescription() const
263 Set FIELD number of components.
265 inline void FIELD_::setNumberOfComponents(int NumberOfComponents)
267 _numberOfComponents=NumberOfComponents;
270 Get FIELD number of components.
272 inline int FIELD_::getNumberOfComponents() const
274 return _numberOfComponents ;
277 Set FIELD number of values.
279 It must be the same than in the associated SUPPORT object.
281 inline void FIELD_::setNumberOfValues(int NumberOfValues)
283 _numberOfValues=NumberOfValues;
286 Get FIELD number of value.
288 inline int FIELD_::getNumberOfValues() const
290 return _numberOfValues ;
293 // inline void FIELD_::setComponentType(int *ComponentType)
295 // _componentsTypes=ComponentType ;
297 // inline int * FIELD_::getComponentType() const
299 // return _componentsTypes ;
301 // inline int FIELD_::getComponentTypeI(int i) const
303 // return _componentsTypes[i-1] ;
307 Set FIELD components names.
309 Duplicate the ComponentsNames string array to put components names in
310 FIELD. ComponentsNames size must be equal to number of components.
312 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
314 if (NULL == _componentsNames)
315 _componentsNames = new string[_numberOfComponents] ;
316 for (int i=0; i<_numberOfComponents; i++)
317 _componentsNames[i]=ComponentsNames[i] ;
320 Set FIELD i^th component name.
322 i must be >=1 and <= number of components.
324 inline void FIELD_::setComponentName(int i, const string ComponentName)
326 _componentsNames[i-1]=ComponentName ;
329 Get a reference to the string array which contain the components names.
331 This Array size is equal to number of components
333 inline const string * FIELD_::getComponentsNames() const
335 return _componentsNames ;
338 Get the name of the i^th component.
340 inline string FIELD_::getComponentName(int i) const
342 return _componentsNames[i-1] ;
345 Set FIELD components descriptions.
347 Duplicate the ComponentsDescriptions string array to put components
348 descriptions in FIELD.
349 ComponentsDescriptions size must be equal to number of components.
351 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
353 if (NULL == _componentsDescriptions)
354 _componentsDescriptions = new string[_numberOfComponents] ;
355 for (int i=0; i<_numberOfComponents; i++)
356 _componentsDescriptions[i]=ComponentsDescriptions[i] ;
359 Set FIELD i^th component description.
361 i must be >=1 and <= number of components.
363 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
365 _componentsDescriptions[i-1]=ComponentDescription ;
368 Get a reference to the string array which contain the components descriptions.
370 This Array size is equal to number of components
372 inline const string * FIELD_::getComponentsDescriptions() const
374 return _componentsDescriptions ;
377 Get the description of the i^th component.
379 inline string FIELD_::getComponentDescription(int i) const
381 return _componentsDescriptions[i-1];
386 Set FIELD components UNIT.
388 Duplicate the ComponentsUnits UNIT array to put components
390 ComponentsUnits size must be equal to number of components.
392 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
394 if (NULL == _componentsUnits)
395 _componentsUnits = new UNIT[_numberOfComponents] ;
396 for (int i=0; i<_numberOfComponents; i++)
397 _componentsUnits[i]=ComponentsUnits[i] ;
400 Get a reference to the UNIT array which contain the components units.
402 This Array size is equal to number of components
404 inline const UNIT * FIELD_::getComponentsUnits() const
406 return _componentsUnits ;
409 Get the UNIT of the i^th component.
411 inline const UNIT * FIELD_::getComponentUnit(int i) const
413 return &_componentsUnits[i-1] ;
416 Set FIELD components unit.
418 Duplicate the MEDComponentsUnits string array to put components
420 MEDComponentsUnits size must be equal to number of components.
423 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
425 if (NULL == _MEDComponentsUnits)
426 _MEDComponentsUnits = new string[_numberOfComponents] ;
427 for (int i=0; i<_numberOfComponents; i++)
428 _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
431 Set FIELD i^th component unit.
433 i must be >=1 and <= number of components.
435 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
437 _MEDComponentsUnits[i-1]=MEDComponentUnit ;
440 Get a reference to the string array which contain the components units.
442 This Array size is equal to number of components
444 inline const string * FIELD_::getMEDComponentsUnits() const
446 return _MEDComponentsUnits ;
449 Get the string for unit of the i^th component.
451 inline string FIELD_::getMEDComponentUnit(int i) const
453 return _MEDComponentsUnits[i-1] ;
456 Set the iteration number where FIELD has been calculated.
458 inline void FIELD_::setIterationNumber(int IterationNumber)
460 _iterationNumber=IterationNumber;
463 Get the iteration number where FIELD has been calculated.
465 inline int FIELD_::getIterationNumber() const
467 return _iterationNumber ;
470 Set the time (in second) where FIELD has been calculated.
472 inline void FIELD_::setTime(double Time)
477 Get the time (in second) where FIELD has been calculated.
479 inline double FIELD_::getTime() const
484 Set the order number where FIELD has been calculated.
486 It corresponds to internal iteration during one time step.
488 inline void FIELD_::setOrderNumber(int OrderNumber)
490 _orderNumber=OrderNumber ;
493 Get the order number where FIELD has been calculated.
495 inline int FIELD_::getOrderNumber() const
497 return _orderNumber ;
500 Get a reference to the SUPPORT object associated to FIELD.
502 inline const SUPPORT * FIELD_::getSupport() const
507 Set the reference to the SUPPORT object associated to FIELD.
509 Reference is not duplicate, so it must not be deleted.
511 inline void FIELD_::setSupport(const SUPPORT * support)
516 Get the FIELD med value type (MED_INT32 or MED_REEL64).
518 inline med_type_champ FIELD_::getValueType () const
523 Set the FIELD med value type (MED_INT32 or MED_REEL64).
525 inline void FIELD_::setValueType (const med_type_champ ValueType)
527 _valueType = ValueType ;
530 /////////////////////////
531 // END OF CLASS FIELD_ //
532 /////////////////////////
536 This template class contains informations related with a FIELD :
537 - Values of the field
542 template <class T> class FIELD : public FIELD_
545 // ------- Drivers Management Part
548 //-----------------------//
550 //-----------------------//
553 virtual GENDRIVER * run(const string & fileName, FIELD<T> * ptrFIELD) const = 0 ;
556 //-------------------------------------------------------//
557 template <class T2> class INSTANCE_DE : public INSTANCE
558 //-------------------------------------------------------//
561 GENDRIVER * run(const string & fileName, FIELD<T> * ptrFIELD) const
563 return new T2(fileName,ptrFIELD);
567 // static INSTANCE_DE<MED_FIELD_RDONLY_DRIVER> inst_med_rdonly ;
568 static INSTANCE_DE<MED_FIELD_RDWR_DRIVER<T> > inst_med ;
569 static INSTANCE_DE<VTK_FIELD_DRIVER<T> > inst_vtk ;
571 static const INSTANCE * const instances[] ;
573 // ------ End of Drivers Management Part
575 // array of value of type T
576 MEDARRAY<T> *_value ;
579 void _operation(const FIELD& m,const FIELD& n, const medModeSwitch mode, char* Op);
580 void _operationInitialize(const FIELD& m,const FIELD& n, char* Op);
581 void _add_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
582 void _sub_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
583 void _mul_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
584 void _div_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
587 FIELD & operator=(const FIELD &m); // A FAIRE
590 FIELD(const FIELD &m);
591 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
592 FIELD(const SUPPORT * Support, driverTypes driverType,
593 const string & fileName="", const string & fieldName="",
594 const int iterationNumber = -1, const int orderNumber = -1)
595 throw (MEDEXCEPTION);
598 const FIELD operator+(const FIELD& m) const;
599 const FIELD operator-(const FIELD& m) const;
600 const FIELD operator*(const FIELD& m) const;
601 const FIELD operator/(const FIELD& m) const;
602 const FIELD operator-() const;
603 FIELD& operator+=(const FIELD& m);
604 FIELD& operator-=(const FIELD& m);
605 FIELD& operator*=(const FIELD& m);
606 FIELD& operator/=(const FIELD& m);
607 static FIELD* add(const FIELD& m, const FIELD& n);
608 static FIELD* sub(const FIELD& m, const FIELD& n);
609 static FIELD* mul(const FIELD& m, const FIELD& n);
610 static FIELD* div(const FIELD& m, const FIELD& n);
611 double normMax() const throw (MEDEXCEPTION);
612 double norm2() const throw (MEDEXCEPTION);
613 void applyLin(T a, T b);
614 template <T T_function(T)> void applyFunc();
615 static FIELD* scalarProduct(const FIELD& m, const FIELD& n);
616 double normL2(int component, const FIELD<double> * p_field_volume=NULL) const;
617 double normL2(const FIELD<double> * p_field_volume=NULL) const;
618 double normL1(int component, const FIELD<double> * p_field_volume=NULL) const;
619 double normL1(const FIELD<double> * p_field_volume=NULL) const;
621 friend class MED_FIELD_RDONLY_DRIVER<T>;
622 friend class MED_FIELD_WRONLY_DRIVER<T>;
624 friend class VTK_FIELD_DRIVER<T>;
625 //friend class MED_FIELD_RDWR_DRIVER <T>;
628 void rmDriver(int index=0);
629 int addDriver(driverTypes driverType,
630 const string & fileName="Default File Name.med",
631 const string & driverFieldName="Default Field Name") ;
632 int addDriver(GENDRIVER & driver);
634 void allocValue(const int NumberOfComponents);
635 void allocValue(const int NumberOfComponents, const int LengthValue);
639 inline void read(int index=0);
640 inline void read(const GENDRIVER & genDriver);
641 inline void write(int index=0, const string & driverName = "");
642 inline void write(const GENDRIVER &);
644 inline void writeAppend(int index=0, const string & driverName = "");
645 inline void writeAppend(const GENDRIVER &);
647 inline void setValue(MEDARRAY<T> *Value);
649 inline MEDARRAY<T>* getvalue() const;
650 inline const T* getValue(medModeSwitch Mode) const;
651 inline const T* getValueI(medModeSwitch Mode,int i) const;
652 inline T getValueIJ(int i,int j) const;
654 inline void setValue(medModeSwitch mode, T* value);
655 inline void setValueI(medModeSwitch mode, int i, T* value);
656 inline void setValueIJ(int i, int j, T value);
659 This fonction feeds the FIELD<double> private attributs _value with the
660 volume of each cells belonging to the argument Support. The field has to be
661 initialised via the constructor FIELD<double>(const SUPPORT * , const int )
662 with Support as SUPPORT argument, 1 has the number of components, and Support
663 has be a SUPPORT on 3D cells. This initialisation could be done by the empty
664 constructor followed by a setSupport and setNumberOfComponents call but it has
665 be followed by a setValueType(MED_REEL64) call.
667 void getVolume() const throw (MEDEXCEPTION) ;
669 This fonction feeds the FIELD<double> private attributs _value with the
670 area of each cells (or faces) belonging to the attribut _support. The field
671 has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
672 const int ) with 1 has the number of components, and _support has be a
673 SUPPORT on 2D cells or 3D faces. This initialisation could be done by the
674 empty constructor followed by a setSupport and setNumberOfComponents call but
675 it has be followed by a setValueType(MED_REEL64) call.
677 void getArea() const throw (MEDEXCEPTION) ;
679 This fonction feeds the FIELD<double> private attributs _value with the
680 length of each segments belonging to the attribut _support. The field has
681 to be initialised via the constructor FIELD<double>(const SUPPORT * ,
682 const int ) with 1 has the number of components, and _support has be a
683 SUPPORT on 3D edges or 2D faces. This initialisation could be done by the
684 empty constructor followed by a setSupport and setNumberOfComponents call but
685 it has be followed by a setValueType(MED_REEL64) call.
687 void getLength() const throw (MEDEXCEPTION) ;
689 This fonction feeds the FIELD<double> private attributs _value with the
690 normal vector of each faces belonging to the attribut _support. The field
691 has to be initialised via the constructor FIELD<double>(const SUPPORT * ,
692 const int ) with the space dimension has the number of components, and
693 _support has be a SUPPORT on 3D or 2D faces. This initialisation could be done
694 by the empty constructor followed by a setSupport and setNumberOfComponents
695 call but it has be followed by a setValueType(MED_REEL64) call.
697 void getNormal() const throw (MEDEXCEPTION) ;
699 This fonction feeds the FIELD<double> private attributs _value with the
700 barycenter of each faces or cells or edges belonging to the attribut _support.
701 The field has to be initialised via the constructor
702 FIELD<double>(const SUPPORT * ,const int ) with the space dimension has the
703 number of components, and _support has be a SUPPORT on 3D cells or 2D faces.
704 This initialisation could be done by the empty constructor followed by a
705 setSupport and setNumberOfComponents call but it has be followed by a
706 setValueType(MED_REEL64) call.
708 void getBarycenter() const throw (MEDEXCEPTION) ;
712 // --------------------
713 // Implemented Methods
714 // --------------------
715 using namespace MEDMEM;
718 Constructor with no parameter, most of the attribut members are set to NULL.
720 template <class T> FIELD<T>::FIELD():
721 _value((MEDARRAY<T>*)NULL)
723 MESSAGE("Constructeur FIELD sans parametre");
727 Constructor with parameters such that all attrribut are set but the _value
728 attribut is allocated but not set.
730 template <class T> FIELD<T>::FIELD(const SUPPORT * Support,
731 const int NumberOfComponents, const medModeSwitch Mode) throw (MEDEXCEPTION) :
732 FIELD_(Support, NumberOfComponents)
734 BEGIN_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
738 _numberOfValues = Support->getNumberOfElements(MED_ALL_ELEMENTS);
740 catch (MEDEXCEPTION &ex) {
741 MESSAGE("No value defined ! ("<<ex.what()<<")");
742 _value = (MEDARRAY<T>*)NULL ;
744 MESSAGE("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
745 if (0<_numberOfValues) {
746 _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues,Mode);
749 _value = (MEDARRAY<T>*)NULL ;
751 END_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
758 template <class T> void FIELD<T>::init ()
765 template <class T> FIELD<T>::FIELD(const FIELD & m):
768 MESSAGE("Constructeur FIELD de recopie");
769 if (m._value != NULL)
771 // copie only default !
772 _value = new MEDARRAY<T>(* m._value,false);
775 _value = (MEDARRAY<T> *) NULL;
776 //_drivers = m._drivers;
784 template <class T> FIELD<T> & FIELD<T>::operator=(const FIELD &m)
786 MESSAGE("Appel de FIELD<T>::operator=") ;
787 // operator= on FIELD_
793 Overload addition operator.
794 This operation is authorized only for compatible fields that have the same support.
795 The compatibility checking includes equality tests of the folowing data members:/n
797 - _numberOfComponents
800 - _MEDComponentsUnits.
802 The data members of the returned field are initialized, based on the first field, except for the name,
803 which is the combination of the two field's names and the operator.
804 Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> /n
805 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
806 When using python, this operator calls the copy constructor in any case.
807 The user has to be aware that when using operator + in associatives expressions like
808 <tt> a = b + c + d +e; </tt> /n
809 no optimisation is performed : the evaluation of last expression requires the construction of
813 const FIELD<T> FIELD<T>::operator+(const FIELD & m) const
815 BEGIN_OF("FIELD<T>::operator+(const FIELD & m)");
816 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
818 // Select mode : avoid if possible any calculation of other mode for fields m or *this
820 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
821 mode=m.getvalue()->getMode();
823 mode=this->getvalue()->getMode();
825 // Creation of the result - memory is allocated by FIELD constructor
826 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
827 //result._operation(*this,m,mode,"+"); // perform Atribute's initialization & addition
828 result._operationInitialize(*this,m,"+"); // perform Atribute's initialization
829 result._add_in_place(*this,m,mode); // perform addition
831 END_OF("FIELD<T>::operator+(const FIELD & m)");
835 /*! Overloaded Operator +=
836 * Operations are directly performed in the first field's array.
837 * This operation is authorized only for compatible fields that have the same support.
840 FIELD<T>& FIELD<T>::operator+=(const FIELD & m)
842 BEGIN_OF("FIELD<T>::operator+=(const FIELD & m)");
843 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
845 // We choose to keep *this mode, even if it may cost a re-calculation for m
846 medModeSwitch mode=this->getvalue()->getMode();
847 const T* value1=m.getValue(mode); // get pointers to the values we are adding
849 // get a non const pointer to the inside array of values and perform operation
850 T * value=const_cast<T *> (getValue(mode));
851 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
852 const T* endV=value+size; // pointer to the end of value
853 for(;value!=endV; value1++,value++)
855 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
856 this->getvalue()->clearOtherMode();
857 END_OF("FIELD<T>::operator+=(const FIELD & m)");
862 /*! Addition of fields. Static member function.
863 * The function return a pointer to a new created field that holds the addition.
864 * Data members are checked for compatibility and initialized.
865 * The user is in charge of memory deallocation.
868 FIELD<T>* FIELD<T>::add(const FIELD& m, const FIELD& n)
870 BEGIN_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
871 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
873 // Select mode : avoid if possible any calculation of other mode for fields m or *this
875 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
876 mode=m.getvalue()->getMode();
878 mode=n.getvalue()->getMode();
880 // Creation of a new field
881 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
882 result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
883 result->_add_in_place(m,n,mode); // perform addition
885 END_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
890 Overload substraction operator.
891 This operation is authorized only for compatible fields that have the same support.
892 The compatibility checking includes equality tests of the folowing data members:/n
894 - _numberOfComponents
897 - _MEDComponentsUnits.
899 The data members of the returned field are initialized, based on the first field, except for the name,
900 which is the combination of the two field's names and the operator.
901 Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> /n
902 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
903 When using python, this operator calls the copy constructor in any case.
904 The user has to be aware that when using operator - in associatives expressions like
905 <tt> a = b - c - d -e; </tt> /n
906 no optimisation is performed : the evaluation of last expression requires the construction of
910 const FIELD<T> FIELD<T>::operator-(const FIELD & m) const
912 BEGIN_OF("FIELD<T>::operator-(const FIELD & m)");
913 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
915 // Select mode : avoid if possible any calculation of other mode for fields m or *this
917 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
918 mode=m.getvalue()->getMode();
920 mode=this->getvalue()->getMode();
922 // Creation of the result - memory is allocated by FIELD constructor
923 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
924 //result._operation(*this,m,mode,"-"); // perform Atribute's initialization & substraction
925 result._operationInitialize(*this,m,"-"); // perform Atribute's initialization
926 result._sub_in_place(*this,m,mode); // perform substracion
928 END_OF("FIELD<T>::operator-(const FIELD & m)");
933 const FIELD<T> FIELD<T>::operator-() const
935 BEGIN_OF("FIELD<T>::operator-()");
937 medModeSwitch mode=this->getvalue()->getMode();
938 // Creation of the result - memory is allocated by FIELD constructor
939 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
940 // Atribute's initialization
941 result.setName("- "+getName());
942 result.setComponentsNames(getComponentsNames());
943 // not yet implemented setComponentType(getComponentType());
944 result.setComponentsDescriptions(getComponentsDescriptions());
945 result.setMEDComponentsUnits(getMEDComponentsUnits());
946 result.setComponentsUnits(getComponentsUnits());
947 result.setIterationNumber(getIterationNumber());
948 result.setTime(getTime());
949 result.setOrderNumber(getOrderNumber());
950 result.setValueType(getValueType());
952 const T* value1=getValue(mode);
953 // get a non const pointer to the inside array of values and perform operation
954 T * value=const_cast<T *> (result.getValue(mode));
955 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
956 const T* endV=value+size; // pointer to the end of value
958 for(;value!=endV; value1++,value++)
960 END_OF("FIELD<T>::operator-=(const FIELD & m)");
964 /*! Overloaded Operator -=
965 * Operations are directly performed in the first field's array.
966 * This operation is authorized only for compatible fields that have the same support.
969 FIELD<T>& FIELD<T>::operator-=(const FIELD & m)
971 BEGIN_OF("FIELD<T>::operator-=(const FIELD & m)");
972 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
974 // We choose to keep *this mode, even if it may cost a re-calculation for m
975 medModeSwitch mode=this->getvalue()->getMode();
976 const T* value1=m.getValue(mode); // get pointers to the values we are adding
978 // get a non const pointer to the inside array of values and perform operation
979 T * value=const_cast<T *> (getValue(mode));
980 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
981 const T* endV=value+size; // pointer to the end of value
983 for(;value!=endV; value1++,value++)
985 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
986 this->getvalue()->clearOtherMode();
987 END_OF("FIELD<T>::operator-=(const FIELD & m)");
992 /*! Substraction of fields. Static member function.
993 * The function return a pointer to a new created field that holds the substraction.
994 * Data members are checked for compatibility and initialized.
995 * The user is in charge of memory deallocation.
998 FIELD<T>* FIELD<T>::sub(const FIELD& m, const FIELD& n)
1000 BEGIN_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
1001 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1003 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1005 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1006 mode=m.getvalue()->getMode();
1008 mode=n.getvalue()->getMode();
1010 // Creation of a new field
1011 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1012 result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1013 result->_sub_in_place(m,n,mode); // perform substraction
1015 END_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
1020 Overload multiplication operator.
1021 This operation is authorized only for compatible fields that have the same support.
1022 The compatibility checking includes equality tests of the folowing data members:/n
1024 - _numberOfComponents
1027 - _MEDComponentsUnits.
1029 The data members of the returned field are initialized, based on the first field, except for the name,
1030 which is the combination of the two field's names and the operator.
1031 Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> /n
1032 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1033 When using python, this operator calls the copy constructor in any case.
1034 The user has to be aware that when using operator * in associatives expressions like
1035 <tt> a = b * c * d *e; </tt> /n
1036 no optimisation is performed : the evaluation of last expression requires the construction of
1040 const FIELD<T> FIELD<T>::operator*(const FIELD & m) const
1042 BEGIN_OF("FIELD<T>::operator*(const FIELD & m)");
1043 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1045 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1047 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
1048 mode=m.getvalue()->getMode();
1050 mode=this->getvalue()->getMode();
1052 // Creation of the result - memory is allocated by FIELD constructor
1053 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
1054 //result._operation(*this,m,mode,"*"); // perform Atribute's initialization & multiplication
1055 result._operationInitialize(*this,m,"*"); // perform Atribute's initialization
1056 result._mul_in_place(*this,m,mode); // perform multiplication
1058 END_OF("FIELD<T>::operator*(const FIELD & m)");
1062 /*! Overloaded Operator *=
1063 * Operations are directly performed in the first field's array.
1064 * This operation is authorized only for compatible fields that have the same support.
1067 FIELD<T>& FIELD<T>::operator*=(const FIELD & m)
1069 BEGIN_OF("FIELD<T>::operator*=(const FIELD & m)");
1070 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1072 // We choose to keep *this mode, even if it may cost a re-calculation for m
1073 medModeSwitch mode=this->getvalue()->getMode();
1074 const T* value1=m.getValue(mode); // get pointers to the values we are adding
1076 // get a non const pointer to the inside array of values and perform operation
1077 T * value=const_cast<T *> (getValue(mode));
1078 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1079 const T* endV=value+size; // pointer to the end of value
1081 for(;value!=endV; value1++,value++)
1083 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1084 this->getvalue()->clearOtherMode();
1085 END_OF("FIELD<T>::operator*=(const FIELD & m)");
1090 /*! Multiplication of fields. Static member function.
1091 * The function return a pointer to a new created field that holds the multiplication.
1092 * Data members are checked for compatibility and initialized.
1093 * The user is in charge of memory deallocation.
1096 FIELD<T>* FIELD<T>::mul(const FIELD& m, const FIELD& n)
1098 BEGIN_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1099 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1101 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1103 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1104 mode=m.getvalue()->getMode();
1106 mode=n.getvalue()->getMode();
1108 // Creation of a new field
1109 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1110 result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1111 result->_mul_in_place(m,n,mode); // perform multiplication
1113 END_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
1119 Overload division operator.
1120 This operation is authorized only for compatible fields that have the same support.
1121 The compatibility checking includes equality tests of the folowing data members:/n
1123 - _numberOfComponents
1126 - _MEDComponentsUnits.
1128 The data members of the returned field are initialized, based on the first field, except for the name,
1129 which is the combination of the two field's names and the operator.
1130 Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> /n
1131 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1132 When using python, this operator calls the copy constructor in any case.
1133 The user has to be aware that when using operator / in associatives expressions like
1134 <tt> a = b / c / d /e; </tt> /n
1135 no optimisation is performed : the evaluation of last expression requires the construction of
1139 const FIELD<T> FIELD<T>::operator/(const FIELD & m) const
1141 BEGIN_OF("FIELD<T>::operator/(const FIELD & m)");
1142 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1144 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1146 if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
1147 mode=m.getvalue()->getMode();
1149 mode=this->getvalue()->getMode();
1151 // Creation of the result - memory is allocated by FIELD constructor
1152 FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
1153 //result._operation(*this,m,mode,"/"); // perform Atribute's initialization & division
1154 result._operationInitialize(*this,m,"/"); // perform Atribute's initialization
1155 result._div_in_place(*this,m,mode); // perform division
1157 END_OF("FIELD<T>::operator/(const FIELD & m)");
1162 /*! Overloaded Operator /=
1163 * Operations are directly performed in the first field's array.
1164 * This operation is authorized only for compatible fields that have the same support.
1167 FIELD<T>& FIELD<T>::operator/=(const FIELD & m)
1169 BEGIN_OF("FIELD<T>::operator/=(const FIELD & m)");
1170 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1172 // We choose to keep *this mode, even if it may cost a re-calculation for m
1173 medModeSwitch mode=this->getvalue()->getMode();
1174 const T* value1=m.getValue(mode); // get pointers to the values we are adding
1176 // get a non const pointer to the inside array of values and perform operation
1177 T * value=const_cast<T *> (getValue(mode));
1178 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1179 const T* endV=value+size; // pointer to the end of value
1181 for(;value!=endV; value1++,value++)
1183 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1184 this->getvalue()->clearOtherMode();
1185 END_OF("FIELD<T>::operator/=(const FIELD & m)");
1190 /*! Division of fields. Static member function.
1191 * The function return a pointer to a new created field that holds the division.
1192 * Data members are checked for compatibility and initialized.
1193 * The user is in charge of memory deallocation.
1196 FIELD<T>* FIELD<T>::div(const FIELD& m, const FIELD& n)
1198 BEGIN_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1199 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1201 // Select mode : avoid if possible any calculation of other mode for fields m or *this
1203 if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
1204 mode=m.getvalue()->getMode();
1206 mode=n.getvalue()->getMode();
1208 // Creation of a new field
1209 FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
1210 result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1211 result->_div_in_place(m,n,mode); // perform division
1213 END_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
1220 This internal method initialize the members of a new field created to hold the result of the operation Op .
1221 Initialization is based on the first field, except for the name, which is the combination of the two field's names
1226 void FIELD<T>::_operationInitialize(const FIELD& m,const FIELD& n, char* Op)
1228 MESSAGE("Appel methode interne _add" << Op);
1230 // Atribute's initialization (copy of the first field's attributes)
1231 // Other data members (_support, _numberOfValues) are initialized in the field's constr.
1232 setName(m.getName()+" "+Op+" "+n.getName());
1233 setComponentsNames(m.getComponentsNames());
1234 // not yet implemented setComponentType(m.getComponentType());
1235 setComponentsDescriptions(m.getComponentsDescriptions());
1236 setMEDComponentsUnits(m.getMEDComponentsUnits());
1238 // The following data member may differ from field m to n.
1239 // The initialization is done based on the first field.
1240 setComponentsUnits(m.getComponentsUnits());
1241 setIterationNumber(m.getIterationNumber());
1242 setTime(m.getTime());
1243 setOrderNumber(m.getOrderNumber());
1244 setValueType(m.getValueType());
1250 Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
1251 This method is applied to a just created field with medModeSwitch mode.
1252 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1257 void FIELD<T>::_add_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1259 // get pointers to the values we are adding
1260 const T* value1=m.getValue(mode);
1261 const T* value2=n.getValue(mode);
1262 // get a non const pointer to the inside array of values and perform operation
1263 T * value=const_cast<T *> (getValue(mode));
1265 const int size=getNumberOfValues()*getNumberOfComponents();
1267 const T* endV1=value1+size;
1268 for(;value1!=endV1; value1++,value2++,value++)
1269 *value=(*value1)+(*value2);
1274 Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
1275 This method is applied to a just created field with medModeSwitch mode.
1276 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1281 void FIELD<T>::_sub_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1283 // get pointers to the values we are adding
1284 const T* value1=m.getValue(mode);
1285 const T* value2=n.getValue(mode);
1286 // get a non const pointer to the inside array of values and perform operation
1287 T * value=const_cast<T *> (getValue(mode));
1289 const int size=getNumberOfValues()*getNumberOfComponents();
1291 const T* endV1=value1+size;
1292 for(;value1!=endV1; value1++,value2++,value++)
1293 *value=(*value1)-(*value2);
1298 Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
1299 This method is applied to a just created field with medModeSwitch mode.
1300 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1305 void FIELD<T>::_mul_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1307 // get pointers to the values we are adding
1308 const T* value1=m.getValue(mode);
1309 const T* value2=n.getValue(mode);
1310 // get a non const pointer to the inside array of values and perform operation
1311 T * value=const_cast<T *> (getValue(mode));
1313 const int size=getNumberOfValues()*getNumberOfComponents();
1315 const T* endV1=value1+size;
1316 for(;value1!=endV1; value1++,value2++,value++)
1317 *value=(*value1)*(*value2);
1322 Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
1323 This method is applied to a just created field with medModeSwitch mode.
1324 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1329 void FIELD<T>::_div_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
1331 // get pointers to the values we are adding
1332 const T* value1=m.getValue(mode);
1333 const T* value2=n.getValue(mode);
1334 // get a non const pointer to the inside array of values and perform operation
1335 T * value=const_cast<T *> (getValue(mode));
1337 const int size=getNumberOfValues()*getNumberOfComponents();
1339 const T* endV1=value1+size;
1340 for(;value1!=endV1; value1++,value2++,value++)
1341 *value=(*value1)/(*value2);
1346 template <class T> double FIELD<T>::normMax() const throw (MEDEXCEPTION)
1348 const T* value=getValue(getvalue()->getMode()); // get pointer to the values
1349 const int size=getNumberOfValues()*getNumberOfComponents();
1350 if (size <= 0) // Size of array has to be strictly positive
1353 diagnosis="FIELD<T>::normMax() : cannot compute the norm of "+getName()+
1354 " : it size is non positive!";
1355 throw MEDEXCEPTION(diagnosis.c_str());
1357 const T* lastvalue=value+size; // get pointer just after last value
1358 const T* pMax=value; // pointer to the max value
1359 const T* pMin=value; // pointer to the min value
1361 // get pointers to the max & min value of array
1362 while ( ++value != lastvalue )
1364 if ( *pMin > *value )
1366 if ( *pMax < *value )
1370 T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
1371 T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
1373 return Max>Min ? static_cast<double>(Max) : static_cast<double>(Min);
1376 /*! Return Euclidien norm
1378 template <class T> double FIELD<T>::norm2() const throw (MEDEXCEPTION)
1380 const T* value=this->getValue(this->getvalue()->getMode()); // get const pointer to the values
1381 const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
1382 if (size <= 0) // Size of array has to be strictly positive
1385 diagnosis="FIELD<T>::norm2() : cannot compute the norm of "+getName()+
1386 " : it size is non positive!";
1387 throw MEDEXCEPTION(diagnosis.c_str());
1389 const T* lastvalue=value+size; // point just after last value
1391 T result((T)0); // init
1392 for( ; value!=lastvalue ; ++value)
1393 result += (*value) * (*value);
1395 return std::sqrt(static_cast<double> (result));
1399 /*! Apply to each (scalar) field component the template parameter T_function,
1400 * which is a pointer to function.
1401 * Since the pointer is known at compile time, the function is inlined into the inner loop!
1402 * calculation is done "in place".
1405 * \code myField.applyFunc<std::sqrt>(); // apply sqare root function \endcode
1406 * \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
1408 template <class T> template <T T_function(T)>
1409 void FIELD<T>::applyFunc()
1411 medModeSwitch mode=getvalue()->getMode();
1413 // get a non const pointer to the inside array of values and perform operation
1414 T * value=const_cast<T *> (getValue(mode));
1415 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1417 if (size>0) // for a negative size, there is nothing to do
1419 const T* lastvalue=value+size; // pointer to the end of value
1420 for(;value!=lastvalue; ++value) // apply linear transformation
1421 *value = T_function(*value);
1422 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1423 getvalue()->clearOtherMode();
1428 /*! Apply to each (scalar) field component the linear function x -> ax+b.
1429 * calculation is done "in place".
1431 template <class T> void FIELD<T>::applyLin(T a, T b)
1433 medModeSwitch mode=getvalue()->getMode();
1435 // get a non const pointer to the inside array of values and perform operation in place
1436 T * value=const_cast<T *> (getValue(mode));
1437 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1439 if (size>0) // for a negative size, there is nothing to do
1441 const T* lastvalue=value+size; // pointer to the end of value
1442 for(;value!=lastvalue; ++value) // apply linear transformation
1443 *value = a*(*value)+b;
1444 // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
1445 getvalue()->clearOtherMode();
1451 * Return a pointer to a new field that holds the scalar product. Static member function.
1452 * This operation is authorized only for compatible fields that have the same support.
1453 * The compatibility checking includes equality tests of the folowing data members:/n
1455 * - _numberOfComponents
1457 * - _componentsTypes
1458 * - _MEDComponentsUnits.
1459 * Data members are initialized.
1460 * The new field point to the same support and has one component.
1461 * Each value of it is the scalar product of the two argument's fields.
1462 * The user is in charge of memory deallocation.
1464 template <class T> FIELD<T>* FIELD<T>::scalarProduct(const FIELD & m, const FIELD & n)
1466 FIELD_::_checkFieldCompatibility( m, n); // may throw exception
1467 // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
1468 const medModeSwitch mode=MED_FULL_INTERLACE;
1470 const int numberOfElements=m.getNumberOfValues(); // strictly positive
1471 const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
1473 // Creation & init of a the result field on the same suppot, with one component
1474 FIELD<T>* result = new FIELD<T>(m.getSupport(),1,mode);
1475 result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
1476 result->setIterationNumber(m.getIterationNumber());
1477 result->setTime(m.getTime());
1478 result->setOrderNumber(m.getOrderNumber());
1479 result->setValueType(m.getValueType());
1481 const T* value1=m.getValue(mode); // get const pointer to the values
1482 const T* value2=n.getValue(mode); // get const pointer to the values
1483 // get a non const pointer to the inside array of values and perform operation
1484 T * value=const_cast<T *> (result->getValue(mode));
1486 const T* lastvalue=value+numberOfElements; // pointing just after last value of result
1487 for ( ; value!=lastvalue ; ++value ) // loop on all elements
1489 *value=(T)0; // initialize value
1490 const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
1491 for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
1492 *value += (*value1) * (*value2);
1497 /*! Return L2 Norm of the field's component.
1498 * Cannot be applied to a field with a support on nodes.
1499 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1501 template <class T> double FIELD<T>::normL2(int component, const FIELD<double> * p_field_volume) const
1503 _checkNormCompatibility(p_field_volume); // may throw exception
1504 if ( component<1 || component>getNumberOfComponents() )
1505 throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
1507 const FIELD<double> * p_field_size=p_field_volume;
1508 if(!p_field_volume) // if the user don't supply the volume
1509 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1511 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1512 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1513 const T* value=getValueI( MED_NO_INTERLACE, component); // get pointer to the component's values
1514 const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
1516 double integrale=0.0;
1518 for (; value!=lastvalue ; ++value ,++vol)
1520 integrale += static_cast<double>((*value) * (*value)) * (*vol);
1524 if(!p_field_volume) // if the user didn't supply the volume
1525 delete p_field_size; // delete temporary volume field
1527 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1529 return integrale/totVol;
1532 /*! Return L2 Norm of the field.
1533 * Cannot be applied to a field with a support on nodes.
1534 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1536 template <class T> double FIELD<T>::normL2(const FIELD<double> * p_field_volume) const
1538 _checkNormCompatibility(p_field_volume); // may throw exception
1539 const FIELD<double> * p_field_size=p_field_volume;
1540 if(!p_field_volume) // if the user don't supply the volume
1541 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1543 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1544 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1545 const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
1546 const T* value=getValue( MED_NO_INTERLACE); // get pointer to the field's values
1549 const double* p_vol=vol;
1550 for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
1553 double integrale=0.0;
1554 for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
1555 for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
1556 integrale += static_cast<double>((*value) * (*value)) * (*p_vol);
1558 if(!p_field_volume) // if the user didn't supply the volume
1559 delete p_field_size; // delete temporary volume field
1561 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1563 return integrale/totVol;
1566 /*! Return L1 Norm of the field's component.
1567 * Cannot be applied to a field with a support on nodes.
1568 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1570 template <class T> double FIELD<T>::normL1(int component, const FIELD<double> * p_field_volume) const
1572 _checkNormCompatibility(p_field_volume); // may throw exception
1573 if ( component<1 || component>getNumberOfComponents() )
1574 throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
1576 const FIELD<double> * p_field_size=p_field_volume;
1577 if(!p_field_volume) // if the user don't supply the volume
1578 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1580 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1581 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1582 const T* value=getValueI( MED_NO_INTERLACE, component); // get pointer to the component's values
1583 const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
1585 double integrale=0.0;
1587 for (; value!=lastvalue ; ++value ,++vol)
1589 integrale += std::abs( static_cast<double>(*value) ) * (*vol);
1593 if(!p_field_volume) // if the user didn't supply the volume
1594 delete p_field_size; // delete temporary volume field
1596 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1598 return integrale/totVol;
1601 /*! Return L1 Norm of the field.
1602 * Cannot be applied to a field with a support on nodes.
1603 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
1605 template <class T> double FIELD<T>::normL1(const FIELD<double> * p_field_volume) const
1607 _checkNormCompatibility(p_field_volume); // may throw exception
1608 const FIELD<double> * p_field_size=p_field_volume;
1609 if(!p_field_volume) // if the user don't supply the volume
1610 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
1612 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
1613 const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
1614 const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
1615 const T* value=getValue( MED_NO_INTERLACE); // get pointer to the field's values
1618 const double* p_vol=vol;
1619 for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
1622 double integrale=0.0;
1623 for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
1624 for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
1625 integrale += std::abs( static_cast<double>(*value) ) * (*p_vol);
1627 if(!p_field_volume) // if the user didn't supply the volume
1628 delete p_field_size; // delete temporary volume field
1630 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
1632 return integrale/totVol;
1639 Constructor with parameters; the object is set via a file and its associated
1640 driver. For the moment only the MED_DRIVER is considered and if the last two
1641 argument (iterationNumber and orderNumber) are not set; their default value
1642 is -1. If the field fieldDriverName with the iteration number
1643 iterationNumber and the order number orderNumber does not exist in the file
1644 fieldDriverName; the constructor raises an exception.
1646 template <class T> FIELD<T>::FIELD(const SUPPORT * Support,
1647 driverTypes driverType,
1648 const string & fileName/*=""*/,
1649 const string & fieldDriverName/*=""*/,
1650 const int iterationNumber,
1651 const int orderNumber)
1652 throw (MEDEXCEPTION)
1654 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) : ";
1663 _value = (MEDARRAY<T>*)NULL;
1665 _iterationNumber = iterationNumber;
1667 _orderNumber = orderNumber;
1673 MED_FIELD_RDONLY_DRIVER<T> myDriver(fileName,this);
1674 myDriver.setFieldName(fieldDriverName);
1675 current = addDriver(myDriver);
1678 // current = addDriver(driverType,fileName,fieldDriverName);
1679 // switch(_drivers[current]->getAccessMode() ) {
1680 // case MED_WRONLY : {
1681 // MESSAGE("FIELD<T>::FIELD(driverTypes driverType, .....) : driverType must have a MED_RDONLY or MED_RDWR accessMode");
1682 // rmDriver(current);
1689 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Driver type unknown : can't create driver!"));
1694 _drivers[current]->open();
1695 _drivers[current]->read();
1696 _drivers[current]->close();
1704 template <class T> FIELD<T>::~FIELD()
1706 BEGIN_OF(" Destructeur FIELD<T>::~FIELD()");
1708 if (_value) delete _value;
1709 END_OF(" Destructeur FIELD<T>::~FIELD()");
1715 template <class T> void FIELD<T>::allocValue(const int NumberOfComponents)
1717 const char* LOC = "FIELD<T>::allocValue(const int NumberOfComponents)" ;
1720 _numberOfComponents = NumberOfComponents ;
1721 if (_componentsTypes == NULL)
1722 _componentsTypes = new int[NumberOfComponents] ;
1723 if (_componentsNames == NULL)
1724 _componentsNames = new string[NumberOfComponents];
1725 if (_componentsDescriptions == NULL)
1726 _componentsDescriptions = new string[NumberOfComponents];
1727 if (_componentsUnits == NULL)
1728 _componentsUnits = new UNIT[NumberOfComponents];
1729 if (_MEDComponentsUnits == NULL)
1730 _MEDComponentsUnits = new string[NumberOfComponents];
1731 for (int i=0;i<NumberOfComponents;i++) {
1732 _componentsTypes[i] = 0 ;
1736 _numberOfValues = _support->getNumberOfElements(MED_ALL_ELEMENTS);
1737 MESSAGE(LOC <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
1739 _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
1743 catch (MEDEXCEPTION &ex) {
1744 MESSAGE("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
1745 _value = (MEDARRAY<T>*)NULL ;
1749 END_OF("void FIELD<T>::allocValue(const int NumberOfComponents)");
1755 template <class T> void FIELD<T>::allocValue(const int NumberOfComponents, const int LengthValue)
1757 BEGIN_OF("void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)");
1759 _numberOfComponents = NumberOfComponents ;
1760 if (_componentsTypes == NULL)
1761 _componentsTypes = new int[NumberOfComponents] ;
1762 if (_componentsNames == NULL)
1763 _componentsNames = new string[NumberOfComponents];
1764 if (_componentsDescriptions == NULL)
1765 _componentsDescriptions = new string[NumberOfComponents];
1766 if (_componentsUnits == NULL)
1767 _componentsUnits = new UNIT[NumberOfComponents];
1768 if (_MEDComponentsUnits == NULL)
1769 _MEDComponentsUnits = new string[NumberOfComponents];
1770 for (int i=0;i<NumberOfComponents;i++) {
1771 _componentsTypes[i] = 0 ;
1774 MESSAGE("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
1775 _numberOfValues = LengthValue ;
1776 _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
1780 END_OF("void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)");
1786 template <class T> void FIELD<T>::deallocValue()
1788 BEGIN_OF("void FIELD<T>::deallocValue()");
1789 _numberOfValues = 0 ;
1790 _numberOfComponents = 0 ;
1794 END_OF("void FIELD<T>::deallocValue()");
1797 // -----------------
1799 // -----------------
1802 template <class T> FIELD<T>::INSTANCE_DE<MED_FIELD_RDWR_DRIVER<T> > FIELD<T>::inst_med ;
1804 template <class T> FIELD<T>::INSTANCE_DE<VTK_FIELD_DRIVER<T> > FIELD<T>::inst_vtk ;
1806 template <class T> const typename FIELD<T>::INSTANCE * const FIELD<T>::instances[] = { &FIELD<T>::inst_med, &FIELD<T>::inst_vtk} ;
1810 Create the specified driver and return its index reference to path to
1811 read or write methods.
1814 template <class T> int FIELD<T>::addDriver(driverTypes driverType,
1815 const string & fileName/*="Default File Name.med"*/,
1816 const string & driverName/*="Default Field Name"*/)
1818 const char * LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Field Name\") : ";
1822 int itDriver = (int) NO_DRIVER;
1828 SCRUTE(instances[driverType]);
1833 itDriver = (int) driverType ;
1842 case GIBI_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"));
1848 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"));
1853 if (itDriver == ((int) NO_DRIVER))
1854 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"));
1856 driver = instances[itDriver]->run(fileName, this) ;
1858 _drivers.push_back(driver);
1860 current = _drivers.size()-1;
1862 _drivers[current]->setFieldName(driverName);
1871 Duplicate the given driver and return its index reference to path to
1872 read or write methods.
1874 template <class T> inline int FIELD<T>::addDriver (GENDRIVER & driver )
1876 const char * LOC = "FIELD<T>::addDriver(GENDRIVER &) : ";
1879 // duplicate driver to delete it with destructor !
1880 GENDRIVER * newDriver = driver.copy() ;
1882 _drivers.push_back(newDriver);
1883 return _drivers.size() -1 ;
1889 Remove the driver referenced by its index.
1891 template <class T> void FIELD<T>::rmDriver (int index/*=0*/)
1893 const char * LOC = "FIELD<T>::rmDriver (int index=0): ";
1896 if ( _drivers[index] ) {
1897 //_drivers.erase(&_drivers[index]);
1899 MESSAGE ("detruire");
1902 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1903 << "The <index given is invalid, index must be between 0 and |"
1912 Read FIELD in the file specified in the driver given by its index.
1914 template <class T> inline void FIELD<T>::read(int index/*=0*/)
1916 const char * LOC = "FIELD<T>::read(int index=0) : ";
1919 if ( _drivers[index] ) {
1920 _drivers[index]->open();
1921 _drivers[index]->read();
1922 _drivers[index]->close();
1925 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1926 << "The index given is invalid, index must be between 0 and |"
1934 Write FIELD in the file specified in the driver given by its index.
1936 template <class T> inline void FIELD<T>::write(int index/*=0*/, const string & driverName /*= ""*/)
1938 const char * LOC = "FIELD<T>::write(int index=0, const string & driverName = \"\") : ";
1941 if( _drivers[index] ) {
1942 _drivers[index]->open();
1943 if (driverName != "") _drivers[index]->setFieldName(driverName);
1944 _drivers[index]->write();
1945 _drivers[index]->close();
1948 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1949 << "The index given is invalid, index must be between 0 and |"
1957 Write FIELD in the file specified in the driver given by its index. Use this
1958 method for ASCII drivers (e.g. VTK_DRIVER)
1960 template <class T> inline void FIELD<T>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
1962 const char * LOC = "FIELD<T>::write(int index=0, const string & driverName = \"\") : ";
1965 if( _drivers[index] ) {
1966 _drivers[index]->openAppend();
1967 if (driverName != "") _drivers[index]->setFieldName(driverName);
1968 _drivers[index]->writeAppend();
1969 _drivers[index]->close();
1972 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
1973 << "The index given is invalid, index must be between 0 and |"
1982 Write FIELD with the driver which is equal to the given driver.
1986 template <class T> inline void FIELD<T>::write(const GENDRIVER & genDriver)
1988 const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
1991 for (unsigned int index=0; index < _drivers.size(); index++ )
1992 if ( *_drivers[index] == genDriver ) {
1993 _drivers[index]->open();
1994 _drivers[index]->write();
1995 _drivers[index]->close();
2004 Write FIELD with the driver which is equal to the given driver.
2006 Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
2008 template <class T> inline void FIELD<T>::writeAppend(const GENDRIVER & genDriver)
2010 const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
2013 for (unsigned int index=0; index < _drivers.size(); index++ )
2014 if ( *_drivers[index] == genDriver ) {
2015 _drivers[index]->openAppend();
2016 _drivers[index]->writeAppend();
2017 _drivers[index]->close();
2026 Read FIELD with the driver which is equal to the given driver.
2030 template <class T> inline void FIELD<T>::read(const GENDRIVER & genDriver)
2032 const char * LOC = " FIELD<T>::read(const GENDRIVER &) : ";
2035 for (unsigned int index=0; index < _drivers.size(); index++ )
2036 if ( *_drivers[index] == genDriver ) {
2037 _drivers[index]->open();
2038 _drivers[index]->read();
2039 _drivers[index]->close();
2048 Destroy the MEDARRAY<T> in FIELD and put the new one without copy.
2051 template <class T> inline void FIELD<T>::setValue(MEDARRAY<T> *Value)
2053 if (NULL != _value) delete _value ;
2059 Return a reference to the MEDARRAY<T> in FIELD.
2062 template <class T> inline MEDARRAY<T>* FIELD<T>::getvalue() const
2068 Return a reference to values array to read them.
2070 template <class T> inline const T* FIELD<T>::getValue(medModeSwitch Mode) const
2072 return _value->get(Mode) ;
2076 Return a reference to i^{th} row or column - component - (depend on Mode value)
2077 of FIELD values array.
2079 template <class T> inline const T* FIELD<T>::getValueI(medModeSwitch Mode,int i) const
2081 if ( Mode == MED_FULL_INTERLACE )
2083 return _value->getRow(i) ;
2085 ASSERT ( Mode == MED_NO_INTERLACE);
2086 return _value->getColumn(i);
2090 Return the value of i^{th} element and j^{th} component.
2092 template <class T> inline T FIELD<T>::getValueIJ(int i,int j) const
2094 return _value->getIJ(i,j) ;
2098 Copy new values array in FIELD according to the given mode.
2100 Array must have right size. If not results are unpredicable.
2102 template <class T> inline void FIELD<T>::setValue(medModeSwitch mode, T* value)
2104 _value->set(mode,value);
2108 Update values array in FIELD with the given ones according to specified mode.
2110 template <class T> inline void FIELD<T>::setValueI(medModeSwitch mode, int i, T* value)
2113 if (MED_FULL_INTERLACE == mode)
2114 _value->setI(i,value);
2115 else if (MED_NO_INTERLACE == mode)
2116 _value->setJ(i,value);
2118 throw MEDEXCEPTION(LOCALIZED("FIELD<T>::setValueI : bad medModeSwitch")) ;
2122 Set the value of i^{th} element and j^{th} component with the given one.
2124 template <class T> inline void FIELD<T>::setValueIJ(int i, int j, T value)
2126 _value->setIJ(i,j,value);
2134 Fill values array with volume values.
2136 template <class T> void FIELD<T>::getVolume() const throw (MEDEXCEPTION)
2138 const char * LOC = "FIELD<double>::getVolume() const : ";
2141 // The field has to be initilised by a non empty support and a
2142 // number of components = 1 and its value type has to be set to MED_REEL64
2143 // (ie a FIELD<double>)
2145 if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_REEL64))
2146 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"));
2152 Fill values array with area values.
2154 template <class T> void FIELD<T>::getArea() const throw (MEDEXCEPTION)
2156 const char * LOC = "FIELD<double>::getArea() const : ";
2159 // The field has to be initilised by a non empty support and a
2160 // number of components = 1 and its value type has to be set to MED_REEL64
2161 // (ie a FIELD<double>)
2163 if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_REEL64))
2164 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"));
2170 Fill values array with length values.
2172 template <class T> void FIELD<T>::getLength() const throw (MEDEXCEPTION)
2174 const char * LOC = "FIELD<double>::getLength() const : ";
2177 // The field has to be initilised by a non empty support and a
2178 // number of components = 1 and its value type has to be set to MED_REEL64
2179 // (ie a FIELD<double>)
2181 if ((_support == (SUPPORT *) NULL) || (_numberOfComponents != 1) || (_valueType != MED_REEL64))
2182 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"));
2188 Fill values array with normal values.
2190 template <class T> void FIELD<T>::getNormal() const throw (MEDEXCEPTION)
2192 const char * LOC = "FIELD<double>::getNormal() const : ";
2195 // The field has to be initilised by a non empty support and a
2196 // number of components = 1 and its value type has to be set to MED_REEL64
2197 // (ie a FIELD<double>)
2199 if (_support == (SUPPORT *) NULL)
2200 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"));
2202 int dim_space = _support->getMesh()->getSpaceDimension();
2204 if ((_numberOfComponents != dim_space) || (_valueType != MED_REEL64))
2205 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"));
2211 Fill values array with barycenter values.
2213 template <class T> void FIELD<T>::getBarycenter() const throw (MEDEXCEPTION)
2215 const char * LOC = "FIELD<double>::getBarycenter() const : ";
2218 // The field has to be initilised by a non empty support and a number of
2219 //components = space dimension and its value type has to be set to MED_REEL64
2220 // (ie a FIELD<double>)
2222 if (_support == (SUPPORT *) NULL)
2223 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"));
2225 int dim_space = _support->getMesh()->getSpaceDimension();
2227 if ((_numberOfComponents != dim_space) || (_valueType != MED_REEL64))
2228 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"));
2233 #endif /* FIELD_HXX */