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