Salome HOME
API modification : input string have type :
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingFieldDouble.cxx
1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay (CEA/DEN)
20
21 #include "MEDCouplingFieldDouble.hxx"
22 #include "MEDCouplingFieldTemplate.hxx"
23 #include "MEDCouplingUMesh.hxx"
24 #include "MEDCouplingTimeDiscretization.hxx"
25 #include "MEDCouplingFieldDiscretization.hxx"
26 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
27 #include "MEDCouplingNatureOfField.hxx"
28
29 #include "InterpKernelAutoPtr.hxx"
30
31 #include <sstream>
32 #include <limits>
33 #include <algorithm>
34 #include <functional>
35
36 using namespace ParaMEDMEM;
37
38
39 /*!
40  * Creates a new MEDCouplingFieldDouble, of given spatial type and time discretization.
41  * For more info, see \ref MEDCouplingFirstSteps3.
42  * \param [in] type - the type of spatial discretization of the created field, one of
43  *        (\ref ParaMEDMEM::ON_CELLS "ON_CELLS", 
44  *         \ref ParaMEDMEM::ON_NODES "ON_NODES",
45  *         \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT", 
46  *         \ref ParaMEDMEM::ON_GAUSS_NE "ON_GAUSS_NE",
47  *         \ref ParaMEDMEM::ON_NODES_KR "ON_NODES_KR").
48  * \param [in] td - the type of time discretization of the created field, one of
49  *        (\ref ParaMEDMEM::NO_TIME "NO_TIME", 
50  *         \ref ParaMEDMEM::ONE_TIME "ONE_TIME", 
51  *         \ref ParaMEDMEM::LINEAR_TIME "LINEAR_TIME", 
52  *         \ref ParaMEDMEM::CONST_ON_TIME_INTERVAL "CONST_ON_TIME_INTERVAL").
53  * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
54  *         caller is to delete this field using decrRef() as it is no more needed. 
55  */
56 MEDCouplingFieldDouble* MEDCouplingFieldDouble::New(TypeOfField type, TypeOfTimeDiscretization td)
57 {
58   return new MEDCouplingFieldDouble(type,td);
59 }
60
61 /*!
62  * Creates a new MEDCouplingFieldDouble, of a given time discretization and with a
63  * spatial type and supporting mesh copied from a given 
64  * \ref MEDCouplingFieldTemplatesPage "field template".
65  * For more info, see \ref MEDCouplingFirstSteps3.
66  * \warning This method does not deeply copy neither the mesh nor the spatial
67  * discretization. Only a shallow copy (reference) is done for the mesh and the spatial
68  * discretization!
69  * \param [in] ft - the \ref MEDCouplingFieldTemplatesPage "field template" defining
70  *        the spatial discretization and the supporting mesh.
71  * \param [in] td - the type of time discretization of the created field, one of
72  *        (\ref ParaMEDMEM::NO_TIME "NO_TIME", 
73  *         \ref ParaMEDMEM::ONE_TIME "ONE_TIME", 
74  *         \ref ParaMEDMEM::LINEAR_TIME "LINEAR_TIME", 
75  *         \ref ParaMEDMEM::CONST_ON_TIME_INTERVAL "CONST_ON_TIME_INTERVAL").
76  * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
77  *         caller is to delete this field using decrRef() as it is no more needed. 
78  */
79 MEDCouplingFieldDouble *MEDCouplingFieldDouble::New(const MEDCouplingFieldTemplate& ft, TypeOfTimeDiscretization td)
80 {
81   return new MEDCouplingFieldDouble(ft,td);
82 }
83
84 /*!
85  * Sets a time \a unit of \a this field. For more info, see \ref MEDCouplingFirstSteps3.
86  * \param [in] unit \a unit (string) in which time is measured.
87  */
88 void MEDCouplingFieldDouble::setTimeUnit(const std::string& unit)
89 {
90   _time_discr->setTimeUnit(unit);
91 }
92
93 /*!
94  * Returns a time unit of \a this field.
95  * \return a string describing units in which time is measured.
96  */
97 std::string MEDCouplingFieldDouble::getTimeUnit() const
98 {
99   return _time_discr->getTimeUnit();
100 }
101
102 /*!
103  * This method if possible the time information (time unit, time iteration, time unit and time value) with its support
104  * that is to say its mesh.
105  * 
106  * \throw  If \c this->_mesh is null an exception will be thrown. An exception will also be throw if the spatial discretization is
107  *         NO_TIME.
108  */
109 void MEDCouplingFieldDouble::synchronizeTimeWithSupport()
110 {
111   _time_discr->synchronizeTimeWith(_mesh);
112 }
113
114 /*!
115  * Returns a new MEDCouplingFieldDouble which is a copy of \a this one. The data
116  * of \a this field is copied either deep or shallow depending on \a recDeepCpy
117  * parameter. But the underlying mesh is always shallow copied.
118  * Data that can be copied either deeply or shallow are:
119  * - \ref MEDCouplingTemporalDisc "temporal discretization" data that holds array(s)
120  * of field values,
121  * - \ref MEDCouplingSpatialDisc "a spatial discretization".
122  * 
123  * \c clone(false) is rather dedicated for advanced users that want to limit the amount
124  * of memory. It allows the user to perform methods like operator+(), operator*()
125  * etc. with \a this and the returned field. If the user wants to duplicate deeply the
126  * underlying mesh he should call cloneWithMesh() method or deepCpy() instead. 
127  * \warning The underlying \b mesh of the returned field is **always the same**
128  *         (pointer) as \a this one **whatever the value** of \a recDeepCpy parameter.
129  *  \param [in] recDeepCpy - if \c true, the copy of the underlying data arrays is
130  *         deep, else all data arrays of \a this field are shared by the new field.
131  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
132  *         caller is to delete this field using decrRef() as it is no more needed.
133  * \sa cloneWithMesh()
134  */
135 MEDCouplingFieldDouble *MEDCouplingFieldDouble::clone(bool recDeepCpy) const
136 {
137   return new MEDCouplingFieldDouble(*this,recDeepCpy);
138 }
139
140 /*!
141  * Returns a new MEDCouplingFieldDouble which is a copy of \a this one. The data
142  * of \a this field is copied either deep or shallow depending on \a recDeepCpy
143  * parameter. But the underlying mesh is always deep copied.
144  * Data that can be copied either deeply or shallow are:
145  * - \ref MEDCouplingTemporalDisc "temporal discretization" data that holds array(s)
146  * of field values,
147  * - \ref MEDCouplingSpatialDisc "a spatial discretization".
148  * 
149  * This method behaves exactly like clone() except that here the underlying **mesh is
150  * always deeply duplicated**, whatever the value \a recDeepCpy parameter.
151  * The result of \c cloneWithMesh(true) is exactly the same as that of deepCpy().
152  * So the resulting field can not be used together with \a this one in the methods
153  * like operator+(), operator*() etc. To avoid deep copying the underlying mesh,
154  * the user can call clone().
155  *  \param [in] recDeepCpy - if \c true, the copy of the underlying data arrays is
156  *         deep, else all data arrays of \a this field are shared by the new field.
157  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
158  *         caller is to delete this field using decrRef() as it is no more needed.
159  * \sa clone()
160  */
161 MEDCouplingFieldDouble *MEDCouplingFieldDouble::cloneWithMesh(bool recDeepCpy) const
162 {
163   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(recDeepCpy);
164   if(_mesh)
165     {
166       MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> mCpy=_mesh->deepCpy();
167       ret->setMesh(mCpy);
168     }
169   return ret.retn();
170 }
171
172 /*!
173  * Returns a new MEDCouplingFieldDouble which is a deep copy of \a this one **including
174  * the mesh**.
175  * The result of this method is exactly the same as that of \c cloneWithMesh(true).
176  * So the resulting field can not be used together with \a this one in the methods
177  * like operator+(), operator*() etc. To avoid deep copying the underlying mesh,
178  * the user can call clone().
179  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
180  *         caller is to delete this field using decrRef() as it is no more needed.
181  * \sa cloneWithMesh()
182  */
183 MEDCouplingFieldDouble *MEDCouplingFieldDouble::deepCpy() const
184 {
185   return cloneWithMesh(true);
186 }
187
188 /*!
189  * Creates a new MEDCouplingFieldDouble of given
190  * \ref MEDCouplingTemporalDisc "temporal discretization". The result field either
191  * shares the data array(s) with \a this field, or holds a deep copy of it, depending on
192  * \a deepCopy parameter. But the underlying \b mesh is always **shallow copied**.
193  * \param [in] td - the type of time discretization of the created field, one of
194  *        (\ref ParaMEDMEM::NO_TIME "NO_TIME", 
195  *         \ref ParaMEDMEM::ONE_TIME "ONE_TIME", 
196  *         \ref ParaMEDMEM::LINEAR_TIME "LINEAR_TIME", 
197  *         \ref ParaMEDMEM::CONST_ON_TIME_INTERVAL "CONST_ON_TIME_INTERVAL").
198  * \param [in] deepCopy - if \c true, the copy of the underlying data arrays is
199  *         deep, else all data arrays of \a this field are shared by the new field.
200  * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
201  *         caller is to delete this field using decrRef() as it is no more needed. 
202  * 
203  * \ref cpp_mcfielddouble_buildNewTimeReprFromThis "Here is a C++ example."<br>
204  * \ref py_mcfielddouble_buildNewTimeReprFromThis "Here is a Python example."
205  * \sa clone()
206  */
207 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildNewTimeReprFromThis(TypeOfTimeDiscretization td, bool deepCopy) const
208 {
209   MEDCouplingTimeDiscretization *tdo=_time_discr->buildNewTimeReprFromThis(td,deepCopy);
210   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretization> disc;
211   if(_type)
212     disc=_type->clone();
213   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),tdo,disc.retn());
214   ret->setMesh(getMesh());
215   ret->setName(getName());
216   ret->setDescription(getDescription());
217   return ret.retn();
218 }
219
220 /*!
221  * This method converts a field on nodes (\a this) to a cell field (returned field). The convertion is a \b non \b conservative remapping !
222  * This method is useful only for users that need a fast convertion from node to cell spatial discretization. The algorithm applied is simply to attach
223  * to each cell the average of values on nodes constituting this cell.
224  *
225  * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
226  *         caller is to delete this field using decrRef() as it is no more needed. The returned field will share the same mesh object object than those in \a this.
227  * \throw If \a this spatial discretization is empty or not ON_NODES.
228  * \throw If \a this is not coherent (see MEDCouplingFieldDouble::checkCoherency).
229  * 
230  * \warning This method is a \b non \b conservative method of remapping from node spatial discretization to cell spatial discretization.
231  * If a conservative method of interpolation is required ParaMEDMEM::MEDCouplingRemapper class should be used instead with "P1P0" method.
232  */
233 MEDCouplingFieldDouble *MEDCouplingFieldDouble::nodeToCellDiscretization() const
234 {
235   checkCoherency();
236   TypeOfField tf(getTypeOfField());
237   if(tf!=ON_NODES)
238     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::nodeToCellDiscretization : this field is expected to be on ON_NODES !");
239   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(clone(false));
240   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretizationP0> nsp(new MEDCouplingFieldDiscretizationP0);
241   ret->setDiscretization(nsp);
242   const MEDCouplingMesh *m(getMesh());//m is non empty thanks to checkCoherency call
243   int nbCells(m->getNumberOfCells());
244   std::vector<DataArrayDouble *> arrs(getArrays());
245   std::size_t sz(arrs.size());
246   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > outArrsSafe(sz); std::vector<DataArrayDouble *> outArrs(sz);
247   for(std::size_t j=0;j<sz;j++)
248     {
249       int nbCompo(arrs[j]->getNumberOfComponents());
250       outArrsSafe[j]=DataArrayDouble::New(); outArrsSafe[j]->alloc(nbCells,nbCompo);
251       outArrsSafe[j]->copyStringInfoFrom(*arrs[j]);
252       outArrs[j]=outArrsSafe[j];
253       double *pt(outArrsSafe[j]->getPointer());
254       const double *srcPt(arrs[j]->begin());
255       for(int i=0;i<nbCells;i++,pt+=nbCompo)
256         {
257           std::vector<int> nodeIds;
258           m->getNodeIdsOfCell(i,nodeIds);
259           std::fill(pt,pt+nbCompo,0.);
260           std::size_t nbNodesInCell(nodeIds.size());
261           for(std::size_t k=0;k<nbNodesInCell;k++)
262             std::transform(srcPt+nodeIds[k]*nbCompo,srcPt+(nodeIds[k]+1)*nbCompo,pt,pt,std::plus<double>());
263           if(nbNodesInCell!=0)
264             std::transform(pt,pt+nbCompo,pt,std::bind2nd(std::multiplies<double>(),1./((double)nbNodesInCell)));
265           else
266             {
267               std::ostringstream oss; oss << "MEDCouplingFieldDouble::nodeToCellDiscretization : Cell id #" << i << " has been detected to have no nodes !";
268               throw INTERP_KERNEL::Exception(oss.str().c_str());
269             }
270         }
271     }
272   ret->setArrays(outArrs);
273   return ret.retn();
274 }
275
276 /*!
277  * This method converts a field on cell (\a this) to a node field (returned field). The convertion is a \b non \b conservative remapping !
278  * This method is useful only for users that need a fast convertion from cell to node spatial discretization. The algorithm applied is simply to attach
279  * to each node the average of values on cell sharing this node. If \a this lies on a mesh having orphan nodes the values applied on them will be NaN (division by 0.).
280  *
281  * \return MEDCouplingFieldDouble* - a new instance of MEDCouplingFieldDouble. The
282  *         caller is to delete this field using decrRef() as it is no more needed. The returned field will share the same mesh object object than those in \a this.
283  * \throw If \a this spatial discretization is empty or not ON_CELLS.
284  * \throw If \a this is not coherent (see MEDCouplingFieldDouble::checkCoherency).
285  *
286  * \warning This method is a \b non \b conservative method of remapping from cell spatial discretization to node spatial discretization.
287  * If a conservative method of interpolation is required ParaMEDMEM::MEDCouplingRemapper class should be used instead with "P0P1" method.
288  */
289 MEDCouplingFieldDouble *MEDCouplingFieldDouble::cellToNodeDiscretization() const
290 {
291   checkCoherency();
292   TypeOfField tf(getTypeOfField());
293   if(tf!=ON_CELLS)
294     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::cellToNodeDiscretization : this field is expected to be on ON_CELLS !");
295   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(clone(false));
296   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretizationP1> nsp(new MEDCouplingFieldDiscretizationP1);
297   ret->setDiscretization(nsp);
298   const MEDCouplingMesh *m(getMesh());//m is non empty thanks to checkCoherency call
299   int nbCells(m->getNumberOfCells()),nbNodes(m->getNumberOfNodes());
300   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> rn(DataArrayInt::New()),rni(DataArrayInt::New());
301   m->getReverseNodalConnectivity(rn,rni);
302   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> rni2(rni->deltaShiftIndex());
303   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> rni3(rni2->convertToDblArr()); rni2=0;
304   std::vector<DataArrayDouble *> arrs(getArrays());
305   std::size_t sz(arrs.size());
306   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > outArrsSafe(sz); std::vector<DataArrayDouble *> outArrs(sz);
307   for(std::size_t j=0;j<sz;j++)
308     {
309       int nbCompo(arrs[j]->getNumberOfComponents());
310       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp(arrs[j]->selectByTupleIdSafe(rn->begin(),rn->end()));
311       outArrsSafe[j]=(tmp->accumulatePerChunck(rni->begin(),rni->end())); tmp=0;
312       outArrsSafe[j]->divideEqual(rni3);
313       outArrsSafe[j]->copyStringInfoFrom(*arrs[j]);
314       outArrs[j]=outArrsSafe[j];
315     }
316   ret->setArrays(outArrs);
317   return ret.retn();
318 }
319
320 /*!
321  * Copies tiny info (component names, name and description) from an \a other field to
322  * \a this one.
323  * \warning The underlying mesh is not renamed (for safety reason).
324  *  \param [in] other - the field to copy the tiny info from.
325  *  \throw If \a this->getNumberOfComponents() != \a other->getNumberOfComponents()
326  */
327 void MEDCouplingFieldDouble::copyTinyStringsFrom(const MEDCouplingField *other)
328 {
329   MEDCouplingField::copyTinyStringsFrom(other);
330   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
331   if(otherC)
332     {
333       _time_discr->copyTinyStringsFrom(*otherC->_time_discr);
334     }
335 }
336
337 /*!
338  * Copies only times, order and iteration from an \a other field to
339  * \a this one. The underlying mesh is not impacted by this method.
340  * Arrays are not impacted neither.
341  *  \param [in] other - the field to tiny attributes from.
342  *  \throw If \a this->getNumberOfComponents() != \a other->getNumberOfComponents()
343  */
344 void MEDCouplingFieldDouble::copyTinyAttrFrom(const MEDCouplingFieldDouble *other)
345 {
346   if(other)
347     {
348       _time_discr->copyTinyAttrFrom(*other->_time_discr);
349     }
350   
351 }
352
353 void MEDCouplingFieldDouble::copyAllTinyAttrFrom(const MEDCouplingFieldDouble *other)
354 {
355   copyTinyStringsFrom(other);
356   copyTinyAttrFrom(other);
357 }
358
359 /*!
360  * Returns a string describing \a this field. This string is outputted by \c print
361  * Python command. The string includes info on
362  * - name,
363  * - description,
364  * - \ref MEDCouplingSpatialDisc "spatial discretization",
365  * - \ref MEDCouplingTemporalDisc "time discretization",
366  * - \ref NatureOfField,
367  * - components,
368  * - mesh.
369  *
370  *  \return std::string - the string describing \a this field.
371  */
372 std::string MEDCouplingFieldDouble::simpleRepr() const
373 {
374   std::ostringstream ret;
375   ret << "FieldDouble with name : \"" << getName() << "\"\n";
376   ret << "Description of field is : \"" << getDescription() << "\"\n";
377   if(_type)
378     { ret << "FieldDouble space discretization is : " << _type->getStringRepr() << "\n"; }
379   else
380     { ret << "FieldDouble has no spatial discretization !\n"; }
381   if(_time_discr)
382     { ret << "FieldDouble time discretization is : " << _time_discr->getStringRepr() << "\n"; }
383   else
384     { ret << "FieldDouble has no time discretization !\n"; }
385   ret << "FieldDouble nature of field is : \"" << MEDCouplingNatureOfField::GetReprNoThrow(_nature) << "\"\n";
386   if(getArray())
387     {
388       if(getArray()->isAllocated())
389         {
390           int nbOfCompo=getArray()->getNumberOfComponents();
391           ret << "FieldDouble default array has " << nbOfCompo << " components and " << getArray()->getNumberOfTuples() << " tuples.\n";
392           ret << "FieldDouble default array has following info on components : ";
393           for(int i=0;i<nbOfCompo;i++)
394             ret << "\"" << getArray()->getInfoOnComponent(i) << "\" ";
395           ret << "\n";
396         }
397       else
398         {
399           ret << "Array set but not allocated !\n";
400         }
401     }
402   if(_mesh)
403     ret << "Mesh support information :\n__________________________\n" << _mesh->simpleRepr();
404   else
405     ret << "Mesh support information : No mesh set !\n";
406   return ret.str();
407 }
408
409 /*!
410  * Returns a string describing \a this field. The string includes info on
411  * - name,
412  * - description,
413  * - \ref MEDCouplingSpatialDisc "spatial discretization",
414  * - \ref MEDCouplingTemporalDisc "time discretization",
415  * - components,
416  * - mesh,
417  * - contents of data arrays.
418  *
419  *  \return std::string - the string describing \a this field.
420  */
421 std::string MEDCouplingFieldDouble::advancedRepr() const
422 {
423   std::ostringstream ret;
424   ret << "FieldDouble with name : \"" << getName() << "\"\n";
425   ret << "Description of field is : \"" << getDescription() << "\"\n";
426   if(_type)
427     { ret << "FieldDouble space discretization is : " << _type->getStringRepr() << "\n"; }
428   else
429     { ret << "FieldDouble has no space discretization set !\n"; }
430   if(_time_discr)
431     { ret << "FieldDouble time discretization is : " << _time_discr->getStringRepr() << "\n"; }
432   else
433     { ret << "FieldDouble has no time discretization set !\n"; }
434   if(getArray())
435     ret << "FieldDouble default array has " << getArray()->getNumberOfComponents() << " components and " << getArray()->getNumberOfTuples() << " tuples.\n";
436   if(_mesh)
437     ret << "Mesh support information :\n__________________________\n" << _mesh->advancedRepr();
438   else
439     ret << "Mesh support information : No mesh set !\n";
440   std::vector<DataArrayDouble *> arrays;
441   _time_discr->getArrays(arrays);
442   int arrayId=0;
443   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++,arrayId++)
444     {
445       ret << "Array #" << arrayId << " :\n__________\n";
446       if(*iter)
447         (*iter)->reprWithoutNameStream(ret);
448       else
449         ret << "Array empty !";
450       ret << "\n";
451     }
452   return ret.str();
453 }
454
455 void MEDCouplingFieldDouble::writeVTK(const std::string& fileName, bool isBinary) const
456 {
457   std::vector<const MEDCouplingFieldDouble *> fs(1,this);
458   MEDCouplingFieldDouble::WriteVTK(fileName,fs,isBinary);
459 }
460
461 bool MEDCouplingFieldDouble::isEqualIfNotWhy(const MEDCouplingField *other, double meshPrec, double valsPrec, std::string& reason) const
462 {
463   if(!other)
464     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::isEqualIfNotWhy : other instance is NULL !");
465   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
466   if(!otherC)
467     {
468       reason="field given in input is not castable in MEDCouplingFieldDouble !";
469       return false;
470     }
471   if(!MEDCouplingField::isEqualIfNotWhy(other,meshPrec,valsPrec,reason))
472     return false;
473   if(!_time_discr->isEqualIfNotWhy(otherC->_time_discr,valsPrec,reason))
474     {
475       reason.insert(0,"In FieldDouble time discretizations differ :");
476       return false;
477     }
478   return true;
479 }
480
481 /*!
482  * Checks equality of \a this and \a other field. Only numeric data is considered,
483  * i.e. names, description etc are not compared.
484  *  \param [in] other - the field to compare with.
485  *  \param [in] meshPrec - a precision used to compare node coordinates of meshes.
486  *  \param [in] valsPrec - a precision used to compare data arrays of the two fields.
487  *  \return bool - \c true if the two fields are equal, \c false else.
488  *  \throw If \a other == NULL.
489  *  \throw If the spatial discretization of \a this field is NULL.
490  */
491 bool MEDCouplingFieldDouble::isEqualWithoutConsideringStr(const MEDCouplingField *other, double meshPrec, double valsPrec) const
492 {
493   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
494   if(!otherC)
495     return false;
496   if(!MEDCouplingField::isEqualWithoutConsideringStr(other,meshPrec,valsPrec))
497     return false;
498   if(!_time_discr->isEqualWithoutConsideringStr(otherC->_time_discr,valsPrec))
499     return false;
500   return true;
501 }
502
503 /*!
504  * This method states if \a this and 'other' are compatibles each other before performing any treatment.
505  * This method is good for methods like : mergeFields.
506  * This method is not very demanding compared to areStrictlyCompatible that is better for operation on fields.
507  */
508 bool MEDCouplingFieldDouble::areCompatibleForMerge(const MEDCouplingField *other) const
509 {
510   if(!MEDCouplingField::areCompatibleForMerge(other))
511     return false;
512   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
513   if(!otherC)
514     return false;
515   if(!_time_discr->areCompatible(otherC->_time_discr))
516     return false;
517   return true;
518 }
519
520 /*!
521  * This method is more strict than MEDCouplingField::areCompatibleForMerge method.
522  * This method is used for operation on fields to operate a first check before attempting operation.
523  */
524 bool MEDCouplingFieldDouble::areStrictlyCompatible(const MEDCouplingField *other) const
525 {
526   std::string tmp;
527   if(!MEDCouplingField::areStrictlyCompatible(other))
528     return false;
529   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
530   if(!otherC)
531     return false;
532   if(!_time_discr->areStrictlyCompatible(otherC->_time_discr,tmp))
533     return false;
534   return true;
535 }
536
537 /*!
538  * Method with same principle than MEDCouplingFieldDouble::areStrictlyCompatible method except that
539  * number of components between \a this and 'other' can be different here (for operator*).
540  */
541 bool MEDCouplingFieldDouble::areCompatibleForMul(const MEDCouplingField *other) const
542 {
543   if(!MEDCouplingField::areStrictlyCompatible(other))
544     return false;
545   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
546   if(!otherC)
547     return false;
548   if(!_time_discr->areStrictlyCompatibleForMul(otherC->_time_discr))
549     return false;
550   return true;
551 }
552
553 /*!
554  * Method with same principle than MEDCouplingFieldDouble::areStrictlyCompatible method except that
555  * number of components between \a this and 'other' can be different here (for operator/).
556  */
557 bool MEDCouplingFieldDouble::areCompatibleForDiv(const MEDCouplingField *other) const
558 {
559   if(!MEDCouplingField::areStrictlyCompatible(other))
560     return false;
561   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
562   if(!otherC)
563     return false;
564   if(!_time_discr->areStrictlyCompatibleForDiv(otherC->_time_discr))
565     return false;
566   return true;
567 }
568
569 /*!
570  * This method is invocated before any attempt of melding. This method is very close to areStrictlyCompatible,
571  * except that \a this and other can have different number of components.
572  */
573 bool MEDCouplingFieldDouble::areCompatibleForMeld(const MEDCouplingFieldDouble *other) const
574 {
575   if(!MEDCouplingField::areStrictlyCompatible(other))
576     return false;
577   if(!_time_discr->areCompatibleForMeld(other->_time_discr))
578     return false;
579   return true;
580 }
581
582 /*!
583  * Permutes values of \a this field according to a given permutation array for cells
584  * renumbering. The underlying mesh is deeply copied and its cells are also permuted. 
585  * The number of cells remains the same; for that the permutation array \a old2NewBg
586  * should not contain equal ids.
587  * ** Warning, this method modifies the mesh aggreagated by \a this (by performing a deep copy ) **.
588  *
589  *  \param [in] old2NewBg - the permutation array in "Old to New" mode. Its length is
590  *         to be equal to \a this->getMesh()->getNumberOfCells().
591  *  \param [in] check - if \c true, \a old2NewBg is transformed to a new permutation
592  *         array, so that its maximal cell id to correspond to (be less than) the number
593  *         of cells in mesh. This new array is then used for the renumbering. If \a 
594  *         check == \c false, \a old2NewBg is used as is, that is less secure as validity 
595  *         of ids in \a old2NewBg is not checked.
596  *  \throw If the mesh is not set.
597  *  \throw If the spatial discretization of \a this field is NULL.
598  *  \throw If \a check == \c true and \a old2NewBg contains equal ids.
599  *  \throw If mesh nature does not allow renumbering (e.g. structured mesh).
600  * 
601  *  \ref cpp_mcfielddouble_renumberCells "Here is a C++ example".<br>
602  *  \ref  py_mcfielddouble_renumberCells "Here is a Python example".
603  */
604 void MEDCouplingFieldDouble::renumberCells(const int *old2NewBg, bool check)
605 {
606   renumberCellsWithoutMesh(old2NewBg,check);
607   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=_mesh->deepCpy();
608   m->renumberCells(old2NewBg,check);
609   setMesh(m);
610   updateTime();
611 }
612
613 /*!
614  * Permutes values of \a this field according to a given permutation array for cells
615  * renumbering. The underlying mesh is \b not permuted. 
616  * The number of cells remains the same; for that the permutation array \a old2NewBg
617  * should not contain equal ids.
618  * This method performs a part of job of renumberCells(). The reasonable use of this
619  * method is only for multi-field instances lying on the same mesh to avoid a
620  * systematic duplication and renumbering of _mesh attribute. 
621  * \warning Use this method with a lot of care!
622  *  \param [in] old2NewBg - the permutation array in "Old to New" mode. Its length is
623  *         to be equal to \a this->getMesh()->getNumberOfCells().
624  *  \param [in] check - if \c true, \a old2NewBg is transformed to a new permutation
625  *         array, so that its maximal cell id to correspond to (be less than) the number
626  *         of cells in mesh. This new array is then used for the renumbering. If \a 
627  *         check == \c false, \a old2NewBg is used as is, that is less secure as validity 
628  *         of ids in \a old2NewBg is not checked.
629  *  \throw If the mesh is not set.
630  *  \throw If the spatial discretization of \a this field is NULL.
631  *  \throw If \a check == \c true and \a old2NewBg contains equal ids.
632  *  \throw If mesh nature does not allow renumbering (e.g. structured mesh).
633  */
634 void MEDCouplingFieldDouble::renumberCellsWithoutMesh(const int *old2NewBg, bool check)
635 {
636    if(!_mesh)
637      throw INTERP_KERNEL::Exception("Expecting a defined mesh to be able to operate a renumbering !");
638    if(!((const MEDCouplingFieldDiscretization *)_type))
639      throw INTERP_KERNEL::Exception("Expecting a spatial discretization to be able to operate a renumbering !");
640   //
641   _type->renumberCells(old2NewBg,check);
642   std::vector<DataArrayDouble *> arrays;
643   _time_discr->getArrays(arrays);
644   std::vector<DataArray *> arrays2(arrays.size()); std::copy(arrays.begin(),arrays.end(),arrays2.begin());
645   _type->renumberArraysForCell(_mesh,arrays2,old2NewBg,check);
646   //
647   updateTime();
648 }
649
650 /*!
651  * Permutes values of \a this field according to a given permutation array for node
652  * renumbering. The underlying mesh is deeply copied and its nodes are also permuted. 
653  * The number of nodes can change, contrary to renumberCells().
654  *  \param [in] old2NewBg - the permutation array in "Old to New" mode. Its length is
655  *         to be equal to \a this->getMesh()->getNumberOfNodes().
656  *  \param [in] eps - a precision used to compare field values at merged nodes. If
657  *         the values differ more than \a eps, an exception is thrown.
658  *  \throw If the mesh is not set.
659  *  \throw If the spatial discretization of \a this field is NULL.
660  *  \throw If \a check == \c true and \a old2NewBg contains equal ids.
661  *  \throw If mesh nature does not allow renumbering (e.g. structured mesh).
662  *  \throw If values at merged nodes deffer more than \a eps.
663  * 
664  *  \ref cpp_mcfielddouble_renumberNodes "Here is a C++ example".<br>
665  *  \ref  py_mcfielddouble_renumberNodes "Here is a Python example".
666  */
667 void MEDCouplingFieldDouble::renumberNodes(const int *old2NewBg, double eps)
668 {
669   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
670   if(!meshC)
671     throw INTERP_KERNEL::Exception("Invalid mesh to apply renumberNodes on it !");
672   int nbOfNodes=meshC->getNumberOfNodes();
673   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
674   int newNbOfNodes=*std::max_element(old2NewBg,old2NewBg+nbOfNodes)+1;
675   renumberNodesWithoutMesh(old2NewBg,newNbOfNodes,eps);
676   meshC2->renumberNodes(old2NewBg,newNbOfNodes);
677   setMesh(meshC2);
678 }
679
680 /*!
681  * Permutes values of \a this field according to a given permutation array for nodes
682  * renumbering. The underlying mesh is \b not permuted. 
683  * The number of nodes can change, contrary to renumberCells().
684  * A given epsilon specifies a threshold of error in case of two nodes are merged but
685  * the difference of values on these nodes are higher than \a eps.
686  * This method performs a part of job of renumberNodes(), excluding node renumbering
687  * in mesh. The reasonable use of this
688  * method is only for multi-field instances lying on the same mesh to avoid a
689  * systematic duplication and renumbering of _mesh attribute. 
690  * \warning Use this method with a lot of care!
691  * \warning In case of an exception thrown, the contents of the data array can be
692  *         partially modified until the exception occurs. 
693  *  \param [in] old2NewBg - the permutation array in "Old to New" mode. Its length is
694  *         to be equal to \a this->getMesh()->getNumberOfNodes().
695  *  \param [in] newNbOfNodes - a number of nodes in the mesh after renumbering.
696  *  \param [in] eps - a precision used to compare field values at merged nodes. If
697  *         the values differ more than \a eps, an exception is thrown.
698  *  \throw If the mesh is not set.
699  *  \throw If the spatial discretization of \a this field is NULL.
700  *  \throw If values at merged nodes deffer more than \a eps.
701  */
702 void MEDCouplingFieldDouble::renumberNodesWithoutMesh(const int *old2NewBg, int newNbOfNodes, double eps)
703 {
704   if(!((const MEDCouplingFieldDiscretization *)_type))
705     throw INTERP_KERNEL::Exception("Expecting a spatial discretization to be able to operate a renumbering !");
706   std::vector<DataArrayDouble *> arrays;
707   _time_discr->getArrays(arrays);
708   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
709     if(*iter)
710       _type->renumberValuesOnNodes(eps,old2NewBg,newNbOfNodes,*iter);
711 }
712
713 /*!
714  * Returns all tuple ids of \a this scalar field that fit the range [\a vmin,
715  * \a vmax]. This method calls DataArrayDouble::getIdsInRange().
716  *  \param [in] vmin - a lower boundary of the range. Tuples with values less than \a
717  *         vmin are not included in the result array.
718  *  \param [in] vmax - an upper boundary of the range. Tuples with values more than \a
719  *         vmax are not included in the result array.
720  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of selected
721  *          tuples. The caller is to delete this array using decrRef() as it is no
722  *          more needed.
723  *  \throw If the data array is not set.
724  *  \throw If \a this->getNumberOfComponents() != 1.
725  */
726 DataArrayInt *MEDCouplingFieldDouble::getIdsInRange(double vmin, double vmax) const
727 {
728   if(getArray()==0)
729     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getIdsInRange : no default array set !");
730   return getArray()->getIdsInRange(vmin,vmax);
731 }
732
733 /*!
734  * Builds a newly created field, that the caller will have the responsability to deal with (decrRef()).
735  * This method makes the assumption that the field is correctly defined when this method is called, no check of this will be done.
736  * This method returns a restriction of \a this so that only tuples with ids specified in \a part will be contained in the returned field.
737  * Parameter \a part specifies **cell ids whatever the spatial discretization of this** (
738  * \ref ParaMEDMEM::ON_CELLS "ON_CELLS", 
739  * \ref ParaMEDMEM::ON_NODES "ON_NODES",
740  * \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT", 
741  * \ref ParaMEDMEM::ON_GAUSS_NE "ON_GAUSS_NE",
742  * \ref ParaMEDMEM::ON_NODES_KR "ON_NODES_KR").
743  *
744  * For example, \a this is a field on cells lying on a mesh that have 10 cells, \a part contains following cell ids [3,7,6].
745  * Then the returned field will lie on mesh having 3 cells and the returned field will contain 3 tuples.<br>
746  * Tuple #0 of the result field will refer to the cell #0 of returned mesh. The cell #0 of returned mesh will be equal to the cell #3 of \a this->getMesh().<br>
747  * Tuple #1 of the result field will refer to the cell #1 of returned mesh. The cell #1 of returned mesh will be equal to the cell #7 of \a this->getMesh().<br>
748  * Tuple #2 of the result field will refer to the cell #2 of returned mesh. The cell #2 of returned mesh will be equal to the cell #6 of \a this->getMesh().
749  *
750  * Let, for example, \a this be a field on nodes lying on a mesh that have 10 cells and 11 nodes, and \a part contains following cellIds [3,7,6].
751  * Thus \a this currently contains 11 tuples. If the restriction of mesh to 3 cells leads to a mesh with 6 nodes, then the returned field
752  * will contain 6 tuples and \a this field will lie on this restricted mesh. 
753  *
754  *  \param [in] part - an array of cell ids to include to the result field.
755  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The caller is to delete this field using decrRef() as it is no more needed.
756  *
757  *  \ref cpp_mcfielddouble_subpart1 "Here is a C++ example".<br>
758  *  \ref  py_mcfielddouble_subpart1 "Here is a Python example".
759  *  \sa MEDCouplingFieldDouble::buildSubPartRange
760  */
761
762 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildSubPart(const DataArrayInt *part) const
763 {
764   if(part==0)
765     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::buildSubPart : not empty array must be passed to this method !");
766   return buildSubPart(part->begin(),part->end());
767 }
768
769 /*!
770  * Builds a newly created field, that the caller will have the responsability to deal with.
771  * \n This method makes the assumption that \a this field is correctly defined when this method is called (\a this->checkCoherency() returns without any exception thrown), **no check of this will be done**.
772  * \n This method returns a restriction of \a this so that only tuple ids specified in [ \a partBg , \a partEnd ) will be contained in the returned field.
773  * \n Parameter [\a partBg, \a partEnd ) specifies **cell ids whatever the spatial discretization** of \a this (
774  * \ref ParaMEDMEM::ON_CELLS "ON_CELLS", 
775  * \ref ParaMEDMEM::ON_NODES "ON_NODES",
776  * \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT", 
777  * \ref ParaMEDMEM::ON_GAUSS_NE "ON_GAUSS_NE",
778  * \ref ParaMEDMEM::ON_NODES_KR "ON_NODES_KR").
779  *
780  * For example, \a this is a field on cells lying on a mesh that have 10 cells, \a partBg contains the following cell ids [3,7,6].
781  * Then the returned field will lie on mesh having 3 cells and will contain 3 tuples.
782  *- Tuple #0 of the result field will refer to the cell #0 of returned mesh. The cell #0 of returned mesh will be equal to the cell #3 of \a this->getMesh().
783  *- Tuple #1 of the result field will refer to the cell #1 of returned mesh. The cell #1 of returned mesh will be equal to the cell #7 of \a this->getMesh().
784  *- Tuple #2 of the result field will refer to the cell #2 of returned mesh. The cell #2 of returned mesh will be equal to the cell #6 of \a this->getMesh().
785  *
786  * Let, for example, \a this be a field on nodes lying on a mesh that have 10 cells and 11 nodes, and \a partBg contains following cellIds [3,7,6].
787  * Thus \a this currently contains 11 tuples. If the restriction of mesh to 3 cells leads to a mesh with 6 nodes, then the returned field
788  * will contain 6 tuples and \a this field will lie on this restricted mesh. 
789  *
790  * \param [in] partBg - start (included) of input range of cell ids to select [ \a partBg, \a partEnd )
791  * \param [in] partEnd - end (not included) of input range of cell ids to select [ \a partBg, \a partEnd )
792  * \return a newly allocated field the caller should deal with.
793  * 
794  * \throw if there is presence of an invalid cell id in [ \a partBg, \a partEnd ) regarding the number of cells of \a this->getMesh().
795  *
796  * \ref cpp_mcfielddouble_subpart1 "Here a C++ example."<br>
797  * \ref py_mcfielddouble_subpart1 "Here a Python example."
798  * \sa ParaMEDMEM::MEDCouplingFieldDouble::buildSubPart(const DataArrayInt *) const, MEDCouplingFieldDouble::buildSubPartRange
799  */
800 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildSubPart(const int *partBg, const int *partEnd) const
801 {
802   if(!((const MEDCouplingFieldDiscretization *)_type))
803     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::buildSubPart : Expecting a not NULL spatial discretization !");
804   DataArrayInt *arrSelect;
805   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=_type->buildSubMeshData(_mesh,partBg,partEnd,arrSelect);
806   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arrSelect2(arrSelect);
807   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(false);//quick shallow copy.
808   const MEDCouplingFieldDiscretization *disc=getDiscretization();
809   if(disc)
810     ret->setDiscretization(MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretization>(disc->clonePart(partBg,partEnd)));
811   ret->setMesh(m);
812   std::vector<DataArrayDouble *> arrays;
813   _time_discr->getArrays(arrays);
814   std::vector<DataArrayDouble *> arrs;
815   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > arrsSafe;
816   const int *arrSelBg=arrSelect->begin();
817   const int *arrSelEnd=arrSelect->end();
818   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
819     {
820       DataArrayDouble *arr=0;
821       if(*iter)
822         arr=(*iter)->selectByTupleIdSafe(arrSelBg,arrSelEnd);
823       arrs.push_back(arr); arrsSafe.push_back(arr);
824     }
825   ret->_time_discr->setArrays(arrs,0);
826   return ret.retn();
827 }
828
829 /*!
830  * This method is equivalent to MEDCouplingFieldDouble::buildSubPart, the only difference is that the input range of cell ids is
831  * given using a range given \a begin, \a end and \a step to optimize the part computation.
832  * 
833  * \sa MEDCouplingFieldDouble::buildSubPart
834  */
835 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildSubPartRange(int begin, int end, int step) const
836 {
837   if(!((const MEDCouplingFieldDiscretization *)_type))
838     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::buildSubPart : Expecting a not NULL spatial discretization !");
839   DataArrayInt *arrSelect;
840   int beginOut,endOut,stepOut;
841   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=_type->buildSubMeshDataRange(_mesh,begin,end,step,beginOut,endOut,stepOut,arrSelect);
842   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arrSelect2(arrSelect);
843   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(false);//quick shallow copy.
844   const MEDCouplingFieldDiscretization *disc=getDiscretization();
845   if(disc)
846     ret->setDiscretization(MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretization>(disc->clonePartRange(begin,end,step)));
847   ret->setMesh(m);
848   std::vector<DataArrayDouble *> arrays;
849   _time_discr->getArrays(arrays);
850   std::vector<DataArrayDouble *> arrs;
851   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > arrsSafe;
852   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
853     {
854       DataArrayDouble *arr=0;
855       if(*iter)
856         {
857           if(arrSelect)
858             {
859               const int *arrSelBg=arrSelect->begin();
860               const int *arrSelEnd=arrSelect->end();
861               arr=(*iter)->selectByTupleIdSafe(arrSelBg,arrSelEnd);
862             }
863           else
864             arr=(*iter)->selectByTupleId2(beginOut,endOut,stepOut);
865         }
866       arrs.push_back(arr); arrsSafe.push_back(arr);
867     }
868   ret->_time_discr->setArrays(arrs,0);
869   return ret.retn();
870 }
871
872 /*!
873  * Returns a type of \ref MEDCouplingTemporalDisc "time discretization" of \a this field.
874  *  \return ParaMEDMEM::TypeOfTimeDiscretization - an enum item describing the time
875  *          discretization type.
876  */
877 TypeOfTimeDiscretization MEDCouplingFieldDouble::getTimeDiscretization() const
878 {
879   return _time_discr->getEnum();
880 }
881
882 MEDCouplingFieldDouble::MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td):MEDCouplingField(type),
883                                                                                               _time_discr(MEDCouplingTimeDiscretization::New(td))
884 {
885 }
886
887 /*!
888  * ** WARINING : This method do not deeply copy neither mesh nor spatial discretization. Only a shallow copy (reference) is done for mesh and spatial discretization ! **
889  */
890 MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldTemplate& ft, TypeOfTimeDiscretization td):MEDCouplingField(ft,false),
891                                                                                                                 _time_discr(MEDCouplingTimeDiscretization::New(td))
892 {
893 }
894
895 MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCopy):MEDCouplingField(other,deepCopy),
896                                                                                                    _time_discr(other._time_discr->performCpy(deepCopy))
897 {
898 }
899
900 MEDCouplingFieldDouble::MEDCouplingFieldDouble(NatureOfField n, MEDCouplingTimeDiscretization *td, MEDCouplingFieldDiscretization *type):MEDCouplingField(type,n),_time_discr(td)
901 {
902 }
903
904 MEDCouplingFieldDouble::~MEDCouplingFieldDouble()
905 {
906   delete _time_discr;
907 }
908
909 /*!
910  * Checks if \a this field is correctly defined, else an exception is thrown.
911  *  \throw If the mesh is not set.
912  *  \throw If the data array is not set.
913  *  \throw If the spatial discretization of \a this field is NULL.
914  *  \throw If \a this->getTimeTolerance() < 0.
915  *  \throw If the temporal discretization data is incorrect.
916  *  \throw If mesh data does not correspond to field data.
917  */
918 void MEDCouplingFieldDouble::checkCoherency() const
919 {
920   if(!_mesh)
921     throw INTERP_KERNEL::Exception("Field invalid because no mesh specified !");
922   if(!((const MEDCouplingFieldDiscretization *)_type))
923     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::checkCoherency : no spatial discretization !");
924   _time_discr->checkCoherency();
925   _type->checkCoherencyBetween(_mesh,getArray());
926 }
927
928 /*!
929  * Accumulate values of a given component of \a this field.
930  *  \param [in] compId - the index of the component of interest.
931  *  \return double - a sum value of *compId*-th component.
932  *  \throw If the data array is not set.
933  *  \throw If \a the condition ( 0 <= \a compId < \a this->getNumberOfComponents() ) is
934  *         not respected.
935  */
936 double MEDCouplingFieldDouble::accumulate(int compId) const
937 {
938   if(getArray()==0)
939     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::accumulate : no default array defined !");
940   return getArray()->accumulate(compId);
941 }
942
943 /*!
944  * Accumulates values of each component of \a this array.
945  *  \param [out] res - an array of length \a this->getNumberOfComponents(), allocated 
946  *         by the caller, that is filled by this method with sum value for each
947  *         component.
948  *  \throw If the data array is not set.
949  */
950 void MEDCouplingFieldDouble::accumulate(double *res) const
951 {
952   if(getArray()==0)
953     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::accumulate : no default array defined !");
954   getArray()->accumulate(res);
955 }
956
957 /*!
958  * Returns the maximal value within \a this scalar field. Values of all arrays stored
959  * in \a this->_time_discr are checked.
960  *  \return double - the maximal value among all values of \a this field.
961  *  \throw If \a this->getNumberOfComponents() != 1
962  *  \throw If the data array is not set.
963  *  \throw If there is an empty data array in \a this field.
964  */
965 double MEDCouplingFieldDouble::getMaxValue() const
966 {
967   std::vector<DataArrayDouble *> arrays;
968   _time_discr->getArrays(arrays);
969   double ret=-std::numeric_limits<double>::max();
970   bool isExistingArr=false;
971   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
972     {
973       if(*iter)
974         {
975           isExistingArr=true;
976           int loc;
977           ret=std::max(ret,(*iter)->getMaxValue(loc));
978         }
979     }
980   if(!isExistingArr)
981     throw INTERP_KERNEL::Exception("getMaxValue : No arrays defined !");
982   return ret;
983 }
984
985 /*!
986  * Returns the maximal value and all its locations within \a this scalar field.
987  * Only the first of available data arrays is checked.
988  *  \param [out] tupleIds - a new instance of DataArrayInt containg indices of
989  *               tuples holding the maximal value. The caller is to delete it using
990  *               decrRef() as it is no more needed.
991  *  \return double - the maximal value among all values of the first array of \a this filed.
992  *  \throw If \a this->getNumberOfComponents() != 1.
993  *  \throw If there is an empty data array in \a this field.
994  */
995 double MEDCouplingFieldDouble::getMaxValue2(DataArrayInt*& tupleIds) const
996 {
997   std::vector<DataArrayDouble *> arrays;
998   _time_discr->getArrays(arrays);
999   double ret=-std::numeric_limits<double>::max();
1000   bool isExistingArr=false;
1001   tupleIds=0;
1002   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret1;
1003   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1004     {
1005       if(*iter)
1006         {
1007           isExistingArr=true;
1008           DataArrayInt *tmp;
1009           ret=std::max(ret,(*iter)->getMaxValue2(tmp));
1010           MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmpSafe(tmp);
1011           if(!((const DataArrayInt *)ret1))
1012             ret1=tmpSafe;
1013         }
1014     }
1015   if(!isExistingArr)
1016     throw INTERP_KERNEL::Exception("getMaxValue2 : No arrays defined !");
1017   tupleIds=ret1.retn();
1018   return ret;
1019 }
1020
1021 /*!
1022  * Returns the minimal value within \a this scalar field. Values of all arrays stored
1023  * in \a this->_time_discr are checked.
1024  *  \return double - the minimal value among all values of \a this field.
1025  *  \throw If \a this->getNumberOfComponents() != 1
1026  *  \throw If the data array is not set.
1027  *  \throw If there is an empty data array in \a this field.
1028  */
1029 double MEDCouplingFieldDouble::getMinValue() const
1030 {
1031   std::vector<DataArrayDouble *> arrays;
1032   _time_discr->getArrays(arrays);
1033   double ret=std::numeric_limits<double>::max();
1034   bool isExistingArr=false;
1035   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1036     {
1037       if(*iter)
1038         {
1039           isExistingArr=true;
1040           int loc;
1041           ret=std::min(ret,(*iter)->getMinValue(loc));
1042         }
1043     }
1044   if(!isExistingArr)
1045     throw INTERP_KERNEL::Exception("getMinValue : No arrays defined !");
1046   return ret;
1047 }
1048
1049 /*!
1050  * Returns the minimal value and all its locations within \a this scalar field.
1051  * Only the first of available data arrays is checked.
1052  *  \param [out] tupleIds - a new instance of DataArrayInt containg indices of
1053  *               tuples holding the minimal value. The caller is to delete it using
1054  *               decrRef() as it is no more needed.
1055  *  \return double - the minimal value among all values of the first array of \a this filed.
1056  *  \throw If \a this->getNumberOfComponents() != 1.
1057  *  \throw If there is an empty data array in \a this field.
1058  */
1059 double MEDCouplingFieldDouble::getMinValue2(DataArrayInt*& tupleIds) const
1060 {
1061   std::vector<DataArrayDouble *> arrays;
1062   _time_discr->getArrays(arrays);
1063   double ret=-std::numeric_limits<double>::max();
1064   bool isExistingArr=false;
1065   tupleIds=0;
1066   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret1;
1067   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1068     {
1069       if(*iter)
1070         {
1071           isExistingArr=true;
1072           DataArrayInt *tmp;
1073           ret=std::max(ret,(*iter)->getMinValue2(tmp));
1074           MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmpSafe(tmp);
1075           if(!((const DataArrayInt *)ret1))
1076             ret1=tmpSafe;
1077         }
1078     }
1079   if(!isExistingArr)
1080     throw INTERP_KERNEL::Exception("getMinValue2 : No arrays defined !");
1081   tupleIds=ret1.retn();
1082   return ret;
1083 }
1084
1085 /*!
1086  * Returns the average value of \a this scalar field.
1087  *  \return double - the average value over all values of the data array.
1088  *  \throw If \a this->getNumberOfComponents() != 1
1089  *  \throw If the data array is not set or it is empty.
1090  */
1091 double MEDCouplingFieldDouble::getAverageValue() const
1092 {
1093   if(getArray()==0)
1094     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getAverageValue : no default array defined !");
1095   return getArray()->getAverageValue();
1096 }
1097
1098 /*!
1099  * This method returns the euclidean norm of \a this field.
1100  * \f[
1101  * \sqrt{\sum_{0 \leq i < nbOfEntity}val[i]*val[i]}
1102  * \f]
1103  *  \throw If the data array is not set.
1104  */
1105 double MEDCouplingFieldDouble::norm2() const
1106 {
1107   if(getArray()==0)
1108     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::norm2 : no default array defined !");
1109   return getArray()->norm2();
1110 }
1111
1112 /*!
1113  * This method returns the max norm of \a this field.
1114  * \f[
1115  * \max_{0 \leq i < nbOfEntity}{abs(val[i])}
1116  * \f]
1117  *  \throw If the data array is not set.
1118  */
1119 double MEDCouplingFieldDouble::normMax() const
1120 {
1121   if(getArray()==0)
1122     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::normMax : no default array defined !");
1123   return getArray()->normMax();
1124 }
1125
1126 /*!
1127  * Computes sums of values of each component of \a this field wighted with
1128  * values returned by buildMeasureField().  
1129  *  \param [out] res - pointer to an array of result sum values, of size at least \a
1130  *         this->getNumberOfComponents(), that is to be allocated by the caller.
1131  *  \param [in] isWAbs - if \c true (default), \c abs() is applied to the weighs computed by
1132  *         buildMeasureField() that makes this method slower. If a user is sure that all
1133  *         cells of the underlying mesh have correct orientation, he can put \a isWAbs ==
1134  *         \c false that speeds up this method.
1135  *  \throw If the mesh is not set.
1136  *  \throw If the data array is not set.
1137  */
1138 void MEDCouplingFieldDouble::getWeightedAverageValue(double *res, bool isWAbs) const
1139 {
1140   if(getArray()==0)
1141     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getWeightedAverageValue : no default array defined !");
1142   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> w=buildMeasureField(isWAbs);
1143   double deno=w->getArray()->accumulate(0);
1144   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=getArray()->deepCpy();
1145   arr->multiplyEqual(w->getArray());
1146   std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::multiplies<double>(),1./deno));
1147   arr->accumulate(res);
1148 }
1149
1150 /*!
1151  * Computes a sum of values of a given component of \a this field wighted with
1152  * values returned by buildMeasureField().
1153  *  \param [in] compId - an index of the component of interest.
1154  *  \param [in] isWAbs - if \c true (default), \c abs() is applied to the weighs computed by
1155  *         buildMeasureField() that makes this method slower. If a user is sure that all
1156  *         cells of the underlying mesh have correct orientation, he can put \a isWAbs ==
1157  *         \c false that speeds up this method.
1158  *  \throw If the mesh is not set.
1159  *  \throw If the data array is not set.
1160  *  \throw If \a compId is not valid.
1161            A valid range is ( 0 <= \a compId < \a this->getNumberOfComponents() ).
1162  */
1163 double MEDCouplingFieldDouble::getWeightedAverageValue(int compId, bool isWAbs) const
1164 {
1165   int nbComps=getArray()->getNumberOfComponents();
1166   if(compId<0 || compId>=nbComps)
1167     {
1168       std::ostringstream oss; oss << "MEDCouplingFieldDouble::getWeightedAverageValue : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
1169       throw INTERP_KERNEL::Exception(oss.str().c_str());
1170     }
1171   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
1172   getWeightedAverageValue(res,isWAbs);
1173   return res[compId];
1174 }
1175
1176 /*!
1177  * Returns the \c normL1 of values of a given component of \a this field:
1178  * \f[
1179  * \frac{\sum_{0 \leq i < nbOfEntity}|val[i]*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}
1180  * \f]
1181  *  \param [in] compId - an index of the component of interest.
1182  *  \throw If the mesh is not set.
1183  *  \throw If the spatial discretization of \a this field is NULL.
1184  *  \throw If \a compId is not valid.
1185            A valid range is ( 0 <= \a compId < \a this->getNumberOfComponents() ).
1186  */
1187 double MEDCouplingFieldDouble::normL1(int compId) const
1188 {
1189   if(!_mesh)
1190     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1 !");
1191   if(!((const MEDCouplingFieldDiscretization *)_type))
1192     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL1 !");
1193   int nbComps=getArray()->getNumberOfComponents();
1194   if(compId<0 || compId>=nbComps)
1195     {
1196       std::ostringstream oss; oss << "MEDCouplingFieldDouble::normL1 : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
1197       throw INTERP_KERNEL::Exception(oss.str().c_str());
1198     }
1199   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
1200   _type->normL1(_mesh,getArray(),res);
1201   return res[compId];
1202 }
1203
1204 /*!
1205  * Returns the \c normL1 of values of each component of \a this field:
1206  * \f[
1207  * \frac{\sum_{0 \leq i < nbOfEntity}|val[i]*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}
1208  * \f]
1209  *  \param [out] res - pointer to an array of result values, of size at least \a
1210  *         this->getNumberOfComponents(), that is to be allocated by the caller.
1211  *  \throw If the mesh is not set.
1212  *  \throw If the spatial discretization of \a this field is NULL.
1213  */
1214 void MEDCouplingFieldDouble::normL1(double *res) const
1215 {
1216   if(!_mesh)
1217     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1");
1218   if(!((const MEDCouplingFieldDiscretization *)_type))
1219     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL1 !");
1220   _type->normL1(_mesh,getArray(),res);
1221 }
1222
1223 /*!
1224  * Returns the \c normL2 of values of a given component of \a this field:
1225  * \f[
1226  * \sqrt{\frac{\sum_{0 \leq i < nbOfEntity}|val[i]^{2}*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}}
1227  * \f]
1228  *  \param [in] compId - an index of the component of interest.
1229  *  \throw If the mesh is not set.
1230  *  \throw If the spatial discretization of \a this field is NULL.
1231  *  \throw If \a compId is not valid.
1232            A valid range is ( 0 <= \a compId < \a this->getNumberOfComponents() ).
1233  */
1234 double MEDCouplingFieldDouble::normL2(int compId) const
1235 {
1236   if(!_mesh)
1237     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL2");
1238   if(!((const MEDCouplingFieldDiscretization *)_type))
1239     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL2 !");
1240   int nbComps=getArray()->getNumberOfComponents();
1241   if(compId<0 || compId>=nbComps)
1242     {
1243       std::ostringstream oss; oss << "MEDCouplingFieldDouble::normL2 : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
1244       throw INTERP_KERNEL::Exception(oss.str().c_str());
1245     }
1246   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
1247   _type->normL2(_mesh,getArray(),res);
1248   return res[compId];
1249 }
1250
1251 /*!
1252  * Returns the \c normL2 of values of each component of \a this field:
1253  * \f[
1254  * \sqrt{\frac{\sum_{0 \leq i < nbOfEntity}|val[i]^{2}*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}}
1255  * \f]
1256  *  \param [out] res - pointer to an array of result values, of size at least \a
1257  *         this->getNumberOfComponents(), that is to be allocated by the caller.
1258  *  \throw If the mesh is not set.
1259  *  \throw If the spatial discretization of \a this field is NULL.
1260  */
1261 void MEDCouplingFieldDouble::normL2(double *res) const
1262 {
1263   if(!_mesh)
1264     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL2");
1265   if(!((const MEDCouplingFieldDiscretization *)_type))
1266     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL2 !");
1267   _type->normL2(_mesh,getArray(),res);
1268 }
1269
1270 /*!
1271  * Computes a sum of values of a given component of \a this field multiplied by
1272  * values returned by buildMeasureField().
1273  * This method is useful to check the conservativity of interpolation method.
1274  *  \param [in] compId - an index of the component of interest.
1275  *  \param [in] isWAbs - if \c true (default), \c abs() is applied to the weighs computed by
1276  *         buildMeasureField() that makes this method slower. If a user is sure that all
1277  *         cells of the underlying mesh have correct orientation, he can put \a isWAbs ==
1278  *         \c false that speeds up this method.
1279  *  \throw If the mesh is not set.
1280  *  \throw If the data array is not set.
1281  *  \throw If \a compId is not valid.
1282            A valid range is ( 0 <= \a compId < \a this->getNumberOfComponents() ).
1283  */
1284 double MEDCouplingFieldDouble::integral(int compId, bool isWAbs) const
1285 {
1286   if(!_mesh)
1287     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform integral");
1288   if(!((const MEDCouplingFieldDiscretization *)_type))
1289     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform integral !");
1290   int nbComps=getArray()->getNumberOfComponents();
1291   if(compId<0 || compId>=nbComps)
1292     {
1293       std::ostringstream oss; oss << "MEDCouplingFieldDouble::integral : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
1294       throw INTERP_KERNEL::Exception(oss.str().c_str());
1295     }
1296   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
1297   _type->integral(_mesh,getArray(),isWAbs,res);
1298   return res[compId];
1299 }
1300
1301 /*!
1302  * Computes a sum of values of each component of \a this field multiplied by
1303  * values returned by buildMeasureField().
1304  * This method is useful to check the conservativity of interpolation method.
1305  *  \param [in] isWAbs - if \c true (default), \c abs() is applied to the weighs computed by
1306  *         buildMeasureField() that makes this method slower. If a user is sure that all
1307  *         cells of the underlying mesh have correct orientation, he can put \a isWAbs ==
1308  *         \c false that speeds up this method.
1309  *  \param [out] res - pointer to an array of result sum values, of size at least \a
1310  *         this->getNumberOfComponents(), that is to be allocated by the caller.
1311  *  \throw If the mesh is not set.
1312  *  \throw If the data array is not set.
1313  *  \throw If the spatial discretization of \a this field is NULL.
1314  */
1315 void MEDCouplingFieldDouble::integral(bool isWAbs, double *res) const
1316 {
1317   if(!_mesh)
1318     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform integral2");
1319   if(!((const MEDCouplingFieldDiscretization *)_type))
1320     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform integral2 !");
1321   _type->integral(_mesh,getArray(),isWAbs,res);
1322 }
1323
1324 /*!
1325  * Returns a value at a given cell of a structured mesh. The cell is specified by its
1326  * (i,j,k) index.
1327  *  \param [in] i - a index of node coordinates array along X axis. The cell is
1328  *         located between the i-th and ( i + 1 )-th nodes along X axis.
1329  *  \param [in] j - a index of node coordinates array along Y axis. The cell is
1330  *         located between the j-th and ( j + 1 )-th nodes along Y axis.
1331  *  \param [in] k - a index of node coordinates array along Z axis. The cell is
1332  *         located between the k-th and ( k + 1 )-th nodes along Z axis.
1333  *  \param [out] res - pointer to an array returning a feild value, of size at least
1334  *         \a this->getNumberOfComponents(), that is to be allocated by the caller.
1335  *  \throw If the spatial discretization of \a this field is NULL.
1336  *  \throw If the mesh is not set.
1337  *  \throw If the mesh is not a structured one.
1338  *
1339  *  \ref cpp_mcfielddouble_getValueOnPos "Here is a C++ example".<br>
1340  *  \ref  py_mcfielddouble_getValueOnPos "Here is a Python example".
1341  */
1342 void MEDCouplingFieldDouble::getValueOnPos(int i, int j, int k, double *res) const
1343 {
1344   const DataArrayDouble *arr=_time_discr->getArray();
1345   if(!_mesh)
1346     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOnPos");
1347   if(!((const MEDCouplingFieldDiscretization *)_type))
1348     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOnPos !");
1349   _type->getValueOnPos(arr,_mesh,i,j,k,res);
1350 }
1351
1352 /*!
1353  * Returns a value of \a this at a given point using spatial discretization.
1354  *  \param [in] spaceLoc - the point of interest.
1355  *  \param [out] res - pointer to an array returning a feild value, of size at least
1356  *         \a this->getNumberOfComponents(), that is to be allocated by the caller.
1357  *  \throw If the spatial discretization of \a this field is NULL.
1358  *  \throw If the mesh is not set.
1359  *  \throw If \a spaceLoc is out of the spatial discretization.
1360  *
1361  *  \ref cpp_mcfielddouble_getValueOn "Here is a C++ example".<br>
1362  *  \ref  py_mcfielddouble_getValueOn "Here is a Python example".
1363  */
1364 void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const
1365 {
1366   const DataArrayDouble *arr=_time_discr->getArray();
1367   if(!_mesh)
1368     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOn");
1369   if(!((const MEDCouplingFieldDiscretization *)_type))
1370     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOnPos !");
1371   _type->getValueOn(arr,_mesh,spaceLoc,res);
1372 }
1373
1374 /*!
1375  * Returns values of \a this at given points using spatial discretization.
1376  *  \param [in] spaceLoc - coordinates of points of interest in full-interlace
1377  *          mode. This array is to be of size ( \a nbOfPoints * \a this->getNumberOfComponents() ).
1378  *  \param [in] nbOfPoints - number of points of interest.
1379  *  \return DataArrayDouble * - a new instance of DataArrayDouble holding field
1380  *         values relating to the input points. This array is of size \a nbOfPoints
1381  *         tuples per \a this->getNumberOfComponents() components. The caller is to 
1382  *         delete this array using decrRef() as it is no more needed.
1383  *  \throw If the spatial discretization of \a this field is NULL.
1384  *  \throw If the mesh is not set.
1385  *  \throw If any point in \a spaceLoc is out of the spatial discretization.
1386  *
1387  *  \ref cpp_mcfielddouble_getValueOnMulti "Here is a C++ example".<br>
1388  *  \ref  py_mcfielddouble_getValueOnMulti "Here is a Python example".
1389  */
1390 DataArrayDouble *MEDCouplingFieldDouble::getValueOnMulti(const double *spaceLoc, int nbOfPoints) const
1391 {
1392   const DataArrayDouble *arr=_time_discr->getArray();
1393   if(!_mesh)
1394     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOnMulti");
1395   if(!((const MEDCouplingFieldDiscretization *)_type))
1396     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOnMulti !");
1397   return _type->getValueOnMulti(arr,_mesh,spaceLoc,nbOfPoints);
1398 }
1399
1400 /*!
1401  * Returns a value of \a this field at a given point at a given time using spatial discretization.
1402  * If the time is not covered by \a this->_time_discr, an exception is thrown.
1403  *  \param [in] spaceLoc - the point of interest.
1404  *  \param [in] time - the time of interest.
1405  *  \param [out] res - pointer to an array returning a feild value, of size at least
1406  *         \a this->getNumberOfComponents(), that is to be allocated by the caller.
1407  *  \throw If the spatial discretization of \a this field is NULL.
1408  *  \throw If the mesh is not set.
1409  *  \throw If \a spaceLoc is out of the spatial discretization.
1410  *  \throw If \a time is not covered by \a this->_time_discr.
1411  *
1412  *  \ref cpp_mcfielddouble_getValueOn_time "Here is a C++ example".<br>
1413  *  \ref  py_mcfielddouble_getValueOn_time "Here is a Python example".
1414  */
1415 void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, double *res) const
1416 {
1417   std::vector< const DataArrayDouble *> arrs=_time_discr->getArraysForTime(time);
1418   if(!_mesh)
1419     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOn");
1420   if(!((const MEDCouplingFieldDiscretization *)_type))
1421     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOn !");
1422   std::vector<double> res2;
1423   for(std::vector< const DataArrayDouble *>::const_iterator iter=arrs.begin();iter!=arrs.end();iter++)
1424     {
1425       int sz=(int)res2.size();
1426       res2.resize(sz+(*iter)->getNumberOfComponents());
1427       _type->getValueOn(*iter,_mesh,spaceLoc,&res2[sz]);
1428     }
1429   _time_discr->getValueForTime(time,res2,res);
1430 }
1431
1432 /*!
1433  * Apply a liner function to a given component of \a this field, so that
1434  * a component value <em>(x)</em> becomes \f$ a * x + b \f$.
1435  *  \param [in] a - the first coefficient of the function.
1436  *  \param [in] b - the second coefficient of the function.
1437  *  \param [in] compoId - the index of component to modify.
1438  *  \throw If the data array(s) is(are) not set.
1439  */
1440 void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId)
1441 {
1442   _time_discr->applyLin(a,b,compoId);
1443 }
1444
1445 /*!
1446  * This method sets \a this to a uniform scalar field with one component.
1447  * All tuples will have the same value 'value'.
1448  * An exception is thrown if no underlying mesh is defined.
1449  */
1450 MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator=(double value) throw(INTERP_KERNEL::Exception)
1451 {
1452   if(!_mesh)
1453     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::operator= : no mesh defined !");
1454   if(!((const MEDCouplingFieldDiscretization *)_type))
1455     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform operator = !");
1456   int nbOfTuple=_type->getNumberOfTuples(_mesh);
1457   _time_discr->setOrCreateUniformValueOnAllComponents(nbOfTuple,value);
1458   return *this;
1459 }
1460
1461 /*!
1462  * Creates data array(s) of \a this field by using a C function for value generation.
1463  *  \param [in] nbOfComp - the number of components for \a this field to have.
1464  *  \param [in] func - the function used to compute values of \a this field.
1465  *         This function is to compute a field value basing on coordinates of value
1466  *         location point.
1467  *  \throw If the mesh is not set.
1468  *  \throw If \a func returns \c false.
1469  *  \throw If the spatial discretization of \a this field is NULL.
1470  *
1471  *  \ref cpp_mcfielddouble_fillFromAnalytic_c_func "Here is a C++ example".
1472  */
1473 void MEDCouplingFieldDouble::fillFromAnalytic(int nbOfComp, FunctionToEvaluate func)
1474 {
1475   if(!_mesh)
1476     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic : no mesh defined !");
1477   if(!((const MEDCouplingFieldDiscretization *)_type))
1478     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic !");
1479   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1480   _time_discr->fillFromAnalytic(loc,nbOfComp,func);
1481 }
1482
1483 /*!
1484  * Creates data array(s) of \a this field by using a function for value generation.<br>
1485  * The function is applied to coordinates of value location points. For example, if
1486  * \a this field is on cells, the function is applied to cell barycenters.
1487  * For more info on supported expressions that can be used in the function, see \ref
1488  * MEDCouplingArrayApplyFuncExpr. <br>
1489  * The function can include arbitrary named variables
1490  * (e.g. "x","y" or "va44") to refer to components of point coordinates. Names of
1491  * variables are sorted in \b alphabetical \b order to associate a variable name with a
1492  * component. For example, in the expression "2*x+z", "x" stands for the component #0
1493  * and "z" stands for the component #1 (\b not #2)!<br>
1494  * In a general case, a value resulting from the function evaluation is assigned to all
1495  * components of a field value. But there is a possibility to have its own expression for
1496  * each component within one function. For this purpose, there are predefined variable
1497  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
1498  * the component #0 etc). A factor of such a variable is added to the
1499  * corresponding component only.<br>
1500  * For example, \a nbOfComp == 4, coordinates of a 3D point are (1.,3.,7.), then
1501  *   - "2*x + z"               produces (5.,5.,5.,5.)
1502  *   - "2*x + 0*y + z"         produces (9.,9.,9.,9.)
1503  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,4.)
1504  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,4.)
1505  *
1506  *  \param [in] nbOfComp - the number of components for \a this field to have.
1507  *  \param [in] func - the function used to compute values of \a this field.
1508  *         This function is used to compute a field value basing on coordinates of value
1509  *         location point. For example, if \a this field is on cells, the function
1510  *         is applied to cell barycenters.
1511  *  \throw If the mesh is not set.
1512  *  \throw If the spatial discretization of \a this field is NULL.
1513  *  \throw If computing \a func fails.
1514  *
1515  *  \ref cpp_mcfielddouble_fillFromAnalytic "Here is a C++ example".<br>
1516  *  \ref  py_mcfielddouble_fillFromAnalytic "Here is a Python example".
1517  */
1518 void MEDCouplingFieldDouble::fillFromAnalytic(int nbOfComp, const std::string& func)
1519 {
1520   if(!_mesh)
1521     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic : no mesh defined !");
1522   if(!((const MEDCouplingFieldDiscretization *)_type))
1523     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic !");
1524   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1525   _time_discr->fillFromAnalytic(loc,nbOfComp,func);
1526 }
1527
1528 /*!
1529  * Creates data array(s) of \a this field by using a function for value generation.<br>
1530  * The function is applied to coordinates of value location points. For example, if
1531  * \a this field is on cells, the function is applied to cell barycenters.<br>
1532  * This method differs from
1533  * \ref ParaMEDMEM::MEDCouplingFieldDouble::fillFromAnalytic(int nbOfComp, const std::string& func) "fillFromAnalytic()"
1534  * by the way how variable
1535  * names, used in the function, are associated with components of coordinates of field
1536  * location points; here, a variable name corresponding to a component is retrieved from
1537  * a corresponding node coordinates array (where it is set via
1538  * DataArrayDouble::setInfoOnComponent()).<br>
1539  * For more info on supported expressions that can be used in the function, see \ref
1540  * MEDCouplingArrayApplyFuncExpr. <br> 
1541  * In a general case, a value resulting from the function evaluation is assigned to all
1542  * components of a field value. But there is a possibility to have its own expression for
1543  * each component within one function. For this purpose, there are predefined variable
1544  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
1545  * the component #0 etc). A factor of such a variable is added to the
1546  * corresponding component only.<br>
1547  * For example, \a nbOfComp == 4, names of spatial components are "x", "y" and "z",
1548  * coordinates of a 3D point are (1.,3.,7.), then
1549  *   - "2*x + z"               produces (9.,9.,9.,9.)
1550  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,8.)
1551  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,8.)
1552  *
1553  *  \param [in] nbOfComp - the number of components for \a this field to have.
1554  *  \param [in] func - the function used to compute values of \a this field.
1555  *         This function is used to compute a field value basing on coordinates of value
1556  *         location point. For example, if \a this field is on cells, the function
1557  *         is applied to cell barycenters.
1558  *  \throw If the mesh is not set.
1559  *  \throw If the spatial discretization of \a this field is NULL.
1560  *  \throw If computing \a func fails.
1561  *
1562  *  \ref cpp_mcfielddouble_fillFromAnalytic2 "Here is a C++ example".<br>
1563  *  \ref  py_mcfielddouble_fillFromAnalytic2 "Here is a Python example".
1564  */
1565 void MEDCouplingFieldDouble::fillFromAnalytic2(int nbOfComp, const std::string& func)
1566 {
1567   if(!_mesh)
1568     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic2 : no mesh defined !");
1569   if(!((const MEDCouplingFieldDiscretization *)_type))
1570     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic2 !");
1571   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1572   _time_discr->fillFromAnalytic2(loc,nbOfComp,func);
1573 }
1574
1575 /*!
1576  * Creates data array(s) of \a this field by using a function for value generation.<br>
1577  * The function is applied to coordinates of value location points. For example, if
1578  * \a this field is on cells, the function is applied to cell barycenters.<br>
1579  * This method differs from
1580  * \ref ParaMEDMEM::MEDCouplingFieldDouble::fillFromAnalytic(int nbOfComp, const std::string& func) "fillFromAnalytic()"
1581  * by the way how variable
1582  * names, used in the function, are associated with components of coordinates of field
1583  * location points; here, a component index of a variable is defined by a
1584  * rank of the variable within the input array \a varsOrder.<br>
1585  * For more info on supported expressions that can be used in the function, see \ref
1586  * MEDCouplingArrayApplyFuncExpr.
1587  * In a general case, a value resulting from the function evaluation is assigned to all
1588  * components of a field value. But there is a possibility to have its own expression for
1589  * each component within one function. For this purpose, there are predefined variable
1590  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
1591  * the component #0 etc). A factor of such a variable is added to the
1592  * corresponding component only.<br>
1593  * For example, \a nbOfComp == 4, names of
1594  * spatial components are given in \a varsOrder: ["x", "y","z"], coordinates of a
1595  * 3D point are (1.,3.,7.), then
1596  *   - "2*x + z"               produces (9.,9.,9.,9.)
1597  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,8.)
1598  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,8.)
1599  *
1600  *  \param [in] nbOfComp - the number of components for \a this field to have.
1601  *  \param [in] func - the function used to compute values of \a this field.
1602  *         This function is used to compute a field value basing on coordinates of value
1603  *         location point. For example, if \a this field is on cells, the function
1604  *         is applied to cell barycenters.
1605  *  \throw If the mesh is not set.
1606  *  \throw If the spatial discretization of \a this field is NULL.
1607  *  \throw If computing \a func fails.
1608  *
1609  *  \ref cpp_mcfielddouble_fillFromAnalytic3 "Here is a C++ example".<br>
1610  *  \ref  py_mcfielddouble_fillFromAnalytic3 "Here is a Python example".
1611  */
1612 void MEDCouplingFieldDouble::fillFromAnalytic3(int nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func)
1613 {
1614   if(!_mesh)
1615     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic2 : no mesh defined !");
1616   if(!((const MEDCouplingFieldDiscretization *)_type))
1617     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic3 !");
1618   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1619   _time_discr->fillFromAnalytic3(loc,nbOfComp,varsOrder,func);
1620 }
1621
1622 /*!
1623  * Modifies values of \a this field by applying a C function to each tuple of all
1624  * data arrays.
1625  *  \param [in] nbOfComp - the number of components for \a this field to have.
1626  *  \param [in] func - the function used to compute values of \a this field.
1627  *         This function is to compute a field value basing on a current field value.
1628  *  \throw If \a func returns \c false.
1629  *
1630  *  \ref cpp_mcfielddouble_applyFunc_c_func "Here is a C++ example".
1631  */
1632 void MEDCouplingFieldDouble::applyFunc(int nbOfComp, FunctionToEvaluate func)
1633 {
1634   _time_discr->applyFunc(nbOfComp,func);
1635 }
1636
1637 /*!
1638  * Fill \a this field with a given value.<br>
1639  * This method is a specialization of other overloaded methods. When \a nbOfComp == 1
1640  * this method is equivalent to ParaMEDMEM::MEDCouplingFieldDouble::operator=().
1641  *  \param [in] nbOfComp - the number of components for \a this field to have.
1642  *  \param [in] val - the value to assign to every atomic value of \a this field.
1643  *  \throw If the spatial discretization of \a this field is NULL.
1644  *  \throw If the mesh is not set.
1645  *
1646  *  \ref cpp_mcfielddouble_applyFunc_val "Here is a C++ example".<br>
1647  *  \ref  py_mcfielddouble_applyFunc_val "Here is a Python example".
1648  */
1649 void MEDCouplingFieldDouble::applyFunc(int nbOfComp, double val)
1650 {
1651   if(!_mesh)
1652     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::applyFunc : no mesh defined !");
1653   if(!((const MEDCouplingFieldDiscretization *)_type))
1654     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform applyFunc !");
1655   int nbOfTuple=_type->getNumberOfTuples(_mesh);
1656   _time_discr->setUniformValue(nbOfTuple,nbOfComp,val);
1657 }
1658
1659 /*!
1660  * Modifies values of \a this field by applying a function to each tuple of all
1661  * data arrays.
1662  * For more info on supported expressions that can be used in the function, see \ref
1663  * MEDCouplingArrayApplyFuncExpr. <br>
1664  * The function can include arbitrary named variables
1665  * (e.g. "x","y" or "va44") to refer to components of a field value. Names of
1666  * variables are sorted in \b alphabetical \b order to associate a variable name with a
1667  * component. For example, in the expression "2*x+z", "x" stands for the component #0
1668  * and "z" stands for the component #1 (\b not #2)!<br>
1669  * In a general case, a value resulting from the function evaluation is assigned to all
1670  * components of a field value. But there is a possibility to have its own expression for
1671  * each component within one function. For this purpose, there are predefined variable
1672  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
1673  * the component #0 etc). A factor of such a variable is added to the
1674  * corresponding component only.<br>
1675  * For example, \a nbOfComp == 4, components of a field value are (1.,3.,7.), then
1676  *   - "2*x + z"               produces (5.,5.,5.,5.)
1677  *   - "2*x + 0*y + z"         produces (9.,9.,9.,9.)
1678  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,4.)
1679  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,4.)
1680  *
1681  *  \param [in] nbOfComp - the number of components for \a this field to have.
1682  *  \param [in] func - the function used to compute values of \a this field.
1683  *         This function is to compute a field value basing on a current field value.
1684  *  \throw If computing \a func fails.
1685  *
1686  *  \ref cpp_mcfielddouble_applyFunc "Here is a C++ example".<br>
1687  *  \ref  py_mcfielddouble_applyFunc "Here is a Python example".
1688  */
1689 void MEDCouplingFieldDouble::applyFunc(int nbOfComp, const std::string& func)
1690 {
1691   _time_discr->applyFunc(nbOfComp,func);
1692 }
1693
1694
1695 /*!
1696  * Modifies values of \a this field by applying a function to each tuple of all
1697  * data arrays.
1698  * For more info on supported expressions that can be used in the function, see \ref
1699  * MEDCouplingArrayApplyFuncExpr. <br>
1700  * This method differs from
1701  * \ref ParaMEDMEM::MEDCouplingFieldDouble::applyFunc(int nbOfComp, const std::string& func) "applyFunc()"
1702  * by the way how variable
1703  * names, used in the function, are associated with components of field values;
1704  * here, a variable name corresponding to a component is retrieved from
1705  * component information of an array (where it is set via
1706  * DataArrayDouble::setInfoOnComponent()).<br>
1707  * In a general case, a value resulting from the function evaluation is assigned to all
1708  * components of a field value. But there is a possibility to have its own expression for
1709  * each component within one function. For this purpose, there are predefined variable
1710  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
1711  * the component #0 etc). A factor of such a variable is added to the
1712  * corresponding component only.<br>
1713  * For example, \a nbOfComp == 4, components of a field value are (1.,3.,7.), then
1714  *   - "2*x + z"               produces (5.,5.,5.,5.)
1715  *   - "2*x + 0*y + z"         produces (9.,9.,9.,9.)
1716  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,4.)
1717  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,4.)
1718  *
1719  *  \param [in] nbOfComp - the number of components for \a this field to have.
1720  *  \param [in] func - the function used to compute values of \a this field.
1721  *         This function is to compute a new field value basing on a current field value.
1722  *  \throw If computing \a func fails.
1723  *
1724  *  \ref cpp_mcfielddouble_applyFunc2 "Here is a C++ example".<br>
1725  *  \ref  py_mcfielddouble_applyFunc2 "Here is a Python example".
1726  */
1727 void MEDCouplingFieldDouble::applyFunc2(int nbOfComp, const std::string& func)
1728 {
1729   _time_discr->applyFunc2(nbOfComp,func);
1730 }
1731
1732 /*!
1733  * Modifies values of \a this field by applying a function to each tuple of all
1734  * data arrays.
1735  * This method differs from
1736  * \ref ParaMEDMEM::MEDCouplingFieldDouble::applyFunc(int nbOfComp, const std::string& func) "applyFunc()"
1737  * by the way how variable
1738  * names, used in the function, are associated with components of field values;
1739  * here, a component index of a variable is defined by a
1740  * rank of the variable within the input array \a varsOrder.<br>
1741  * For more info on supported expressions that can be used in the function, see \ref
1742  * MEDCouplingArrayApplyFuncExpr.
1743  * In a general case, a value resulting from the function evaluation is assigned to all
1744  * components of a field value. But there is a possibility to have its own expression for
1745  * each component within one function. For this purpose, there are predefined variable
1746  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
1747  * the component #0 etc). A factor of such a variable is added to the
1748  * corresponding component only.<br>
1749  * For example, \a nbOfComp == 4, names of
1750  * components are given in \a varsOrder: ["x", "y","z"], components of a
1751  * 3D vector are (1.,3.,7.), then
1752  *   - "2*x + z"               produces (9.,9.,9.,9.)
1753  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,8.)
1754  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,8.)
1755  *
1756  *  \param [in] nbOfComp - the number of components for \a this field to have.
1757  *  \param [in] func - the function used to compute values of \a this field.
1758  *         This function is to compute a new field value basing on a current field value.
1759  *  \throw If computing \a func fails.
1760  *
1761  *  \ref cpp_mcfielddouble_applyFunc3 "Here is a C++ example".<br>
1762  *  \ref  py_mcfielddouble_applyFunc3 "Here is a Python example".
1763  */
1764 void MEDCouplingFieldDouble::applyFunc3(int nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func)
1765 {
1766   _time_discr->applyFunc3(nbOfComp,varsOrder,func);
1767 }
1768
1769 /*!
1770  * Modifies values of \a this field by applying a function to each atomic value of all
1771  * data arrays. The function computes a new single value basing on an old single value.
1772  * For more info on supported expressions that can be used in the function, see \ref
1773  * MEDCouplingArrayApplyFuncExpr. <br>
1774  * The function can include **only one** arbitrary named variable
1775  * (e.g. "x","y" or "va44") to refer to a field atomic value. <br>
1776  * In a general case, a value resulting from the function evaluation is assigned to 
1777  * a single field value. But there is a possibility to have its own expression for
1778  * each component within one function. For this purpose, there are predefined variable
1779  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
1780  * the component #0 etc). A factor of such a variable is added to the
1781  * corresponding component only.<br>
1782  * For example, components of a field value are (1.,3.,7.), then
1783  *   - "2*x - 1"               produces (1.,5.,13.)
1784  *   - "2*x*IVec + (x+3)*KVec" produces (2.,0.,10.)
1785  *   - "2*x*IVec + (x+3)*KVec + 1" produces (3.,1.,11.)
1786  *
1787  *  \param [in] func - the function used to compute values of \a this field.
1788  *         This function is to compute a field value basing on a current field value.
1789  *  \throw If computing \a func fails.
1790  *
1791  *  \ref cpp_mcfielddouble_applyFunc_same_nb_comp "Here is a C++ example".<br>
1792  *  \ref  py_mcfielddouble_applyFunc_same_nb_comp "Here is a Python example".
1793  */
1794 void MEDCouplingFieldDouble::applyFunc(const std::string& func)
1795 {
1796   _time_discr->applyFunc(func);
1797 }
1798
1799 /*!
1800  * Applyies the function specified by the string repr 'func' on each tuples on all arrays contained in _time_discr.
1801  * The field will contain exactly the same number of components after the call.
1802  * Use is not warranted for the moment !
1803  */
1804 void MEDCouplingFieldDouble::applyFuncFast32(const std::string& func)
1805 {
1806   _time_discr->applyFuncFast32(func);
1807 }
1808
1809 /*!
1810  * Applyies the function specified by the string repr 'func' on each tuples on all arrays contained in _time_discr.
1811  * The field will contain exactly the same number of components after the call.
1812  * Use is not warranted for the moment !
1813  */
1814 void MEDCouplingFieldDouble::applyFuncFast64(const std::string& func)
1815 {
1816   _time_discr->applyFuncFast64(func);
1817 }
1818
1819 /*!
1820  * Returns number of components in the data array. For more info on the data arrays,
1821  * see \ref MEDCouplingArrayPage.
1822  *  \return int - the number of components in the data array.
1823  *  \throw If the data array is not set.
1824  */
1825 int MEDCouplingFieldDouble::getNumberOfComponents() const
1826 {
1827   if(getArray()==0)
1828     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getNumberOfComponents : No array specified !");
1829   return getArray()->getNumberOfComponents();
1830 }
1831
1832 /*!
1833  * Returns number of tuples in \a this field, that depends on 
1834  * - the number of entities in the underlying mesh
1835  * - \ref MEDCouplingSpatialDisc "spatial discretization" of \a this field (e.g. number
1836  * of Gauss points if \a this->getTypeOfField() == 
1837  * \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT").
1838  *
1839  * The returned value does **not depend** on the number of tuples in the data array
1840  * (which has to be equal to the returned value), \b contrary to
1841  * getNumberOfComponents() and getNumberOfValues() that retrieve information from the
1842  * data array.
1843  * \warning No checkCoherency() is done here.
1844  * For more info on the data arrays, see \ref MEDCouplingArrayPage.
1845  *  \return int - the number of tuples.
1846  *  \throw If the mesh is not set.
1847  *  \throw If the spatial discretization of \a this field is NULL.
1848  *  \throw If the spatial discretization is not fully defined.
1849  */
1850 int MEDCouplingFieldDouble::getNumberOfTuples() const
1851 {
1852   if(!_mesh)
1853     throw INTERP_KERNEL::Exception("Impossible to retrieve number of tuples because no mesh specified !");
1854   if(!((const MEDCouplingFieldDiscretization *)_type))
1855     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getNumberOfTuples !");
1856   return _type->getNumberOfTuples(_mesh);
1857 }
1858
1859 /*!
1860  * Returns number of atomic double values in the data array of \a this field.
1861  * For more info on the data arrays, see \ref MEDCouplingArrayPage.
1862  *  \return int - (number of tuples) * (number of components) of the
1863  *  data array.
1864  *  \throw If the data array is not set.
1865  */
1866 int MEDCouplingFieldDouble::getNumberOfValues() const
1867 {
1868   if(getArray()==0)
1869     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getNumberOfValues : No array specified !");
1870   return getArray()->getNbOfElems();
1871 }
1872
1873 /*!
1874  * Sets own modification time by the most recently modified element of data (the mesh,
1875  * the data array etc). For more info, see \ref MEDCouplingTimeLabelPage.
1876  */
1877 void MEDCouplingFieldDouble::updateTime() const
1878 {
1879   MEDCouplingField::updateTime();
1880   updateTimeWith(*_time_discr);
1881 }
1882
1883 std::size_t MEDCouplingFieldDouble::getHeapMemorySizeWithoutChildren() const
1884 {
1885   return MEDCouplingField::getHeapMemorySizeWithoutChildren();
1886 }
1887
1888 std::vector<const BigMemoryObject *> MEDCouplingFieldDouble::getDirectChildren() const
1889 {
1890   std::vector<const BigMemoryObject *> ret(MEDCouplingField::getDirectChildren());
1891   if(_time_discr)
1892     {
1893       std::vector<const BigMemoryObject *> ret2(_time_discr->getDirectChildren());
1894       ret.insert(ret.end(),ret2.begin(),ret2.end());
1895     }
1896   return ret;
1897 }
1898
1899 /*!
1900  * Sets \ref NatureOfField.
1901  *  \param [in] nat - an item of enum ParaMEDMEM::NatureOfField.
1902  */
1903 void MEDCouplingFieldDouble::setNature(NatureOfField nat)
1904 {
1905   MEDCouplingField::setNature(nat);
1906   if(_type)
1907     _type->checkCompatibilityWithNature(nat);
1908 }
1909
1910 /*!
1911  * This method synchronizes time information (time, iteration, order, time unit) regarding the information in \c this->_mesh.
1912  * \throw If no mesh is set in this. Or if \a this is not compatible with time setting (typically NO_TIME)
1913  */
1914 void MEDCouplingFieldDouble::synchronizeTimeWithMesh()
1915 {
1916   if(!_mesh)
1917     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::synchronizeTimeWithMesh : no mesh set in this !");
1918   int it=-1,ordr=-1;
1919   double val=_mesh->getTime(it,ordr);
1920   std::string timeUnit(_mesh->getTimeUnit());
1921   setTime(val,it,ordr);
1922   setTimeUnit(timeUnit);
1923 }
1924
1925 /*!
1926  * Returns a value of \a this field of type either
1927  * \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT" or
1928  * \ref ParaMEDMEM::ON_GAUSS_NE "ON_GAUSS_NE".
1929  *  \param [in] cellId - an id of cell of interest.
1930  *  \param [in] nodeIdInCell - a node index within the cell.
1931  *  \param [in] compoId - an index of component.
1932  *  \return double - the field value corresponding to the specified parameters.
1933  *  \throw If the data array is not set.
1934  *  \throw If the mesh is not set.
1935  *  \throw If the spatial discretization of \a this field is NULL.
1936  *  \throw If \a this field if of type other than 
1937  *         \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT" or
1938  *         \ref ParaMEDMEM::ON_GAUSS_NE "ON_GAUSS_NE".
1939  */
1940 double MEDCouplingFieldDouble::getIJK(int cellId, int nodeIdInCell, int compoId) const
1941 {
1942   if(!((const MEDCouplingFieldDiscretization *)_type))
1943     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getIJK !");
1944   return _type->getIJK(_mesh,getArray(),cellId,nodeIdInCell,compoId);
1945 }
1946
1947 /*!
1948  * Sets the data array. 
1949  *  \param [in] array - the data array holding values of \a this field. It's size
1950  *         should correspond to the mesh and
1951  *         \ref MEDCouplingSpatialDisc "spatial discretization" of \a this field
1952  *         (see getNumberOfTuples()), but this size is not checked here.
1953  */
1954 void MEDCouplingFieldDouble::setArray(DataArrayDouble *array)
1955 {
1956   _time_discr->setArray(array,this);
1957 }
1958
1959 /*!
1960  * Sets the data array holding values corresponding to an end of a time interval
1961  * for which \a this field is defined.
1962  *  \param [in] array - the data array holding values of \a this field. It's size
1963  *         should correspond to the mesh and
1964  *         \ref MEDCouplingSpatialDisc "spatial discretization" of \a this field
1965  *         (see getNumberOfTuples()), but this size is not checked here.
1966  */
1967 void MEDCouplingFieldDouble::setEndArray(DataArrayDouble *array)
1968 {
1969   _time_discr->setEndArray(array,this);
1970 }
1971
1972 /*!
1973  * Sets all data arrays needed to define the field values.
1974  *  \param [in] arrs - a vector of DataArrayDouble's holding values of \a this
1975  *         field. Size of each array should correspond to the mesh and
1976  *         \ref MEDCouplingSpatialDisc "spatial discretization" of \a this field
1977  *         (see getNumberOfTuples()), but this size is not checked here.
1978  *  \throw If number of arrays in \a arrs does not correspond to type of
1979  *         \ref MEDCouplingTemporalDisc "temporal discretization" of \a this field.
1980  */
1981 void MEDCouplingFieldDouble::setArrays(const std::vector<DataArrayDouble *>& arrs)
1982 {
1983   _time_discr->setArrays(arrs,this);
1984 }
1985
1986 void MEDCouplingFieldDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
1987 {
1988   tinyInfo.clear();
1989   _time_discr->getTinySerializationStrInformation(tinyInfo);
1990   tinyInfo.push_back(_name);
1991   tinyInfo.push_back(_desc);
1992   tinyInfo.push_back(getTimeUnit());
1993 }
1994
1995 /*!
1996  * This method retrieves some critical values to resize and prepare remote instance.
1997  * The first two elements returned in tinyInfo correspond to the parameters to give in constructor.
1998  * @param tinyInfo out parameter resized correctly after the call. The length of this vector is tiny.
1999  */
2000 void MEDCouplingFieldDouble::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
2001 {
2002   if(!((const MEDCouplingFieldDiscretization *)_type))
2003     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getTinySerializationIntInformation !");
2004   tinyInfo.clear();
2005   tinyInfo.push_back((int)_type->getEnum());
2006   tinyInfo.push_back((int)_time_discr->getEnum());
2007   tinyInfo.push_back((int)_nature);
2008   _time_discr->getTinySerializationIntInformation(tinyInfo);
2009   std::vector<int> tinyInfo2;
2010   _type->getTinySerializationIntInformation(tinyInfo2);
2011   tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
2012   tinyInfo.push_back((int)tinyInfo2.size());
2013 }
2014
2015 /*!
2016  * This method retrieves some critical values to resize and prepare remote instance.
2017  * @param tinyInfo out parameter resized correctly after the call. The length of this vector is tiny.
2018  */
2019 void MEDCouplingFieldDouble::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
2020 {
2021   if(!((const MEDCouplingFieldDiscretization *)_type))
2022     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getTinySerializationDbleInformation !");
2023   tinyInfo.clear();
2024   _time_discr->getTinySerializationDbleInformation(tinyInfo);
2025   std::vector<double> tinyInfo2;
2026   _type->getTinySerializationDbleInformation(tinyInfo2);
2027   tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
2028   tinyInfo.push_back((int)tinyInfo2.size());//very bad, lack of time to improve it
2029 }
2030
2031 /*!
2032  * This method has to be called to the new instance filled by CORBA, MPI, File...
2033  * @param tinyInfoI is the value retrieves from distant result of getTinySerializationIntInformation on source instance to be copied.
2034  * @param dataInt out parameter. If not null the pointer is already owned by \a this after the call of this method. In this case no decrRef must be applied.
2035  * @param arrays out parameter is a vector resized to the right size. The pointers in the vector is already owned by \a this after the call of this method.
2036  *               No decrRef must be applied to every instances in returned vector.
2037  */
2038 void MEDCouplingFieldDouble::resizeForUnserialization(const std::vector<int>& tinyInfoI, DataArrayInt *&dataInt, std::vector<DataArrayDouble *>& arrays)
2039 {
2040   if(!((const MEDCouplingFieldDiscretization *)_type))
2041     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform resizeForUnserialization !");
2042   dataInt=0;
2043   std::vector<int> tinyInfoITmp(tinyInfoI);
2044   int sz=tinyInfoITmp.back();
2045   tinyInfoITmp.pop_back();
2046   std::vector<int> tinyInfoITmp2(tinyInfoITmp.begin(),tinyInfoITmp.end()-sz);
2047   std::vector<int> tinyInfoI2(tinyInfoITmp2.begin()+3,tinyInfoITmp2.end());
2048   _time_discr->resizeForUnserialization(tinyInfoI2,arrays);
2049   std::vector<int> tinyInfoITmp3(tinyInfoITmp.end()-sz,tinyInfoITmp.end());
2050   _type->resizeForUnserialization(tinyInfoITmp3,dataInt);
2051 }
2052
2053 void MEDCouplingFieldDouble::finishUnserialization(const std::vector<int>& tinyInfoI, const std::vector<double>& tinyInfoD, const std::vector<std::string>& tinyInfoS)
2054 {
2055   if(!((const MEDCouplingFieldDiscretization *)_type))
2056     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform finishUnserialization !");
2057   std::vector<int> tinyInfoI2(tinyInfoI.begin()+3,tinyInfoI.end());
2058   //
2059   std::vector<double> tmp(tinyInfoD);
2060   int sz=(int)tinyInfoD.back();//very bad, lack of time to improve it
2061   tmp.pop_back();
2062   std::vector<double> tmp1(tmp.begin(),tmp.end()-sz);
2063   std::vector<double> tmp2(tmp.end()-sz,tmp.end());
2064   //
2065   _time_discr->finishUnserialization(tinyInfoI2,tmp1,tinyInfoS);
2066   _nature=(NatureOfField)tinyInfoI[2];
2067   _type->finishUnserialization(tmp2);
2068   int nbOfElemS=(int)tinyInfoS.size();
2069   _name=tinyInfoS[nbOfElemS-3];
2070   _desc=tinyInfoS[nbOfElemS-2];
2071   setTimeUnit(tinyInfoS[nbOfElemS-1]);
2072 }
2073
2074 /*!
2075  * Contrary to MEDCouplingPointSet class the returned arrays are \b not the responsabilities of the caller.
2076  * The values returned must be consulted only in readonly mode.
2077  */
2078 void MEDCouplingFieldDouble::serialize(DataArrayInt *&dataInt, std::vector<DataArrayDouble *>& arrays) const
2079 {
2080   if(!((const MEDCouplingFieldDiscretization *)_type))
2081     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform serialize !");
2082   _time_discr->getArrays(arrays);
2083   _type->getSerializationIntArray(dataInt);
2084 }
2085
2086 /*!
2087  * Tries to set an \a other mesh as the support of \a this field. An attempt fails, if
2088  * a current and the \a other meshes are different with use of specified equality
2089  * criteria, and then an exception is thrown.
2090  *  \param [in] other - the mesh to use as the field support if this mesh can be
2091  *         considered equal to the current mesh.
2092  *  \param [in] levOfCheck - defines equality criteria used for mesh comparison. For
2093  *         it's meaning explanation, see MEDCouplingMesh::checkGeoEquivalWith() which
2094  *         is used for mesh comparison.
2095  *  \param [in] precOnMesh - a precision used to compare nodes of the two meshes.
2096  *         It is used as \a prec parameter of MEDCouplingMesh::checkGeoEquivalWith().
2097  *  \param [in] eps - a precision used at node renumbering (if needed) to compare field
2098  *         values at merged nodes. If the values differ more than \a eps, an
2099  *         exception is thrown.
2100  *  \throw If the mesh is not set.
2101  *  \throw If \a other == NULL.
2102  *  \throw If any of the meshes is not well defined.
2103  *  \throw If the two meshes do not match.
2104  *  \throw If field values at merged nodes (if any) deffer more than \a eps.
2105  *
2106  *  \ref cpp_mcfielddouble_changeUnderlyingMesh "Here is a C++ example".<br>
2107  *  \ref  py_mcfielddouble_changeUnderlyingMesh "Here is a Python example".
2108  */
2109 void MEDCouplingFieldDouble::changeUnderlyingMesh(const MEDCouplingMesh *other, int levOfCheck, double precOnMesh, double eps)
2110 {
2111   if(_mesh==0 || other==0)
2112     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::changeUnderlyingMesh : is expected to operate on not null meshes !");
2113   DataArrayInt *cellCor=0,*nodeCor=0;
2114   other->checkGeoEquivalWith(_mesh,levOfCheck,precOnMesh,cellCor,nodeCor);
2115   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellCor2(cellCor),nodeCor2(nodeCor);
2116   if(cellCor)
2117     renumberCellsWithoutMesh(cellCor->getConstPointer(),false);
2118   if(nodeCor)
2119     renumberNodesWithoutMesh(nodeCor->getConstPointer(),nodeCor->getMaxValueInArray()+1,eps);
2120   setMesh(const_cast<MEDCouplingMesh *>(other));
2121 }
2122
2123 /*!
2124  * Subtracts another field from \a this one in case when the two fields have different
2125  * supporting meshes. The subtraction is performed provided that the tho meshes can be
2126  * considered equal with use of specified equality criteria, else an exception is thrown.
2127  * If the meshes match, the mesh of \a f is set to \a this field (\a this is permuted if 
2128  * necessary) and field values are subtracted. No interpolation is done here, only an
2129  * analysis of two underlying mesh is done to see if the meshes are geometrically
2130  * equivalent.<br>
2131  * The job of this method consists in calling
2132  * \a this->changeUnderlyingMesh() with \a f->getMesh() as the first parameter, and then
2133  * \a this -= \a f.<br>
2134  * This method requires that \a f and \a this are coherent (checkCoherency()) and that \a f
2135  * and \a this are coherent for a merge.<br>
2136  * "DM" in the method name stands for "different meshes".
2137  *  \param [in] f - the field to subtract from this.
2138  *  \param [in] levOfCheck - defines equality criteria used for mesh comparison. For
2139  *         it's meaning explanation, see MEDCouplingMesh::checkGeoEquivalWith() which
2140  *         is used for mesh comparison.
2141  *  \param [in] precOnMesh - a precision used to compare nodes of the two meshes.
2142  *         It is used as \a prec parameter of MEDCouplingMesh::checkGeoEquivalWith().
2143  *  \param [in] eps - a precision used at node renumbering (if needed) to compare field
2144  *         values at merged nodes. If the values differ more than \a eps, an
2145  *         exception is thrown.
2146  *  \throw If \a f == NULL.
2147  *  \throw If any of the meshes is not set or is not well defined.
2148  *  \throw If the two meshes do not match.
2149  *  \throw If the two fields are not coherent for merge.
2150  *  \throw If field values at merged nodes (if any) deffer more than \a eps.
2151  *
2152  *  \ref cpp_mcfielddouble_substractInPlaceDM "Here is a C++ example".<br>
2153  *  \ref  py_mcfielddouble_substractInPlaceDM "Here is a Python example".
2154  *  \sa changeUnderlyingMesh().
2155  */
2156 void MEDCouplingFieldDouble::substractInPlaceDM(const MEDCouplingFieldDouble *f, int levOfCheck, double precOnMesh, double eps)
2157 {
2158   checkCoherency();
2159   if(!f)
2160     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::substractInPlaceDM : input field is NULL !");
2161   f->checkCoherency();
2162   if(!areCompatibleForMerge(f))
2163     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::substractInPlaceDM : Fields are not compatible ; unable to apply mergeFields on them !");
2164   changeUnderlyingMesh(f->getMesh(),levOfCheck,precOnMesh,eps);
2165   operator-=(*f);
2166 }
2167
2168 /*!
2169  * Merges coincident nodes of the underlying mesh. If some nodes are coincident, the
2170  * underlying mesh is replaced by a new mesh instance where the coincident nodes are merged.
2171  *  \param [in] eps - a precision used to compare nodes of the two meshes.
2172  *  \param [in] epsOnVals - a precision used to compare field
2173  *         values at merged nodes. If the values differ more than \a epsOnVals, an
2174  *         exception is thrown.
2175  *  \return bool - \c true if some nodes have been merged and hence \a this field lies
2176  *         on another mesh.
2177  *  \throw If the mesh is of type not inheriting from MEDCouplingPointSet.
2178  *  \throw If the mesh is not well defined.
2179  *  \throw If the spatial discretization of \a this field is NULL.
2180  *  \throw If the data array is not set.
2181  *  \throw If field values at merged nodes (if any) deffer more than \a epsOnVals.
2182  */
2183 bool MEDCouplingFieldDouble::mergeNodes(double eps, double epsOnVals)
2184 {
2185   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
2186   if(!meshC)
2187     throw INTERP_KERNEL::Exception("Invalid support mesh to apply mergeNodes on it : must be a MEDCouplingPointSet one !");
2188   if(!((const MEDCouplingFieldDiscretization *)_type))
2189     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform mergeNodes !");
2190   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
2191   bool ret;
2192   int ret2;
2193   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->mergeNodes(eps,ret,ret2);
2194   if(!ret)//no nodes have been merged.
2195     return ret;
2196   std::vector<DataArrayDouble *> arrays;
2197   _time_discr->getArrays(arrays);
2198   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
2199     if(*iter)
2200       _type->renumberValuesOnNodes(epsOnVals,arr->getConstPointer(),meshC2->getNumberOfNodes(),*iter);
2201   setMesh(meshC2);
2202   return true;
2203 }
2204
2205 /*!
2206  * Merges coincident nodes of the underlying mesh. If some nodes are coincident, the
2207  * underlying mesh is replaced by a new mesh instance where the coincident nodes are
2208  * merged.<br>
2209  * In contrast to mergeNodes(), location of merged nodes is changed to be at their barycenter.
2210  *  \param [in] eps - a precision used to compare nodes of the two meshes.
2211  *  \param [in] epsOnVals - a precision used to compare field
2212  *         values at merged nodes. If the values differ more than \a epsOnVals, an
2213  *         exception is thrown.
2214  *  \return bool - \c true if some nodes have been merged and hence \a this field lies
2215  *         on another mesh.
2216  *  \throw If the mesh is of type not inheriting from MEDCouplingPointSet.
2217  *  \throw If the mesh is not well defined.
2218  *  \throw If the spatial discretization of \a this field is NULL.
2219  *  \throw If the data array is not set.
2220  *  \throw If field values at merged nodes (if any) deffer more than \a epsOnVals.
2221  */
2222 bool MEDCouplingFieldDouble::mergeNodes2(double eps, double epsOnVals)
2223 {
2224   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
2225   if(!meshC)
2226     throw INTERP_KERNEL::Exception("Invalid support mesh to apply mergeNodes on it : must be a MEDCouplingPointSet one !");
2227   if(!((const MEDCouplingFieldDiscretization *)_type))
2228     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform mergeNodes2 !");
2229   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
2230   bool ret;
2231   int ret2;
2232   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->mergeNodes2(eps,ret,ret2);
2233   if(!ret)//no nodes have been merged.
2234     return ret;
2235   std::vector<DataArrayDouble *> arrays;
2236   _time_discr->getArrays(arrays);
2237   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
2238     if(*iter)
2239       _type->renumberValuesOnNodes(epsOnVals,arr->getConstPointer(),meshC2->getNumberOfNodes(),*iter);
2240   setMesh(meshC2);
2241   return true;
2242 }
2243
2244 /*!
2245  * Removes from the underlying mesh nodes not used in any cell. If some nodes are
2246  * removed, the underlying mesh is replaced by a new mesh instance where the unused
2247  * nodes are removed.<br>
2248  *  \param [in] epsOnVals - a precision used to compare field
2249  *         values at merged nodes. If the values differ more than \a epsOnVals, an
2250  *         exception is thrown.
2251  *  \return bool - \c true if some nodes have been removed and hence \a this field lies
2252  *         on another mesh.
2253  *  \throw If the mesh is of type not inheriting from MEDCouplingPointSet.
2254  *  \throw If the mesh is not well defined.
2255  *  \throw If the spatial discretization of \a this field is NULL.
2256  *  \throw If the data array is not set.
2257  *  \throw If field values at merged nodes (if any) deffer more than \a epsOnVals.
2258  */
2259 bool MEDCouplingFieldDouble::zipCoords(double epsOnVals)
2260 {
2261   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
2262   if(!meshC)
2263     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::zipCoords : Invalid support mesh to apply zipCoords on it : must be a MEDCouplingPointSet one !");
2264   if(!((const MEDCouplingFieldDiscretization *)_type))
2265     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform zipCoords !");
2266   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
2267   int oldNbOfNodes=meshC2->getNumberOfNodes();
2268   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->zipCoordsTraducer();
2269   if(meshC2->getNumberOfNodes()!=oldNbOfNodes)
2270     {
2271       std::vector<DataArrayDouble *> arrays;
2272       _time_discr->getArrays(arrays);
2273       for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
2274         if(*iter)
2275           _type->renumberValuesOnNodes(epsOnVals,arr->getConstPointer(),meshC2->getNumberOfNodes(),*iter);
2276       setMesh(meshC2);
2277       return true;
2278     }
2279   return false;
2280 }
2281
2282 /*!
2283  * Removes duplicates of cells from the understanding mesh. If some cells are
2284  * removed, the underlying mesh is replaced by a new mesh instance where the cells
2285  * duplicates are removed.<br>
2286  *  \param [in] compType - specifies a cell comparison technique. Meaning of its
2287  *          valid values [0,1,2] is explained in the description of
2288  *          MEDCouplingPointSet::zipConnectivityTraducer() which is called by this method.
2289  *  \param [in] epsOnVals - a precision used to compare field
2290  *         values at merged cells. If the values differ more than \a epsOnVals, an
2291  *         exception is thrown.
2292  *  \return bool - \c true if some cells have been removed and hence \a this field lies
2293  *         on another mesh.
2294  *  \throw If the mesh is not an instance of MEDCouplingUMesh.
2295  *  \throw If the mesh is not well defined.
2296  *  \throw If the spatial discretization of \a this field is NULL.
2297  *  \throw If the data array is not set.
2298  *  \throw If field values at merged cells (if any) deffer more than \a epsOnVals.
2299  */
2300 bool MEDCouplingFieldDouble::zipConnectivity(int compType, double epsOnVals)
2301 {
2302   const MEDCouplingUMesh *meshC=dynamic_cast<const MEDCouplingUMesh *>(_mesh);
2303   if(!meshC)
2304     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::zipConnectivity : Invalid support mesh to apply zipCoords on it : must be a MEDCouplingPointSet one !");
2305   if(!((const MEDCouplingFieldDiscretization *)_type))
2306     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform zipConnectivity !");
2307   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> meshC2((MEDCouplingUMesh *)meshC->deepCpy());
2308   int oldNbOfCells=meshC2->getNumberOfCells();
2309   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->zipConnectivityTraducer(compType);
2310   if(meshC2->getNumberOfCells()!=oldNbOfCells)
2311     {
2312       std::vector<DataArrayDouble *> arrays;
2313       _time_discr->getArrays(arrays);
2314       for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
2315         if(*iter)
2316           _type->renumberValuesOnCells(epsOnVals,meshC,arr->getConstPointer(),meshC2->getNumberOfCells(),*iter);
2317       setMesh(meshC2);
2318       return true;
2319     }
2320   return false;
2321 }
2322
2323 /*!
2324  * This method calls MEDCouplingUMesh::buildSlice3D method. So this method makes the assumption that underlying mesh exists.
2325  * For the moment, this method is implemented for fields on cells.
2326  * 
2327  * \return a newly allocated field double containing the result that the user should deallocate.
2328  */
2329 MEDCouplingFieldDouble *MEDCouplingFieldDouble::extractSlice3D(const double *origin, const double *vec, double eps) const
2330 {
2331   const MEDCouplingMesh *mesh=getMesh();
2332   if(!mesh)
2333     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::extractSlice3D : underlying mesh is null !");
2334   if(getTypeOfField()!=ON_CELLS)
2335     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::extractSlice3D : only implemented for fields on cells !");
2336   const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh(mesh->buildUnstructured());
2337   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(false);
2338   ret->setMesh(umesh);
2339   DataArrayInt *cellIds=0;
2340   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mesh2=umesh->buildSlice3D(origin,vec,eps,cellIds);
2341   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIds2=cellIds;
2342   ret->setMesh(mesh2);
2343   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tupleIds=computeTupleIdsToSelectFromCellIds(cellIds->begin(),cellIds->end());
2344   std::vector<DataArrayDouble *> arrays;
2345   _time_discr->getArrays(arrays);
2346   int i=0;
2347   std::vector<DataArrayDouble *> newArr(arrays.size());
2348   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > newArr2(arrays.size());
2349   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++,i++)
2350     {
2351       if(*iter)
2352         {
2353           newArr2[i]=(*iter)->selectByTupleIdSafe(cellIds->begin(),cellIds->end());
2354           newArr[i]=newArr2[i];
2355         }
2356     }
2357   ret->setArrays(newArr);
2358   return ret.retn();
2359 }
2360
2361 /*!
2362  * Divides every cell of the underlying mesh into simplices (triangles in 2D and
2363  * tetrahedra in 3D). If some cells are divided, the underlying mesh is replaced by a new
2364  * mesh instance containing the simplices.<br> 
2365  *  \param [in] policy - specifies a pattern used for splitting. For its description, see
2366  *          MEDCouplingUMesh::simplexize().
2367  *  \return bool - \c true if some cells have been divided and hence \a this field lies
2368  *         on another mesh.
2369  *  \throw If \a policy has an invalid value. For valid values, see the description of 
2370  *         MEDCouplingUMesh::simplexize().
2371  *  \throw If MEDCouplingMesh::simplexize() is not applicable to the underlying mesh.
2372  *  \throw If the mesh is not well defined.
2373  *  \throw If the spatial discretization of \a this field is NULL.
2374  *  \throw If the data array is not set.
2375  */
2376 bool MEDCouplingFieldDouble::simplexize(int policy)
2377 {
2378   if(!_mesh)
2379     throw INTERP_KERNEL::Exception("No underlying mesh on this field to perform simplexize !");
2380   if(!((const MEDCouplingFieldDiscretization *)_type))
2381     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform simplexize !");
2382   int oldNbOfCells=_mesh->getNumberOfCells();
2383   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> meshC2(_mesh->deepCpy());
2384   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->simplexize(policy);
2385   int newNbOfCells=meshC2->getNumberOfCells();
2386   if(oldNbOfCells==newNbOfCells)
2387     return false;
2388   std::vector<DataArrayDouble *> arrays;
2389   _time_discr->getArrays(arrays);
2390   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
2391     if(*iter)
2392       _type->renumberValuesOnCellsR(_mesh,arr->getConstPointer(),arr->getNbOfElems(),*iter);
2393   setMesh(meshC2);
2394   return true;
2395 }
2396
2397 /*!
2398  * Creates a new MEDCouplingFieldDouble filled with the doubly contracted product of
2399  * every tensor of \a this 6-componental field.
2400  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, whose
2401  *          each tuple is calculated from the tuple <em>(t)</em> of \a this field as
2402  *          follows: \f$ t[0]^2+t[1]^2+t[2]^2+2*t[3]^2+2*t[4]^2+2*t[5]^2\f$. 
2403  *          This new field lies on the same mesh as \a this one. The caller is to delete
2404  *          this field using decrRef() as it is no more needed.
2405  *  \throw If \a this->getNumberOfComponents() != 6.
2406  *  \throw If the spatial discretization of \a this field is NULL.
2407  */
2408 MEDCouplingFieldDouble *MEDCouplingFieldDouble::doublyContractedProduct() const
2409 {
2410   if(!((const MEDCouplingFieldDiscretization *)_type))
2411     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform doublyContractedProduct !");
2412   MEDCouplingTimeDiscretization *td=_time_discr->doublyContractedProduct();
2413   td->copyTinyAttrFrom(*_time_discr);
2414   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2415   ret->setName("DoublyContractedProduct");
2416   ret->setMesh(getMesh());
2417   return ret.retn();
2418 }
2419
2420 /*!
2421  * Creates a new MEDCouplingFieldDouble filled with the determinant of a square
2422  * matrix defined by every tuple of \a this field, having either 4, 6 or 9 components.
2423  * The case of 6 components corresponds to that of the upper triangular matrix. 
2424  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, whose
2425  *          each tuple is the determinant of matrix of the corresponding tuple of \a this 
2426  *          field. This new field lies on the same mesh as \a this one. The caller is to 
2427  *          delete this field using decrRef() as it is no more needed.
2428  *  \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2429  *  \throw If the spatial discretization of \a this field is NULL.
2430  */
2431 MEDCouplingFieldDouble *MEDCouplingFieldDouble::determinant() const
2432 {
2433   if(!((const MEDCouplingFieldDiscretization *)_type))
2434     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform determinant !");
2435   MEDCouplingTimeDiscretization *td=_time_discr->determinant();
2436   td->copyTinyAttrFrom(*_time_discr);
2437   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2438   ret->setName("Determinant");
2439   ret->setMesh(getMesh());
2440   return ret.retn();
2441 }
2442
2443
2444 /*!
2445  * Creates a new MEDCouplingFieldDouble with 3 components filled with 3 eigenvalues of
2446  * an upper triangular matrix defined by every tuple of \a this 6-componental field.
2447  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, 
2448  *          having 3 components, whose each tuple contains the eigenvalues of the matrix of
2449  *          corresponding tuple of \a this field. This new field lies on the same mesh as
2450  *          \a this one. The caller is to delete this field using decrRef() as it is no
2451  *          more needed.  
2452  *  \throw If \a this->getNumberOfComponents() != 6.
2453  *  \throw If the spatial discretization of \a this field is NULL.
2454  */
2455 MEDCouplingFieldDouble *MEDCouplingFieldDouble::eigenValues() const
2456 {
2457   if(!((const MEDCouplingFieldDiscretization *)_type))
2458     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform eigenValues !");
2459   MEDCouplingTimeDiscretization *td=_time_discr->eigenValues();
2460   td->copyTinyAttrFrom(*_time_discr);
2461   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2462   ret->setName("EigenValues");
2463   ret->setMesh(getMesh());
2464   return ret.retn();
2465 }
2466
2467 /*!
2468  * Creates a new MEDCouplingFieldDouble with 9 components filled with 3 eigenvectors of
2469  * an upper triangular matrix defined by every tuple of \a this 6-componental field.
2470  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, 
2471  *          having 9 components, whose each tuple contains the eigenvectors of the matrix of
2472  *          corresponding tuple of \a this field. This new field lies on the same mesh as
2473  *          \a this one. The caller is to delete this field using decrRef() as it is no
2474  *          more needed.  
2475  *  \throw If \a this->getNumberOfComponents() != 6.
2476  *  \throw If the spatial discretization of \a this field is NULL.
2477  */
2478 MEDCouplingFieldDouble *MEDCouplingFieldDouble::eigenVectors() const
2479 {
2480   if(!((const MEDCouplingFieldDiscretization *)_type))
2481     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform eigenVectors !");
2482   MEDCouplingTimeDiscretization *td=_time_discr->eigenVectors();
2483   td->copyTinyAttrFrom(*_time_discr);
2484   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2485   ret->setName("EigenVectors");
2486   ret->setMesh(getMesh());
2487   return ret.retn();
2488 }
2489
2490 /*!
2491  * Creates a new MEDCouplingFieldDouble filled with the inverse matrix of
2492  * a matrix defined by every tuple of \a this field having either 4, 6 or 9
2493  * components. The case of 6 components corresponds to that of the upper triangular
2494  * matrix.
2495  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, 
2496  *          having the same number of components as \a this one, whose each tuple
2497  *          contains the inverse matrix of the matrix of corresponding tuple of \a this
2498  *          field. This new field lies on the same mesh as \a this one. The caller is to
2499  *          delete this field using decrRef() as it is no more needed.  
2500  *  \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2501  *  \throw If the spatial discretization of \a this field is NULL.
2502  */
2503 MEDCouplingFieldDouble *MEDCouplingFieldDouble::inverse() const
2504 {
2505   if(!((const MEDCouplingFieldDiscretization *)_type))
2506     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform inverse !");
2507   MEDCouplingTimeDiscretization *td=_time_discr->inverse();
2508   td->copyTinyAttrFrom(*_time_discr);
2509   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2510   ret->setName("Inversion");
2511   ret->setMesh(getMesh());
2512   return ret.retn();
2513 }
2514
2515 /*!
2516  * Creates a new MEDCouplingFieldDouble filled with the trace of
2517  * a matrix defined by every tuple of \a this field having either 4, 6 or 9
2518  * components. The case of 6 components corresponds to that of the upper triangular
2519  * matrix.
2520  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, 
2521  *          having 1 component, whose each tuple is the trace of the matrix of
2522  *          corresponding tuple of \a this field.
2523  *          This new field lies on the same mesh as \a this one. The caller is to
2524  *          delete this field using decrRef() as it is no more needed.  
2525  *  \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2526  *  \throw If the spatial discretization of \a this field is NULL.
2527  */
2528 MEDCouplingFieldDouble *MEDCouplingFieldDouble::trace() const
2529 {
2530   if(!((const MEDCouplingFieldDiscretization *)_type))
2531     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform trace !");
2532   MEDCouplingTimeDiscretization *td=_time_discr->trace();
2533   td->copyTinyAttrFrom(*_time_discr);
2534   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2535   ret->setName("Trace");
2536   ret->setMesh(getMesh());
2537   return ret.retn();
2538 }
2539
2540 /*!
2541  * Creates a new MEDCouplingFieldDouble filled with the stress deviator tensor of
2542  * a stress tensor defined by every tuple of \a this 6-componental field.
2543  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, 
2544  *          having same number of components and tuples as \a this field,
2545  *          whose each tuple contains the stress deviator tensor of the stress tensor of
2546  *          corresponding tuple of \a this field. This new field lies on the same mesh as
2547  *          \a this one. The caller is to delete this field using decrRef() as it is no
2548  *          more needed.  
2549  *  \throw If \a this->getNumberOfComponents() != 6.
2550  *  \throw If the spatial discretization of \a this field is NULL.
2551  */
2552 MEDCouplingFieldDouble *MEDCouplingFieldDouble::deviator() const
2553 {
2554   if(!((const MEDCouplingFieldDiscretization *)_type))
2555     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform deviator !");
2556   MEDCouplingTimeDiscretization *td=_time_discr->deviator();
2557   td->copyTinyAttrFrom(*_time_discr);
2558   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2559   ret->setName("Deviator");
2560   ret->setMesh(getMesh());
2561   return ret.retn();
2562 }
2563
2564 /*!
2565  * Creates a new MEDCouplingFieldDouble filled with the magnitude of
2566  * every vector of \a this field.
2567  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble, 
2568  *          having one component, whose each tuple is the magnitude of the vector
2569  *          of corresponding tuple of \a this field. This new field lies on the
2570  *          same mesh as \a this one. The caller is to
2571  *          delete this field using decrRef() as it is no more needed.  
2572  *  \throw If the spatial discretization of \a this field is NULL.
2573  */
2574 MEDCouplingFieldDouble *MEDCouplingFieldDouble::magnitude() const
2575 {
2576   if(!((const MEDCouplingFieldDiscretization *)_type))
2577     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform magnitude !");
2578   MEDCouplingTimeDiscretization *td=_time_discr->magnitude();
2579   td->copyTinyAttrFrom(*_time_discr);
2580   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2581   ret->setName("Magnitude");
2582   ret->setMesh(getMesh());
2583   return ret.retn();
2584 }
2585
2586 /*!
2587  * Creates a new scalar MEDCouplingFieldDouble filled with the maximal value among
2588  * values of every tuple of \a this field.
2589  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2590  *          This new field lies on the same mesh as \a this one. The caller is to
2591  *          delete this field using decrRef() as it is no more needed.  
2592  *  \throw If the spatial discretization of \a this field is NULL.
2593  */
2594 MEDCouplingFieldDouble *MEDCouplingFieldDouble::maxPerTuple() const
2595 {
2596   if(!((const MEDCouplingFieldDiscretization *)_type))
2597     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform maxPerTuple !");
2598   MEDCouplingTimeDiscretization *td=_time_discr->maxPerTuple();
2599   td->copyTinyAttrFrom(*_time_discr);
2600   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2601   std::ostringstream oss;
2602   oss << "Max_" << getName();
2603   ret->setName(oss.str());
2604   ret->setMesh(getMesh());
2605   return ret.retn();
2606 }
2607
2608 /*!
2609  * Changes number of components in \a this field. If \a newNbOfComp is less
2610  * than \a this->getNumberOfComponents() then each tuple
2611  * is truncated to have \a newNbOfComp components, keeping first components. If \a
2612  * newNbOfComp is more than \a this->getNumberOfComponents() then 
2613  * each tuple is populated with \a dftValue to have \a newNbOfComp components.  
2614  *  \param [in] newNbOfComp - number of components for the new field to have.
2615  *  \param [in] dftValue - value assigned to new values added to \a this field.
2616  *  \throw If \a this is not allocated.
2617  */
2618 void MEDCouplingFieldDouble::changeNbOfComponents(int newNbOfComp, double dftValue)
2619 {
2620   _time_discr->changeNbOfComponents(newNbOfComp,dftValue);
2621 }
2622
2623 /*!
2624  * Creates a new MEDCouplingFieldDouble composed of selected components of \a this field.
2625  * The new MEDCouplingFieldDouble has the same number of tuples but includes components
2626  * specified by \a compoIds parameter. So that getNbOfElems() of the result field
2627  * can be either less, same or more than \a this->getNumberOfValues().
2628  *  \param [in] compoIds - sequence of zero based indices of components to include
2629  *              into the new field.
2630  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble that the caller
2631  *          is to delete using decrRef() as it is no more needed.
2632  *  \throw If a component index (\a i) is not valid: 
2633  *         \a i < 0 || \a i >= \a this->getNumberOfComponents().
2634  */
2635 MEDCouplingFieldDouble *MEDCouplingFieldDouble::keepSelectedComponents(const std::vector<int>& compoIds) const
2636 {
2637   if(!((const MEDCouplingFieldDiscretization *)_type))
2638     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform keepSelectedComponents !");
2639   MEDCouplingTimeDiscretization *td=_time_discr->keepSelectedComponents(compoIds);
2640   td->copyTinyAttrFrom(*_time_discr);
2641   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2642   ret->setName(getName());
2643   ret->setMesh(getMesh());
2644   return ret.retn();
2645 }
2646
2647
2648 /*!
2649  * Copy all components in a specified order from another field.
2650  * The number of tuples in \a this and the other field can be different.
2651  *  \param [in] f - the field to copy data from.
2652  *  \param [in] compoIds - sequence of zero based indices of components, data of which is
2653  *              to be copied.
2654  *  \throw If the two fields have different number of data arrays.
2655  *  \throw If a data array is set in one of fields and is not set in the other.
2656  *  \throw If \a compoIds.size() != \a a->getNumberOfComponents().
2657  *  \throw If \a compoIds[i] < 0 or \a compoIds[i] > \a this->getNumberOfComponents().
2658  */
2659 void MEDCouplingFieldDouble::setSelectedComponents(const MEDCouplingFieldDouble *f, const std::vector<int>& compoIds)
2660 {
2661   _time_discr->setSelectedComponents(f->_time_discr,compoIds);
2662 }
2663
2664 /*!
2665  * Sorts value within every tuple of \a this field.
2666  *  \param [in] asc - if \a true, the values are sorted in ascending order, else,
2667  *              in descending order.
2668  *  \throw If a data array is not allocated.
2669  */
2670 void MEDCouplingFieldDouble::sortPerTuple(bool asc)
2671 {
2672   _time_discr->sortPerTuple(asc);
2673 }
2674
2675 /*!
2676  * Creates a new MEDCouplingFieldDouble by concatenating two given fields.
2677  * Values of
2678  * the first field precede values of the second field within the result field.
2679  *  \param [in] f1 - the first field.
2680  *  \param [in] f2 - the second field.
2681  *  \return MEDCouplingFieldDouble * - the result field. It is a new instance of
2682  *          MEDCouplingFieldDouble. The caller is to delete this mesh using decrRef() 
2683  *          as it is no more needed.
2684  *  \throw If the fields are not compatible for the merge.
2685  *  \throw If the spatial discretization of \a f1 is NULL.
2686  *  \throw If the time discretization of \a f1 is NULL.
2687  *
2688  *  \ref cpp_mcfielddouble_MergeFields "Here is a C++ example".<br>
2689  *  \ref  py_mcfielddouble_MergeFields "Here is a Python example".
2690  */
2691 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MergeFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2692 {
2693   if(!f1->areCompatibleForMerge(f2))
2694     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MergeFields on them !");
2695   const MEDCouplingMesh *m1(f1->getMesh()),*m2(f2->getMesh());
2696   if(!f1->_time_discr)
2697     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MergeFields : no time discr of f1 !");
2698   if(!f1->_type)
2699     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MergeFields : no spatial discr of f1 !");
2700   MEDCouplingTimeDiscretization *td=f1->_time_discr->aggregate(f2->_time_discr);
2701   td->copyTinyAttrFrom(*f1->_time_discr);
2702   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2703   ret->setName(f1->getName());
2704   ret->setDescription(f1->getDescription());
2705   if(m1)
2706     {
2707       MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=m1->mergeMyselfWith(m2);
2708       ret->setMesh(m);
2709     }
2710   return ret.retn();
2711 }
2712
2713 /*!
2714  * Creates a new MEDCouplingFieldDouble by concatenating all given fields.
2715  * Values of the *i*-th field precede values of the (*i*+1)-th field within the result.
2716  * If there is only one field in \a a, a deepCopy() (except time information of mesh and
2717  * field) of the field is returned. 
2718  * Generally speaking the first field in \a a is used to assign tiny attributes of the
2719  * returned field. 
2720  *  \param [in] a - a vector of fields (MEDCouplingFieldDouble) to concatenate.
2721  *  \return MEDCouplingFieldDouble * - the result field. It is a new instance of
2722  *          MEDCouplingFieldDouble. The caller is to delete this mesh using decrRef() 
2723  *          as it is no more needed.
2724  *  \throw If \a a is empty.
2725  *  \throw If the fields are not compatible for the merge.
2726  *
2727  *  \ref cpp_mcfielddouble_MergeFields "Here is a C++ example".<br>
2728  *  \ref  py_mcfielddouble_MergeFields "Here is a Python example".
2729  */
2730 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MergeFields(const std::vector<const MEDCouplingFieldDouble *>& a)
2731 {
2732   if(a.size()<1)
2733     throw INTERP_KERNEL::Exception("FieldDouble::MergeFields : size of array must be >= 1 !");
2734   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > ms(a.size());
2735   std::vector< const MEDCouplingUMesh *> ms2(a.size());
2736   std::vector< const MEDCouplingTimeDiscretization *> tds(a.size());
2737   std::vector<const MEDCouplingFieldDouble *>::const_iterator it=a.begin();
2738   const MEDCouplingFieldDouble *ref=(*it++);
2739   if(!ref)
2740     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MergeFields : presence of NULL instance in first place of input vector !");
2741   for(;it!=a.end();it++)
2742     if(!ref->areCompatibleForMerge(*it))
2743       throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MergeFields on them !");
2744   for(int i=0;i<(int)a.size();i++)
2745     {
2746       if(a[i]->getMesh())
2747         { ms[i]=a[i]->getMesh()->buildUnstructured(); ms2[i]=ms[i]; }
2748       else
2749         { ms[i]=0; ms2[i]=0; }
2750       tds[i]=a[i]->_time_discr;
2751     }
2752   MEDCouplingTimeDiscretization *td=tds[0]->aggregate(tds);
2753   td->copyTinyAttrFrom(*(a[0]->_time_discr));
2754   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(a[0]->getNature(),td,a[0]->_type->clone());
2755   ret->setName(a[0]->getName());
2756   ret->setDescription(a[0]->getDescription());
2757   if(ms2[0])
2758     {
2759       MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=MEDCouplingUMesh::MergeUMeshes(ms2);
2760       m->copyTinyInfoFrom(ms2[0]);
2761       ret->setMesh(m);
2762     }
2763   return ret.retn();
2764 }
2765
2766 /*!
2767  * Creates a new MEDCouplingFieldDouble by concatenating components of two given fields.
2768  * The number of components in the result field is a sum of the number of components of
2769  * given fields. The number of tuples in the result field is same as that of each of given
2770  * arrays.
2771  * Number of tuples in the given fields must be the same.
2772  *  \param [in] f1 - a field to include in the result field.
2773  *  \param [in] f2 - another field to include in the result field.
2774  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2775  *          The caller is to delete this result field using decrRef() as it is no more
2776  *          needed.
2777  *  \throw If the fields are not compatible for a meld (areCompatibleForMeld()).
2778  *  \throw If any of data arrays is not allocated.
2779  *  \throw If \a f1->getNumberOfTuples() != \a f2->getNumberOfTuples()
2780  */
2781 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MeldFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2782 {
2783   if(!f1->areCompatibleForMeld(f2))
2784     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MeldFields on them !");
2785   MEDCouplingTimeDiscretization *td=f1->_time_discr->meld(f2->_time_discr);
2786   td->copyTinyAttrFrom(*f1->_time_discr);
2787   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2788   ret->setMesh(f1->getMesh());
2789   return ret.retn();
2790 }
2791
2792 /*!
2793  * Returns a new MEDCouplingFieldDouble containing a dot product of two given fields, 
2794  * so that the i-th tuple of the result field is a sum of products of j-th components of
2795  * i-th tuples of given fields (\f$ f_i = \sum_{j=1}^n f1_j * f2_j \f$). 
2796  * Number of tuples and components in the given fields must be the same.
2797  *  \param [in] f1 - a given field.
2798  *  \param [in] f2 - another given field.
2799  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2800  *          The caller is to delete this result field using decrRef() as it is no more
2801  *          needed.
2802  *  \throw If either \a f1 or \a f2 is NULL.
2803  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
2804  *         differ not only in values.
2805  */
2806 MEDCouplingFieldDouble *MEDCouplingFieldDouble::DotFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2807 {
2808   if(!f1)
2809     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::DotFields : input field is NULL !");
2810   if(!f1->areStrictlyCompatible(f2))
2811     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply DotFields on them !");
2812   MEDCouplingTimeDiscretization *td=f1->_time_discr->dot(f2->_time_discr);
2813   td->copyTinyAttrFrom(*f1->_time_discr);
2814   MEDCouplingFieldDouble *ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2815   ret->setMesh(f1->getMesh());
2816   return ret;
2817 }
2818
2819 /*!
2820  * Returns a new MEDCouplingFieldDouble containing a cross product of two given fields, 
2821  * so that
2822  * the i-th tuple of the result field is a 3D vector which is a cross
2823  * product of two vectors defined by the i-th tuples of given fields.
2824  * Number of tuples in the given fields must be the same.
2825  * Number of components in the given fields must be 3.
2826  *  \param [in] f1 - a given field.
2827  *  \param [in] f2 - another given field.
2828  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2829  *          The caller is to delete this result field using decrRef() as it is no more
2830  *          needed.
2831  *  \throw If either \a f1 or \a f2 is NULL.
2832  *  \throw If \a f1->getNumberOfComponents() != 3
2833  *  \throw If \a f2->getNumberOfComponents() != 3
2834  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
2835  *         differ not only in values.
2836  */
2837 MEDCouplingFieldDouble *MEDCouplingFieldDouble::CrossProductFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2838 {
2839   if(!f1)
2840     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::CrossProductFields : input field is NULL !");
2841   if(!f1->areStrictlyCompatible(f2))
2842     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply CrossProductFields on them !");
2843   MEDCouplingTimeDiscretization *td=f1->_time_discr->crossProduct(f2->_time_discr);
2844   td->copyTinyAttrFrom(*f1->_time_discr);
2845   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2846   ret->setMesh(f1->getMesh());
2847   return ret.retn();
2848 }
2849
2850 /*!
2851  * Returns a new MEDCouplingFieldDouble containing maximal values of two given fields.
2852  * Number of tuples and components in the given fields must be the same.
2853  *  \param [in] f1 - a field to compare values with another one.
2854  *  \param [in] f2 - another field to compare values with the first one.
2855  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2856  *          The caller is to delete this result field using decrRef() as it is no more
2857  *          needed.
2858  *  \throw If either \a f1 or \a f2 is NULL.
2859  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
2860  *         differ not only in values.
2861  *
2862  *  \ref cpp_mcfielddouble_MaxFields "Here is a C++ example".<br>
2863  *  \ref  py_mcfielddouble_MaxFields "Here is a Python example".
2864  */
2865 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MaxFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2866 {
2867   if(!f1)
2868     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MaxFields : input field is NULL !");
2869   if(!f1->areStrictlyCompatible(f2))
2870     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MaxFields on them !");
2871   MEDCouplingTimeDiscretization *td=f1->_time_discr->max(f2->_time_discr);
2872   td->copyTinyAttrFrom(*f1->_time_discr);
2873   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2874   ret->setMesh(f1->getMesh());
2875   return ret.retn();
2876 }
2877
2878 /*!
2879  * Returns a new MEDCouplingFieldDouble containing minimal values of two given fields.
2880  * Number of tuples and components in the given fields must be the same.
2881  *  \param [in] f1 - a field to compare values with another one.
2882  *  \param [in] f2 - another field to compare values with the first one.
2883  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2884  *          The caller is to delete this result field using decrRef() as it is no more
2885  *          needed.
2886  *  \throw If either \a f1 or \a f2 is NULL.
2887  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
2888  *         differ not only in values.
2889  *
2890  *  \ref cpp_mcfielddouble_MaxFields "Here is a C++ example".<br>
2891  *  \ref  py_mcfielddouble_MaxFields "Here is a Python example".
2892  */
2893 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MinFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2894 {
2895   if(!f1)
2896     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MinFields : input field is NULL !");
2897   if(!f1->areStrictlyCompatible(f2))
2898     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MinFields on them !");
2899   MEDCouplingTimeDiscretization *td=f1->_time_discr->min(f2->_time_discr);
2900   td->copyTinyAttrFrom(*f1->_time_discr);
2901   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2902   ret->setMesh(f1->getMesh());
2903   return ret.retn();
2904 }
2905
2906 /*!
2907  * Returns a copy of \a this field in which sign of all values is reversed.
2908  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble
2909  *         containing the same number of tuples and components as \a this field. 
2910  *         The caller is to delete this result field using decrRef() as it is no more
2911  *         needed. 
2912  *  \throw If the spatial discretization of \a this field is NULL.
2913  *  \throw If a data array is not allocated.
2914  */
2915 MEDCouplingFieldDouble *MEDCouplingFieldDouble::negate() const
2916 {
2917   if(!((const MEDCouplingFieldDiscretization *)_type))
2918     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform negate !");
2919   MEDCouplingTimeDiscretization *td=_time_discr->negate();
2920   td->copyTinyAttrFrom(*_time_discr);
2921   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
2922   ret->setMesh(getMesh());
2923   return ret.retn();
2924 }
2925
2926 /*!
2927  * Returns a new MEDCouplingFieldDouble containing sum values of corresponding values of
2928  * two given fields ( _f_ [ i, j ] = _f1_ [ i, j ] + _f2_ [ i, j ] ).
2929  * Number of tuples and components in the given fields must be the same.
2930  *  \param [in] f1 - a field to sum up.
2931  *  \param [in] f2 - another field to sum up.
2932  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2933  *          The caller is to delete this result field using decrRef() as it is no more
2934  *          needed.
2935  *  \throw If either \a f1 or \a f2 is NULL.
2936  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
2937  *         differ not only in values.
2938  */
2939 MEDCouplingFieldDouble *MEDCouplingFieldDouble::AddFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2940 {
2941   if(!f1)
2942     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::AddFields : input field is NULL !");
2943   if(!f1->areStrictlyCompatible(f2))
2944     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply AddFields on them !");
2945   MEDCouplingTimeDiscretization *td=f1->_time_discr->add(f2->_time_discr);
2946   td->copyTinyAttrFrom(*f1->_time_discr);
2947   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2948   ret->setMesh(f1->getMesh());
2949   return ret.retn();
2950 }
2951
2952 /*!
2953  * Adds values of another MEDCouplingFieldDouble to values of \a this one
2954  * ( _this_ [ i, j ] += _other_ [ i, j ] ) using DataArrayDouble::addEqual().
2955  * The two fields must have same number of tuples, components and same underlying mesh.
2956  *  \param [in] other - the field to add to \a this one.
2957  *  \return const MEDCouplingFieldDouble & - a reference to \a this field.
2958  *  \throw If \a other is NULL.
2959  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
2960  *         differ not only in values.
2961  */
2962 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator+=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
2963 {
2964   if(!areStrictlyCompatible(&other))
2965     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply += on them !");
2966   _time_discr->addEqual(other._time_discr);
2967   return *this;
2968 }
2969
2970 /*!
2971  * Returns a new MEDCouplingFieldDouble containing subtraction of corresponding values of
2972  * two given fields ( _f_ [ i, j ] = _f1_ [ i, j ] - _f2_ [ i, j ] ).
2973  * Number of tuples and components in the given fields must be the same.
2974  *  \param [in] f1 - a field to subtract from.
2975  *  \param [in] f2 - a field to subtract.
2976  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
2977  *          The caller is to delete this result field using decrRef() as it is no more
2978  *          needed.
2979  *  \throw If either \a f1 or \a f2 is NULL.
2980  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
2981  *         differ not only in values.
2982  */
2983 MEDCouplingFieldDouble *MEDCouplingFieldDouble::SubstractFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
2984 {
2985   if(!f1)
2986     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::SubstractFields : input field is NULL !");
2987   if(!f1->areStrictlyCompatible(f2))
2988     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply SubstractFields on them !");
2989   MEDCouplingTimeDiscretization *td=f1->_time_discr->substract(f2->_time_discr);
2990   td->copyTinyAttrFrom(*f1->_time_discr);
2991   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2992   ret->setMesh(f1->getMesh());
2993   return ret.retn();
2994 }
2995
2996 /*!
2997  * Subtract values of another MEDCouplingFieldDouble from values of \a this one
2998  * ( _this_ [ i, j ] -= _other_ [ i, j ] ) using DataArrayDouble::substractEqual().
2999  * The two fields must have same number of tuples, components and same underlying mesh.
3000  *  \param [in] other - the field to subtract from \a this one.
3001  *  \return const MEDCouplingFieldDouble & - a reference to \a this field.
3002  *  \throw If \a other is NULL.
3003  *  \throw If the fields are not strictly compatible (areStrictlyCompatible()), i.e. they
3004  *         differ not only in values.
3005  */
3006 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator-=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
3007 {
3008   if(!areStrictlyCompatible(&other))
3009     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply -= on them !");
3010   _time_discr->substractEqual(other._time_discr);
3011   return *this;
3012 }
3013
3014 /*!
3015  * Returns a new MEDCouplingFieldDouble containing product values of
3016  * two given fields. There are 2 valid cases.
3017  * 1.  The fields have same number of tuples and components. Then each value of
3018  *   the result field (_f_) is a product of the corresponding values of _f1_ and
3019  *   _f2_, i.e. _f_ [ i, j ] = _f1_ [ i, j ] * _f2_ [ i, j ].
3020  * 2.  The fields have same number of tuples and one field, say _f2_, has one
3021  *   component. Then
3022  *   _f_ [ i, j ] = _f1_ [ i, j ] * _f2_ [ i, 0 ].
3023  *
3024  * The two fields must have same number of tuples and same underlying mesh.
3025  *  \param [in] f1 - a factor field.
3026  *  \param [in] f2 - another factor field.
3027  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
3028  *          The caller is to delete this result field using decrRef() as it is no more
3029  *          needed.
3030  *  \throw If either \a f1 or \a f2 is NULL.
3031  *  \throw If the fields are not compatible for production (areCompatibleForMul()),
3032  *         i.e. they differ not only in values and possibly number of components.
3033  */
3034 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MultiplyFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
3035 {
3036   if(!f1)
3037     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MultiplyFields : input field is NULL !");
3038   if(!f1->areCompatibleForMul(f2))
3039     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MultiplyFields on them !");
3040   MEDCouplingTimeDiscretization *td=f1->_time_discr->multiply(f2->_time_discr);
3041   td->copyTinyAttrFrom(*f1->_time_discr);
3042   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
3043   ret->setMesh(f1->getMesh());
3044   return ret.retn();
3045 }
3046
3047 /*!
3048  * Multiply values of another MEDCouplingFieldDouble to values of \a this one
3049  * using DataArrayDouble::multiplyEqual().
3050  * The two fields must have same number of tuples and same underlying mesh.
3051  * There are 2 valid cases.
3052  * 1.  The fields have same number of components. Then each value of
3053  *   \a other is multiplied to the corresponding value of \a this field, i.e.
3054  *   _this_ [ i, j ] *= _other_ [ i, j ].
3055  * 2. The _other_ field has one component. Then
3056  *   _this_ [ i, j ] *= _other_ [ i, 0 ].
3057  *
3058  * The two fields must have same number of tuples and same underlying mesh.
3059  *  \param [in] other - an field to multiply to \a this one.
3060  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
3061  *          The caller is to delete this result field using decrRef() as it is no more
3062  *          needed.
3063  *  \throw If \a other is NULL.
3064  *  \throw If the fields are not strictly compatible for production
3065  *         (areCompatibleForMul()),
3066  *         i.e. they differ not only in values and possibly in number of components.
3067  */
3068 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator*=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
3069 {
3070   if(!areCompatibleForMul(&other))
3071     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply *= on them !");
3072   _time_discr->multiplyEqual(other._time_discr);
3073   return *this;
3074 }
3075
3076 /*!
3077  * Returns a new MEDCouplingFieldDouble containing division of two given fields.
3078  * There are 2 valid cases.
3079  * 1.  The fields have same number of tuples and components. Then each value of
3080  *   the result field (_f_) is a division of the corresponding values of \a f1 and
3081  *   \a f2, i.e. _f_ [ i, j ] = _f1_ [ i, j ] / _f2_ [ i, j ].
3082  * 2.  The fields have same number of tuples and _f2_ has one component. Then
3083  *   _f_ [ i, j ] = _f1_ [ i, j ] / _f2_ [ i, 0 ].
3084  *
3085  *  \param [in] f1 - a numerator field.
3086  *  \param [in] f2 - a denominator field.
3087  *  \return MEDCouplingFieldDouble * - the new instance of MEDCouplingFieldDouble.
3088  *          The caller is to delete this result field using decrRef() as it is no more
3089  *          needed.
3090  *  \throw If either \a f1 or \a f2 is NULL.
3091  *  \throw If the fields are not compatible for division (areCompatibleForDiv()),
3092  *         i.e. they differ not only in values and possibly in number of components.
3093  */
3094 MEDCouplingFieldDouble *MEDCouplingFieldDouble::DivideFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
3095 {
3096   if(!f1)
3097     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::DivideFields : input field is NULL !");
3098   if(!f1->areCompatibleForDiv(f2))
3099     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply DivideFields on them !");
3100   MEDCouplingTimeDiscretization *td=f1->_time_discr->divide(f2->_time_discr);
3101   td->copyTinyAttrFrom(*f1->_time_discr);
3102   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
3103   ret->setMesh(f1->getMesh());
3104   return ret.retn();
3105 }
3106
3107 /*!
3108  * Divide values of \a this field by values of another MEDCouplingFieldDouble
3109  * using DataArrayDouble::divideEqual().
3110  * The two fields must have same number of tuples and same underlying mesh.
3111  * There are 2 valid cases.
3112  * 1.  The fields have same number of components. Then each value of
3113  *    \a this field is divided by the corresponding value of \a other one, i.e.
3114  *   _this_ [ i, j ] /= _other_ [ i, j ].
3115  * 2.  The \a other field has one component. Then
3116  *   _this_ [ i, j ] /= _other_ [ i, 0 ].
3117  *
3118  *  \warning No check of division by zero is performed!
3119  *  \param [in] other - an field to divide \a this one by.
3120  *  \throw If \a other is NULL.
3121  *  \throw If the fields are not compatible for division (areCompatibleForDiv()),
3122  *         i.e. they differ not only in values and possibly in number of components.
3123  */
3124 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator/=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
3125 {
3126   if(!areCompatibleForDiv(&other))
3127     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply /= on them !");
3128   _time_discr->divideEqual(other._time_discr);
3129   return *this;
3130 }
3131
3132 /*!
3133  * Directly called by MEDCouplingFieldDouble::operator^.
3134  * 
3135  * \sa MEDCouplingFieldDouble::operator^
3136  */
3137 MEDCouplingFieldDouble *MEDCouplingFieldDouble::PowFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2)
3138 {
3139   if(!f1)
3140     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::PowFields : input field is NULL !");
3141   if(!f1->areCompatibleForMul(f2))
3142     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply PowFields on them !");
3143   MEDCouplingTimeDiscretization *td=f1->_time_discr->pow(f2->_time_discr);
3144   td->copyTinyAttrFrom(*f1->_time_discr);
3145   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
3146   ret->setMesh(f1->getMesh());
3147   return ret.retn();
3148 }
3149
3150 /*!
3151  * Directly call MEDCouplingFieldDouble::PowFields static method.
3152  * 
3153  * \sa MEDCouplingFieldDouble::PowFields
3154  */
3155 MEDCouplingFieldDouble *MEDCouplingFieldDouble::operator^(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception)
3156 {
3157   return PowFields(this,&other);
3158 }
3159
3160 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator^=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
3161 {
3162   if(!areCompatibleForDiv(&other))
3163     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply /= on them !");
3164   _time_discr->powEqual(other._time_discr);
3165   return *this;
3166 }
3167
3168 /*!
3169  * Writes the field series \a fs and the mesh the fields lie on in the VTK file \a fileName.
3170  * If \a fs is empty no file is written.
3171  * The result file is valid provided that no exception is thrown.
3172  * \warning All the fields must be named and lie on the same non NULL mesh.
3173  *  \param [in] fileName - the name of a VTK file to write in.
3174  *  \param [in] fs - the fields to write.
3175  *  \param [in] isBinary - specifies the VTK format of the written file. By default true (Binary mode)
3176  *  \throw If \a fs[ 0 ] == NULL.
3177  *  \throw If the fields lie not on the same mesh.
3178  *  \throw If the mesh is not set.
3179  *  \throw If any of the fields has no name.
3180  *
3181  *  \ref cpp_mcfielddouble_WriteVTK "Here is a C++ example".<br>
3182  *  \ref  py_mcfielddouble_WriteVTK "Here is a Python example".
3183  */
3184 void MEDCouplingFieldDouble::WriteVTK(const std::string& fileName, const std::vector<const MEDCouplingFieldDouble *>& fs, bool isBinary)
3185 {
3186   if(fs.empty())
3187     return;
3188   std::size_t nfs=fs.size();
3189   if(!fs[0])
3190     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : 1st instance of field is NULL !");
3191   const MEDCouplingMesh *m=fs[0]->getMesh();
3192   if(!m)
3193     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : 1st instance of field lies on NULL mesh !");
3194   for(std::size_t i=1;i<nfs;i++)
3195     if(fs[i]->getMesh()!=m)
3196       throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : Fields are not lying on a same mesh ! Expected by VTK ! MEDCouplingFieldDouble::setMesh or MEDCouplingFieldDouble::changeUnderlyingMesh can help to that.");
3197   if(!m)
3198     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : Fields are lying on a same mesh but it is empty !");
3199   MEDCouplingAutoRefCountObjectPtr<DataArrayByte> byteArr;
3200   if(isBinary)
3201     { byteArr=DataArrayByte::New(); byteArr->alloc(0,1); }
3202   std::ostringstream coss,noss;
3203   for(std::size_t i=0;i<nfs;i++)
3204     {
3205       const MEDCouplingFieldDouble *cur=fs[i];
3206       std::string name(cur->getName());
3207       if(name.empty())
3208         {
3209           std::ostringstream oss; oss << "MEDCouplingFieldDouble::WriteVTK : Field in pos #" << i << " has no name !";
3210           throw INTERP_KERNEL::Exception(oss.str().c_str());
3211         }
3212       TypeOfField typ=cur->getTypeOfField();
3213       if(typ==ON_CELLS)
3214         cur->getArray()->writeVTK(coss,8,cur->getName(),byteArr);
3215       else if(typ==ON_NODES)
3216         cur->getArray()->writeVTK(noss,8,cur->getName(),byteArr);
3217       else
3218         throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : only node and cell fields supported for the moment !");
3219     }
3220   m->writeVTKAdvanced(fileName,coss.str(),noss.str(),byteArr);
3221 }
3222
3223 void MEDCouplingFieldDouble::reprQuickOverview(std::ostream& stream) const
3224 {
3225   stream << "MEDCouplingFieldDouble C++ instance at " << this << ". Name : \"" << _name << "\"." << std::endl;
3226   const char *nat=0;
3227   try
3228     {
3229       nat=MEDCouplingNatureOfField::GetRepr(_nature);
3230       stream << "Nature of field : " << nat << ".\n";
3231     }
3232   catch(INTERP_KERNEL::Exception& /*e*/)
3233     {  }
3234   const MEDCouplingFieldDiscretization *fd(_type);
3235   if(!fd)
3236     stream << "No spatial discretization set !";
3237   else
3238     fd->reprQuickOverview(stream);
3239   stream << std::endl;
3240   if(!_mesh)
3241     stream << "\nNo mesh support defined !";
3242   else
3243     {
3244       std::ostringstream oss;
3245       _mesh->reprQuickOverview(oss);
3246       std::string tmp(oss.str());
3247       stream << "\nMesh info : " << tmp.substr(0,tmp.find('\n'));
3248     }
3249   if(_time_discr)
3250     {
3251       const DataArrayDouble *arr=_time_discr->getArray();
3252       if(arr)
3253         {
3254           stream << "\n\nArray info : ";
3255           arr->reprQuickOverview(stream);
3256         }
3257       else
3258         {
3259           stream << "\n\nNo data array set !";
3260         }
3261     }
3262 }