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