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