src/MEDCoupling/MEDCouplingMemArray.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF
   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 DataArrayDouble *DataArrayDouble::operatePerTuple(std::function<double(const double *bg, const double *endd)> func) const
2341 {
2342   checkAllocated();
2343   std::size_t nbOfComp(getNumberOfComponents());
2344   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2345   mcIdType nbOfTuple(getNumberOfTuples());
2346   ret->alloc(nbOfTuple,1);
2347   const double *src=getConstPointer();
2348   double *dest=ret->getPointer();
2349   for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=nbOfComp)
2350     *dest=func(src,src+nbOfComp);
2351   return ret.retn();
2352 }
2353
2354 /*!
2355  * Computes the maximal value within every tuple of \a this array.
2356  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2357  *          same number of tuples as \a this array and one component.
2358  *          The caller is to delete this result array using decrRef() as it is no more
2359  *          needed.
2360  *  \throw If \a this is not allocated.
2361  *  \sa DataArrayDouble::maxPerTupleWithCompoId, DataArrayDouble::minPerTuple
2362  */
2363 DataArrayDouble *DataArrayDouble::maxPerTuple() const
2364 {
2365   return this->operatePerTuple([](const double *bg, const double *endd) { return *std::max_element(bg,endd); });
2366 }
2367
2368 /*!
2369  * Computes the minimal value within every tuple of \a this array.
2370  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2371  *          same number of tuples as \a this array and one component.
2372  *          The caller is to delete this result array using decrRef() as it is no more
2373  *          needed.
2374  *  \throw If \a this is not allocated.
2375  *  \sa DataArrayDouble::maxPerTuple
2376  */
2377 DataArrayDouble *DataArrayDouble::minPerTuple() const
2378 {
2379   return this->operatePerTuple([](const double *bg, const double *endd) { return *std::min_element(bg,endd); });
2380 }
2381
2382 /*!
2383  * Computes the maximal value within every tuple of \a this array and it returns the first component
2384  * id for each tuple that corresponds to the maximal value within the tuple.
2385  *
2386  *  \param [out] compoIdOfMaxPerTuple - the new new instance of DataArrayInt containing the
2387  *          same number of tuples and only one component.
2388  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2389  *          same number of tuples as \a this array and one component.
2390  *          The caller is to delete this result array using decrRef() as it is no more
2391  *          needed.
2392  *  \throw If \a this is not allocated.
2393  *  \sa DataArrayDouble::maxPerTuple
2394  */
2395 DataArrayDouble *DataArrayDouble::maxPerTupleWithCompoId(DataArrayIdType* &compoIdOfMaxPerTuple) const
2396 {
2397   checkAllocated();
2398   std::size_t nbOfComp(getNumberOfComponents());
2399   MCAuto<DataArrayDouble> ret0=DataArrayDouble::New();
2400   MCAuto<DataArrayIdType> ret1=DataArrayIdType::New();
2401   mcIdType nbOfTuple=getNumberOfTuples();
2402   ret0->alloc(nbOfTuple,1); ret1->alloc(nbOfTuple,1);
2403   const double *src=getConstPointer();
2404   double *dest=ret0->getPointer(); mcIdType *dest1=ret1->getPointer();
2405   for(mcIdType i=0;i<nbOfTuple;i++,dest++,dest1++,src+=nbOfComp)
2406     {
2407       const double *loc=std::max_element(src,src+nbOfComp);
2408       *dest=*loc;
2409       *dest1=ToIdType(std::distance(src,loc));
2410     }
2411   compoIdOfMaxPerTuple=ret1.retn();
2412   return ret0.retn();
2413 }
2414
2415 /*!
2416  * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c this->getNumberOfTuples() tuples.
2417  * \n This returned array contains the euclidian distance for each tuple in \a this.
2418  * \n So the returned array can be seen as a dense symmetrical matrix whose diagonal elements are equal to 0.
2419  * \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)
2420  *
2421  * \warning use this method with care because it can leads to big amount of consumed memory !
2422  *
2423  * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2424  *
2425  * \throw If \a this is not allocated.
2426  *
2427  * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrixWith
2428  */
2429 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrix() const
2430 {
2431   checkAllocated();
2432   std::size_t nbOfComp(getNumberOfComponents());
2433   mcIdType nbOfTuples(getNumberOfTuples());
2434   const double *inData=getConstPointer();
2435   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2436   ret->alloc(nbOfTuples*nbOfTuples,1);
2437   double *outData=ret->getPointer();
2438   for(mcIdType i=0;i<nbOfTuples;i++)
2439     {
2440       outData[i*nbOfTuples+i]=0.;
2441       for(mcIdType j=i+1;j<nbOfTuples;j++)
2442         {
2443           double dist=0.;
2444           for(std::size_t k=0;k<nbOfComp;k++)
2445             { double delta=inData[i*nbOfComp+k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2446           dist=sqrt(dist);
2447           outData[i*nbOfTuples+j]=dist;
2448           outData[j*nbOfTuples+i]=dist;
2449         }
2450     }
2451   return ret.retn();
2452 }
2453
2454 /*!
2455  * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c other->getNumberOfTuples() tuples.
2456  * \n This returned array contains the euclidian distance for each tuple in \a other with each tuple in \a this.
2457  * \n So the returned array can be seen as a dense rectangular matrix with \c other->getNumberOfTuples() rows and \c this->getNumberOfTuples() columns.
2458  * \n Output rectangular matrix is sorted along rows.
2459  * \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)
2460  *
2461  * \warning use this method with care because it can leads to big amount of consumed memory !
2462  *
2463  * \param [in] other DataArrayDouble instance having same number of components than \a this.
2464  * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2465  *
2466  * \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.
2467  *
2468  * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrix
2469  */
2470 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const
2471 {
2472   if(!other)
2473     throw INTERP_KERNEL::Exception("DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : input parameter is null !");
2474   checkAllocated();
2475   other->checkAllocated();
2476   std::size_t nbOfComp(getNumberOfComponents());
2477   std::size_t otherNbOfComp(other->getNumberOfComponents());
2478   if(nbOfComp!=otherNbOfComp)
2479     {
2480       std::ostringstream oss; oss << "DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : this nb of compo=" << nbOfComp << " and other nb of compo=" << otherNbOfComp << ". It should match !";
2481       throw INTERP_KERNEL::Exception(oss.str().c_str());
2482     }
2483   mcIdType nbOfTuples(getNumberOfTuples());
2484   mcIdType otherNbOfTuples(other->getNumberOfTuples());
2485   const double *inData=getConstPointer();
2486   const double *inDataOther=other->getConstPointer();
2487   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2488   ret->alloc(otherNbOfTuples*nbOfTuples,1);
2489   double *outData=ret->getPointer();
2490   for(mcIdType i=0;i<otherNbOfTuples;i++,inDataOther+=nbOfComp)
2491     {
2492       for(mcIdType j=0;j<nbOfTuples;j++)
2493         {
2494           double dist=0.;
2495           for(std::size_t k=0;k<nbOfComp;k++)
2496             { double delta=inDataOther[k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2497           dist=sqrt(dist);
2498           outData[i*nbOfTuples+j]=dist;
2499         }
2500     }
2501   return ret.retn();
2502 }
2503
2504 /*!
2505  * This method expects that \a this stores 3 tuples containing 2 components each.
2506  * Each of this tuples represent a point into 2D space.
2507  * 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).
2508  * 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[.
2509  *
2510  *  \throw If \a this is not allocated.
2511  *  \throw If \a this has not 3 tuples of 2 components
2512  *  \throw If tuples/points in \a this are aligned
2513  */
2514 void DataArrayDouble::asArcOfCircle(double center[2], double& radius, double& ang) const
2515 {
2516   checkAllocated();
2517   INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(1e-14);
2518   if(getNumberOfTuples()!=3 && getNumberOfComponents()!=2)
2519     throw INTERP_KERNEL::Exception("DataArrayDouble::asArcCircle : this method expects");
2520   const double *pt(begin());
2521   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]));
2522   {
2523     INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeLin> e1(new INTERP_KERNEL::EdgeLin(n0,n2)),e2(new INTERP_KERNEL::EdgeLin(n2,n1));
2524     INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
2525     bool colinearity(inters.areColinears());
2526     if(colinearity)
2527       throw INTERP_KERNEL::Exception("DataArrayDouble::asArcOfCircle : 3 points in this have been detected as colinear !");
2528   }
2529   INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeArcCircle> ret(new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1));
2530   const double *c(ret->getCenter());
2531   center[0]=c[0]; center[1]=c[1];
2532   radius=ret->getRadius();
2533   ang=ret->getAngle();
2534 }
2535
2536 /*!
2537  * Sorts value within every tuple of \a this array.
2538  *  \param [in] asc - if \a true, the values are sorted in ascending order, else,
2539  *              in descending order.
2540  *  \throw If \a this is not allocated.
2541  */
2542 void DataArrayDouble::sortPerTuple(bool asc)
2543 {
2544   checkAllocated();
2545   double *pt=getPointer();
2546   mcIdType nbOfTuple(getNumberOfTuples());
2547   std::size_t nbOfComp(getNumberOfComponents());
2548   if(asc)
2549     for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2550       std::sort(pt,pt+nbOfComp);
2551   else
2552     for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2553       std::sort(pt,pt+nbOfComp,std::greater<double>());
2554   declareAsNew();
2555 }
2556
2557 /*!
2558  * Modify all elements of \a this array, so that
2559  * an element _x_ becomes \f$ numerator / x \f$.
2560  *  \warning If an exception is thrown because of presence of 0.0 element in \a this
2561  *           array, all elements processed before detection of the zero element remain
2562  *           modified.
2563  *  \param [in] numerator - the numerator used to modify array elements.
2564  *  \throw If \a this is not allocated.
2565  *  \throw If there is an element equal to 0.0 in \a this array.
2566  */
2567 void DataArrayDouble::applyInv(double numerator)
2568 {
2569   checkAllocated();
2570   double *ptr=getPointer();
2571   std::size_t nbOfElems=getNbOfElems();
2572   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2573     {
2574       if(std::abs(*ptr)>std::numeric_limits<double>::min())
2575         {
2576           *ptr=numerator/(*ptr);
2577         }
2578       else
2579         {
2580           std::ostringstream oss; oss << "DataArrayDouble::applyInv : presence of null value in tuple #" << i/getNumberOfComponents() << " component #" << i%getNumberOfComponents();
2581           oss << " !";
2582           throw INTERP_KERNEL::Exception(oss.str().c_str());
2583         }
2584     }
2585   declareAsNew();
2586 }
2587
2588 /*!
2589  * Modify all elements of \a this array, so that
2590  * an element _x_ becomes <em> val ^ x </em>. Contrary to DataArrayInt::applyPow
2591  * all values in \a this have to be >= 0 if val is \b not integer.
2592  *  \param [in] val - the value used to apply pow on all array elements.
2593  *  \throw If \a this is not allocated.
2594  *  \warning If an exception is thrown because of presence of 0 element in \a this
2595  *           array and \a val is \b not integer, all elements processed before detection of the zero element remain
2596  *           modified.
2597  */
2598 void DataArrayDouble::applyPow(double val)
2599 {
2600   checkAllocated();
2601   double *ptr=getPointer();
2602   std::size_t nbOfElems=getNbOfElems();
2603   int val2=(int)val;
2604   bool isInt=((double)val2)==val;
2605   if(!isInt)
2606     {
2607       for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2608         {
2609           if(*ptr>=0)
2610             *ptr=pow(*ptr,val);
2611           else
2612             {
2613               std::ostringstream oss; oss << "DataArrayDouble::applyPow (double) : At elem # " << i << " value is " << *ptr << " ! must be >=0. !";
2614               throw INTERP_KERNEL::Exception(oss.str().c_str());
2615             }
2616         }
2617     }
2618   else
2619     {
2620       for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2621         *ptr=pow(*ptr,val2);
2622     }
2623   declareAsNew();
2624 }
2625
2626 /*!
2627  * Modify all elements of \a this array, so that
2628  * an element _x_ becomes \f$ val ^ x \f$.
2629  *  \param [in] val - the value used to apply pow on all array elements.
2630  *  \throw If \a this is not allocated.
2631  *  \throw If \a val < 0.
2632  *  \warning If an exception is thrown because of presence of 0 element in \a this
2633  *           array, all elements processed before detection of the zero element remain
2634  *           modified.
2635  */
2636 void DataArrayDouble::applyRPow(double val)
2637 {
2638   checkAllocated();
2639   if(val<0.)
2640     throw INTERP_KERNEL::Exception("DataArrayDouble::applyRPow : the input value has to be >= 0 !");
2641   double *ptr=getPointer();
2642   std::size_t nbOfElems=getNbOfElems();
2643   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2644     *ptr=pow(val,*ptr);
2645   declareAsNew();
2646 }
2647
2648 /*!
2649  * Returns a new DataArrayDouble created from \a this one by applying \a
2650  * FunctionToEvaluate to every tuple of \a this array. Textual data is not copied.
2651  * For more info see \ref MEDCouplingArrayApplyFunc
2652  *  \param [in] nbOfComp - number of components in the result array.
2653  *  \param [in] func - the \a FunctionToEvaluate declared as
2654  *              \c bool (*\a func)(\c const \c double *\a pos, \c double *\a res),
2655  *              where \a pos points to the first component of a tuple of \a this array
2656  *              and \a res points to the first component of a tuple of the result array.
2657  *              Note that length (number of components) of \a pos can differ from
2658  *              that of \a res.
2659  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2660  *          same number of tuples as \a this array.
2661  *          The caller is to delete this result array using decrRef() as it is no more
2662  *          needed.
2663  *  \throw If \a this is not allocated.
2664  *  \throw If \a func returns \a false.
2665  */
2666 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, FunctionToEvaluate func) const
2667 {
2668   checkAllocated();
2669   DataArrayDouble *newArr=DataArrayDouble::New();
2670   mcIdType nbOfTuples(getNumberOfTuples());
2671   std::size_t oldNbOfComp(getNumberOfComponents());
2672   newArr->alloc(nbOfTuples,nbOfComp);
2673   const double *ptr=getConstPointer();
2674   double *ptrToFill=newArr->getPointer();
2675   for(mcIdType i=0;i<nbOfTuples;i++)
2676     {
2677       if(!func(ptr+i*oldNbOfComp,ptrToFill+i*nbOfComp))
2678         {
2679           std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2680           std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2681           oss << ") : Evaluation of function failed !";
2682           newArr->decrRef();
2683           throw INTERP_KERNEL::Exception(oss.str().c_str());
2684         }
2685     }
2686   return newArr;
2687 }
2688
2689 /*!
2690  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2691  * tuple of \a this array. Textual data is not copied.
2692  * For more info see \ref MEDCouplingArrayApplyFunc1.
2693  *  \param [in] nbOfComp - number of components in the result array.
2694  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2695  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2696  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2697  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2698  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2699  *          same number of tuples as \a this array and \a nbOfComp components.
2700  *          The caller is to delete this result array using decrRef() as it is no more
2701  *          needed.
2702  *  \throw If \a this is not allocated.
2703  *  \throw If computing \a func fails.
2704  */
2705 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2706 {
2707   INTERP_KERNEL::ExprParser expr(func);
2708   expr.parse();
2709   std::set<std::string> vars;
2710   expr.getTrueSetOfVars(vars);
2711   std::vector<std::string> varsV(vars.begin(),vars.end());
2712   return applyFuncNamedCompo(nbOfComp,varsV,func,isSafe);
2713 }
2714
2715 /*!
2716  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2717  * tuple of \a this array. Textual data is not copied. This method works by tuples (whatever its size).
2718  * If \a this is a one component array, call applyFuncOnThis instead that performs the same work faster.
2719  *
2720  * For more info see \ref MEDCouplingArrayApplyFunc0.
2721  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2722  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2723  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2724  *                       If false the computation is carried on without any notification. When false the evaluation is a little faster.
2725  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2726  *          same number of tuples and components as \a this array.
2727  *          The caller is to delete this result array using decrRef() as it is no more
2728  *          needed.
2729  *  \sa applyFuncOnThis
2730  *  \throw If \a this is not allocated.
2731  *  \throw If computing \a func fails.
2732  */
2733 DataArrayDouble *DataArrayDouble::applyFunc(const std::string& func, bool isSafe) const
2734 {
2735   std::size_t nbOfComp(getNumberOfComponents());
2736   if(nbOfComp<=0)
2737     throw INTERP_KERNEL::Exception("DataArrayDouble::applyFunc : output number of component must be > 0 !");
2738   checkAllocated();
2739   mcIdType nbOfTuples(getNumberOfTuples());
2740   MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2741   newArr->alloc(nbOfTuples,nbOfComp);
2742   INTERP_KERNEL::ExprParser expr(func);
2743   expr.parse();
2744   std::set<std::string> vars;
2745   expr.getTrueSetOfVars(vars);
2746   if(vars.size()>1)
2747     {
2748       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 : ";
2749       std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2750       throw INTERP_KERNEL::Exception(oss.str().c_str());
2751     }
2752   if(vars.empty())
2753     {
2754       expr.prepareFastEvaluator();
2755       newArr->rearrange(1);
2756       newArr->fillWithValue(expr.evaluateDouble());
2757       newArr->rearrange(nbOfComp);
2758       return newArr.retn();
2759     }
2760   std::vector<std::string> vars2(vars.begin(),vars.end());
2761   double buff,*ptrToFill(newArr->getPointer());
2762   const double *ptr(begin());
2763   std::vector<double> stck;
2764   expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2765   expr.prepareFastEvaluator();
2766   if(!isSafe)
2767     {
2768       for(mcIdType i=0;i<nbOfTuples;i++)
2769         {
2770           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2771             {
2772               buff=*ptr;
2773               expr.evaluateDoubleInternal(stck);
2774               *ptrToFill=stck.back();
2775               stck.pop_back();
2776             }
2777         }
2778     }
2779   else
2780     {
2781       for(mcIdType i=0;i<nbOfTuples;i++)
2782         {
2783           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2784             {
2785               buff=*ptr;
2786               try
2787               {
2788                   expr.evaluateDoubleInternalSafe(stck);
2789               }
2790               catch(INTERP_KERNEL::Exception& e)
2791               {
2792                   std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2793                   oss << buff;
2794                   oss << ") : Evaluation of function failed !" << e.what();
2795                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2796               }
2797               *ptrToFill=stck.back();
2798               stck.pop_back();
2799             }
2800         }
2801     }
2802   return newArr.retn();
2803 }
2804
2805 /*!
2806  * This method is a non const method that modify the array in \a this.
2807  * This method only works on one component array. It means that function \a func must
2808  * contain at most one variable.
2809  * This method is a specialization of applyFunc method with one parameter on one component array.
2810  *
2811  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2812  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2813  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2814  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2815  *
2816  * \sa applyFunc
2817  */
2818 void DataArrayDouble::applyFuncOnThis(const std::string& func, bool isSafe)
2819 {
2820   std::size_t nbOfComp(getNumberOfComponents());
2821   if(nbOfComp<=0)
2822     throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncOnThis : output number of component must be > 0 !");
2823   checkAllocated();
2824   mcIdType nbOfTuples(getNumberOfTuples());
2825   INTERP_KERNEL::ExprParser expr(func);
2826   expr.parse();
2827   std::set<std::string> vars;
2828   expr.getTrueSetOfVars(vars);
2829   if(vars.size()>1)
2830     {
2831       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 : ";
2832       std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2833       throw INTERP_KERNEL::Exception(oss.str().c_str());
2834     }
2835   if(vars.empty())
2836     {
2837       expr.prepareFastEvaluator();
2838       std::vector<std::string> compInfo(getInfoOnComponents());
2839       rearrange(1);
2840       fillWithValue(expr.evaluateDouble());
2841       rearrange(nbOfComp);
2842       setInfoOnComponents(compInfo);
2843       return ;
2844     }
2845   std::vector<std::string> vars2(vars.begin(),vars.end());
2846   double buff,*ptrToFill(getPointer());
2847   const double *ptr(begin());
2848   std::vector<double> stck;
2849   expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2850   expr.prepareFastEvaluator();
2851   if(!isSafe)
2852     {
2853       for(mcIdType i=0;i<nbOfTuples;i++)
2854         {
2855           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2856             {
2857               buff=*ptr;
2858               expr.evaluateDoubleInternal(stck);
2859               *ptrToFill=stck.back();
2860               stck.pop_back();
2861             }
2862         }
2863     }
2864   else
2865     {
2866       for(mcIdType i=0;i<nbOfTuples;i++)
2867         {
2868           for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2869             {
2870               buff=*ptr;
2871               try
2872               {
2873                   expr.evaluateDoubleInternalSafe(stck);
2874               }
2875               catch(INTERP_KERNEL::Exception& e)
2876               {
2877                   std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2878                   oss << buff;
2879                   oss << ") : Evaluation of function failed !" << e.what();
2880                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2881               }
2882               *ptrToFill=stck.back();
2883               stck.pop_back();
2884             }
2885         }
2886     }
2887 }
2888
2889 /*!
2890  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2891  * tuple of \a this array. Textual data is not copied.
2892  * For more info see \ref MEDCouplingArrayApplyFunc2.
2893  *  \param [in] nbOfComp - number of components in the result array.
2894  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2895  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2896  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2897  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2898  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2899  *          same number of tuples as \a this array.
2900  *          The caller is to delete this result array using decrRef() as it is no more
2901  *          needed.
2902  *  \throw If \a this is not allocated.
2903  *  \throw If \a func contains vars that are not in \a this->getInfoOnComponent().
2904  *  \throw If computing \a func fails.
2905  */
2906 DataArrayDouble *DataArrayDouble::applyFuncCompo(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2907 {
2908   return applyFuncNamedCompo(nbOfComp,getVarsOnComponent(),func,isSafe);
2909 }
2910
2911 /*!
2912  * Returns a new DataArrayDouble created from \a this one by applying a function to every
2913  * tuple of \a this array. Textual data is not copied.
2914  * For more info see \ref MEDCouplingArrayApplyFunc3.
2915  *  \param [in] nbOfComp - number of components in the result array.
2916  *  \param [in] varsOrder - sequence of vars defining their order.
2917  *  \param [in] func - the expression defining how to transform a tuple of \a this array.
2918  *              Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2919  *  \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2920  *              If false the computation is carried on without any notification. When false the evaluation is a little faster.
2921  *  \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2922  *          same number of tuples as \a this array.
2923  *          The caller is to delete this result array using decrRef() as it is no more
2924  *          needed.
2925  *  \throw If \a this is not allocated.
2926  *  \throw If \a func contains vars not in \a varsOrder.
2927  *  \throw If computing \a func fails.
2928  */
2929 DataArrayDouble *DataArrayDouble::applyFuncNamedCompo(std::size_t nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func, bool isSafe) const
2930 {
2931   if(nbOfComp<=0)
2932     throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncNamedCompo : output number of component must be > 0 !");
2933   std::vector<std::string> varsOrder2(varsOrder);
2934   std::size_t oldNbOfComp(getNumberOfComponents());
2935   for(std::size_t i=varsOrder.size();i<oldNbOfComp;i++)
2936     varsOrder2.push_back(std::string());
2937   checkAllocated();
2938   mcIdType nbOfTuples(getNumberOfTuples());
2939   INTERP_KERNEL::ExprParser expr(func);
2940   expr.parse();
2941   std::set<std::string> vars;
2942   expr.getTrueSetOfVars(vars);
2943   if(vars.size()>oldNbOfComp)
2944     {
2945       std::ostringstream oss; oss << "The field has " << oldNbOfComp << " components and there are ";
2946       oss << vars.size() << " variables : ";
2947       std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2948       throw INTERP_KERNEL::Exception(oss.str().c_str());
2949     }
2950   MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2951   newArr->alloc(nbOfTuples,nbOfComp);
2952   INTERP_KERNEL::AutoPtr<double> buff(new double[oldNbOfComp]);
2953   double *buffPtr(buff),*ptrToFill;
2954   std::vector<double> stck;
2955   for(std::size_t iComp=0;iComp<nbOfComp;iComp++)
2956     {
2957       expr.prepareExprEvaluationDouble(varsOrder2,(int)oldNbOfComp,(int)nbOfComp,(int)iComp,buffPtr,buffPtr+oldNbOfComp);
2958       expr.prepareFastEvaluator();
2959       const double *ptr(getConstPointer());
2960       ptrToFill=newArr->getPointer()+iComp;
2961       if(!isSafe)
2962         {
2963           for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2964             {
2965               std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2966               expr.evaluateDoubleInternal(stck);
2967               *ptrToFill=stck.back();
2968               stck.pop_back();
2969             }
2970         }
2971       else
2972         {
2973           for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2974             {
2975               std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2976               try
2977               {
2978                   expr.evaluateDoubleInternalSafe(stck);
2979                   *ptrToFill=stck.back();
2980                   stck.pop_back();
2981               }
2982               catch(INTERP_KERNEL::Exception& e)
2983               {
2984                   std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2985                   std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2986                   oss << ") : Evaluation of function failed !" << e.what();
2987                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2988               }
2989             }
2990         }
2991     }
2992   return newArr.retn();
2993 }
2994
2995 void DataArrayDouble::applyFuncFast32(const std::string& func)
2996 {
2997   checkAllocated();
2998   INTERP_KERNEL::ExprParser expr(func);
2999   expr.parse();
3000   char *funcStr=expr.compileX86();
3001   MYFUNCPTR funcPtr;
3002   *((void **)&funcPtr)=funcStr;//he he...
3003   //
3004   double *ptr=getPointer();
3005   std::size_t nbOfComp=getNumberOfComponents();
3006   mcIdType nbOfTuples=getNumberOfTuples();
3007   std::size_t nbOfElems=nbOfTuples*nbOfComp;
3008   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
3009     *ptr=funcPtr(*ptr);
3010   declareAsNew();
3011 }
3012
3013 void DataArrayDouble::applyFuncFast64(const std::string& func)
3014 {
3015   checkAllocated();
3016   INTERP_KERNEL::ExprParser expr(func);
3017   expr.parse();
3018   char *funcStr=expr.compileX86_64();
3019   MYFUNCPTR funcPtr;
3020   *((void **)&funcPtr)=funcStr;//he he...
3021   //
3022   double *ptr=getPointer();
3023   std::size_t nbOfComp=getNumberOfComponents();
3024   mcIdType nbOfTuples=getNumberOfTuples();
3025   std::size_t nbOfElems=nbOfTuples*nbOfComp;
3026   for(std::size_t i=0;i<nbOfElems;i++,ptr++)
3027     *ptr=funcPtr(*ptr);
3028   declareAsNew();
3029 }
3030
3031 /*!
3032  * \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.
3033  */
3034 MCAuto<DataArrayDouble> DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const
3035 {
3036   checkAllocated();
3037   if(getNumberOfComponents()!=3)
3038     throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !");
3039   mcIdType nbTuples(getNumberOfTuples());
3040   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3041   ret->alloc(nbTuples,3);
3042   Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer());
3043   return ret;
3044 }
3045
3046 DataArrayDoubleIterator *DataArrayDouble::iterator()
3047 {
3048   return new DataArrayDoubleIterator(this);
3049 }
3050
3051 /*!
3052  * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3053  * array whose values are within a given range. Textual data is not copied.
3054  *  \param [in] vmin - a lowest acceptable value (included).
3055  *  \param [in] vmax - a greatest acceptable value (included).
3056  *  \return DataArrayInt * - the new instance of DataArrayInt.
3057  *          The caller is to delete this result array using decrRef() as it is no more
3058  *          needed.
3059  *  \throw If \a this->getNumberOfComponents() != 1.
3060  *
3061  *  \sa DataArrayDouble::findIdsNotInRange
3062  *
3063  *  \if ENABLE_EXAMPLES
3064  *  \ref cpp_mcdataarraydouble_getidsinrange "Here is a C++ example".<br>
3065  *  \ref py_mcdataarraydouble_getidsinrange "Here is a Python example".
3066  *  \endif
3067  */
3068 DataArrayIdType *DataArrayDouble::findIdsInRange(double vmin, double vmax) const
3069 {
3070   checkAllocated();
3071   if(getNumberOfComponents()!=1)
3072     throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsInRange : this must have exactly one component !");
3073   const double *cptr(begin());
3074   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3075   mcIdType nbOfTuples(getNumberOfTuples());
3076   for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3077     if(*cptr>=vmin && *cptr<=vmax)
3078       ret->pushBackSilent(i);
3079   return ret.retn();
3080 }
3081
3082 /*!
3083  * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3084  * array whose values are not within a given range. Textual data is not copied.
3085  *  \param [in] vmin - a lowest not acceptable value (excluded).
3086  *  \param [in] vmax - a greatest not acceptable value (excluded).
3087  *  \return DataArrayInt * - the new instance of DataArrayInt.
3088  *          The caller is to delete this result array using decrRef() as it is no more
3089  *          needed.
3090  *  \throw If \a this->getNumberOfComponents() != 1.
3091  *
3092  *  \sa DataArrayDouble::findIdsInRange
3093  */
3094 DataArrayIdType *DataArrayDouble::findIdsNotInRange(double vmin, double vmax) const
3095 {
3096   checkAllocated();
3097   if(getNumberOfComponents()!=1)
3098     throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsNotInRange : this must have exactly one component !");
3099   const double *cptr(begin());
3100   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3101   mcIdType nbOfTuples(getNumberOfTuples());
3102   for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3103     if(*cptr<vmin || *cptr>vmax)
3104       ret->pushBackSilent(i);
3105   return ret.retn();
3106 }
3107
3108 /*!
3109  * Returns a new DataArrayDouble by concatenating two given arrays, so that (1) the number
3110  * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3111  * the number of component in the result array is same as that of each of given arrays.
3112  * Info on components is copied from the first of the given arrays. Number of components
3113  * in the given arrays must be  the same.
3114  *  \param [in] a1 - an array to include in the result array.
3115  *  \param [in] a2 - another array to include in the result array.
3116  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3117  *          The caller is to delete this result array using decrRef() as it is no more
3118  *          needed.
3119  *  \throw If both \a a1 and \a a2 are NULL.
3120  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents().
3121  */
3122 DataArrayDouble *DataArrayDouble::Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2)
3123 {
3124   std::vector<const DataArrayDouble *> tmp(2);
3125   tmp[0]=a1; tmp[1]=a2;
3126   return Aggregate(tmp);
3127 }
3128
3129 /*!
3130  * Returns a new DataArrayDouble by concatenating all given arrays, so that (1) the number
3131  * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3132  * the number of component in the result array is same as that of each of given arrays.
3133  * Info on components is copied from the first of the given arrays. Number of components
3134  * in the given arrays must be  the same.
3135  * If the number of non null of elements in \a arr is equal to one the returned object is a copy of it
3136  * not the object itself.
3137  *  \param [in] arr - a sequence of arrays to include in the result array.
3138  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3139  *          The caller is to delete this result array using decrRef() as it is no more
3140  *          needed.
3141  *  \throw If all arrays within \a arr are NULL.
3142  *  \throw If getNumberOfComponents() of arrays within \a arr.
3143  */
3144 DataArrayDouble *DataArrayDouble::Aggregate(const std::vector<const DataArrayDouble *>& arr)
3145 {
3146   std::vector<const DataArrayDouble *> a;
3147   for(std::vector<const DataArrayDouble *>::const_iterator it4=arr.begin();it4!=arr.end();it4++)
3148     if(*it4)
3149       a.push_back(*it4);
3150   if(a.empty())
3151     throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : input list must contain at least one NON EMPTY DataArrayDouble !");
3152   std::vector<const DataArrayDouble *>::const_iterator it=a.begin();
3153   std::size_t nbOfComp((*it)->getNumberOfComponents());
3154   mcIdType nbt=(*it++)->getNumberOfTuples();
3155   for(mcIdType i=1;it!=a.end();it++,i++)
3156     {
3157       if((*it)->getNumberOfComponents()!=nbOfComp)
3158         throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : Nb of components mismatch for array aggregation !");
3159       nbt+=(*it)->getNumberOfTuples();
3160     }
3161   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
3162   ret->alloc(nbt,nbOfComp);
3163   double *pt=ret->getPointer();
3164   for(it=a.begin();it!=a.end();it++)
3165     pt=std::copy((*it)->getConstPointer(),(*it)->getConstPointer()+(*it)->getNbOfElems(),pt);
3166   ret->copyStringInfoFrom(*(a[0]));
3167   return ret.retn();
3168 }
3169
3170 /*!
3171  * Returns a new DataArrayDouble containing a dot product of two given arrays, so that
3172  * the i-th tuple of the result array is a sum of products of j-th components of i-th
3173  * tuples of given arrays (\f$ a_i = \sum_{j=1}^n a1_j * a2_j \f$).
3174  * Info on components and name is copied from the first of the given arrays.
3175  * Number of tuples and components in the given arrays must be the same.
3176  *  \param [in] a1 - a given array.
3177  *  \param [in] a2 - another given array.
3178  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3179  *          The caller is to delete this result array using decrRef() as it is no more
3180  *          needed.
3181  *  \throw If either \a a1 or \a a2 is NULL.
3182  *  \throw If any given array is not allocated.
3183  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3184  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3185  */
3186 DataArrayDouble *DataArrayDouble::Dot(const DataArrayDouble *a1, const DataArrayDouble *a2)
3187 {
3188   if(!a1 || !a2)
3189     throw INTERP_KERNEL::Exception("DataArrayDouble::Dot : input DataArrayDouble instance is NULL !");
3190   a1->checkAllocated();
3191   a2->checkAllocated();
3192   std::size_t nbOfComp(a1->getNumberOfComponents());
3193   if(nbOfComp!=a2->getNumberOfComponents())
3194     throw INTERP_KERNEL::Exception("Nb of components mismatch for array Dot !");
3195   mcIdType nbOfTuple(a1->getNumberOfTuples());
3196   if(nbOfTuple!=a2->getNumberOfTuples())
3197     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Dot !");
3198   DataArrayDouble *ret=DataArrayDouble::New();
3199   ret->alloc(nbOfTuple,1);
3200   double *retPtr=ret->getPointer();
3201   const double *a1Ptr=a1->begin(),*a2Ptr(a2->begin());
3202   for(mcIdType i=0;i<nbOfTuple;i++)
3203     {
3204       double sum=0.;
3205       for(std::size_t j=0;j<nbOfComp;j++)
3206         sum+=a1Ptr[i*nbOfComp+j]*a2Ptr[i*nbOfComp+j];
3207       retPtr[i]=sum;
3208     }
3209   ret->setInfoOnComponent(0,a1->getInfoOnComponent(0));
3210   ret->setName(a1->getName());
3211   return ret;
3212 }
3213
3214 /*!
3215  * Returns a new DataArrayDouble containing a cross product of two given arrays, so that
3216  * the i-th tuple of the result array contains 3 components of a vector which is a cross
3217  * product of two vectors defined by the i-th tuples of given arrays.
3218  * Info on components is copied from the first of the given arrays.
3219  * Number of tuples in the given arrays must be the same.
3220  * Number of components in the given arrays must be 3.
3221  *  \param [in] a1 - a given array.
3222  *  \param [in] a2 - another given array.
3223  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3224  *          The caller is to delete this result array using decrRef() as it is no more
3225  *          needed.
3226  *  \throw If either \a a1 or \a a2 is NULL.
3227  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3228  *  \throw If \a a1->getNumberOfComponents() != 3
3229  *  \throw If \a a2->getNumberOfComponents() != 3
3230  */
3231 DataArrayDouble *DataArrayDouble::CrossProduct(const DataArrayDouble *a1, const DataArrayDouble *a2)
3232 {
3233   if(!a1 || !a2)
3234     throw INTERP_KERNEL::Exception("DataArrayDouble::CrossProduct : input DataArrayDouble instance is NULL !");
3235   std::size_t nbOfComp(a1->getNumberOfComponents());
3236   if(nbOfComp!=a2->getNumberOfComponents())
3237     throw INTERP_KERNEL::Exception("Nb of components mismatch for array crossProduct !");
3238   if(nbOfComp!=3)
3239     throw INTERP_KERNEL::Exception("Nb of components must be equal to 3 for array crossProduct !");
3240   mcIdType nbOfTuple(a1->getNumberOfTuples());
3241   if(nbOfTuple!=a2->getNumberOfTuples())
3242     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array crossProduct !");
3243   DataArrayDouble *ret=DataArrayDouble::New();
3244   ret->alloc(nbOfTuple,3);
3245   double *retPtr=ret->getPointer();
3246   const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3247   for(mcIdType i=0;i<nbOfTuple;i++)
3248     {
3249       retPtr[3*i]=a1Ptr[3*i+1]*a2Ptr[3*i+2]-a1Ptr[3*i+2]*a2Ptr[3*i+1];
3250       retPtr[3*i+1]=a1Ptr[3*i+2]*a2Ptr[3*i]-a1Ptr[3*i]*a2Ptr[3*i+2];
3251       retPtr[3*i+2]=a1Ptr[3*i]*a2Ptr[3*i+1]-a1Ptr[3*i+1]*a2Ptr[3*i];
3252     }
3253   ret->copyStringInfoFrom(*a1);
3254   return ret;
3255 }
3256
3257 /*!
3258  * Returns a new DataArrayDouble containing maximal values of two given arrays.
3259  * Info on components is copied from the first of the given arrays.
3260  * Number of tuples and components in the given arrays must be the same.
3261  *  \param [in] a1 - an array to compare values with another one.
3262  *  \param [in] a2 - another array to compare values with the first one.
3263  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3264  *          The caller is to delete this result array using decrRef() as it is no more
3265  *          needed.
3266  *  \throw If either \a a1 or \a a2 is NULL.
3267  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3268  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3269  */
3270 DataArrayDouble *DataArrayDouble::Max(const DataArrayDouble *a1, const DataArrayDouble *a2)
3271 {
3272   if(!a1 || !a2)
3273     throw INTERP_KERNEL::Exception("DataArrayDouble::Max : input DataArrayDouble instance is NULL !");
3274   std::size_t nbOfComp(a1->getNumberOfComponents());
3275   if(nbOfComp!=a2->getNumberOfComponents())
3276     throw INTERP_KERNEL::Exception("Nb of components mismatch for array Max !");
3277   mcIdType nbOfTuple(a1->getNumberOfTuples());
3278   if(nbOfTuple!=a2->getNumberOfTuples())
3279     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Max !");
3280   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3281   ret->alloc(nbOfTuple,nbOfComp);
3282   double *retPtr(ret->getPointer());
3283   const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3284   std::size_t nbElem(nbOfTuple*nbOfComp);
3285   for(std::size_t i=0;i<nbElem;i++)
3286     retPtr[i]=std::max(a1Ptr[i],a2Ptr[i]);
3287   ret->copyStringInfoFrom(*a1);
3288   return ret.retn();
3289 }
3290
3291 /*!
3292  * Returns a new DataArrayDouble containing minimal values of two given arrays.
3293  * Info on components is copied from the first of the given arrays.
3294  * Number of tuples and components in the given arrays must be the same.
3295  *  \param [in] a1 - an array to compare values with another one.
3296  *  \param [in] a2 - another array to compare values with the first one.
3297  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3298  *          The caller is to delete this result array using decrRef() as it is no more
3299  *          needed.
3300  *  \throw If either \a a1 or \a a2 is NULL.
3301  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3302  *  \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3303  */
3304 DataArrayDouble *DataArrayDouble::Min(const DataArrayDouble *a1, const DataArrayDouble *a2)
3305 {
3306   if(!a1 || !a2)
3307     throw INTERP_KERNEL::Exception("DataArrayDouble::Min : input DataArrayDouble instance is NULL !");
3308   std::size_t nbOfComp(a1->getNumberOfComponents());
3309   if(nbOfComp!=a2->getNumberOfComponents())
3310     throw INTERP_KERNEL::Exception("Nb of components mismatch for array min !");
3311   mcIdType nbOfTuple(a1->getNumberOfTuples());
3312   if(nbOfTuple!=a2->getNumberOfTuples())
3313     throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array min !");
3314   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3315   ret->alloc(nbOfTuple,nbOfComp);
3316   double *retPtr(ret->getPointer());
3317   const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3318   std::size_t nbElem(nbOfTuple*nbOfComp);
3319   for(std::size_t i=0;i<nbElem;i++)
3320     retPtr[i]=std::min(a1Ptr[i],a2Ptr[i]);
3321   ret->copyStringInfoFrom(*a1);
3322   return ret.retn();
3323 }
3324
3325 /*!
3326  * Returns a new DataArrayDouble that is the result of pow of two given arrays. There are 3
3327  * valid cases.
3328  *
3329  *  \param [in] a1 - an array to pow up.
3330  *  \param [in] a2 - another array to sum up.
3331  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
3332  *          The caller is to delete this result array using decrRef() as it is no more
3333  *          needed.
3334  *  \throw If either \a a1 or \a a2 is NULL.
3335  *  \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3336  *  \throw If \a a1->getNumberOfComponents() != 1 or \a a2->getNumberOfComponents() != 1.
3337  *  \throw If there is a negative value in \a a1.
3338  */
3339 DataArrayDouble *DataArrayDouble::Pow(const DataArrayDouble *a1, const DataArrayDouble *a2)
3340 {
3341   if(!a1 || !a2)
3342     throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : at least one of input instances is null !");
3343   mcIdType nbOfTuple=a1->getNumberOfTuples();
3344   mcIdType nbOfTuple2=a2->getNumberOfTuples();
3345   std::size_t nbOfComp=a1->getNumberOfComponents();
3346   std::size_t nbOfComp2=a2->getNumberOfComponents();
3347   if(nbOfTuple!=nbOfTuple2)
3348     throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of tuples mismatches !");
3349   if(nbOfComp!=1 || nbOfComp2!=1)
3350     throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of components of both arrays must be equal to 1 !");
3351   MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(nbOfTuple,1);
3352   const double *ptr1(a1->begin()),*ptr2(a2->begin());
3353   double *ptr=ret->getPointer();
3354   for(mcIdType i=0;i<nbOfTuple;i++,ptr1++,ptr2++,ptr++)
3355     {
3356       if(*ptr1>=0)
3357         {
3358           *ptr=pow(*ptr1,*ptr2);
3359         }
3360       else
3361         {
3362           std::ostringstream oss; oss << "DataArrayDouble::Pow : on tuple #" << i << " of a1 value is < 0 (" << *ptr1 << ") !";
3363           throw INTERP_KERNEL::Exception(oss.str().c_str());
3364         }
3365     }
3366   return ret.retn();
3367 }
3368
3369 /*!
3370  * Apply pow on values of another DataArrayDouble to values of \a this one.
3371  *
3372  *  \param [in] other - an array to pow to \a this one.
3373  *  \throw If \a other is NULL.
3374  *  \throw If \a this->getNumberOfTuples() != \a other->getNumberOfTuples()
3375  *  \throw If \a this->getNumberOfComponents() != 1 or \a other->getNumberOfComponents() != 1
3376  *  \throw If there is a negative value in \a this.
3377  */
3378 void DataArrayDouble::powEqual(const DataArrayDouble *other)
3379 {
3380   if(!other)
3381     throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : input instance is null !");
3382   mcIdType nbOfTuple=getNumberOfTuples();
3383   mcIdType nbOfTuple2=other->getNumberOfTuples();
3384   std::size_t nbOfComp=getNumberOfComponents();
3385   std::size_t nbOfComp2=other->getNumberOfComponents();
3386   if(nbOfTuple!=nbOfTuple2)
3387     throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of tuples mismatches !");
3388   if(nbOfComp!=1 || nbOfComp2!=1)
3389     throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of components of both arrays must be equal to 1 !");
3390   double *ptr=getPointer();
3391   const double *ptrc=other->begin();
3392   for(mcIdType i=0;i<nbOfTuple;i++,ptrc++,ptr++)
3393     {
3394       if(*ptr>=0)
3395         *ptr=pow(*ptr,*ptrc);
3396       else
3397         {
3398           std::ostringstream oss; oss << "DataArrayDouble::powEqual : on tuple #" << i << " of this value is < 0 (" << *ptr << ") !";
3399           throw INTERP_KERNEL::Exception(oss.str().c_str());
3400         }
3401     }
3402   declareAsNew();
3403 }
3404
3405 /*!
3406  * This method is \b NOT wrapped into python because it can be useful only for performance reasons in C++ context.
3407  * All values in \a this must be 0. or 1. within eps error. 0 means false, 1 means true.
3408  * If an another value than 0 or 1 appear (within eps precision) an INTERP_KERNEL::Exception will be thrown.
3409  *
3410  * \throw if \a this is not allocated.
3411  * \throw if \a this has not exactly one component.
3412  */
3413 std::vector<bool> DataArrayDouble::toVectorOfBool(double eps) const
3414 {
3415   checkAllocated();
3416   if(getNumberOfComponents()!=1)
3417     throw INTERP_KERNEL::Exception("DataArrayDouble::toVectorOfBool : must be applied on single component array !");
3418   mcIdType nbt(getNumberOfTuples());
3419   std::vector<bool> ret(nbt);
3420   const double *pt(begin());
3421   for(mcIdType i=0;i<nbt;i++)
3422     {
3423       if(fabs(pt[i])<eps)
3424         ret[i]=false;
3425       else if(fabs(pt[i]-1.)<eps)
3426         ret[i]=true;
3427       else
3428         {
3429           std::ostringstream oss; oss << "DataArrayDouble::toVectorOfBool : the tuple #" << i << " has value " << pt[i] << " is invalid ! must be 0. or 1. !";
3430           throw INTERP_KERNEL::Exception(oss.str().c_str());
3431         }
3432     }
3433   return ret;
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::getTinySerializationIntInformation(std::vector<mcIdType>& tinyInfo) const
3441 {
3442   tinyInfo.resize(2);
3443   if(isAllocated())
3444     {
3445       tinyInfo[0]=getNumberOfTuples();
3446       tinyInfo[1]=ToIdType(getNumberOfComponents());
3447     }
3448   else
3449     {
3450       tinyInfo[0]=-1;
3451       tinyInfo[1]=-1;
3452     }
3453 }
3454
3455 /*!
3456  * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3457  * Server side.
3458  */
3459 void DataArrayDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
3460 {
3461   if(isAllocated())
3462     {
3463       std::size_t nbOfCompo(getNumberOfComponents());
3464       tinyInfo.resize(nbOfCompo+1);
3465       tinyInfo[0]=getName();
3466       for(std::size_t i=0;i<nbOfCompo;i++)
3467         tinyInfo[i+1]=getInfoOnComponent(i);
3468     }
3469   else
3470     {
3471       tinyInfo.resize(1);
3472       tinyInfo[0]=getName();
3473     }
3474 }
3475
3476 /*!
3477  * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3478  * This method returns if a feeding is needed.
3479  */
3480 bool DataArrayDouble::resizeForUnserialization(const std::vector<mcIdType>& tinyInfoI)
3481 {
3482   mcIdType nbOfTuple=tinyInfoI[0];
3483   mcIdType nbOfComp=tinyInfoI[1];
3484   if(nbOfTuple!=-1 || nbOfComp!=-1)
3485     {
3486       alloc(nbOfTuple,nbOfComp);
3487       return true;
3488     }
3489   return false;
3490 }
3491
3492 /*!
3493  * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3494  */
3495 void DataArrayDouble::finishUnserialization(const std::vector<mcIdType>& tinyInfoI, const std::vector<std::string>& tinyInfoS)
3496 {
3497   setName(tinyInfoS[0]);
3498   if(isAllocated())
3499     {
3500       std::size_t nbOfCompo(getNumberOfComponents());
3501       for(std::size_t i=0;i<nbOfCompo;i++)
3502         setInfoOnComponent(i,tinyInfoS[i+1]);
3503     }
3504 }
3505
3506 /*!
3507  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn
3508  * around an axe ( \a center, \a vect) and with angle \a angle.
3509  */
3510 void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3511 {
3512   if(!center || !vect)
3513     throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !");
3514   double sina(sin(angle));
3515   double cosa(cos(angle));
3516   double vectorNorm[3];
3517   double matrix[9];
3518   double matrixTmp[9];
3519   double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
3520   if(norm<std::numeric_limits<double>::min())
3521     throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !");
3522   std::transform(vect,vect+3,vectorNorm,std::bind(std::multiplies<double>(),std::placeholders::_1,1/norm));
3523   //rotation matrix computation
3524   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;
3525   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
3526   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
3527   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
3528   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind(std::multiplies<double>(),std::placeholders::_1,1-cosa));
3529   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3530   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
3531   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
3532   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
3533   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind(std::multiplies<double>(),std::placeholders::_1,sina));
3534   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3535   //rotation matrix computed.
3536   double tmp[3];
3537   for(mcIdType i=0; i<nbNodes; i++)
3538     {
3539       std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<double>());
3540       coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
3541       coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
3542       coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
3543     }
3544 }
3545
3546 void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3547 {
3548   double matrix[9],matrix2[9],matrix3[9];
3549   double vect[3],crossVect[3];
3550   INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3551   crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3552   crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3553   crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3554   double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3555   matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni;
3556   matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni;
3557   matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni;
3558   matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv;
3559   matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc;
3560   matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni;
3561   for(mcIdType i=0;i<3;i++)
3562     for(mcIdType j=0;j<3;j++)
3563       {
3564         double val(0.);
3565         for(mcIdType k=0;k<3;k++)
3566           val+=matrix[3*i+k]*matrix2[3*k+j];
3567         matrix3[3*i+j]=val;
3568       }
3569   //rotation matrix computed.
3570   double tmp[3];
3571   for(mcIdType i=0; i<nbNodes; i++)
3572     {
3573       std::transform(coordsIn+i*3,coordsIn+(i+1)*3,point,tmp,std::minus<double>());
3574       coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0];
3575       coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1];
3576       coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2];
3577     }
3578 }
3579
3580 void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9])
3581 {
3582   double vect[3],crossVect[3];
3583   INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3584   crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3585   crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3586   crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3587   double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3588   baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv;
3589   baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc;
3590   baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni;
3591 }
3592
3593 /*!
3594  * \param [in] seg2 : coordinates of input seg2 expected to have spacedim==2
3595  * \param [in] tri3 : coordinates of input tri3 also expected to have spacedim==2
3596  * \param [out] coeffs : the result of integration normalized to 1. along \a seg2 inside tri3 sorted by the node id of \a tri3
3597  * \param [out] length : the length of seg2. That is too say the length of integration
3598  */
3599 void DataArrayDouble::ComputeIntegralOfSeg2IntoTri3(const double seg2[4], const double tri3[6], double coeffs[3], double& length)
3600 {
3601   length=INTERP_KERNEL::norme_vecteur(seg2,seg2+2);
3602   double mid[2];
3603   INTERP_KERNEL::mid_of_seg2(seg2,seg2+2,mid);
3604   INTERP_KERNEL::barycentric_coords<2>(tri3,mid,coeffs); // integral along seg2 is equal to value at the center of SEG2 !
3605 }
3606
3607 /*!
3608  * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords
3609  * around the center point \a center and with angle \a angle.
3610  */
3611 void DataArrayDouble::Rotate2DAlg(const double *center, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3612 {
3613   double cosa=cos(angle);
3614   double sina=sin(angle);
3615   double matrix[4];
3616   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
3617   double tmp[2];
3618   for(mcIdType i=0; i<nbNodes; i++)
3619     {
3620       std::transform(coordsIn+i*2,coordsIn+(i+1)*2,center,tmp,std::minus<double>());
3621       coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
3622       coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
3623     }
3624 }
3625
3626 DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):DataArrayIterator<double>(da)
3627 {
3628 }
3629
3630 DataArrayDoubleTuple::DataArrayDoubleTuple(double *pt, std::size_t nbOfComp):DataArrayTuple<double>(pt,nbOfComp)
3631 {
3632 }
3633
3634
3635 std::string DataArrayDoubleTuple::repr() const
3636 {
3637   std::ostringstream oss; oss.precision(17); oss << "(";
3638   for(std::size_t i=0;i<_nb_of_compo-1;i++)
3639     oss << _pt[i] << ", ";
3640   oss << _pt[_nb_of_compo-1] << ")";
3641   return oss.str();
3642 }
3643
3644 double DataArrayDoubleTuple::doubleValue() const
3645 {
3646   return this->zeValue();
3647 }
3648
3649 /*!
3650  * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayDouble::decrRef.
3651  * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayDouble::useArray with ownership set to \b false.
3652  * 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
3653  * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3654  */
3655 DataArrayDouble *DataArrayDoubleTuple::buildDADouble(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3656 {
3657   return this->buildDA(nbOfTuples,nbOfCompo);
3658 }
3659
3660 /*!
3661  * Returns a full copy of \a this. For more info on copying data arrays see
3662  * \ref MEDCouplingArrayBasicsCopyDeep.
3663  *  \return DataArrayInt * - a new instance of DataArrayInt.
3664  */
3665 DataArrayInt32 *DataArrayInt32::deepCopy() const
3666 {
3667   return new DataArrayInt32(*this);
3668 }
3669
3670 MCAuto<DataArrayInt64> DataArrayInt32::convertToInt64Arr() const
3671 {
3672   this->checkAllocated();
3673   MCAuto<DataArrayInt64> ret(DataArrayInt64::New());
3674   ret->alloc(this->getNumberOfTuples(),this->getNumberOfComponents());
3675   ret->copyStringInfoFrom(*this);
3676   const std::int32_t *pt(this->begin());
3677   std::for_each(ret->getPointer(),ret->getPointer()+ret->getNbOfElems(),[&pt](std::int64_t& val) { val = std::int64_t(*pt++); });
3678   return ret;
3679 }
3680
3681 DataArrayInt32Iterator *DataArrayInt32::iterator()
3682 {
3683   return new DataArrayInt32Iterator(this);
3684 }
3685
3686
3687 DataArrayInt32Iterator::DataArrayInt32Iterator(DataArrayInt32 *da):DataArrayIterator<Int32>(da)
3688 {
3689 }
3690
3691 DataArrayInt32Tuple::DataArrayInt32Tuple(Int32 *pt, std::size_t nbOfComp):DataArrayTuple<Int32>(pt,nbOfComp)
3692 {
3693 }
3694
3695 std::string DataArrayInt32Tuple::repr() const
3696 {
3697   std::ostringstream oss; oss << "(";
3698   for(std::size_t i=0;i<_nb_of_compo-1;i++)
3699     oss << _pt[i] << ", ";
3700   oss << _pt[_nb_of_compo-1] << ")";
3701   return oss.str();
3702 }
3703
3704 Int32 DataArrayInt32Tuple::intValue() const
3705 {
3706   return this->zeValue();
3707 }
3708
3709 /*!
3710  * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayInt::decrRef.
3711  * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayInt::useArray with ownership set to \b false.
3712  * 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
3713  * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3714  */
3715 DataArrayInt32 *DataArrayInt32Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3716 {
3717   return this->buildDA(nbOfTuples,nbOfCompo);
3718 }
3719
3720 MCAuto<DataArrayInt32> DataArrayInt64::convertToInt32Arr() const
3721 {
3722   this->checkAllocated();
3723   MCAuto<DataArrayInt32> ret(DataArrayInt32::New());
3724   ret->alloc(this->getNumberOfTuples(),this->getNumberOfComponents());
3725   ret->copyStringInfoFrom(*this);
3726   const std::int64_t *pt(this->begin());
3727   std::for_each(ret->getPointer(),ret->getPointer()+ret->getNbOfElems(),[&pt](std::int32_t& val) { val = std::int32_t(*pt++); });
3728   return ret;
3729 }
3730
3731 DataArrayInt64Iterator *DataArrayInt64::iterator()
3732 {
3733   return new DataArrayInt64Iterator(this);
3734 }
3735
3736
3737 DataArrayInt64Iterator::DataArrayInt64Iterator(DataArrayInt64 *da):DataArrayIterator<Int64>(da)
3738 {
3739 }
3740
3741 DataArrayInt64Tuple::DataArrayInt64Tuple(Int64 *pt, std::size_t nbOfComp):DataArrayTuple<Int64>(pt,nbOfComp)
3742 {
3743 }
3744
3745 std::string DataArrayInt64Tuple::repr() const
3746 {
3747   std::ostringstream oss; oss << "(";
3748   for(std::size_t i=0;i<_nb_of_compo-1;i++)
3749     oss << _pt[i] << ", ";
3750   oss << _pt[_nb_of_compo-1] << ")";
3751   return oss.str();
3752 }
3753
3754 Int64 DataArrayInt64Tuple::intValue() const
3755 {
3756   return this->zeValue();
3757 }
3758
3759 DataArrayInt64 *DataArrayInt64Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3760 {
3761   return this->buildDA(nbOfTuples,nbOfCompo);
3762 }
3763
3764
3765 DataArrayInt64 *DataArrayInt64::deepCopy() const
3766 {
3767   return new DataArrayInt64(*this);
3768 }