src/MEDCoupling/MEDCouplingFieldDouble.cxx

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