src/MEDCoupling/MEDCouplingMemArray.cxx

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