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