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