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