Salome HOME
Merge from V6_main 12/04/2013
[modules/med.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  * Creates a new instance of MEDCouplingFieldDouble of given type. The caller is responsable for the returned field.
40  *
41  * \param [in] type type of spatial discretization of a created field (\ref ParaMEDMEM::ON_CELLS "ON_CELLS", \ref ParaMEDMEM::ON_NODES "ON_NODES", \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT", \ref ParaMEDMEM::ON_GAUSS_NE "ON_GAUSS_NE", \ref ParaMEDMEM::ON_NODES_KR "ON_NODES_KR").
42  * \param [in] td type of time discretization of a created field (\ref ParaMEDMEM::NO_TIME "NO_TIME", \ref ParaMEDMEM::ONE_TIME "ONE_TIME", \ref ParaMEDMEM::LINEAR_TIME "LINEAR_TIME", \ref ParaMEDMEM::CONST_ON_TIME_INTERVAL "CONST_ON_TIME_INTERVAL").
43  
44  * \return a newly allocated field the caller should deal with.
45  */
46 MEDCouplingFieldDouble *MEDCouplingFieldDouble::New(TypeOfField type, TypeOfTimeDiscretization td)
47 {
48   return new MEDCouplingFieldDouble(type,td);
49 }
50
51 /*!
52  * Creates a new instance of MEDCouplingFieldDouble of given type. The caller is responsable for the returned field.
53  * ** WARINING : This method do not deeply copy neither mesh nor spatial discretization. Only a shallow copy (reference) is done for mesh and spatial discretization ! **
54  *
55  * \param [in] ft \ref MEDCouplingFieldTemplatesPage "field template" defining its spatial discretization and supporting mesh.
56  * \param [in] td type of time discretization of a created field (\ref ParaMEDMEM::NO_TIME "NO_TIME", \ref ParaMEDMEM::ONE_TIME "ONE_TIME", \ref ParaMEDMEM::LINEAR_TIME "LINEAR_TIME", \ref ParaMEDMEM::CONST_ON_TIME_INTERVAL "CONST_ON_TIME_INTERVAL")
57  
58  * \return a newly allocated field the caller should deal with.
59  */
60 MEDCouplingFieldDouble *MEDCouplingFieldDouble::New(const MEDCouplingFieldTemplate& ft, TypeOfTimeDiscretization td)
61 {
62   return new MEDCouplingFieldDouble(ft,td);
63 }
64
65 /*!
66  * Sets time \a unit of \a this field.
67  *
68  * \param [in] unit \a unit (string) in which time is measured.
69  */
70 void MEDCouplingFieldDouble::setTimeUnit(const char *unit)
71 {
72   _time_discr->setTimeUnit(unit);
73 }
74
75 /*!
76  * Returns a time unit of \a this field.
77  *
78  * \return a string describing units in which time is measured.
79  */
80 const char *MEDCouplingFieldDouble::getTimeUnit() const
81 {
82   return _time_discr->getTimeUnit();
83 }
84
85 /*!
86  * This method if possible the time information (time unit, time iteration, time unit and time value) with its support
87  * that is to say its mesh.
88  * 
89  * \throw  If \c this->_mesh is null an exception will be thrown. An exception will also be throw if the spatial discretization is
90  *         NO_TIME.
91  */
92 void MEDCouplingFieldDouble::synchronizeTimeWithSupport() throw(INTERP_KERNEL::Exception)
93 {
94   _time_discr->synchronizeTimeWith(_mesh);
95 }
96
97 /*!
98  * This method performs a copy of \a this **without any copy of the underlying mesh** ( see warning section of this method).
99  * The copy of arrays is deep if \b recDeepCpy equals to true, no copy of arrays is done if \b recDeepCpy equals to false.
100  *
101  * \c clone(false) is rather dedicated for advanced users that want to limit the amount of memory.
102  * 
103  * It allows the user to perform methods
104  * MEDCouplingFieldDouble::AddFields, MEDCouplingFieldDouble::MultiplyFields with \a this and the returned field.
105  * 
106  * \warning The \b underlying \b mesh of the returned field is \b always the same (same pointer) than \a this **whatever the value** of \a recDeepCpy parameter.
107  * If the user wants to duplicated deeply the underlying mesh he should call MEDCouplingFieldDouble::cloneWithMesh method or MEDCouplingFieldDouble::deepCpy instead.
108  *
109  * \param [in] recDeepCpy specifies if underlying arrays in \a this should be copied or only attached to the returned field.
110  * \return a newly allocated MEDCouplingFieldDouble instance that the caller should deal with.
111  * \sa ParaMEDMEM::MEDCouplingFieldDouble::cloneWithMesh(bool recDeepCpy) const
112  */
113 MEDCouplingFieldDouble *MEDCouplingFieldDouble::clone(bool recDeepCpy) const
114 {
115   return new MEDCouplingFieldDouble(*this,recDeepCpy);
116 }
117
118 /*!
119  * This method behaves exactly like MEDCouplingFieldDouble::clone method **except that here the underlying mesh is systematically **
120  * (whatever the value of the input parameter \a recDeepCpy) **deeply duplicated**.
121  *
122  * The result of \c cloneWithMesh(true) is exactly the same than calling \ref MEDCouplingFieldDouble::deepCpy "deepCpy".
123  * 
124  * So the resulting field of this call cannot be called with \a this with the following methods MEDCouplingFieldDouble::AddFields, MEDCouplingFieldDouble::MultiplyFields ...
125  * To avoid to deep copy the underlying mesh the user should call MEDCouplingFieldDouble::clone method instead.
126
127  * \param [in] recDeepCpy specifies if underlying arrays in \a this should be copied or only attached to the returned field.
128  * \return a newly allocated MEDCouplingFieldDouble instance that the caller should deal with.
129  * \sa ParaMEDMEM::MEDCouplingFieldDouble::clone(bool recDeepCpy) const
130  */
131 MEDCouplingFieldDouble *MEDCouplingFieldDouble::cloneWithMesh(bool recDeepCpy) const
132 {
133   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(recDeepCpy);
134   if(_mesh)
135     {
136       MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> mCpy=_mesh->deepCpy();
137       ret->setMesh(mCpy);
138     }
139   return ret.retn();
140 }
141
142 /*!
143  * This method performs a deepCpy of \a this (**mesh included**)!
144  * So the resulting field of this call cannot be called with \a this with following methods MEDCouplingFieldDouble::AddFields, MEDCouplingFieldDouble::MultiplyFields ...
145  * To avoid deep copying the underlying mesh the user should call MEDCouplingFieldDouble::clone method instead.
146  * This method is exactly equivalent to MEDCouplingFieldDouble::cloneWithMesh called with parameter true.
147  *
148  * \return a newly allocated MEDCouplingFieldDouble instance that the caller should deal with.
149  * \sa ParaMEDMEM::MEDCouplingFieldDouble::cloneWithMesh(bool recDeepCpy) const
150  */
151 MEDCouplingFieldDouble *MEDCouplingFieldDouble::deepCpy() const
152 {
153   return cloneWithMesh(true);
154 }
155
156 /*!
157  * TODOC
158  * 
159  * \param [in] td type of time discretization of a created field (\ref ParaMEDMEM::NO_TIME "NO_TIME", \ref ParaMEDMEM::ONE_TIME "ONE_TIME", \ref ParaMEDMEM::LINEAR_TIME "LINEAR_TIME", \ref ParaMEDMEM::CONST_ON_TIME_INTERVAL "CONST_ON_TIME_INTERVAL").
160  * \param [in] deepCopy specifies if underlying arrays in \a this should be copied or only attached to the returned field.
161  * \return a newly allocated MEDCouplingFieldDouble instance that the caller should deal with.
162  *
163  * \ref cpp_mcfielddouble_buildnewtimereprfromthis "Here a C++ example."
164  
165  * \ref py_mcfielddouble_buildnewtimereprfromthis "Here a Python example."
166  * \sa ParaMEDMEM::MEDCouplingFieldDouble::clone(bool recDeepCpy) const
167  */
168 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildNewTimeReprFromThis(TypeOfTimeDiscretization td, bool deepCopy) const
169 {
170   MEDCouplingTimeDiscretization *tdo=_time_discr->buildNewTimeReprFromThis(td,deepCopy);
171   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretization> disc;
172   if(_type)
173     disc=_type->clone();
174   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),tdo,disc.retn());
175   ret->setMesh(getMesh());
176   ret->setName(getName());
177   ret->setDescription(getDescription());
178   return ret.retn();
179 }
180
181 /*!
182  * Copy tiny info (component names, name, description) but warning the underlying mesh is not renamed (for safety reason).
183  */
184 void MEDCouplingFieldDouble::copyTinyStringsFrom(const MEDCouplingField *other) throw(INTERP_KERNEL::Exception)
185 {
186   MEDCouplingField::copyTinyStringsFrom(other);
187   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
188   if(otherC)
189     {
190       _time_discr->copyTinyStringsFrom(*otherC->_time_discr);
191     }
192 }
193
194 /*!
195  * Copy only times, order, iteration from other. The underlying mesh is not impacted by this method.
196  * Arrays are not impacted too.
197  */
198 void MEDCouplingFieldDouble::copyTinyAttrFrom(const MEDCouplingFieldDouble *other) throw(INTERP_KERNEL::Exception)
199 {
200   if(other)
201     {
202       _time_discr->copyTinyAttrFrom(*other->_time_discr);
203     }
204   
205 }
206
207 void MEDCouplingFieldDouble::copyAllTinyAttrFrom(const MEDCouplingFieldDouble *other) throw(INTERP_KERNEL::Exception)
208 {
209   copyTinyStringsFrom(other);
210   copyTinyAttrFrom(other);
211 }
212
213 std::string MEDCouplingFieldDouble::simpleRepr() const
214 {
215   std::ostringstream ret;
216   ret << "FieldDouble with name : \"" << getName() << "\"\n";
217   ret << "Description of field is : \"" << getDescription() << "\"\n";
218   if(_type)
219     { ret << "FieldDouble space discretization is : " << _type->getStringRepr() << "\n"; }
220   else
221     { ret << "FieldDouble has no spatial discretization !\n"; }
222   if(_time_discr)
223     { ret << "FieldDouble time discretization is : " << _time_discr->getStringRepr() << "\n"; }
224   else
225     { ret << "FieldDouble has no time discretization !\n"; }
226   ret << "FieldDouble nature of field is : \"" << MEDCouplingNatureOfField::GetReprNoThrow(_nature) << "\"\n";
227   if(getArray())
228     {
229       if(getArray()->isAllocated())
230         {
231           int nbOfCompo=getArray()->getNumberOfComponents();
232           ret << "FieldDouble default array has " << nbOfCompo << " components and " << getArray()->getNumberOfTuples() << " tuples.\n";
233           ret << "FieldDouble default array has following info on components : ";
234           for(int i=0;i<nbOfCompo;i++)
235             ret << "\"" << getArray()->getInfoOnComponent(i) << "\" ";
236           ret << "\n";
237         }
238       else
239         {
240           ret << "Array set but not allocated !\n";
241         }
242     }
243   if(_mesh)
244     ret << "Mesh support information :\n__________________________\n" << _mesh->simpleRepr();
245   else
246     ret << "Mesh support information : No mesh set !\n";
247   return ret.str();
248 }
249
250 std::string MEDCouplingFieldDouble::advancedRepr() const
251 {
252   std::ostringstream ret;
253   ret << "FieldDouble with name : \"" << getName() << "\"\n";
254   ret << "Description of field is : \"" << getDescription() << "\"\n";
255   if(_type)
256     { ret << "FieldDouble space discretization is : " << _type->getStringRepr() << "\n"; }
257   else
258     { ret << "FieldDouble has no space discretization set !\n"; }
259   if(_time_discr)
260     { ret << "FieldDouble time discretization is : " << _time_discr->getStringRepr() << "\n"; }
261   else
262     { ret << "FieldDouble has no time discretization set !\n"; }
263   if(getArray())
264     ret << "FieldDouble default array has " << getArray()->getNumberOfComponents() << " components and " << getArray()->getNumberOfTuples() << " tuples.\n";
265   if(_mesh)
266     ret << "Mesh support information :\n__________________________\n" << _mesh->advancedRepr();
267   else
268     ret << "Mesh support information : No mesh set !\n";
269   std::vector<DataArrayDouble *> arrays;
270   _time_discr->getArrays(arrays);
271   int arrayId=0;
272   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++,arrayId++)
273     {
274       ret << "Array #" << arrayId << " :\n__________\n";
275       if(*iter)
276         (*iter)->reprWithoutNameStream(ret);
277       else
278         ret << "Array empty !";
279       ret << "\n";
280     }
281   return ret.str();
282 }
283
284 void MEDCouplingFieldDouble::writeVTK(const char *fileName) const throw(INTERP_KERNEL::Exception)
285 {
286   std::vector<const MEDCouplingFieldDouble *> fs(1,this);
287   MEDCouplingFieldDouble::WriteVTK(fileName,fs);
288 }
289
290 bool MEDCouplingFieldDouble::isEqualIfNotWhy(const MEDCouplingField *other, double meshPrec, double valsPrec, std::string& reason) const throw(INTERP_KERNEL::Exception)
291 {
292   if(!other)
293     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::isEqualIfNotWhy : other instance is NULL !");
294   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
295   if(!otherC)
296     {
297       reason="field given in input is not castable in MEDCouplingFieldDouble !";
298       return false;
299     }
300   if(!MEDCouplingField::isEqualIfNotWhy(other,meshPrec,valsPrec,reason))
301     return false;
302   if(!_time_discr->isEqualIfNotWhy(otherC->_time_discr,valsPrec,reason))
303     {
304       reason.insert(0,"In FieldDouble time discretizations differ :");
305       return false;
306     }
307   return true;
308 }
309
310 bool MEDCouplingFieldDouble::isEqualWithoutConsideringStr(const MEDCouplingField *other, double meshPrec, double valsPrec) const
311 {
312   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
313   if(!otherC)
314     return false;
315   if(!MEDCouplingField::isEqualWithoutConsideringStr(other,meshPrec,valsPrec))
316     return false;
317   if(!_time_discr->isEqualWithoutConsideringStr(otherC->_time_discr,valsPrec))
318     return false;
319   return true;
320 }
321
322 /*!
323  * This method states if \a this and 'other' are compatibles each other before performing any treatment.
324  * This method is good for methods like : mergeFields.
325  * This method is not very demanding compared to areStrictlyCompatible that is better for operation on fields.
326  */
327 bool MEDCouplingFieldDouble::areCompatibleForMerge(const MEDCouplingField *other) const
328 {
329   if(!MEDCouplingField::areCompatibleForMerge(other))
330     return false;
331   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
332   if(!otherC)
333     return false;
334   if(!_time_discr->areCompatible(otherC->_time_discr))
335     return false;
336   return true;
337 }
338
339 /*!
340  * This method is more strict than MEDCouplingField::areCompatibleForMerge method.
341  * This method is used for operation on fields to operate a first check before attempting operation.
342  */
343 bool MEDCouplingFieldDouble::areStrictlyCompatible(const MEDCouplingField *other) const
344 {
345   std::string tmp;
346   if(!MEDCouplingField::areStrictlyCompatible(other))
347     return false;
348   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
349   if(!otherC)
350     return false;
351   if(!_time_discr->areStrictlyCompatible(otherC->_time_discr,tmp))
352     return false;
353   return true;
354 }
355
356 /*!
357  * Method with same principle than MEDCouplingFieldDouble::areStrictlyCompatible method except that
358  * number of components between \a this and 'other' can be different here (for operator*).
359  */
360 bool MEDCouplingFieldDouble::areCompatibleForMul(const MEDCouplingField *other) const
361 {
362   if(!MEDCouplingField::areStrictlyCompatible(other))
363     return false;
364   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
365   if(!otherC)
366     return false;
367   if(!_time_discr->areStrictlyCompatibleForMul(otherC->_time_discr))
368     return false;
369   return true;
370 }
371
372 /*!
373  * Method with same principle than MEDCouplingFieldDouble::areStrictlyCompatible method except that
374  * number of components between \a this and 'other' can be different here (for operator/).
375  */
376 bool MEDCouplingFieldDouble::areCompatibleForDiv(const MEDCouplingField *other) const
377 {
378   if(!MEDCouplingField::areStrictlyCompatible(other))
379     return false;
380   const MEDCouplingFieldDouble *otherC=dynamic_cast<const MEDCouplingFieldDouble *>(other);
381   if(!otherC)
382     return false;
383   if(!_time_discr->areStrictlyCompatibleForDiv(otherC->_time_discr))
384     return false;
385   return true;
386 }
387
388 /*!
389  * This method is invocated before any attempt of melding. This method is very close to areStrictlyCompatible,
390  * except that \a this and other can have different number of components.
391  */
392 bool MEDCouplingFieldDouble::areCompatibleForMeld(const MEDCouplingFieldDouble *other) const
393 {
394   if(!MEDCouplingField::areStrictlyCompatible(other))
395     return false;
396   if(!_time_discr->areCompatibleForMeld(other->_time_discr))
397     return false;
398   return true;
399 }
400
401 /*!
402  * This method performs a clone of mesh and a renumbering of underlying cells of it. The number of cells remains the same.
403  * The values of field are impacted in consequence to have the same geometrical field.
404  */
405 void MEDCouplingFieldDouble::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
406 {
407   renumberCellsWithoutMesh(old2NewBg,check);
408   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=_mesh->deepCpy();
409   m->renumberCells(old2NewBg,check);
410   setMesh(m);
411   updateTime();
412 }
413
414 /*!
415  * \b WARNING : use this method with lot of care !
416  * This method performs half job of MEDCouplingFieldDouble::renumberCells. That is to say no permutation of cells is done on underlying mesh.
417  * That is to say, the field content is changed by this method. The reason of this method is only for multi-field instances lying on the same mesh to
418  * avoid a systematic duplication and renumbering of _mesh attribute.
419  */
420 void MEDCouplingFieldDouble::renumberCellsWithoutMesh(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
421 {
422    if(!_mesh)
423      throw INTERP_KERNEL::Exception("Expecting a defined mesh to be able to operate a renumbering !");
424    if(!((const MEDCouplingFieldDiscretization *)_type))
425      throw INTERP_KERNEL::Exception("Expecting a spatial discretization to be able to operate a renumbering !");
426   //
427   _type->renumberCells(old2NewBg,check);
428   std::vector<DataArrayDouble *> arrays;
429   _time_discr->getArrays(arrays);
430   _type->renumberArraysForCell(_mesh,arrays,old2NewBg,check);
431   //
432   updateTime();
433 }
434
435 /*!
436  * This method performs a clone of mesh and a renumbering of underlying nodes of it. The number of nodes remains not compulsory the same as renumberCells method.
437  * The values of field are impacted in consequence to have the same geometrical field.
438  * 
439  * \sa MEDCouplingFieldDouble::renumberNodesWithoutMesh
440  */
441 void MEDCouplingFieldDouble::renumberNodes(const int *old2NewBg, double eps) throw(INTERP_KERNEL::Exception)
442 {
443   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
444   if(!meshC)
445     throw INTERP_KERNEL::Exception("Invalid mesh to apply renumberNodes on it !");
446   int nbOfNodes=meshC->getNumberOfNodes();
447   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
448   int newNbOfNodes=*std::max_element(old2NewBg,old2NewBg+nbOfNodes)+1;
449   renumberNodesWithoutMesh(old2NewBg,newNbOfNodes,eps);
450   meshC2->renumberNodes(old2NewBg,newNbOfNodes);
451   setMesh(meshC2);
452 }
453
454 /*!
455  * \b WARNING : use this method with lot of care !
456  * ** WARNING : in case of throw the content in array can be partially modified until the exception raises **
457  * This method performs half job of MEDCouplingFieldDouble::renumberNodes. That is to say no permutation of nodes is done on underlying mesh.
458  * That is to say, the field content is changed by this method. As the API suggests, this method can performs the half job of nodes contraction.
459  * That's why an epsilon is given to specify a threshold of error in case of two nodes are merged but the difference of values on these nodes are higher than \a eps.
460  */
461 void MEDCouplingFieldDouble::renumberNodesWithoutMesh(const int *old2NewBg, int newNbOfNodes, double eps) throw(INTERP_KERNEL::Exception)
462 {
463   if(!((const MEDCouplingFieldDiscretization *)_type))
464     throw INTERP_KERNEL::Exception("Expecting a spatial discretization to be able to operate a renumbering !");
465   std::vector<DataArrayDouble *> arrays;
466   _time_discr->getArrays(arrays);
467   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
468     if(*iter)
469       _type->renumberValuesOnNodes(eps,old2NewBg,newNbOfNodes,*iter);
470 }
471
472 /*!
473  * This method makes the assumption that the default array is set. If not an exception will be thrown.
474  * This method is usable only if the default array has exactly one component. If not an exception will be thrown too.
475  * This method returns all tuples ids that fit the range [vmin,vmax].
476  * The caller has the responsability of the returned DataArrayInt.
477  */
478 DataArrayInt *MEDCouplingFieldDouble::getIdsInRange(double vmin, double vmax) const throw(INTERP_KERNEL::Exception)
479 {
480   if(getArray()==0)
481     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getIdsInRange : no default array set !");
482   return getArray()->getIdsInRange(vmin,vmax);
483 }
484
485 /*!
486  * Builds a newly created field, that the caller will have the responsability to deal with (decrRef).
487  * This method makes the assumption that the field is correctly defined when this method is called, no check of this will be done.
488  * This method returns a restriction of \a this so that only tuples id specified in 'part' will be contained in returned field.
489  * Parameter 'part' specifies \b cell \b ids \b whatever \b the \b spatial \b discretization of \a this (ON_CELLS, ON_NODES, ON_GAUSS_PT, ON_GAUSS_NE)
490  *
491  * If \a this is a field on cell lying on a mesh that have 10 cells. If part contains following cellIds [3,7,6].
492  * In this case the returned field will lie on mesh having 3 cells and the returned field will contain 3 tuples.
493  * Tuple#0 of return field will refer to the cell#0 of returned mesh. The cell #0 of returned mesh will be equal to the cell#3 of 'this->getMesh()'
494  * Tuple#1 of return field will refer to the cell#1 of returned mesh. The cell #1 of returned mesh will be equal to the cell#7 of 'this->getMesh()'
495  * Tuple#2 of return field will refer to the cell#2 of returned mesh. The cell #2 of returned mesh will be equal to the cell#6 of 'this->getMesh()'
496  *
497  * If \a this is field on node lying on a mesh that have 10 cells and 11 nodes for example. If part contains following cellIds [3,7,6].
498  * \a this is currently contains 11 tuples. If the restriction of mesh to 3 cells leads to a mesh with 6 nodes, the returned field,
499  * will contain 6 tuples and this field will lie on this restricted mesh.
500  *
501  * \sa MEDCouplingFieldDouble::buildSubPartRange
502  */
503 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildSubPart(const DataArrayInt *part) const throw(INTERP_KERNEL::Exception)
504 {
505   if(part==0)
506     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::buildSubPart : not empty array must be passed to this method !");
507   return buildSubPart(part->begin(),part->end());
508 }
509
510 /*!
511  * Builds a newly created field, that the caller will have the responsability to deal with.
512  * \n This method makes the assumption that the field \a this is correctly defined when this method is called (\c this->checkCoherency() returns without any exception thrown), **no check of this will be done**.
513  * \n This method returns a restriction of \a this so that only tuples id specified in [ \a partBg , \a partEnd ) will be contained in returned field. 
514  * \n Parameter [\a partBg, \a partEnd ) specifies \b cell \b ids \b whatever \b the \b spatial \b discretization of \a this
515  * (\ref ParaMEDMEM::ON_CELLS "ON_CELLS", \ref ParaMEDMEM::ON_NODES "ON_CELLS", \ref ParaMEDMEM::ON_GAUSS_PT "ON_GAUSS_PT", \ref ParaMEDMEM::ON_GAUSS_NE "ON_GAUSS_NE")
516  *
517  * If \a this is a field on cell lying on a mesh that have 10 cells. If part contains following cellIds [3,7,6].
518  * In this case the returned field will lie on mesh having 3 cells and the returned field will contain 3 tuples.
519  *
520  *- Tuple#0 of return field will refer to the cell#0 of returned mesh. The cell #0 of returned mesh will be equal to the cell#3 of \c this->getMesh()
521  *- Tuple#1 of return field will refer to the cell#1 of returned mesh. The cell #1 of returned mesh will be equal to the cell#7 of \c this->getMesh()
522  *- Tuple#2 of return field will refer to the cell#2 of returned mesh. The cell #2 of returned mesh will be equal to the cell#6 of \c this->getMesh()
523  *
524  * If \a this is field on node lying on a mesh that have 10 cells and 11 nodes for example. So \a this is currently contains 11 tuples.
525  * \n If part contains following cellIds [3,7,6].
526  * \n If the restriction of mesh to 3 cells leads to a mesh with 6 nodes, the returned field,
527  * will contain 6 tuples (and same number of components \c this->getArray()->getNumberOfComponents() ) and this field will lie on this restricted mesh.
528  *
529  * \param [in] partBg start (included) of input range cell ids to select [ \a partBg, \a partEnd )
530  * \param [in] partEnd end (not included) of input range cell ids to select [ \a partBg, \a partEnd )
531  * \return a newly allocated field the caller should deal with.
532  * 
533  * \throw if there is presence of an invalid cell id in [ \a partBg, \a partEnd ) regarding the number of cells of \c this->getMesh()
534  *
535  * \ref cpp_mcfielddouble_subpart1 "Here a C++ example."
536  *
537  * \ref py_mcfielddouble_subpart1 "Here a Python example."
538  * \sa ParaMEDMEM::MEDCouplingFieldDouble::buildSubPart(const DataArrayInt *) const, MEDCouplingFieldDouble::buildSubPartRange
539  */
540 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildSubPart(const int *partBg, const int *partEnd) const throw(INTERP_KERNEL::Exception)
541 {
542   if(!((const MEDCouplingFieldDiscretization *)_type))
543     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::buildSubPart : Expecting a not NULL spatial discretization !");
544   DataArrayInt *arrSelect;
545   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=_type->buildSubMeshData(_mesh,partBg,partEnd,arrSelect);
546   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arrSelect2(arrSelect);
547   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(false);//quick shallow copy.
548   const MEDCouplingFieldDiscretization *disc=getDiscretization();
549   if(disc)
550     ret->setDiscretization(MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretization>(disc->clonePart(partBg,partEnd)));
551   ret->setMesh(m);
552   std::vector<DataArrayDouble *> arrays;
553   _time_discr->getArrays(arrays);
554   std::vector<DataArrayDouble *> arrs;
555   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > arrsSafe;
556   const int *arrSelBg=arrSelect->begin();
557   const int *arrSelEnd=arrSelect->end();
558   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
559     {
560       DataArrayDouble *arr=0;
561       if(*iter)
562         arr=(*iter)->selectByTupleIdSafe(arrSelBg,arrSelEnd);
563       arrs.push_back(arr); arrsSafe.push_back(arr);
564     }
565   ret->_time_discr->setArrays(arrs,0);
566   return ret.retn();
567 }
568
569 /*!
570  * This method is equivalent to MEDCouplingFieldDouble::buildSubPart, the only difference is that the input range of cell ids is
571  * given using a range given \a begin, \a end and \a step to optimize the part computation.
572  * 
573  * \sa MEDCouplingFieldDouble::buildSubPart
574  */
575 MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildSubPartRange(int begin, int end, int step) const throw(INTERP_KERNEL::Exception)
576 {
577   if(!((const MEDCouplingFieldDiscretization *)_type))
578     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::buildSubPart : Expecting a not NULL spatial discretization !");
579   DataArrayInt *arrSelect;
580   int beginOut,endOut,stepOut;
581   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=_type->buildSubMeshDataRange(_mesh,begin,end,step,beginOut,endOut,stepOut,arrSelect);
582   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arrSelect2(arrSelect);
583   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(false);//quick shallow copy.
584   const MEDCouplingFieldDiscretization *disc=getDiscretization();
585   if(disc)
586     ret->setDiscretization(MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretization>(disc->clonePartRange(begin,end,step)));
587   ret->setMesh(m);
588   std::vector<DataArrayDouble *> arrays;
589   _time_discr->getArrays(arrays);
590   std::vector<DataArrayDouble *> arrs;
591   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > arrsSafe;
592   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
593     {
594       DataArrayDouble *arr=0;
595       if(*iter)
596         {
597           if(arrSelect)
598             {
599               const int *arrSelBg=arrSelect->begin();
600               const int *arrSelEnd=arrSelect->end();
601               arr=(*iter)->selectByTupleIdSafe(arrSelBg,arrSelEnd);
602             }
603           else
604             arr=(*iter)->selectByTupleId2(beginOut,endOut,stepOut);
605         }
606       arrs.push_back(arr); arrsSafe.push_back(arr);
607     }
608   ret->_time_discr->setArrays(arrs,0);
609   return ret.retn();
610 }
611
612 TypeOfTimeDiscretization MEDCouplingFieldDouble::getTimeDiscretization() const
613 {
614   return _time_discr->getEnum();
615 }
616
617 MEDCouplingFieldDouble::MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td):MEDCouplingField(type),
618                                                                                               _time_discr(MEDCouplingTimeDiscretization::New(td))
619 {
620 }
621
622 /*!
623  * ** WARINING : This method do not deeply copy neither mesh nor spatial discretization. Only a shallow copy (reference) is done for mesh and spatial discretization ! **
624  */
625 MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldTemplate& ft, TypeOfTimeDiscretization td):MEDCouplingField(ft,false),
626                                                                                                                 _time_discr(MEDCouplingTimeDiscretization::New(td))
627 {
628 }
629
630 MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCopy):MEDCouplingField(other,deepCopy),
631                                                                                                    _time_discr(other._time_discr->performCpy(deepCopy))
632 {
633 }
634
635 MEDCouplingFieldDouble::MEDCouplingFieldDouble(NatureOfField n, MEDCouplingTimeDiscretization *td, MEDCouplingFieldDiscretization *type):MEDCouplingField(type,n),_time_discr(td)
636 {
637 }
638
639 MEDCouplingFieldDouble::~MEDCouplingFieldDouble()
640 {
641   delete _time_discr;
642 }
643
644 void MEDCouplingFieldDouble::checkCoherency() const throw(INTERP_KERNEL::Exception)
645 {
646   if(!_mesh)
647     throw INTERP_KERNEL::Exception("Field invalid because no mesh specified !");
648   if(!((const MEDCouplingFieldDiscretization *)_type))
649     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::checkCoherency : no spatial discretization !");
650   _time_discr->checkCoherency();
651   _type->checkCoherencyBetween(_mesh,getArray());
652 }
653
654 /*!
655  * Returns the accumulation (the sum) of comId_th component of each tuples of \b default and \b only \b default array.
656  */
657 double MEDCouplingFieldDouble::accumulate(int compId) const
658 {
659   if(getArray()==0)
660     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::accumulate : no default array defined !");
661   return getArray()->accumulate(compId);
662 }
663
664 /*!
665  * Returns the accumulation (the sum) of all tuples of \b default and \b only default array.
666  * The res is expected to be of size getNumberOfComponents().
667  */
668 void MEDCouplingFieldDouble::accumulate(double *res) const
669 {
670   if(getArray()==0)
671     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::accumulate : no default array defined !");
672   getArray()->accumulate(res);
673 }
674
675 /*!
676  * This method returns the max value in \a this. \a this is expected to be a field with exactly \b one component. If not an exception will be thrown.
677  * To getMaxValue on vector field applyFunc is needed before. This method looks only on all arrays stored in 'this->_time_discr'.
678  * If no arrays exists, an exception will be thrown.
679  */
680 double MEDCouplingFieldDouble::getMaxValue() const throw(INTERP_KERNEL::Exception)
681 {
682   std::vector<DataArrayDouble *> arrays;
683   _time_discr->getArrays(arrays);
684   double ret=-std::numeric_limits<double>::max();
685   bool isExistingArr=false;
686   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
687     {
688       if(*iter)
689         {
690           isExistingArr=true;
691           int loc;
692           ret=std::max(ret,(*iter)->getMaxValue(loc));
693         }
694     }
695   if(!isExistingArr)
696     throw INTERP_KERNEL::Exception("getMaxValue : No arrays defined !");
697   return ret;
698 }
699
700 /*!
701  * This method is an extension of ParaMEDMEM::MEDCouplingFieldDouble::getMaxValue method because the returned 
702  * value is the same but this method also returns to you a tupleIds object which the caller have the responsibility
703  * to deal with. The main difference is that the returned tupleIds is those corresponding the first set array.
704  * If you have more than one array set (in LINEAR_TIME instance for example) only the first not null array will be used
705  * to compute tupleIds.
706  */
707 double MEDCouplingFieldDouble::getMaxValue2(DataArrayInt*& tupleIds) const throw(INTERP_KERNEL::Exception)
708 {
709   std::vector<DataArrayDouble *> arrays;
710   _time_discr->getArrays(arrays);
711   double ret=-std::numeric_limits<double>::max();
712   bool isExistingArr=false;
713   tupleIds=0;
714   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret1;
715   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
716     {
717       if(*iter)
718         {
719           isExistingArr=true;
720           DataArrayInt *tmp;
721           ret=std::max(ret,(*iter)->getMaxValue2(tmp));
722           MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmpSafe(tmp);
723           if(!((const DataArrayInt *)ret1))
724             ret1=tmpSafe;
725         }
726     }
727   if(!isExistingArr)
728     throw INTERP_KERNEL::Exception("getMaxValue2 : No arrays defined !");
729   tupleIds=ret1.retn();
730   return ret;
731 }
732
733 /*!
734  * This method returns the min value in \a this. \a this is expected to be a field with exactly \b one component. If not an exception will be thrown.
735  * To getMinValue on vector field applyFunc is needed before. This method looks only on all arrays stored in 'this->_time_discr'.
736  * If no arrays exists, an exception will be thrown.
737  */
738 double MEDCouplingFieldDouble::getMinValue() const throw(INTERP_KERNEL::Exception)
739 {
740   std::vector<DataArrayDouble *> arrays;
741   _time_discr->getArrays(arrays);
742   double ret=std::numeric_limits<double>::max();
743   bool isExistingArr=false;
744   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
745     {
746       if(*iter)
747         {
748           isExistingArr=true;
749           int loc;
750           ret=std::min(ret,(*iter)->getMinValue(loc));
751         }
752     }
753   if(!isExistingArr)
754     throw INTERP_KERNEL::Exception("getMinValue : No arrays defined !");
755   return ret;
756 }
757
758 /*!
759  * This method is an extension of ParaMEDMEM::MEDCouplingFieldDouble::getMinValue method because the returned 
760  * value is the same but this method also returns to you a tupleIds object which the caller have the responsibility
761  * to deal with. The main difference is that the returned tupleIds is those corresponding the first set array.
762  * If you have more than one array set (in LINEAR_TIME instance for example) only the first not null array will be used
763  * to compute tupleIds.
764  */
765 double MEDCouplingFieldDouble::getMinValue2(DataArrayInt*& tupleIds) const throw(INTERP_KERNEL::Exception)
766 {
767   std::vector<DataArrayDouble *> arrays;
768   _time_discr->getArrays(arrays);
769   double ret=-std::numeric_limits<double>::max();
770   bool isExistingArr=false;
771   tupleIds=0;
772   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret1;
773   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
774     {
775       if(*iter)
776         {
777           isExistingArr=true;
778           DataArrayInt *tmp;
779           ret=std::max(ret,(*iter)->getMinValue2(tmp));
780           MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmpSafe(tmp);
781           if(!((const DataArrayInt *)ret1))
782             ret1=tmpSafe;
783         }
784     }
785   if(!isExistingArr)
786     throw INTERP_KERNEL::Exception("getMinValue2 : No arrays defined !");
787   tupleIds=ret1.retn();
788   return ret;
789 }
790
791 /*!
792  * This method returns the average value in \a this. \a this is expected to be a field with exactly \b one component. If not an exception will be thrown.
793  * To getAverageValue on vector field applyFunc is needed before. This method looks only \b default array \b and \b only \b default.
794  * If default array does not exist, an exception will be thrown.
795  */
796 double MEDCouplingFieldDouble::getAverageValue() const throw(INTERP_KERNEL::Exception)
797 {
798   if(getArray()==0)
799     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getAverageValue : no default array defined !");
800   return getArray()->getAverageValue();
801 }
802
803 /*!
804  * This method returns the euclidean norm of \a this.
805  * \f[
806  * \sqrt{\sum_{0 \leq i < nbOfEntity}val[i]*val[i]}
807  * \f]
808  * If default array does not exist, an exception will be thrown.
809  */
810 double MEDCouplingFieldDouble::norm2() const throw(INTERP_KERNEL::Exception)
811 {
812   if(getArray()==0)
813     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::norm2 : no default array defined !");
814   return getArray()->norm2();
815 }
816
817 /*!
818  * This method returns the max norm of \a this.
819  * \f[
820  * \max_{0 \leq i < nbOfEntity}{abs(val[i])}
821  * \f]
822  * If default array does not exist, an exception will be thrown.
823  */
824 double MEDCouplingFieldDouble::normMax() const throw(INTERP_KERNEL::Exception)
825 {
826   if(getArray()==0)
827     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::normMax : no default array defined !");
828   return getArray()->normMax();
829 }
830
831 /*!
832  * This method returns the average value in \a this weighted by ParaMEDMEM::MEDCouplingField::buildMeasureField.
833  * \a this is expected to be a field with exactly \b one component. If not an exception will be thrown.
834  * To getAverageValue on vector field applyFunc is needed before. This method looks only \b default array \b and \b only \b default.
835  * If default array does not exist, an exception will be thrown.
836  * 
837  * \param [out] res the location where the result will be stored. \a res is expected to be a location with \c this->getNumberOfComponents() places available.
838  * \param [in] isWAbs specifies if abs is applied on measure on underlying mesh before performing computation. For a user already sure that all cells of its underlying mesh
839  *                    are all well oriented this parameter can be set to false to be 'faster'. By default this parameter is true.
840  */
841 void MEDCouplingFieldDouble::getWeightedAverageValue(double *res, bool isWAbs) const throw(INTERP_KERNEL::Exception)
842 {
843   if(getArray()==0)
844     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getWeightedAverageValue : no default array defined !");
845   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> w=buildMeasureField(isWAbs);
846   double deno=w->getArray()->accumulate(0);
847   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=getArray()->deepCpy();
848   arr->multiplyEqual(w->getArray());
849   std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::multiplies<double>(),1./deno));
850   arr->accumulate(res);
851 }
852
853 /*!
854  * This method returns the average value in \a this weighted by ParaMEDMEM::MEDCouplingField::buildMeasureField.
855  * \a this is expected to be a field with exactly \b one component. If not an exception will be thrown.
856  * To getAverageValue on vector field applyFunc is needed before. This method looks only \b default array \b and \b only \b default.
857  * If default array does not exist, an exception will be thrown.
858  * 
859  * \param [in] compId The component id that should be in [0, \c this->getNumberOfComponents() ). If not an INTERP_KERNEL::Exception will be thrown.
860  * \param [in] isWAbs specifies if abs is applied on measure on underlying mesh before performing computation. For a user already sure that all cells of its underlying mesh
861  *                    are all well oriented this parameter can be set to false to be 'faster'. By default this parameter is true in C++ not in python (overloading confusion).
862  */
863 double MEDCouplingFieldDouble::getWeightedAverageValue(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception)
864 {
865   int nbComps=getArray()->getNumberOfComponents();
866   if(compId<0 || compId>=nbComps)
867     {
868       std::ostringstream oss; oss << "MEDCouplingFieldDouble::getWeightedAverageValue : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
869       throw INTERP_KERNEL::Exception(oss.str().c_str());
870     }
871   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
872   getWeightedAverageValue(res,isWAbs);
873   return res[compId];
874 }
875
876 /*!
877  * Returns the normL1 of current field on compId component :
878  * \f[
879  * \frac{\sum_{0 \leq i < nbOfEntity}|val[i]*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}
880  * \f]
881  * If compId>=nbOfComponent an exception is thrown.
882  */
883 double MEDCouplingFieldDouble::normL1(int compId) const throw(INTERP_KERNEL::Exception)
884 {
885   if(!_mesh)
886     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1 !");
887   if(!((const MEDCouplingFieldDiscretization *)_type))
888     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL1 !");
889   int nbComps=getArray()->getNumberOfComponents();
890   if(compId<0 || compId>=nbComps)
891     {
892       std::ostringstream oss; oss << "MEDCouplingFieldDouble::normL1 : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
893       throw INTERP_KERNEL::Exception(oss.str().c_str());
894     }
895   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
896   _type->normL1(_mesh,getArray(),res);
897   return res[compId];
898 }
899
900 /*!
901  * Returns the normL1 of current field on each components :
902  * \f[
903  * \frac{\sum_{0 \leq i < nbOfEntity}|val[i]*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}
904  * \f]
905  * The res is expected to be of size getNumberOfComponents().
906  */
907 void MEDCouplingFieldDouble::normL1(double *res) const throw(INTERP_KERNEL::Exception)
908 {
909   if(!_mesh)
910     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1");
911   if(!((const MEDCouplingFieldDiscretization *)_type))
912     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL1 !");
913   _type->normL1(_mesh,getArray(),res);
914 }
915
916 /*!
917  * Returns the normL2 of current field on compId component :
918  * \f[
919  * \sqrt{\frac{\sum_{0 \leq i < nbOfEntity}|val[i]^{2}*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}}
920  * \f]
921  * If compId>=nbOfComponent an exception is thrown.
922  */
923 double MEDCouplingFieldDouble::normL2(int compId) const throw(INTERP_KERNEL::Exception)
924 {
925   if(!_mesh)
926     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL2");
927   if(!((const MEDCouplingFieldDiscretization *)_type))
928     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL2 !");
929   int nbComps=getArray()->getNumberOfComponents();
930   if(compId<0 || compId>=nbComps)
931     {
932       std::ostringstream oss; oss << "MEDCouplingFieldDouble::normL2 : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
933       throw INTERP_KERNEL::Exception(oss.str().c_str());
934     }
935   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
936   _type->normL2(_mesh,getArray(),res);
937   return res[compId];
938 }
939
940 /*!
941  * Returns the normL2 of current field on each components :
942  * \f[
943  * \sqrt{\frac{\sum_{0 \leq i < nbOfEntity}|val[i]^{2}*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}}
944  * \f]
945  * The res is expected to be of size getNumberOfComponents().
946  */
947 void MEDCouplingFieldDouble::normL2(double *res) const throw(INTERP_KERNEL::Exception)
948 {
949   if(!_mesh)
950     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL2");
951   if(!((const MEDCouplingFieldDiscretization *)_type))
952     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform normL2 !");
953   _type->normL2(_mesh,getArray(),res);
954 }
955
956 /*!
957  * Returns the accumulation (the sum) of comId_th component of each tuples weigthed by the field
958  * returns by getWeightingField relative of the _type of field of default array.
959  * This method is useful to check the conservativity of interpolation method.
960  */
961 double MEDCouplingFieldDouble::integral(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception)
962 {
963   if(!_mesh)
964     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform integral");
965   if(!((const MEDCouplingFieldDiscretization *)_type))
966     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform integral !");
967   int nbComps=getArray()->getNumberOfComponents();
968   if(compId<0 || compId>=nbComps)
969     {
970       std::ostringstream oss; oss << "MEDCouplingFieldDouble::integral : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !";
971       throw INTERP_KERNEL::Exception(oss.str().c_str());
972     }
973   INTERP_KERNEL::AutoPtr<double> res=new double[nbComps];
974   _type->integral(_mesh,getArray(),isWAbs,res);
975   return res[compId];
976 }
977
978 /*!
979  * Returns the accumulation (the sum) of each tuples weigthed by the field
980  * returns by getWeightingField relative of the _type of field of default array.
981  * This method is useful to check the conservativity of interpolation method.
982  */
983 void MEDCouplingFieldDouble::integral(bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
984 {
985   if(!_mesh)
986     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform integral2");
987   if(!((const MEDCouplingFieldDiscretization *)_type))
988     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform integral2 !");
989   _type->integral(_mesh,getArray(),isWAbs,res);
990 }
991
992 /*!
993  * This method is reserved for field lying on structured mesh spatial support. It returns the value of cell localized by (i,j,k)
994  * If spatial support is not structured mesh an exception will be thrown.
995  * @param res out array expected to be equal to size getNumberOfComponents()
996  */
997 void MEDCouplingFieldDouble::getValueOnPos(int i, int j, int k, double *res) const throw(INTERP_KERNEL::Exception)
998 {
999   const DataArrayDouble *arr=_time_discr->getArray();
1000   if(!_mesh)
1001     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOnPos");
1002   if(!((const MEDCouplingFieldDiscretization *)_type))
1003     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOnPos !");
1004   _type->getValueOnPos(arr,_mesh,i,j,k,res);
1005 }
1006
1007 /*!
1008  * Returns value of \a this on default time of point 'spaceLoc' using spatial discretization.
1009  * If 'point' is outside the spatial discretization of this an exception will be thrown.
1010  */
1011 void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception)
1012 {
1013   const DataArrayDouble *arr=_time_discr->getArray();
1014   if(!_mesh)
1015     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOn");
1016   if(!((const MEDCouplingFieldDiscretization *)_type))
1017     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOnPos !");
1018   _type->getValueOn(arr,_mesh,spaceLoc,res);
1019 }
1020
1021 /*!
1022  * Returns a newly allocated array with 'nbOfPoints' tuples and nb of components equal to 'this->getNumberOfComponents()'.
1023  */
1024 DataArrayDouble *MEDCouplingFieldDouble::getValueOnMulti(const double *spaceLoc, int nbOfPoints) const throw(INTERP_KERNEL::Exception)
1025 {
1026   const DataArrayDouble *arr=_time_discr->getArray();
1027   if(!_mesh)
1028     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOnMulti");
1029   if(!((const MEDCouplingFieldDiscretization *)_type))
1030     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOnMulti !");
1031   return _type->getValueOnMulti(arr,_mesh,spaceLoc,nbOfPoints);
1032 }
1033
1034 /*!
1035  * Returns value of \a this on time 'time' of point 'spaceLoc' using spatial discretization.
1036  * If 'time' is not covered by this->_time_discr an exception will be thrown.
1037  * If 'point' is outside the spatial discretization of this an exception will be thrown.
1038  */
1039 void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception)
1040 {
1041   std::vector< const DataArrayDouble *> arrs=_time_discr->getArraysForTime(time);
1042   if(!_mesh)
1043     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform getValueOn");
1044   if(!((const MEDCouplingFieldDiscretization *)_type))
1045     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getValueOn !");
1046   std::vector<double> res2;
1047   for(std::vector< const DataArrayDouble *>::const_iterator iter=arrs.begin();iter!=arrs.end();iter++)
1048     {
1049       int sz=(int)res2.size();
1050       res2.resize(sz+(*iter)->getNumberOfComponents());
1051       _type->getValueOn(*iter,_mesh,spaceLoc,&res2[sz]);
1052     }
1053   _time_discr->getValueForTime(time,res2,res);
1054 }
1055
1056 /*!
1057  * Applies a*x+b on 'compoId'th component of each cell.
1058  */
1059 void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId)
1060 {
1061   _time_discr->applyLin(a,b,compoId);
1062 }
1063
1064 /*!
1065  * This method sets \a this to a uniform scalar field with one component.
1066  * All tuples will have the same value 'value'.
1067  * An exception is thrown if no underlying mesh is defined.
1068  */
1069 MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator=(double value) throw(INTERP_KERNEL::Exception)
1070 {
1071   if(!_mesh)
1072     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::operator= : no mesh defined !");
1073   if(!((const MEDCouplingFieldDiscretization *)_type))
1074     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform operator = !");
1075   int nbOfTuple=_type->getNumberOfTuples(_mesh);
1076   _time_discr->setOrCreateUniformValueOnAllComponents(nbOfTuple,value);
1077   return *this;
1078 }
1079
1080 /*!
1081  * This method is very similar to this one MEDCouplingMesh::fillFromAnalytic.
1082  * See MEDCouplingMesh::fillFromAnalytic method doc to have more details.
1083  * The main difference is that the field as been started to be constructed here.
1084  * An exception is thrown if no underlying mesh is set before the call of this method.
1085  */
1086 void MEDCouplingFieldDouble::fillFromAnalytic(int nbOfComp, FunctionToEvaluate func) throw(INTERP_KERNEL::Exception)
1087 {
1088   if(!_mesh)
1089     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic : no mesh defined !");
1090   if(!((const MEDCouplingFieldDiscretization *)_type))
1091     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic !");
1092   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1093   _time_discr->fillFromAnalytic(loc,nbOfComp,func);
1094 }
1095
1096 /*!
1097  * This method is very similar to this one MEDCouplingMesh::fillFromAnalytic.
1098  * See MEDCouplingMesh::fillFromAnalytic method doc to have more details.
1099  * The main difference is that the field as been started to be constructed here.
1100  * An exception is thrown if no underlying mesh is set before the call of this method.
1101  */
1102 void MEDCouplingFieldDouble::fillFromAnalytic(int nbOfComp, const char *func) throw(INTERP_KERNEL::Exception)
1103 {
1104   if(!_mesh)
1105     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic : no mesh defined !");
1106   if(!((const MEDCouplingFieldDiscretization *)_type))
1107     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic !");
1108   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1109   _time_discr->fillFromAnalytic(loc,nbOfComp,func);
1110 }
1111
1112 /*!
1113  * This method is very similar to this one MEDCouplingMesh::fillFromAnalytic2.
1114  * See MEDCouplingMesh::fillFromAnalytic method doc to have more details.
1115  * The main difference is that the field as been started to be constructed here.
1116  * An exception is throw if no underlying mesh is set before the call of this method.
1117  */
1118 void MEDCouplingFieldDouble::fillFromAnalytic2(int nbOfComp, const char *func) throw(INTERP_KERNEL::Exception)
1119 {
1120   if(!_mesh)
1121     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic2 : no mesh defined !");
1122   if(!((const MEDCouplingFieldDiscretization *)_type))
1123     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic2 !");
1124   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1125   _time_discr->fillFromAnalytic2(loc,nbOfComp,func);
1126 }
1127
1128 /*!
1129  * This method is very similar to this one MEDCouplingMesh::fillFromAnalytic3.
1130  * See MEDCouplingMesh::fillFromAnalytic method doc to have more details.
1131  * The main difference is that the field as been started to be constructed here.
1132  * An exception is thrown if no underlying mesh is set before the call of this method.
1133  */
1134 void MEDCouplingFieldDouble::fillFromAnalytic3(int nbOfComp, const std::vector<std::string>& varsOrder, const char *func) throw(INTERP_KERNEL::Exception)
1135 {
1136   if(!_mesh)
1137     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::fillFromAnalytic2 : no mesh defined !");
1138   if(!((const MEDCouplingFieldDiscretization *)_type))
1139     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform fillFromAnalytic3 !");
1140   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> loc=_type->getLocalizationOfDiscValues(_mesh);
1141   _time_discr->fillFromAnalytic3(loc,nbOfComp,varsOrder,func);
1142 }
1143
1144 /*!
1145  * Applyies the function specified by pointer 'func' on each tuples on all arrays contained in _time_discr.
1146  * If '*func' returns false during one evaluation an exception will be thrown.
1147  */
1148 void MEDCouplingFieldDouble::applyFunc(int nbOfComp, FunctionToEvaluate func)
1149 {
1150   _time_discr->applyFunc(nbOfComp,func);
1151 }
1152
1153 /*!
1154  * This method is a specialization of other overloaded methods. When 'nbOfComp' equals 1 this method is equivalent to
1155  * ParaMEDMEM::MEDCouplingFieldDouble::operator=.
1156  */
1157 void MEDCouplingFieldDouble::applyFunc(int nbOfComp, double val)
1158 {
1159   if(!_mesh)
1160     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::applyFunc : no mesh defined !");
1161   if(!((const MEDCouplingFieldDiscretization *)_type))
1162     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform applyFunc !");
1163   int nbOfTuple=_type->getNumberOfTuples(_mesh);
1164   _time_discr->setUniformValue(nbOfTuple,nbOfComp,val);
1165 }
1166
1167 /*!
1168  * Applyies the function specified by the string repr 'func' on each tuples on all arrays contained in _time_discr.
1169  * If '*func' fails in evaluation during one evaluation an exception will be thrown.
1170  * The field will contain 'nbOfComp' components after the call.
1171  */
1172 void MEDCouplingFieldDouble::applyFunc(int nbOfComp, const char *func) throw(INTERP_KERNEL::Exception)
1173 {
1174   _time_discr->applyFunc(nbOfComp,func);
1175 }
1176
1177 /*!
1178  * This method is equivalent to MEDCouplingFieldDouble::applyFunc, except that here components info are used to determine variables position in 'func'.
1179  * If there is vars detected in 'func' that is not in an info on components an exception will be thrown.
1180  */
1181 void MEDCouplingFieldDouble::applyFunc2(int nbOfComp, const char *func) throw(INTERP_KERNEL::Exception)
1182 {
1183   _time_discr->applyFunc2(nbOfComp,func);
1184 }
1185
1186 /*!
1187  * This method is equivalent to MEDCouplingFieldDouble::applyFunc, except that here 'varsOrder' is used to determine variables position in 'func'.
1188  * If there is vars detected in 'func' that is not in 'varsOrder' an exception will be thrown.
1189  */
1190 void MEDCouplingFieldDouble::applyFunc3(int nbOfComp, const std::vector<std::string>& varsOrder, const char *func) throw(INTERP_KERNEL::Exception)
1191 {
1192   _time_discr->applyFunc3(nbOfComp,varsOrder,func);
1193 }
1194
1195 /*!
1196  * Applyies the function specified by the string repr 'func' on each tuples on all arrays contained in _time_discr.
1197  * If '*func' fails in evaluation during one evaluation an exception will be thrown.
1198  * The field will contain exactly the same number of components after the call.
1199  */
1200 void MEDCouplingFieldDouble::applyFunc(const char *func) throw(INTERP_KERNEL::Exception)
1201 {
1202   _time_discr->applyFunc(func);
1203 }
1204
1205 /*!
1206  * Applyies the function specified by the string repr 'func' on each tuples on all arrays contained in _time_discr.
1207  * The field will contain exactly the same number of components after the call.
1208  * Use is not warranted for the moment !
1209  */
1210 void MEDCouplingFieldDouble::applyFuncFast32(const char *func) throw(INTERP_KERNEL::Exception)
1211 {
1212   _time_discr->applyFuncFast32(func);
1213 }
1214
1215 /*!
1216  * Applyies the function specified by the string repr 'func' on each tuples on all arrays contained in _time_discr.
1217  * The field will contain exactly the same number of components after the call.
1218  * Use is not warranted for the moment !
1219  */
1220 void MEDCouplingFieldDouble::applyFuncFast64(const char *func) throw(INTERP_KERNEL::Exception)
1221 {
1222   _time_discr->applyFuncFast64(func);
1223 }
1224
1225 /*!
1226  * This method makes the assumption that the default array has been set before.
1227  * If not an exception will be sent.
1228  * If default array set, the number of components will be sent.
1229  */
1230 int MEDCouplingFieldDouble::getNumberOfComponents() const throw(INTERP_KERNEL::Exception)
1231 {
1232   if(getArray()==0)
1233     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getNumberOfComponents : No array specified !");
1234   return getArray()->getNumberOfComponents();
1235 }
1236
1237 /*!
1238  * This method makes the assumption that _mesh has be set before the call of this method and description of gauss
1239  * localizations in case of Gauss field. If not an exception will sent.
1240  * \b Contrary to MEDCouplingFieldDouble::getNumberOfComponents and MEDCouplingFieldDouble::getNumberOfValues is
1241  * \b not aware of the presence of the default array.
1242  * \b WARNING \b no coherency check is done here. MEDCouplingFieldDouble::checkCoherency method should be called to check that !
1243  */
1244 int MEDCouplingFieldDouble::getNumberOfTuples() const throw(INTERP_KERNEL::Exception)
1245 {
1246   if(!_mesh)
1247     throw INTERP_KERNEL::Exception("Impossible to retrieve number of tuples because no mesh specified !");
1248   if(!((const MEDCouplingFieldDiscretization *)_type))
1249     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getNumberOfTuples !");
1250   return _type->getNumberOfTuples(_mesh);
1251 }
1252
1253 /*!
1254  * This method makes the assumption that the default array has been set before.
1255  * If not an exception will be sent.
1256  * If default array set, the number of values present in the default array will be sent.
1257  */
1258 int MEDCouplingFieldDouble::getNumberOfValues() const throw(INTERP_KERNEL::Exception)
1259 {
1260   if(getArray()==0)
1261     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getNumberOfValues : No array specified !");
1262   return getArray()->getNbOfElems();
1263 }
1264
1265 void MEDCouplingFieldDouble::updateTime() const
1266 {
1267   MEDCouplingField::updateTime();
1268   updateTimeWith(*_time_discr);
1269 }
1270
1271 std::size_t MEDCouplingFieldDouble::getHeapMemorySize() const
1272 {
1273   std::size_t ret=0;
1274   if(_time_discr)
1275     ret+=_time_discr->getHeapMemorySize();
1276   return MEDCouplingField::getHeapMemorySize()+ret;
1277 }
1278
1279 void MEDCouplingFieldDouble::setNature(NatureOfField nat) throw(INTERP_KERNEL::Exception)
1280 {
1281   MEDCouplingField::setNature(nat);
1282   if(_type)
1283     _type->checkCompatibilityWithNature(nat);
1284 }
1285
1286 /*!
1287  * This method synchronizes time information (time, iteration, order, time unit) regarding the information in \c this->_mesh.
1288  * \throw If no mesh is set in this. Or if \a this is not compatible with time setting (typically NO_TIME)
1289  */
1290 void MEDCouplingFieldDouble::synchronizeTimeWithMesh() throw(INTERP_KERNEL::Exception)
1291 {
1292   if(!_mesh)
1293     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::synchronizeTimeWithMesh : no mesh set in this !");
1294   int it=-1,ordr=-1;
1295   double val=_mesh->getTime(it,ordr);
1296   std::string timeUnit(_mesh->getTimeUnit());
1297   setTime(val,it,ordr);
1298   setTimeUnit(timeUnit.c_str());
1299 }
1300
1301 double MEDCouplingFieldDouble::getIJK(int cellId, int nodeIdInCell, int compoId) const
1302 {
1303   if(!((const MEDCouplingFieldDiscretization *)_type))
1304     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getIJK !");
1305   return _type->getIJK(_mesh,getArray(),cellId,nodeIdInCell,compoId);
1306 }
1307
1308 void MEDCouplingFieldDouble::setArray(DataArrayDouble *array)
1309 {
1310   _time_discr->setArray(array,this);
1311 }
1312
1313 void MEDCouplingFieldDouble::setEndArray(DataArrayDouble *array)
1314 {
1315   _time_discr->setEndArray(array,this);
1316 }
1317
1318 void MEDCouplingFieldDouble::setArrays(const std::vector<DataArrayDouble *>& arrs) throw(INTERP_KERNEL::Exception)
1319 {
1320   _time_discr->setArrays(arrs,this);
1321 }
1322
1323 void MEDCouplingFieldDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
1324 {
1325   tinyInfo.clear();
1326   _time_discr->getTinySerializationStrInformation(tinyInfo);
1327   tinyInfo.push_back(_name);
1328   tinyInfo.push_back(_desc);
1329   tinyInfo.push_back(getTimeUnit());
1330 }
1331
1332 /*!
1333  * This method retrieves some critical values to resize and prepare remote instance.
1334  * The first two elements returned in tinyInfo correspond to the parameters to give in constructor.
1335  * @param tinyInfo out parameter resized correctly after the call. The length of this vector is tiny.
1336  */
1337 void MEDCouplingFieldDouble::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
1338 {
1339   if(!((const MEDCouplingFieldDiscretization *)_type))
1340     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getTinySerializationIntInformation !");
1341   tinyInfo.clear();
1342   tinyInfo.push_back((int)_type->getEnum());
1343   tinyInfo.push_back((int)_time_discr->getEnum());
1344   tinyInfo.push_back((int)_nature);
1345   _time_discr->getTinySerializationIntInformation(tinyInfo);
1346   std::vector<int> tinyInfo2;
1347   _type->getTinySerializationIntInformation(tinyInfo2);
1348   tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
1349   tinyInfo.push_back((int)tinyInfo2.size());
1350 }
1351
1352 /*!
1353  * This method retrieves some critical values to resize and prepare remote instance.
1354  * @param tinyInfo out parameter resized correctly after the call. The length of this vector is tiny.
1355  */
1356 void MEDCouplingFieldDouble::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
1357 {
1358   if(!((const MEDCouplingFieldDiscretization *)_type))
1359     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform getTinySerializationDbleInformation !");
1360   tinyInfo.clear();
1361   _time_discr->getTinySerializationDbleInformation(tinyInfo);
1362   std::vector<double> tinyInfo2;
1363   _type->getTinySerializationDbleInformation(tinyInfo2);
1364   tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
1365   tinyInfo.push_back((int)tinyInfo2.size());//very bad, lack of time to improve it
1366 }
1367
1368 /*!
1369  * This method has to be called to the new instance filled by CORBA, MPI, File...
1370  * @param tinyInfoI is the value retrieves from distant result of getTinySerializationIntInformation on source instance to be copied.
1371  * @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.
1372  * @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.
1373  *               No decrRef must be applied to every instances in returned vector.
1374  */
1375 void MEDCouplingFieldDouble::resizeForUnserialization(const std::vector<int>& tinyInfoI, DataArrayInt *&dataInt, std::vector<DataArrayDouble *>& arrays)
1376 {
1377   if(!((const MEDCouplingFieldDiscretization *)_type))
1378     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform resizeForUnserialization !");
1379   dataInt=0;
1380   std::vector<int> tinyInfoITmp(tinyInfoI);
1381   int sz=tinyInfoITmp.back();
1382   tinyInfoITmp.pop_back();
1383   std::vector<int> tinyInfoITmp2(tinyInfoITmp.begin(),tinyInfoITmp.end()-sz);
1384   std::vector<int> tinyInfoI2(tinyInfoITmp2.begin()+3,tinyInfoITmp2.end());
1385   _time_discr->resizeForUnserialization(tinyInfoI2,arrays);
1386   std::vector<int> tinyInfoITmp3(tinyInfoITmp.end()-sz,tinyInfoITmp.end());
1387   _type->resizeForUnserialization(tinyInfoITmp3,dataInt);
1388 }
1389
1390 void MEDCouplingFieldDouble::finishUnserialization(const std::vector<int>& tinyInfoI, const std::vector<double>& tinyInfoD, const std::vector<std::string>& tinyInfoS)
1391 {
1392   if(!((const MEDCouplingFieldDiscretization *)_type))
1393     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform finishUnserialization !");
1394   std::vector<int> tinyInfoI2(tinyInfoI.begin()+3,tinyInfoI.end());
1395   //
1396   std::vector<double> tmp(tinyInfoD);
1397   int sz=(int)tinyInfoD.back();//very bad, lack of time to improve it
1398   tmp.pop_back();
1399   std::vector<double> tmp1(tmp.begin(),tmp.end()-sz);
1400   std::vector<double> tmp2(tmp.end()-sz,tmp.end());
1401   //
1402   _time_discr->finishUnserialization(tinyInfoI2,tmp1,tinyInfoS);
1403   _nature=(NatureOfField)tinyInfoI[2];
1404   _type->finishUnserialization(tmp2);
1405   int nbOfElemS=(int)tinyInfoS.size();
1406   _name=tinyInfoS[nbOfElemS-3];
1407   _desc=tinyInfoS[nbOfElemS-2];
1408   setTimeUnit(tinyInfoS[nbOfElemS-1].c_str());
1409 }
1410
1411 /*!
1412  * Contrary to MEDCouplingPointSet class the returned arrays are \b not the responsabilities of the caller.
1413  * The values returned must be consulted only in readonly mode.
1414  */
1415 void MEDCouplingFieldDouble::serialize(DataArrayInt *&dataInt, std::vector<DataArrayDouble *>& arrays) const
1416 {
1417   if(!((const MEDCouplingFieldDiscretization *)_type))
1418     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform serialize !");
1419   _time_discr->getArrays(arrays);
1420   _type->getSerializationIntArray(dataInt);
1421 }
1422
1423 /*!
1424  * This method tries to to change the mesh support of \a this following the parameter 'levOfCheck' and 'precOnMesh'.
1425  * Semantic of 'levOfCheck' is explained in MEDCouplingMesh::checkGeoEquivalWith method. This method is used to perform the job.
1426  * If this->_mesh is not defined or other an exeption will be throw.
1427  */
1428 void MEDCouplingFieldDouble::changeUnderlyingMesh(const MEDCouplingMesh *other, int levOfCheck, double precOnMesh, double eps) throw(INTERP_KERNEL::Exception)
1429 {
1430   if(_mesh==0 || other==0)
1431     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::changeUnderlyingMesh : is expected to operate on not null meshes !");
1432   DataArrayInt *cellCor=0,*nodeCor=0;
1433   other->checkGeoEquivalWith(_mesh,levOfCheck,precOnMesh,cellCor,nodeCor);
1434   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellCor2(cellCor),nodeCor2(nodeCor);
1435   if(cellCor)
1436     renumberCellsWithoutMesh(cellCor->getConstPointer(),false);
1437   if(nodeCor)
1438     renumberNodesWithoutMesh(nodeCor->getConstPointer(),nodeCor->getMaxValueInArray()+1,eps);
1439   setMesh(const_cast<MEDCouplingMesh *>(other));
1440 }
1441
1442 /*!
1443  * This method is an extension of MEDCouplingFieldDouble::operator-=. It allows a user to operate a difference of 2 fields (\a this and 'f') even if they do not share same meshes.
1444  * No interpolation will be done here only an analyze of two underlying mesh will be done to see if the meshes are geometrically equivalent. If yes, the eventual renumbering will be done and operator-= applyed after.
1445  * This method requires that 'f' and \a this are coherent (check coherency) and that 'f' and \a this would be coherent for a merge.
1446  * Semantic of 'levOfCheck' is explained in MEDCouplingMesh::checkGeoEquivalWith method.
1447  * \param [in] precOnMesh precision for the mesh comparison between \c this->getMesh() and \c f->getMesh()
1448  * \param [in] eps the precision on values that can appear in case of 
1449  */
1450 void MEDCouplingFieldDouble::substractInPlaceDM(const MEDCouplingFieldDouble *f, int levOfCheck, double precOnMesh, double eps) throw(INTERP_KERNEL::Exception)
1451 {
1452   checkCoherency();
1453   if(!f)
1454     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::substractInPlaceDM : input field is NULL !");
1455   f->checkCoherency();
1456   if(!areCompatibleForMerge(f))
1457     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::substractInPlaceDM : Fields are not compatible ; unable to apply mergeFields on them !");
1458   changeUnderlyingMesh(f->getMesh(),levOfCheck,precOnMesh,eps);
1459   operator-=(*f);
1460 }
1461
1462 /*!
1463  * Merge nodes of underlying mesh. In case of some node will be merged the underlying mesh instance will change.
1464  * The first 'eps' stands for geometric approximation. The second 'epsOnVals' is for epsilon on values in case of node merging.
1465  * If 2 nodes distant from less than 'eps' and with value different with more than 'epsOnVals' an exception will be thrown.
1466  */
1467 bool MEDCouplingFieldDouble::mergeNodes(double eps, double epsOnVals) throw(INTERP_KERNEL::Exception)
1468 {
1469   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
1470   if(!meshC)
1471     throw INTERP_KERNEL::Exception("Invalid support mesh to apply mergeNodes on it : must be a MEDCouplingPointSet one !");
1472   if(!((const MEDCouplingFieldDiscretization *)_type))
1473     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform mergeNodes !");
1474   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
1475   bool ret;
1476   int ret2;
1477   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->mergeNodes(eps,ret,ret2);
1478   if(!ret)//no nodes have been merged.
1479     return ret;
1480   std::vector<DataArrayDouble *> arrays;
1481   _time_discr->getArrays(arrays);
1482   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1483     if(*iter)
1484       _type->renumberValuesOnNodes(epsOnVals,arr->getConstPointer(),meshC2->getNumberOfNodes(),*iter);
1485   setMesh(meshC2);
1486   return true;
1487 }
1488
1489 /*!
1490  * Merge nodes with (barycenter computation) of underlying mesh. In case of some node will be merged the underlying mesh instance will change.
1491  * The first 'eps' stands for geometric approximation. The second 'epsOnVals' is for epsilon on values in case of node merging.
1492  * If 2 nodes distant from less than 'eps' and with value different with more than 'epsOnVals' an exception will be thrown.
1493  */
1494 bool MEDCouplingFieldDouble::mergeNodes2(double eps, double epsOnVals) throw(INTERP_KERNEL::Exception)
1495 {
1496   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
1497   if(!meshC)
1498     throw INTERP_KERNEL::Exception("Invalid support mesh to apply mergeNodes on it : must be a MEDCouplingPointSet one !");
1499   if(!((const MEDCouplingFieldDiscretization *)_type))
1500     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform mergeNodes2 !");
1501   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
1502   bool ret;
1503   int ret2;
1504   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->mergeNodes2(eps,ret,ret2);
1505   if(!ret)//no nodes have been merged.
1506     return ret;
1507   std::vector<DataArrayDouble *> arrays;
1508   _time_discr->getArrays(arrays);
1509   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1510     if(*iter)
1511       _type->renumberValuesOnNodes(epsOnVals,arr->getConstPointer(),meshC2->getNumberOfNodes(),*iter);
1512   setMesh(meshC2);
1513   return true;
1514 }
1515
1516 /*!
1517  * This method applyies ParaMEDMEM::MEDCouplingPointSet::zipCoords method on 'this->_mesh' that should be set and of type ParaMEDMEM::MEDCouplingPointSet.
1518  * If some nodes have disappeared true is returned.
1519  * 'epsOnVals' stands for epsilon in case of merge of cells. This value is used as tolerance in case the corresponding values differ.
1520  */
1521 bool MEDCouplingFieldDouble::zipCoords(double epsOnVals) throw(INTERP_KERNEL::Exception)
1522 {
1523   const MEDCouplingPointSet *meshC=dynamic_cast<const MEDCouplingPointSet *>(_mesh);
1524   if(!meshC)
1525     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::zipCoords : Invalid support mesh to apply zipCoords on it : must be a MEDCouplingPointSet one !");
1526   if(!((const MEDCouplingFieldDiscretization *)_type))
1527     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform zipCoords !");
1528   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> meshC2((MEDCouplingPointSet *)meshC->deepCpy());
1529   int oldNbOfNodes=meshC2->getNumberOfNodes();
1530   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->zipCoordsTraducer();
1531   if(meshC2->getNumberOfNodes()!=oldNbOfNodes)
1532     {
1533       std::vector<DataArrayDouble *> arrays;
1534       _time_discr->getArrays(arrays);
1535       for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1536         if(*iter)
1537           _type->renumberValuesOnNodes(epsOnVals,arr->getConstPointer(),meshC2->getNumberOfNodes(),*iter);
1538       setMesh(meshC2);
1539       return true;
1540     }
1541   return false;
1542 }
1543
1544 /*!
1545  * This method applyies ParaMEDMEM::MEDCouplingUMesh::zipConnectivityTraducer on 'this->_mesh' that should be set and of type ParaMEDMEM::MEDCouplingUMesh.
1546  * The semantic of 'compType' is given in ParaMEDMEM::MEDCouplingUMesh::zipConnectivityTraducer method.
1547  * 'epsOnVals' stands for epsilon in case of merge of cells. This value is used as tolerance in case the corresponding values differ.
1548  */
1549 bool MEDCouplingFieldDouble::zipConnectivity(int compType, double epsOnVals) throw(INTERP_KERNEL::Exception)
1550 {
1551   const MEDCouplingUMesh *meshC=dynamic_cast<const MEDCouplingUMesh *>(_mesh);
1552   if(!meshC)
1553     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::zipConnectivity : Invalid support mesh to apply zipCoords on it : must be a MEDCouplingPointSet one !");
1554   if(!((const MEDCouplingFieldDiscretization *)_type))
1555     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform zipConnectivity !");
1556   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> meshC2((MEDCouplingUMesh *)meshC->deepCpy());
1557   int oldNbOfCells=meshC2->getNumberOfCells();
1558   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->zipConnectivityTraducer(compType);
1559   if(meshC2->getNumberOfCells()!=oldNbOfCells)
1560     {
1561       std::vector<DataArrayDouble *> arrays;
1562       _time_discr->getArrays(arrays);
1563       for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1564         if(*iter)
1565           _type->renumberValuesOnCells(epsOnVals,meshC,arr->getConstPointer(),meshC2->getNumberOfCells(),*iter);
1566       setMesh(meshC2);
1567       return true;
1568     }
1569   return false;
1570 }
1571
1572 /*!
1573  * This method calls MEDCouplingUMesh::buildSlice3D method. So this method makes the assumption that underlying mesh exists.
1574  * For the moment, this method is implemented for fields on cells.
1575  * 
1576  * \return a newly allocated field double containing the result that the user should deallocate.
1577  */
1578 MEDCouplingFieldDouble *MEDCouplingFieldDouble::extractSlice3D(const double *origin, const double *vec, double eps) const throw(INTERP_KERNEL::Exception)
1579 {
1580   const MEDCouplingMesh *mesh=getMesh();
1581   if(!mesh)
1582     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::extractSlice3D : underlying mesh is null !");
1583   if(getTypeOfField()!=ON_CELLS)
1584     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::extractSlice3D : only implemented for fields on cells !");
1585   const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh(mesh->buildUnstructured());
1586   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=clone(false);
1587   ret->setMesh(umesh);
1588   DataArrayInt *cellIds=0;
1589   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mesh2=umesh->buildSlice3D(origin,vec,eps,cellIds);
1590   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIds2=cellIds;
1591   ret->setMesh(mesh2);
1592   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tupleIds=computeTupleIdsToSelectFromCellIds(cellIds->begin(),cellIds->end());
1593   std::vector<DataArrayDouble *> arrays;
1594   _time_discr->getArrays(arrays);
1595   int i=0;
1596   std::vector<DataArrayDouble *> newArr(arrays.size());
1597   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> > newArr2(arrays.size());
1598   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++,i++)
1599     {
1600       if(*iter)
1601         {
1602           newArr2[i]=(*iter)->selectByTupleIdSafe(cellIds->begin(),cellIds->end());
1603           newArr[i]=newArr2[i];
1604         }
1605     }
1606   ret->setArrays(newArr);
1607   return ret.retn();
1608 }
1609
1610 /*!
1611  * This method applyies ParaMEDMEM::MEDCouplingUMesh::simplexize on 'this->_mesh'.
1612  * The semantic of 'policy' is given in ParaMEDMEM::MEDCouplingUMesh::simplexize method.
1613  */
1614 bool MEDCouplingFieldDouble::simplexize(int policy) throw(INTERP_KERNEL::Exception)
1615 {
1616   if(!_mesh)
1617     throw INTERP_KERNEL::Exception("No underlying mesh on this field to perform simplexize !");
1618   if(!((const MEDCouplingFieldDiscretization *)_type))
1619     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform simplexize !");
1620   int oldNbOfCells=_mesh->getNumberOfCells();
1621   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> meshC2(_mesh->deepCpy());
1622   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arr=meshC2->simplexize(policy);
1623   int newNbOfCells=meshC2->getNumberOfCells();
1624   if(oldNbOfCells==newNbOfCells)
1625     return false;
1626   std::vector<DataArrayDouble *> arrays;
1627   _time_discr->getArrays(arrays);
1628   for(std::vector<DataArrayDouble *>::const_iterator iter=arrays.begin();iter!=arrays.end();iter++)
1629     if(*iter)
1630       _type->renumberValuesOnCellsR(_mesh,arr->getConstPointer(),arr->getNbOfElems(),*iter);
1631   setMesh(meshC2);
1632   return true;
1633 }
1634
1635 MEDCouplingFieldDouble *MEDCouplingFieldDouble::doublyContractedProduct() const throw(INTERP_KERNEL::Exception)
1636 {
1637   if(!((const MEDCouplingFieldDiscretization *)_type))
1638     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform doublyContractedProduct !");
1639   MEDCouplingTimeDiscretization *td=_time_discr->doublyContractedProduct();
1640   td->copyTinyAttrFrom(*_time_discr);
1641   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1642   ret->setName("DoublyContractedProduct");
1643   ret->setMesh(getMesh());
1644   return ret.retn();
1645 }
1646
1647 MEDCouplingFieldDouble *MEDCouplingFieldDouble::determinant() const throw(INTERP_KERNEL::Exception)
1648 {
1649   if(!((const MEDCouplingFieldDiscretization *)_type))
1650     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform determinant !");
1651   MEDCouplingTimeDiscretization *td=_time_discr->determinant();
1652   td->copyTinyAttrFrom(*_time_discr);
1653   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1654   ret->setName("Determinant");
1655   ret->setMesh(getMesh());
1656   return ret.retn();
1657 }
1658
1659 MEDCouplingFieldDouble *MEDCouplingFieldDouble::eigenValues() const throw(INTERP_KERNEL::Exception)
1660 {
1661   if(!((const MEDCouplingFieldDiscretization *)_type))
1662     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform eigenValues !");
1663   MEDCouplingTimeDiscretization *td=_time_discr->eigenValues();
1664   td->copyTinyAttrFrom(*_time_discr);
1665   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1666   ret->setName("EigenValues");
1667   ret->setMesh(getMesh());
1668   return ret.retn();
1669 }
1670
1671 MEDCouplingFieldDouble *MEDCouplingFieldDouble::eigenVectors() const throw(INTERP_KERNEL::Exception)
1672 {
1673   if(!((const MEDCouplingFieldDiscretization *)_type))
1674     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform eigenVectors !");
1675   MEDCouplingTimeDiscretization *td=_time_discr->eigenVectors();
1676   td->copyTinyAttrFrom(*_time_discr);
1677   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1678   ret->setName("EigenVectors");
1679   ret->setMesh(getMesh());
1680   return ret.retn();
1681 }
1682
1683 MEDCouplingFieldDouble *MEDCouplingFieldDouble::inverse() const throw(INTERP_KERNEL::Exception)
1684 {
1685   if(!((const MEDCouplingFieldDiscretization *)_type))
1686     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform inverse !");
1687   MEDCouplingTimeDiscretization *td=_time_discr->inverse();
1688   td->copyTinyAttrFrom(*_time_discr);
1689   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1690   ret->setName("Inversion");
1691   ret->setMesh(getMesh());
1692   return ret.retn();
1693 }
1694
1695 MEDCouplingFieldDouble *MEDCouplingFieldDouble::trace() const throw(INTERP_KERNEL::Exception)
1696 {
1697   if(!((const MEDCouplingFieldDiscretization *)_type))
1698     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform trace !");
1699   MEDCouplingTimeDiscretization *td=_time_discr->trace();
1700   td->copyTinyAttrFrom(*_time_discr);
1701   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1702   ret->setName("Trace");
1703   ret->setMesh(getMesh());
1704   return ret.retn();
1705 }
1706
1707 MEDCouplingFieldDouble *MEDCouplingFieldDouble::deviator() const throw(INTERP_KERNEL::Exception)
1708 {
1709   if(!((const MEDCouplingFieldDiscretization *)_type))
1710     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform deviator !");
1711   MEDCouplingTimeDiscretization *td=_time_discr->deviator();
1712   td->copyTinyAttrFrom(*_time_discr);
1713   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1714   ret->setName("Deviator");
1715   ret->setMesh(getMesh());
1716   return ret.retn();
1717 }
1718
1719 MEDCouplingFieldDouble *MEDCouplingFieldDouble::magnitude() const throw(INTERP_KERNEL::Exception)
1720 {
1721   if(!((const MEDCouplingFieldDiscretization *)_type))
1722     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform magnitude !");
1723   MEDCouplingTimeDiscretization *td=_time_discr->magnitude();
1724   td->copyTinyAttrFrom(*_time_discr);
1725   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1726   ret->setName("Magnitude");
1727   ret->setMesh(getMesh());
1728   return ret.retn();
1729 }
1730
1731 MEDCouplingFieldDouble *MEDCouplingFieldDouble::maxPerTuple() const throw(INTERP_KERNEL::Exception)
1732 {
1733   if(!((const MEDCouplingFieldDiscretization *)_type))
1734     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform maxPerTuple !");
1735   MEDCouplingTimeDiscretization *td=_time_discr->maxPerTuple();
1736   td->copyTinyAttrFrom(*_time_discr);
1737   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1738   std::ostringstream oss;
1739   oss << "Max_" << getName();
1740   ret->setName(oss.str().c_str());
1741   ret->setMesh(getMesh());
1742   return ret.retn();
1743 }
1744
1745 void MEDCouplingFieldDouble::changeNbOfComponents(int newNbOfComp, double dftValue) throw(INTERP_KERNEL::Exception)
1746 {
1747   _time_discr->changeNbOfComponents(newNbOfComp,dftValue);
1748 }
1749
1750 MEDCouplingFieldDouble *MEDCouplingFieldDouble::keepSelectedComponents(const std::vector<int>& compoIds) const throw(INTERP_KERNEL::Exception)
1751 {
1752   if(!((const MEDCouplingFieldDiscretization *)_type))
1753     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform keepSelectedComponents !");
1754   MEDCouplingTimeDiscretization *td=_time_discr->keepSelectedComponents(compoIds);
1755   td->copyTinyAttrFrom(*_time_discr);
1756   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1757   ret->setName(getName());
1758   ret->setMesh(getMesh());
1759   return ret.retn();
1760 }
1761
1762 void MEDCouplingFieldDouble::setSelectedComponents(const MEDCouplingFieldDouble *f, const std::vector<int>& compoIds) throw(INTERP_KERNEL::Exception)
1763 {
1764   _time_discr->setSelectedComponents(f->_time_discr,compoIds);
1765 }
1766
1767 void MEDCouplingFieldDouble::sortPerTuple(bool asc) throw(INTERP_KERNEL::Exception)
1768 {
1769   _time_discr->sortPerTuple(asc);
1770 }
1771
1772 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MergeFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1773 {
1774   if(!f1->areCompatibleForMerge(f2))
1775     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MergeFields on them !");
1776   const MEDCouplingMesh *m1=f1->getMesh();
1777   const MEDCouplingMesh *m2=f2->getMesh();
1778   if(!m1)
1779     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MergeFields : no underlying mesh of f1 !");
1780   if(!f1->_time_discr)
1781     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MergeFields : no time discr of f1 !");
1782   if(!f1->_type)
1783     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MergeFields : no spatial discr of f1 !");
1784   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> m=m1->mergeMyselfWith(m2);
1785   MEDCouplingTimeDiscretization *td=f1->_time_discr->aggregate(f2->_time_discr);
1786   td->copyTinyAttrFrom(*f1->_time_discr);
1787   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1788   ret->setMesh(m);
1789   ret->setName(f1->getName());
1790   ret->setDescription(f1->getDescription());
1791   return ret.retn();
1792 }
1793
1794 /*!
1795  * This method returns a newly created field that is the union of all fields in input array 'a'.
1796  * This method expects that 'a' is non empty. If not an exception will be thrown.
1797  * If there is only one field in 'a' a deepCopy (except time information of mesh and field) of the unique field instance in 'a' will be returned.
1798  * Generally speaking the first instance field in 'a' will be used to assign tiny attributes of returned field.
1799  */
1800 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MergeFields(const std::vector<const MEDCouplingFieldDouble *>& a) throw(INTERP_KERNEL::Exception)
1801 {
1802   if(a.size()<1)
1803     throw INTERP_KERNEL::Exception("FieldDouble::MergeFields : size of array must be >= 1 !");
1804   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > ms(a.size());
1805   std::vector< const MEDCouplingUMesh *> ms2(a.size());
1806   std::vector< const MEDCouplingTimeDiscretization *> tds(a.size());
1807   std::vector<const MEDCouplingFieldDouble *>::const_iterator it=a.begin();
1808   const MEDCouplingFieldDouble *ref=(*it++);
1809   if(!ref)
1810     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MergeFields : presence of NULL instance in first place of input vector !");
1811   for(;it!=a.end();it++)
1812     if(!ref->areCompatibleForMerge(*it))
1813       throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MergeFields on them !");
1814   for(int i=0;i<(int)a.size();i++)
1815     {
1816       if(!a[i]->getMesh())
1817         throw INTERP_KERNEL::Exception("MergeFields : A field as no underlying mesh !");
1818       ms[i]=a[i]->getMesh()->buildUnstructured();
1819       ms2[i]=ms[i];
1820       tds[i]=a[i]->_time_discr;
1821     }
1822   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=MEDCouplingUMesh::MergeUMeshes(ms2);
1823   m->setName(ms2[0]->getName()); m->setDescription(ms2[0]->getDescription());
1824   MEDCouplingTimeDiscretization *td=tds[0]->aggregate(tds);
1825   td->copyTinyAttrFrom(*(a[0]->_time_discr));
1826   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(a[0]->getNature(),td,a[0]->_type->clone());
1827   ret->setMesh(m);
1828   ret->setName(a[0]->getName());
1829   ret->setDescription(a[0]->getDescription());
1830   return ret.retn();
1831 }
1832
1833 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MeldFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1834 {
1835   if(!f1->areCompatibleForMeld(f2))
1836     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MeldFields on them !");
1837   MEDCouplingTimeDiscretization *td=f1->_time_discr->meld(f2->_time_discr);
1838   td->copyTinyAttrFrom(*f1->_time_discr);
1839   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1840   ret->setMesh(f1->getMesh());
1841   return ret.retn();
1842 }
1843
1844 MEDCouplingFieldDouble *MEDCouplingFieldDouble::DotFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1845 {
1846   if(!f1)
1847     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::DotFields : input field is NULL !");
1848   if(!f1->areStrictlyCompatible(f2))
1849     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply DotFields on them !");
1850   MEDCouplingTimeDiscretization *td=f1->_time_discr->dot(f2->_time_discr);
1851   td->copyTinyAttrFrom(*f1->_time_discr);
1852   MEDCouplingFieldDouble *ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1853   ret->setMesh(f1->getMesh());
1854   return ret;
1855 }
1856
1857 MEDCouplingFieldDouble *MEDCouplingFieldDouble::CrossProductFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1858 {
1859   if(!f1)
1860     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::CrossProductFields : input field is NULL !");
1861   if(!f1->areStrictlyCompatible(f2))
1862     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply CrossProductFields on them !");
1863   MEDCouplingTimeDiscretization *td=f1->_time_discr->crossProduct(f2->_time_discr);
1864   td->copyTinyAttrFrom(*f1->_time_discr);
1865   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1866   ret->setMesh(f1->getMesh());
1867   return ret.retn();
1868 }
1869
1870 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MaxFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1871 {
1872   if(!f1)
1873     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MaxFields : input field is NULL !");
1874   if(!f1->areStrictlyCompatible(f2))
1875     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MaxFields on them !");
1876   MEDCouplingTimeDiscretization *td=f1->_time_discr->max(f2->_time_discr);
1877   td->copyTinyAttrFrom(*f1->_time_discr);
1878   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1879   ret->setMesh(f1->getMesh());
1880   return ret.retn();
1881 }
1882
1883 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MinFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1884 {
1885   if(!f1)
1886     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MinFields : input field is NULL !");
1887   if(!f1->areStrictlyCompatible(f2))
1888     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MinFields on them !");
1889   MEDCouplingTimeDiscretization *td=f1->_time_discr->min(f2->_time_discr);
1890   td->copyTinyAttrFrom(*f1->_time_discr);
1891   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1892   ret->setMesh(f1->getMesh());
1893   return ret.retn();
1894 }
1895
1896 MEDCouplingFieldDouble *MEDCouplingFieldDouble::negate() const throw(INTERP_KERNEL::Exception)
1897 {
1898   if(!((const MEDCouplingFieldDiscretization *)_type))
1899     throw INTERP_KERNEL::Exception("No spatial discretization underlying this field to perform negate !");
1900   MEDCouplingTimeDiscretization *td=_time_discr->negate();
1901   td->copyTinyAttrFrom(*_time_discr);
1902   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(getNature(),td,_type->clone());
1903   ret->setMesh(getMesh());
1904   return ret.retn();
1905 }
1906
1907 MEDCouplingFieldDouble *MEDCouplingFieldDouble::AddFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1908 {
1909   if(!f1)
1910     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::AddFields : input field is NULL !");
1911   if(!f1->areStrictlyCompatible(f2))
1912     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply AddFields on them !");
1913   MEDCouplingTimeDiscretization *td=f1->_time_discr->add(f2->_time_discr);
1914   td->copyTinyAttrFrom(*f1->_time_discr);
1915   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1916   ret->setMesh(f1->getMesh());
1917   return ret.retn();
1918 }
1919
1920 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator+=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
1921 {
1922   if(!areStrictlyCompatible(&other))
1923     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply += on them !");
1924   _time_discr->addEqual(other._time_discr);
1925   return *this;
1926 }
1927
1928 MEDCouplingFieldDouble *MEDCouplingFieldDouble::SubstractFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1929 {
1930   if(!f1)
1931     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::SubstractFields : input field is NULL !");
1932   if(!f1->areStrictlyCompatible(f2))
1933     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply SubstractFields on them !");
1934   MEDCouplingTimeDiscretization *td=f1->_time_discr->substract(f2->_time_discr);
1935   td->copyTinyAttrFrom(*f1->_time_discr);
1936   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1937   ret->setMesh(f1->getMesh());
1938   return ret.retn();
1939 }
1940
1941 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator-=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
1942 {
1943   if(!areStrictlyCompatible(&other))
1944     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply -= on them !");
1945   _time_discr->substractEqual(other._time_discr);
1946   return *this;
1947 }
1948
1949 MEDCouplingFieldDouble *MEDCouplingFieldDouble::MultiplyFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1950 {
1951   if(!f1)
1952     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::MultiplyFields : input field is NULL !");
1953   if(!f1->areCompatibleForMul(f2))
1954     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply MultiplyFields on them !");
1955   MEDCouplingTimeDiscretization *td=f1->_time_discr->multiply(f2->_time_discr);
1956   td->copyTinyAttrFrom(*f1->_time_discr);
1957   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1958   ret->setMesh(f1->getMesh());
1959   return ret.retn();
1960 }
1961
1962 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator*=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
1963 {
1964   if(!areCompatibleForMul(&other))
1965     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply *= on them !");
1966   _time_discr->multiplyEqual(other._time_discr);
1967   return *this;
1968 }
1969
1970 MEDCouplingFieldDouble *MEDCouplingFieldDouble::DivideFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1971 {
1972   if(!f1)
1973     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::DivideFields : input field is NULL !");
1974   if(!f1->areCompatibleForDiv(f2))
1975     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply DivideFields on them !");
1976   MEDCouplingTimeDiscretization *td=f1->_time_discr->divide(f2->_time_discr);
1977   td->copyTinyAttrFrom(*f1->_time_discr);
1978   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
1979   ret->setMesh(f1->getMesh());
1980   return ret.retn();
1981 }
1982
1983 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator/=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
1984 {
1985   if(!areCompatibleForDiv(&other))
1986     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply /= on them !");
1987   _time_discr->divideEqual(other._time_discr);
1988   return *this;
1989 }
1990
1991 /*!
1992  * Directly called by MEDCouplingFieldDouble::operator^.
1993  * 
1994  * \sa MEDCouplingFieldDouble::operator^
1995  */
1996 MEDCouplingFieldDouble *MEDCouplingFieldDouble::PowFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception)
1997 {
1998   if(!f1)
1999     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::PowFields : input field is NULL !");
2000   if(!f1->areCompatibleForMul(f2))
2001     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply PowFields on them !");
2002   MEDCouplingTimeDiscretization *td=f1->_time_discr->pow(f2->_time_discr);
2003   td->copyTinyAttrFrom(*f1->_time_discr);
2004   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=new MEDCouplingFieldDouble(f1->getNature(),td,f1->_type->clone());
2005   ret->setMesh(f1->getMesh());
2006   return ret.retn();
2007 }
2008
2009 /*!
2010  * Directly call MEDCouplingFieldDouble::PowFields static method.
2011  * 
2012  * \sa MEDCouplingFieldDouble::PowFields
2013  */
2014 MEDCouplingFieldDouble *MEDCouplingFieldDouble::operator^(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception)
2015 {
2016   return PowFields(this,&other);
2017 }
2018
2019 const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator^=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception)
2020 {
2021   if(!areCompatibleForDiv(&other))
2022     throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply /= on them !");
2023   _time_discr->powEqual(other._time_discr);
2024   return *this;
2025 }
2026
2027 /*!
2028  * This method writes the field series 'fs' in the VTK file 'fileName'.
2029  * If 'fs' is empty no file is written. If fields lies on more than one mesh an exception will be thrown and no file will be written too.
2030  * If the single mesh is empty an exception will be thrown.
2031  * Finally there is a field in 'fs' with no name an exception will be thrown too.
2032  */
2033 void MEDCouplingFieldDouble::WriteVTK(const char *fileName, const std::vector<const MEDCouplingFieldDouble *>& fs) throw(INTERP_KERNEL::Exception)
2034 {
2035   if(fs.empty())
2036     return;
2037   std::size_t nfs=fs.size();
2038   if(!fs[0])
2039     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : 1st instance of field is NULL !");
2040   const MEDCouplingMesh *m=fs[0]->getMesh();
2041   if(!m)
2042     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : 1st instance of field lies on NULL mesh !");
2043   for(std::size_t i=1;i<nfs;i++)
2044     if(fs[i]->getMesh()!=m)
2045       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.");
2046   if(!m)
2047     throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::WriteVTK : Fields are lying on a same mesh but it is empty !");
2048   std::ostringstream coss,noss;
2049   for(std::size_t i=0;i<nfs;i++)
2050     {
2051       const MEDCouplingFieldDouble *cur=fs[i];
2052       std::string name(cur->getName());
2053       if(name.empty())
2054         {
2055           std::ostringstream oss; oss << "MEDCouplingFieldDouble::WriteVTK : Field in pos #" << i << " has no name !";
2056           throw INTERP_KERNEL::Exception(oss.str().c_str());
2057         }
2058       TypeOfField typ=cur->getTypeOfField();
2059       if(typ==ON_CELLS)
2060         cur->getArray()->writeVTK(coss,8,cur->getName());
2061       else if(typ==ON_NODES)
2062         cur->getArray()->writeVTK(noss,8,cur->getName());
2063     }
2064   m->writeVTKAdvanced(fileName,coss.str(),noss.str());
2065 }
2066
2067 void MEDCouplingFieldDouble::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
2068 {
2069   stream << "MEDCouplingFieldDouble C++ instance at " << this << ". Name : \"" << _name << "\"." << std::endl;
2070   const char *nat=0;
2071   try
2072     {
2073       nat=MEDCouplingNatureOfField::GetRepr(_nature);
2074       stream << "Nature of field : " << nat << ".\n";
2075     }
2076   catch(INTERP_KERNEL::Exception& e)
2077     {  }
2078   const MEDCouplingFieldDiscretization *fd(_type);
2079   if(!fd)
2080     stream << "No spatial discretization set !";
2081   else
2082     fd->reprQuickOverview(stream);
2083   stream << std::endl;
2084   if(!_mesh)
2085     stream << "\nNo mesh support defined !";
2086   else
2087     {
2088       std::ostringstream oss;
2089       _mesh->reprQuickOverview(oss);
2090       std::string tmp(oss.str());
2091       stream << "\nMesh info : " << tmp.substr(0,tmp.find('\n'));
2092     }
2093   if(_time_discr)
2094     {
2095       const DataArrayDouble *arr=_time_discr->getArray();
2096       if(arr)
2097         {
2098           stream << "\n\nArray info : ";
2099           arr->reprQuickOverview(stream);
2100         }
2101       else
2102         {
2103           stream << "\n\nNo data array set !";
2104         }
2105     }
2106 }