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