1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
31 #include "MEDMEM_Utilities.hxx"
32 #include "MEDMEM_MedVersion.hxx"
33 #include "MEDMEM_Exception.hxx"
34 #include "MEDMEM_define.hxx"
35 #include "MEDMEM_Support.hxx"
36 #include "MEDMEM_Unit.hxx"
37 #include "MEDMEM_nArray.hxx"
38 #include "MEDMEM_GenDriver.hxx"
39 #include "MEDMEM_RCBase.hxx"
40 #include "MEDMEM_ArrayInterface.hxx"
41 #include "MEDMEM_SetInterlacingType.hxx"
42 #include "MEDMEM_FieldForward.hxx"
43 #include "MEDMEM_GaussLocalization.hxx"
44 #include "InterpKernelGaussCoords.hxx"
45 #include "PointLocator.hxx"
62 struct MinMax<double> {
63 static double getMin() { return DBL_MIN; }
64 static double getMax() { return DBL_MAX; }
69 static int getMin() { return INT_MIN; }
70 static int getMax() { return INT_MAX; }
73 template < typename T > struct SET_VALUE_TYPE {
74 static const MED_EN::med_type_champ _valueType = MED_EN::MED_UNDEFINED_TYPE;};
75 template < > struct SET_VALUE_TYPE<double> {
76 static const MED_EN::med_type_champ _valueType = MED_EN::MED_REEL64; };
77 template < > struct SET_VALUE_TYPE<int> {
78 static const MED_EN::med_type_champ _valueType = MED_EN::MED_INT32; };
80 /*!\defgroup FIELD_io Reading and writing files
82 Fields can be read or written to/from MED files.
86 For reading a field a typical use consists in :
87 - reading the mesh associated on which the field lies
88 - read the field, specifying its time step and order number
92 //reading mesh from file
93 MESH mesh(MED_DRIVER, "file.med", "my_Mesh");
94 //reading the field from the file
95 FIELD<double> field(group,MED_DRIVER,"file.med","my_Field",1,1,&mesh);
98 It is also possible to read a field without specifying its support. In this case, the field constructor
99 creates a support with no link to the initial mesh:
101 FIELD<double> field(MED_DRIVER, "file.med", "myField",1,1);
102 SUPPORT* support= field->getSupport();
105 See also \ref FIELD_constructors
109 When it comes to write fields, it is enough to call write() method.
110 A typical use will be :
113 mesh.write(MED_DRIVER, "myResultFile.med");
114 field.write(MED_DRIVER, "myResultFile.med");
117 \defgroup FIELD_constructors
119 The different field constructors correspond to the two main
120 ways a field is used :
121 - either it is read from a file to be consulted,
122 - or it can be created from scratch with a link to a support on which the values will be built.
124 \defgroup FIELD_algo Numerical operations on fields
125 This section groups together the different operators that enable the user to
126 treat the FIELD objects as high-level numerical arrays, giving operators for
127 numerical treatment (overloading of basic operators, algorithms, etc...)
129 \defgroup FIELD_getset Basic Get/Set operations
131 This sections groups together the basic operations
132 that describe access to all the elements constitutive of the description of the field :
134 - time iteration number(compulsory),
135 - inner loop iteration number(compulsory),
137 - description(optional),
138 - number of components(compulsory),
139 - components names(optional),
140 - components description(optional).
142 Some of these items are compulsory because they are essential to the field in order to define
143 its structure or to be identified inside a MED file during the write process. The other ones
144 are there for additional information and can be overlooked if not necessary.
146 When creating a field by reading a file, all the parameters are set according to the file
147 data and can be consulted via the get methods. When creating a file from scratch, the
148 name and number of components are set by the constructor, but the other items have to be
149 set via the setXXX methods.
151 \defgroup FIELD_gauss Gauss points
153 In MED, it is possible to declare a Gauss model that
154 describes the location of Gauss points in a reference cell.
155 This Gauss model definition is contained in the
156 \a GAUSS_LOCALIZATION class. A \a GAUSS_LOCALIZATION object
157 is associated to a field and to a type.
159 It is not permitted to define a Gauss model in a polygonal
160 or polyhedric element.
162 The Gauss model can be :
163 - loaded from a MED file,
164 - written to a MED file,
165 - used to define a FIELD with multiple localizations per element.
167 \section gauss_constructors Constructing a Gauss Model
169 A Gauss model can be constructed with the following constructor :
170 \param locName defines a name associated with the gauss model
171 \param typeGeo names the type to which the Gauss model is assocaited
172 \param nGauss defines the number of Gauss points
173 \param cooRef defines an array giving the coordinates of the nodes of the reference element (dimension : spaceDimension * number of nodes for type \a typeGeo)
174 \param cooGauss defines an array giving the coordinates of the nodes of the Gauss points (dimension : spaceDimension * \a nGauss
176 \param wg weights associated with each Gauss point (dimension : \a nGauss)
178 Example : in 2D, a Gauss model definition for a triangle
179 would be written as :
182 string locname("gauss model");
183 double cooRef[6] ={0.0, 0.0, 1.0, 0.0, 0.0, 1.0};
184 double cooGauss[6]={0.2, 0.2, 0.8, 0.1, 0.1, 0.8};
185 double wg[3]={0.3334, 0.3334, 0.3334};
186 GAUSS_LOCALIZATION model(locname,
199 This class contains all the informations related with a template class FIELD :
200 - Components descriptions
201 - Time step description
202 - Location of the values (a SUPPORT class)
205 class MEDMEM_EXPORT FIELD_ : public RCBASE // GENERIC POINTER TO a template <class T, class INTERLACING_TAG> class FIELD
223 string _description ;
226 Pointer to the support the field deals with.
229 const SUPPORT * _support ;
233 Number of field's components.
236 int _numberOfComponents ;
239 Number of field's values.
240 doesn't take care of _numberOfComponents
241 and number of Gauss points.
244 int _numberOfValues ;
248 Array of size _numberOfComponents. \n
249 (constant, scalar, vector, tensor)\n
250 We could use an array of integer to store
251 numbers of values: \n
253 - space dimension for vector,\n
254 - space dimension square for tensor.\n
255 So numbers of values per entities would be
256 sum of _componentsTypes array.
258 Not implemented yet! All type are scalar !
261 vector<int> _componentsTypes ;
264 Array of size _numberOfComponents
265 storing components names if any.
268 vector<string> _componentsNames;
271 Array of size _numberOfComponents
272 storing components descriptions if any.
275 vector<string> _componentsDescriptions;
278 Array of size _numberOfComponents
279 storing components units if any.
282 vector<UNIT> _componentsUnits;
285 Array of size _numberOfComponents
286 storing components units if any.
289 vector<string> _MEDComponentsUnits;
292 Iteration number of the field.
295 int _iterationNumber ;
304 Order number of the field.
310 At the initialization step of the field using the constructors; this attribute,
311 the value type (integer or real) , is set automatically. There is a get method
312 but not a set method for this attribute.
315 MED_EN::med_type_champ _valueType ;
318 At the initialization step of the field using the constructors; this attribute,
319 the interlacing type (full interlace or no interlace field value storage), is set
320 automatically. There is a get method but not a set method for this attribute.
323 MED_EN::medModeSwitch _interlacingType;
325 vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
326 static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
327 static void _deepCheckFieldCompatibility(const FIELD_& m, const FIELD_& n, bool checkUnit=true) throw (MEDEXCEPTION);
328 void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL,
329 const bool nodalAllowed = false) const throw (MEDEXCEPTION);
330 FIELD<double>* _getFieldSize(const SUPPORT *subSupport=NULL) const;
340 FIELD_(const SUPPORT * Support, const int NumberOfComponents);
345 FIELD_(const FIELD_ &m);
354 friend class MED_MED_RDONLY_DRIVER22;
355 friend class MED_MED_WRONLY_DRIVER22;
356 friend class MED_MED_RDWR_DRIVER22;
357 friend class VTK_MED_DRIVER;
359 FIELD_& operator=(const FIELD_ &m);
361 virtual void rmDriver(int index=0);
369 /*! Creates a driver for reading/writing fields in a file.
370 \param driverType specifies the file type (MED_DRIVER, VTK_DRIVER)
371 \param fileName name of the output file
372 \param driverFieldName name of the field
373 \param access specifies whether the file is opened for read, write or both.
376 virtual int addDriver(driverTypes driverType,
377 const string & fileName="Default File Name.med",
378 const string & driverFieldName="Default Field Nam",
379 MED_EN::med_mode_acces access=MED_EN::RDWR) ;
381 virtual int addDriver( GENDRIVER & driver);
382 virtual void read (driverTypes driverType, const std::string & fileName);
383 virtual void read (const GENDRIVER &);
384 virtual void read(int index=0);
385 virtual void openAppend( void );
386 virtual void write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
387 virtual void write(driverTypes driverType,
388 const std::string & fileName,
389 MED_EN::med_mode_acces medMode=MED_EN::RDWR);
391 /*! Triggers the writing of the field with respect to the driver handle
392 \a index given by \a addDriver(...) method. */
393 virtual void write(int index=0);
394 /*!\if MEDMEM_ug @} \endif */
396 virtual void writeAppend(const GENDRIVER &);
397 virtual void writeAppend(int index=0, const string & driverName="");
399 inline void setName(const string Name);
400 inline string getName() const;
401 inline void setDescription(const string Description);
402 inline string getDescription() const;
403 inline const SUPPORT * getSupport() const;
404 inline void setSupport(const SUPPORT * support);
405 inline void setNumberOfComponents(const int NumberOfComponents);
406 inline int getNumberOfComponents() const;
407 inline void setNumberOfValues(const int NumberOfValues);
408 inline int getNumberOfValues() const;
409 inline void setComponentsNames(const string * ComponentsNames);
410 inline void setComponentName(int i, const string ComponentName);
411 inline const string * getComponentsNames() const;
412 inline string getComponentName(int i) const;
413 inline void setComponentsDescriptions(const string * ComponentsDescriptions);
414 inline void setComponentDescription(int i, const string ComponentDescription);
415 inline const string * getComponentsDescriptions() const;
416 inline string getComponentDescription(int i) const;
418 // provisoire : en attendant de regler le probleme des unites !
419 inline void setComponentsUnits(const UNIT * ComponentsUnits);
420 inline const UNIT * getComponentsUnits() const;
421 inline const UNIT * getComponentUnit(int i) const;
422 inline void setMEDComponentsUnits(const string * MEDComponentsUnits);
423 inline void setMEDComponentUnit(int i, const string MEDComponentUnit);
424 inline const string * getMEDComponentsUnits() const;
425 inline string getMEDComponentUnit(int i) const;
427 inline void setIterationNumber(int IterationNumber);
428 inline int getIterationNumber() const;
429 inline void setTime(double Time);
430 inline double getTime() const;
431 inline void setOrderNumber(int OrderNumber);
432 inline int getOrderNumber() const;
434 inline MED_EN::med_type_champ getValueType () const;
435 inline MED_EN::medModeSwitch getInterlacingType() const;
436 virtual inline bool getGaussPresence() const throw (MEDEXCEPTION);
438 void copyGlobalInfo(const FIELD_& m);
445 \addtogroup FIELD_getset
450 Sets FIELD name. The length should not exceed MED_TAILLE_NOM
451 as defined in Med (i.e. 32 characters).
453 inline void FIELD_::setName(const string Name)
460 inline string FIELD_::getName() const
465 Sets FIELD description. The length should not exceed MED_TAILLE_DESC as defined in Med (i.e. 200 characters).
467 inline void FIELD_::setDescription(const string Description)
469 _description=Description;
472 Gets FIELD description.
474 inline string FIELD_::getDescription() const
479 Sets FIELD number of components.
481 inline void FIELD_::setNumberOfComponents(const int NumberOfComponents)
483 _numberOfComponents=NumberOfComponents;
484 _componentsTypes.resize(_numberOfComponents);
485 _componentsNames.resize(_numberOfComponents);
486 _componentsDescriptions.resize(_numberOfComponents);
487 _componentsUnits.resize(_numberOfComponents);
488 _MEDComponentsUnits.resize(_numberOfComponents);
491 Gets FIELD number of components.
493 inline int FIELD_::getNumberOfComponents() const
495 return _numberOfComponents ;
498 Sets FIELD number of values.
500 It must be the same than in the associated SUPPORT object.
502 inline void FIELD_::setNumberOfValues(const int NumberOfValues)
504 _numberOfValues=NumberOfValues;
507 Gets FIELD number of value.
509 inline int FIELD_::getNumberOfValues() const
511 return _numberOfValues ;
515 Sets FIELD components names.
517 Duplicates the ComponentsNames string array to put components names in
518 FIELD. ComponentsNames size must be equal to number of components.
520 inline void FIELD_::setComponentsNames(const string * ComponentsNames)
522 _componentsNames.resize(_numberOfComponents);
523 for (int i=0; i<_numberOfComponents; i++)
524 _componentsNames[i]=ComponentsNames[i] ;
527 Sets FIELD i^th component name.
529 i must be >=1 and <= number of components.
532 inline void FIELD_::setComponentName(int i, const string ComponentName)
534 const char * LOC = " FIELD_::setComponentName() : ";
536 if( i<1 || i>_numberOfComponents )
537 throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
539 _componentsNames[i-1]=ComponentName ;
542 Gets a reference to the string array which contain the components names.
544 This Array size is equal to number of components
546 inline const string * FIELD_::getComponentsNames() const
548 return &(_componentsNames[0]) ;
551 Gets the name of the i^th component.
554 inline string FIELD_::getComponentName(int i) const
556 const char * LOC = " FIELD_::getComponentName() : ";
558 if( i<1 || i>_numberOfComponents )
559 throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
561 return _componentsNames[i-1] ;
564 Sets FIELD components descriptions.
566 Duplicates the ComponentsDescriptions string array to put components
567 descriptions in FIELD.
568 ComponentsDescriptions size must be equal to number of components.
570 inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
572 _componentsDescriptions.resize(_numberOfComponents);
573 for (int i=0; i<_numberOfComponents; i++)
574 _componentsDescriptions[i]=ComponentsDescriptions[i] ;
577 Sets FIELD i^th component description.
579 i must be >=1 and <= number of components.
582 inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
584 const char * LOC = " FIELD_::setComponentDescription() : ";
586 if( i<1 || i>_numberOfComponents )
587 throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
589 _componentsDescriptions[i-1]=ComponentDescription ;
592 Gets a reference to the string array which contain the components descriptions.
594 This Array size is equal to number of components
596 inline const string * FIELD_::getComponentsDescriptions() const
598 return &(_componentsDescriptions[0]);
601 Gets the description of the i^th component.
604 inline string FIELD_::getComponentDescription(int i) const
606 const char * LOC = " FIELD_::setComponentDescription() : ";
608 if( i<1 || i>_numberOfComponents )
609 throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
611 return _componentsDescriptions[i-1];
616 Sets FIELD components UNIT.
618 Duplicates the ComponentsUnits UNIT array to put components
620 ComponentsUnits size must be equal to number of components.
623 inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
625 _componentsUnits.resize(_numberOfComponents);
626 for (int i=0; i<_numberOfComponents; i++)
627 _componentsUnits[i]=ComponentsUnits[i] ;
630 Gets a reference to the UNIT array which contain the components units.
632 This array size is equal to number of components
635 inline const UNIT * FIELD_::getComponentsUnits() const
637 return &(_componentsUnits[0]);
640 Gets the UNIT of the i^th component.
643 inline const UNIT * FIELD_::getComponentUnit(int i) const
645 const char * LOC = " FIELD_::getComponentUnit() : ";
647 if( i<1 || i>_numberOfComponents )
648 throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
650 return &_componentsUnits[i-1] ;
653 Sets FIELD components unit.
655 Duplicates the MEDComponentsUnits string array to put components
657 MEDComponentsUnits size must be equal to number of components.
660 inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
662 _MEDComponentsUnits.resize(_numberOfComponents);
663 for (int i=0; i<_numberOfComponents; i++)
664 _MEDComponentsUnits[i]=MEDComponentsUnits[i] ;
667 Sets FIELD i^th component unit.
669 i must be >=1 and <= number of components.
672 inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
674 const char * LOC = " FIELD_::setMEDComponentUnit() : ";
676 if( i<1 || i>_numberOfComponents )
677 throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
679 _MEDComponentsUnits[i-1]=MEDComponentUnit ;
682 Gets a reference to the string array which contain the components units.
684 This array size is equal to number of components
686 inline const string * FIELD_::getMEDComponentsUnits() const
688 return &(_MEDComponentsUnits[0]);
691 Gets the string for unit of the i^th component.
694 inline string FIELD_::getMEDComponentUnit(int i) const
696 const char * LOC = " FIELD_::getMEDComponentUnit() : ";
698 if( i<1 || i>_numberOfComponents )
699 throw MEDEXCEPTION(STRING(LOC)<<" invalid index" );
701 return _MEDComponentsUnits[i-1] ;
704 Sets the iteration number where FIELD has been calculated.
706 inline void FIELD_::setIterationNumber(int IterationNumber)
708 _iterationNumber=IterationNumber;
711 Gets the iteration number where FIELD has been calculated.
713 inline int FIELD_::getIterationNumber() const
715 return _iterationNumber ;
718 Sets the time when FIELD has been calculated.
720 inline void FIELD_::setTime(double Time)
725 Gets the time when FIELD has been calculated.
727 inline double FIELD_::getTime() const
732 Sets the order number where FIELD has been calculated.
734 It corresponds to internal iteration during one time step.
736 inline void FIELD_::setOrderNumber(int OrderNumber)
738 _orderNumber=OrderNumber ;
741 Gets the order number where FIELD has been calculated.
743 inline int FIELD_::getOrderNumber() const
745 return _orderNumber ;
748 Gets a reference to the SUPPORT object associated to FIELD.
750 inline const SUPPORT * FIELD_::getSupport() const
755 Sets the reference to the SUPPORT object associated to FIELD.
757 Reference is not duplicate, so it must not be deleted.
759 inline void FIELD_::setSupport(const SUPPORT * support)
761 //A.G. Addings for RC
762 if(_support!=support)
765 _support->removeReference();
768 _support->addReference();
772 Gets the FIELD med value type (MED_INT32 or MED_REEL64).
774 inline MED_EN::med_type_champ FIELD_::getValueType () const
780 Gets the FIELD med interlacing type (MED_FULL_INTERLACE or MED_NO_INTERLACE).
782 inline MED_EN::medModeSwitch FIELD_::getInterlacingType () const
784 return _interlacingType ;
786 /*!\if MEDMEM_ug @} \endif*/
789 \addtogroup FIELD_gauss
794 Determines whether the field stores several Gauss points per element.
796 inline bool FIELD_::getGaussPresence() const throw (MEDEXCEPTION)
798 const char * LOC = "FIELD_::getGaussPresence() : ";
799 throw MEDEXCEPTION(STRING(LOC) << " This FIELD_ doesn't rely on a FIELD<T>" );
802 /*!\if MEDMEM_ug @} \endif*/
804 } //End namespace MEDMEM
806 /////////////////////////
807 // END OF CLASS FIELD_ //
808 /////////////////////////
812 This template class contains informations related with a FIELD :
813 - Values of the field, their type (real or integer), the storage mode (full interlace,
814 no interlace or no interlace by type).
821 template<class T2> class MED_FIELD_RDONLY_DRIVER;
822 template<class T2> class MED_FIELD_WRONLY_DRIVER;
823 template<class T2> class VTK_FIELD_DRIVER;
827 class INTERLACING_TAG
828 > class FIELD : public FIELD_
832 typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array ArrayNoGauss;
833 typedef typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array ArrayGauss;
834 typedef typename MEDMEM_ArrayInterface<T,NoInterlace,NoGauss>::Array ArrayNo;
835 typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array ArrayFull;
836 typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,NoGauss>::Array ArrayNoByType;
837 typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,Gauss>::Array ArrayNoByTypeGauss;
838 typedef MEDMEM_Array_ Array;
839 typedef T ElementType;
840 typedef INTERLACING_TAG InterlacingTag;
841 typedef map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> locMap;
843 // array of value of type T
846 // MESH, to be used for field reading from a file (if desired to link
847 // to existing support instead of new support creation for the field)
854 map<MED_EN::medGeometryElement,GAUSS_LOCALIZATION_*> _gaussModel; //A changer quand les drivers seront template de l'entrelacement
856 static T _scalarForPow;
860 void _operation(const FIELD& m,const FIELD& n, const char* Op);
861 void _operationInitialize(const FIELD& m,const FIELD& n, const char* Op);
862 void _add_in_place(const FIELD& m,const FIELD& n);
863 void _sub_in_place(const FIELD& m,const FIELD& n);
864 void _mul_in_place(const FIELD& m,const FIELD& n);
865 void _div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION);
868 FIELD(const FIELD &m);
869 FIELD(const SUPPORT * Support, const int NumberOfComponents) throw (MEDEXCEPTION);
870 FIELD(driverTypes driverType,
871 const string & fileName, const string & fieldDriverName,
872 const int iterationNumber=-1, const int orderNumber=-1,
874 throw (MEDEXCEPTION);
875 FIELD(const SUPPORT * Support, driverTypes driverType,
876 const string & fileName="", const string & fieldName="",
877 const int iterationNumber = -1, const int orderNumber = -1)
878 throw (MEDEXCEPTION);
882 FIELD & operator=(const FIELD &m);
883 FIELD & operator=(T value);
884 FIELD *operator+(const FIELD& m) const;
885 FIELD *operator-(const FIELD& m) const;
886 FIELD *operator*(const FIELD& m) const;
887 FIELD *operator/(const FIELD& m) const;
888 FIELD *operator-() const;
889 FIELD& operator+=(const FIELD& m);
890 FIELD& operator-=(const FIELD& m);
891 FIELD& operator*=(const FIELD& m);
892 FIELD& operator/=(const FIELD& m);
894 void applyLin(T a, T b, int icomp);
895 static FIELD* add(const FIELD& m, const FIELD& n);
896 static FIELD* addDeep(const FIELD& m, const FIELD& n);
897 static FIELD* sub(const FIELD& m, const FIELD& n);
898 static FIELD* subDeep(const FIELD& m, const FIELD& n);
899 static FIELD* mul(const FIELD& m, const FIELD& n);
900 static FIELD* mulDeep(const FIELD& m, const FIELD& n);
901 static FIELD* div(const FIELD& m, const FIELD& n);
902 static FIELD* divDeep(const FIELD& m, const FIELD& n);
903 double normMax() const throw (MEDEXCEPTION);
905 //------- TDG and BS addings
907 void getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION);
908 vector<int> getHistogram(int &nbint) throw (MEDEXCEPTION);
909 FIELD<double>* buildGradient() const throw (MEDEXCEPTION);
910 FIELD<double>* buildNorm2Field() const throw (MEDEXCEPTION);
912 //-------------------
914 double norm2() const throw (MEDEXCEPTION);
915 void applyLin(T a, T b);
916 template <T T_function(T)> void applyFunc();
917 void applyPow(T scalar);
918 static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
919 double normL2(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
920 double normL2(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
921 double normL1(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
922 double normL1(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
923 double integral(const SUPPORT *subSupport=NULL) const throw (MEDEXCEPTION);
924 FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
926 friend class MED_FIELD_RDONLY_DRIVER<T>;
927 friend class MED_FIELD_WRONLY_DRIVER<T>;
928 friend class VTK_FIELD_DRIVER<T>;
931 void rmDriver(int index=0);
932 int addDriver(driverTypes driverType,
933 const string & fileName="Default File Name.med",
934 const string & driverFieldName="Default Field Name",
935 MED_EN::med_mode_acces access=MED_EN::RDWR) ;
937 int addDriver(GENDRIVER & driver);
939 void allocValue(const int NumberOfComponents);
940 void allocValue(const int NumberOfComponents, const int LengthValue);
944 inline void read(int index=0);
945 inline void read(const GENDRIVER & genDriver);
946 inline void read(driverTypes driverType, const std::string& filename);
947 inline void write(int index=0);
948 inline void write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
949 inline void write(driverTypes driverType,
950 const std::string& filename,
951 MED_EN::med_mode_acces medMode=MED_EN::RDWR);
953 inline void writeAppend(int index=0, const string & driverName = "");
954 inline void writeAppend(const GENDRIVER &);
956 inline MEDMEM_Array_ * getArray() const throw (MEDEXCEPTION);
957 inline ArrayGauss * getArrayGauss() const throw (MEDEXCEPTION);
958 inline ArrayNoGauss * getArrayNoGauss() const throw (MEDEXCEPTION);
959 inline bool getGaussPresence() const throw (MEDEXCEPTION);
961 inline int getValueLength() const throw (MEDEXCEPTION);
964 \addtogroup FIELD_value
967 /*! Returns a pointer to the value array.*/
968 inline const T* getValue() const throw (MEDEXCEPTION);
969 inline const T* getRow(int i) const throw (MEDEXCEPTION);
970 inline const T* getColumn(int j) const throw (MEDEXCEPTION);
972 Returns the value of \f$ i^{th} \f$ element and \f$ j^{th}\f$ component.
973 This method only works with fields having no particular Gauss point
974 definition (i.e., fields having one value per element).
975 This method makes the retrieval of the value independent from the
976 interlacing pattern, but it is slower than the complete retrieval
977 obtained by the \b getValue() method.
980 inline T getValueIJ(int i,int j) const throw (MEDEXCEPTION);
983 Returns the \f$ j^{th}\f$ component of \f$ k^{th}\f$ Gauss points of \f$ i^{th}\f$ value.
984 This method is compatible with elements having more than one Gauss point.
985 This method makes the retrieval of the value independent from the
986 interlacing pattern, but it is slower than the complete retrieval
987 obtained by the \b getValue() method.
989 inline T getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION);
991 inline int getValueByTypeLength(int t) const throw (MEDEXCEPTION);
992 inline const T* getValueByType(int t) const throw (MEDEXCEPTION);
993 inline T getValueIJByType(int i,int j,int t) const throw (MEDEXCEPTION);
994 inline T getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION);
997 The following example describes the creation of a FIELD.
999 \example FIELDcreate.cxx
1001 \if MEDMEM_ug @} \endif */
1003 bool getValueOnElement(int eltIdInSup,T* retValues) const throw (MEDEXCEPTION);
1004 void getValueOnPoint(const double* coords, double* output) const throw (MEDEXCEPTION);
1005 void getValueOnPoints(int nb_points, const double* coords, double* output) const throw (MEDEXCEPTION);
1007 const int getNumberOfGeometricTypes() const throw (MEDEXCEPTION);
1008 const GAUSS_LOCALIZATION<INTERLACING_TAG> & getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1009 const GAUSS_LOCALIZATION<INTERLACING_TAG> * getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1010 const GAUSS_LOCALIZATION_* getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1011 void setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG> & gaussloc);
1012 void setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc);
1013 const int * getNumberOfGaussPoints() const throw (MEDEXCEPTION);
1014 const int getNumberOfGaussPoints( MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
1015 const int getNbGaussI(int i) const throw (MEDEXCEPTION);
1016 const int * getNumberOfElements() const throw (MEDEXCEPTION);
1017 const MED_EN::medGeometryElement * getGeometricTypes() const throw (MEDEXCEPTION);
1018 bool isOnAllElements() const throw (MEDEXCEPTION);
1020 inline void setArray(MEDMEM_Array_ *value) throw (MEDEXCEPTION);
1022 FIELD<double, FullInterlace>* getGaussPointsCoordinates() const throw (MEDEXCEPTION);
1025 \addtogroup FIELD_value
1030 This method makes it possible to have the field pointing to
1031 an existing value array. The ordering of the elements in the value array must
1032 conform to the MEDMEM ordering (I,K,J) : the outer loop is on the elements,
1033 the intermediate loop is on the Gauss points, the inner loop is on
1036 inline void setValue( T* value) throw (MEDEXCEPTION);
1037 inline void setRow( int i, T* value) throw (MEDEXCEPTION);
1038 inline void setColumn( int i, T* value) throw (MEDEXCEPTION);
1040 Sets the value of \f$ i^{th} \f$ element and \f$ j^{th}\f$ component with \a value.
1042 inline void setValueIJ(int i, int j, T value) throw (MEDEXCEPTION);
1043 /*! \if MEDMEM_ug @} \endif */
1044 inline void setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION);
1045 inline void setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION);
1046 inline void setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION);
1048 typedef void (*myFuncType)(const double *,T*);
1049 void fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION);
1050 typedef void (*myFuncType2)(const T *,T*);
1051 FIELD<T,INTERLACING_TAG> *execFunc(int nbOfComponents, myFuncType2 f) throw (MEDEXCEPTION);
1055 #include "MEDMEM_DriverFactory.hxx"
1059 template <class T,class INTERLACING_TAG> T FIELD<T, INTERLACING_TAG>::_scalarForPow=1;
1061 // --------------------
1062 // Implemented Methods
1063 // --------------------
1066 Constructor with no parameter, most of the attribut members are set to NULL.
1068 template <class T, class INTERLACING_TAG>
1069 FIELD<T, INTERLACING_TAG>::FIELD():FIELD_()
1071 MESSAGE_MED("Constructeur FIELD sans parametre");
1073 //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
1074 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE);
1075 FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
1077 //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
1078 ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE);
1079 FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
1081 _value = ( ArrayNoGauss * ) NULL;
1083 _mesh = ( MESH* ) NULL;
1087 \addtogroup FIELD_constructors FIELD<T> constructors
1092 Constructor that allocates the value array with the dimensions provided by
1093 \a NumberOfComponents and the dimension of \a Support. The value array is
1094 allocated but not initialized.
1095 This constructor does not allow the creation of fields with Gauss points.
1096 \param Support support on which the field lies
1097 \param NumberOfComponents number of components of the variable stored. For instance,
1098 it will be 3 for a (vx,vy,vz) vector.
1101 FIELD<double> field (support, 3);
1102 int nbelem = support->getNumberOfElements(MED_ALL_ELEMENTS);
1103 for (int i=1; i<=nbelem; i++)
1105 field->setValueIJ(i,j,0.0);
1108 template <class T, class INTERLACING_TAG>
1109 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
1110 const int NumberOfComponents) throw (MEDEXCEPTION) :
1111 FIELD_(Support, NumberOfComponents),_value(NULL)
1113 const char* LOC = "FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)";
1117 //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
1118 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
1119 FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
1121 //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
1122 ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
1123 FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
1127 // becarefull about the numbre of gauss point
1128 _numberOfValues = Support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
1130 #if defined(_DEBUG_) || defined(_DEBUG)
1131 catch (MEDEXCEPTION &ex)
1133 catch (MEDEXCEPTION )
1136 MESSAGE_MED("No value defined ! ("<<ex.what()<<")");
1138 MESSAGE_MED("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
1139 if ( _numberOfValues > 0 )
1141 if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
1143 const int * nbelgeo = Support->getNumberOfElements();
1144 vector<int> nbelgeoc( Support->getNumberOfTypes() + 1 );
1146 for ( int t = 1; t < (int)nbelgeoc.size(); ++t )
1147 nbelgeoc[t] = nbelgeoc[t-1] + nbelgeo[t-1];
1148 _value = new ArrayNoByType (_numberOfComponents,_numberOfValues,
1149 Support->getNumberOfTypes(), &nbelgeoc[0]);
1153 _value = new ArrayNoGauss (_numberOfComponents,_numberOfValues);
1157 _mesh = ( MESH* ) NULL;
1168 template <class T, class INTERLACING_TAG> void FIELD<T, INTERLACING_TAG>::init ()
1175 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::FIELD(const FIELD & m):
1178 MESSAGE_MED("Constructeur FIELD de recopie");
1180 // RECOPIE PROFONDE <> de l'operateur= Rmq from EF
1181 if (m._value != NULL)
1183 if ( m.getGaussPresence() )
1184 _value = new ArrayGauss( *(static_cast< ArrayGauss * > (m._value) ) ,false);
1186 _value = new ArrayNoGauss( *(static_cast< ArrayNoGauss * > (m._value)) ,false);
1189 _value = (ArrayNoGauss *) NULL;
1190 locMap::const_iterator it;
1191 for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
1192 _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
1193 new GAUSS_LOCALIZATION<INTERLACING_TAG>(
1194 *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
1197 _valueType = m._valueType;
1198 _interlacingType = m._interlacingType;
1199 //drivers = m._drivers;
1202 _mesh->addReference();
1210 template <class T, class INTERLACING_TAG>
1211 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(const FIELD &m)
1213 MESSAGE_MED("Appel de FIELD<T>::operator=") ;
1214 if ( this == &m) return *this;
1216 // copy values array
1217 FIELD_::operator=(m); // Driver are ignored & ?copie su pointeur de Support?
1219 _value = m._value; //PROBLEME RECOPIE DES POINTEURS PAS COHERENT AVEC LE
1220 //CONSTRUCTEUR PAR RECOPIE
1221 //CF :Commentaire dans MEDMEM_Array
1222 locMap::const_iterator it;
1223 for ( it = m._gaussModel.begin();it != m._gaussModel.end(); it++ )
1224 _gaussModel[static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second)->getType()]=
1225 new GAUSS_LOCALIZATION<INTERLACING_TAG>(
1226 *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ( (*it).second )
1229 _valueType = m._valueType;
1230 _interlacingType = m._interlacingType;
1234 _mesh->removeReference();
1237 _mesh->addReference();
1243 Initializes all the field values to \a value
1245 template <class T, class INTERLACING_TAG>
1246 FIELD<T, INTERLACING_TAG> & FIELD<T, INTERLACING_TAG>::operator=(T value)
1248 MESSAGE_MED("Appel de FIELD<T>::operator= T") ;
1249 int size=getNumberOfComponents()*getNumberOfValues();
1250 T* ptr= const_cast<T*>( getValue());
1251 for (int i=0; i< size; i++)
1257 /*!\addtogroup FIELD_algo
1262 Overload addition operator.
1263 This operation is authorized only for compatible fields that have the same support.
1264 The compatibility checking includes equality tests of the folowing data members:\n
1266 - _numberOfComponents
1269 - _MEDComponentsUnits.
1271 The data members of the returned field are initialized, based on the first field, except for the name,
1272 which is the combination of the two field's names and the operator.
1273 Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> \n
1274 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1275 When using python, this operator calls the copy constructor in any case.
1276 The user has to be aware that when using operator + in associatives expressions like
1277 <tt> a = b + c + d +e; </tt> \n
1278 no optimisation is performed : the evaluation of last expression requires the construction of
1281 template <class T, class INTERLACING_TAG>
1282 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator+(const FIELD & m) const
1284 const char* LOC = "FIELD<T>::operator+(const FIELD & m)";
1286 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1288 // Creation of the result - memory is allocated by FIELD constructor
1289 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1290 result->_operationInitialize(*this,m,"+"); // perform Atribute's initialization
1291 result->_add_in_place(*this,m); // perform addition
1297 /*! Overloaded Operator +=
1298 * Operations are directly performed in the first field's array.
1299 * This operation is authorized only for compatible fields that have the same support.
1301 template <class T, class INTERLACING_TAG>
1302 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator+=(const FIELD & m)
1304 const char* LOC = "FIELD<T>::operator+=(const FIELD & m)";
1306 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1308 const T* value1=m.getValue(); // get pointers to the values we are adding
1310 // get a non const pointer to the inside array of values and perform operation
1311 T * value=const_cast<T *> (getValue());
1312 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1313 const T* endV=value+size; // pointer to the end of value
1314 for(;value!=endV; value1++,value++)
1321 /*! Addition of fields. Static member function.
1322 * The function return a pointer to a new created field that holds the addition.
1323 * Data members are checked for compatibility and initialized.
1324 * The user is in charge of memory deallocation.
1326 template <class T, class INTERLACING_TAG>
1327 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::add(const FIELD& m, const FIELD& n)
1329 const char* LOC = "FIELD<T>::add(const FIELD & m, const FIELD& n)";
1331 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1333 // Creation of a new field
1334 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1335 m.getNumberOfComponents());
1336 result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1337 result->_add_in_place(m,n); // perform addition
1343 /*! Same as add method except that field check is deeper.
1345 template <class T, class INTERLACING_TAG>
1346 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::addDeep(const FIELD& m, const FIELD& n)
1348 const char* LOC = "FIELD<T>::addDeep(const FIELD & m, const FIELD& n)";
1350 FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1352 // Creation of a new field
1353 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1354 m.getNumberOfComponents());
1355 result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
1356 result->_add_in_place(m,n); // perform addition
1363 Overload substraction operator.
1364 This operation is authorized only for compatible fields that have the same support.
1365 The compatibility checking includes equality tests of the folowing data members:\n
1367 - _numberOfComponents
1370 - _MEDComponentsUnits.
1372 The data members of the returned field are initialized, based on the first field, except for the name,
1373 which is the combination of the two field's names and the operator.
1374 Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> \n
1375 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1376 When using python, this operator calls the copy constructor in any case.
1377 The user has to be aware that when using operator - in associatives expressions like
1378 <tt> a = b - c - d -e; </tt> \n
1379 no optimisation is performed : the evaluation of last expression requires the construction of
1382 template <class T, class INTERLACING_TAG>
1383 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator-(const FIELD & m) const
1385 const char* LOC = "FIELD<T>::operator-(const FIELD & m)";
1387 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1389 // Creation of the result - memory is allocated by FIELD constructor
1390 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1391 result->_operationInitialize(*this,m,"-"); // perform Atribute's initialization
1392 result->_sub_in_place(*this,m); // perform substracion
1398 template <class T, class INTERLACING_TAG>
1399 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator-() const
1401 const char* LOC = "FIELD<T>::operator-()";
1404 // Creation of the result - memory is allocated by FIELD constructor
1405 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1406 // Atribute's initialization
1407 result->setName("- "+getName());
1408 result->setComponentsNames(getComponentsNames());
1409 // not yet implemented setComponentType(getComponentType());
1410 result->setComponentsDescriptions(getComponentsDescriptions());
1411 result->setMEDComponentsUnits(getMEDComponentsUnits());
1412 result->setComponentsUnits(getComponentsUnits());
1413 result->setIterationNumber(getIterationNumber());
1414 result->setTime(getTime());
1415 result->setOrderNumber(getOrderNumber());
1417 const T* value1=getValue();
1418 // get a non const pointer to the inside array of values and perform operation
1419 T * value=const_cast<T *> (result->getValue());
1420 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1421 const T* endV=value+size; // pointer to the end of value
1423 for(;value!=endV; value1++,value++)
1424 *value = -(*value1);
1429 /*! Overloaded Operator -=
1430 * Operations are directly performed in the first field's array.
1431 * This operation is authorized only for compatible fields that have the same support.
1433 template <class T, class INTERLACING_TAG>
1434 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator-=(const FIELD & m)
1436 const char* LOC = "FIELD<T>::operator-=(const FIELD & m)";
1438 FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
1440 const T* value1=m.getValue();
1442 // get a non const pointer to the inside array of values and perform operation
1443 T * value=const_cast<T *> (getValue());
1444 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1445 const T* endV=value+size; // pointer to the end of value
1447 for(;value!=endV; value1++,value++)
1456 /*! Apply to a given field component the linear function x -> ax+b.
1457 * calculation is done "in place".
1459 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b, int icomp)
1461 // get a non const pointer to the inside array of values and perform operation in place
1462 T * value=const_cast<T *> (getValue());
1464 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1466 if (size>0) // for a negative size, there is nothing to do
1469 const T* lastvalue=value+size; // pointer to the end of value
1471 for(;value!=lastvalue; value+=getNumberOfComponents()) // apply linear transformation
1472 *value = a*(*value)+b;
1476 /*! Substraction of fields. Static member function.
1477 * The function return a pointer to a new created field that holds the substraction.
1478 * Data members are checked for compatibility and initialized.
1479 * The user is in charge of memory deallocation.
1481 template <class T, class INTERLACING_TAG>
1482 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::sub(const FIELD& m, const FIELD& n)
1484 const char* LOC = "FIELD<T>::sub(const FIELD & m, const FIELD& n)";
1486 FIELD_::_checkFieldCompatibility(m, n); // may throw exception
1488 // Creation of a new field
1489 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1490 m.getNumberOfComponents());
1491 result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1492 result->_sub_in_place(m,n); // perform substraction
1498 /*! Same as sub method except that field check is deeper.
1500 template <class T, class INTERLACING_TAG>
1501 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::subDeep(const FIELD& m, const FIELD& n)
1503 const char* LOC = "FIELD<T>::subDeep(const FIELD & m, const FIELD& n)";
1505 FIELD_::_deepCheckFieldCompatibility(m, n); // may throw exception
1507 // Creation of a new field
1508 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1509 m.getNumberOfComponents());
1510 result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
1511 result->_sub_in_place(m,n); // perform substraction
1518 Overload multiplication operator.
1519 This operation is authorized only for compatible fields that have the same support.
1520 The compatibility checking includes equality tests of the folowing data members:\n
1522 - _numberOfComponents
1525 - _MEDComponentsUnits.
1527 The data members of the returned field are initialized, based on the first field, except for the name,
1528 which is the combination of the two field's names and the operator.
1529 Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> \n
1530 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1531 When using python, this operator calls the copy constructor in any case.
1532 The user has to be aware that when using operator * in associatives expressions like
1533 <tt> a = b * c * d *e; </tt> \n
1534 no optimisation is performed : the evaluation of last expression requires the construction of
1537 template <class T, class INTERLACING_TAG>
1538 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator*(const FIELD & m) const
1540 const char* LOC = "FIELD<T>::operator*(const FIELD & m)";
1542 FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1544 // Creation of the result - memory is allocated by FIELD constructor
1545 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1546 result->_operationInitialize(*this,m,"*"); // perform Atribute's initialization
1547 result->_mul_in_place(*this,m); // perform multiplication
1553 /*! Overloaded Operator *=
1554 * Operations are directly performed in the first field's array.
1555 * This operation is authorized only for compatible fields that have the same support.
1557 template <class T, class INTERLACING_TAG>
1558 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator*=(const FIELD & m)
1560 const char* LOC = "FIELD<T>::operator*=(const FIELD & m)";
1562 FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1564 const T* value1=m.getValue();
1566 // get a non const pointer to the inside array of values and perform operation
1567 T * value=const_cast<T *> (getValue());
1568 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1569 const T* endV=value+size; // pointer to the end of value
1571 for(;value!=endV; value1++,value++)
1579 /*! Multiplication of fields. Static member function.
1580 * The function return a pointer to a new created field that holds the multiplication.
1581 * Data members are checked for compatibility and initialized.
1582 * The user is in charge of memory deallocation.
1584 template <class T, class INTERLACING_TAG>
1585 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mul(const FIELD& m, const FIELD& n)
1587 const char* LOC = "FIELD<T>::mul(const FIELD & m, const FIELD& n)";
1589 FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1591 // Creation of a new field
1592 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1593 m.getNumberOfComponents());
1594 result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1595 result->_mul_in_place(m,n); // perform multiplication
1601 /*! Same as mul method except that field check is deeper.
1603 template <class T, class INTERLACING_TAG>
1604 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::mulDeep(const FIELD& m, const FIELD& n)
1606 const char* LOC = "FIELD<T>::mulDeep(const FIELD & m, const FIELD& n)";
1608 FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1610 // Creation of a new field
1611 FIELD<T, INTERLACING_TAG>* result = new FIELD<T,INTERLACING_TAG>(m.getSupport(),
1612 m.getNumberOfComponents());
1613 result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
1614 result->_mul_in_place(m,n); // perform multiplication
1621 Overload division operator.
1622 This operation is authorized only for compatible fields that have the same support.
1623 The compatibility checking includes equality tests of the folowing data members:\n
1625 - _numberOfComponents
1628 - _MEDComponentsUnits.
1630 The data members of the returned field are initialized, based on the first field, except for the name,
1631 which is the combination of the two field's names and the operator.
1632 Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> \n
1633 In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
1634 When using python, this operator calls the copy constructor in any case.
1635 The user has to be aware that when using operator / in associatives expressions like
1636 <tt> a = b / c / d /e; </tt> \n
1637 no optimisation is performed : the evaluation of last expression requires the construction of
1640 template <class T, class INTERLACING_TAG>
1641 FIELD<T, INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::operator/(const FIELD & m) const
1643 const char* LOC = "FIELD<T>::operator/(const FIELD & m)";
1645 FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1647 // Creation of the result - memory is allocated by FIELD constructor
1648 FIELD<T, INTERLACING_TAG> *result=new FIELD<T, INTERLACING_TAG>(this->getSupport(),this->getNumberOfComponents());
1651 result->_operationInitialize(*this,m,"/"); // perform Atribute's initialization
1652 result->_div_in_place(*this,m); // perform division
1654 catch(MEDEXCEPTION& e)
1656 result->removeReference();
1665 /*! Overloaded Operator /=
1666 * Operations are directly performed in the first field's array.
1667 * This operation is authorized only for compatible fields that have the same support.
1669 template <class T, class INTERLACING_TAG>
1670 FIELD<T, INTERLACING_TAG>& FIELD<T, INTERLACING_TAG>::operator/=(const FIELD & m)
1672 const char* LOC = "FIELD<T>::operator/=(const FIELD & m)";
1674 FIELD_::_checkFieldCompatibility(*this, m, false); // may throw exception
1676 const T* value1=m.getValue(); // get pointers to the values we are adding
1678 // get a non const pointer to the inside array of values and perform operation
1679 T * value=const_cast<T *> (getValue());
1680 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
1681 const T* endV=value+size; // pointer to the end of value
1683 for(;value!=endV; value1++,value++)
1691 /*! Division of fields. Static member function.
1692 * The function return a pointer to a new created field that holds the division.
1693 * Data members are checked for compatibility and initialized.
1694 * The user is in charge of memory deallocation.
1696 template <class T, class INTERLACING_TAG>
1697 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::div(const FIELD& m, const FIELD& n)
1699 const char* LOC = "FIELD<T>::div(const FIELD & m, const FIELD& n)";
1701 FIELD_::_checkFieldCompatibility(m, n, false); // may throw exception
1703 // Creation of a new field
1704 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1705 m.getNumberOfComponents());
1708 result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1709 result->_div_in_place(m,n); // perform division
1711 catch(MEDEXCEPTION& e)
1713 result->removeReference();
1720 /*! Same as div method except that field check is deeper.
1722 template <class T,class INTERLACING_TAG>
1723 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::divDeep(const FIELD& m, const FIELD& n)
1725 const char* LOC = "FIELD<T>::divDeep(const FIELD & m, const FIELD& n)";
1727 FIELD_::_deepCheckFieldCompatibility(m, n, false); // may throw exception
1729 // Creation of a new field
1730 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),
1731 m.getNumberOfComponents());
1734 result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
1735 result->_div_in_place(m,n); // perform division
1737 catch(MEDEXCEPTION& e)
1739 result->removeReference();
1753 This internal method initialize the members of a new field created to hold the result of the operation Op .
1754 Initialization is based on the first field, except for the name, which is the combination of the two field's names
1758 template <class T, class INTERLACING_TAG>
1759 void FIELD<T, INTERLACING_TAG>::_operationInitialize(const FIELD& m,const FIELD& n, const char* Op)
1761 MESSAGE_MED("Appel methode interne " << Op);
1763 // Atribute's initialization (copy of the first field's attributes)
1764 // Other data members (_support, _numberOfValues) are initialized in the field's constr.
1765 setName(m.getName()+" "+Op+" "+n.getName());
1766 setComponentsNames(m.getComponentsNames());
1767 setComponentsDescriptions(m.getComponentsDescriptions());
1768 setMEDComponentsUnits(m.getMEDComponentsUnits());
1770 // The following data member may differ from field m to n.
1771 // The initialization is done based on the first field.
1773 setComponentsUnits(m.getComponentsUnits());
1775 setIterationNumber(m.getIterationNumber());
1776 setTime(m.getTime());
1777 setOrderNumber(m.getOrderNumber());
1783 Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
1784 This method is applied to a just created field with medModeSwitch mode.
1785 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1789 template <class T, class INTERLACING_TAG>
1790 void FIELD<T, INTERLACING_TAG>::_add_in_place(const FIELD& m,const FIELD& n)
1792 // get pointers to the values we are adding
1793 const T* value1=m.getValue();
1794 const T* value2=n.getValue();
1795 // get a non const pointer to the inside array of values and perform operation
1796 T * value=const_cast<T *> (getValue());
1798 const int size=getNumberOfValues()*getNumberOfComponents();
1800 const T* endV1=value1+size;
1801 for(;value1!=endV1; value1++,value2++,value++)
1802 *value=(*value1)+(*value2);
1807 Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
1808 This method is applied to a just created field with medModeSwitch mode.
1809 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1813 template <class T, class INTERLACING_TAG>
1814 void FIELD<T, INTERLACING_TAG>::_sub_in_place(const FIELD& m,const FIELD& n)
1816 // get pointers to the values we are adding
1817 const T* value1=m.getValue();
1818 const T* value2=n.getValue();
1819 // get a non const pointer to the inside array of values and perform operation
1820 T * value=const_cast<T *> (getValue());
1822 const int size=getNumberOfValues()*getNumberOfComponents();
1824 const T* endV1=value1+size;
1825 for(;value1!=endV1; value1++,value2++,value++)
1826 *value=(*value1)-(*value2);
1831 Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
1832 This method is applied to a just created field with medModeSwitch mode.
1833 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1837 template <class T, class INTERLACING_TAG>
1838 void FIELD<T, INTERLACING_TAG>::_mul_in_place(const FIELD& m,const FIELD& n)
1840 // get pointers to the values we are adding
1841 const T* value1=m.getValue();
1842 const T* value2=n.getValue();
1843 // get a non const pointer to the inside array of values and perform operation
1844 T * value=const_cast<T *> (getValue());
1846 const int size=getNumberOfValues()*getNumberOfComponents();
1848 const T* endV1=value1+size;
1849 for(;value1!=endV1; value1++,value2++,value++)
1850 *value=(*value1)*(*value2);
1855 Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
1856 This method is applied to a just created field with medModeSwitch mode.
1857 For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
1861 template <class T, class INTERLACING_TAG>
1862 void FIELD<T, INTERLACING_TAG>::_div_in_place(const FIELD& m,const FIELD& n) throw (MEDEXCEPTION)
1864 // get pointers to the values we are adding
1865 const T* value1=m.getValue();
1866 const T* value2=n.getValue();
1867 // get a non const pointer to the inside array of values and perform operation
1868 T * value=const_cast<T *> (getValue());
1870 const int size=getNumberOfValues()*getNumberOfComponents();
1872 const T* endV1=value1+size;
1873 for(;value1!=endV1; value1++,value2++,value++){
1874 if ( *value2 == 0 ) { // FAIRE PLUTOT UN TRY CATCH Rmq from EF
1876 diagnosis="FIELD<T,INTERLACING_TAG>::_div_in_place(...) : Divide by zero !";
1877 throw MEDEXCEPTION(diagnosis.c_str());
1879 *value=(*value1)/(*value2);
1884 \addtogroup FIELD_algo
1888 /*! Return maximum of all absolute values contained in the array (all elements and all components are browsed).
1890 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::normMax() const throw (MEDEXCEPTION)
1892 const T* value=getValue(); // get pointer to the values
1893 const int size=getNumberOfValues()*getNumberOfComponents();
1894 if (size <= 0) // Size of array has to be strictly positive
1897 diagnosis="FIELD<T,INTERLACIN_TAG>::normMax() : cannot compute the norm of "+getName()+
1898 " : it size is non positive!";
1899 throw MEDEXCEPTION(diagnosis.c_str());
1901 const T* lastvalue=value+size; // get pointer just after last value
1902 const T* pMax=value; // pointer to the max value
1903 const T* pMin=value; // pointer to the min value
1905 // get pointers to the max & min value of array
1906 while ( ++value != lastvalue )
1908 if ( *pMin > *value )
1910 if ( *pMax < *value )
1914 T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
1915 T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
1917 return Max>Min ? double(Max) : double(Min);
1920 /*! Return Euclidian norm for all elements of the array.
1922 template <class T, class INTERLACIN_TAG> double FIELD<T, INTERLACIN_TAG>::norm2() const throw (MEDEXCEPTION)
1924 const T* value=this->getValue(); // get const pointer to the values
1925 const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
1926 if (size <= 0) // Size of array has to be strictly positive
1929 diagnosis="FIELD<T,INTERLACIN_TAG>::norm2() : cannot compute the norm of "+getName()+
1930 " : it size is non positive!";
1931 throw MEDEXCEPTION(diagnosis.c_str());
1933 const T* lastvalue=value+size; // point just after last value
1935 T result((T)0); // init
1936 for( ; value!=lastvalue ; ++value)
1937 result += (*value) * (*value);
1939 return std::sqrt(double(result));
1943 //------------- TDG and BS addings
1945 /*! Return Extrema of field
1947 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::getMinMax(T &vmin, T &vmax) throw (MEDEXCEPTION)
1949 const T* value=getValue(); // get pointer to the values
1950 const int size=getNumberOfValues()*getNumberOfComponents();
1951 const T* lastvalue=value+size; // point just after last value
1953 if (size <= 0){ // Size of array has to be strictly positive
1956 diagnosis="FIELD<T,INTERLACIN_TAG>::getMinMax() : cannot compute the extremums of "+getName()+
1957 " : its size is non positive!";
1958 throw MEDEXCEPTION(diagnosis.c_str());
1962 vmax=MinMax<T>::getMin(); // init a max value
1963 vmin=MinMax<T>::getMax(); // init a min value
1965 for( ; value!=lastvalue ; ++value){
1966 if ( vmin > *value )
1968 if ( vmax < *value )
1982 /*! Return Histogram of field
1984 template <class T, class INTERLACIN_TAG> vector<int> FIELD<T, INTERLACIN_TAG>::getHistogram(int &nbint) throw (MEDEXCEPTION)
1986 const T* value=getValue(); // get pointer to the values
1987 const int size=getNumberOfValues()*getNumberOfComponents();
1988 const T* lastvalue=value+size; // point just after last value
1990 if (size <= 0){ // Size of array has to be strictly positive
1993 diagnosis="FIELD<T,INTERLACIN_TAG>::getHistogram() : cannot compute the histogram of "+getName()+
1994 " : it size is non positive!";
1995 throw MEDEXCEPTION(diagnosis.c_str());
1998 vector<int> Histogram(nbint) ;
2002 for( j=0 ; j!=nbint ; j++) Histogram[j]=0 ;
2004 getMinMax(vmin,vmax);
2005 for( ; value!=lastvalue ; ++value){
2006 if(*value==vmax) j = nbint-1;
2007 else j = (int)(((double)nbint * (*value-vmin))/(vmax-vmin));
2015 /*! Return vectorial gradient field
2017 template <class T, class INTERLACIN_TAG>
2018 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildGradient() const throw (MEDEXCEPTION)
2020 const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildGradient() : ";
2023 // space dimension of input mesh
2024 int spaceDim = getSupport()->getMesh()->getSpaceDimension();
2025 double *x = new double[spaceDim];
2027 FIELD<double, FullInterlace>* Gradient =
2028 new FIELD<double, FullInterlace>(getSupport(),spaceDim);
2030 string name("gradient of ");
2032 Gradient->setName(name);
2033 string descr("gradient of ");
2034 descr += getDescription();
2035 Gradient->setDescription(descr);
2037 if( _numberOfComponents > 1 ){
2040 throw MEDEXCEPTION("gradient calculation only on scalar field");
2043 for(int i=1;i<=spaceDim;i++){
2044 string nameC("gradient of ");
2046 Gradient->setComponentName(i,nameC);
2047 Gradient->setComponentDescription(i,"gradient");
2048 string MEDComponentUnit = getMEDComponentUnit(1)+getSupport()->getMesh()->getCoordinatesUnits()[i-1];
2049 Gradient->setMEDComponentUnit(i,MEDComponentUnit);
2052 Gradient->setIterationNumber(getIterationNumber());
2053 Gradient->setOrderNumber(getOrderNumber());
2054 Gradient->setTime(getTime());
2056 // typ of entity on what is field
2057 MED_EN::medEntityMesh typ = getSupport()->getEntity();
2063 const double *coord;
2071 // read connectivity array to have the list of nodes contained by an element
2072 C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,typ,MED_ALL_ELEMENTS);
2073 iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,typ);
2074 // calculate reverse connectivity to have the list of elements which contains node i
2075 revC = getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,typ);
2076 indC = getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,typ);
2077 // coordinates of each node
2078 coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
2079 // number of elements
2080 NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2081 // barycenter field of elements
2082 FIELD<double, FullInterlace>* barycenter = getSupport()->getMesh()->getBarycenter(getSupport());
2084 // calculate gradient vector for each element i
2085 for (int i = 1; i < NumberOf + 1; i++) {
2087 // listElements contains elements which contains a node of element i
2088 set <int> listElements;
2089 set <int>::iterator elemIt;
2090 listElements.clear();
2092 // for each node j of element i
2093 for (int ij = iC[i-1]; ij < iC[i]; ij++) {
2095 for (int k = indC[j-1]; k < indC[j]; k++) {
2096 // c element contains node j
2098 // we put the elements in set
2100 listElements.insert(c);
2103 // coordinates of barycentre of element i in space of dimension spaceDim
2104 for (int j = 0; j < spaceDim; j++)
2105 x[j] = barycenter->getValueIJ(i,j+1);
2107 for (int j = 0; j < spaceDim; j++) {
2108 // value of field of element i
2109 double val = getValueIJ(i,1);
2111 // calculate gradient for each neighbor element
2112 for (elemIt = listElements.begin(); elemIt != listElements.end(); elemIt++) {
2115 for (int l = 0; l < spaceDim; l++) {
2116 // coordinate of barycenter of element elem
2117 double xx = barycenter->getValueIJ(elem, l+1);
2118 d2 += (x[l]-xx) * (x[l]-xx);
2120 grad += (barycenter->getValueIJ(elem,j+1)-x[j])*(getValueIJ(elem,1)-val)/sqrt(d2);
2122 if (listElements.size() != 0) grad /= listElements.size();
2123 Gradient->setValueIJ(i,j+1,grad);
2130 // read connectivity array to have the list of nodes contained by an element
2131 C = getSupport()->getMesh()->getConnectivity(MED_FULL_INTERLACE,MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
2132 iC = getSupport()->getMesh()->getConnectivityIndex(MED_NODAL,MED_CELL);
2133 // calculate reverse connectivity to have the list of elements which contains node i
2134 revC=getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
2135 indC=getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
2136 // coordinates of each node
2137 coord = getSupport()->getMesh()->getCoordinates(MED_FULL_INTERLACE);
2139 // calculate gradient for each node
2140 NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2141 for (int i=1; i<NumberOf+1; i++){
2142 // listNodes contains nodes neigbor of node i
2143 set <int> listNodes;
2144 set <int>::iterator nodeIt ;
2146 for(int j=indC[i-1];j<indC[i];j++){
2147 // c element contains node i
2149 // we put the nodes of c element in set
2150 for(int k=iC[c-1];k<iC[c];k++)
2152 listNodes.insert(C[k-1]);
2154 // coordinates of node i in space of dimension spaceDim
2155 for(int j=0;j<spaceDim;j++)
2156 x[j] = coord[(i-1)*spaceDim+j];
2158 for(int j=0;j<spaceDim;j++){
2160 double val = getValueIJ(i,1);
2162 // calculate gradient for each neighbor node
2163 for(nodeIt=listNodes.begin();nodeIt!=listNodes.end();nodeIt++){
2166 for(int l=0;l<spaceDim;l++){
2167 double xx = coord[(node-1)*spaceDim+l];
2168 d2 += (x[l]-xx) * (x[l]-xx);
2170 grad += (coord[(node-1)*spaceDim+j]-x[j])*(getValueIJ(node,1)-val)/sqrt(d2);
2172 if(listNodes.size() != 0) grad /= listNodes.size();
2173 Gradient->setValueIJ(i,j+1,grad);
2177 case MED_ALL_ENTITIES:
2180 throw MEDEXCEPTION("gradient calculation not yet implemented on all elements");
2190 /*! Return scalar norm2 field
2192 template <class T, class INTERLACIN_TAG>
2193 FIELD<double, FullInterlace>* FIELD<T, INTERLACIN_TAG>::buildNorm2Field() const throw (MEDEXCEPTION)
2195 const char * LOC = "FIELD<T, INTERLACIN_TAG>::buildNorm2Field() : ";
2198 FIELD<double, FullInterlace>* Norm2Field =
2199 new FIELD<double, FullInterlace>(getSupport(),1);
2201 string name("norm2 of ");
2203 Norm2Field->setName(name);
2204 string descr("norm2 of ");
2205 descr += getDescription();
2206 Norm2Field->setDescription(descr);
2208 string nameC("norm2 of ");
2210 Norm2Field->setComponentName(1,nameC);
2211 Norm2Field->setComponentDescription(1,"norm2");
2212 string MEDComponentUnit = getMEDComponentUnit(1);
2213 Norm2Field->setMEDComponentUnit(1,MEDComponentUnit);
2215 Norm2Field->setIterationNumber(getIterationNumber());
2216 Norm2Field->setOrderNumber(getOrderNumber());
2217 Norm2Field->setTime(getTime());
2219 // calculate nom2 for each element
2220 int NumberOf = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
2221 for (int i=1; i<NumberOf+1; i++){
2223 for(int j=1;j<=getNumberOfComponents();j++)
2224 norm2 += getValueIJ(i,j)*getValueIJ(i,j);
2225 Norm2Field->setValueIJ(i,1,sqrt(norm2));
2233 /*! Apply to each (scalar) field component the template parameter T_function,
2234 * which is a pointer to function.
2235 * Since the pointer is known at compile time, the function is inlined into the inner loop!
2236 * calculation is done "in place".
2239 * \code myField.applyFunc<std::sqrt>(); // apply sqare root function \endcode
2240 * \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
2242 template <class T, class INTERLACIN_TAG> template <T T_function(T)>
2243 void FIELD<T, INTERLACIN_TAG>::applyFunc()
2245 // get a non const pointer to the inside array of values and perform operation
2246 T * value=const_cast<T *> (getValue());
2247 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
2249 if (size>0) // for a negative size, there is nothing to do
2251 const T* lastvalue=value+size; // pointer to the end of value
2252 for(;value!=lastvalue; ++value) // apply linear transformation
2253 *value = T_function(*value);
2257 template <class T, class INTERLACIN_TAG> T FIELD<T, INTERLACIN_TAG>::pow(T x)
2259 return (T)::pow((double)x,FIELD<T, INTERLACIN_TAG>::_scalarForPow);
2262 /*! Apply to each (scalar) field component the math function pow.
2263 * calculation is done "in place".
2266 * \code myField.applyFunc<std::sqrt>(); // apply sqare root function \endcode
2267 * \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
2269 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyPow(T scalar)
2271 FIELD<T, INTERLACIN_TAG>::_scalarForPow=scalar;
2272 applyFunc<FIELD<T, INTERLACIN_TAG>::pow>();
2275 /*! Apply to each (scalar) field component the linear function x -> ax+b.
2276 * calculation is done "in place".
2278 template <class T, class INTERLACIN_TAG> void FIELD<T, INTERLACIN_TAG>::applyLin(T a, T b)
2280 // get a non const pointer to the inside array of values and perform operation in place
2281 T * value=const_cast<T *> (getValue());
2282 const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
2284 if (size>0) // for a negative size, there is nothing to do
2286 const T* lastvalue=value+size; // pointer to the end of value
2287 for(;value!=lastvalue; ++value) // apply linear transformation
2288 *value = a*(*value)+b;
2294 * Return a pointer to a new field that holds the scalar product. Static member function.
2295 * This operation is authorized only for compatible fields that have the same support.
2296 * The compatibility checking includes equality tests of the folowing data members:\n
2298 * - _numberOfComponents
2300 * - _componentsTypes
2301 * - _MEDComponentsUnits.
2302 * Data members are initialized.
2303 * The new field point to the same support and has one component.
2304 * Each value of it is the scalar product of the two argument's fields.
2305 * The user is in charge of memory deallocation.
2307 template <class T, class INTERLACING_TAG>
2308 FIELD<T, INTERLACING_TAG>*
2309 FIELD<T, INTERLACING_TAG>::scalarProduct(const FIELD & m, const FIELD & n, bool deepCheck)
2312 FIELD_::_checkFieldCompatibility( m, n, false); // may throw exception
2314 FIELD_::_deepCheckFieldCompatibility(m, n, false);
2316 // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
2317 // result type imply INTERLACING_TAG=FullInterlace for m & n
2319 const int numberOfElements=m.getNumberOfValues(); // strictly positive
2320 const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
2322 // Creation & init of a the result field on the same support, with one component
2323 // You have to be careful about the interlacing mode, because in the computation step,
2324 // it seems to assume the the interlacing mode is the FullInterlacing
2326 FIELD<T, INTERLACING_TAG>* result = new FIELD<T, INTERLACING_TAG>(m.getSupport(),1);
2327 result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
2328 result->setIterationNumber(m.getIterationNumber());
2329 result->setTime(m.getTime());
2330 result->setOrderNumber(m.getOrderNumber());
2332 const T* value1=m.getValue(); // get const pointer to the values
2333 const T* value2=n.getValue(); // get const pointer to the values
2334 // get a non const pointer to the inside array of values and perform operation
2335 T * value=const_cast<T *> (result->getValue());
2337 const T* lastvalue=value+numberOfElements; // pointing just after last value of result
2338 for ( ; value!=lastvalue ; ++value ) // loop on all elements
2340 *value=(T)0; // initialize value
2341 const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
2342 for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
2343 *value += (*value1) * (*value2);
2348 /*! Return L2 Norm of the field's component.
2349 * Cannot be applied to a field with a support on nodes.
2350 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2351 * For the nodal field, p_field_volume must be for all cells even if the field is partial.
2353 template <class T, class INTERLACING_TAG>
2354 double FIELD<T, INTERLACING_TAG>::normL2(int component,
2355 const FIELD<double, FullInterlace> * p_field_volume) const
2357 _checkNormCompatibility(p_field_volume, /*nodalAllowed=*/true); // may throw exception
2358 if ( component<1 || component>getNumberOfComponents() )
2359 throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
2361 const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2362 if(!p_field_volume) // if the user don't supply the volume
2363 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
2365 p_field_size->addReference();
2366 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2367 const double* vol=p_field_size->getValue();
2368 // Il n'est vraiment pas optimal de mixer des champs dans des modes d'entrelacement
2369 // different juste pour le calcul
2372 double integrale=0.0;
2375 if ( getSupport()->getEntity() == MED_NODE ) // issue 20120: [CEA 206] normL2 on NODE field
2377 //Most frequently the FIELD is on the whole mesh and
2378 // there is no need in optimizing iterations from supporting nodes-> back to cells,
2379 // so we iterate just on all cells
2380 const MESH * mesh = getSupport()->getMesh()->convertInMESH();
2381 const int nbCells = mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
2382 const int *C = mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
2383 const int *iC = mesh->getConnectivityIndex(MED_NODAL,MED_CELL);
2384 for (int i = 0; i < nbCells; ++i, ++vol) {
2385 // calculate integral on current element as average summ of values on all it's nodes
2386 double curCellValue = 0;
2387 try { // we expect exception with partial fields for nodes w/o values
2388 for (int ij = iC[i]; ij < iC[i+1]; ij++) {
2390 curCellValue += getValueIJ( node, component );
2393 catch ( MEDEXCEPTION ) {
2396 int nbNodes = iC[i+1]-iC[i];
2397 curCellValue /= nbNodes;
2398 integrale += (curCellValue * curCellValue) * std::abs(*vol);
2399 totVol+=std::abs(*vol);
2401 mesh->removeReference();
2405 if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2406 const T* value = getValue();
2407 value = value + (component-1) * getNumberOfValues();
2408 const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
2409 for (; value!=lastvalue ; ++value ,++vol) {
2410 integrale += double((*value) * (*value)) * std::abs(*vol);
2411 totVol+=std::abs(*vol);
2414 else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2415 ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2416 for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2417 integrale += anArray->getIJ(i,component) * anArray->getIJ(i,component) * std::abs(*vol);
2418 totVol+=std::abs(*vol);
2421 else { // FULL_INTERLACE
2422 ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2423 for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2424 integrale += anArray->getIJ(i,component) * anArray->getIJ(i,component) * std::abs(*vol);
2425 totVol+=std::abs(*vol);
2431 p_field_size->removeReference(); // delete temporary volume field
2434 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2435 return integrale/totVol;
2438 /*! Return L2 Norm of the field.
2439 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2440 * For the nodal field, p_field_volume must be for all cells even if the field is partial.
2442 template <class T, class INTERLACING_TAG>
2443 double FIELD<T, INTERLACING_TAG>::normL2(const FIELD<double, FullInterlace> * p_field_volume) const
2445 _checkNormCompatibility(p_field_volume, /*nodalAllowed=*/true); // may throw exception
2446 const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2447 if(!p_field_volume) // if the user don't supply the volume
2448 p_field_size=_getFieldSize(); // we calculate the volume
2450 p_field_size->addReference();
2451 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2452 const double* vol=p_field_size->getValue();
2453 const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
2456 double integrale=0.0;
2459 if ( getSupport()->getEntity() == MED_NODE ) // issue 20120: [CEA 206] normL2 on NODE field
2461 //Most frequently the FIELD is on the whole mesh and
2462 // there is no need in optimizing iterations from supporting nodes-> back to cells,
2463 // so we iterate just on all cells
2464 const MESH * mesh = getSupport()->getMesh()->convertInMESH();
2465 const int nbCells = mesh->getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS);
2466 const int *C = mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS);
2467 const int *iC = mesh->getConnectivityIndex(MED_NODAL,MED_CELL);
2468 int nbComp = getNumberOfComponents();
2469 for (int i = 0; i < nbCells; ++i, ++vol) {
2470 // calculate integral on current element as average summ of values on all it's nodes
2471 int nbNodes = iC[i+1]-iC[i];
2472 vector< double > curCellValue( nbComp, 0 );
2473 try { // we expect exception with partial fields for nodes w/o values
2474 for (int ij = iC[i]; ij < iC[i+1]; ij++) {
2476 for ( int j = 0; j < nbComp; ++j )
2477 curCellValue[ j ] += getValueIJ( node, j+1 ) / nbNodes;
2480 catch ( MEDEXCEPTION ) {
2484 for ( int j = 0; j < nbComp; ++j ) {
2485 integrale += (curCellValue[j] * curCellValue[j]) * std::abs(*vol);
2487 totVol+=std::abs(*vol);
2489 mesh->removeReference();
2490 if ( nbCells > 0 && totVol == 0.)
2491 throw MEDEXCEPTION("can't compute sobolev norm : "
2492 "none of elements has values on all it's nodes");
2496 const double* p_vol=vol;
2497 for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2498 totVol+=std::abs(*p_vol);
2500 if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2501 const T* value = getValue();
2502 for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
2503 for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
2504 integrale += (*value) * (*value) * std::abs(*p_vol);
2508 else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2509 ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2510 for (int j=1; j<=anArray->getDim(); j++) {
2512 for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2513 integrale += anArray->getIJ(i,j) * anArray->getIJ(i,j) * std::abs(*p_vol);
2517 else { // FULL_INTERLACE
2518 ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2519 for (int j=1; j<=anArray->getDim(); j++) {
2521 for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2522 integrale += anArray->getIJ(i,j) * anArray->getIJ(i,j) * std::abs(*p_vol);
2528 p_field_size->removeReference(); // delete temporary volume field
2531 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2532 return integrale/totVol;
2535 /*! Return L1 Norm of the field's component.
2536 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2538 template <class T, class INTERLACING_TAG>
2539 double FIELD<T, INTERLACING_TAG>::normL1(int component,
2540 const FIELD<double, FullInterlace > * p_field_volume) const
2542 _checkNormCompatibility(p_field_volume); // may throw exception
2543 if ( component<1 || component>getNumberOfComponents() )
2544 throw MEDEXCEPTION(STRING("FIELD<T,INTERLACING_TAG>::normL1() : The component argument should be between 1 and the number of components"));
2546 const FIELD<double,FullInterlace> * p_field_size=p_field_volume;
2547 if(!p_field_volume) // if the user don't supply the volume
2548 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implÃÃÅ mentation dans mesh]
2550 p_field_size->addReference();
2551 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2552 const double* vol = p_field_size->getValue();
2554 double integrale=0.0;
2557 if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2558 const T* value = getValue();
2559 const T* lastvalue = value + getNumberOfValues(); // pointing just after the end of column
2560 for (; value!=lastvalue ; ++value ,++vol) {
2561 integrale += std::abs( *value * *vol );
2562 totVol+=std::abs(*vol);
2565 else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2566 ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2567 for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2568 integrale += std::abs( anArray->getIJ(i,component) * (*vol));
2569 totVol+=std::abs(*vol);
2572 else { // FULL_INTERLACE
2573 ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2574 for (int i=1; i <= anArray->getNbElem() ; i++, ++vol ) {
2575 integrale += std::abs( anArray->getIJ(i,component) * *vol);
2576 totVol+=std::abs(*vol);
2581 p_field_size->removeReference(); // delete temporary volume field
2583 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2584 return integrale/totVol;
2587 /*! Return L1 Norm of the field.
2588 * Cannot be applied to a field with a support on nodes.
2589 * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
2591 template <class T, class INTERLACING_TAG>
2592 double FIELD<T, INTERLACING_TAG>::normL1(const FIELD<double, FullInterlace> * p_field_volume) const
2594 _checkNormCompatibility(p_field_volume); // may throw exception
2595 const FIELD<double, FullInterlace> * p_field_size=p_field_volume;
2596 if(!p_field_volume) // if the user don't supply the volume
2597 p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
2599 p_field_size->addReference();
2600 // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
2601 const double* vol = p_field_size->getValue();
2602 const double* lastvol = vol+getNumberOfValues(); // pointing just after the end of vol
2604 double integrale=0.0;
2606 const double* p_vol=vol;
2607 for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
2608 totVol+=std::abs(*p_vol);
2610 if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE ) {
2611 const T* value = getValue();
2612 for (int i=1; i<=getNumberOfComponents(); ++i) { // compute integral on all components
2613 for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol) {
2614 integrale += std::abs( *value * *p_vol );
2618 else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
2619 ArrayNoByType* anArray = dynamic_cast< ArrayNoByType * > ( getArrayNoGauss() );
2620 for (int j=1; j<=anArray->getDim(); j++) {
2622 for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2623 integrale += std::abs( anArray->getIJ(i,j) * *p_vol );
2627 else { // FULL_INTERLACE
2628 ArrayFull* anArray = dynamic_cast< ArrayFull * > ( getArrayNoGauss() );
2629 for (int j=1; j<=anArray->getDim(); j++) {
2631 for (p_vol=vol; i<=anArray->getNbElem() || p_vol!=lastvol; i++, ++p_vol ) {
2632 integrale += std::abs( anArray->getIJ(i,j) * *p_vol );
2637 p_field_size->removeReference();
2639 throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
2640 return integrale/totVol;
2644 * \brief Return integral of the field.
2645 * \param subSupport - optional part of a field to consider.
2646 * \retval double - value of integral
2648 template <class T, class INTERLACING_TAG>
2649 double FIELD<T, INTERLACING_TAG>::integral(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
2651 const char* LOC = "FIELD<>::integral(subSupport): ";
2653 double integrale = 0;
2655 if (!subSupport ) subSupport = _support;
2657 // check feasibility
2658 if ( getGaussPresence() )
2659 throw MEDEXCEPTION(STRING(LOC)<<"Gauss numbers greater than one are not yet implemented!");
2660 if ( subSupport->getEntity() != _support->getEntity())
2661 throw MEDEXCEPTION(STRING(LOC)<<"Different support entity of this field and subSupport");
2662 if ( subSupport->getEntity() == MED_EN::MED_NODE )
2663 throw MEDEXCEPTION(STRING(LOC)<<"Integral of nodal field not yet supported");
2666 const int nbElems = subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2667 const bool subOnAll = ( subSupport->isOnAllElements() );
2668 const bool myOnAll = ( _support->isOnAllElements() );
2669 const int* subNums = !subOnAll ? subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS) : 0;
2670 const int* myNums = !myOnAll ? _support->getNumber(MED_EN::MED_ALL_ELEMENTS) : 0;
2671 if ( !subOnAll && !subNums )
2672 throw MEDEXCEPTION(STRING(LOC)<<"Invalid support: no element numbers");
2673 if ( !myOnAll && !myNums )
2674 throw MEDEXCEPTION(STRING(LOC)<<"Invalid field support: no element numbers");
2675 if ( subOnAll && !myOnAll )
2676 return integral(NULL);
2678 // get size of elements
2679 const FIELD<double, FullInterlace> * cellSize=_getFieldSize(subSupport);
2680 const double* size = cellSize->getValue();
2681 const double* lastSize = size + nbElems; // pointing just after the end of size
2683 const T* value = getValue();
2685 // calculate integrale
2686 if ( (subOnAll && _support->isOnAllElements()) || subSupport == _support )
2688 const double* p_vol;
2689 if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2691 for (int j=1; j<=getNumberOfComponents(); ++j)
2692 for ( p_vol=size; p_vol != lastSize; ++value ,++p_vol)
2693 integrale += std::abs( *value * *p_vol );
2695 else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
2697 typename ArrayNoByType::InterlacingPolicy* indexer =
2698 dynamic_cast< typename ArrayNoByType::InterlacingPolicy * > ( getArrayNoGauss() );
2699 for (int i, j=1; j<=getNumberOfComponents(); j++)
2700 for (i = 1, p_vol=size; p_vol!=lastSize; i++, ++p_vol )
2701 integrale += std::abs( value[indexer->getIndex(i,j)] * *p_vol );
2703 else // FULL_INTERLACE
2705 for ( p_vol=size; p_vol != lastSize; ++p_vol)
2706 for (int j=0; j<getNumberOfComponents(); ++j, ++value)
2707 integrale += std::abs( *value * *p_vol );
2712 // find index for each element of subSupport
2713 PointerOf<int> index;
2714 if ( _support->isOnAllElements() )
2716 index.set( subNums );
2718 else // find common of two partial supports
2720 // hope that numbers are in increasing order
2721 index.set( nbElems );
2722 for (int ii = 0; ii < nbElems; ii++)
2724 bool allNumsFound = true;
2725 int i = 0, iSub = 0;
2726 for ( ; iSub < nbElems; ++iSub )
2728 while ( i < getNumberOfValues() && subNums[iSub] > myNums[i] )
2730 if (i == getNumberOfValues() /*subNums[iSub] > myNums[i]*/) // no more myNums
2732 index[iSub] = 0; // no such number in myNums
2735 else if ( subNums[iSub] == myNums[i] ) // elem number found
2736 index[iSub] = ++i; // -- index counts from 1
2737 else // subNums[iSub] < myNums[i]
2738 allNumsFound = (index[iSub] = 0); // no such number in myNums
2740 if ( iSub != nbElems || !allNumsFound )
2742 // check if numbers are in increasing order
2743 bool increasingOrder = true;
2744 for ( iSub = 1; iSub < nbElems && increasingOrder; ++iSub )
2745 increasingOrder = ( subNums[iSub-1] < subNums[iSub] );
2746 for ( i = 1; i < getNumberOfValues() && increasingOrder; ++i )
2747 increasingOrder = ( myNums[i-1] < myNums[i] );
2749 if ( !increasingOrder )
2750 for ( iSub = 0; iSub < nbElems; ++iSub )
2753 index[iSub] = _support->getValIndFromGlobalNumber( subNums[iSub] );
2755 catch (MEDEXCEPTION)
2763 if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE )
2765 for (int j=0; j<getNumberOfComponents(); ++j)
2767 value = getValue() + j * getNumberOfValues();
2768 for ( int i = 0; i < nbElems; ++i )
2770 integrale += std::abs( value[ index[i]-1 ] * size[i] );
2773 else if ( getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE )
2775 typename ArrayNoByType::InterlacingPolicy* indexer =
2776 dynamic_cast< typename ArrayNoByType::InterlacingPolicy * > ( getArrayNoGauss() );
2777 for (int j=1; j<=getNumberOfComponents(); j++)
2778 for ( int i = 0; i < nbElems; ++i )
2780 integrale += std::abs( value[indexer->getIndex(index[i],j)] * size[i] );
2782 else // FULL_INTERLACE
2784 const int dim = getNumberOfComponents();
2785 for ( int i = 0; i < nbElems; ++i )
2787 for (int j=0; j<dim; ++j)
2788 integrale += std::abs( value[ dim*(index[i]-1) + j] * size[i] );
2791 cellSize->removeReference();
2795 /*! Return a new field (to deallocate with delete) lying on subSupport that is included by
2796 * this->_support with corresponding values extracting from this->_value.
2798 template <class T, class INTERLACING_TAG>
2799 FIELD<T, INTERLACING_TAG>* FIELD<T, INTERLACING_TAG>::extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION)
2801 if(!subSupport->belongsTo(*_support))
2802 throw MEDEXCEPTION("FIELD<T>::extract : subSupport not included in this->_support !");
2803 if(_support->isOnAllElements() && subSupport->isOnAllElements())
2804 return new FIELD<T, INTERLACING_TAG>(*this);
2806 FIELD<T, INTERLACING_TAG> *ret = new FIELD<T, INTERLACING_TAG>(subSupport,
2807 _numberOfComponents);
2810 throw MEDEXCEPTION("FIELD<T>::extract : invalid support detected !");
2812 T* valuesToSet=(T*)ret->getValue();
2814 int nbOfEltsSub=subSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2815 const int *eltsSub=subSupport->getNumber(MED_EN::MED_ALL_ELEMENTS);
2816 T* tempVals=new T[_numberOfComponents];
2817 for(int i=0;i<nbOfEltsSub;i++)
2819 if(!getValueOnElement(eltsSub[i],tempVals))
2820 throw MEDEXCEPTION("Problem in belongsTo function !!!");
2821 for(int j=0;j<_numberOfComponents;j++)
2822 valuesToSet[i*_numberOfComponents+j]=tempVals[j];
2826 ret->copyGlobalInfo(*this);
2834 \addtogroup FIELD_io
2838 Constructor with parameters; the object is set via a file and its associated
2839 driver. For the moment only the MED_DRIVER is considered and if the last two
2840 argument (iterationNumber and orderNumber) are not set; their default value
2841 is -1. If the field fieldDriverName with the iteration number
2842 iterationNumber and the order number orderNumber does not exist in the file
2843 fieldDriverName; the constructor raises an exception.
2845 template <class T, class INTERLACING_TAG>
2846 FIELD<T, INTERLACING_TAG>::FIELD(const SUPPORT * Support,
2847 driverTypes driverType,
2848 const string & fileName/*=""*/,
2849 const string & fieldDriverName/*=""*/,
2850 const int iterationNumber,
2851 const int orderNumber) throw (MEDEXCEPTION)
2853 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) : ";
2860 _mesh = ( MESH* ) NULL;
2862 //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2863 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2864 FIELD_::_valueType=SET_VALUE_TYPE<T>::_valueType;
2866 //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2867 ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2868 FIELD_::_interlacingType=SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2871 //A.G. Addings for RC
2873 _support->addReference();
2874 // OCC 10/03/2006 -- According to the rules defined with help of
2875 // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2876 // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2877 // FIELD template - MSVC++ 2003 compiler generated an error here.
2878 // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2881 _iterationNumber = iterationNumber;
2883 _orderNumber = orderNumber;
2885 current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
2887 _drivers[current]->open();
2888 _drivers[current]->read();
2889 _drivers[current]->close();
2895 If the mesh argument is not initialized or passed NULL,
2896 this constructor, at least, allows to create a FIELD without creating any
2897 SUPPORT then without having to load a MESH object, a support is created. It
2898 provides the meshName related mesh but doesn't not set a mesh in the created
2900 If the passed mesh contains corresponding support, this support will be used
2901 for the field. This support will be found in mesh by name of one of profiles,
2902 on which the FIELD lays in MED-file. This has sense for the case, then MED-file
2903 was created by MEDMEM, and so name of profile contains name of corresponding support.
2905 template <class T, class INTERLACING_TAG>
2906 FIELD<T,INTERLACING_TAG>::FIELD(driverTypes driverType,
2907 const string & fileName,
2908 const string & fieldDriverName,
2909 const int iterationNumber,
2910 const int orderNumber,
2912 throw (MEDEXCEPTION) :FIELD_()
2915 const char* LOC = "FIELD<T,INTERLACING_TAG>::FIELD( driverTypes driverType, const string & fileName, string & fieldDriverName, int iterationNumber, int orderNumber) : ";
2922 _mesh->addReference();
2924 //INITIALISATION DE _valueType DS LE CONSTRUCTEUR DE FIELD_
2925 ASSERT_MED(FIELD_::_valueType == MED_EN::MED_UNDEFINED_TYPE)
2926 FIELD_::_valueType = SET_VALUE_TYPE<T>::_valueType;
2928 //INITIALISATION DE _interlacingType DS LE CONSTRUCTEUR DE FIELD_
2929 ASSERT_MED(FIELD_::_interlacingType == MED_EN::MED_UNDEFINED_INTERLACE)
2930 FIELD_::_interlacingType = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
2932 _support = (SUPPORT *) NULL;
2933 // OCC 10/03/2006 -- According to the rules defined with help of
2934 // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
2935 // MEDMEM_Array<> template using INTERLACING_TAG parameter of
2936 // FIELD template - MSVC++ 2003 compiler generated an error here.
2937 // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
2940 _iterationNumber = iterationNumber;
2942 _orderNumber = orderNumber;
2944 current = addDriver(driverType,fileName,fieldDriverName,MED_EN::RDONLY);
2946 _drivers[current]->open();
2947 _drivers[current]->read();
2948 _drivers[current]->close();
2959 template <class T, class INTERLACING_TAG> FIELD<T, INTERLACING_TAG>::~FIELD()
2961 const char* LOC = " Destructeur FIELD<T, INTERLACING_TAG>::~FIELD()";
2964 if (_value) delete _value; _value=0;
2965 locMap::const_iterator it;
2966 for ( it = _gaussModel.begin();it != _gaussModel.end(); it++ )
2967 delete (*it).second;
2968 _gaussModel.clear();
2970 _mesh->removeReference();
2978 template <class T, class INTERLACING_TAG>
2979 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)
2981 const char* LOC = "FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents)";
2984 _numberOfComponents = NumberOfComponents ;
2985 _componentsTypes.resize(NumberOfComponents);
2986 _componentsNames.resize(NumberOfComponents);
2987 _componentsDescriptions.resize(NumberOfComponents);
2988 _componentsUnits.resize(NumberOfComponents);
2989 _MEDComponentsUnits.resize(NumberOfComponents);
2990 for (int i=0;i<NumberOfComponents;i++) {
2991 _componentsTypes[i] = 0 ;
2995 // becarefull about the number of gauss point
2996 _numberOfValues = _support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
2997 MESSAGE_MED(PREFIX_MED <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
2999 //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
3000 _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
3004 catch (MEDEXCEPTION &ex) {
3005 MESSAGE_MED("No value defined, problem with NumberOfComponents (and may be _support) size of MEDARRAY<T>::_value !");
3006 // OCC 10/03/2006 -- According to the rules defined with help of
3007 // MEDMEM_IntrelacingTraits class, it is not allowed to instantiate
3008 // MEDMEM_Array<> template using INTERLACING_TAG parameter of
3009 // FIELD template - MSVC++ 2003 compiler generated an error here.
3010 // _value = (MEDMEM_Array<T, INTERLACING_TAG> *) NULL;
3021 template <class T, class INTERLACING_TAG>
3022 void FIELD<T, INTERLACING_TAG>::allocValue(const int NumberOfComponents,
3023 const int LengthValue)
3025 const char* LOC = "void FIELD<T>::allocValue(const int NumberOfComponents,const int LengthValue)";
3028 _numberOfComponents = NumberOfComponents ;
3029 _componentsTypes.resize(NumberOfComponents);
3030 _componentsNames.resize(NumberOfComponents);
3031 _componentsDescriptions.resize(NumberOfComponents);
3032 _componentsUnits.resize(NumberOfComponents);
3033 _MEDComponentsUnits.resize(NumberOfComponents);
3034 for (int i=0;i<NumberOfComponents;i++) {
3035 _componentsTypes[i] = 0 ;
3038 MESSAGE_MED("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
3039 _numberOfValues = LengthValue ;
3041 //EF : A modifier lors de l'intégration de la classe de localisation des points de gauss
3042 _value = new ArrayNoGauss(_numberOfComponents,_numberOfValues);
3053 template <class T, class INTERLACING_TAG>
3054 void FIELD<T, INTERLACING_TAG>::deallocValue()
3056 const char* LOC = "void FIELD<T, INTERLACING_TAG>::deallocValue()";
3058 _numberOfValues = 0 ;
3059 _numberOfComponents = 0 ;
3060 if (_value != NULL) {
3072 \addtogroup FIELD_io
3075 // -----------------
3077 // -----------------
3080 Creates the specified driver and return its index reference to path to
3081 read or write methods.
3083 \param driverType specifies the file type (MED_DRIVER or VTK_DRIVER)
3084 \param fileName name of the output file
3085 \param driverName name of the field
3086 \param access access type (read, write or both)
3090 template <class T, class INTERLACING_TAG>
3091 int FIELD<T, INTERLACING_TAG>::addDriver(driverTypes driverType,
3092 const string & fileName/*="Default File Name.med"*/,
3093 const string & driverName/*="Default Field Name"*/,
3094 MED_EN::med_mode_acces access)
3096 //jfa tmp (as last argument has no default value):const char * LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName=\"Default File Name.med\",const string & driverName=\"Default Field Name\",MED_EN::med_mode_acces access) : ";
3100 const char* LOC = "FIELD<T>::addDriver(driverTypes driverType, const string & fileName,const string & driverName,MED_EN::med_mode_acces access) :";
3103 SCRUTE_MED(driverType);
3105 driver = DRIVERFACTORY::buildDriverForField(driverType,fileName,this,access);
3107 _drivers.push_back(driver);
3109 int current = _drivers.size()-1;
3111 _drivers[current]->setFieldName(driverName);
3119 Read FIELD in the file specified in the driver given by its index.
3121 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(int index/*=0*/)
3123 const char * LOC = "FIELD<T, INTERLACING_TAG>::read(int index=0) : ";
3126 if ( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3127 _drivers[index]->open();
3130 _drivers[index]->read();
3132 catch ( const MEDEXCEPTION& ex )
3134 _drivers[index]->close();
3137 _drivers[index]->close();
3140 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3141 << "The index given is invalid, index must be between 0 and |"
3149 Read FIELD with the driver. Additional information (name etc.) to select a field
3150 must be set to the field.
3152 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(const GENDRIVER & driver )
3154 const char* LOC = " FIELD<T, INTERLACING_TAG>::read(const GENDRIVER &) : ";
3157 // For the case where driver does not know about me since it has been created through
3158 // constructor witout parameters: create newDriver knowing me and get missing data
3159 // from driver using merge()
3160 std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
3161 driver.getFileName(),
3163 newDriver->merge( driver );
3170 catch ( const MEDEXCEPTION& ex )
3181 Read FIELD with driver of the given type. Additional information (name etc.) to select a field
3182 must be set to the field.
3184 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::read(driverTypes driverType, const std::string& filename)
3186 const char* LOC = " FIELD<T, INTERLACING_TAG>::read(driverTypes driverType, const std::string& filename) : ";
3189 std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driverType, filename,
3196 catch ( const MEDEXCEPTION& ex )
3206 /*! \if MEDMEM_ug @} \endif */
3209 Duplicates the given driver and return its index reference to path to
3210 read or write methods.
3212 template <class T, class INTERLACING_TAG>
3213 inline int FIELD<T, INTERLACING_TAG>::addDriver (GENDRIVER & driver )
3217 const char* LOC = "FIELD<T, INTERLACING_TAG>::addDriver(GENDRIVER &) : ";
3220 GENDRIVER * newDriver =
3221 DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
3222 driver.getFileName(), this,
3223 driver.getAccessMode());
3224 _drivers.push_back(newDriver);
3226 current = _drivers.size()-1;
3227 SCRUTE_MED(current);
3228 driver.setId(current);
3230 newDriver->merge( driver );
3231 newDriver->setId(current);
3237 Remove the driver referenced by its index.
3239 template <class T, class INTERLACING_TAG>
3240 void FIELD<T, INTERLACING_TAG>::rmDriver (int index/*=0*/)
3242 const char * LOC = "FIELD<T, INTERLACING_TAG>::rmDriver (int index=0): ";
3245 if ( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3246 MESSAGE_MED ("detruire");
3249 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3250 << "The <index given is invalid, index must be between 0 and |"
3259 \addtogroup FIELD_io
3264 Writes FIELD in the file specified by the driver handle \a index.
3269 // Attaching the friver to file "output.med", meshname "Mesh"
3270 int driver_handle = mesh.addDriver(MED_DRIVER, "output.med", "Mesh");
3271 // Writing the content of mesh to the file
3272 mesh.write(driver_handle);
3275 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(int index/*=0*/)
3277 const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
3280 if( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3281 _drivers[index]->open();
3284 _drivers[index]->write();
3286 catch ( const MEDEXCEPTION& ex )
3288 _drivers[index]->close();
3291 _drivers[index]->close();
3294 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3295 << "The index given is invalid, index must be between 0 and |"
3302 Write FIELD with the given driver.
3304 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(const GENDRIVER & driver, MED_EN::med_mode_acces medMode/*=MED_EN::RDWR*/)
3306 const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
3309 // For the case where driver does not know about me since it has been created through
3310 // constructor witout parameters: create newDriver knowing me and get missing data
3311 // from driver using merge()
3312 std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driver.getDriverType(),
3313 driver.getFileName(),
3315 newDriver->merge( driver );
3316 if ( newDriver->getDriverType() == MED_DRIVER )
3317 newDriver->setAccessMode( MED_EN::med_mode_acces( getMedAccessMode( medMode ) ));
3324 catch ( const MEDEXCEPTION& ex )
3335 Write FIELD with driver of the given type.
3337 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::write(driverTypes driverType, const std::string& filename, MED_EN::med_mode_acces medMode/*=MED_EN::RDWR*/)
3339 const char* LOC = " FIELD<T, INTERLACING_TAG>::write(driverTypes driverType, const std::string& filename) : ";
3342 std::auto_ptr<GENDRIVER> newDriver( DRIVERFACTORY::buildDriverForField(driverType, filename,
3344 if ( newDriver->getDriverType() == MED_DRIVER )
3345 newDriver->setAccessMode( MED_EN::med_mode_acces( getMedAccessMode( medMode ) ));
3352 catch ( const MEDEXCEPTION& ex )
3362 /*! \if MEDMEM_ug @} \endif */
3364 Write FIELD in the file specified in the driver given by its index. Use this
3365 method for ASCII drivers (e.g. VTK_DRIVER)
3367 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
3369 const char * LOC = "FIELD<T,INTERLACING_TAG>::write(int index=0, const string & driverName = \"\") : ";
3372 if( index>=0 && index<(int)_drivers.size() && _drivers[index] ) {
3373 _drivers[index]->openAppend();
3374 if (driverName != "") _drivers[index]->setFieldName(driverName);
3377 _drivers[index]->writeAppend();
3379 catch ( const MEDEXCEPTION& ex )
3381 _drivers[index]->close();
3384 _drivers[index]->close();
3387 throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
3388 << "The index given is invalid, index must be between 0 and |"
3397 Write FIELD with the driver which is equal to the given driver.
3399 Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
3401 template <class T, class INTERLACING_TAG> inline void FIELD<T, INTERLACING_TAG>::writeAppend(const GENDRIVER & genDriver)
3403 const char* LOC = " FIELD<T, INTERLACING_TAG>::write(const GENDRIVER &) : ";
3406 for (unsigned int index=0; index < _drivers.size(); index++ )
3407 if ( *_drivers[index] == genDriver ) {
3408 _drivers[index]->openAppend();
3411 _drivers[index]->writeAppend();
3413 catch ( const MEDEXCEPTION& ex )
3415 _drivers[index]->close();
3418 _drivers[index]->close();
3426 Fills in already allocated retValues array the values related to eltIdInSup.
3427 If the element does not exist in this->_support false is returned, true otherwise.
3429 template <class T, class INTERLACING_TAG>
3430 bool FIELD<T, INTERLACING_TAG>::getValueOnElement(int eltIdInSup,T* retValues)
3431 const throw (MEDEXCEPTION)
3436 if(_support->isOnAllElements())
3438 int nbOfEltsThis=_support->getMesh()->getNumberOfElements(_support->getEntity(),MED_EN::MED_ALL_ELEMENTS);
3439 if(eltIdInSup>nbOfEltsThis)
3441 const T* valsThis=getValue();
3442 for(int j=0;j<_numberOfComponents;j++)
3443 retValues[j]=valsThis[(eltIdInSup-1)*_numberOfComponents+j];
3448 int nbOfEltsThis=_support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
3449 const int *eltsThis=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
3452 for(iThis=0;iThis<nbOfEltsThis && !found;)
3453 if(eltsThis[iThis]==eltIdInSup)
3459 const T* valsThis=getValue();
3460 for(int j=0;j<_numberOfComponents;j++)
3461 retValues[j]=valsThis[iThis*_numberOfComponents+j];
3467 * \brief Retrieve value in a given point
3468 * \param coords - point coordinates
3469 * \param output - output buffer
3471 template <class T, class INTERLACING_TAG>
3472 void FIELD<T, INTERLACING_TAG>::getValueOnPoint(const double* coords, double* output) const throw (MEDEXCEPTION)
3474 getValueOnPoints(1, coords, output);
3478 * \brief Retrieve values in given points
3479 * \param nb_points - number of points
3480 * \param coords - point coordinates
3481 * \param output - output buffer
3483 template <class T, class INTERLACING_TAG>
3484 void FIELD<T, INTERLACING_TAG>::getValueOnPoints(int nb_points, const double* coords, double* output) const throw (MEDEXCEPTION)
3486 const char* LOC = " FIELD<T, INTERLACING_TAG>::getValueOnPoints(int nb_points, const double* coords, double* output) : ";
3487 // check operation feasibility
3488 if ( !getSupport() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL Support"));
3489 if ( !getSupport()->getMesh() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL Mesh"));
3490 if ( !_value ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"NULL _value"));
3491 if ( getGaussPresence() ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Not implemeneted for Gauss points"));
3493 MED_EN::medEntityMesh entity = getSupport()->getEntity();
3494 if ( entity != MED_EN::MED_CELL &&
3495 entity != MED_EN::MED_NODE )
3496 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support must be on CELLs or NODEs"));
3498 // initialize output value
3499 for ( int j = 0; j < nb_points*getNumberOfComponents(); ++j )
3502 const MEDMEM::MESH* mesh = getSupport()->getMesh()->convertInMESH();
3503 MEDMEM::AutoDeref derefMesh( mesh );
3505 const double* point = coords;
3506 double* value = output;
3508 if ( entity == MED_EN::MED_CELL )
3510 MEDMEM::PointLocator pLocator (*mesh);
3511 for ( int i = 0; i < nb_points; ++i)
3513 // find the cell enclosing the point
3514 std::list<int> cellIds = pLocator.locate( point );
3515 int nbCells = cellIds.size();
3517 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Point is out of mesh"));
3520 std::list<int>::iterator iCell = cellIds.begin();
3521 for ( ; iCell != cellIds.end(); ++iCell )
3522 for ( int j = 1; j <= getNumberOfComponents(); ++j )
3523 value[j-1] += getValueIJ( *iCell, j ) / nbCells;
3526 point += mesh->getSpaceDimension();
3527 value += getNumberOfComponents();
3530 else // MED_EN::MED_NODE
3532 const double * allCoords = mesh->getCoordinates( MED_EN::MED_FULL_INTERLACE );
3534 MEDMEM::PointLocatorInSimplex pLocator (*mesh);
3535 for ( int i = 0; i < nb_points; ++i)
3537 // find nodes of the simplex enclosing the point
3538 std::list<int> nodeIds = pLocator.locate( point );
3539 int nbNodes = nodeIds.size();
3541 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Point is out of mesh"));
3542 if ( nbNodes != mesh->getMeshDimension() + 1 )
3543 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid nb of points of simplex: "<<nbNodes));
3545 // get coordinates of simplex nodes
3546 std::vector<const double*> nodeCoords( nbNodes );
3547 std::list<int>::iterator iNode = nodeIds.begin();
3549 for ( ; n < nbNodes; ++iNode, ++n )
3550 nodeCoords[n] = allCoords + (*iNode-1) * mesh->getSpaceDimension();
3552 // compute wegths of simplex nodes
3554 pLocator.getNodeWightsInSimplex( nodeCoords, coords, nodeWgt );
3557 for ( n = 0, iNode = nodeIds.begin(); iNode != nodeIds.end(); ++iNode, ++n )
3558 for ( int j = 1; j <= getNumberOfComponents(); ++j )
3559 value[j-1] += getValueIJ( *iNode, j ) * nodeWgt[ n ];
3562 point += mesh->getSpaceDimension();
3563 value += getNumberOfComponents();
3570 Return the coordinates of the gauss points
3571 The returned field has SPACEDIM components
3573 template <class T, class INTERLACING_TAG>
3574 FIELD<double, FullInterlace>* FIELD<T, INTERLACING_TAG>::getGaussPointsCoordinates() const
3575 throw (MEDEXCEPTION)
3577 const char * LOC = "FIELD::getGaussPointsCoordinates() : ";
3581 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3583 const MEDMEM::MESH* mesh = getSupport()->getMesh()->convertInMESH();
3584 MEDMEM::AutoDeref derefMesh( mesh );
3586 const double * coord = mesh->getCoordinates(MED_FULL_INTERLACE);
3587 int spaceDim = mesh->getSpaceDimension();
3589 //Init calculator of the gauss point coordinates
3590 INTERP_KERNEL::GaussCoords calculator;
3591 locMap::const_iterator it;
3593 int nb_type = getSupport()->getNumberOfTypes();
3594 int length_values = getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
3595 const medGeometryElement* types = getSupport()->getTypes();
3596 medEntityMesh support_entity = getSupport()->getEntity();
3597 bool isOnAll = getSupport()->isOnAllElements();
3599 const int* global_connectivity = 0;
3600 const GAUSS_LOCALIZATION<INTERLACING_TAG>* gaussLock = NULL;
3602 typedef typename MEDMEM_ArrayInterface<double,INTERLACING_TAG,NoGauss>::Array ArrayCoord;
3603 typedef typename MEDMEM_ArrayInterface<double,INTERLACING_TAG,Gauss>::Array TArrayGauss;
3605 vector<int> nbelgeoc, nbgaussgeo;
3607 nbelgeoc.resize(nb_type+1, 0);
3608 nbgaussgeo.resize(nb_type+1, 0);
3610 for ( int iType = 0 ; iType < nb_type ; iType++ ) {
3612 medGeometryElement elem_type = types[iType] ;
3613 if(elem_type == MED_EN::MED_POLYGON && elem_type == MED_EN::MED_POLYHEDRA )
3614 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Bad cell type : "<<MED_EN::geoNames[elem_type]<<" !!! "));
3616 it = _gaussModel.find(elem_type);
3618 if(it == _gaussModel.end())
3619 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Gauss localization not defined for "<<MED_EN::geoNames[elem_type]<<" type!!! "));
3620 gaussLock = static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second);
3622 ArrayCoord coord = gaussLock->getGsCoo();
3623 double* gaussCoord = new double[coord.getNbElem()*coord.getDim()];
3625 for( int i = 1 ; i <= coord.getNbElem() ; i++ ) {
3626 for( int j = 1 ; j <= coord.getDim() ; j++ ) {
3627 gaussCoord[idx++] = coord.getIJ(i,j);
3632 ArrayCoord ref = gaussLock->getRefCoo();
3633 double* refCoord = new double[ref.getNbElem()*ref.getDim()];
3634 for( int i = 1 ; i <= ref.getNbElem() ; i++ ) {
3635 for( int j = 1 ; j <= ref.getDim() ; j++ ) {
3636 refCoord[idx++] = ref.getIJ(i,j);
3640 INTERP_KERNEL::NormalizedCellType normType;
3642 case MED_EN::MED_SEG2 : normType = INTERP_KERNEL::NORM_SEG2;break;
3643 case MED_EN::MED_SEG3 : normType = INTERP_KERNEL::NORM_SEG3;break;
3644 default : normType = (INTERP_KERNEL::NormalizedCellType) ((((unsigned long)elem_type/100-2)*10) + ((unsigned long)elem_type%100));
3648 calculator.addGaussInfo(normType,
3651 gaussLock->getNbGauss(),
3655 // Preapre Info for the gauss array
3656 nbelgeoc [ iType+1 ] = nbelgeoc[ iType ] + getSupport()->getNumberOfElements(elem_type);
3657 nbgaussgeo [ iType+1 ] = gaussLock->getNbGauss();
3659 delete [] gaussCoord;
3663 FIELD<double, FullInterlace>* gpCoord =
3664 new FIELD<double, FullInterlace>(getSupport(),spaceDim);
3665 gpCoord->setName("Gauss Points Coordinates");
3666 gpCoord->setDescription("Gauss Points Coordinates");
3668 for(int dimId = 1 ; dimId <= spaceDim; dimId++) {
3671 gpCoord->setComponentName(dimId,"X");
3672 gpCoord->setComponentDescription(dimId,"X coordinate of the gauss point");
3675 gpCoord->setComponentName(dimId,"Y");
3676 gpCoord->setComponentDescription(dimId,"Y coordinate of the gauss point");
3679 gpCoord->setComponentName(dimId,"Z");
3680 gpCoord->setComponentDescription(dimId,"Z coordinate of the gauss point");
3684 gpCoord->setMEDComponentUnit(dimId, mesh->getCoordinatesUnits()[dimId-1]);
3687 gpCoord->setIterationNumber(getIterationNumber());
3688 gpCoord->setOrderNumber(getOrderNumber());
3689 gpCoord->setTime(getTime());
3691 TArrayGauss *arrayGauss = new TArrayGauss(spaceDim, length_values,
3692 nb_type, & nbelgeoc[0], & nbgaussgeo[0]);
3693 gpCoord->setArray(arrayGauss);
3698 //Calculation of the coordinates
3700 for ( int i = 0 ; i < nb_type ; i++ ) {
3702 medGeometryElement type = types[i] ;
3703 INTERP_KERNEL::NormalizedCellType normType;
3705 case MED_EN::MED_SEG2 : normType = INTERP_KERNEL::NORM_SEG2;break;
3706 case MED_EN::MED_SEG3 : normType = INTERP_KERNEL::NORM_SEG3;break;
3707 default : normType = (INTERP_KERNEL::NormalizedCellType) ((((unsigned long)type/100-2)*10) + ((unsigned long)type%100));
3711 it = _gaussModel.find(type);
3713 gaussLock = static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> * > ((*it).second);
3714 int nb_entity_type = getSupport()->getNumberOfElements(type);
3718 global_connectivity = mesh->getConnectivity(MED_NODAL,support_entity,type);
3721 const int * supp_number = getSupport()->getNumber(type);
3722 const int * connectivity = mesh->getConnectivity(MED_NODAL,support_entity,MED_ALL_ELEMENTS);
3723 const int * connectivityIndex = mesh->getConnectivityIndex(MED_NODAL,support_entity);
3724 int * global_connectivity_tmp = new int[(type%100)*nb_entity_type];
3726 for (int k_type = 0; k_type<nb_entity_type; k_type++) {
3727 for (int j_ent = 0; j_ent<(type%100); j_ent++) {
3728 global_connectivity_tmp[k_type*(type%100)+j_ent] = connectivity[connectivityIndex[supp_number[k_type]-1]+j_ent-1];
3731 global_connectivity = global_connectivity_tmp;
3734 int nbNodes = (type%100);
3735 double* gCoord = NULL;
3738 for ( int elem = 0; elem < nb_entity_type; elem++ ) {
3739 int elem_index = nbNodes*elem;
3740 Ni = new int[nbNodes];
3741 for( int idx = 0 ; idx < nbNodes; idx++ ) {
3742 Ni[idx] = global_connectivity[ elem_index+idx ] - 1;
3745 gCoord = calculator.calculateCoords(normType,
3749 int resultIndex = 0;
3750 for( int k = 0; k < gaussLock->getNbGauss(); k++ ) {
3751 for( int dimId = 1; dimId <= spaceDim; dimId++ ) {
3752 gpCoord->setValueIJK(index,dimId,(k+1),gCoord[resultIndex]);
3760 if (!isOnAll && type != MED_EN::MED_POLYHEDRA && type != MED_EN::MED_POLYGON) {
3761 delete [] global_connectivity ;
3770 Destroy the MEDARRAY<T> in FIELD and put the new one without copy.
3773 template <class T, class INTERLACING_TAG>
3774 inline void FIELD<T, INTERLACING_TAG>::setArray(MEDMEM_Array_ * Value)
3775 throw (MEDEXCEPTION)
3777 if (NULL != _value) delete _value ;
3783 Return a reference to the MEDARRAY<T, INTERLACING_TAG> in FIELD.
3786 template <class T, class INTERLACING_TAG>
3787 inline MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() const throw (MEDEXCEPTION)
3789 const char* LOC = "MEDMEM_Array_ * FIELD<T, INTERLACING_TAG>::getArray() : ";
3794 template <class T,class INTERLACING_TAG> inline
3795 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,Gauss>::Array *
3796 FIELD<T, INTERLACING_TAG>::getArrayGauss() const throw (MEDEXCEPTION)
3798 const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayGauss() : ";
3801 if ( getGaussPresence() )
3802 return static_cast<ArrayGauss *> (_value);
3804 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
3805 "The field has no Gauss Point"));
3811 template <class T,class INTERLACING_TAG> inline
3812 typename MEDMEM_ArrayInterface<T,INTERLACING_TAG,NoGauss>::Array *
3813 FIELD<T, INTERLACING_TAG>::getArrayNoGauss() const throw (MEDEXCEPTION)
3815 const char * LOC = "FIELD<T, INTERLACING_TAG>::getArrayNoGauss() : ";
3818 if ( ! getGaussPresence() )
3819 return static_cast < ArrayNoGauss * > (_value);
3821 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<
3822 "The field has Gauss Point"));
3828 template <class T,class INTERLACING_TAG> inline bool
3829 FIELD<T, INTERLACING_TAG>::getGaussPresence() const throw (MEDEXCEPTION)
3832 return _value->getGaussPresence();
3834 throw MEDEXCEPTION("FIELD<T, INTERLACING_TAG>::getGaussPresence() const : Can't call getGaussPresence on a null _value");
3838 Return the actual length of the reference to values array returned by getValue.
3839 Take care of number of components and number of Gauss points by geometric type
3841 template <class T, class INTERLACING_TAG>
3842 inline int FIELD<T, INTERLACING_TAG>::getValueLength() const
3843 throw (MEDEXCEPTION)
3845 if ( getGaussPresence() )
3846 return static_cast<ArrayGauss *>(_value)->getArraySize() ;
3848 return static_cast<ArrayNoGauss *>(_value)->getArraySize() ;
3853 \defgroup FIELD_value Field values
3855 These methods are provided for accessing the values of a field.
3856 There are two ways to do so : one consists in using accessors
3857 that retrieve elements or group of elements from the entire field. Typical use is
3859 FIELD field(MED_DRIVER, "result.med","Pressure");
3860 double P0=field.getValueIJ(1,1);
3863 Another way is to retrieve the pointer to the array that contains the
3864 variable values. In this case, the user should be aware of the interlacing mode
3865 so that no mistakes are made when retrieving the values.
3868 FIELD field(MED_DRIVER, "result.med","Pressure");
3869 double* ptrP=field.getValue();
3877 Returns a reference to values array to read them.
3879 template <class T, class INTERLACIN_TAG>
3880 inline const T* FIELD<T, INTERLACIN_TAG>::getValue() const throw (MEDEXCEPTION)
3882 const char* LOC = "FIELD<T, INTERLACING_TAG>::getValue() : ";
3884 if ( getGaussPresence() )
3885 return static_cast<ArrayGauss *>(_value)->getPtr() ;
3887 return static_cast<ArrayNoGauss *>(_value)->getPtr() ;
3890 Returns a reference to \f$ i^{th} \f$ row
3891 of FIELD values array.
3892 If a faster accessor is intended you may use getArray() once,
3893 then MEDMEM_Array accessors.
3894 Be careful if field support is not on all elements getRow
3895 use support->getValIndFromGlobalNumber(i).
3897 template <class T,class INTERLACING_TAG> inline
3899 FIELD<T,INTERLACING_TAG>::getRow(int i) const throw (MEDEXCEPTION)
3901 const char* LOC; LOC = "FIELD<T,INTERLACING_TAG>::getRow(int i) : ";
3905 valIndex = _support->getValIndFromGlobalNumber(i);
3907 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3909 if ( getGaussPresence() )
3910 return static_cast<ArrayGauss *>(_value)->getRow(valIndex) ;
3912 return static_cast<ArrayNoGauss *>(_value)->getRow(valIndex) ;
3916 Returns a reference to $j^{th}$ column
3917 of FIELD values array.
3919 template <class T,class INTERLACING_TAG> inline const T*
3920 FIELD<T,INTERLACING_TAG>::getColumn(int j) const throw (MEDEXCEPTION)
3922 if ( getGaussPresence() )
3923 return static_cast<ArrayGauss *>(_value)->getColumn(j) ;
3925 return static_cast<ArrayNoGauss *>(_value)->getColumn(j) ;
3929 Returns the value of $i^{th}$ element and $j^{th}$ component.
3930 This method only works with fields having no particular Gauss point
3931 definition (i.e., fields having one value per element).
3932 This method makes the retrieval of the value independent from the
3933 interlacing pattern, but it is slower than the complete retrieval
3934 obtained by the \b getValue() method.
3936 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJ(int i,int j) const throw (MEDEXCEPTION)
3938 const char * LOC = "getValueIJ(..)";
3941 valIndex = _support->getValIndFromGlobalNumber(i);
3943 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3945 if ( getGaussPresence() )
3946 return static_cast<ArrayGauss *>(_value)->getIJ(valIndex,j) ;
3948 return static_cast<ArrayNoGauss *>(_value)->getIJ(valIndex,j) ;
3952 Returns the $j^{th}$ component of $k^{th}$ Gauss points of $i^{th}$ value.
3953 This method is compatible with elements having more than one Gauss point.
3954 This method makes the retrieval of the value independent from the
3955 interlacing pattern, but it is slower than the complete retrieval
3956 obtained by the \b getValue() method.
3958 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION)
3960 const char * LOC = "getValueIJK(..)";
3963 valIndex = _support->getValIndFromGlobalNumber(i);
3965 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
3967 if ( getGaussPresence() )
3968 return static_cast<ArrayGauss *>(_value)->getIJK(valIndex,j,k) ;
3970 return static_cast<ArrayNoGauss *>(_value)->getIJK(valIndex,j,k) ;
3972 /*! \if MEDMEM_ug @} \endif */
3975 Return number of values of a geomertic type in NoInterlaceByType mode
3977 template <class T, class INTERLACIN_TAG>
3978 inline int FIELD<T, INTERLACIN_TAG>::getValueByTypeLength(int t) const throw (MEDEXCEPTION)
3980 const char * LOC ="getValueByTypeLength() : ";
3981 if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
3982 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
3984 if ( getGaussPresence() ) {
3985 ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
3986 if ( t < 1 || t > array->getNbGeoType() )
3987 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
3988 return array->getLengthOfType( t );
3991 ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
3992 if ( t < 1 || t > array->getNbGeoType() )
3993 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Invalid type: "<< t ));
3994 return array->getLengthOfType( t );
3999 Return a reference to values array to read them.
4001 template <class T, class INTERLACIN_TAG>
4002 inline const T* FIELD<T, INTERLACIN_TAG>::getValueByType(int t) const throw (MEDEXCEPTION)
4004 if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4005 throw MEDEXCEPTION(LOCALIZED("getValueByType() : not MED_NO_INTERLACE_BY_TYPE field" ));
4007 if ( getGaussPresence() ) {
4008 ArrayNoByTypeGauss* array = static_cast<ArrayNoByTypeGauss *>(_value);
4009 return array->getPtr() + array->getIndex( t );
4012 ArrayNoByType* array = static_cast<ArrayNoByType *>(_value);
4013 return array->getPtr() + array->getIndex( t );
4018 Return the value of i^{th} element in indicated type t and j^{th} component.
4020 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJByType(int i,int j, int t) const throw (MEDEXCEPTION)
4022 const char * LOC = "getValueIJByType(..)";
4023 if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4024 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4026 if ( getGaussPresence() )
4027 return static_cast<ArrayNoByTypeGauss *>(_value)->getIJByType(i,j,t) ;
4029 return static_cast<ArrayNoByType *>(_value)->getIJByType(i,j,t) ;
4033 Return the j^{th} component of k^{th} gauss points of i^{th} value with type t.
4035 template <class T,class INTERLACING_TAG> inline T FIELD<T,INTERLACING_TAG>::getValueIJKByType(int i,int j,int k,int t) const throw (MEDEXCEPTION)
4037 const char * LOC = "getValueIJKByType(..)";
4038 if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4039 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4041 if ( getGaussPresence() )
4042 return static_cast<ArrayNoByTypeGauss *>(_value)->getIJKByType(i,j,k,t) ;
4044 return static_cast<ArrayNoByType *>(_value)->getIJKByType(i,j,k,t) ;
4048 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGeometricTypes() const throw (MEDEXCEPTION)
4050 const char * LOC = "getNumberOfGeometricTypes(..)";
4053 return _support->getNumberOfTypes();
4055 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4060 \addtogroup FIELD_gauss
4065 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> &
4066 FIELD<T,INTERLACING_TAG>::getGaussLocalization(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
4068 const char * LOC ="getGaussLocalization(MED_EN::medGeometryElement geomElement) : ";
4069 const GAUSS_LOCALIZATION_ * locPtr=0;
4071 locMap::const_iterator it;
4072 if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4073 locPtr = (*it).second;
4074 return *static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
4077 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
4081 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION<INTERLACING_TAG> *
4082 FIELD<T,INTERLACING_TAG>::getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
4084 const char * LOC ="getGaussLocalizationPtr(MED_EN::medGeometryElement geomElement) : ";
4085 const GAUSS_LOCALIZATION_ * locPtr=0;
4087 locMap::const_iterator it;
4088 if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4089 locPtr = (*it).second;
4090 return static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr);
4093 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type" ));
4098 * \brief Return GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
4100 template <class T,class INTERLACING_TAG> const GAUSS_LOCALIZATION_ *
4101 FIELD<T,INTERLACING_TAG>::getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION)
4103 const char * LOC ="getGaussLocalizationRoot(MED_EN::medGeometryElement geomElement) : ";
4105 locMap::const_iterator it;
4106 if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4107 return (*it).second;
4110 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Can't find any GaussLocalization on this geometric type: "<< geomElement ));
4115 * \brief Take onership of GAUSS_LOCALIZATION_* whose interlacing type may differ from one of the field
4117 template <class T,class INTERLACING_TAG> void
4118 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, GAUSS_LOCALIZATION_* gaussloc)
4120 locMap::iterator it = _gaussModel.find(geomElement);
4121 if ( it != _gaussModel.end() ) {
4123 it->second = gaussloc;
4126 _gaussModel[ geomElement ] = gaussloc;
4131 template <class T,class INTERLACING_TAG> void
4132 FIELD<T,INTERLACING_TAG>::setGaussLocalization(MED_EN::medGeometryElement geomElement, const GAUSS_LOCALIZATION<INTERLACING_TAG>& gaussloc)
4134 locMap::iterator it = _gaussModel.find(geomElement);
4135 if ( it != _gaussModel.end() ) {
4137 it->second = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
4140 _gaussModel[ geomElement ] = new GAUSS_LOCALIZATION<INTERLACING_TAG> (gaussloc);
4145 Returns number of Gauss points for this medGeometryElement.
4148 if there is no GAUSS_LOCALIZATION having this medGeometryElement but
4149 the medGeometryElement exist in the SUPPORT, getNumberOfGaussPoints
4152 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) const
4153 throw (MEDEXCEPTION)
4155 const char * LOC ="getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
4156 const GAUSS_LOCALIZATION_ * locPtr=0;
4158 locMap::const_iterator it;
4159 if ( ( it = _gaussModel.find(geomElement)) != _gaussModel.end() ) {
4160 locPtr = (*it).second;
4161 return static_cast<const GAUSS_LOCALIZATION<INTERLACING_TAG> *>(locPtr)->getNbGauss();
4166 if ( _support->getNumberOfElements(geomElement) ) return 1;
4167 } catch ( MEDEXCEPTION & ex) {
4168 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<< "GeometricType not found !" )) ;
4171 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4173 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Should never execute this!" ));
4178 Returns number of Gauss points for each geometric type.
4181 if there is no gauss points whatever the geometric type is
4182 it returns an exception. (renvoyer un tableau de 1 ?)
4184 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfGaussPoints() const
4185 throw (MEDEXCEPTION)
4187 const char * LOC ="const int * getNumberOfGaussPoints(MED_EN::medGeometryElement geomElement) : ";
4190 if ( getGaussPresence() ) {
4191 return static_cast<ArrayGauss *>(_value)->getNbGaussGeo()+1;
4193 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"value hasn't Gauss points " ));
4196 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Value not defined" ));
4200 Returns number of Gauss points for element n°i.
4201 The i index is a global index (take care of previous element
4202 on different geometric type).
4204 template <class T,class INTERLACING_TAG> const int FIELD<T,INTERLACING_TAG>::getNbGaussI(int i) const throw (MEDEXCEPTION)
4206 const char * LOC = "getNbGaussI(..)";
4210 valIndex = _support->getValIndFromGlobalNumber(i);
4212 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4215 if ( getGaussPresence() )
4216 return static_cast<ArrayGauss *>(_value)->getNbGauss(valIndex) ;
4218 return static_cast<ArrayNoGauss *>(_value)->getNbGauss(valIndex) ;
4220 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_value not defined" ));
4227 template <class T,class INTERLACING_TAG> const int * FIELD<T,INTERLACING_TAG>::getNumberOfElements() const throw (MEDEXCEPTION)
4229 const char * LOC = "getNumberOfElements(..)";
4232 return _support->getNumberOfElements();
4234 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4238 template <class T,class INTERLACING_TAG> const MED_EN::medGeometryElement * FIELD<T,INTERLACING_TAG>::getGeometricTypes() const throw (MEDEXCEPTION)
4240 const char * LOC = "getGeometricTypes(..)";
4243 return _support->getTypes();
4245 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4248 template <class T,class INTERLACING_TAG> bool FIELD<T,INTERLACING_TAG>::isOnAllElements() const throw (MEDEXCEPTION)
4250 const char * LOC = "isOnAllElements(..)";
4253 return _support->isOnAllElements();
4255 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not defined" ));
4261 \addtogroup FIELD_value
4266 Copy new values array in FIELD according to the given mode.
4268 Array must have right size. If not results are unpredicable.
4269 In MED_FULL_INTERLACE mode, values are stored elementwise in X1,Y1,Z1,X2,Y2,Z2.. order.
4270 In MED_NO_INTERLACE mode, values are stored componentwise in X1,X2,X3,...,Y1,Y2,Y3,... order.
4272 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValue( T* value) throw (MEDEXCEPTION)
4274 if ( getGaussPresence() )
4275 static_cast<ArrayGauss *>(_value)->setPtr(value) ;
4277 static_cast<ArrayNoGauss *>(_value)->setPtr(value) ;
4281 Update values array in the j^{th} row of FIELD values array with the given ones and
4282 according to specified mode.
4284 template <class T,class INTERLACING_TAG>
4285 inline void FIELD<T,INTERLACING_TAG>::setRow( int i, T* value) throw (MEDEXCEPTION)
4287 const char * LOC = "FIELD<T,INTERLACING_TAG>::setRow(int i, T* value) : ";
4290 valIndex = _support->getValIndFromGlobalNumber(i);
4292 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
4294 if ( getGaussPresence() )
4295 static_cast<ArrayGauss *>(_value)->setRow(valIndex, value) ;
4297 static_cast<ArrayNoGauss *>(_value)->setRow(valIndex, value) ;
4301 Update values array in the $j^{th}$ column of FIELD values array with the given ones and
4302 according to specified mode.
4304 template <class T,class INTERLACING_TAG>
4305 inline void FIELD<T,INTERLACING_TAG>::setColumn( int j, T* value) throw (MEDEXCEPTION)
4307 if ( getGaussPresence() )
4308 static_cast<ArrayGauss *>(_value)->setColumn(j, value) ;
4310 static_cast<ArrayNoGauss *>(_value)->setColumn(j, value) ;
4314 Sets the value of i^{th} element and j^{th} component with the given one.
4316 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) throw (MEDEXCEPTION)
4318 const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
4321 valIndex = _support->getValIndFromGlobalNumber(i);
4323 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
4325 if ( getGaussPresence() )
4326 static_cast<ArrayGauss *>(_value)->setIJ(valIndex,j,value) ;
4328 static_cast<ArrayNoGauss *>(_value)->setIJ(valIndex,j,value) ;
4332 Set the value of i^{th} element, j^{th} component and k^{th} gauss point with the given one.
4334 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJK(int i, int j, int k, T value) throw (MEDEXCEPTION)
4336 const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJ(int i, int j, T value) : ";
4339 valIndex = _support->getValIndFromGlobalNumber(i);
4341 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Support not define |" ));
4343 if ( getGaussPresence() )
4344 static_cast<ArrayGauss *>(_value)->setIJK(valIndex,j,k,value) ;
4346 static_cast<ArrayNoGauss *>(_value)->setIJK(valIndex,j,k,value) ;
4350 Set the value of i^{th} element and j^{th} component with the given one.
4352 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) throw (MEDEXCEPTION)
4354 const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJByType(int i, int j, int t, T value) : ";
4355 if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4356 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4358 if ( getGaussPresence() )
4359 return static_cast<ArrayNoByTypeGauss *>(_value)->setIJByType(i,j,t,value) ;
4361 return static_cast<ArrayNoByType *>(_value)->setIJByType(i,j,t,value) ;
4365 Set the value of component of k^{th} gauss points of i^{th} element and j^{th} component with the given one.
4367 template <class T,class INTERLACING_TAG> inline void FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int k, int t, T value) throw (MEDEXCEPTION)
4369 const char * LOC = "FIELD<T,INTERLACING_TAG>::setValueIJKByType(int i, int j, int t, int k, T value) : ";
4370 if ( getInterlacingType() != MED_EN::MED_NO_INTERLACE_BY_TYPE )
4371 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"not MED_NO_INTERLACE_BY_TYPE field" ));
4373 if ( getGaussPresence() )
4374 return static_cast<ArrayNoByTypeGauss *>(_value)->setIJKByType(i,j,k,t,value) ;
4376 return static_cast<ArrayNoByType *>(_value)->setIJKByType(i,j,k,t,value) ;
4378 /*! \if MEDMEM_ug @} \endif */
4385 Fill array by using T_Analytic.
4386 WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
4387 Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
4389 template <class T, class INTERLACING_TAG>
4390 void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION)
4392 const char * LOC = "void FIELD<T, INTERLACING_TAG>::fillFromAnalytic(myFuncType f) : ";
4394 if (_support == (SUPPORT *) NULL)
4395 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"No Support defined."));
4397 const GMESH * mesh = _support->getMesh();
4398 int spaceDim = mesh->getSpaceDimension();
4399 const double * coord;
4401 const double * bary;
4402 FIELD<double,FullInterlace> * barycenterField=0;
4404 double ** xyz=new double* [spaceDim];
4405 bool deallocateXyz=false;
4406 if(_support->getEntity()==MED_EN::MED_NODE)
4408 const MESH * unstructured = _support->getMesh()->convertInMESH();
4409 if (_support->isOnAllElements())
4411 coord=unstructured->getCoordinates(MED_EN::MED_NO_INTERLACE);
4412 for(i=0; i<spaceDim; i++)
4413 xyz[i]=(double *)coord+i*_numberOfValues;
4417 coord = unstructured->getCoordinates(MED_EN::MED_FULL_INTERLACE);
4418 const int * nodesNumber=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
4419 for(i=0; i<spaceDim; i++)
4420 xyz[i]=new double[_numberOfValues];
4422 for(i=0;i<_numberOfValues;i++)
4424 for(j=0;j<spaceDim;j++)
4425 xyz[j][i]=coord[(nodesNumber[i]-1)*spaceDim+j];
4428 unstructured->removeReference();
4432 barycenterField = mesh->getBarycenter(_support);
4433 bary = barycenterField->getValue();
4434 for(i=0; i<spaceDim; i++)
4435 xyz[i]=new double[_numberOfValues];
4437 for(i=0;i<_numberOfValues;i++) {
4438 for(j=0;j<spaceDim;j++)
4439 xyz[j][i]=bary[i*spaceDim+j];
4442 T* valsToSet=(T*)getValue();
4443 double *temp=new double[spaceDim];
4444 for(i=0;i<_numberOfValues;i++)
4446 for(j=0;j<spaceDim;j++)
4448 f(temp,valsToSet+i*_numberOfComponents);
4452 delete barycenterField;
4454 for(j=0;j<spaceDim;j++)
4460 Execute a function on _values on 'this' and put the result on a newly created field that has to be deallocated.
4461 WARNING : "this" must have allocated its array by setting this->_support and this->_numberOfComponents properly.
4462 Typically you should use it on a field built with constructor FIELD<T>::FIELD<T>(SUPPORT *,int nbOfComponents)
4464 template <class T, class INTERLACING_TAG>
4465 FIELD<T,INTERLACING_TAG> *FIELD<T, INTERLACING_TAG>::execFunc(int nbOfComponents,
4466 myFuncType2 f) throw (MEDEXCEPTION)
4468 FIELD<T,INTERLACING_TAG> *ret=new FIELD<T,INTERLACING_TAG>(_support,nbOfComponents);
4469 const T* valsInput=getValue();
4470 T* valsOutPut=(T*)ret->getValue();
4472 for(i=0;i<_numberOfValues;i++)
4473 f(valsInput+i*_numberOfComponents,valsOutPut+i*nbOfComponents);
4477 }//End namespace MEDMEM
4479 #endif /* FIELD_HXX */