src/MEDCoupling/MEDCouplingMemArray.cxx

   1 // Copyright (C) 2007-2020  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, or (at your option) any later version.
   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 (EDF R&D)
  20
  21 #include "MEDCouplingMemArray.txx"
  22
  23 #include "BBTree.txx"
  24 #include "GenMathFormulae.hxx"
  25 #include "InterpKernelAutoPtr.hxx"
  26 #include "InterpKernelExprParser.hxx"
  27
  28 #include "InterpKernelAutoPtr.hxx"
  29 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  30 #include "InterpKernelAutoPtr.hxx"
  31 #include "InterpKernelGeo2DNode.hxx"
  32 #include "InterpKernelGeo2DEdgeLin.hxx"
  33
  34 #include <set>
  35 #include <cmath>
  36 #include <limits>
  37 #include <numeric>
  38 #include <algorithm>
  39 #include <functional>
  40
  41 typedef double (*MYFUNCPTR)(double);
  42
  43 using namespace MEDCoupling;
  44
  45 template class MEDCOUPLING_EXPORT MEDCoupling::MemArray<mcIdType>;
  46 template class MEDCOUPLING_EXPORT MEDCoupling::MemArray<double>;
  47 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplate<mcIdType>;
  48 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplate<double>;
  49 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateClassic<Int32>;
  50 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateClassic<Int64>;
  51 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateClassic<double>;
  52 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateFP<double>;
  53 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayIterator<double>;
  54 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayIterator<mcIdType>;
  55 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscrete<Int32>;
  56 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscreteSigned<Int32>;
  57 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscrete<Int64>;
  58 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscreteSigned<Int64>;
  59 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTuple<mcIdType>;
  60 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTuple<double>;
  61 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTuple<float>;
  62
  63 void MEDCoupling::DACheckNbOfTuplesAndComp(const DataArray *da, mcIdType nbOfTuples, std::size_t nbOfCompo, const std::string& msg)
  64 {
  65   if(!da)
  66     throw INTERP_KERNEL::Exception("DACheckNbOfTuplesAndComp : null input object !");
  67   da->checkNbOfTuplesAndComp(nbOfTuples,nbOfCompo,msg);
  68 }
  69
  70 template<mcIdType SPACEDIM>
  71 void DataArrayDouble::findCommonTuplesAlg(const double *bbox, mcIdType nbNodes, mcIdType limitNodeId, double prec, DataArrayIdType *c, DataArrayIdType *cI) const
  72 {
  73   const double *coordsPtr=getConstPointer();
  74   BBTreePts<SPACEDIM,mcIdType> myTree(bbox,0,0,nbNodes,prec);
  75   std::vector<bool> isDone(nbNodes);
  76   for(mcIdType i=0;i<nbNodes;i++)
  77     {
  78       if(!isDone[i])
  79         {
  80           std::vector<mcIdType> intersectingElems;
  81           myTree.getElementsAroundPoint(coordsPtr+i*SPACEDIM,intersectingElems);
  82           if(intersectingElems.size()>1)
  83             {
  84               std::vector<mcIdType> commonNodes;
  85               for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
  86                 if(*it!=i)
  87                   if(*it>=limitNodeId)
  88                     {
  89                       commonNodes.push_back(*it);
  90                       isDone[*it]=true;
  91                     }
  92               if(!commonNodes.empty())
  93                 {
  94                   cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size())+1);
  95                   c->pushBackSilent(i);
  96                   c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
  97                 }
  98             }
  99         }
 100     }
 101 }
 102
 103 template<mcIdType SPACEDIM>
 104 void DataArrayDouble::FindTupleIdsNearTuplesAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, const double *pos, mcIdType nbOfTuples, double eps,
 105                                                 DataArrayIdType *c, DataArrayIdType *cI)
 106 {
 107   for(mcIdType i=0;i<nbOfTuples;i++)
 108     {
 109       std::vector<mcIdType> intersectingElems;
 110       myTree.getElementsAroundPoint(pos+i*SPACEDIM,intersectingElems);
 111       std::vector<mcIdType> commonNodes;
 112       for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
 113         commonNodes.push_back(*it);
 114       cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size()));
 115       c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
 116     }
 117 }
 118
 119 template<mcIdType SPACEDIM>
 120 void DataArrayDouble::FindClosestTupleIdAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, double dist, const double *pos, mcIdType nbOfTuples, const double *thisPt, mcIdType thisNbOfTuples, mcIdType *res)
 121 {
 122   double distOpt = std::max(dist, std::numeric_limits<double>::epsilon());
 123   const double *p(pos);
 124   mcIdType *r(res);
 125   for(mcIdType i=0;i<nbOfTuples;i++,p+=SPACEDIM,r++)
 126     {
 127       while(true)
 128         {
 129           mcIdType elem=-1;
 130           double ret=myTree.getElementsAroundPoint2(p,distOpt,elem);
 131           if(ret!=std::numeric_limits<double>::max())
 132             {
 133               distOpt=std::max(ret,1e-4);
 134               *r=elem;
 135               break;
 136             }
 137           else
 138             { distOpt=2*distOpt; continue; }
 139         }
 140     }
 141 }
 142
 143 mcIdType DataArray::EffectiveCircPerm(mcIdType nbOfShift, mcIdType nbOfTuples)
 144 {
 145   if(nbOfTuples<=0)
 146     throw INTERP_KERNEL::Exception("DataArray::EffectiveCircPerm : number of tuples is expected to be > 0 !");
 147   if(nbOfShift>=0)
 148     {
 149       return nbOfShift%nbOfTuples;
 150     }
 151   else
 152     {
 153       mcIdType tmp(-nbOfShift);
 154       tmp=tmp%nbOfTuples;
 155       return nbOfTuples-tmp;
 156     }
 157 }
 158
 159 std::size_t DataArray::getHeapMemorySizeWithoutChildren() const
 160 {
 161   std::size_t sz1=_name.capacity();
 162   std::size_t sz2=_info_on_compo.capacity();
 163   std::size_t sz3=0;
 164   for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
 165     sz3+=(*it).capacity();
 166   return sz1+sz2+sz3;
 167 }
 168
 169 std::vector<const BigMemoryObject *> DataArray::getDirectChildrenWithNull() const
 170 {
 171   return std::vector<const BigMemoryObject *>();
 172 }
 173
 174 /*!
 175  * Sets the attribute \a _name of \a this array.
 176  * See \ref MEDCouplingArrayBasicsName "DataArrays infos" for more information.
 177  *  \param [in] name - new array name
 178  */
 179 void DataArray::setName(const std::string& name)
 180 {
 181   _name=name;
 182 }
 183
 184 /*!
 185  * Copies textual data from an \a other DataArray. The copied data are
 186  * - the name attribute,
 187  * - the information of components.
 188  *
 189  * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
 190  *
 191  *  \param [in] other - another instance of DataArray to copy the textual data from.
 192  *  \throw If number of components of \a this array differs from that of the \a other.
 193  */
 194 void DataArray::copyStringInfoFrom(const DataArray& other)
 195 {
 196   if(_info_on_compo.size()!=other._info_on_compo.size())
 197     throw INTERP_KERNEL::Exception("Size of arrays mismatches on copyStringInfoFrom !");
 198   _name=other._name;
 199   _info_on_compo=other._info_on_compo;
 200 }
 201
 202 void DataArray::copyPartOfStringInfoFrom(const DataArray& other, const std::vector<std::size_t>& compoIds)
 203 {
 204   std::size_t nbOfCompoOth=other.getNumberOfComponents();
 205   std::size_t newNbOfCompo=compoIds.size();
 206   for(std::size_t i=0;i<newNbOfCompo;i++)
 207     if(compoIds[i]>=nbOfCompoOth)
 208       {
 209         std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompoOth << ")";
 210         throw INTERP_KERNEL::Exception(oss.str().c_str());
 211       }
 212   for(std::size_t i=0;i<newNbOfCompo;i++)
 213     setInfoOnComponent(i,other.getInfoOnComponent(compoIds[i]));
 214 }
 215
 216 void DataArray::copyPartOfStringInfoFrom2(const std::vector<std::size_t>& compoIds, const DataArray& other)
 217 {
 218   if(compoIds.size()!=other.getNumberOfComponents())
 219     throw INTERP_KERNEL::Exception("Given compoIds has not the same size as number of components of given array !");
 220   std::size_t partOfCompoToSet=compoIds.size();
 221   std::size_t nbOfCompo=getNumberOfComponents();
 222   for(std::size_t i=0;i<partOfCompoToSet;i++)
 223     if(compoIds[i]>=nbOfCompo)
 224       {
 225         std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompo << ")";
 226         throw INTERP_KERNEL::Exception(oss.str().c_str());
 227       }
 228   for(std::size_t i=0;i<partOfCompoToSet;i++)
 229     setInfoOnComponent(compoIds[i],other.getInfoOnComponent(i));
 230 }
 231
 232 bool DataArray::areInfoEqualsIfNotWhy(const DataArray& other, std::string& reason) const
 233 {
 234   std::ostringstream oss;
 235   if(_name!=other._name)
 236     {
 237       oss << "Names DataArray mismatch : this name=\"" << _name << " other name=\"" << other._name << "\" !";
 238       reason=oss.str();
 239       return false;
 240     }
 241   if(_info_on_compo!=other._info_on_compo)
 242     {
 243       oss << "Components DataArray mismatch : \nThis components=";
 244       for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
 245         oss << "\"" << *it << "\",";
 246       oss << "\nOther components=";
 247       for(std::vector<std::string>::const_iterator it=other._info_on_compo.begin();it!=other._info_on_compo.end();it++)
 248         oss << "\"" << *it << "\",";
 249       reason=oss.str();
 250       return false;
 251     }
 252   return true;
 253 }
 254
 255 /*!
 256  * Compares textual information of \a this DataArray with that of an \a other one.
 257  * The compared data are
 258  * - the name attribute,
 259  * - the information of components.
 260  *
 261  * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
 262  *  \param [in] other - another instance of DataArray to compare the textual data of.
 263  *  \return bool - \a true if the textual information is same, \a false else.
 264  */
 265 bool DataArray::areInfoEquals(const DataArray& other) const
 266 {
 267   std::string tmp;
 268   return areInfoEqualsIfNotWhy(other,tmp);
 269 }
 270
 271 void DataArray::reprWithoutNameStream(std::ostream& stream) const
 272 {
 273   stream << "Number of components : "<< getNumberOfComponents() << "\n";
 274   stream << "Info of these components : ";
 275   for(std::vector<std::string>::const_iterator iter=_info_on_compo.begin();iter!=_info_on_compo.end();iter++)
 276     stream << "\"" << *iter << "\"   ";
 277   stream << "\n";
 278 }
 279
 280 std::string DataArray::cppRepr(const std::string& varName) const
 281 {
 282   std::ostringstream ret;
 283   reprCppStream(varName,ret);
 284   return ret.str();
 285 }
 286
 287 /*!
 288  * Sets information on all components. To know more on format of this information
 289  * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 290  *  \param [in] info - a vector of strings.
 291  *  \throw If size of \a info differs from the number of components of \a this.
 292  */
 293 void DataArray::setInfoOnComponents(const std::vector<std::string>& info)
 294 {
 295   if(getNumberOfComponents()!=info.size())
 296     {
 297       std::ostringstream oss; oss << "DataArray::setInfoOnComponents : input is of size " << info.size() << " whereas number of components is equal to " << getNumberOfComponents() << " !";
 298       throw INTERP_KERNEL::Exception(oss.str().c_str());
 299     }
 300   _info_on_compo=info;
 301 }
 302
 303 /*!
 304  * This method is only a dispatcher towards DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3 depending on the true
 305  * type of \a this and \a aBase.
 306  *
 307  * \throw If \a aBase and \a this do not have the same type.
 308  *
 309  * \sa DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3.
 310  */
 311 void DataArray::setPartOfValuesBase3(const DataArray *aBase, const mcIdType *bgTuples, const mcIdType *endTuples, mcIdType bgComp, mcIdType endComp, mcIdType stepComp, bool strictCompoCompare)
 312 {
 313   if(!aBase)
 314     throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object is NULL !");
 315   DataArrayDouble *this1(dynamic_cast<DataArrayDouble *>(this));
 316   DataArrayIdType *this2(dynamic_cast<DataArrayIdType *>(this));
 317   DataArrayChar *this3(dynamic_cast<DataArrayChar *>(this));
 318   const DataArrayDouble *a1(dynamic_cast<const DataArrayDouble *>(aBase));
 319   const DataArrayIdType *a2(dynamic_cast<const DataArrayIdType *>(aBase));
 320   const DataArrayChar *a3(dynamic_cast<const DataArrayChar *>(aBase));
 321   if(this1 && a1)
 322     {
 323       this1->setPartOfValues3(a1,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
 324       return ;
 325     }
 326   if(this2 && a2)
 327     {
 328       this2->setPartOfValues3(a2,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
 329       return ;
 330     }
 331   if(this3 && a3)
 332     {
 333       this3->setPartOfValues3(a3,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
 334       return ;
 335     }
 336   throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object and this do not have the same type !");
 337 }
 338
 339 std::vector<std::string> DataArray::getVarsOnComponent() const
 340 {
 341   std::size_t nbOfCompo=_info_on_compo.size();
 342   std::vector<std::string> ret(nbOfCompo);
 343   for(std::size_t i=0;i<nbOfCompo;i++)
 344     ret[i]=getVarOnComponent(i);
 345   return ret;
 346 }
 347
 348 std::vector<std::string> DataArray::getUnitsOnComponent() const
 349 {
 350   std::size_t nbOfCompo=_info_on_compo.size();
 351   std::vector<std::string> ret(nbOfCompo);
 352   for(std::size_t i=0;i<nbOfCompo;i++)
 353     ret[i]=getUnitOnComponent(i);
 354   return ret;
 355 }
 356
 357 /*!
 358  * Returns information on a component specified by an index.
 359  * To know more on format of this information
 360  * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 361  *  \param [in] i - the index (zero based) of the component of interest.
 362  *  \return std::string - a string containing the information on \a i-th component.
 363  *  \throw If \a i is not a valid component index.
 364  */
 365 std::string DataArray::getInfoOnComponent(std::size_t i) const
 366 {
 367   if(i<_info_on_compo.size())
 368     return _info_on_compo[i];
 369   else
 370     {
 371       std::ostringstream oss; oss << "DataArray::getInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
 372       throw INTERP_KERNEL::Exception(oss.str().c_str());
 373     }
 374 }
 375
 376 /*!
 377  * Returns the var part of the full information of the \a i-th component.
 378  * For example, if \c getInfoOnComponent(0) returns "SIGXY [N/m^2]", then
 379  * \c getVarOnComponent(0) returns "SIGXY".
 380  * If a unit part of information is not detected by presence of
 381  * two square brackets, then the full information is returned.
 382  * To read more about the component information format, see
 383  * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 384  *  \param [in] i - the index (zero based) of the component of interest.
 385  *  \return std::string - a string containing the var information, or the full info.
 386  *  \throw If \a i is not a valid component index.
 387  */
 388 std::string DataArray::getVarOnComponent(std::size_t i) const
 389 {
 390   if(i<_info_on_compo.size())
 391     {
 392       return GetVarNameFromInfo(_info_on_compo[i]);
 393     }
 394   else
 395     {
 396       std::ostringstream oss; oss << "DataArray::getVarOnComponent : Specified component id is out of range  (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
 397       throw INTERP_KERNEL::Exception(oss.str().c_str());
 398     }
 399 }
 400
 401 /*!
 402  * Returns the unit part of the full information of the \a i-th component.
 403  * For example, if \c getInfoOnComponent(0) returns "SIGXY [ N/m^2]", then
 404  * \c getUnitOnComponent(0) returns " N/m^2".
 405  * If a unit part of information is not detected by presence of
 406  * two square brackets, then an empty string is returned.
 407  * To read more about the component information format, see
 408  * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 409  *  \param [in] i - the index (zero based) of the component of interest.
 410  *  \return std::string - a string containing the unit information, if any, or "".
 411  *  \throw If \a i is not a valid component index.
 412  */
 413 std::string DataArray::getUnitOnComponent(std::size_t i) const
 414 {
 415   if(i<_info_on_compo.size())
 416     {
 417       return GetUnitFromInfo(_info_on_compo[i]);
 418     }
 419   else
 420     {
 421       std::ostringstream oss; oss << "DataArray::getUnitOnComponent : Specified component id is out of range  (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
 422       throw INTERP_KERNEL::Exception(oss.str().c_str());
 423     }
 424 }
 425
 426 /*!
 427  * Returns the var part of the full component information.
 428  * For example, if \a info == "SIGXY [N/m^2]", then this method returns "SIGXY".
 429  * If a unit part of information is not detected by presence of
 430  * two square brackets, then the whole \a info is returned.
 431  * To read more about the component information format, see
 432  * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 433  *  \param [in] info - the full component information.
 434  *  \return std::string - a string containing only var information, or the \a info.
 435  */
 436 std::string DataArray::GetVarNameFromInfo(const std::string& info)
 437 {
 438   std::size_t p1=info.find_last_of('[');
 439   std::size_t p2=info.find_last_of(']');
 440   if(p1==std::string::npos || p2==std::string::npos)
 441     return info;
 442   if(p1>p2)
 443     return info;
 444   if(p1==0)
 445     return std::string();
 446   std::size_t p3=info.find_last_not_of(' ',p1-1);
 447   return info.substr(0,p3+1);
 448 }
 449
 450 /*!
 451  * Returns the unit part of the full component information.
 452  * For example, if \a info == "SIGXY [ N/m^2]", then this method returns " N/m^2".
 453  * If a unit part of information is not detected by presence of
 454  * two square brackets, then an empty string is returned.
 455  * To read more about the component information format, see
 456  * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 457  *  \param [in] info - the full component information.
 458  *  \return std::string - a string containing only unit information, if any, or "".
 459  */
 460 std::string DataArray::GetUnitFromInfo(const std::string& info)
 461 {
 462   std::size_t p1=info.find_last_of('[');
 463   std::size_t p2=info.find_last_of(']');
 464   if(p1==std::string::npos || p2==std::string::npos)
 465     return std::string();
 466   if(p1>p2)
 467     return std::string();
 468   return info.substr(p1+1,p2-p1-1);
 469 }
 470
 471 /*!
 472  * This method put in info format the result of the merge of \a var and \a unit.
 473  * The standard format for that is "var [unit]".
 474  * Inversely you can retrieve the var part or the unit part of info string using resp. GetVarNameFromInfo and GetUnitFromInfo.
 475  */
 476 std::string DataArray::BuildInfoFromVarAndUnit(const std::string& var, const std::string& unit)
 477 {
 478   std::ostringstream oss;
 479   oss << var << " [" << unit << "]";
 480   return oss.str();
 481 }
 482
 483 std::string DataArray::GetAxisTypeRepr(MEDCouplingAxisType at)
 484 {
 485   switch(at)
 486     {
 487     case AX_CART:
 488       return std::string("AX_CART");
 489     case AX_CYL:
 490       return std::string("AX_CYL");
 491     case AX_SPHER:
 492       return std::string("AX_SPHER");
 493     default:
 494       throw INTERP_KERNEL::Exception("DataArray::GetAxisTypeRepr : unrecognized axis type enum !");
 495     }
 496 }
 497
 498 /*!
 499  * Returns a new DataArray by concatenating all given arrays, so that (1) the number
 500  * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
 501  * the number of component in the result array is same as that of each of given arrays.
 502  * Info on components is copied from the first of the given arrays. Number of components
 503  * in the given arrays must be  the same.
 504  *  \param [in] arrs - a sequence of arrays to include in the result array. All arrays must have the same type.
 505  *  \return DataArray * - the new instance of DataArray (that can be either DataArrayInt, DataArrayDouble, DataArrayChar).
 506  *          The caller is to delete this result array using decrRef() as it is no more
 507  *          needed.
 508  *  \throw If all arrays within \a arrs are NULL.
 509  *  \throw If all not null arrays in \a arrs have not the same type.
 510  *  \throw If getNumberOfComponents() of arrays within \a arrs.
 511  */
 512 DataArray *DataArray::Aggregate(const std::vector<const DataArray *>& arrs)
 513 {
 514   std::vector<const DataArray *> arr2;
 515   for(std::vector<const DataArray *>::const_iterator it=arrs.begin();it!=arrs.end();it++)
 516     if(*it)
 517       arr2.push_back(*it);
 518   if(arr2.empty())
 519     throw INTERP_KERNEL::Exception("DataArray::Aggregate : only null instance in input vector !");
 520   std::vector<const DataArrayDouble *> arrd;
 521   std::vector<const DataArrayIdType *> arri;
 522   std::vector<const DataArrayChar *> arrc;
 523   for(std::vector<const DataArray *>::const_iterator it=arr2.begin();it!=arr2.end();it++)
 524     {
 525       const DataArrayDouble *a=dynamic_cast<const DataArrayDouble *>(*it);
 526       if(a)
 527         { arrd.push_back(a); continue; }
 528       const DataArrayIdType *b=dynamic_cast<const DataArrayIdType *>(*it);
 529       if(b)
 530         { arri.push_back(b); continue; }
 531       const DataArrayChar *c=dynamic_cast<const DataArrayChar *>(*it);
 532       if(c)
 533         { arrc.push_back(c); continue; }
 534       throw INTERP_KERNEL::Exception("DataArray::Aggregate : presence of not null instance in inuput that is not in [DataArrayDouble, DataArrayInt, DataArrayChar] !");
 535     }
 536   if(arr2.size()==arrd.size())
 537     return DataArrayDouble::Aggregate(arrd);
 538   if(arr2.size()==arri.size())
 539     return DataArrayIdType::Aggregate(arri);
 540   if(arr2.size()==arrc.size())
 541     return DataArrayChar::Aggregate(arrc);
 542   throw INTERP_KERNEL::Exception("DataArray::Aggregate : all input arrays must have the same type !");
 543 }
 544
 545 /*!
 546  * Sets information on a component specified by an index.
 547  * To know more on format of this information
 548  * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 549  *  \warning Don't pass NULL as \a info!
 550  *  \param [in] i - the index (zero based) of the component of interest.
 551  *  \param [in] info - the string containing the information.
 552  *  \throw If \a i is not a valid component index.
 553  */
 554 void DataArray::setInfoOnComponent(std::size_t i, const std::string& info)
 555 {
 556   if(i<_info_on_compo.size())
 557     _info_on_compo[i]=info;
 558   else
 559     {
 560       std::ostringstream oss; oss << "DataArray::setInfoOnComponent : Specified component id is out of range  (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
 561       throw INTERP_KERNEL::Exception(oss.str().c_str());
 562     }
 563 }
 564
 565 /*!
 566  * Sets information on all components. This method can change number of components
 567  * at certain conditions; if the conditions are not respected, an exception is thrown.
 568  * The number of components can be changed in \a this only if \a this is not allocated.
 569  * The condition of number of components must not be changed.
 570  *
 571  * To know more on format of the component information see
 572  * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
 573  *  \param [in] info - a vector of component infos.
 574  *  \throw If \a this->getNumberOfComponents() != \a info.size() && \a this->isAllocated()
 575  */
 576 void DataArray::setInfoAndChangeNbOfCompo(const std::vector<std::string>& info)
 577 {
 578   if(getNumberOfComponents()!=info.size())
 579     {
 580       if(!isAllocated())
 581         _info_on_compo=info;
 582       else
 583         {
 584           std::ostringstream oss; oss << "DataArray::setInfoAndChangeNbOfCompo : input is of size " << info.size() << " whereas number of components is equal to " << getNumberOfComponents() << "  and this is already allocated !";
 585           throw INTERP_KERNEL::Exception(oss.str().c_str());
 586         }
 587     }
 588   else
 589     _info_on_compo=info;
 590 }
 591
 592 void DataArray::checkNbOfTuples(mcIdType nbOfTuples, const std::string& msg) const
 593 {
 594   if(getNumberOfTuples()!=nbOfTuples)
 595     {
 596       std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " <<  nbOfTuples << " having " << getNumberOfTuples() << " !";
 597       throw INTERP_KERNEL::Exception(oss.str().c_str());
 598     }
 599 }
 600
 601 void DataArray::checkNbOfComps(std::size_t nbOfCompo, const std::string& msg) const
 602 {
 603   if (getNumberOfComponents()!=nbOfCompo)
 604     {
 605       std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << nbOfCompo << " having " << getNumberOfComponents() << " !";
 606       throw INTERP_KERNEL::Exception(oss.str().c_str());
 607     }
 608 }
 609
 610 void DataArray::checkNbOfElems(mcIdType nbOfElems, const std::string& msg) const
 611 {
 612   if(getNbOfElems()!=nbOfElems)
 613     {
 614       std::ostringstream oss; oss << msg << " : mismatch number of elems : Expected " << nbOfElems << " having " << getNbOfElems() << " !";
 615       throw INTERP_KERNEL::Exception(oss.str().c_str());
 616     }
 617 }
 618
 619 void DataArray::checkNbOfTuplesAndComp(const DataArray& other, const std::string& msg) const
 620 {
 621   if(getNumberOfTuples()!=other.getNumberOfTuples())
 622     {
 623       std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " <<  other.getNumberOfTuples() << " having " << getNumberOfTuples() << " !";
 624       throw INTERP_KERNEL::Exception(oss.str().c_str());
 625     }
 626   if(getNumberOfComponents()!=other.getNumberOfComponents())
 627     {
 628       std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << other.getNumberOfComponents() << " having " << getNumberOfComponents() << " !";
 629       throw INTERP_KERNEL::Exception(oss.str().c_str());
 630     }
 631 }
 632
 633 void DataArray::checkNbOfTuplesAndComp(mcIdType nbOfTuples, std::size_t nbOfCompo, const std::string& msg) const
 634 {
 635   checkNbOfTuples(nbOfTuples,msg);
 636   checkNbOfComps(nbOfCompo,msg);
 637 }
 638
 639 /*!
 640  * Simply this method checks that \b value is in [0,\b ref).
 641  */
 642 void DataArray::CheckValueInRange(mcIdType ref, mcIdType value, const std::string& msg)
 643 {
 644   if(value<0 || value>=ref)
 645     {
 646       std::ostringstream oss; oss << "DataArray::CheckValueInRange : " << msg  << " ! Expected in range [0," << ref << "[ having " << value << " !";
 647       throw INTERP_KERNEL::Exception(oss.str().c_str());
 648     }
 649 }
 650
 651 /*!
 652  * This method checks that [\b start, \b end) is compliant with ref length \b value.
 653  * typically start in [0,\b value) and end in [0,\b value). If value==start and start==end, it is supported.
 654  */
 655 void DataArray::CheckValueInRangeEx(mcIdType value, mcIdType start, mcIdType end, const std::string& msg)
 656 {
 657   if(start<0 || start>=value)
 658     {
 659       if(value!=start || end!=start)
 660         {
 661           std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg  << " ! Expected start " << start << " of input range, in [0," << value << "[ !";
 662           throw INTERP_KERNEL::Exception(oss.str().c_str());
 663         }
 664     }
 665   if(end<0 || end>value)
 666     {
 667       std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg  << " ! Expected end " << end << " of input range, in [0," << value << "] !";
 668       throw INTERP_KERNEL::Exception(oss.str().c_str());
 669     }
 670 }
 671
 672 void DataArray::CheckClosingParInRange(mcIdType ref, mcIdType value, const std::string& msg)
 673 {
 674   if(value<0 || value>ref)
 675     {
 676       std::ostringstream oss; oss << "DataArray::CheckClosingParInRange : " << msg  << " ! Expected input range in [0," << ref << "] having closing open parenthesis " << value << " !";
 677       throw INTERP_KERNEL::Exception(oss.str().c_str());
 678     }
 679 }
 680
 681 /*!
 682  * This method is useful to slice work among a pool of threads or processes. \a begin, \a end \a step is the input whole slice of work to perform,
 683  * typically it is a whole slice of tuples of DataArray or cells, nodes of a mesh...
 684  *
 685  * The input \a sliceId should be an id in [0, \a nbOfSlices) that specifies the slice of work.
 686  *
 687  * \param [in] start - the start of the input slice of the whole work to perform split into slices.
 688  * \param [in] stop - the stop of the input slice of the whole work to perform split into slices.
 689  * \param [in] step - the step (that can be <0) of the input slice of the whole work to perform split into slices.
 690  * \param [in] sliceId - the slice id considered
 691  * \param [in] nbOfSlices - the number of slices (typically the number of cores on which the work is expected to be sliced)
 692  * \param [out] startSlice - the start of the slice considered
 693  * \param [out] stopSlice - the stop of the slice consided
 694  *
 695  * \throw If \a step == 0
 696  * \throw If \a nbOfSlices not > 0
 697  * \throw If \a sliceId not in [0,nbOfSlices)
 698  */
 699 void DataArray::GetSlice(mcIdType start, mcIdType stop, mcIdType step, mcIdType sliceId, mcIdType nbOfSlices, mcIdType& startSlice, mcIdType& stopSlice)
 700 {
 701   DataArrayTools<mcIdType>::GetSlice(start, stop, step, sliceId, nbOfSlices, startSlice, stopSlice);
 702 }
 703
 704 mcIdType DataArray::GetNumberOfItemGivenBES(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
 705 {
 706   return DataArrayTools<mcIdType>::GetNumberOfItemGivenBES(begin, end, step, msg);
 707 }
 708
 709 mcIdType DataArray::GetNumberOfItemGivenBESRelative(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
 710 {
 711   return DataArrayTools<mcIdType>::GetNumberOfItemGivenBESRelative(begin, end, step, msg);
 712 }
 713
 714 mcIdType DataArray::GetPosOfItemGivenBESRelativeNoThrow(mcIdType value, mcIdType begin, mcIdType end, mcIdType step)
 715 {
 716   return DataArrayTools<mcIdType>::GetPosOfItemGivenBESRelativeNoThrow(value, begin, end, step);
 717 }
 718
 719 /*!
 720  * Returns a new instance of DataArrayDouble. The caller is to delete this array
 721  * using decrRef() as it is no more needed.
 722  */
 723 DataArrayDouble *DataArrayDouble::New()
 724 {
 725   return new DataArrayDouble;
 726 }
 727
 728 /*!
 729  * Returns the only one value in \a this, if and only if number of elements
 730  * (nb of tuples * nb of components) is equal to 1, and that \a this is allocated.
 731  *  \return double - the sole value stored in \a this array.
 732  *  \throw If at least one of conditions stated above is not fulfilled.
 733  */
 734 double DataArrayDouble::doubleValue() const
 735 {
 736   if(isAllocated())
 737     {
 738       if(getNbOfElems()==1)
 739         {
 740           return *getConstPointer();
 741         }
 742       else
 743         throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is allocated but number of elements is not equal to 1 !");
 744     }
 745   else
 746     throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is not allocated !");
 747 }
 748
 749 /*!
 750  * Returns a full copy of \a this. For more info on copying data arrays see
 751  * \ref MEDCouplingArrayBasicsCopyDeep.
 752  *  \return DataArrayDouble * - a new instance of DataArrayDouble. The caller is to
 753  *          delete this array using decrRef() as it is no more needed.
 754  */
 755 DataArrayDouble *DataArrayDouble::deepCopy() const
 756 {
 757   return new DataArrayDouble(*this);
 758 }
 759
 760 /*!
 761  * Checks that \a this array is consistently **increasing** or **decreasing** in value,
 762  * with at least absolute difference value of |\a eps| at each step.
 763  * If not an exception is thrown.
 764  *  \param [in] increasing - if \a true, the array values should be increasing.
 765  *  \param [in] eps - minimal absolute difference between the neighbor values at which
 766  *                    the values are considered different.
 767  *  \throw If sequence of values is not strictly monotonic in agreement with \a
 768  *         increasing arg.
 769  *  \throw If \a this->getNumberOfComponents() != 1.
 770  *  \throw If \a this is not allocated.
 771  */
 772 void DataArrayDouble::checkMonotonic(bool increasing, double eps) const
 773 {
 774   if(!isMonotonic(increasing,eps))
 775     {
 776       if (increasing)
 777         throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not INCREASING monotonic !");
 778       else
 779         throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not DECREASING monotonic !");
 780     }
 781 }
 782
 783 /*!
 784  * Checks that \a this array is consistently **increasing** or **decreasing** in value,
 785  * with at least absolute difference value of |\a eps| at each step.
 786  *  \param [in] increasing - if \a true, array values should be increasing.
 787  *  \param [in] eps - minimal absolute difference between the neighbor values at which
 788  *                    the values are considered different.
 789  *  \return bool - \a true if values change in accordance with \a increasing arg.
 790  *  \throw If \a this->getNumberOfComponents() != 1.
 791  *  \throw If \a this is not allocated.
 792  */
 793 bool DataArrayDouble::isMonotonic(bool increasing, double eps) const
 794 {
 795   checkAllocated();
 796   if(getNumberOfComponents()!=1)
 797     throw INTERP_KERNEL::Exception("DataArrayDouble::isMonotonic : only supported with 'this' array with ONE component !");
 798   mcIdType nbOfElements(getNumberOfTuples());
 799   const double *ptr=getConstPointer();
 800   if(nbOfElements==0)
 801     return true;
 802   double ref=ptr[0];
 803   double absEps=fabs(eps);
 804   if(increasing)
 805     {
 806       for(mcIdType i=1;i<nbOfElements;i++)
 807         {
 808           if(ptr[i]<(ref+absEps))
 809             return false;
 810           ref=ptr[i];
 811         }
 812       return true;
 813     }
 814   else
 815     {
 816       for(mcIdType i=1;i<nbOfElements;i++)
 817         {
 818           if(ptr[i]>(ref-absEps))
 819             return false;
 820           ref=ptr[i];
 821         }
 822       return true;
 823     }
 824 }
 825
 826 void DataArrayDouble::writeVTK(std::ostream& ofs, mcIdType indent, const std::string& nameInFile, DataArrayByte *byteArr) const
 827 {
 828   static const char SPACE[4]={' ',' ',' ',' '};
 829   checkAllocated();
 830   std::string idt(indent,' ');
 831   ofs.precision(17);
 832   ofs << idt << "<DataArray type=\"Float32\" Name=\"" << nameInFile << "\" NumberOfComponents=\"" << getNumberOfComponents() << "\"";
 833   //
 834   bool areAllEmpty(true);
 835   for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
 836     if(!(*it).empty())
 837       areAllEmpty=false;
 838   if(!areAllEmpty)
 839     for(std::size_t i=0;i<_info_on_compo.size();i++)
 840       ofs << " ComponentName" << i << "=\"" << _info_on_compo[i] << "\"";
 841   //
 842   if(byteArr)
 843     {
 844       ofs << " format=\"appended\" offset=\"" << byteArr->getNumberOfTuples() << "\">";
 845       INTERP_KERNEL::AutoPtr<float> tmp(new float[getNbOfElems()]);
 846       float *pt(tmp);
 847       // to make Visual C++ happy : instead of std::copy(begin(),end(),(float *)tmp);
 848       for(const double *src=begin();src!=end();src++,pt++)
 849         *pt=float(*src);
 850       const char *data(reinterpret_cast<const char *>((float *)tmp));
 851       std::size_t sz(getNbOfElems()*sizeof(float));
 852       byteArr->insertAtTheEnd(data,data+sz);
 853       byteArr->insertAtTheEnd(SPACE,SPACE+4);
 854     }
 855   else
 856     {
 857       ofs << " RangeMin=\"" << getMinValueInArray() << "\" RangeMax=\"" << getMaxValueInArray() << "\" format=\"ascii\">\n" << idt;
 858       std::copy(begin(),end(),std::ostream_iterator<double>(ofs," "));
 859     }
 860   ofs << std::endl << idt << "</DataArray>\n";
 861 }
 862
 863 void DataArrayDouble::reprCppStream(const std::string& varName, std::ostream& stream) const
 864 {
 865   mcIdType nbTuples=getNumberOfTuples();
 866   std::size_t nbComp=getNumberOfComponents();
 867   const double *data(getConstPointer());
 868   stream.precision(17);
 869   stream << "DataArrayDouble *" << varName << "=DataArrayDouble::New();" << std::endl;
 870   if(nbTuples*nbComp>=1)
 871     {
 872       stream << "const double " << varName << "Data[" << nbTuples*nbComp << "]={";
 873       std::copy(data,data+nbTuples*nbComp-1,std::ostream_iterator<double>(stream,","));
 874       stream << data[nbTuples*nbComp-1] << "};" << std::endl;
 875       stream << varName << "->useArray(" << varName << "Data,false,CPP_DEALLOC," << nbTuples << "," << nbComp << ");" << std::endl;
 876     }
 877   else
 878     stream << varName << "->alloc(" << nbTuples << "," << nbComp << ");" << std::endl;
 879   stream << varName << "->setName(\"" << getName() << "\");" << std::endl;
 880 }
 881
 882 /*!
 883  * Method that gives a quick overvien of \a this for python.
 884  */
 885 void DataArrayDouble::reprQuickOverview(std::ostream& stream) const
 886 {
 887   static const std::size_t MAX_NB_OF_BYTE_IN_REPR=300;
 888   stream << "DataArrayDouble C++ instance at " << this << ". ";
 889   if(isAllocated())
 890     {
 891       std::size_t nbOfCompo(_info_on_compo.size());
 892       if(nbOfCompo>=1)
 893         {
 894           mcIdType nbOfTuples(getNumberOfTuples());
 895           stream << "Number of tuples : " << nbOfTuples << ". Number of components : " << nbOfCompo << "." << std::endl;
 896           reprQuickOverviewData(stream,MAX_NB_OF_BYTE_IN_REPR);
 897         }
 898       else
 899         stream << "Number of components : 0.";
 900     }
 901   else
 902     stream << "*** No data allocated ****";
 903 }
 904
 905 void DataArrayDouble::reprQuickOverviewData(std::ostream& stream, std::size_t maxNbOfByteInRepr) const
 906 {
 907   const double *data=begin();
 908   mcIdType nbOfTuples(getNumberOfTuples());
 909   std::size_t nbOfCompo(_info_on_compo.size());
 910   std::ostringstream oss2; oss2 << "[";
 911   oss2.precision(17);
 912   std::string oss2Str(oss2.str());
 913   bool isFinished=true;
 914   for(mcIdType i=0;i<nbOfTuples && isFinished;i++)
 915     {
 916       if(nbOfCompo>1)
 917         {
 918           oss2 << "(";
 919           for(std::size_t j=0;j<nbOfCompo;j++,data++)
 920             {
 921               oss2 << *data;
 922               if(j!=nbOfCompo-1) oss2 << ", ";
 923             }
 924           oss2 << ")";
 925         }
 926       else
 927         oss2 << *data++;
 928       if(i!=nbOfTuples-1) oss2 << ", ";
 929       std::string oss3Str(oss2.str());
 930       if(oss3Str.length()<maxNbOfByteInRepr)
 931         oss2Str=oss3Str;
 932       else
 933         isFinished=false;
 934     }
 935   stream << oss2Str;
 936   if(!isFinished)
 937     stream << "... ";
 938   stream << "]";
 939 }
 940
 941 /*!
 942  * Equivalent to DataArrayDouble::isEqual except that if false the reason of
 943  * mismatch is given.
 944  *
 945  * \param [in] other the instance to be compared with \a this
 946  * \param [in] prec the precision to compare numeric data of the arrays.
 947  * \param [out] reason In case of inequality returns the reason.
 948  * \sa DataArrayDouble::isEqual
 949  */
 950 bool DataArrayDouble::isEqualIfNotWhy(const DataArrayDouble& other, double prec, std::string& reason) const
 951 {
 952   if(!areInfoEqualsIfNotWhy(other,reason))
 953     return false;
 954   return _mem.isEqual(other._mem,prec,reason);
 955 }
 956
 957 /*!
 958  * Checks if \a this and another DataArrayDouble are fully equal. For more info see
 959  * \ref MEDCouplingArrayBasicsCompare.
 960  *  \param [in] other - an instance of DataArrayDouble to compare with \a this one.
 961  *  \param [in] prec - precision value to compare numeric data of the arrays.
 962  *  \return bool - \a true if the two arrays are equal, \a false else.
 963  */
 964 bool DataArrayDouble::isEqual(const DataArrayDouble& other, double prec) const
 965 {
 966   std::string tmp;
 967   return isEqualIfNotWhy(other,prec,tmp);
 968 }
 969
 970 /*!
 971  * Checks if values of \a this and another DataArrayDouble are equal. For more info see
 972  * \ref MEDCouplingArrayBasicsCompare.
 973  *  \param [in] other - an instance of DataArrayDouble to compare with \a this one.
 974  *  \param [in] prec - precision value to compare numeric data of the arrays.
 975  *  \return bool - \a true if the values of two arrays are equal, \a false else.
 976  */
 977 bool DataArrayDouble::isEqualWithoutConsideringStr(const DataArrayDouble& other, double prec) const
 978 {
 979   std::string tmp;
 980   return _mem.isEqual(other._mem,prec,tmp);
 981 }
 982
 983 /*!
 984  * This method checks that all tuples in \a other are in \a this.
 985  * If true, the output param \a tupleIds contains the tuples ids of \a this that correspond to tupes in \a this.
 986  * For each i in [ 0 , other->getNumberOfTuples() ) tuple #i in \a other is equal ( regarding input precision \a prec ) to tuple tupleIds[i] in \a this.
 987  *
 988  * \param [in] other - the array having the same number of components than \a this.
 989  * \param [out] tupleIds - the tuple ids containing the same number of tuples than \a other has.
 990  * \sa DataArrayDouble::findCommonTuples
 991  */
 992 bool DataArrayDouble::areIncludedInMe(const DataArrayDouble *other, double prec, DataArrayIdType *&tupleIds) const
 993 {
 994   if(!other)
 995     throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : input array is NULL !");
 996   checkAllocated(); other->checkAllocated();
 997   if(getNumberOfComponents()!=other->getNumberOfComponents())
 998     throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : the number of components does not match !");
 999   MCAuto<DataArrayDouble> a=DataArrayDouble::Aggregate(this,other);
1000   DataArrayIdType *c=0,*ci=0;
1001   a->findCommonTuples(prec,getNumberOfTuples(),c,ci);
1002   MCAuto<DataArrayIdType> cSafe(c),ciSafe(ci);
1003   mcIdType newNbOfTuples=-1;
1004   MCAuto<DataArrayIdType> ids=DataArrayIdType::ConvertIndexArrayToO2N(a->getNumberOfTuples(),c->begin(),ci->begin(),ci->end(),newNbOfTuples);
1005   MCAuto<DataArrayIdType> ret1=ids->selectByTupleIdSafeSlice(getNumberOfTuples(),a->getNumberOfTuples(),1);
1006   tupleIds=ret1.retn();
1007   return newNbOfTuples==getNumberOfTuples();
1008 }
1009
1010 /*!
1011  * Searches for tuples coincident within \a prec tolerance. Each tuple is considered
1012  * as coordinates of a point in getNumberOfComponents()-dimensional space. The
1013  * distance separating two points is computed with the infinite norm.
1014  *
1015  * Indices of coincident tuples are stored in output arrays.
1016  * A pair of arrays (\a comm, \a commIndex) is called "Surjective Format 2".
1017  *
1018  * This method is typically used by MEDCouplingPointSet::findCommonNodes() and
1019  * MEDCouplingUMesh::mergeNodes().
1020  *  \param [in] prec - minimal absolute distance between two tuples (infinite norm) at which they are
1021  *              considered not coincident.
1022  *  \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1023  *              tuples have id strictly lower than \a limitTupleId then they are not returned.
1024  *  \param [out] comm - the array holding ids (== indices) of coincident tuples.
1025  *               \a comm->getNumberOfComponents() == 1.
1026  *               \a comm->getNumberOfTuples() == \a commIndex->back().
1027  *  \param [out] commIndex - the array dividing all indices stored in \a comm into
1028  *               groups of (indices of) coincident tuples. Its every value is a tuple
1029  *               index where a next group of tuples begins. For example the second
1030  *               group of tuples in \a comm is described by following range of indices:
1031  *               [ \a commIndex[1], \a commIndex[2] ). \a commIndex->getNumberOfTuples()-1
1032  *               gives the number of groups of coincident tuples.
1033  *  \throw If \a this is not allocated.
1034  *  \throw If the number of components is not in [1,2,3,4].
1035  *
1036  *  \if ENABLE_EXAMPLES
1037  *  \ref cpp_mcdataarraydouble_findcommontuples "Here is a C++ example".
1038  *
1039  *  \ref py_mcdataarraydouble_findcommontuples  "Here is a Python example".
1040  *  \endif
1041  *  \sa DataArrayInt::ConvertIndexArrayToO2N(), DataArrayDouble::areIncludedInMe
1042  */
1043 void DataArrayDouble::findCommonTuples(double prec, mcIdType limitTupleId, DataArrayIdType *&comm, DataArrayIdType *&commIndex) const
1044 {
1045   checkAllocated();
1046   std::size_t nbOfCompo=getNumberOfComponents();
1047   if ((nbOfCompo<1) || (nbOfCompo>4)) //test before work
1048     throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2, 3 or 4.");
1049
1050   mcIdType nbOfTuples(getNumberOfTuples());
1051   //
1052   MCAuto<DataArrayIdType> c(DataArrayIdType::New()),cI(DataArrayIdType::New()); c->alloc(0,1); cI->pushBackSilent(0);
1053   switch(nbOfCompo)
1054   {
1055     case 4:
1056       findCommonTuplesAlg<4>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1057       break;
1058     case 3:
1059       findCommonTuplesAlg<3>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1060       break;
1061     case 2:
1062       findCommonTuplesAlg<2>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1063       break;
1064     case 1:
1065       findCommonTuplesAlg<1>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1066       break;
1067     default:
1068       throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : nb of components managed are 1,2,3 and 4 ! not implemented for other number of components !");
1069   }
1070   comm=c.retn();
1071   commIndex=cI.retn();
1072 }
1073
1074 /*!
1075  * This methods returns the minimal distance between the two set of points \a this and \a other.
1076  * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1077  * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1078  *
1079  * \param [out] thisTupleId the tuple id in \a this corresponding to the returned minimal distance
1080  * \param [out] otherTupleId the tuple id in \a other corresponding to the returned minimal distance
1081  * \return the minimal distance between the two set of points \a this and \a other.
1082  * \sa DataArrayDouble::findClosestTupleId
1083  */
1084 double DataArrayDouble::minimalDistanceTo(const DataArrayDouble *other, mcIdType& thisTupleId, mcIdType& otherTupleId) const
1085 {
1086   MCAuto<DataArrayIdType> part1=findClosestTupleId(other);
1087   std::size_t nbOfCompo=getNumberOfComponents();
1088   mcIdType otherNbTuples=other->getNumberOfTuples();
1089   const double *thisPt(begin()),*otherPt(other->begin());
1090   const mcIdType *part1Pt(part1->begin());
1091   double ret=std::numeric_limits<double>::max();
1092   for(mcIdType i=0;i<otherNbTuples;i++,part1Pt++,otherPt+=nbOfCompo)
1093     {
1094       double tmp(0.);
1095       for(std::size_t j=0;j<nbOfCompo;j++)
1096         tmp+=(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j])*(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j]);
1097       if(tmp<ret)
1098         { ret=tmp; thisTupleId=*part1Pt; otherTupleId=i; }
1099     }
1100   return sqrt(ret);
1101 }
1102
1103 /*!
1104  * This methods returns for each tuple in \a other which tuple in \a this is the closest.
1105  * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1106  * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1107  *
1108  * \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components.
1109  * \sa DataArrayDouble::minimalDistanceTo
1110  */
1111 DataArrayIdType *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const
1112 {
1113   if(!other)
1114     throw INTERP_KERNEL::Exception("DataArrayDouble::findClosestTupleId : other instance is NULL !");
1115   checkAllocated(); other->checkAllocated();
1116   std::size_t nbOfCompo(getNumberOfComponents());
1117   if(nbOfCompo!=other->getNumberOfComponents())
1118     {
1119       std::ostringstream oss; oss << "DataArrayDouble::findClosestTupleId : number of components in this is " << nbOfCompo;
1120       oss << ", whereas number of components in other is " << other->getNumberOfComponents() << "! Should be equal !";
1121       throw INTERP_KERNEL::Exception(oss.str().c_str());
1122     }
1123   mcIdType nbOfTuples(other->getNumberOfTuples());
1124   mcIdType thisNbOfTuples(getNumberOfTuples());
1125   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(nbOfTuples,1);
1126   double bounds[6];
1127   getMinMaxPerComponent(bounds);
1128   switch(nbOfCompo)
1129   {
1130     case 3:
1131       {
1132         double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4]));
1133         double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta);
1134         double characSize=pow((delta*delta*delta)/((double)thisNbOfTuples),1./3.);
1135         BBTreePts<3,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1136         FindClosestTupleIdAlg<3>(myTree,3.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1137         break;
1138       }
1139     case 2:
1140       {
1141         double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2]));
1142         double delta=std::max(xDelta,yDelta);
1143         double characSize=sqrt(delta/(double)thisNbOfTuples);
1144         BBTreePts<2,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1145         FindClosestTupleIdAlg<2>(myTree,2.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1146         break;
1147       }
1148     case 1:
1149       {
1150         double characSize=fabs(bounds[1]-bounds[0])/FromIdType<double>(thisNbOfTuples);
1151         BBTreePts<1,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1152         FindClosestTupleIdAlg<1>(myTree,1.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1153         break;
1154       }
1155     default:
1156       throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for findClosestTupleId. Must be 1, 2 or 3.");
1157   }
1158   return ret.retn();
1159 }
1160
1161 /*!
1162  * This method expects that \a this and \a otherBBoxFrmt arrays are bounding box arrays ( as the output of MEDCouplingPointSet::getBoundingBoxForBBTree method ).
1163  * This method will return a DataArrayInt array having the same number of tuples than \a this. This returned array tells for each cell in \a this
1164  * how many bounding boxes in \a otherBBoxFrmt.
1165  * So, this method expects that \a this and \a otherBBoxFrmt have the same number of components.
1166  *
1167  * \param [in] otherBBoxFrmt - It is an array .
1168  * \param [in] eps - the absolute precision of the detection. when eps < 0 the bboxes are enlarged so more interactions are detected. Inversely when > 0 the bboxes are stretched.
1169  * \sa MEDCouplingPointSet::getBoundingBoxForBBTree
1170  * \throw If \a this and \a otherBBoxFrmt have not the same number of components.
1171  * \throw If \a this and \a otherBBoxFrmt number of components is not even (BBox format).
1172  */
1173 DataArrayIdType *DataArrayDouble::computeNbOfInteractionsWith(const DataArrayDouble *otherBBoxFrmt, double eps) const
1174 {
1175   if(!otherBBoxFrmt)
1176     throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : input array is NULL !");
1177   if(!isAllocated() || !otherBBoxFrmt->isAllocated())
1178     throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : this and input array must be allocated !");
1179   std::size_t nbOfComp(getNumberOfComponents());
1180   mcIdType nbOfTuples(getNumberOfTuples());
1181   if(nbOfComp!=otherBBoxFrmt->getNumberOfComponents())
1182     {
1183       std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : this number of components (" << nbOfComp << ") must be equal to the number of components of input array (" << otherBBoxFrmt->getNumberOfComponents() << ") !";
1184       throw INTERP_KERNEL::Exception(oss.str().c_str());
1185     }
1186   if(nbOfComp%2!=0)
1187     {
1188       std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : Number of components (" << nbOfComp << ") is not even ! It should be to be compatible with bbox format !";
1189       throw INTERP_KERNEL::Exception(oss.str().c_str());
1190     }
1191   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfTuples,1);
1192   const double *thisBBPtr(begin());
1193   mcIdType *retPtr(ret->getPointer());
1194   switch(nbOfComp/2)
1195   {
1196     case 3:
1197       {
1198         BBTree<3,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1199         for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1200           *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1201         break;
1202       }
1203     case 2:
1204       {
1205         BBTree<2,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1206         for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1207           *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1208         break;
1209       }
1210     case 1:
1211       {
1212         BBTree<1,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1213         for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1214           *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1215         break;
1216       }
1217     default:
1218       throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : space dimension supported are [1,2,3] !");
1219   }
1220
1221   return ret.retn();
1222 }
1223
1224 /*!
1225  * Returns a copy of \a this array by excluding coincident tuples. Each tuple is
1226  * considered as coordinates of a point in getNumberOfComponents()-dimensional
1227  * space. The distance between tuples is computed using norm2. If several tuples are
1228  * not far each from other than \a prec, only one of them remains in the result
1229  * array. The order of tuples in the result array is same as in \a this one except
1230  * that coincident tuples are excluded.
1231  *  \param [in] prec - minimal absolute distance between two tuples at which they are
1232  *              considered not coincident.
1233  *  \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1234  *              tuples have id strictly lower than \a limitTupleId then they are not excluded.
1235  *  \return DataArrayDouble * - the new instance of DataArrayDouble that the caller
1236  *          is to delete using decrRef() as it is no more needed.
1237  *  \throw If \a this is not allocated.
1238  *  \throw If the number of components is not in [1,2,3,4].
1239  *
1240  *  \if ENABLE_EXAMPLES
1241  *  \ref py_mcdataarraydouble_getdifferentvalues "Here is a Python example".
1242  *  \endif
1243  */
1244 DataArrayDouble *DataArrayDouble::getDifferentValues(double prec, mcIdType limitTupleId) const
1245 {
1246   checkAllocated();
1247   DataArrayIdType *c0=0,*cI0=0;
1248   findCommonTuples(prec,limitTupleId,c0,cI0);
1249   MCAuto<DataArrayIdType> c(c0),cI(cI0);
1250   mcIdType newNbOfTuples=-1;
1251   MCAuto<DataArrayIdType> o2n=DataArrayIdType::ConvertIndexArrayToO2N(getNumberOfTuples(),c0->begin(),cI0->begin(),cI0->end(),newNbOfTuples);
1252   return renumberAndReduce(o2n->getConstPointer(),newNbOfTuples);
1253 }
1254
1255 /*!
1256  * Copy all components in a specified order from another DataArrayDouble.
1257  * Both numerical and textual data is copied. The number of tuples in \a this and
1258  * the other array can be different.
1259  *  \param [in] a - the array to copy data from.
1260  *  \param [in] compoIds - sequence of zero based indices of components, data of which is
1261  *              to be copied.
1262  *  \throw If \a a is NULL.
1263  *  \throw If \a compoIds.size() != \a a->getNumberOfComponents().
1264  *  \throw If \a compoIds[i] < 0 or \a compoIds[i] > \a this->getNumberOfComponents().
1265  *
1266  *  \if ENABLE_EXAMPLES
1267  *  \ref py_mcdataarraydouble_setselectedcomponents "Here is a Python example".
1268  *  \endif
1269  */
1270 void DataArrayDouble::setSelectedComponents(const DataArrayDouble *a, const std::vector<std::size_t>& compoIds)
1271 {
1272   if(!a)
1273     throw INTERP_KERNEL::Exception("DataArrayDouble::setSelectedComponents : input DataArrayDouble is NULL !");
1274   checkAllocated();
1275   copyPartOfStringInfoFrom2(compoIds,*a);
1276   std::size_t partOfCompoSz=compoIds.size();
1277   std::size_t nbOfCompo=getNumberOfComponents();
1278   mcIdType nbOfTuples=std::min(getNumberOfTuples(),a->getNumberOfTuples());
1279   const double *ac=a->getConstPointer();
1280   double *nc=getPointer();
1281   for(mcIdType i=0;i<nbOfTuples;i++)
1282     for(std::size_t j=0;j<partOfCompoSz;j++,ac++)
1283       nc[nbOfCompo*i+compoIds[j]]=*ac;
1284 }
1285 /*!
1286  * Checks if 0.0 value is present in \a this array. If it is the case, an exception
1287  * is thrown.
1288  * \throw If zero is found in \a this array.
1289  */
1290 void DataArrayDouble::checkNoNullValues() const
1291 {
1292   const double *tmp=getConstPointer();
1293   mcIdType nbOfElems=getNbOfElems();
1294   const double *where=std::find(tmp,tmp+nbOfElems,0.);
1295   if(where!=tmp+nbOfElems)
1296     throw INTERP_KERNEL::Exception("A value 0.0 have been detected !");
1297 }
1298
1299 /*!
1300  * Computes minimal and maximal value in each component. An output array is filled
1301  * with \c 2 * \a this->getNumberOfComponents() values, so the caller is to allocate
1302  * enough memory before calling this method.
1303  *  \param [out] bounds - array of size at least 2 *\a this->getNumberOfComponents().
1304  *               It is filled as follows:<br>
1305  *               \a bounds[0] = \c min_of_component_0 <br>
1306  *               \a bounds[1] = \c max_of_component_0 <br>
1307  *               \a bounds[2] = \c min_of_component_1 <br>
1308  *               \a bounds[3] = \c max_of_component_1 <br>
1309  *               ...
1310  */
1311 void DataArrayDouble::getMinMaxPerComponent(double *bounds) const
1312 {
1313   checkAllocated();
1314   std::size_t dim=getNumberOfComponents();
1315   for (std::size_t idim=0; idim<dim; idim++)
1316     {
1317       bounds[idim*2]=std::numeric_limits<double>::max();
1318       bounds[idim*2+1]=-std::numeric_limits<double>::max();
1319     }
1320   const double *ptr=getConstPointer();
1321   mcIdType nbOfTuples=getNumberOfTuples();
1322   for(mcIdType i=0;i<nbOfTuples;i++)
1323     {
1324       for(std::size_t idim=0;idim<dim;idim++)
1325         {
1326           if(bounds[idim*2]>ptr[i*dim+idim])
1327             {
1328               bounds[idim*2]=ptr[i*dim+idim];
1329             }
1330           if(bounds[idim*2+1]<ptr[i*dim+idim])
1331             {
1332               bounds[idim*2+1]=ptr[i*dim+idim];
1333             }
1334         }
1335     }
1336 }
1337
1338 /*!
1339  * This method retrieves a newly allocated DataArrayDouble instance having same number of tuples than \a this and twice number of components than \a this
1340  * to store both the min and max per component of each tuples.
1341  * \param [in] epsilon the width of the bbox (identical in each direction) - 0.0 by default
1342  *
1343  * \return a newly created DataArrayDouble instance having \c this->getNumberOfTuples() tuples and 2 * \c this->getNumberOfComponent() components
1344  *
1345  * \throw If \a this is not allocated yet.
1346  */
1347 DataArrayDouble *DataArrayDouble::computeBBoxPerTuple(double epsilon) const
1348 {
1349   checkAllocated();
1350   const double *dataPtr=getConstPointer();
1351   std::size_t nbOfCompo=getNumberOfComponents();
1352   mcIdType nbTuples=getNumberOfTuples();
1353   MCAuto<DataArrayDouble> bbox=DataArrayDouble::New();
1354   bbox->alloc(nbTuples,2*nbOfCompo);
1355   double *bboxPtr=bbox->getPointer();
1356   for(mcIdType i=0;i<nbTuples;i++)
1357     {
1358       for(std::size_t j=0;j<nbOfCompo;j++)
1359         {
1360           bboxPtr[2*nbOfCompo*i+2*j]=dataPtr[nbOfCompo*i+j]-epsilon;
1361           bboxPtr[2*nbOfCompo*i+2*j+1]=dataPtr[nbOfCompo*i+j]+epsilon;
1362         }
1363     }
1364   return bbox.retn();
1365 }
1366
1367 /*!
1368  * For each tuples **t** in \a other, this method retrieves tuples in \a this that are equal to **t**.
1369  * Two tuples are considered equal if the euclidian distance between the two tuples is lower than \a eps.
1370  *
1371  * \param [in] other a DataArrayDouble having same number of components than \a this.
1372  * \param [in] eps absolute precision representing distance (using infinite norm) between 2 tuples behind which 2 tuples are considered equal.
1373  * \param [out] c will contain the set of tuple ids in \a this that are equal to to the tuple ids in \a other contiguously.
1374  *             \a cI allows to extract information in \a c.
1375  * \param [out] cI is an indirection array that allows to extract the data contained in \a c.
1376  *
1377  * \throw In case of:
1378  *  - \a this is not allocated
1379  *  - \a other is not allocated or null
1380  *  - \a this and \a other do not have the same number of components
1381  *  - if number of components of \a this is not in [1,2,3]
1382  *
1383  * \sa MEDCouplingPointSet::getNodeIdsNearPoints, DataArrayDouble::getDifferentValues
1384  */
1385 void DataArrayDouble::computeTupleIdsNearTuples(const DataArrayDouble *other, double eps, DataArrayIdType *& c, DataArrayIdType *& cI) const
1386 {
1387   if(!other)
1388     throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : input pointer other is null !");
1389   checkAllocated();
1390   other->checkAllocated();
1391   std::size_t nbOfCompo=getNumberOfComponents();
1392   std::size_t otherNbOfCompo=other->getNumberOfComponents();
1393   if(nbOfCompo!=otherNbOfCompo)
1394     throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : number of components should be equal between this and other !");
1395   mcIdType nbOfTuplesOther=other->getNumberOfTuples();
1396   mcIdType nbOfTuples=getNumberOfTuples();
1397   MCAuto<DataArrayIdType> cArr(DataArrayIdType::New()),cIArr(DataArrayIdType::New()); cArr->alloc(0,1); cIArr->pushBackSilent(0);
1398   switch(nbOfCompo)
1399   {
1400     case 3:
1401       {
1402         BBTreePts<3,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1403         FindTupleIdsNearTuplesAlg<3>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1404         break;
1405       }
1406     case 2:
1407       {
1408         BBTreePts<2,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1409         FindTupleIdsNearTuplesAlg<2>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1410         break;
1411       }
1412     case 1:
1413       {
1414         BBTreePts<1,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1415         FindTupleIdsNearTuplesAlg<1>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1416         break;
1417       }
1418     default:
1419       throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for computeTupleIdsNearTuples. Must be 1, 2 or 3.");
1420   }
1421   c=cArr.retn(); cI=cIArr.retn();
1422 }
1423
1424 /*!
1425  * This method recenter tuples in \b this in order to be centered at the origin to benefit about the advantages of maximal precision to be around the box
1426  * around origin of 'radius' 1.
1427  *
1428  * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
1429  */
1430 void DataArrayDouble::recenterForMaxPrecision(double eps)
1431 {
1432   checkAllocated();
1433   std::size_t dim=getNumberOfComponents();
1434   std::vector<double> bounds(2*dim);
1435   getMinMaxPerComponent(&bounds[0]);
1436   for(std::size_t i=0;i<dim;i++)
1437     {
1438       double delta=bounds[2*i+1]-bounds[2*i];
1439       double offset=(bounds[2*i]+bounds[2*i+1])/2.;
1440       if(delta>eps)
1441         applyLin(1./delta,-offset/delta,i);
1442       else
1443         applyLin(1.,-offset,i);
1444     }
1445 }
1446
1447 /*!
1448  * Returns the maximal value and all its locations within \a this one-dimensional array.
1449  *  \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1450  *               tuples holding the maximal value. The caller is to delete it using
1451  *               decrRef() as it is no more needed.
1452  *  \return double - the maximal value among all values of \a this array.
1453  *  \throw If \a this->getNumberOfComponents() != 1
1454  *  \throw If \a this->getNumberOfTuples() < 1
1455  */
1456 double DataArrayDouble::getMaxValue2(DataArrayIdType*& tupleIds) const
1457 {
1458   mcIdType tmp;
1459   tupleIds=0;
1460   double ret=getMaxValue(tmp);
1461   tupleIds=findIdsInRange(ret,ret);
1462   return ret;
1463 }
1464
1465 /*!
1466  * Returns the minimal value and all its locations within \a this one-dimensional array.
1467  *  \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1468  *               tuples holding the minimal value. The caller is to delete it using
1469  *               decrRef() as it is no more needed.
1470  *  \return double - the minimal value among all values of \a this array.
1471  *  \throw If \a this->getNumberOfComponents() != 1
1472  *  \throw If \a this->getNumberOfTuples() < 1
1473  */
1474 double DataArrayDouble::getMinValue2(DataArrayIdType*& tupleIds) const
1475 {
1476   mcIdType tmp;
1477   tupleIds=0;
1478   double ret=getMinValue(tmp);
1479   tupleIds=findIdsInRange(ret,ret);
1480   return ret;
1481 }
1482
1483 /*!
1484  * This method returns the number of values in \a this that are equals ( within an absolute precision of \a eps ) to input parameter \a value.
1485  * This method only works for single component array.
1486  *
1487  * \return a value in [ 0, \c this->getNumberOfTuples() )
1488  *
1489  * \throw If \a this is not allocated
1490  *
1491  */
1492 mcIdType DataArrayDouble::count(double value, double eps) const
1493 {
1494   mcIdType ret=0;
1495   checkAllocated();
1496   if(getNumberOfComponents()!=1)
1497     throw INTERP_KERNEL::Exception("DataArrayDouble::count : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1498   const double *vals=begin();
1499   mcIdType nbOfTuples=getNumberOfTuples();
1500   for(mcIdType i=0;i<nbOfTuples;i++,vals++)
1501     if(fabs(*vals-value)<=eps)
1502       ret++;
1503   return ret;
1504 }
1505
1506 /*!
1507  * Returns the average value of \a this one-dimensional array.
1508  *  \return double - the average value over all values of \a this array.
1509  *  \throw If \a this->getNumberOfComponents() != 1
1510  *  \throw If \a this->getNumberOfTuples() < 1
1511  */
1512 double DataArrayDouble::getAverageValue() const
1513 {
1514   if(getNumberOfComponents()!=1)
1515     throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1516   mcIdType nbOfTuples(getNumberOfTuples());
1517   if(nbOfTuples<=0)
1518     throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : array exists but number of tuples must be > 0 !");
1519   const double *vals=getConstPointer();
1520   double ret=std::accumulate(vals,vals+nbOfTuples,0.);
1521   return ret/FromIdType<double>(nbOfTuples);
1522 }
1523
1524 /*!
1525  * Returns the Euclidean norm of the vector defined by \a this array.
1526  *  \return double - the value of the Euclidean norm, i.e.
1527  *          the square root of the inner product of vector.
1528  *  \throw If \a this is not allocated.
1529  */
1530 double DataArrayDouble::norm2() const
1531 {
1532   checkAllocated();
1533   double ret=0.;
1534   std::size_t nbOfElems=getNbOfElems();
1535   const double *pt=getConstPointer();
1536   for(std::size_t i=0;i<nbOfElems;i++,pt++)
1537     ret+=(*pt)*(*pt);
1538   return sqrt(ret);
1539 }
1540
1541 /*!
1542  * Returns the maximum norm of the vector defined by \a this array.
1543  * This method works even if the number of components is different from one.
1544  * If the number of elements in \a this is 0, -1. is returned.
1545  *  \return double - the value of the maximum norm, i.e.
1546  *          the maximal absolute value among values of \a this array (whatever its number of components).
1547  *  \throw If \a this is not allocated.
1548  */
1549 double DataArrayDouble::normMax() const
1550 {
1551   checkAllocated();
1552   double ret(-1.);
1553   std::size_t nbOfElems(getNbOfElems());
1554   const double *pt(getConstPointer());
1555   for(std::size_t i=0;i<nbOfElems;i++,pt++)
1556     {
1557       double val(std::abs(*pt));
1558       if(val>ret)
1559         ret=val;
1560     }
1561   return ret;
1562 }
1563
1564 /*!
1565  * Returns the maximum norm of for each component of \a this array.
1566  * If the number of elements in \a this is 0, -1. is returned.
1567 *  \param [out] res - pointer to an array of result values, of size at least \a
1568  *         this->getNumberOfComponents(), that is to be allocated by the caller.
1569  *  \throw If \a this is not allocated.
1570  */
1571 void DataArrayDouble::normMaxPerComponent(double * res) const
1572 {
1573   checkAllocated();
1574   mcIdType nbOfTuples(getNumberOfTuples());
1575   std::size_t nbOfCompos(getNumberOfComponents());
1576   std::fill(res, res+nbOfCompos, -1.0);
1577   const double *pt(getConstPointer());
1578   for(mcIdType i=0;i<nbOfTuples;i++)
1579     for (std::size_t j=0; j<nbOfCompos; j++, pt++)
1580       {
1581         double val(std::abs(*pt));
1582         if(val>res[j])
1583           res[j]=val;
1584       }
1585 }
1586
1587
1588 /*!
1589  * Returns the minimum norm (absolute value) of the vector defined by \a this array.
1590  * This method works even if the number of components is different from one.
1591  * If the number of elements in \a this is 0, std::numeric_limits<double>::max() is returned.
1592  *  \return double - the value of the minimum norm, i.e.
1593  *          the minimal absolute value among values of \a this array (whatever its number of components).
1594  *  \throw If \a this is not allocated.
1595  */
1596 double DataArrayDouble::normMin() const
1597 {
1598   checkAllocated();
1599   double ret(std::numeric_limits<double>::max());
1600   std::size_t nbOfElems(getNbOfElems());
1601   const double *pt(getConstPointer());
1602   for(std::size_t i=0;i<nbOfElems;i++,pt++)
1603     {
1604       double val(std::abs(*pt));
1605       if(val<ret)
1606         ret=val;
1607     }
1608   return ret;
1609 }
1610
1611 /*!
1612  * Accumulates values of each component of \a this array.
1613  *  \param [out] res - an array of length \a this->getNumberOfComponents(), allocated
1614  *         by the caller, that is filled by this method with sum value for each
1615  *         component.
1616  *  \throw If \a this is not allocated.
1617  */
1618 void DataArrayDouble::accumulate(double *res) const
1619 {
1620   checkAllocated();
1621   const double *ptr=getConstPointer();
1622   mcIdType nbTuple(getNumberOfTuples());
1623   std::size_t nbComps(getNumberOfComponents());
1624   std::fill(res,res+nbComps,0.);
1625   for(mcIdType i=0;i<nbTuple;i++)
1626     std::transform(ptr+i*nbComps,ptr+(i+1)*nbComps,res,res,std::plus<double>());
1627 }
1628
1629 /*!
1630  * This method returns the min distance from an external tuple defined by [ \a tupleBg , \a tupleEnd ) to \a this and
1631  * the first tuple in \a this that matches the returned distance. If there is no tuples in \a this an exception will be thrown.
1632  *
1633  *
1634  * \a this is expected to be allocated and expected to have a number of components equal to the distance from \a tupleBg to
1635  * \a tupleEnd. If not an exception will be thrown.
1636  *
1637  * \param [in] tupleBg start pointer (included) of input external tuple
1638  * \param [in] tupleEnd end pointer (not included) of input external tuple
1639  * \param [out] tupleId the tuple id in \a this that matches the min of distance between \a this and input external tuple
1640  * \return the min distance.
1641  * \sa MEDCouplingUMesh::distanceToPoint
1642  */
1643 double DataArrayDouble::distanceToTuple(const double *tupleBg, const double *tupleEnd, mcIdType& tupleId) const
1644 {
1645   checkAllocated();
1646   mcIdType nbTuple(getNumberOfTuples());
1647   std::size_t nbComps(getNumberOfComponents());
1648   if(nbComps!=(std::size_t)std::distance(tupleBg,tupleEnd))
1649     { std::ostringstream oss; oss << "DataArrayDouble::distanceToTuple : size of input tuple is " << std::distance(tupleBg,tupleEnd) << " should be equal to the number of components in this : " << nbComps << " !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
1650   if(nbTuple==0)
1651     throw INTERP_KERNEL::Exception("DataArrayDouble::distanceToTuple : no tuple in this ! No distance to compute !");
1652   double ret0=std::numeric_limits<double>::max();
1653   tupleId=-1;
1654   const double *work=getConstPointer();
1655   for(mcIdType i=0;i<nbTuple;i++)
1656     {
1657       double val=0.;
1658       for(std::size_t j=0;j<nbComps;j++,work++)
1659         val+=(*work-tupleBg[j])*((*work-tupleBg[j]));
1660       if(val>=ret0)
1661         continue;
1662       else
1663         { ret0=val; tupleId=i; }
1664     }
1665   return sqrt(ret0);
1666 }
1667
1668 /*!
1669  * Accumulate values of the given component of \a this array.
1670  *  \param [in] compId - the index of the component of interest.
1671  *  \return double - a sum value of \a compId-th component.
1672  *  \throw If \a this is not allocated.
1673  *  \throw If \a the condition ( 0 <= \a compId < \a this->getNumberOfComponents() ) is
1674  *         not respected.
1675  */
1676 double DataArrayDouble::accumulate(std::size_t compId) const
1677 {
1678   checkAllocated();
1679   const double *ptr=getConstPointer();
1680   mcIdType nbTuple(getNumberOfTuples());
1681   std::size_t nbComps(getNumberOfComponents());
1682   if(compId>=nbComps)
1683     throw INTERP_KERNEL::Exception("DataArrayDouble::accumulate : Invalid compId specified : No such nb of components !");
1684   double ret=0.;
1685   for(mcIdType i=0;i<nbTuple;i++)
1686     ret+=ptr[i*nbComps+compId];
1687   return ret;
1688 }
1689
1690 /*!
1691  * This method accumulate using addition tuples in \a this using input index array [ \a bgOfIndex, \a endOfIndex ).
1692  * The returned array will have same number of components than \a this and number of tuples equal to
1693  * \c std::distance(bgOfIndex,endOfIndex) \b minus \b one.
1694  *
1695  * The input index array is expected to be ascendingly sorted in which the all referenced ids should be in [0, \c this->getNumberOfTuples).
1696  * This method is quite useful for users that need to put a field on cells to field on nodes on the same mesh without a need of conservation.
1697  *
1698  * \param [in] bgOfIndex - begin (included) of the input index array.
1699  * \param [in] endOfIndex - end (excluded) of the input index array.
1700  * \return DataArrayDouble * - the new instance having the same number of components than \a this.
1701  *
1702  * \throw If bgOfIndex or end is NULL.
1703  * \throw If input index array is not ascendingly sorted.
1704  * \throw If there is an id in [ \a bgOfIndex, \a endOfIndex ) not in [0, \c this->getNumberOfTuples).
1705  * \throw If std::distance(bgOfIndex,endOfIndex)==0.
1706  */
1707 DataArrayDouble *DataArrayDouble::accumulatePerChunck(const mcIdType *bgOfIndex, const mcIdType *endOfIndex) const
1708 {
1709   if(!bgOfIndex || !endOfIndex)
1710     throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : input pointer NULL !");
1711   checkAllocated();
1712   std::size_t nbCompo(getNumberOfComponents());
1713   mcIdType nbOfTuples(getNumberOfTuples());
1714   std::size_t sz=std::distance(bgOfIndex,endOfIndex);
1715   if(sz<1)
1716     throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : invalid size of input index array !");
1717   sz--;
1718   MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(sz,nbCompo);
1719   const mcIdType *w=bgOfIndex;
1720   if(*w<0 || *w>=nbOfTuples)
1721     throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : The first element of the input index not in [0,nbOfTuples) !");
1722   const double *srcPt=begin()+(*w)*nbCompo;
1723   double *tmp=ret->getPointer();
1724   for(std::size_t i=0;i<sz;i++,tmp+=nbCompo,w++)
1725     {
1726       std::fill(tmp,tmp+nbCompo,0.);
1727       if(w[1]>=w[0])
1728         {
1729           for(mcIdType j=w[0];j<w[1];j++,srcPt+=nbCompo)
1730             {
1731               if(j>=0 && j<nbOfTuples)
1732                 std::transform(srcPt,srcPt+nbCompo,tmp,tmp,std::plus<double>());
1733               else
1734                 {
1735                   std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array points to id " << j << " should be in [0," << nbOfTuples << ") !";
1736                   throw INTERP_KERNEL::Exception(oss.str().c_str());
1737                 }
1738             }
1739         }
1740       else
1741         {
1742           std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array is not in ascendingly sorted.";
1743           throw INTERP_KERNEL::Exception(oss.str().c_str());
1744         }
1745     }
1746   ret->copyStringInfoFrom(*this);
1747   return ret.retn();
1748 }
1749
1750 /*!
1751  * This method is close to numpy cumSum except that number of element is equal to \a this->getNumberOfTuples()+1. First element of DataArray returned is equal to 0.
1752  * This method expects that \a this as only one component. The returned array will have \a this->getNumberOfTuples()+1 tuple with also one component.
1753  * The ith element of returned array is equal to the sum of elements in \a this with rank strictly lower than i.
1754  *
1755  * \return DataArrayDouble - A newly built array containing cum sum of \a this.
1756  */
1757 MCAuto<DataArrayDouble> DataArrayDouble::cumSum() const
1758 {
1759   checkAllocated();
1760   checkNbOfComps(1,"DataArrayDouble::cumSum : this is expected to be single component");
1761   mcIdType nbOfTuple(getNumberOfTuples());
1762   MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuple+1,1);
1763   double *ptr(ret->getPointer());
1764   ptr[0]=0.;
1765   const double *thisPtr(begin());
1766   for(mcIdType i=0;i<nbOfTuple;i++)
1767     ptr[i+1]=ptr[i]+thisPtr[i];
1768   return ret;
1769 }
1770
1771 /*!
1772  * Converts each 2D point defined by the tuple of \a this array from the Polar to the
1773  * Cartesian coordinate system. The two components of the tuple of \a this array are
1774  * considered to contain (1) radius and (2) angle of the point in the Polar CS.
1775  *  \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1776  *          contains X and Y coordinates of the point in the Cartesian CS. The caller
1777  *          is to delete this array using decrRef() as it is no more needed. The array
1778  *          does not contain any textual info on components.
1779  *  \throw If \a this->getNumberOfComponents() != 2.
1780  * \sa fromCartToPolar
1781  */
1782 DataArrayDouble *DataArrayDouble::fromPolarToCart() const
1783 {
1784   checkAllocated();
1785   std::size_t nbOfComp(getNumberOfComponents());
1786   if(nbOfComp!=2)
1787     throw INTERP_KERNEL::Exception("DataArrayDouble::fromPolarToCart : must be an array with exactly 2 components !");
1788   mcIdType nbOfTuple(getNumberOfTuples());
1789   DataArrayDouble *ret(DataArrayDouble::New());
1790   ret->alloc(nbOfTuple,2);
1791   double *w(ret->getPointer());
1792   const double *wIn(getConstPointer());
1793   for(mcIdType i=0;i<nbOfTuple;i++,w+=2,wIn+=2)
1794     {
1795       w[0]=wIn[0]*cos(wIn[1]);
1796       w[1]=wIn[0]*sin(wIn[1]);
1797     }
1798   return ret;
1799 }
1800
1801 /*!
1802  * Converts each 3D point defined by the tuple of \a this array from the Cylindrical to
1803  * the Cartesian coordinate system. The three components of the tuple of \a this array
1804  * are considered to contain (1) radius, (2) azimuth and (3) altitude of the point in
1805  * the Cylindrical CS.
1806  *  \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1807  *          contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1808  *          on the third component is copied from \a this array. The caller
1809  *          is to delete this array using decrRef() as it is no more needed.
1810  *  \throw If \a this->getNumberOfComponents() != 3.
1811  * \sa fromCartToCyl
1812  */
1813 DataArrayDouble *DataArrayDouble::fromCylToCart() const
1814 {
1815   checkAllocated();
1816   std::size_t nbOfComp(getNumberOfComponents());
1817   if(nbOfComp!=3)
1818     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCylToCart : must be an array with exactly 3 components !");
1819   mcIdType nbOfTuple(getNumberOfTuples());
1820   DataArrayDouble *ret(DataArrayDouble::New());
1821   ret->alloc(getNumberOfTuples(),3);
1822   double *w(ret->getPointer());
1823   const double *wIn(getConstPointer());
1824   for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1825     {
1826       w[0]=wIn[0]*cos(wIn[1]);
1827       w[1]=wIn[0]*sin(wIn[1]);
1828       w[2]=wIn[2];
1829     }
1830   ret->setInfoOnComponent(2,getInfoOnComponent(2));
1831   return ret;
1832 }
1833
1834 /*!
1835  * Converts each 3D point defined by the tuple of \a this array from the Spherical to
1836  * the Cartesian coordinate system. The three components of the tuple of \a this array
1837  * are considered to contain (1) radius, (2) polar angle and (3) azimuthal angle of the
1838  * point in the Cylindrical CS.
1839  *  \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1840  *          contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1841  *          on the third component is copied from \a this array. The caller
1842  *          is to delete this array using decrRef() as it is no more needed.
1843  *  \throw If \a this->getNumberOfComponents() != 3.
1844  * \sa fromCartToSpher
1845  */
1846 DataArrayDouble *DataArrayDouble::fromSpherToCart() const
1847 {
1848   checkAllocated();
1849   std::size_t nbOfComp(getNumberOfComponents());
1850   if(nbOfComp!=3)
1851     throw INTERP_KERNEL::Exception("DataArrayDouble::fromSpherToCart : must be an array with exactly 3 components !");
1852   mcIdType nbOfTuple(getNumberOfTuples());
1853   DataArrayDouble *ret(DataArrayDouble::New());
1854   ret->alloc(getNumberOfTuples(),3);
1855   double *w(ret->getPointer());
1856   const double *wIn(getConstPointer());
1857   for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1858     {
1859       w[0]=wIn[0]*cos(wIn[2])*sin(wIn[1]);
1860       w[1]=wIn[0]*sin(wIn[2])*sin(wIn[1]);
1861       w[2]=wIn[0]*cos(wIn[1]);
1862     }
1863   return ret;
1864 }
1865
1866 /*!
1867  * This method returns a new array containing the same number of tuples than \a this. To do this, this method needs \a at parameter to specify the convention of \a this.
1868  * All the tuples of the returned array will be in cartesian sense. So if \a at equals to AX_CART the returned array is basically a deep copy of \a this.
1869  * If \a at equals to AX_CYL the returned array will be the result of operation cylindric to cartesian of \a this...
1870  *
1871  * \param [in] atOfThis - The axis type of \a this.
1872  * \return DataArrayDouble * - the new instance of DataArrayDouble (that must be dealed by caller) containing the result of the cartesianizification of \a this.
1873  */
1874 DataArrayDouble *DataArrayDouble::cartesianize(MEDCouplingAxisType atOfThis) const
1875 {
1876   checkAllocated();
1877   std::size_t nbOfComp(getNumberOfComponents());
1878   MCAuto<DataArrayDouble> ret;
1879   switch(atOfThis)
1880     {
1881     case AX_CART:
1882       ret=deepCopy();
1883       break;
1884     case AX_CYL:
1885       if(nbOfComp==3)
1886         {
1887           ret=fromCylToCart();
1888           break;
1889         }
1890       if(nbOfComp==2)
1891         {
1892           ret=fromPolarToCart();
1893           break;
1894         }
1895       else
1896         throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1897     case AX_SPHER:
1898       if(nbOfComp==3)
1899         {
1900           ret=fromSpherToCart();
1901           break;
1902         }
1903       if(nbOfComp==2)
1904         {
1905           ret=fromPolarToCart();
1906           break;
1907         }
1908       else
1909         throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1910     default:
1911       throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : not recognized axis type ! Only AX_CART, AX_CYL and AX_SPHER supported !");
1912     }
1913   ret->copyStringInfoFrom(*this);
1914   return ret.retn();
1915 }
1916
1917 /*!
1918  * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to polar.
1919  * This method expects that \a this has exactly 2 components.
1920  * \sa fromPolarToCart
1921  */
1922 DataArrayDouble *DataArrayDouble::fromCartToPolar() const
1923 {
1924   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1925   checkAllocated();
1926   std::size_t nbOfComp(getNumberOfComponents());
1927   mcIdType nbTuples(getNumberOfTuples());
1928   if(nbOfComp!=2)
1929     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToPolar : must be an array with exactly 2 components !");
1930   ret->alloc(nbTuples,2);
1931   double *retPtr(ret->getPointer());
1932   const double *ptr(begin());
1933   for(mcIdType i=0;i<nbTuples;i++,ptr+=2,retPtr+=2)
1934     {
1935       retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1936       retPtr[1]=atan2(ptr[1],ptr[0]);
1937     }
1938   return ret.retn();
1939 }
1940
1941 /*!
1942  * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical.
1943  * This method expects that \a this has exactly 3 components.
1944  * \sa fromCylToCart
1945  */
1946 DataArrayDouble *DataArrayDouble::fromCartToCyl() const
1947 {
1948   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1949   checkAllocated();
1950   std::size_t nbOfComp(getNumberOfComponents());
1951   mcIdType nbTuples(getNumberOfTuples());
1952   if(nbOfComp!=3)
1953     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCyl : must be an array with exactly 3 components !");
1954   ret->alloc(nbTuples,3);
1955   double *retPtr(ret->getPointer());
1956   const double *ptr(begin());
1957   for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1958     {
1959       retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1960       retPtr[1]=atan2(ptr[1],ptr[0]);
1961       retPtr[2]=ptr[2];
1962     }
1963   return ret.retn();
1964 }
1965
1966 /*!
1967  * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to spherical coordinates.
1968  * \sa fromSpherToCart
1969  */
1970 DataArrayDouble *DataArrayDouble::fromCartToSpher() const
1971 {
1972   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1973   checkAllocated();
1974   std::size_t nbOfComp(getNumberOfComponents());
1975   mcIdType nbTuples(getNumberOfTuples());
1976   if(nbOfComp!=3)
1977     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToSpher : must be an array with exactly 3 components !");
1978   ret->alloc(nbTuples,3);
1979   double *retPtr(ret->getPointer());
1980   const double *ptr(begin());
1981   for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1982     {
1983       retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]+ptr[2]*ptr[2]);
1984       retPtr[1]=acos(ptr[2]/retPtr[0]);
1985       retPtr[2]=atan2(ptr[1],ptr[0]);
1986     }
1987   return ret.retn();
1988 }
1989
1990 /*!
1991  * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical relative to the given \a center and a \a vector.
1992  * This method expects that \a this has exactly 3 components.
1993  * \sa MEDCouplingFieldDouble::computeVectorFieldCyl
1994  */
1995 DataArrayDouble *DataArrayDouble::fromCartToCylGiven(const DataArrayDouble *coords, const double center[3], const double vect[3]) const
1996 {
1997   if(!coords)
1998     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : input coords are NULL !");
1999   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
2000   checkAllocated(); coords->checkAllocated();
2001   std::size_t nbOfComp(getNumberOfComponents());
2002   mcIdType nbTuples(getNumberOfTuples());
2003   if(nbOfComp!=3)
2004     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : must be an array with exactly 3 components !");
2005   if(coords->getNumberOfComponents()!=3)
2006     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have exactly 3 components !");
2007   if(coords->getNumberOfTuples()!=nbTuples)
2008     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have the same number of tuples !");
2009   ret->alloc(nbTuples,nbOfComp);
2010   double magOfVect(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
2011   if(magOfVect<1e-12)
2012     throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : magnitude of vect is too low !");
2013   double Ur[3],Uteta[3],Uz[3],*retPtr(ret->getPointer());
2014   const double *coo(coords->begin()),*vectField(begin());
2015   std::transform(vect,vect+3,Uz,std::bind2nd(std::multiplies<double>(),1./magOfVect));
2016   for(mcIdType i=0;i<nbTuples;i++,vectField+=3,retPtr+=3,coo+=3)
2017     {
2018       std::transform(coo,coo+3,center,Ur,std::minus<double>());
2019       Uteta[0]=Uz[1]*Ur[2]-Uz[2]*Ur[1]; Uteta[1]=Uz[2]*Ur[0]-Uz[0]*Ur[2]; Uteta[2]=Uz[0]*Ur[1]-Uz[1]*Ur[0];
2020       double magOfTeta(sqrt(Uteta[0]*Uteta[0]+Uteta[1]*Uteta[1]+Uteta[2]*Uteta[2]));
2021       std::transform(Uteta,Uteta+3,Uteta,std::bind2nd(std::multiplies<double>(),1./magOfTeta));
2022       Ur[0]=Uteta[1]*Uz[2]-Uteta[2]*Uz[1]; Ur[1]=Uteta[2]*Uz[0]-Uteta[0]*Uz[2]; Ur[2]=Uteta[0]*Uz[1]-Uteta[1]*Uz[0];
2023       retPtr[0]=Ur[0]*vectField[0]+Ur[1]*vectField[1]+Ur[2]*vectField[2];
2024       retPtr[1]=Uteta[0]*vectField[0]+Uteta[1]*vectField[1]+Uteta[2]*vectField[2];
2025       retPtr[2]=Uz[0]*vectField[0]+Uz[1]*vectField[1]+Uz[2]*vectField[2];
2026     }
2027   ret->copyStringInfoFrom(*this);
2028   return ret.retn();
2029 }
2030
2031 /*!
2032  * Computes the doubly contracted product of every tensor defined by the tuple of \a this
2033  * array containing 6 components.
2034  *  \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2035  *          is calculated from the tuple <em>(t)</em> of \a this array as follows:
2036  *          \f$ t[0]^2+t[1]^2+t[2]^2+2*t[3]^2+2*t[4]^2+2*t[5]^2\f$.
2037  *         The caller is to delete this result array using decrRef() as it is no more needed.
2038  *  \throw If \a this->getNumberOfComponents() != 6.
2039  */
2040 DataArrayDouble *DataArrayDouble::doublyContractedProduct() const
2041 {
2042   checkAllocated();
2043   std::size_t nbOfComp(getNumberOfComponents());
2044   if(nbOfComp!=6)
2045     throw INTERP_KERNEL::Exception("DataArrayDouble::doublyContractedProduct : must be an array with exactly 6 components !");
2046   DataArrayDouble *ret=DataArrayDouble::New();
2047   mcIdType nbOfTuple=getNumberOfTuples();
2048   ret->alloc(nbOfTuple,1);
2049   const double *src=getConstPointer();
2050   double *dest=ret->getPointer();
2051   for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2052     *dest=src[0]*src[0]+src[1]*src[1]+src[2]*src[2]+2.*src[3]*src[3]+2.*src[4]*src[4]+2.*src[5]*src[5];
2053   return ret;
2054 }
2055
2056 /*!
2057  * Computes the determinant of every square matrix defined by the tuple of \a this
2058  * array, which contains either 4, 6 or 9 components. The case of 6 components
2059  * corresponds to that of the upper triangular matrix.
2060  *  \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2061  *          is the determinant of matrix of the corresponding tuple of \a this array.
2062  *          The caller is to delete this result array using decrRef() as it is no more
2063  *          needed.
2064  *  \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2065  */
2066 DataArrayDouble *DataArrayDouble::determinant() const
2067 {
2068   checkAllocated();
2069   DataArrayDouble *ret=DataArrayDouble::New();
2070   mcIdType nbOfTuple=getNumberOfTuples();
2071   ret->alloc(nbOfTuple,1);
2072   const double *src=getConstPointer();
2073   double *dest=ret->getPointer();
2074   switch(getNumberOfComponents())
2075   {
2076     case 6:
2077       for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2078         *dest=src[0]*src[1]*src[2]+2.*src[4]*src[5]*src[3]-src[0]*src[4]*src[4]-src[2]*src[3]*src[3]-src[1]*src[5]*src[5];
2079       return ret;
2080     case 4:
2081       for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2082         *dest=src[0]*src[3]-src[1]*src[2];
2083       return ret;
2084     case 9:
2085       for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2086         *dest=src[0]*src[4]*src[8]+src[1]*src[5]*src[6]+src[2]*src[3]*src[7]-src[0]*src[5]*src[7]-src[1]*src[3]*src[8]-src[2]*src[4]*src[6];
2087       return ret;
2088     default:
2089       ret->decrRef();
2090       throw INTERP_KERNEL::Exception("DataArrayDouble::determinant : Invalid number of components ! must be in 4,6,9 !");
2091   }
2092 }
2093
2094 /*!
2095  * Computes 3 eigenvalues of every upper triangular matrix defined by the tuple of
2096  * \a this array, which contains 6 components.
2097  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing 3
2098  *          components, whose each tuple contains the eigenvalues of the matrix of
2099  *          corresponding tuple of \a this array.
2100  *          The caller is to delete this result array using decrRef() as it is no more
2101  *          needed.
2102  *  \throw If \a this->getNumberOfComponents() != 6.
2103  */
2104 DataArrayDouble *DataArrayDouble::eigenValues() const
2105 {
2106   checkAllocated();
2107   std::size_t nbOfComp=getNumberOfComponents();
2108   if(nbOfComp!=6)
2109     throw INTERP_KERNEL::Exception("DataArrayDouble::eigenValues : must be an array with exactly 6 components !");
2110   DataArrayDouble *ret=DataArrayDouble::New();
2111   mcIdType nbOfTuple=getNumberOfTuples();
2112   ret->alloc(nbOfTuple,3);
2113   const double *src=getConstPointer();
2114   double *dest=ret->getPointer();
2115   for(mcIdType i=0;i<nbOfTuple;i++,dest+=3,src+=6)
2116     INTERP_KERNEL::computeEigenValues6(src,dest);
2117   return ret;
2118 }
2119
2120 /*!
2121  * Computes 3 eigenvectors of every upper triangular matrix defined by the tuple of
2122  * \a this array, which contains 6 components.
2123  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing 9
2124  *          components, whose each tuple contains 3 eigenvectors of the matrix of
2125  *          corresponding tuple of \a this array.
2126  *          The caller is to delete this result array using decrRef() as it is no more
2127  *          needed.
2128  *  \throw If \a this->getNumberOfComponents() != 6.
2129  */
2130 DataArrayDouble *DataArrayDouble::eigenVectors() const
2131 {
2132   checkAllocated();
2133   std::size_t nbOfComp=getNumberOfComponents();
2134   if(nbOfComp!=6)
2135     throw INTERP_KERNEL::Exception("DataArrayDouble::eigenVectors : must be an array with exactly 6 components !");
2136   DataArrayDouble *ret=DataArrayDouble::New();
2137   mcIdType nbOfTuple=getNumberOfTuples();
2138   ret->alloc(nbOfTuple,9);
2139   const double *src=getConstPointer();
2140   double *dest=ret->getPointer();
2141   for(mcIdType i=0;i<nbOfTuple;i++,src+=6)
2142     {
2143       double tmp[3];
2144       INTERP_KERNEL::computeEigenValues6(src,tmp);
2145       for(mcIdType j=0;j<3;j++,dest+=3)
2146         INTERP_KERNEL::computeEigenVectorForEigenValue6(src,tmp[j],1e-12,dest);
2147     }
2148   return ret;
2149 }
2150
2151 /*!
2152  * Computes the inverse matrix of every matrix defined by the tuple of \a this
2153  * array, which contains either 4, 6 or 9 components. The case of 6 components
2154  * corresponds to that of the upper triangular matrix.
2155  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2156  *          same number of components as \a this one, whose each tuple is the inverse
2157  *          matrix of the matrix of corresponding tuple of \a this array.
2158  *          The caller is to delete this result array using decrRef() as it is no more
2159  *          needed.
2160  *  \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2161  */
2162 DataArrayDouble *DataArrayDouble::inverse() const
2163 {
2164   checkAllocated();
2165   std::size_t nbOfComp=getNumberOfComponents();
2166   if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2167     throw INTERP_KERNEL::Exception("DataArrayDouble::inversion : must be an array with 4,6 or 9 components !");
2168   DataArrayDouble *ret=DataArrayDouble::New();
2169   mcIdType nbOfTuple=getNumberOfTuples();
2170   ret->alloc(nbOfTuple,nbOfComp);
2171   const double *src=getConstPointer();
2172   double *dest=ret->getPointer();
2173   if(nbOfComp==6)
2174     for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2175       {
2176         double det=src[0]*src[1]*src[2]+2.*src[4]*src[5]*src[3]-src[0]*src[4]*src[4]-src[2]*src[3]*src[3]-src[1]*src[5]*src[5];
2177         dest[0]=(src[1]*src[2]-src[4]*src[4])/det;
2178         dest[1]=(src[0]*src[2]-src[5]*src[5])/det;
2179         dest[2]=(src[0]*src[1]-src[3]*src[3])/det;
2180         dest[3]=(src[5]*src[4]-src[3]*src[2])/det;
2181         dest[4]=(src[5]*src[3]-src[0]*src[4])/det;
2182         dest[5]=(src[3]*src[4]-src[1]*src[5])/det;
2183       }
2184   else if(nbOfComp==4)
2185     for(mcIdType i=0;i<nbOfTuple;i++,dest+=4,src+=4)
2186       {
2187         double det=src[0]*src[3]-src[1]*src[2];
2188         dest[0]=src[3]/det;
2189         dest[1]=-src[1]/det;
2190         dest[2]=-src[2]/det;
2191         dest[3]=src[0]/det;
2192       }
2193   else
2194     for(mcIdType i=0;i<nbOfTuple;i++,dest+=9,src+=9)
2195       {
2196         double det=src[0]*src[4]*src[8]+src[1]*src[5]*src[6]+src[2]*src[3]*src[7]-src[0]*src[5]*src[7]-src[1]*src[3]*src[8]-src[2]*src[4]*src[6];
2197         dest[0]=(src[4]*src[8]-src[7]*src[5])/det;
2198         dest[1]=(src[7]*src[2]-src[1]*src[8])/det;
2199         dest[2]=(src[1]*src[5]-src[4]*src[2])/det;
2200         dest[3]=(src[6]*src[5]-src[3]*src[8])/det;
2201         dest[4]=(src[0]*src[8]-src[6]*src[2])/det;
2202         dest[5]=(src[2]*src[3]-src[0]*src[5])/det;
2203         dest[6]=(src[3]*src[7]-src[6]*src[4])/det;
2204         dest[7]=(src[6]*src[1]-src[0]*src[7])/det;
2205         dest[8]=(src[0]*src[4]-src[1]*src[3])/det;
2206       }
2207   return ret;
2208 }
2209
2210 /*!
2211  * Computes the trace of every matrix defined by the tuple of \a this
2212  * array, which contains either 4, 6 or 9 components. The case of 6 components
2213  * corresponds to that of the upper triangular matrix.
2214  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing
2215  *          1 component, whose each tuple is the trace of
2216  *          the matrix of corresponding tuple of \a this array.
2217  *          The caller is to delete this result array using decrRef() as it is no more
2218  *          needed.
2219  *  \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2220  */
2221 DataArrayDouble *DataArrayDouble::trace() const
2222 {
2223   checkAllocated();
2224   std::size_t nbOfComp=getNumberOfComponents();
2225   if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2226     throw INTERP_KERNEL::Exception("DataArrayDouble::trace : must be an array with 4,6 or 9 components !");
2227   DataArrayDouble *ret=DataArrayDouble::New();
2228   mcIdType nbOfTuple=getNumberOfTuples();
2229   ret->alloc(nbOfTuple,1);
2230   const double *src=getConstPointer();
2231   double *dest=ret->getPointer();
2232   if(nbOfComp==6)
2233     for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2234       *dest=src[0]+src[1]+src[2];
2235   else if(nbOfComp==4)
2236     for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2237       *dest=src[0]+src[3];
2238   else
2239     for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2240       *dest=src[0]+src[4]+src[8];
2241   return ret;
2242 }
2243
2244 /*!
2245  * Computes the stress deviator tensor of every stress tensor defined by the tuple of
2246  * \a this array, which contains 6 components.
2247  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2248  *          same number of components and tuples as \a this array.
2249  *          The caller is to delete this result array using decrRef() as it is no more
2250  *          needed.
2251  *  \throw If \a this->getNumberOfComponents() != 6.
2252  */
2253 DataArrayDouble *DataArrayDouble::deviator() const
2254 {
2255   checkAllocated();
2256   std::size_t nbOfComp=getNumberOfComponents();
2257   if(nbOfComp!=6)
2258     throw INTERP_KERNEL::Exception("DataArrayDouble::deviator : must be an array with exactly 6 components !");
2259   DataArrayDouble *ret=DataArrayDouble::New();
2260   mcIdType nbOfTuple=getNumberOfTuples();
2261   ret->alloc(nbOfTuple,6);
2262   const double *src=getConstPointer();
2263   double *dest=ret->getPointer();
2264   for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2265     {
2266       double tr=(src[0]+src[1]+src[2])/3.;
2267       dest[0]=src[0]-tr;
2268       dest[1]=src[1]-tr;
2269       dest[2]=src[2]-tr;
2270       dest[3]=src[3];
2271       dest[4]=src[4];
2272       dest[5]=src[5];
2273     }
2274   return ret;
2275 }
2276
2277 /*!
2278  * Computes the magnitude of every vector defined by the tuple of
2279  * \a this array.
2280  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2281  *          same number of tuples as \a this array and one component.
2282  *          The caller is to delete this result array using decrRef() as it is no more
2283  *          needed.
2284  *  \throw If \a this is not allocated.
2285  */
2286 DataArrayDouble *DataArrayDouble::magnitude() const
2287 {
2288   checkAllocated();
2289   std::size_t nbOfComp=getNumberOfComponents();
2290   DataArrayDouble *ret=DataArrayDouble::New();
2291   mcIdType nbOfTuple=getNumberOfTuples();
2292   ret->alloc(nbOfTuple,1);
2293   const double *src=getConstPointer();
2294   double *dest=ret->getPointer();
2295   for(mcIdType i=0;i<nbOfTuple;i++,dest++)
2296     {
2297       double sum=0.;
2298       for(std::size_t j=0;j<nbOfComp;j++,src++)
2299         sum+=(*src)*(*src);
2300       *dest=sqrt(sum);
2301     }
2302   return ret;
2303 }
2304
2305 /*!
2306  * Computes the maximal value within every tuple of \a this array.
2307  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2308  *          same number of tuples as \a this array and one component.
2309  *          The caller is to delete this result array using decrRef() as it is no more
2310  *          needed.
2311  *  \throw If \a this is not allocated.
2312  *  \sa DataArrayDouble::maxPerTupleWithCompoId
2313  */
2314 DataArrayDouble *DataArrayDouble::maxPerTuple() const
2315 {
2316   checkAllocated();
2317   std::size_t nbOfComp(getNumberOfComponents());
2318   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2319   mcIdType nbOfTuple(getNumberOfTuples());
2320   ret->alloc(nbOfTuple,1);
2321   const double *src=getConstPointer();
2322   double *dest=ret->getPointer();
2323   for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=nbOfComp)
2324     *dest=*std::max_element(src,src+nbOfComp);
2325   return ret.retn();
2326 }
2327
2328 /*!
2329  * Computes the maximal value within every tuple of \a this array and it returns the first component
2330  * id for each tuple that corresponds to the maximal value within the tuple.
2331  *
2332  *  \param [out] compoIdOfMaxPerTuple - the new new instance of DataArrayInt containing the
2333  *          same number of tuples and only one component.
2334  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2335  *          same number of tuples as \a this array and one component.
2336  *          The caller is to delete this result array using decrRef() as it is no more
2337  *          needed.
2338  *  \throw If \a this is not allocated.
2339  *  \sa DataArrayDouble::maxPerTuple
2340  */
2341 DataArrayDouble *DataArrayDouble::maxPerTupleWithCompoId(DataArrayIdType* &compoIdOfMaxPerTuple) const
2342 {
2343   checkAllocated();
2344   std::size_t nbOfComp(getNumberOfComponents());
2345   MCAuto<DataArrayDouble> ret0=DataArrayDouble::New();
2346   MCAuto<DataArrayIdType> ret1=DataArrayIdType::New();
2347   mcIdType nbOfTuple=getNumberOfTuples();
2348   ret0->alloc(nbOfTuple,1); ret1->alloc(nbOfTuple,1);
2349   const double *src=getConstPointer();
2350   double *dest=ret0->getPointer(); mcIdType *dest1=ret1->getPointer();
2351   for(mcIdType i=0;i<nbOfTuple;i++,dest++,dest1++,src+=nbOfComp)
2352     {
2353       const double *loc=std::max_element(src,src+nbOfComp);
2354       *dest=*loc;
2355       *dest1=ToIdType(std::distance(src,loc));
2356     }
2357   compoIdOfMaxPerTuple=ret1.retn();
2358   return ret0.retn();
2359 }
2360
2361 /*!
2362  * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c this->getNumberOfTuples() tuples.
2363  * \n This returned array contains the euclidian distance for each tuple in \a this.
2364  * \n So the returned array can be seen as a dense symmetrical matrix whose diagonal elements are equal to 0.
2365  * \n The returned array has only one component (and **not** \c this->getNumberOfTuples() components to avoid the useless memory consumption due to components info in returned DataArrayDouble)
2366  *
2367  * \warning use this method with care because it can leads to big amount of consumed memory !
2368  *
2369  * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2370  *
2371  * \throw If \a this is not allocated.
2372  *
2373  * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrixWith
2374  */
2375 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrix() const
2376 {
2377   checkAllocated();
2378   std::size_t nbOfComp(getNumberOfComponents());
2379   mcIdType nbOfTuples(getNumberOfTuples());
2380   const double *inData=getConstPointer();
2381   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2382   ret->alloc(nbOfTuples*nbOfTuples,1);
2383   double *outData=ret->getPointer();
2384   for(mcIdType i=0;i<nbOfTuples;i++)
2385     {
2386       outData[i*nbOfTuples+i]=0.;
2387       for(mcIdType j=i+1;j<nbOfTuples;j++)
2388         {
2389           double dist=0.;
2390           for(std::size_t k=0;k<nbOfComp;k++)
2391             { double delta=inData[i*nbOfComp+k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2392           dist=sqrt(dist);
2393           outData[i*nbOfTuples+j]=dist;
2394           outData[j*nbOfTuples+i]=dist;
2395         }
2396     }
2397   return ret.retn();
2398 }
2399
2400 /*!
2401  * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c other->getNumberOfTuples() tuples.
2402  * \n This returned array contains the euclidian distance for each tuple in \a other with each tuple in \a this.
2403  * \n So the returned array can be seen as a dense rectangular matrix with \c other->getNumberOfTuples() rows and \c this->getNumberOfTuples() columns.
2404  * \n Output rectangular matrix is sorted along rows.
2405  * \n The returned array has only one component (and **not** \c this->getNumberOfTuples() components to avoid the useless memory consumption due to components info in returned DataArrayDouble)
2406  *
2407  * \warning use this method with care because it can leads to big amount of consumed memory !
2408  *
2409  * \param [in] other DataArrayDouble instance having same number of components than \a this.
2410  * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2411  *
2412  * \throw If \a this is not allocated, or if \a other is null or if \a other is not allocated, or if number of components of \a other and \a this differs.
2413  *
2414  * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrix
2415  */
2416 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const
2417 {
2418   if(!other)
2419     throw INTERP_KERNEL::Exception("DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : input parameter is null !");
2420   checkAllocated();
2421   other->checkAllocated();
2422   std::size_t nbOfComp(getNumberOfComponents());
2423   std::size_t otherNbOfComp(other->getNumberOfComponents());
2424   if(nbOfComp!=otherNbOfComp)
2425     {
2426       std::ostringstream oss; oss << "DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : this nb of compo=" << nbOfComp << " and other nb of compo=" << otherNbOfComp << ". It should match !";
2427       throw INTERP_KERNEL::Exception(oss.str().c_str());
2428     }
2429   mcIdType nbOfTuples(getNumberOfTuples());
2430   mcIdType otherNbOfTuples(other->getNumberOfTuples());
2431   const double *inData=getConstPointer();
2432   const double *inDataOther=other->getConstPointer();
2433   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2434   ret->alloc(otherNbOfTuples*nbOfTuples,1);
2435   double *outData=ret->getPointer();
2436   for(mcIdType i=0;i<otherNbOfTuples;i++,inDataOther+=nbOfComp)
2437     {
2438       for(mcIdType j=0;j<nbOfTuples;j++)
2439         {
2440           double dist=0.;
2441           for(std::size_t k=0;k<nbOfComp;k++)
2442             { double delta=inDataOther[k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2443           dist=sqrt(dist);
2444           outData[i*nbOfTuples+j]=dist;
2445         }
2446     }
2447   return ret.retn();
2448 }
2449
2450 /*!
2451  * This method expects that \a this stores 3 tuples containing 2 components each.
2452  * Each of this tuples represent a point into 2D space.
2453  * This method tries to find an arc of circle starting from first point (tuple) to 2nd and middle point (tuple) along 3nd and last point (tuple).
2454  * If such arc of circle exists, the corresponding center, radius of circle is returned. And additionnaly the length of arc expressed with an \a ang output variable in ]0,2*pi[.
2455  *
2456  *  \throw If \a this is not allocated.
2457  *  \throw If \a this has not 3 tuples of 2 components
2458  *  \throw If tuples/points in \a this are aligned
2459  */
2460 void DataArrayDouble::asArcOfCircle(double center[2], double& radius, double& ang) const
2461 {
2462   checkAllocated();
2463   INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(1e-14);
2464   if(getNumberOfTuples()!=3 && getNumberOfComponents()!=2)
2465     throw INTERP_KERNEL::Exception("DataArrayDouble::asArcCircle : this method expects");
2466   const double *pt(begin());
2467   MCAuto<INTERP_KERNEL::Node> n0(new INTERP_KERNEL::Node(pt[0],pt[1])),n1(new INTERP_KERNEL::Node(pt[2],pt[3])),n2(new INTERP_KERNEL::Node(pt[4],pt[5]));
2468   {
2469     INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeLin> e1(new INTERP_KERNEL::EdgeLin(n0,n2)),e2(new INTERP_KERNEL::EdgeLin(n2,n1));
2470     INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
2471     bool colinearity(inters.areColinears());
2472     if(colinearity)
2473       throw INTERP_KERNEL::Exception("DataArrayDouble::asArcOfCircle : 3 points in this have been detected as colinear !");
2474   }
2475   INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeArcCircle> ret(new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1));
2476   const double *c(ret->getCenter());
2477   center[0]=c[0]; center[1]=c[1];
2478   radius=ret->getRadius();
2479   ang=ret->getAngle();
2480 }
2481
2482 /*!
2483  * Sorts value within every tuple of \a this array.
2484  *  \param [in] asc - if \a true, the values are sorted in ascending order, else,
2485  *              in descending order.
2486  *  \throw If \a this is not allocated.
2487  */
2488 void DataArrayDouble::sortPerTuple(bool asc)
2489 {
2490   checkAllocated();
2491   double *pt=getPointer();
2492   mcIdType nbOfTuple(getNumberOfTuples());
2493   std::size_t nbOfComp(getNumberOfComponents());
2494   if(asc)
2495     for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2496       std::sort(pt,pt+nbOfComp);
2497   else
2498     for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2499       std::sort(pt,pt+nbOfComp,std::greater<double>());
2500   declareAsNew();
2501 }
2502
2503 /*!
2504  * Modify all elements of \a this array, so that
2505  * an element _x_ becomes \f$ numerator / x \f$.
2506  *  \warning If an exception is thrown because of presence of 0.0 element in \a this
2507  *           array, all elements processed before detection of the zero element remain
2508  *           modified.
2509  *  \param [in] numerator - the numerator used to modify array elements.
2510  *  \throw If \a this is not allocated.
2511  *  \throw If there is an element equal to 0.0 in \a this array.
2512  */
2513 void DataArrayDouble::applyInv(double numerator)
2514 {
2515   checkAllocated();
2516   double *ptr=getPointer();
2517   std::size_t nbOfElems=getNbOfElems();
2518   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2519     {
2520       if(std::abs(*ptr)>std::numeric_limits<double>::min())
2521         {
2522           *ptr=numerator/(*ptr);
2523         }
2524       else
2525         {
2526           std::ostringstream oss; oss << "DataArrayDouble::applyInv : presence of null value in tuple #" << i/getNumberOfComponents() << " component #" << i%getNumberOfComponents();
2527           oss << " !";
2528           throw INTERP_KERNEL::Exception(oss.str().c_str());
2529         }
2530     }
2531   declareAsNew();
2532 }
2533
2534 /*!
2535  * Modify all elements of \a this array, so that
2536  * an element _x_ becomes <em> val ^ x </em>. Contrary to DataArrayInt::applyPow
2537  * all values in \a this have to be >= 0 if val is \b not integer.
2538  *  \param [in] val - the value used to apply pow on all array elements.
2539  *  \throw If \a this is not allocated.
2540  *  \warning If an exception is thrown because of presence of 0 element in \a this
2541  *           array and \a val is \b not integer, all elements processed before detection of the zero element remain
2542  *           modified.
2543  */
2544 void DataArrayDouble::applyPow(double val)
2545 {
2546   checkAllocated();
2547   double *ptr=getPointer();
2548   std::size_t nbOfElems=getNbOfElems();
2549   int val2=(int)val;
2550   bool isInt=((double)val2)==val;
2551   if(!isInt)
2552     {
2553       for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2554         {
2555           if(*ptr>=0)
2556             *ptr=pow(*ptr,val);
2557           else
2558             {
2559               std::ostringstream oss; oss << "DataArrayDouble::applyPow (double) : At elem # " << i << " value is " << *ptr << " ! must be >=0. !";
2560               throw INTERP_KERNEL::Exception(oss.str().c_str());
2561             }
2562         }
2563     }
2564   else
2565     {
2566       for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2567         *ptr=pow(*ptr,val2);
2568     }
2569   declareAsNew();
2570 }
2571
2572 /*!
2573  * Modify all elements of \a this array, so that
2574  * an element _x_ becomes \f$ val ^ x \f$.
2575  *  \param [in] val - the value used to apply pow on all array elements.
2576  *  \throw If \a this is not allocated.
2577  *  \throw If \a val < 0.
2578  *  \warning If an exception is thrown because of presence of 0 element in \a this
2579  *           array, all elements processed before detection of the zero element remain
2580  *           modified.
2581  */
2582 void DataArrayDouble::applyRPow(double val)
2583 {
2584   checkAllocated();
2585   if(val<0.)
2586     throw INTERP_KERNEL::Exception("DataArrayDouble::applyRPow : the input value has to be >= 0 !");
2587   double *ptr=getPointer();
2588   std::size_t nbOfElems=getNbOfElems();
2589   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2590     *ptr=pow(val,*ptr);
2591   declareAsNew();
2592 }
2593
2594 /*!
2595  * Returns a new DataArrayDouble created from \a this one by applying \a
2596  * FunctionToEvaluate to every tuple of \a this array. Textual data is not copied.
2597  * For more info see \ref MEDCouplingArrayApplyFunc
2598  *  \param [in] nbOfComp - number of components in the result array.
2599  *  \param [in] func - the \a FunctionToEvaluate declared as
2600  *              \c bool (*\a func)(\c const \c double *\a pos, \c double *\a res),
2601  *              where \a pos points to the first component of a tuple of \a this array
2602  *              and \a res points to the first component of a tuple of the result array.
2603  *              Note that length (number of components) of \a pos can differ from
2604  *              that of \a res.
2605  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2606  *          same number of tuples as \a this array.
2607  *          The caller is to delete this result array using decrRef() as it is no more
2608  *          needed.
2609  *  \throw If \a this is not allocated.
2610  *  \throw If \a func returns \a false.
2611  */
2612 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, FunctionToEvaluate func) const
2613 {
2614   checkAllocated();
2615   DataArrayDouble *newArr=DataArrayDouble::New();
2616   mcIdType nbOfTuples(getNumberOfTuples());
2617   std::size_t oldNbOfComp(getNumberOfComponents());
2618   newArr->alloc(nbOfTuples,nbOfComp);
2619   const double *ptr=getConstPointer();
2620   double *ptrToFill=newArr->getPointer();
2621   for(mcIdType i=0;i<nbOfTuples;i++)
2622     {
2623       if(!func(ptr+i*oldNbOfComp,ptrToFill+i*nbOfComp))
2624         {
2625           std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2626           std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2627           oss << ") : Evaluation of function failed !";
2628           newArr->decrRef();
2629           throw INTERP_KERNEL::Exception(oss.str().c_str());
2630         }
2631     }
2632   return newArr;
2633 }
2634
2635 /*!
2636  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2637  * tuple of \a this array. Textual data is not copied.
2638  * For more info see \ref MEDCouplingArrayApplyFunc1.
2639  *  \param [in] nbOfComp - number of components in the result array.
2640  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2641  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2642  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2643  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2644  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2645  *          same number of tuples as \a this array and \a nbOfComp components.
2646  *          The caller is to delete this result array using decrRef() as it is no more
2647  *          needed.
2648  *  \throw If \a this is not allocated.
2649  *  \throw If computing \a func fails.
2650  */
2651 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2652 {
2653   INTERP_KERNEL::ExprParser expr(func);
2654   expr.parse();
2655   std::set<std::string> vars;
2656   expr.getTrueSetOfVars(vars);
2657   std::vector<std::string> varsV(vars.begin(),vars.end());
2658   return applyFuncNamedCompo(nbOfComp,varsV,func,isSafe);
2659 }
2660
2661 /*!
2662  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2663  * tuple of \a this array. Textual data is not copied. This method works by tuples (whatever its size).
2664  * If \a this is a one component array, call applyFuncOnThis instead that performs the same work faster.
2665  *
2666  * For more info see \ref MEDCouplingArrayApplyFunc0.
2667  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2668  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2669  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2670  *                       If false the computation is carried on without any notification. When false the evaluation is a little faster.
2671  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2672  *          same number of tuples and components as \a this array.
2673  *          The caller is to delete this result array using decrRef() as it is no more
2674  *          needed.
2675  *  \sa applyFuncOnThis
2676  *  \throw If \a this is not allocated.
2677  *  \throw If computing \a func fails.
2678  */
2679 DataArrayDouble *DataArrayDouble::applyFunc(const std::string& func, bool isSafe) const
2680 {
2681   std::size_t nbOfComp(getNumberOfComponents());
2682   if(nbOfComp<=0)
2683     throw INTERP_KERNEL::Exception("DataArrayDouble::applyFunc : output number of component must be > 0 !");
2684   checkAllocated();
2685   mcIdType nbOfTuples(getNumberOfTuples());
2686   MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2687   newArr->alloc(nbOfTuples,nbOfComp);
2688   INTERP_KERNEL::ExprParser expr(func);
2689   expr.parse();
2690   std::set<std::string> vars;
2691   expr.getTrueSetOfVars(vars);
2692   if(vars.size()>1)
2693     {
2694       std::ostringstream oss; oss << "DataArrayDouble::applyFunc : this method works only with at most one var func expression ! If you need to map comps on variables please use applyFuncCompo or applyFuncNamedCompo instead ! Vars in expr are : ";
2695       std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2696       throw INTERP_KERNEL::Exception(oss.str().c_str());
2697     }
2698   if(vars.empty())
2699     {
2700       expr.prepareFastEvaluator();
2701       newArr->rearrange(1);
2702       newArr->fillWithValue(expr.evaluateDouble());
2703       newArr->rearrange(nbOfComp);
2704       return newArr.retn();
2705     }
2706   std::vector<std::string> vars2(vars.begin(),vars.end());
2707   double buff,*ptrToFill(newArr->getPointer());
2708   const double *ptr(begin());
2709   std::vector<double> stck;
2710   expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2711   expr.prepareFastEvaluator();
2712   if(!isSafe)
2713     {
2714       for(mcIdType i=0;i<nbOfTuples;i++)
2715         {
2716           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2717             {
2718               buff=*ptr;
2719               expr.evaluateDoubleInternal(stck);
2720               *ptrToFill=stck.back();
2721               stck.pop_back();
2722             }
2723         }
2724     }
2725   else
2726     {
2727       for(mcIdType i=0;i<nbOfTuples;i++)
2728         {
2729           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2730             {
2731               buff=*ptr;
2732               try
2733               {
2734                   expr.evaluateDoubleInternalSafe(stck);
2735               }
2736               catch(INTERP_KERNEL::Exception& e)
2737               {
2738                   std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2739                   oss << buff;
2740                   oss << ") : Evaluation of function failed !" << e.what();
2741                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2742               }
2743               *ptrToFill=stck.back();
2744               stck.pop_back();
2745             }
2746         }
2747     }
2748   return newArr.retn();
2749 }
2750
2751 /*!
2752  * This method is a non const method that modify the array in \a this.
2753  * This method only works on one component array. It means that function \a func must
2754  * contain at most one variable.
2755  * This method is a specialization of applyFunc method with one parameter on one component array.
2756  *
2757  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2758  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2759  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2760  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2761  *
2762  * \sa applyFunc
2763  */
2764 void DataArrayDouble::applyFuncOnThis(const std::string& func, bool isSafe)
2765 {
2766   std::size_t nbOfComp(getNumberOfComponents());
2767   if(nbOfComp<=0)
2768     throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncOnThis : output number of component must be > 0 !");
2769   checkAllocated();
2770   mcIdType nbOfTuples(getNumberOfTuples());
2771   INTERP_KERNEL::ExprParser expr(func);
2772   expr.parse();
2773   std::set<std::string> vars;
2774   expr.getTrueSetOfVars(vars);
2775   if(vars.size()>1)
2776     {
2777       std::ostringstream oss; oss << "DataArrayDouble::applyFuncOnThis : this method works only with at most one var func expression ! If you need to map comps on variables please use applyFuncCompo or applyFuncNamedCompo instead ! Vars in expr are : ";
2778       std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2779       throw INTERP_KERNEL::Exception(oss.str().c_str());
2780     }
2781   if(vars.empty())
2782     {
2783       expr.prepareFastEvaluator();
2784       std::vector<std::string> compInfo(getInfoOnComponents());
2785       rearrange(1);
2786       fillWithValue(expr.evaluateDouble());
2787       rearrange(nbOfComp);
2788       setInfoOnComponents(compInfo);
2789       return ;
2790     }
2791   std::vector<std::string> vars2(vars.begin(),vars.end());
2792   double buff,*ptrToFill(getPointer());
2793   const double *ptr(begin());
2794   std::vector<double> stck;
2795   expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2796   expr.prepareFastEvaluator();
2797   if(!isSafe)
2798     {
2799       for(mcIdType i=0;i<nbOfTuples;i++)
2800         {
2801           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2802             {
2803               buff=*ptr;
2804               expr.evaluateDoubleInternal(stck);
2805               *ptrToFill=stck.back();
2806               stck.pop_back();
2807             }
2808         }
2809     }
2810   else
2811     {
2812       for(mcIdType i=0;i<nbOfTuples;i++)
2813         {
2814           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2815             {
2816               buff=*ptr;
2817               try
2818               {
2819                   expr.evaluateDoubleInternalSafe(stck);
2820               }
2821               catch(INTERP_KERNEL::Exception& e)
2822               {
2823                   std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2824                   oss << buff;
2825                   oss << ") : Evaluation of function failed !" << e.what();
2826                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2827               }
2828               *ptrToFill=stck.back();
2829               stck.pop_back();
2830             }
2831         }
2832     }
2833 }
2834
2835 /*!
2836  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2837  * tuple of \a this array. Textual data is not copied.
2838  * For more info see \ref MEDCouplingArrayApplyFunc2.
2839  *  \param [in] nbOfComp - number of components in the result array.
2840  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2841  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2842  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2843  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2844  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2845  *          same number of tuples as \a this array.
2846  *          The caller is to delete this result array using decrRef() as it is no more
2847  *          needed.
2848  *  \throw If \a this is not allocated.
2849  *  \throw If \a func contains vars that are not in \a this->getInfoOnComponent().
2850  *  \throw If computing \a func fails.
2851  */
2852 DataArrayDouble *DataArrayDouble::applyFuncCompo(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2853 {
2854   return applyFuncNamedCompo(nbOfComp,getVarsOnComponent(),func,isSafe);
2855 }
2856
2857 /*!
2858  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2859  * tuple of \a this array. Textual data is not copied.
2860  * For more info see \ref MEDCouplingArrayApplyFunc3.
2861  *  \param [in] nbOfComp - number of components in the result array.
2862  *  \param [in] varsOrder - sequence of vars defining their order.
2863  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2864  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2865  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2866  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2867  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2868  *          same number of tuples as \a this array.
2869  *          The caller is to delete this result array using decrRef() as it is no more
2870  *          needed.
2871  *  \throw If \a this is not allocated.
2872  *  \throw If \a func contains vars not in \a varsOrder.
2873  *  \throw If computing \a func fails.
2874  */
2875 DataArrayDouble *DataArrayDouble::applyFuncNamedCompo(std::size_t nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func, bool isSafe) const
2876 {
2877   if(nbOfComp<=0)
2878     throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncNamedCompo : output number of component must be > 0 !");
2879   std::vector<std::string> varsOrder2(varsOrder);
2880   std::size_t oldNbOfComp(getNumberOfComponents());
2881   for(std::size_t i=varsOrder.size();i<oldNbOfComp;i++)
2882     varsOrder2.push_back(std::string());
2883   checkAllocated();
2884   mcIdType nbOfTuples(getNumberOfTuples());
2885   INTERP_KERNEL::ExprParser expr(func);
2886   expr.parse();
2887   std::set<std::string> vars;
2888   expr.getTrueSetOfVars(vars);
2889   if(vars.size()>oldNbOfComp)
2890     {
2891       std::ostringstream oss; oss << "The field has " << oldNbOfComp << " components and there are ";
2892       oss << vars.size() << " variables : ";
2893       std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2894       throw INTERP_KERNEL::Exception(oss.str().c_str());
2895     }
2896   MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2897   newArr->alloc(nbOfTuples,nbOfComp);
2898   INTERP_KERNEL::AutoPtr<double> buff(new double[oldNbOfComp]);
2899   double *buffPtr(buff),*ptrToFill;
2900   std::vector<double> stck;
2901   for(std::size_t iComp=0;iComp<nbOfComp;iComp++)
2902     {
2903       expr.prepareExprEvaluationDouble(varsOrder2,(int)oldNbOfComp,(int)nbOfComp,(int)iComp,buffPtr,buffPtr+oldNbOfComp);
2904       expr.prepareFastEvaluator();
2905       const double *ptr(getConstPointer());
2906       ptrToFill=newArr->getPointer()+iComp;
2907       if(!isSafe)
2908         {
2909           for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2910             {
2911               std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2912               expr.evaluateDoubleInternal(stck);
2913               *ptrToFill=stck.back();
2914               stck.pop_back();
2915             }
2916         }
2917       else
2918         {
2919           for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2920             {
2921               std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2922               try
2923               {
2924                   expr.evaluateDoubleInternalSafe(stck);
2925                   *ptrToFill=stck.back();
2926                   stck.pop_back();
2927               }
2928               catch(INTERP_KERNEL::Exception& e)
2929               {
2930                   std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2931                   std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2932                   oss << ") : Evaluation of function failed !" << e.what();
2933                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2934               }
2935             }
2936         }
2937     }
2938   return newArr.retn();
2939 }
2940
2941 void DataArrayDouble::applyFuncFast32(const std::string& func)
2942 {
2943   checkAllocated();
2944   INTERP_KERNEL::ExprParser expr(func);
2945   expr.parse();
2946   char *funcStr=expr.compileX86();
2947   MYFUNCPTR funcPtr;
2948   *((void **)&funcPtr)=funcStr;//he he...
2949   //
2950   double *ptr=getPointer();
2951   std::size_t nbOfComp=getNumberOfComponents();
2952   mcIdType nbOfTuples=getNumberOfTuples();
2953   std::size_t nbOfElems=nbOfTuples*nbOfComp;
2954   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2955     *ptr=funcPtr(*ptr);
2956   declareAsNew();
2957 }
2958
2959 void DataArrayDouble::applyFuncFast64(const std::string& func)
2960 {
2961   checkAllocated();
2962   INTERP_KERNEL::ExprParser expr(func);
2963   expr.parse();
2964   char *funcStr=expr.compileX86_64();
2965   MYFUNCPTR funcPtr;
2966   *((void **)&funcPtr)=funcStr;//he he...
2967   //
2968   double *ptr=getPointer();
2969   std::size_t nbOfComp=getNumberOfComponents();
2970   mcIdType nbOfTuples=getNumberOfTuples();
2971   std::size_t nbOfElems=nbOfTuples*nbOfComp;
2972   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2973     *ptr=funcPtr(*ptr);
2974   declareAsNew();
2975 }
2976
2977 /*!
2978  * \return a new object that is the result of the symmetry along 3D plane defined by its normal vector \a normalVector and a point \a point.
2979  */
2980 MCAuto<DataArrayDouble> DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const
2981 {
2982   checkAllocated();
2983   if(getNumberOfComponents()!=3)
2984     throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !");
2985   mcIdType nbTuples(getNumberOfTuples());
2986   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
2987   ret->alloc(nbTuples,3);
2988   Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer());
2989   return ret;
2990 }
2991
2992 DataArrayDoubleIterator *DataArrayDouble::iterator()
2993 {
2994   return new DataArrayDoubleIterator(this);
2995 }
2996
2997 /*!
2998  * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
2999  * array whose values are within a given range. Textual data is not copied.
3000  *  \param [in] vmin - a lowest acceptable value (included).
3001  *  \param [in] vmax - a greatest acceptable value (included).
3002  *  \return DataArrayInt * - the new instance of DataArrayInt.
3003  *          The caller is to delete this result array using decrRef() as it is no more
3004  *          needed.
3005  *  \throw If \a this->getNumberOfComponents() != 1.
3006  *
3007  *  \sa DataArrayDouble::findIdsNotInRange
3008  *
3009  *  \if ENABLE_EXAMPLES
3010  *  \ref cpp_mcdataarraydouble_getidsinrange "Here is a C++ example".<br>
3011  *  \ref py_mcdataarraydouble_getidsinrange "Here is a Python example".
3012  *  \endif
3013  */
3014 DataArrayIdType *DataArrayDouble::findIdsInRange(double vmin, double vmax) const
3015 {
3016   checkAllocated();
3017   if(getNumberOfComponents()!=1)
3018     throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsInRange : this must have exactly one component !");
3019   const double *cptr(begin());
3020   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3021   mcIdType nbOfTuples(getNumberOfTuples());
3022   for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3023     if(*cptr>=vmin && *cptr<=vmax)
3024       ret->pushBackSilent(i);
3025   return ret.retn();
3026 }
3027
3028 /*!
3029  * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3030  * array whose values are not within a given range. Textual data is not copied.
3031  *  \param [in] vmin - a lowest not acceptable value (excluded).
3032  *  \param [in] vmax - a greatest not acceptable value (excluded).
3033  *  \return DataArrayInt * - the new instance of DataArrayInt.
3034  *          The caller is to delete this result array using decrRef() as it is no more
3035  *          needed.
3036  *  \throw If \a this->getNumberOfComponents() != 1.
3037  *
3038  *  \sa DataArrayDouble::findIdsInRange
3039  */
3040 DataArrayIdType *DataArrayDouble::findIdsNotInRange(double vmin, double vmax) const
3041 {
3042   checkAllocated();
3043   if(getNumberOfComponents()!=1)
3044     throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsNotInRange : this must have exactly one component !");
3045   const double *cptr(begin());
3046   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3047   mcIdType nbOfTuples(getNumberOfTuples());
3048   for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3049     if(*cptr<vmin || *cptr>vmax)
3050       ret->pushBackSilent(i);
3051   return ret.retn();
3052 }
3053
3054 /*!
3055  * Returns a new DataArrayDouble by concatenating two given arrays, so that (1) the number
3056  * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3057  * the number of component in the result array is same as that of each of given arrays.
3058  * Info on components is copied from the first of the given arrays. Number of components
3059  * in the given arrays must be  the same.
3060  *  \param [in] a1 - an array to include in the result array.
3061  *  \param [in] a2 - another array to include in the result array.
3062  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3063  *          The caller is to delete this result array using decrRef() as it is no more
3064  *          needed.
3065  *  \throw If both \a a1 and \a a2 are NULL.
3066  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents().
3067  */
3068 DataArrayDouble *DataArrayDouble::Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2)
3069 {
3070   std::vector<const DataArrayDouble *> tmp(2);
3071   tmp[0]=a1; tmp[1]=a2;
3072   return Aggregate(tmp);
3073 }
3074
3075 /*!
3076  * Returns a new DataArrayDouble by concatenating all given arrays, so that (1) the number
3077  * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3078  * the number of component in the result array is same as that of each of given arrays.
3079  * Info on components is copied from the first of the given arrays. Number of components
3080  * in the given arrays must be  the same.
3081  * If the number of non null of elements in \a arr is equal to one the returned object is a copy of it
3082  * not the object itself.
3083  *  \param [in] arr - a sequence of arrays to include in the result array.
3084  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3085  *          The caller is to delete this result array using decrRef() as it is no more
3086  *          needed.
3087  *  \throw If all arrays within \a arr are NULL.
3088  *  \throw If getNumberOfComponents() of arrays within \a arr.
3089  */
3090 DataArrayDouble *DataArrayDouble::Aggregate(const std::vector<const DataArrayDouble *>& arr)
3091 {
3092   std::vector<const DataArrayDouble *> a;
3093   for(std::vector<const DataArrayDouble *>::const_iterator it4=arr.begin();it4!=arr.end();it4++)
3094     if(*it4)
3095       a.push_back(*it4);
3096   if(a.empty())
3097     throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : input list must contain at least one NON EMPTY DataArrayDouble !");
3098   std::vector<const DataArrayDouble *>::const_iterator it=a.begin();
3099   std::size_t nbOfComp((*it)->getNumberOfComponents());
3100   mcIdType nbt=(*it++)->getNumberOfTuples();
3101   for(mcIdType i=1;it!=a.end();it++,i++)
3102     {
3103       if((*it)->getNumberOfComponents()!=nbOfComp)
3104         throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : Nb of components mismatch for array aggregation !");
3105       nbt+=(*it)->getNumberOfTuples();
3106     }
3107   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
3108   ret->alloc(nbt,nbOfComp);
3109   double *pt=ret->getPointer();
3110   for(it=a.begin();it!=a.end();it++)
3111     pt=std::copy((*it)->getConstPointer(),(*it)->getConstPointer()+(*it)->getNbOfElems(),pt);
3112   ret->copyStringInfoFrom(*(a[0]));
3113   return ret.retn();
3114 }
3115
3116 /*!
3117  * Returns a new DataArrayDouble containing a dot product of two given arrays, so that
3118  * the i-th tuple of the result array is a sum of products of j-th components of i-th
3119  * tuples of given arrays (\f$ a_i = \sum_{j=1}^n a1_j * a2_j \f$).
3120  * Info on components and name is copied from the first of the given arrays.
3121  * Number of tuples and components in the given arrays must be the same.
3122  *  \param [in] a1 - a given array.
3123  *  \param [in] a2 - another given array.
3124  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3125  *          The caller is to delete this result array using decrRef() as it is no more
3126  *          needed.
3127  *  \throw If either \a a1 or \a a2 is NULL.
3128  *  \throw If any given array is not allocated.
3129  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3130  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3131  */
3132 DataArrayDouble *DataArrayDouble::Dot(const DataArrayDouble *a1, const DataArrayDouble *a2)
3133 {
3134   if(!a1 || !a2)
3135     throw INTERP_KERNEL::Exception("DataArrayDouble::Dot : input DataArrayDouble instance is NULL !");
3136   a1->checkAllocated();
3137   a2->checkAllocated();
3138   std::size_t nbOfComp(a1->getNumberOfComponents());
3139   if(nbOfComp!=a2->getNumberOfComponents())
3140     throw INTERP_KERNEL::Exception("Nb of components mismatch for array Dot !");
3141   mcIdType nbOfTuple(a1->getNumberOfTuples());
3142   if(nbOfTuple!=a2->getNumberOfTuples())
3143     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Dot !");
3144   DataArrayDouble *ret=DataArrayDouble::New();
3145   ret->alloc(nbOfTuple,1);
3146   double *retPtr=ret->getPointer();
3147   const double *a1Ptr=a1->begin(),*a2Ptr(a2->begin());
3148   for(mcIdType i=0;i<nbOfTuple;i++)
3149     {
3150       double sum=0.;
3151       for(std::size_t j=0;j<nbOfComp;j++)
3152         sum+=a1Ptr[i*nbOfComp+j]*a2Ptr[i*nbOfComp+j];
3153       retPtr[i]=sum;
3154     }
3155   ret->setInfoOnComponent(0,a1->getInfoOnComponent(0));
3156   ret->setName(a1->getName());
3157   return ret;
3158 }
3159
3160 /*!
3161  * Returns a new DataArrayDouble containing a cross product of two given arrays, so that
3162  * the i-th tuple of the result array contains 3 components of a vector which is a cross
3163  * product of two vectors defined by the i-th tuples of given arrays.
3164  * Info on components is copied from the first of the given arrays.
3165  * Number of tuples in the given arrays must be the same.
3166  * Number of components in the given arrays must be 3.
3167  *  \param [in] a1 - a given array.
3168  *  \param [in] a2 - another given array.
3169  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3170  *          The caller is to delete this result array using decrRef() as it is no more
3171  *          needed.
3172  *  \throw If either \a a1 or \a a2 is NULL.
3173  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3174  *  \throw If \a a1->getNumberOfComponents() != 3
3175  *  \throw If \a a2->getNumberOfComponents() != 3
3176  */
3177 DataArrayDouble *DataArrayDouble::CrossProduct(const DataArrayDouble *a1, const DataArrayDouble *a2)
3178 {
3179   if(!a1 || !a2)
3180     throw INTERP_KERNEL::Exception("DataArrayDouble::CrossProduct : input DataArrayDouble instance is NULL !");
3181   std::size_t nbOfComp(a1->getNumberOfComponents());
3182   if(nbOfComp!=a2->getNumberOfComponents())
3183     throw INTERP_KERNEL::Exception("Nb of components mismatch for array crossProduct !");
3184   if(nbOfComp!=3)
3185     throw INTERP_KERNEL::Exception("Nb of components must be equal to 3 for array crossProduct !");
3186   mcIdType nbOfTuple(a1->getNumberOfTuples());
3187   if(nbOfTuple!=a2->getNumberOfTuples())
3188     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array crossProduct !");
3189   DataArrayDouble *ret=DataArrayDouble::New();
3190   ret->alloc(nbOfTuple,3);
3191   double *retPtr=ret->getPointer();
3192   const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3193   for(mcIdType i=0;i<nbOfTuple;i++)
3194     {
3195       retPtr[3*i]=a1Ptr[3*i+1]*a2Ptr[3*i+2]-a1Ptr[3*i+2]*a2Ptr[3*i+1];
3196       retPtr[3*i+1]=a1Ptr[3*i+2]*a2Ptr[3*i]-a1Ptr[3*i]*a2Ptr[3*i+2];
3197       retPtr[3*i+2]=a1Ptr[3*i]*a2Ptr[3*i+1]-a1Ptr[3*i+1]*a2Ptr[3*i];
3198     }
3199   ret->copyStringInfoFrom(*a1);
3200   return ret;
3201 }
3202
3203 /*!
3204  * Returns a new DataArrayDouble containing maximal values of two given arrays.
3205  * Info on components is copied from the first of the given arrays.
3206  * Number of tuples and components in the given arrays must be the same.
3207  *  \param [in] a1 - an array to compare values with another one.
3208  *  \param [in] a2 - another array to compare values with the first one.
3209  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3210  *          The caller is to delete this result array using decrRef() as it is no more
3211  *          needed.
3212  *  \throw If either \a a1 or \a a2 is NULL.
3213  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3214  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3215  */
3216 DataArrayDouble *DataArrayDouble::Max(const DataArrayDouble *a1, const DataArrayDouble *a2)
3217 {
3218   if(!a1 || !a2)
3219     throw INTERP_KERNEL::Exception("DataArrayDouble::Max : input DataArrayDouble instance is NULL !");
3220   std::size_t nbOfComp(a1->getNumberOfComponents());
3221   if(nbOfComp!=a2->getNumberOfComponents())
3222     throw INTERP_KERNEL::Exception("Nb of components mismatch for array Max !");
3223   mcIdType nbOfTuple(a1->getNumberOfTuples());
3224   if(nbOfTuple!=a2->getNumberOfTuples())
3225     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Max !");
3226   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3227   ret->alloc(nbOfTuple,nbOfComp);
3228   double *retPtr(ret->getPointer());
3229   const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3230   std::size_t nbElem(nbOfTuple*nbOfComp);
3231   for(std::size_t i=0;i<nbElem;i++)
3232     retPtr[i]=std::max(a1Ptr[i],a2Ptr[i]);
3233   ret->copyStringInfoFrom(*a1);
3234   return ret.retn();
3235 }
3236
3237 /*!
3238  * Returns a new DataArrayDouble containing minimal values of two given arrays.
3239  * Info on components is copied from the first of the given arrays.
3240  * Number of tuples and components in the given arrays must be the same.
3241  *  \param [in] a1 - an array to compare values with another one.
3242  *  \param [in] a2 - another array to compare values with the first one.
3243  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3244  *          The caller is to delete this result array using decrRef() as it is no more
3245  *          needed.
3246  *  \throw If either \a a1 or \a a2 is NULL.
3247  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3248  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3249  */
3250 DataArrayDouble *DataArrayDouble::Min(const DataArrayDouble *a1, const DataArrayDouble *a2)
3251 {
3252   if(!a1 || !a2)
3253     throw INTERP_KERNEL::Exception("DataArrayDouble::Min : input DataArrayDouble instance is NULL !");
3254   std::size_t nbOfComp(a1->getNumberOfComponents());
3255   if(nbOfComp!=a2->getNumberOfComponents())
3256     throw INTERP_KERNEL::Exception("Nb of components mismatch for array min !");
3257   mcIdType nbOfTuple(a1->getNumberOfTuples());
3258   if(nbOfTuple!=a2->getNumberOfTuples())
3259     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array min !");
3260   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3261   ret->alloc(nbOfTuple,nbOfComp);
3262   double *retPtr(ret->getPointer());
3263   const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3264   std::size_t nbElem(nbOfTuple*nbOfComp);
3265   for(std::size_t i=0;i<nbElem;i++)
3266     retPtr[i]=std::min(a1Ptr[i],a2Ptr[i]);
3267   ret->copyStringInfoFrom(*a1);
3268   return ret.retn();
3269 }
3270
3271 /*!
3272  * Returns a new DataArrayDouble that is the result of pow of two given arrays. There are 3
3273  * valid cases.
3274  *
3275  *  \param [in] a1 - an array to pow up.
3276  *  \param [in] a2 - another array to sum up.
3277  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3278  *          The caller is to delete this result array using decrRef() as it is no more
3279  *          needed.
3280  *  \throw If either \a a1 or \a a2 is NULL.
3281  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3282  *  \throw If \a a1->getNumberOfComponents() != 1 or \a a2->getNumberOfComponents() != 1.
3283  *  \throw If there is a negative value in \a a1.
3284  */
3285 DataArrayDouble *DataArrayDouble::Pow(const DataArrayDouble *a1, const DataArrayDouble *a2)
3286 {
3287   if(!a1 || !a2)
3288     throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : at least one of input instances is null !");
3289   mcIdType nbOfTuple=a1->getNumberOfTuples();
3290   mcIdType nbOfTuple2=a2->getNumberOfTuples();
3291   std::size_t nbOfComp=a1->getNumberOfComponents();
3292   std::size_t nbOfComp2=a2->getNumberOfComponents();
3293   if(nbOfTuple!=nbOfTuple2)
3294     throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of tuples mismatches !");
3295   if(nbOfComp!=1 || nbOfComp2!=1)
3296     throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of components of both arrays must be equal to 1 !");
3297   MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(nbOfTuple,1);
3298   const double *ptr1(a1->begin()),*ptr2(a2->begin());
3299   double *ptr=ret->getPointer();
3300   for(mcIdType i=0;i<nbOfTuple;i++,ptr1++,ptr2++,ptr++)
3301     {
3302       if(*ptr1>=0)
3303         {
3304           *ptr=pow(*ptr1,*ptr2);
3305         }
3306       else
3307         {
3308           std::ostringstream oss; oss << "DataArrayDouble::Pow : on tuple #" << i << " of a1 value is < 0 (" << *ptr1 << ") !";
3309           throw INTERP_KERNEL::Exception(oss.str().c_str());
3310         }
3311     }
3312   return ret.retn();
3313 }
3314
3315 /*!
3316  * Apply pow on values of another DataArrayDouble to values of \a this one.
3317  *
3318  *  \param [in] other - an array to pow to \a this one.
3319  *  \throw If \a other is NULL.
3320  *  \throw If \a this->getNumberOfTuples() != \a other->getNumberOfTuples()
3321  *  \throw If \a this->getNumberOfComponents() != 1 or \a other->getNumberOfComponents() != 1
3322  *  \throw If there is a negative value in \a this.
3323  */
3324 void DataArrayDouble::powEqual(const DataArrayDouble *other)
3325 {
3326   if(!other)
3327     throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : input instance is null !");
3328   mcIdType nbOfTuple=getNumberOfTuples();
3329   mcIdType nbOfTuple2=other->getNumberOfTuples();
3330   std::size_t nbOfComp=getNumberOfComponents();
3331   std::size_t nbOfComp2=other->getNumberOfComponents();
3332   if(nbOfTuple!=nbOfTuple2)
3333     throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of tuples mismatches !");
3334   if(nbOfComp!=1 || nbOfComp2!=1)
3335     throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of components of both arrays must be equal to 1 !");
3336   double *ptr=getPointer();
3337   const double *ptrc=other->begin();
3338   for(mcIdType i=0;i<nbOfTuple;i++,ptrc++,ptr++)
3339     {
3340       if(*ptr>=0)
3341         *ptr=pow(*ptr,*ptrc);
3342       else
3343         {
3344           std::ostringstream oss; oss << "DataArrayDouble::powEqual : on tuple #" << i << " of this value is < 0 (" << *ptr << ") !";
3345           throw INTERP_KERNEL::Exception(oss.str().c_str());
3346         }
3347     }
3348   declareAsNew();
3349 }
3350
3351 /*!
3352  * This method is \b NOT wrapped into python because it can be useful only for performance reasons in C++ context.
3353  * All values in \a this must be 0. or 1. within eps error. 0 means false, 1 means true.
3354  * If an another value than 0 or 1 appear (within eps precision) an INTERP_KERNEL::Exception will be thrown.
3355  *
3356  * \throw if \a this is not allocated.
3357  * \throw if \a this has not exactly one component.
3358  */
3359 std::vector<bool> DataArrayDouble::toVectorOfBool(double eps) const
3360 {
3361   checkAllocated();
3362   if(getNumberOfComponents()!=1)
3363     throw INTERP_KERNEL::Exception("DataArrayDouble::toVectorOfBool : must be applied on single component array !");
3364   mcIdType nbt(getNumberOfTuples());
3365   std::vector<bool> ret(nbt);
3366   const double *pt(begin());
3367   for(mcIdType i=0;i<nbt;i++)
3368     {
3369       if(fabs(pt[i])<eps)
3370         ret[i]=false;
3371       else if(fabs(pt[i]-1.)<eps)
3372         ret[i]=true;
3373       else
3374         {
3375           std::ostringstream oss; oss << "DataArrayDouble::toVectorOfBool : the tuple #" << i << " has value " << pt[i] << " is invalid ! must be 0. or 1. !";
3376           throw INTERP_KERNEL::Exception(oss.str().c_str());
3377         }
3378     }
3379   return ret;
3380 }
3381
3382 /*!
3383  * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3384  * Server side.
3385  */
3386 void DataArrayDouble::getTinySerializationIntInformation(std::vector<mcIdType>& tinyInfo) const
3387 {
3388   tinyInfo.resize(2);
3389   if(isAllocated())
3390     {
3391       tinyInfo[0]=getNumberOfTuples();
3392       tinyInfo[1]=ToIdType(getNumberOfComponents());
3393     }
3394   else
3395     {
3396       tinyInfo[0]=-1;
3397       tinyInfo[1]=-1;
3398     }
3399 }
3400
3401 /*!
3402  * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3403  * Server side.
3404  */
3405 void DataArrayDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
3406 {
3407   if(isAllocated())
3408     {
3409       std::size_t nbOfCompo(getNumberOfComponents());
3410       tinyInfo.resize(nbOfCompo+1);
3411       tinyInfo[0]=getName();
3412       for(std::size_t i=0;i<nbOfCompo;i++)
3413         tinyInfo[i+1]=getInfoOnComponent(i);
3414     }
3415   else
3416     {
3417       tinyInfo.resize(1);
3418       tinyInfo[0]=getName();
3419     }
3420 }
3421
3422 /*!
3423  * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3424  * This method returns if a feeding is needed.
3425  */
3426 bool DataArrayDouble::resizeForUnserialization(const std::vector<mcIdType>& tinyInfoI)
3427 {
3428   mcIdType nbOfTuple=tinyInfoI[0];
3429   mcIdType nbOfComp=tinyInfoI[1];
3430   if(nbOfTuple!=-1 || nbOfComp!=-1)
3431     {
3432       alloc(nbOfTuple,nbOfComp);
3433       return true;
3434     }
3435   return false;
3436 }
3437
3438 /*!
3439  * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3440  */
3441 void DataArrayDouble::finishUnserialization(const std::vector<mcIdType>& tinyInfoI, const std::vector<std::string>& tinyInfoS)
3442 {
3443   setName(tinyInfoS[0]);
3444   if(isAllocated())
3445     {
3446       std::size_t nbOfCompo(getNumberOfComponents());
3447       for(std::size_t i=0;i<nbOfCompo;i++)
3448         setInfoOnComponent(i,tinyInfoS[i+1]);
3449     }
3450 }
3451
3452 /*!
3453  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn
3454  * around an axe ( \a center, \a vect) and with angle \a angle.
3455  */
3456 void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3457 {
3458   if(!center || !vect)
3459     throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !");
3460   double sina(sin(angle));
3461   double cosa(cos(angle));
3462   double vectorNorm[3];
3463   double matrix[9];
3464   double matrixTmp[9];
3465   double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
3466   if(norm<std::numeric_limits<double>::min())
3467     throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !");
3468   std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
3469   //rotation matrix computation
3470   matrix[0]=cosa; matrix[1]=0.; matrix[2]=0.; matrix[3]=0.; matrix[4]=cosa; matrix[5]=0.; matrix[6]=0.; matrix[7]=0.; matrix[8]=cosa;
3471   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
3472   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
3473   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
3474   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
3475   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3476   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
3477   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
3478   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
3479   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
3480   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3481   //rotation matrix computed.
3482   double tmp[3];
3483   for(mcIdType i=0; i<nbNodes; i++)
3484     {
3485       std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<double>());
3486       coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
3487       coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
3488       coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
3489     }
3490 }
3491
3492 void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3493 {
3494   double matrix[9],matrix2[9],matrix3[9];
3495   double vect[3],crossVect[3];
3496   INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3497   crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3498   crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3499   crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3500   double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3501   matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni;
3502   matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni;
3503   matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni;
3504   matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv;
3505   matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc;
3506   matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni;
3507   for(mcIdType i=0;i<3;i++)
3508     for(mcIdType j=0;j<3;j++)
3509       {
3510         double val(0.);
3511         for(mcIdType k=0;k<3;k++)
3512           val+=matrix[3*i+k]*matrix2[3*k+j];
3513         matrix3[3*i+j]=val;
3514       }
3515   //rotation matrix computed.
3516   double tmp[3];
3517   for(mcIdType i=0; i<nbNodes; i++)
3518     {
3519       std::transform(coordsIn+i*3,coordsIn+(i+1)*3,point,tmp,std::minus<double>());
3520       coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0];
3521       coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1];
3522       coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2];
3523     }
3524 }
3525
3526 void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9])
3527 {
3528   double vect[3],crossVect[3];
3529   INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3530   crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3531   crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3532   crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3533   double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3534   baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv;
3535   baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc;
3536   baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni;
3537 }
3538
3539 /*!
3540  * \param [in] seg2 : coordinates of input seg2 expected to have spacedim==2
3541  * \param [in] tri3 : coordinates of input tri3 also expected to have spacedim==2
3542  * \param [out] coeffs : the result of integration normalized to 1. along \a seg2 inside tri3 sorted by the node id of \a tri3
3543  * \param [out] length : the length of seg2. That is too say the length of integration
3544  */
3545 void DataArrayDouble::ComputeIntegralOfSeg2IntoTri3(const double seg2[4], const double tri3[6], double coeffs[3], double& length)
3546 {
3547   length=INTERP_KERNEL::norme_vecteur(seg2,seg2+2);
3548   double mid[2];
3549   INTERP_KERNEL::mid_of_seg2(seg2,seg2+2,mid);
3550   INTERP_KERNEL::barycentric_coords<2>(tri3,mid,coeffs); // integral along seg2 is equal to value at the center of SEG2 !
3551 }
3552
3553 /*!
3554  * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords
3555  * around the center point \a center and with angle \a angle.
3556  */
3557 void DataArrayDouble::Rotate2DAlg(const double *center, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3558 {
3559   double cosa=cos(angle);
3560   double sina=sin(angle);
3561   double matrix[4];
3562   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
3563   double tmp[2];
3564   for(mcIdType i=0; i<nbNodes; i++)
3565     {
3566       std::transform(coordsIn+i*2,coordsIn+(i+1)*2,center,tmp,std::minus<double>());
3567       coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
3568       coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
3569     }
3570 }
3571
3572 DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):DataArrayIterator<double>(da)
3573 {
3574 }
3575
3576 DataArrayDoubleTuple::DataArrayDoubleTuple(double *pt, std::size_t nbOfComp):DataArrayTuple<double>(pt,nbOfComp)
3577 {
3578 }
3579
3580
3581 std::string DataArrayDoubleTuple::repr() const
3582 {
3583   std::ostringstream oss; oss.precision(17); oss << "(";
3584   for(std::size_t i=0;i<_nb_of_compo-1;i++)
3585     oss << _pt[i] << ", ";
3586   oss << _pt[_nb_of_compo-1] << ")";
3587   return oss.str();
3588 }
3589
3590 double DataArrayDoubleTuple::doubleValue() const
3591 {
3592   return this->zeValue();
3593 }
3594
3595 /*!
3596  * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayDouble::decrRef.
3597  * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayDouble::useArray with ownership set to \b false.
3598  * This method throws an INTERP_KERNEL::Exception is it is impossible to match sizes of \b this that is too say \b nbOfCompo=this->_nb_of_elem and \bnbOfTuples==1 or
3599  * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3600  */
3601 DataArrayDouble *DataArrayDoubleTuple::buildDADouble(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3602 {
3603   return this->buildDA(nbOfTuples,nbOfCompo);
3604 }
3605
3606 /*!
3607  * Returns a full copy of \a this. For more info on copying data arrays see
3608  * \ref MEDCouplingArrayBasicsCopyDeep.
3609  *  \return DataArrayInt * - a new instance of DataArrayInt.
3610  */
3611 DataArrayInt32 *DataArrayInt32::deepCopy() const
3612 {
3613   return new DataArrayInt32(*this);
3614 }
3615
3616 DataArrayInt32Iterator *DataArrayInt32::iterator()
3617 {
3618   return new DataArrayInt32Iterator(this);
3619 }
3620
3621
3622 DataArrayInt32Iterator::DataArrayInt32Iterator(DataArrayInt32 *da):DataArrayIterator<Int32>(da)
3623 {
3624 }
3625
3626 DataArrayInt32Tuple::DataArrayInt32Tuple(Int32 *pt, std::size_t nbOfComp):DataArrayTuple<Int32>(pt,nbOfComp)
3627 {
3628 }
3629
3630 std::string DataArrayInt32Tuple::repr() const
3631 {
3632   std::ostringstream oss; oss << "(";
3633   for(std::size_t i=0;i<_nb_of_compo-1;i++)
3634     oss << _pt[i] << ", ";
3635   oss << _pt[_nb_of_compo-1] << ")";
3636   return oss.str();
3637 }
3638
3639 Int32 DataArrayInt32Tuple::intValue() const
3640 {
3641   return this->zeValue();
3642 }
3643
3644 /*!
3645  * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayInt::decrRef.
3646  * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayInt::useArray with ownership set to \b false.
3647  * This method throws an INTERP_KERNEL::Exception is it is impossible to match sizes of \b this that is too say \b nbOfCompo=this->_nb_of_elem and \bnbOfTuples==1 or
3648  * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3649  */
3650 DataArrayInt32 *DataArrayInt32Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3651 {
3652   return this->buildDA(nbOfTuples,nbOfCompo);
3653 }
3654
3655 DataArrayInt64Iterator *DataArrayInt64::iterator()
3656 {
3657   return new DataArrayInt64Iterator(this);
3658 }
3659
3660
3661 DataArrayInt64Iterator::DataArrayInt64Iterator(DataArrayInt64 *da):DataArrayIterator<Int64>(da)
3662 {
3663 }
3664
3665 DataArrayInt64Tuple::DataArrayInt64Tuple(Int64 *pt, std::size_t nbOfComp):DataArrayTuple<Int64>(pt,nbOfComp)
3666 {
3667 }
3668
3669 std::string DataArrayInt64Tuple::repr() const
3670 {
3671   std::ostringstream oss; oss << "(";
3672   for(std::size_t i=0;i<_nb_of_compo-1;i++)
3673     oss << _pt[i] << ", ";
3674   oss << _pt[_nb_of_compo-1] << ")";
3675   return oss.str();
3676 }
3677
3678 Int64 DataArrayInt64Tuple::intValue() const
3679 {
3680   return this->zeValue();
3681 }
3682
3683 DataArrayInt64 *DataArrayInt64Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3684 {
3685   return this->buildDA(nbOfTuples,nbOfCompo);
3686 }
3687
3688
3689 DataArrayInt64 *DataArrayInt64::deepCopy() const
3690 {
3691   return new DataArrayInt64(*this);
3692 }