-// MED MEDMEM : MED files in memory
-//
-// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
-//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
-//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//
-// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
-//
-//
-//
-// File : MEDMEM_Field.hxx
-// Module : MED
-
/*
File Field.hxx
$Header$
#define FIELD_HXX
#include <vector>
+#include <algorithm>
+#include <cmath>
#include "utilities.h"
#include "MEDMEM_Exception.hxx"
#include "MEDMEM_MedFieldDriver.hxx"
#include "MEDMEM_MedMedDriver.hxx"
+#include "MEDMEM_VtkFieldDriver.hxx"
+#include "MEDMEM_VtkMedDriver.hxx"
+
using namespace MED_EN;
/*!
*/
+namespace MEDMEM {
class FIELD_ // GENERIC POINTER TO a template <class T> class FIELD
{
protected:
double _time;
int _orderNumber ;
+ // _valueType should be a static const. Here is an initialization exemple
+ // template < classType T > struct SET_VALUE_TYPE { static const med_type_champ _valueType = 0; }
+ // template < > struct SET_VALUE_TYPE<double> { static const med_type_champ _valueType = MED_EN::MED_REEL64; }
+ // template < > struct SET_VALUE_TYPE<int> { static const med_type_champ _valueType = MED_EN::MED_INT32; }
+ // static const med_type_champ _valueType = SET_VALUE_TYPE <T>::_valueType;
med_type_champ _valueType ;
vector<GENDRIVER *> _drivers; // Storage of the drivers currently in use
+ static void _checkFieldCompatibility(const FIELD_& m, const FIELD_& n ) throw (MEDEXCEPTION);
+ void _checkNormCompatibility(const FIELD<double>* p_field_volume=NULL) const throw (MEDEXCEPTION);
+ FIELD<double>* _getFieldSize() const;
public:
friend class MED_MED_WRONLY_DRIVER;
friend class MED_MED_RDWR_DRIVER;
+ friend class VTK_MED_DRIVER;
+
/*!
Constructor.
*/
/*!
Destructor.
*/
- ~FIELD_();
+ virtual ~FIELD_();
// virtual void setIterationNumber (int IterationNumber);
// virtual void setOrderNumber (int OrderNumber);
virtual int addDriver( GENDRIVER & driver);
virtual void read (const GENDRIVER &);
virtual void read(int index=0);
+ virtual void openAppend( void );
virtual void write(const GENDRIVER &);
virtual void write(int index=0, const string & driverName="");
- inline void setName(string Name);
+ virtual void writeAppend(const GENDRIVER &);
+ virtual void writeAppend(int index=0, const string & driverName="");
+
+ inline void setName(const string Name);
inline string getName() const;
- inline void setDescription(string Description);
+ inline void setDescription(const string Description);
inline string getDescription() const;
inline const SUPPORT * getSupport() const;
- inline void setSupport(SUPPORT * support);
- inline void setNumberOfComponents(int NumberOfComponents);
+ inline void setSupport(const SUPPORT * support);
+ inline void setNumberOfComponents(const int NumberOfComponents);
inline int getNumberOfComponents() const;
- inline void setNumberOfValues(int NumberOfValues);
+ inline void setNumberOfValues(const int NumberOfValues);
inline int getNumberOfValues() const;
// inline void setComponentType(int *ComponentType);
// inline int * getComponentType() const;
// inline int getComponentTypeI(int i) const;
- inline void setComponentsNames(string * ComponentsNames);
- inline void setComponentName(int i, string ComponentName);
+ inline void setComponentsNames(const string * ComponentsNames);
+ inline void setComponentName(int i, const string ComponentName);
inline const string * getComponentsNames() const;
inline string getComponentName(int i) const;
- inline void setComponentsDescriptions(string *ComponentsDescriptions);
- inline void setComponentDescription(int i, string ComponentDescription);
+ inline void setComponentsDescriptions(const string * ComponentsDescriptions);
+ inline void setComponentDescription(int i, const string ComponentDescription);
inline const string * getComponentsDescriptions() const;
inline string getComponentDescription(int i) const;
// provisoire : en attendant de regler le probleme des unites !
- inline void setComponentsUnits(UNIT * ComponentsUnits);
+ inline void setComponentsUnits(const UNIT * ComponentsUnits);
inline const UNIT * getComponentsUnits() const;
inline const UNIT * getComponentUnit(int i) const;
- inline void setMEDComponentsUnits(string * MEDComponentsUnits);
- inline void setMEDComponentUnit(int i, string MEDComponentUnit);
+ inline void setMEDComponentsUnits(const string * MEDComponentsUnits);
+ inline void setMEDComponentUnit(int i, const string MEDComponentUnit);
inline const string * getMEDComponentsUnits() const;
inline string getMEDComponentUnit(int i) const;
inline void setOrderNumber(int OrderNumber);
inline int getOrderNumber() const;
- inline void setValueType (med_type_champ ValueType) ;
+ inline void setValueType (const med_type_champ ValueType) ;
inline med_type_champ getValueType () const;
};
+};
// ---------------------------------
// Implemented Methods : constructor
// ---------------------------------
+using namespace MEDMEM;
// -----------------
// Methodes Inline
/*!
Set FIELD name.
*/
-inline void FIELD_::setName(string Name)
+inline void FIELD_::setName(const string Name)
{
_name=Name;
}
/*!
Set FIELD description.
*/
-inline void FIELD_::setDescription(string Description)
+inline void FIELD_::setDescription(const string Description)
{
_description=Description;
}
Duplicate the ComponentsNames string array to put components names in
FIELD. ComponentsNames size must be equal to number of components.
*/
-inline void FIELD_::setComponentsNames(string * ComponentsNames)
+inline void FIELD_::setComponentsNames(const string * ComponentsNames)
{
if (NULL == _componentsNames)
_componentsNames = new string[_numberOfComponents] ;
i must be >=1 and <= number of components.
*/
-inline void FIELD_::setComponentName(int i, string ComponentName)
+inline void FIELD_::setComponentName(int i, const string ComponentName)
{
_componentsNames[i-1]=ComponentName ;
}
descriptions in FIELD.
ComponentsDescriptions size must be equal to number of components.
*/
-inline void FIELD_::setComponentsDescriptions(string *ComponentsDescriptions)
+inline void FIELD_::setComponentsDescriptions(const string * ComponentsDescriptions)
{
if (NULL == _componentsDescriptions)
_componentsDescriptions = new string[_numberOfComponents] ;
i must be >=1 and <= number of components.
*/
-inline void FIELD_::setComponentDescription(int i,string ComponentDescription)
+inline void FIELD_::setComponentDescription(int i,const string ComponentDescription)
{
_componentsDescriptions[i-1]=ComponentDescription ;
}
units in FIELD.
ComponentsUnits size must be equal to number of components.
*/
-inline void FIELD_::setComponentsUnits(UNIT * ComponentsUnits)
+inline void FIELD_::setComponentsUnits(const UNIT * ComponentsUnits)
{
if (NULL == _componentsUnits)
_componentsUnits = new UNIT[_numberOfComponents] ;
MEDComponentsUnits size must be equal to number of components.
*/
-inline void FIELD_::setMEDComponentsUnits(string * MEDComponentsUnits)
+inline void FIELD_::setMEDComponentsUnits(const string * MEDComponentsUnits)
{
if (NULL == _MEDComponentsUnits)
_MEDComponentsUnits = new string[_numberOfComponents] ;
i must be >=1 and <= number of components.
*/
-inline void FIELD_::setMEDComponentUnit(int i, string MEDComponentUnit)
+inline void FIELD_::setMEDComponentUnit(int i, const string MEDComponentUnit)
{
_MEDComponentsUnits[i-1]=MEDComponentUnit ;
}
Reference is not duplicate, so it must not be deleted.
*/
-inline void FIELD_::setSupport(SUPPORT * support)
+inline void FIELD_::setSupport(const SUPPORT * support)
{
_support = support ;
}
/*!
Set the FIELD med value type (MED_INT32 or MED_REEL64).
*/
-inline void FIELD_::setValueType (med_type_champ ValueType)
+inline void FIELD_::setValueType (const med_type_champ ValueType)
{
_valueType = ValueType ;
}
*/
+namespace MEDMEM {
template <class T> class FIELD : public FIELD_
{
// static INSTANCE_DE<MED_FIELD_RDONLY_DRIVER> inst_med_rdonly ;
static INSTANCE_DE<MED_FIELD_RDWR_DRIVER<T> > inst_med ;
+ static INSTANCE_DE<VTK_FIELD_DRIVER<T> > inst_vtk ;
+
static const INSTANCE * const instances[] ;
// ------ End of Drivers Management Part
MEDARRAY<T> *_value ;
private:
+ void _operation(const FIELD& m,const FIELD& n, const medModeSwitch mode, char* Op);
+ void _operationInitialize(const FIELD& m,const FIELD& n, char* Op);
+ void _add_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
+ void _sub_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
+ void _mul_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
+ void _div_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode);
//setValueType() ;
public:
FIELD();
FIELD(const FIELD &m);
FIELD & operator=(const FIELD &m); // A FAIRE
- FIELD(const SUPPORT * Support, const int NumberOfComponents); // Ajout NB Constructeur FIELD avec allocation de memoire de tous ses attribut
+ FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode=MED_FULL_INTERLACE) throw (MEDEXCEPTION) ; // Ajout NB Constructeur FIELD avec allocation de memoire de tous ses attribut
FIELD(const SUPPORT * Support, driverTypes driverType,
- const string & fileName="", const string & fieldName="");
+ const string & fileName="", const string & fieldName="",
+ const int iterationNumber = -1, const int orderNumber = -1)
+ throw (MEDEXCEPTION);
~FIELD();
+ const FIELD operator+(const FIELD& m) const;
+ const FIELD operator-(const FIELD& m) const;
+ const FIELD operator*(const FIELD& m) const;
+ const FIELD operator/(const FIELD& m) const;
+ const FIELD operator-() const;
+ FIELD& operator+=(const FIELD& m);
+ FIELD& operator-=(const FIELD& m);
+ FIELD& operator*=(const FIELD& m);
+ FIELD& operator/=(const FIELD& m);
+ static FIELD* add(const FIELD& m, const FIELD& n);
+ static FIELD* sub(const FIELD& m, const FIELD& n);
+ static FIELD* mul(const FIELD& m, const FIELD& n);
+ static FIELD* div(const FIELD& m, const FIELD& n);
+ double normMax() const throw (MEDEXCEPTION);
+ double norm2() const throw (MEDEXCEPTION);
+ void applyLin(T a, T b);
+ template <T T_function(T)> void applyFunc();
+ static FIELD* scalarProduct(const FIELD& m, const FIELD& n);
+ double normL2(int component, const FIELD<double> * p_field_volume=NULL) const;
+ double normL2(const FIELD<double> * p_field_volume=NULL) const;
+ double normL1(int component, const FIELD<double> * p_field_volume=NULL) const;
+ double normL1(const FIELD<double> * p_field_volume=NULL) const;
+
friend class MED_FIELD_RDONLY_DRIVER<T>;
friend class MED_FIELD_WRONLY_DRIVER<T>;
+
+ friend class VTK_FIELD_DRIVER<T>;
//friend class MED_FIELD_RDWR_DRIVER <T>;
void init ();
inline void write(int index=0, const string & driverName = "");
inline void write(const GENDRIVER &);
+ inline void writeAppend(int index=0, const string & driverName = "");
+ inline void writeAppend(const GENDRIVER &);
+
inline void setValue(MEDARRAY<T> *Value);
inline MEDARRAY<T>* getvalue() const;
*/
void getBarycenter() const throw (MEDEXCEPTION) ;
};
+};
// --------------------
// Implemented Methods
// --------------------
+using namespace MEDMEM;
/*!
- Constructor.
+ Constructor with no parameter, most of the attribut members are set to NULL.
*/
template <class T> FIELD<T>::FIELD():
_value((MEDARRAY<T>*)NULL)
}
/*!
- Constructor.
+ Constructor with parameters such that all attrribut are set but the _value
+ attribut is allocated but not set.
*/
template <class T> FIELD<T>::FIELD(const SUPPORT * Support,
- const int NumberOfComponents):
+ const int NumberOfComponents, const medModeSwitch Mode) throw (MEDEXCEPTION) :
FIELD_(Support, NumberOfComponents)
{
- BEGIN_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)");
+ BEGIN_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
+ SCRUTE(this);
try {
_numberOfValues = Support->getNumberOfElements(MED_ALL_ELEMENTS);
}
MESSAGE("FIELD : constructeur : "<< _numberOfValues <<" et "<< NumberOfComponents);
if (0<_numberOfValues) {
- _value = new MEDARRAY<T>::MEDARRAY(_numberOfComponents,_numberOfValues);
+ _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues,Mode);
_isRead = true ;
} else
_value = (MEDARRAY<T>*)NULL ;
- END_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents)");
+ END_OF("FIELD<T>::FIELD(const SUPPORT * Support, const int NumberOfComponents, const medModeSwitch Mode)");
}
/*!
template <class T> FIELD<T>::FIELD(const FIELD & m):
FIELD_((FIELD_) m)
{
+ MESSAGE("Constructeur FIELD de recopie");
if (m._value != NULL)
{
// copie only default !
- _value = new MEDARRAY<T>::MEDARRAY(* m._value,false);
+ _value = new MEDARRAY<T>(* m._value,false);
}
else
_value = (MEDARRAY<T> *) NULL;
*/
template <class T> FIELD<T> & FIELD<T>::operator=(const FIELD &m)
{
+ MESSAGE("Appel de FIELD<T>::operator=");
+}
+
+/*!
+ Overload addition operator.
+ This operation is authorized only for compatible fields that have the same support.
+ The compatibility checking includes equality tests of the folowing data members:/n
+ - _support
+ - _numberOfComponents
+ - _numberOfValues
+ - _componentsTypes
+ - _MEDComponentsUnits.
+
+ The data members of the returned field are initialized, based on the first field, except for the name,
+ which is the combination of the two field's names and the operator.
+ Advised Utilisation in C++ : <tt> FIELD<T> c = a + b; </tt> /n
+ In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
+ When using python, this operator calls the copy constructor in any case.
+ The user has to be aware that when using operator + in associatives expressions like
+ <tt> a = b + c + d +e; </tt> /n
+ no optimisation is performed : the evaluation of last expression requires the construction of
+ 3 temporary fields.
+*/
+template <class T>
+const FIELD<T> FIELD<T>::operator+(const FIELD & m) const
+{
+ BEGIN_OF("FIELD<T>::operator+(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=this->getvalue()->getMode();
+
+ // Creation of the result - memory is allocated by FIELD constructor
+ FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
+ //result._operation(*this,m,mode,"+"); // perform Atribute's initialization & addition
+ result._operationInitialize(*this,m,"+"); // perform Atribute's initialization
+ result._add_in_place(*this,m,mode); // perform addition
+
+ END_OF("FIELD<T>::operator+(const FIELD & m)");
+ return result;
+}
+
+/*! Overloaded Operator +=
+ * Operations are directly performed in the first field's array.
+ * This operation is authorized only for compatible fields that have the same support.
+ */
+template <class T>
+FIELD<T>& FIELD<T>::operator+=(const FIELD & m)
+{
+ BEGIN_OF("FIELD<T>::operator+=(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // We choose to keep *this mode, even if it may cost a re-calculation for m
+ medModeSwitch mode=this->getvalue()->getMode();
+ const T* value1=m.getValue(mode); // get pointers to the values we are adding
+
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+ const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
+ const T* endV=value+size; // pointer to the end of value
+ for(;value!=endV; value1++,value++)
+ *value += *value1;
+ // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
+ this->getvalue()->clearOtherMode();
+ END_OF("FIELD<T>::operator+=(const FIELD & m)");
+ return *this;
+}
+
+
+/*! Addition of fields. Static member function.
+ * The function return a pointer to a new created field that holds the addition.
+ * Data members are checked for compatibility and initialized.
+ * The user is in charge of memory deallocation.
+ */
+template <class T>
+FIELD<T>* FIELD<T>::add(const FIELD& m, const FIELD& n)
+{
+ BEGIN_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
+ FIELD_::_checkFieldCompatibility(m, n); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=n.getvalue()->getMode();
+
+ // Creation of a new field
+ FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
+ result->_operationInitialize(m,n,"+"); // perform Atribute's initialization
+ result->_add_in_place(m,n,mode); // perform addition
+
+ END_OF("FIELD<T>::add(const FIELD & m, const FIELD& n)");
+ return result;
+}
+
+/*!
+ Overload substraction operator.
+ This operation is authorized only for compatible fields that have the same support.
+ The compatibility checking includes equality tests of the folowing data members:/n
+ - _support
+ - _numberOfComponents
+ - _numberOfValues
+ - _componentsTypes
+ - _MEDComponentsUnits.
+
+ The data members of the returned field are initialized, based on the first field, except for the name,
+ which is the combination of the two field's names and the operator.
+ Advised Utilisation in C++ : <tt> FIELD<T> c = a - b; </tt> /n
+ In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
+ When using python, this operator calls the copy constructor in any case.
+ The user has to be aware that when using operator - in associatives expressions like
+ <tt> a = b - c - d -e; </tt> /n
+ no optimisation is performed : the evaluation of last expression requires the construction of
+ 3 temporary fields.
+*/
+template <class T>
+const FIELD<T> FIELD<T>::operator-(const FIELD & m) const
+{
+ BEGIN_OF("FIELD<T>::operator-(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=this->getvalue()->getMode();
+
+ // Creation of the result - memory is allocated by FIELD constructor
+ FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
+ //result._operation(*this,m,mode,"-"); // perform Atribute's initialization & substraction
+ result._operationInitialize(*this,m,"-"); // perform Atribute's initialization
+ result._sub_in_place(*this,m,mode); // perform substracion
+
+ END_OF("FIELD<T>::operator-(const FIELD & m)");
+ return result;
+}
+
+template <class T>
+const FIELD<T> FIELD<T>::operator-() const
+{
+ BEGIN_OF("FIELD<T>::operator-()");
+
+ medModeSwitch mode=this->getvalue()->getMode();
+ // Creation of the result - memory is allocated by FIELD constructor
+ FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
+ // Atribute's initialization
+ result.setName("- "+getName());
+ result.setComponentsNames(getComponentsNames());
+ // not yet implemented setComponentType(getComponentType());
+ result.setComponentsDescriptions(getComponentsDescriptions());
+ result.setMEDComponentsUnits(getMEDComponentsUnits());
+ result.setComponentsUnits(getComponentsUnits());
+ result.setIterationNumber(getIterationNumber());
+ result.setTime(getTime());
+ result.setOrderNumber(getOrderNumber());
+ result.setValueType(getValueType());
+
+ const T* value1=getValue(mode);
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (result.getValue(mode));
+ const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
+ const T* endV=value+size; // pointer to the end of value
+
+ for(;value!=endV; value1++,value++)
+ *value = -(*value1);
+ END_OF("FIELD<T>::operator-=(const FIELD & m)");
+ return result;
+}
+
+/*! Overloaded Operator -=
+ * Operations are directly performed in the first field's array.
+ * This operation is authorized only for compatible fields that have the same support.
+ */
+template <class T>
+FIELD<T>& FIELD<T>::operator-=(const FIELD & m)
+{
+ BEGIN_OF("FIELD<T>::operator-=(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // We choose to keep *this mode, even if it may cost a re-calculation for m
+ medModeSwitch mode=this->getvalue()->getMode();
+ const T* value1=m.getValue(mode); // get pointers to the values we are adding
+
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+ const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
+ const T* endV=value+size; // pointer to the end of value
+
+ for(;value!=endV; value1++,value++)
+ *value -= *value1;
+ // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
+ this->getvalue()->clearOtherMode();
+ END_OF("FIELD<T>::operator-=(const FIELD & m)");
+ return *this;
+}
+
+
+/*! Substraction of fields. Static member function.
+ * The function return a pointer to a new created field that holds the substraction.
+ * Data members are checked for compatibility and initialized.
+ * The user is in charge of memory deallocation.
+ */
+template <class T>
+FIELD<T>* FIELD<T>::sub(const FIELD& m, const FIELD& n)
+{
+ BEGIN_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
+ FIELD_::_checkFieldCompatibility(m, n); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=n.getvalue()->getMode();
+
+ // Creation of a new field
+ FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
+ result->_operationInitialize(m,n,"-"); // perform Atribute's initialization
+ result->_sub_in_place(m,n,mode); // perform substraction
+
+ END_OF("FIELD<T>::sub(const FIELD & m, const FIELD& n)");
+ return result;
+}
+
+/*!
+ Overload multiplication operator.
+ This operation is authorized only for compatible fields that have the same support.
+ The compatibility checking includes equality tests of the folowing data members:/n
+ - _support
+ - _numberOfComponents
+ - _numberOfValues
+ - _componentsTypes
+ - _MEDComponentsUnits.
+
+ The data members of the returned field are initialized, based on the first field, except for the name,
+ which is the combination of the two field's names and the operator.
+ Advised Utilisation in C++ : <tt> FIELD<T> c = a * b; </tt> /n
+ In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
+ When using python, this operator calls the copy constructor in any case.
+ The user has to be aware that when using operator * in associatives expressions like
+ <tt> a = b * c * d *e; </tt> /n
+ no optimisation is performed : the evaluation of last expression requires the construction of
+ 3 temporary fields.
+*/
+template <class T>
+const FIELD<T> FIELD<T>::operator*(const FIELD & m) const
+{
+ BEGIN_OF("FIELD<T>::operator*(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=this->getvalue()->getMode();
+
+ // Creation of the result - memory is allocated by FIELD constructor
+ FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
+ //result._operation(*this,m,mode,"*"); // perform Atribute's initialization & multiplication
+ result._operationInitialize(*this,m,"*"); // perform Atribute's initialization
+ result._mul_in_place(*this,m,mode); // perform multiplication
+
+ END_OF("FIELD<T>::operator*(const FIELD & m)");
+ return result;
+}
+
+/*! Overloaded Operator *=
+ * Operations are directly performed in the first field's array.
+ * This operation is authorized only for compatible fields that have the same support.
+ */
+template <class T>
+FIELD<T>& FIELD<T>::operator*=(const FIELD & m)
+{
+ BEGIN_OF("FIELD<T>::operator*=(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // We choose to keep *this mode, even if it may cost a re-calculation for m
+ medModeSwitch mode=this->getvalue()->getMode();
+ const T* value1=m.getValue(mode); // get pointers to the values we are adding
+
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+ const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
+ const T* endV=value+size; // pointer to the end of value
+
+ for(;value!=endV; value1++,value++)
+ *value *= *value1;
+ // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
+ this->getvalue()->clearOtherMode();
+ END_OF("FIELD<T>::operator*=(const FIELD & m)");
+ return *this;
+}
+
+
+/*! Multiplication of fields. Static member function.
+ * The function return a pointer to a new created field that holds the multiplication.
+ * Data members are checked for compatibility and initialized.
+ * The user is in charge of memory deallocation.
+ */
+template <class T>
+FIELD<T>* FIELD<T>::mul(const FIELD& m, const FIELD& n)
+{
+ BEGIN_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
+ FIELD_::_checkFieldCompatibility(m, n); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=n.getvalue()->getMode();
+
+ // Creation of a new field
+ FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
+ result->_operationInitialize(m,n,"*"); // perform Atribute's initialization
+ result->_mul_in_place(m,n,mode); // perform multiplication
+
+ END_OF("FIELD<T>::mul(const FIELD & m, const FIELD& n)");
+ return result;
+}
+
+
+/*!
+ Overload division operator.
+ This operation is authorized only for compatible fields that have the same support.
+ The compatibility checking includes equality tests of the folowing data members:/n
+ - _support
+ - _numberOfComponents
+ - _numberOfValues
+ - _componentsTypes
+ - _MEDComponentsUnits.
+
+ The data members of the returned field are initialized, based on the first field, except for the name,
+ which is the combination of the two field's names and the operator.
+ Advised Utilisation in C++ : <tt> FIELD<T> c = a / b; </tt> /n
+ In this case, the (recent) compilators perform optimisation and don't call the copy constructor.
+ When using python, this operator calls the copy constructor in any case.
+ The user has to be aware that when using operator / in associatives expressions like
+ <tt> a = b / c / d /e; </tt> /n
+ no optimisation is performed : the evaluation of last expression requires the construction of
+ 3 temporary fields.
+*/
+template <class T>
+const FIELD<T> FIELD<T>::operator/(const FIELD & m) const
+{
+ BEGIN_OF("FIELD<T>::operator/(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(this->getvalue()->getMode()==m.getvalue()->getMode() || this->getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=this->getvalue()->getMode();
+
+ // Creation of the result - memory is allocated by FIELD constructor
+ FIELD<T> result(this->getSupport(),this->getNumberOfComponents(),mode);
+ //result._operation(*this,m,mode,"/"); // perform Atribute's initialization & division
+ result._operationInitialize(*this,m,"/"); // perform Atribute's initialization
+ result._div_in_place(*this,m,mode); // perform division
+
+ END_OF("FIELD<T>::operator/(const FIELD & m)");
+ return result;
+}
+
+
+/*! Overloaded Operator /=
+ * Operations are directly performed in the first field's array.
+ * This operation is authorized only for compatible fields that have the same support.
+ */
+template <class T>
+FIELD<T>& FIELD<T>::operator/=(const FIELD & m)
+{
+ BEGIN_OF("FIELD<T>::operator/=(const FIELD & m)");
+ FIELD_::_checkFieldCompatibility(*this, m); // may throw exception
+
+ // We choose to keep *this mode, even if it may cost a re-calculation for m
+ medModeSwitch mode=this->getvalue()->getMode();
+ const T* value1=m.getValue(mode); // get pointers to the values we are adding
+
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+ const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
+ const T* endV=value+size; // pointer to the end of value
+
+ for(;value!=endV; value1++,value++)
+ *value /= *value1;
+ // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
+ this->getvalue()->clearOtherMode();
+ END_OF("FIELD<T>::operator/=(const FIELD & m)");
+ return *this;
+}
+
+
+/*! Division of fields. Static member function.
+ * The function return a pointer to a new created field that holds the division.
+ * Data members are checked for compatibility and initialized.
+ * The user is in charge of memory deallocation.
+ */
+template <class T>
+FIELD<T>* FIELD<T>::div(const FIELD& m, const FIELD& n)
+{
+ BEGIN_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
+ FIELD_::_checkFieldCompatibility(m, n); // may throw exception
+
+ // Select mode : avoid if possible any calculation of other mode for fields m or *this
+ medModeSwitch mode;
+ if(m.getvalue()->getMode()==n.getvalue()->getMode() || n.getvalue()->isOtherCalculated())
+ mode=m.getvalue()->getMode();
+ else
+ mode=n.getvalue()->getMode();
+
+ // Creation of a new field
+ FIELD<T>* result = new FIELD<T>(m.getSupport(),m.getNumberOfComponents(),mode);
+ result->_operationInitialize(m,n,"/"); // perform Atribute's initialization
+ result->_div_in_place(m,n,mode); // perform division
+
+ END_OF("FIELD<T>::div(const FIELD & m, const FIELD& n)");
+ return result;
+}
+
+
+/*!
+ \if developper
+ This internal method initialize the members of a new field created to hold the result of the operation Op .
+ Initialization is based on the first field, except for the name, which is the combination of the two field's names
+ and the operator.
+ \endif
+*/
+template <class T>
+void FIELD<T>::_operationInitialize(const FIELD& m,const FIELD& n, char* Op)
+{
+ MESSAGE("Appel methode interne _add" << Op);
+
+ // Atribute's initialization (copy of the first field's attributes)
+ // Other data members (_support, _numberOfValues) are initialized in the field's constr.
+ setName(m.getName()+" "+Op+" "+n.getName());
+ setComponentsNames(m.getComponentsNames());
+ // not yet implemented setComponentType(m.getComponentType());
+ setComponentsDescriptions(m.getComponentsDescriptions());
+ setMEDComponentsUnits(m.getMEDComponentsUnits());
+
+ // The following data member may differ from field m to n.
+ // The initialization is done based on the first field.
+ setComponentsUnits(m.getComponentsUnits());
+ setIterationNumber(m.getIterationNumber());
+ setTime(m.getTime());
+ setOrderNumber(m.getOrderNumber());
+ setValueType(m.getValueType());
+}
+
+
+/*!
+ \if developper
+ Internal method called by FIELD<T>::operator+ and FIELD<T>::add to perform addition "in place".
+ This method is applied to a just created field with medModeSwitch mode.
+ For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
+ it doesn't exist!
+ \endif
+*/
+template <class T>
+void FIELD<T>::_add_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
+{
+ // get pointers to the values we are adding
+ const T* value1=m.getValue(mode);
+ const T* value2=n.getValue(mode);
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+
+ const int size=getNumberOfValues()*getNumberOfComponents();
+ SCRUTE(size);
+ const T* endV1=value1+size;
+ for(;value1!=endV1; value1++,value2++,value++)
+ *value=(*value1)+(*value2);
+}
+
+/*!
+ \if developper
+ Internal method called by FIELD<T>::operator- and FIELD<T>::sub to perform substraction "in place".
+ This method is applied to a just created field with medModeSwitch mode.
+ For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
+ it doesn't exist!
+ \endif
+*/
+template <class T>
+void FIELD<T>::_sub_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
+{
+ // get pointers to the values we are adding
+ const T* value1=m.getValue(mode);
+ const T* value2=n.getValue(mode);
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+
+ const int size=getNumberOfValues()*getNumberOfComponents();
+ SCRUTE(size);
+ const T* endV1=value1+size;
+ for(;value1!=endV1; value1++,value2++,value++)
+ *value=(*value1)-(*value2);
+}
+
+/*!
+ \if developper
+ Internal method called by FIELD<T>::operator* and FIELD<T>::mul to perform multiplication "in place".
+ This method is applied to a just created field with medModeSwitch mode.
+ For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
+ it doesn't exist!
+ \endif
+*/
+template <class T>
+void FIELD<T>::_mul_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
+{
+ // get pointers to the values we are adding
+ const T* value1=m.getValue(mode);
+ const T* value2=n.getValue(mode);
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+
+ const int size=getNumberOfValues()*getNumberOfComponents();
+ SCRUTE(size);
+ const T* endV1=value1+size;
+ for(;value1!=endV1; value1++,value2++,value++)
+ *value=(*value1)*(*value2);
+}
+
+/*!
+ \if developper
+ Internal method called by FIELD<T>::operator/ and FIELD<T>::div to perform division "in place".
+ This method is applied to a just created field with medModeSwitch mode.
+ For this reason, the alternate mode doesn't need to be set to 0 after performing operation :
+ it doesn't exist!
+ \endif
+*/
+template <class T>
+void FIELD<T>::_div_in_place(const FIELD& m,const FIELD& n, const medModeSwitch mode)
+{
+ // get pointers to the values we are adding
+ const T* value1=m.getValue(mode);
+ const T* value2=n.getValue(mode);
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+
+ const int size=getNumberOfValues()*getNumberOfComponents();
+ SCRUTE(size);
+ const T* endV1=value1+size;
+ for(;value1!=endV1; value1++,value2++,value++)
+ *value=(*value1)/(*value2);
+}
+
+/*! Return Max Norm
+ */
+template <class T> double FIELD<T>::normMax() const throw (MEDEXCEPTION)
+{
+ const T* value=getValue(getvalue()->getMode()); // get pointer to the values
+ const int size=getNumberOfValues()*getNumberOfComponents();
+ if (size <= 0) // Size of array has to be strictly positive
+ {
+ string diagnosis;
+ diagnosis="FIELD<T>::normMax() : cannot compute the norm of "+getName()+
+ " : it size is non positive!";
+ throw MEDEXCEPTION(diagnosis.c_str());
+ }
+ const T* lastvalue=value+size; // get pointer just after last value
+ const T* pMax=value; // pointer to the max value
+ const T* pMin=value; // pointer to the min value
+
+ // get pointers to the max & min value of array
+ while ( ++value != lastvalue )
+ {
+ if ( *pMin > *value )
+ pMin=value;
+ if ( *pMax < *value )
+ pMax=value;
+ }
+
+ T Max= *pMax>(T) 0 ? *pMax : -*pMax; // Max=abs(*pMax)
+ T Min= *pMin>(T) 0 ? *pMin : -*pMin; // Min=abs(*pMin)
+
+ return Max>Min ? static_cast<double>(Max) : static_cast<double>(Min);
+}
+
+/*! Return Euclidien norm
+ */
+template <class T> double FIELD<T>::norm2() const throw (MEDEXCEPTION)
+{
+ const T* value=this->getValue(this->getvalue()->getMode()); // get const pointer to the values
+ const int size=getNumberOfValues()*getNumberOfComponents(); // get size of array
+ if (size <= 0) // Size of array has to be strictly positive
+ {
+ string diagnosis;
+ diagnosis="FIELD<T>::norm2() : cannot compute the norm of "+getName()+
+ " : it size is non positive!";
+ throw MEDEXCEPTION(diagnosis.c_str());
+ }
+ const T* lastvalue=value+size; // point just after last value
+
+ T result((T)0); // init
+ for( ; value!=lastvalue ; ++value)
+ result += (*value) * (*value);
+
+ return std::sqrt(static_cast<double> (result));
+}
+
+
+/*! Apply to each (scalar) field component the template parameter T_function,
+ * which is a pointer to function.
+ * Since the pointer is known at compile time, the function is inlined into the inner loop!
+ * calculation is done "in place".
+ * Use examples :
+ *
+ * \code myField.applyFunc<std::sqrt>(); // apply sqare root function \endcode
+ * \code myField.applyFunc<myFunction>(); // apply your own created function \endcode
+ */
+template <class T> template <T T_function(T)>
+void FIELD<T>::applyFunc()
+{
+ medModeSwitch mode=getvalue()->getMode();
+
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (getValue(mode));
+ const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
+
+ if (size>0) // for a negative size, there is nothing to do
+ {
+ const T* lastvalue=value+size; // pointer to the end of value
+ for(;value!=lastvalue; ++value) // apply linear transformation
+ *value = T_function(*value);
+ // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
+ getvalue()->clearOtherMode();
+ }
}
+
+
+/*! Apply to each (scalar) field component the linear function x -> ax+b.
+ * calculation is done "in place".
+ */
+template <class T> void FIELD<T>::applyLin(T a, T b)
+{
+ medModeSwitch mode=getvalue()->getMode();
+
+ // get a non const pointer to the inside array of values and perform operation in place
+ T * value=const_cast<T *> (getValue(mode));
+ const int size=getNumberOfValues()*getNumberOfComponents(); // size of array
+
+ if (size>0) // for a negative size, there is nothing to do
+ {
+ const T* lastvalue=value+size; // pointer to the end of value
+ for(;value!=lastvalue; ++value) // apply linear transformation
+ *value = a*(*value)+b;
+ // if it exists, the alternate mode needs to be cleared because it is not up-to-date anymore
+ getvalue()->clearOtherMode();
+ }
+}
+
/*!
- Constructor.
+ * Return a pointer to a new field that holds the scalar product. Static member function.
+ * This operation is authorized only for compatible fields that have the same support.
+ * The compatibility checking includes equality tests of the folowing data members:/n
+ * - _support
+ * - _numberOfComponents
+ * - _numberOfValues
+ * - _componentsTypes
+ * - _MEDComponentsUnits.
+ * Data members are initialized.
+ * The new field point to the same support and has one component.
+ * Each value of it is the scalar product of the two argument's fields.
+ * The user is in charge of memory deallocation.
+ */
+template <class T> FIELD<T>* FIELD<T>::scalarProduct(const FIELD & m, const FIELD & n)
+{
+ FIELD_::_checkFieldCompatibility( m, n); // may throw exception
+ // we need a MED_FULL_INTERLACE representation of m & n to compute the scalar product
+ const medModeSwitch mode=MED_FULL_INTERLACE;
+
+ const int numberOfElements=m.getNumberOfValues(); // strictly positive
+ const int NumberOfComponents=m.getNumberOfComponents(); // strictly positive
+
+ // Creation & init of a the result field on the same suppot, with one component
+ FIELD<T>* result = new FIELD<T>(m.getSupport(),1,mode);
+ result->setName( "scalarProduct ( " + m.getName() + " , " + n.getName() + " )" );
+ result->setIterationNumber(m.getIterationNumber());
+ result->setTime(m.getTime());
+ result->setOrderNumber(m.getOrderNumber());
+ result->setValueType(m.getValueType());
+
+ const T* value1=m.getValue(mode); // get const pointer to the values
+ const T* value2=n.getValue(mode); // get const pointer to the values
+ // get a non const pointer to the inside array of values and perform operation
+ T * value=const_cast<T *> (result->getValue(mode));
+
+ const T* lastvalue=value+numberOfElements; // pointing just after last value of result
+ for ( ; value!=lastvalue ; ++value ) // loop on all elements
+ {
+ *value=(T)0; // initialize value
+ const T* endofRow=value1+NumberOfComponents; // pointing just after end of row
+ for ( ; value1 != endofRow; ++value1, ++value2) // computation of dot product
+ *value += (*value1) * (*value2);
+ }
+ return result;
+}
+
+/*! Return L2 Norm of the field's component.
+ * Cannot be applied to a field with a support on nodes.
+ * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
+ */
+template <class T> double FIELD<T>::normL2(int component, const FIELD<double> * p_field_volume) const
+{
+ _checkNormCompatibility(p_field_volume); // may throw exception
+ if ( component<1 || component>getNumberOfComponents() )
+ throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
+
+ const FIELD<double> * p_field_size=p_field_volume;
+ if(!p_field_volume) // if the user don't supply the volume
+ p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
+
+ // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
+ const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
+ const T* value=getValueI( MED_NO_INTERLACE, component); // get pointer to the component's values
+ const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
+
+ double integrale=0.0;
+ double totVol=0.0;
+ for (; value!=lastvalue ; ++value ,++vol)
+ {
+ integrale += static_cast<double>((*value) * (*value)) * (*vol);
+ totVol+=*vol;
+ }
+
+ if(!p_field_volume) // if the user didn't supply the volume
+ delete p_field_size; // delete temporary volume field
+ if( totVol <= 0)
+ throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
+
+ return integrale/totVol;
+}
+
+/*! Return L2 Norm of the field.
+ * Cannot be applied to a field with a support on nodes.
+ * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
+ */
+template <class T> double FIELD<T>::normL2(const FIELD<double> * p_field_volume) const
+{
+ _checkNormCompatibility(p_field_volume); // may throw exception
+ const FIELD<double> * p_field_size=p_field_volume;
+ if(!p_field_volume) // if the user don't supply the volume
+ p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
+
+ // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
+ const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
+ const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
+ const T* value=getValue( MED_NO_INTERLACE); // get pointer to the field's values
+
+ double totVol=0.0;
+ const double* p_vol=vol;
+ for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
+ totVol+=*p_vol;
+
+ double integrale=0.0;
+ for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
+ for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
+ integrale += static_cast<double>((*value) * (*value)) * (*p_vol);
+
+ if(!p_field_volume) // if the user didn't supply the volume
+ delete p_field_size; // delete temporary volume field
+ if( totVol <= 0)
+ throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
+
+ return integrale/totVol;
+}
+
+/*! Return L1 Norm of the field's component.
+ * Cannot be applied to a field with a support on nodes.
+ * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
+ */
+template <class T> double FIELD<T>::normL1(int component, const FIELD<double> * p_field_volume) const
+{
+ _checkNormCompatibility(p_field_volume); // may throw exception
+ if ( component<1 || component>getNumberOfComponents() )
+ throw MEDEXCEPTION(STRING("FIELD<T>::normL2() : The component argument should be between 1 and the number of components"));
+
+ const FIELD<double> * p_field_size=p_field_volume;
+ if(!p_field_volume) // if the user don't supply the volume
+ p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
+
+ // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
+ const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
+ const T* value=getValueI( MED_NO_INTERLACE, component); // get pointer to the component's values
+ const T* lastvalue=value+getNumberOfValues(); // pointing just after the end of column
+
+ double integrale=0.0;
+ double totVol=0.0;
+ for (; value!=lastvalue ; ++value ,++vol)
+ {
+ integrale += std::abs( static_cast<double>(*value) ) * (*vol);
+ totVol+=*vol;
+ }
+
+ if(!p_field_volume) // if the user didn't supply the volume
+ delete p_field_size; // delete temporary volume field
+ if( totVol <= 0)
+ throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
+
+ return integrale/totVol;
+}
+
+/*! Return L1 Norm of the field.
+ * Cannot be applied to a field with a support on nodes.
+ * If the optional p_field_volume argument is furnished, the volume is not re-calculated.
+ */
+template <class T> double FIELD<T>::normL1(const FIELD<double> * p_field_volume) const
+{
+ _checkNormCompatibility(p_field_volume); // may throw exception
+ const FIELD<double> * p_field_size=p_field_volume;
+ if(!p_field_volume) // if the user don't supply the volume
+ p_field_size=_getFieldSize(); // we calculate the volume [PROVISOIRE, en attendant l'implémentation dans mesh]
+
+ // get pointer to the element's volumes. MED_FULL_INTERLACE is the default mode for p_field_size
+ const double* vol=p_field_size->getValue(MED_FULL_INTERLACE);
+ const double* lastvol=vol+getNumberOfValues(); // pointing just after the end of vol
+ const T* value=getValue( MED_NO_INTERLACE); // get pointer to the field's values
+
+ double totVol=0.0;
+ const double* p_vol=vol;
+ for (p_vol=vol; p_vol!=lastvol ; ++p_vol) // calculate total volume
+ totVol+=*p_vol;
+
+ double integrale=0.0;
+ for (int i=1; i<=getNumberOfComponents(); ++i) // compute integral on all components
+ for (p_vol=vol; p_vol!=lastvol ; ++value ,++p_vol)
+ integrale += std::abs( static_cast<double>(*value) ) * (*p_vol);
+
+ if(!p_field_volume) // if the user didn't supply the volume
+ delete p_field_size; // delete temporary volume field
+ if( totVol <= 0)
+ throw MEDEXCEPTION(STRING("cannot compute sobolev norm : volume is not positive!"));
+
+ return integrale/totVol;
+}
+
+
+
+
+/*!
+ Constructor with parameters; the object is set via a file and its associated
+ driver. For the moment only the MED_DRIVER is considered and if the last two
+ argument (iterationNumber and orderNumber) are not set; their default value
+ is -1. If the field fieldDriverName with the iteration number
+ iterationNumber and the order number orderNumber does not exist in the file
+ fieldDriverName; the constructor raises an exception.
*/
template <class T> FIELD<T>::FIELD(const SUPPORT * Support,
driverTypes driverType,
const string & fileName/*=""*/,
- const string & fieldDriverName/*=""*/)
+ const string & fieldDriverName/*=""*/,
+ const int iterationNumber,
+ const int orderNumber)
+ throw (MEDEXCEPTION)
{
- const char * LOC = "template <class T> FIELD<T>::FIELD(const SUPPORT * Support, driverTypes driverType, const string & fileName=\"\", const string & fieldName=\"\") : ";
+ 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) : ";
int current;
_support = Support;
_value = (MEDARRAY<T>*)NULL;
- MED_FIELD_RDONLY_DRIVER<T> myDriver(fileName,this);
- myDriver.setFieldName(fieldDriverName);
- current = addDriver(myDriver);
+ _iterationNumber = iterationNumber;
+ _time = 0.0;
+ _orderNumber = orderNumber;
+
+ switch(driverType)
+ {
+ case MED_DRIVER :
+ {
+ MED_FIELD_RDONLY_DRIVER<T> myDriver(fileName,this);
+ myDriver.setFieldName(fieldDriverName);
+ current = addDriver(myDriver);
+ break;
+ }
// current = addDriver(driverType,fileName,fieldDriverName);
// switch(_drivers[current]->getAccessMode() ) {
// case MED_WRONLY : {
// default : {
// }
// }
+ default :
+ {
+ throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Driver type unknown : can't create driver!"));
+ break;
+ }
+ }
+
_drivers[current]->open();
_drivers[current]->read();
_drivers[current]->close();
template <class T> FIELD<T>::~FIELD()
{
BEGIN_OF(" Destructeur FIELD<T>::~FIELD()");
+ SCRUTE(this);
if (_value) delete _value;
END_OF(" Destructeur FIELD<T>::~FIELD()");
}
_numberOfValues = _support->getNumberOfElements(MED_ALL_ELEMENTS);
MESSAGE(LOC <<" : "<<_numberOfValues <<" et "<< NumberOfComponents);
- _value = new MEDARRAY<T>::MEDARRAY(_numberOfComponents,_numberOfValues);
+ _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
_isRead = true ;
}
MESSAGE("FIELD : constructeur : "<<LengthValue <<" et "<< NumberOfComponents);
_numberOfValues = LengthValue ;
- _value = new MEDARRAY<T>::MEDARRAY(_numberOfComponents,_numberOfValues);
+ _value = new MEDARRAY<T>(_numberOfComponents,_numberOfValues);
_isRead = true ;
SCRUTE(_value);
template <class T> FIELD<T>::INSTANCE_DE<MED_FIELD_RDWR_DRIVER<T> > FIELD<T>::inst_med ;
-template <class T> const FIELD<T>::INSTANCE * const FIELD<T>::instances[] = { &FIELD<T>::inst_med } ;
+
+template <class T> FIELD<T>::INSTANCE_DE<VTK_FIELD_DRIVER<T> > FIELD<T>::inst_vtk ;
+
+template <class T> const typename FIELD<T>::INSTANCE * const FIELD<T>::instances[] = { &FIELD<T>::inst_med, &FIELD<T>::inst_vtk} ;
/*!
Create the specified driver and return its index reference to path to
read or write methods.
*/
+
template <class T> int FIELD<T>::addDriver(driverTypes driverType,
const string & fileName/*="Default File Name.med"*/,
const string & driverName/*="Default Field Name"*/)
GENDRIVER * driver;
int current;
+ int itDriver = (int) NO_DRIVER;
BEGIN_OF(LOC);
- driver = instances[driverType]->run(fileName, this) ;
+ SCRUTE(driverType);
+
+ SCRUTE(instances[driverType]);
+
+ switch(driverType)
+ {
+ case MED_DRIVER : {
+ itDriver = (int) driverType ;
+ break ;
+ }
+
+ case VTK_DRIVER : {
+ itDriver = 1 ;
+ break ;
+ }
+
+ case GIBI_DRIVER : {
+ throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "driverType other than MED_DRIVER and VTK_DRIVER has been specified to the method which is not allowed for the object FIELD"));
+ break;
+ }
+
+ case NO_DRIVER : {
+ throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "driverType other than MED_DRIVER and VTK_DRIVER has been specified to the method which is not allowed for the object FIELD"));
+ break;
+ }
+ }
+
+ if (itDriver == ((int) NO_DRIVER))
+ throw MED_EXCEPTION (LOCALIZED(STRING(LOC)<< "driverType other than MED_DRIVER and VTK_DRIVER has been specified to the method which is not allowed for the object FIELD"));
+
+ driver = instances[itDriver]->run(fileName, this) ;
+
_drivers.push_back(driver);
+
current = _drivers.size()-1;
_drivers[current]->setFieldName(driverName);
- return current;
END_OF(LOC);
+ return current;
}
_drivers.push_back(newDriver);
return _drivers.size() -1 ;
-
+
END_OF(LOC);
};
}
else
throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
- << "The index given is invalid, index must be between 0 and |"
+ << "The <index given is invalid, index must be between 0 and |"
<< _drivers.size()
)
);
END_OF(LOC);
}
+/*!
+ Write FIELD in the file specified in the driver given by its index. Use this
+ method for ASCII drivers (e.g. VTK_DRIVER)
+*/
+template <class T> inline void FIELD<T>::writeAppend(int index/*=0*/, const string & driverName /*= ""*/)
+{
+ const char * LOC = "FIELD<T>::write(int index=0, const string & driverName = \"\") : ";
+ BEGIN_OF(LOC);
+
+ if( _drivers[index] ) {
+ _drivers[index]->openAppend();
+ if (driverName != "") _drivers[index]->setFieldName(driverName);
+ _drivers[index]->writeAppend();
+ _drivers[index]->close();
+ }
+ else
+ throw MED_EXCEPTION ( LOCALIZED( STRING(LOC)
+ << "The index given is invalid, index must be between 0 and |"
+ << _drivers.size()
+ )
+ );
+ END_OF(LOC);
+}
+
/*!
\internal
Write FIELD with the driver which is equal to the given driver.
const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
BEGIN_OF(LOC);
- for (int index=0; index < _drivers.size(); index++ )
+ for (unsigned int index=0; index < _drivers.size(); index++ )
if ( *_drivers[index] == genDriver ) {
_drivers[index]->open();
_drivers[index]->write();
}
+/*!
+ \internal
+ Write FIELD with the driver which is equal to the given driver.
+
+ Use by MED object. Use this method for ASCII drivers (e.g. VTK_DRIVER).
+*/
+template <class T> inline void FIELD<T>::writeAppend(const GENDRIVER & genDriver)
+{
+ const char * LOC = " FIELD<T>::write(const GENDRIVER &) : ";
+ BEGIN_OF(LOC);
+
+ for (unsigned int index=0; index < _drivers.size(); index++ )
+ if ( *_drivers[index] == genDriver ) {
+ _drivers[index]->openAppend();
+ _drivers[index]->writeAppend();
+ _drivers[index]->close();
+ }
+
+ END_OF(LOC);
+
+}
+
/*!
\internal
Read FIELD with the driver which is equal to the given driver.
const char * LOC = " FIELD<T>::read(const GENDRIVER &) : ";
BEGIN_OF(LOC);
- for (int index=0; index < _drivers.size(); index++ )
+ for (unsigned int index=0; index < _drivers.size(); index++ )
if ( *_drivers[index] == genDriver ) {
_drivers[index]->open();
_drivers[index]->read();
