src/MEDCoupling/MEDCouplingMemArray.txx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #ifndef __PARAMEDMEM_MEDCOUPLINGMEMARRAY_TXX__
  22 #define __PARAMEDMEM_MEDCOUPLINGMEMARRAY_TXX__
  23
  24 #include "MEDCouplingMemArray.hxx"
  25 #include "NormalizedUnstructuredMesh.hxx"
  26 #include "InterpKernelException.hxx"
  27 #include "InterpolationUtils.hxx"
  28 #include "MCAuto.hxx"
  29
  30 #include <sstream>
  31 #include <cstdlib>
  32 #include <algorithm>
  33
  34 namespace MEDCoupling
  35 {
  36   template<class T>
  37   void MEDCouplingPointer<T>::setInternal(T *pointer)
  38   {
  39     _internal=pointer;
  40     _external=0;
  41   }
  42
  43   template<class T>
  44   void MEDCouplingPointer<T>::setExternal(const T *pointer)
  45   {
  46     _external=pointer;
  47     _internal=0;
  48   }
  49
  50   template<class T>
  51   MemArray<T>::MemArray(const MemArray<T>& other):_nb_of_elem(0),_nb_of_elem_alloc(0),_ownership(false),_dealloc(0),_param_for_deallocator(0)
  52   {
  53     if(!other._pointer.isNull())
  54       {
  55         _nb_of_elem_alloc=other._nb_of_elem;
  56         T *pointer=(T*)malloc(_nb_of_elem_alloc*sizeof(T));
  57         std::copy(other._pointer.getConstPointer(),other._pointer.getConstPointer()+other._nb_of_elem,pointer);
  58         useArray(pointer,true,C_DEALLOC,other._nb_of_elem);
  59       }
  60   }
  61
  62   template<class T>
  63   void MemArray<T>::useArray(const T *array, bool ownership, DeallocType type, std::size_t nbOfElem)
  64   {
  65     destroy();
  66     _nb_of_elem=nbOfElem;
  67     _nb_of_elem_alloc=nbOfElem;
  68     if(ownership)
  69       _pointer.setInternal(const_cast<T *>(array));
  70     else
  71       _pointer.setExternal(array);
  72     _ownership=ownership;
  73     _dealloc=BuildFromType(type);
  74   }
  75
  76   template<class T>
  77   void MemArray<T>::useExternalArrayWithRWAccess(const T *array, std::size_t nbOfElem)
  78   {
  79     destroy();
  80     _nb_of_elem=nbOfElem;
  81     _nb_of_elem_alloc=nbOfElem;
  82     _pointer.setInternal(const_cast<T *>(array));
  83     _ownership=false;
  84     _dealloc=CPPDeallocator;
  85   }
  86
  87   template<class T>
  88   void MemArray<T>::writeOnPlace(std::size_t id, T element0, const T *others, std::size_t sizeOfOthers)
  89   {
  90     if(id+sizeOfOthers>=_nb_of_elem_alloc)
  91       reserve(2*_nb_of_elem+sizeOfOthers+1);
  92     T *pointer=_pointer.getPointer();
  93     pointer[id]=element0;
  94     std::copy(others,others+sizeOfOthers,pointer+id+1);
  95     _nb_of_elem=std::max<std::size_t>(_nb_of_elem,id+sizeOfOthers+1);
  96   }
  97
  98   template<class T>
  99   template<class InputIterator>
 100   void MemArray<T>::insertAtTheEnd(InputIterator first, InputIterator last)
 101   {
 102     T *pointer=_pointer.getPointer();
 103     while(first!=last)
 104       {
 105         if(_nb_of_elem>=_nb_of_elem_alloc)
 106           {
 107             reserve(_nb_of_elem_alloc>0?2*_nb_of_elem_alloc:1);
 108             pointer=_pointer.getPointer();
 109           }
 110         pointer[_nb_of_elem++]=*first++;
 111       }
 112   }
 113
 114   template<class T>
 115   void MemArray<T>::pushBack(T elem)
 116   {
 117     if(_nb_of_elem>=_nb_of_elem_alloc)
 118       reserve(_nb_of_elem_alloc>0?2*_nb_of_elem_alloc:1);
 119     T *pt=getPointer();
 120     pt[_nb_of_elem++]=elem;
 121   }
 122
 123   template<class T>
 124   T MemArray<T>::popBack()
 125   {
 126     if(_nb_of_elem>0)
 127       {
 128         const T *pt=getConstPointer();
 129         return pt[--_nb_of_elem];
 130       }
 131     throw INTERP_KERNEL::Exception("MemArray::popBack : nothing to pop in array !");
 132   }
 133
 134   template<class T>
 135   void MemArray<T>::pack() const
 136   {
 137     (const_cast<MemArray<T> * >(this))->reserve(_nb_of_elem);
 138   }
 139
 140   template<class T>
 141   bool MemArray<T>::isEqual(const MemArray<T>& other, T prec, std::string& reason) const
 142   {
 143     std::ostringstream oss; oss.precision(15);
 144     if(_nb_of_elem!=other._nb_of_elem)
 145       {
 146         oss << "Number of elements in coarse data of DataArray mismatch : this=" << _nb_of_elem << " other=" << other._nb_of_elem;
 147         reason=oss.str();
 148         return false;
 149       }
 150     const T *pt1=_pointer.getConstPointer();
 151     const T *pt2=other._pointer.getConstPointer();
 152     if(pt1==0 && pt2==0)
 153       return true;
 154     if(pt1==0 || pt2==0)
 155       {
 156         oss << "coarse data pointer is defined for only one DataArray instance !";
 157         reason=oss.str();
 158         return false;
 159       }
 160     if(pt1==pt2)
 161       return true;
 162     for(std::size_t i=0;i<_nb_of_elem;i++)
 163       if(pt1[i]-pt2[i]<-prec || (pt1[i]-pt2[i])>prec)
 164         {
 165           oss << "The content of data differs at pos #" << i << " of coarse data ! this[i]=" << pt1[i] << " other[i]=" << pt2[i];
 166           reason=oss.str();
 167           return false;
 168         }
 169     return true;
 170   }
 171
 172   /*!
 173    * \param [in] sl is typically the number of components
 174    * \return True if a not null pointer is present, False if not.
 175    */
 176   template<class T>
 177   bool MemArray<T>::reprHeader(int sl, std::ostream& stream) const
 178   {
 179     stream << "Number of tuples : ";
 180     if(!_pointer.isNull())
 181       {
 182         if(sl!=0)
 183           stream << _nb_of_elem/sl << std::endl << "Internal memory facts : " << _nb_of_elem << "/" << _nb_of_elem_alloc;
 184         else
 185           stream << "Empty Data";
 186       }
 187     else
 188       stream << "No data";
 189     stream << "\n";
 190     stream << "Data content :\n";
 191     bool ret=!_pointer.isNull();
 192     if(!ret)
 193       stream << "No data !\n";
 194     return ret;
 195   }
 196
 197   /*!
 198    * \param [in] sl is typically the number of components
 199    */
 200   template<class T>
 201   void MemArray<T>::repr(int sl, std::ostream& stream) const
 202   {
 203     if(reprHeader(sl,stream))
 204       {
 205         const T *data=getConstPointer();
 206         if(_nb_of_elem!=0 && sl!=0)
 207           {
 208             std::size_t nbOfTuples=_nb_of_elem/std::abs(sl);
 209             for(std::size_t i=0;i<nbOfTuples;i++)
 210               {
 211                 stream << "Tuple #" << i << " : ";
 212                 std::copy(data,data+sl,std::ostream_iterator<T>(stream," "));
 213                 stream << "\n";
 214                 data+=sl;
 215               }
 216           }
 217         else
 218           stream << "Empty Data\n";
 219       }
 220   }
 221
 222   /*!
 223    * \param [in] sl is typically the number of components
 224    */
 225   template<class T>
 226   void MemArray<T>::reprZip(int sl, std::ostream& stream) const
 227   {
 228     stream << "Number of tuples : ";
 229     if(!_pointer.isNull())
 230       {
 231         if(sl!=0)
 232           stream << _nb_of_elem/sl;
 233         else
 234           stream << "Empty Data";
 235       }
 236     else
 237       stream << "No data";
 238     stream << "\n";
 239     stream << "Data content : ";
 240     const T *data=getConstPointer();
 241     if(!_pointer.isNull())
 242       {
 243         if(_nb_of_elem!=0 && sl!=0)
 244           {
 245             std::size_t nbOfTuples=_nb_of_elem/std::abs(sl);
 246             for(std::size_t i=0;i<nbOfTuples;i++)
 247               {
 248                 stream << "|";
 249                 std::copy(data,data+sl,std::ostream_iterator<T>(stream," "));
 250                 stream << "| ";
 251                 data+=sl;
 252               }
 253             stream << "\n";
 254           }
 255         else
 256           stream << "Empty Data\n";
 257       }
 258     else
 259       stream << "No data !\n";
 260   }
 261
 262   /*!
 263    * \param [in] sl is typically the number of components
 264    */
 265   template<class T>
 266   void MemArray<T>::reprNotTooLong(int sl, std::ostream& stream) const
 267   {
 268     if(reprHeader(sl,stream))
 269       {
 270         const T *data=getConstPointer();
 271         if(_nb_of_elem!=0 && sl!=0)
 272           {
 273             std::size_t nbOfTuples=_nb_of_elem/std::abs(sl);
 274             if(nbOfTuples<=1000)
 275               {
 276                 for(std::size_t i=0;i<nbOfTuples;i++)
 277                   {
 278                     stream << "Tuple #" << i << " : ";
 279                     std::copy(data,data+sl,std::ostream_iterator<T>(stream," "));
 280                     stream << "\n";
 281                     data+=sl;
 282                   }
 283               }
 284             else
 285               {// too much tuples -> print the 3 first tuples and 3 last.
 286                 stream << "Tuple #0 : ";
 287                 std::copy(data,data+sl,std::ostream_iterator<T>(stream," ")); stream << "\n";
 288                 stream << "Tuple #1 : ";
 289                 std::copy(data+sl,data+2*sl,std::ostream_iterator<T>(stream," ")); stream << "\n";
 290                 stream << "Tuple #2 : ";
 291                 std::copy(data+2*sl,data+3*sl,std::ostream_iterator<T>(stream," ")); stream << "\n";
 292                 stream << "...\n";
 293                 stream << "Tuple #" << nbOfTuples-3 << " : ";
 294                 std::copy(data+(nbOfTuples-3)*sl,data+(nbOfTuples-2)*sl,std::ostream_iterator<T>(stream," ")); stream << "\n";
 295                 stream << "Tuple #" << nbOfTuples-2 << " : ";
 296                 std::copy(data+(nbOfTuples-2)*sl,data+(nbOfTuples-1)*sl,std::ostream_iterator<T>(stream," ")); stream << "\n";
 297                 stream << "Tuple #" << nbOfTuples-1 << " : ";
 298                 std::copy(data+(nbOfTuples-1)*sl,data+nbOfTuples*sl,std::ostream_iterator<T>(stream," ")); stream << "\n";
 299               }
 300           }
 301         else
 302           stream << "Empty Data\n";
 303       }
 304   }
 305
 306   template<class T>
 307   void MemArray<T>::fillWithValue(const T& val)
 308   {
 309     T *pt=_pointer.getPointer();
 310     std::fill(pt,pt+_nb_of_elem,val);
 311   }
 312
 313   template<class T>
 314   T *MemArray<T>::fromNoInterlace(int nbOfComp) const
 315   {
 316     if(nbOfComp<1)
 317       throw INTERP_KERNEL::Exception("MemArray<T>::fromNoInterlace : number of components must be > 0 !");
 318     const T *pt=_pointer.getConstPointer();
 319     std::size_t nbOfTuples=_nb_of_elem/nbOfComp;
 320     T *ret=(T*)malloc(_nb_of_elem*sizeof(T));
 321     T *w=ret;
 322     for(std::size_t i=0;i<nbOfTuples;i++)
 323       for(int j=0;j<nbOfComp;j++,w++)
 324         *w=pt[j*nbOfTuples+i];
 325     return ret;
 326   }
 327
 328   template<class T>
 329   T *MemArray<T>::toNoInterlace(int nbOfComp) const
 330   {
 331     if(nbOfComp<1)
 332       throw INTERP_KERNEL::Exception("MemArray<T>::toNoInterlace : number of components must be > 0 !");
 333     const T *pt=_pointer.getConstPointer();
 334     std::size_t nbOfTuples=_nb_of_elem/nbOfComp;
 335     T *ret=(T*)malloc(_nb_of_elem*sizeof(T));
 336     T *w=ret;
 337     for(int i=0;i<nbOfComp;i++)
 338       for(std::size_t j=0;j<nbOfTuples;j++,w++)
 339         *w=pt[j*nbOfComp+i];
 340     return ret;
 341   }
 342
 343   template<class T>
 344   void MemArray<T>::sort(bool asc)
 345   {
 346     T *pt=_pointer.getPointer();
 347     if(asc)
 348       std::sort(pt,pt+_nb_of_elem);
 349     else
 350       {
 351         typename std::reverse_iterator<T *> it1(pt+_nb_of_elem);
 352         typename std::reverse_iterator<T *> it2(pt);
 353         std::sort(it1,it2);
 354       }
 355   }
 356
 357   template<class T>
 358   void MemArray<T>::reverse(int nbOfComp)
 359   {
 360     if(nbOfComp<1)
 361       throw INTERP_KERNEL::Exception("MemArray<T>::reverse : only supported with 'this' array with ONE or more than ONE component !");
 362     T *pt=_pointer.getPointer();
 363     if(nbOfComp==1)
 364       {
 365         std::reverse(pt,pt+_nb_of_elem);
 366         return ;
 367       }
 368     else
 369       {
 370         T *pt2=pt+_nb_of_elem-nbOfComp;
 371         std::size_t nbOfTuples=_nb_of_elem/nbOfComp;
 372         for(std::size_t i=0;i<nbOfTuples/2;i++,pt+=nbOfComp,pt2-=nbOfComp)
 373           {
 374             for(int j=0;j<nbOfComp;j++)
 375               std::swap(pt[j],pt2[j]);
 376           }
 377       }
 378   }
 379
 380   template<class T>
 381   void MemArray<T>::alloc(std::size_t nbOfElements)
 382   {
 383     destroy();
 384     _nb_of_elem=nbOfElements;
 385     _nb_of_elem_alloc=nbOfElements;
 386     _pointer.setInternal((T*)malloc(_nb_of_elem_alloc*sizeof(T)));
 387     _ownership=true;
 388     _dealloc=CDeallocator;
 389   }
 390
 391   /*!
 392    * This method performs systematically an allocation of \a newNbOfElements elements in \a this.
 393    * \a _nb_of_elem and \a _nb_of_elem_alloc will \b NOT be systematically equal (contrary to MemArray<T>::reAlloc method.
 394    * So after the call of this method \a _nb_of_elem will be equal tostd::min<std::size_t>(_nb_of_elem,newNbOfElements) and \a _nb_of_elem_alloc equal to
 395    * \a newNbOfElements. This method is typically used to perform a pushBack to avoid systematic allocations-copy-deallocation.
 396    * So after the call of this method the accessible content is perfectly set.
 397    *
 398    * So this method should not be confused with MemArray<T>::reserve that is close to MemArray<T>::reAlloc but not same.
 399    */
 400   template<class T>
 401   void MemArray<T>::reserve(std::size_t newNbOfElements)
 402   {
 403     if(_nb_of_elem_alloc==newNbOfElements)
 404       return ;
 405     T *pointer=(T*)malloc(newNbOfElements*sizeof(T));
 406     std::copy(_pointer.getConstPointer(),_pointer.getConstPointer()+std::min<std::size_t>(_nb_of_elem,newNbOfElements),pointer);
 407     if(_ownership)
 408       DestroyPointer(const_cast<T *>(_pointer.getConstPointer()),_dealloc,_param_for_deallocator);//Do not use getPointer because in case of _external
 409     _pointer.setInternal(pointer);
 410     _nb_of_elem=std::min<std::size_t>(_nb_of_elem,newNbOfElements);
 411     _nb_of_elem_alloc=newNbOfElements;
 412     _ownership=true;
 413     _dealloc=CDeallocator;
 414     _param_for_deallocator=0;
 415   }
 416
 417   /*!
 418    * This method performs systematically an allocation of \a newNbOfElements elements in \a this.
 419    * \a _nb_of_elem and \a _nb_of_elem_alloc will be equal even if only std::min<std::size_t>(_nb_of_elem,newNbOfElements) come from the .
 420    * The remaing part of the new allocated chunk are available but not set previouly !
 421    *
 422    * So this method should not be confused with MemArray<T>::reserve that is close to MemArray<T>::reAlloc but not same.
 423    */
 424   template<class T>
 425   void MemArray<T>::reAlloc(std::size_t newNbOfElements)
 426   {
 427     if(_nb_of_elem==newNbOfElements)
 428       return ;
 429     T *pointer=(T*)malloc(newNbOfElements*sizeof(T));
 430     std::copy(_pointer.getConstPointer(),_pointer.getConstPointer()+std::min<std::size_t>(_nb_of_elem,newNbOfElements),pointer);
 431     if(_ownership)
 432       DestroyPointer(const_cast<T *>(_pointer.getConstPointer()),_dealloc,_param_for_deallocator);//Do not use getPointer because in case of _external
 433     _pointer.setInternal(pointer);
 434     _nb_of_elem=newNbOfElements;
 435     _nb_of_elem_alloc=newNbOfElements;
 436     _ownership=true;
 437     _dealloc=CDeallocator;
 438     _param_for_deallocator=0;
 439   }
 440
 441   template<class T>
 442   void MemArray<T>::CPPDeallocator(void *pt, void *param)
 443   {
 444     delete [] reinterpret_cast<T*>(pt);
 445   }
 446
 447   template<class T>
 448   void MemArray<T>::CDeallocator(void *pt, void *param)
 449   {
 450     free(pt);
 451   }
 452
 453   template<class T>
 454   typename MemArray<T>::Deallocator MemArray<T>::BuildFromType(DeallocType type)
 455   {
 456     switch(type)
 457     {
 458       case CPP_DEALLOC:
 459         return CPPDeallocator;
 460       case C_DEALLOC:
 461         return CDeallocator;
 462       default:
 463         throw INTERP_KERNEL::Exception("Invalid deallocation requested ! Unrecognized enum DeallocType !");
 464     }
 465   }
 466
 467   template<class T>
 468   void MemArray<T>::DestroyPointer(T *pt, typename MemArray<T>::Deallocator dealloc, void *param)
 469   {
 470     if(dealloc)
 471       dealloc(pt,param);
 472   }
 473
 474   template<class T>
 475   void MemArray<T>::destroy()
 476   {
 477     if(_ownership)
 478       DestroyPointer(const_cast<T *>(_pointer.getConstPointer()),_dealloc,_param_for_deallocator);//Do not use getPointer because in case of _external
 479     _pointer.null();
 480     _ownership=false;
 481     _dealloc=NULL;
 482     _param_for_deallocator=NULL;
 483     _nb_of_elem=0;
 484     _nb_of_elem_alloc=0;
 485   }
 486
 487   template<class T>
 488   MemArray<T> &MemArray<T>::operator=(const MemArray<T>& other)
 489   {
 490     alloc(other._nb_of_elem);
 491     std::copy(other._pointer.getConstPointer(),other._pointer.getConstPointer()+_nb_of_elem,_pointer.getPointer());
 492     return *this;
 493   }
 494
 495   //////////////////////////////////
 496
 497   template<class T>
 498   std::size_t DataArrayTemplate<T>::getHeapMemorySizeWithoutChildren() const
 499   {
 500     std::size_t sz(_mem.getNbOfElemAllocated());
 501     sz*=sizeof(T);
 502     return DataArray::getHeapMemorySizeWithoutChildren()+sz;
 503   }
 504
 505   /*!
 506    * Allocates the raw data in memory. If the memory was already allocated, then it is
 507    * freed and re-allocated. See an example of this method use
 508    * \ref MEDCouplingArraySteps1WC "here".
 509    *  \param [in] nbOfTuple - number of tuples of data to allocate.
 510    *  \param [in] nbOfCompo - number of components of data to allocate.
 511    *  \throw If \a nbOfTuple < 0 or \a nbOfCompo < 0.
 512    */
 513   template<class T>
 514   void DataArrayTemplate<T>::alloc(int nbOfTuple, int nbOfCompo)
 515   {
 516     if(nbOfTuple<0 || nbOfCompo<0)
 517       {
 518         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::alloc : request for negative length of data !";
 519         throw INTERP_KERNEL::Exception(oss.str().c_str());
 520       }
 521     _info_on_compo.resize(nbOfCompo);
 522     _mem.alloc(nbOfCompo*(std::size_t)nbOfTuple);
 523     declareAsNew();
 524   }
 525
 526   /*!
 527    * Sets a C array to be used as raw data of \a this. The previously set info
 528    *  of components is retained and re-sized.
 529    * For more info see \ref MEDCouplingArraySteps1.
 530    *  \param [in] array - the C array to be used as raw data of \a this.
 531    *  \param [in] ownership - if \a true, \a array will be deallocated at destruction of \a this.
 532    *  \param [in] type - specifies how to deallocate \a array. If \a type == MEDCoupling::CPP_DEALLOC,
 533    *                     \c delete [] \c array; will be called. If \a type == MEDCoupling::C_DEALLOC,
 534    *                     \c free(\c array ) will be called.
 535    *  \param [in] nbOfTuple - new number of tuples in \a this.
 536    *  \param [in] nbOfCompo - new number of components in \a this.
 537    */
 538   template<class T>
 539   void DataArrayTemplate<T>::useArray(const T *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo)
 540   {
 541     _info_on_compo.resize(nbOfCompo);
 542     _mem.useArray(array,ownership,type,(std::size_t)nbOfTuple*nbOfCompo);
 543     declareAsNew();
 544   }
 545
 546   template<class T>
 547   void DataArrayTemplate<T>::useExternalArrayWithRWAccess(const T *array, int nbOfTuple, int nbOfCompo)
 548   {
 549     _info_on_compo.resize(nbOfCompo);
 550     _mem.useExternalArrayWithRWAccess(array,(std::size_t)nbOfTuple*nbOfCompo);
 551     declareAsNew();
 552   }
 553
 554   /*!
 555    * Returns a value located at specified tuple and component.
 556    * This method is equivalent to DataArrayTemplate<T>::getIJ() except that validity of
 557    * parameters is checked. So this method is safe but expensive if used to go through
 558    * all values of \a this.
 559    *  \param [in] tupleId - index of tuple of interest.
 560    *  \param [in] compoId - index of component of interest.
 561    *  \return double - value located by \a tupleId and \a compoId.
 562    *  \throw If \a this is not allocated.
 563    *  \throw If condition <em>( 0 <= tupleId < this->getNumberOfTuples() )</em> is violated.
 564    *  \throw If condition <em>( 0 <= compoId < this->getNumberOfComponents() )</em> is violated.
 565    */
 566   template<class T>
 567   T DataArrayTemplate<T>::getIJSafe(int tupleId, int compoId) const
 568   {
 569     checkAllocated();
 570     if(tupleId<0 || tupleId>=getNumberOfTuples())
 571       {
 572         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::getIJSafe : request for tupleId " << tupleId << " should be in [0," << getNumberOfTuples() << ") !";
 573         throw INTERP_KERNEL::Exception(oss.str().c_str());
 574       }
 575     if(compoId<0 || compoId>=getNumberOfComponents())
 576       {
 577         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::getIJSafe : request for compoId " << compoId << " should be in [0," << getNumberOfComponents() << ") !";
 578         throw INTERP_KERNEL::Exception(oss.str().c_str());
 579       }
 580     return _mem[tupleId*_info_on_compo.size()+compoId];
 581   }
 582
 583   /*!
 584    * This method \b do \b not modify content of \a this. It only modify its memory footprint if the allocated memory is to high regarding real data to store.
 585    *
 586    * \sa DataArray::getHeapMemorySizeWithoutChildren, DataArrayTemplate<T>::reserve
 587    */
 588   template<class T>
 589   void DataArrayTemplate<T>::pack() const
 590   {
 591     _mem.pack();
 592   }
 593
 594   /*!
 595    * Checks if raw data is allocated. Read more on the raw data
 596    * in \ref MEDCouplingArrayBasicsTuplesAndCompo "DataArrays infos" for more information.
 597    *  \return bool - \a true if the raw data is allocated, \a false else.
 598    */
 599   template<class T>
 600   bool DataArrayTemplate<T>::isAllocated() const
 601   {
 602     return getConstPointer()!=0;
 603   }
 604
 605   /*!
 606    * Checks if raw data is allocated and throws an exception if it is not the case.
 607    *  \throw If the raw data is not allocated.
 608    */
 609   template<class T>
 610   void DataArrayTemplate<T>::checkAllocated() const
 611   {
 612     if(!isAllocated())
 613       {
 614         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::checkAllocated : Array is defined but not allocated ! Call alloc or setValues method first !";
 615         throw INTERP_KERNEL::Exception(oss.str().c_str());
 616       }
 617   }
 618
 619   /*!
 620    * This method desallocated \a this without modification of informations relative to the components.
 621    * After call of this method, DataArrayDouble::isAllocated will return false.
 622    * If \a this is already not allocated, \a this is let unchanged.
 623    */
 624   template<class T>
 625   void DataArrayTemplate<T>::desallocate()
 626   {
 627     _mem.destroy();
 628   }
 629
 630   /*!
 631    * This method reserve nbOfElems elements in memory ( nbOfElems*8 bytes ) \b without impacting the number of tuples in \a this.
 632    * If \a this has already been allocated, this method checks that \a this has only one component. If not an INTERP_KERNEL::Exception will be thrown.
 633    * If \a this has not already been allocated, number of components is set to one.
 634    * This method allows to reduce number of reallocations on invokation of DataArrayDouble::pushBackSilent and DataArrayDouble::pushBackValsSilent on \a this.
 635    *
 636    * \sa DataArrayDouble::pack, DataArrayDouble::pushBackSilent, DataArrayDouble::pushBackValsSilent
 637    */
 638   template<class T>
 639   void DataArrayTemplate<T>::reserve(std::size_t nbOfElems)
 640   {
 641     int nbCompo(getNumberOfComponents());
 642     if(nbCompo==1)
 643       {
 644         _mem.reserve(nbOfElems);
 645       }
 646     else if(nbCompo==0)
 647       {
 648         _mem.reserve(nbOfElems);
 649         _info_on_compo.resize(1);
 650       }
 651     else
 652       {
 653         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::reserve : not available for DataArrayDouble with number of components different than 1 !";
 654         throw INTERP_KERNEL::Exception(oss.str().c_str());
 655       }
 656   }
 657
 658   /*!
 659    * This method adds at the end of \a this the single value \a val. This method do \b not update its time label to avoid useless incrementation
 660    * of counter. So the caller is expected to call TimeLabel::declareAsNew on \a this at the end of the push session.
 661    *
 662    * \param [in] val the value to be added in \a this
 663    * \throw If \a this has already been allocated with number of components different from one.
 664    * \sa DataArrayDouble::pushBackValsSilent
 665    */
 666   template<class T>
 667   void DataArrayTemplate<T>::pushBackSilent(T val)
 668   {
 669     int nbCompo(getNumberOfComponents());
 670     if(nbCompo==1)
 671       _mem.pushBack(val);
 672     else if(nbCompo==0)
 673       {
 674         _info_on_compo.resize(1);
 675         _mem.pushBack(val);
 676       }
 677     else
 678       {
 679         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::pushBackSilent : not available for DataArrayDouble with number of components different than 1 !";
 680         throw INTERP_KERNEL::Exception(oss.str().c_str());
 681       }
 682   }
 683
 684   /*!
 685    * This method adds at the end of \a this a serie of values [\c valsBg,\c valsEnd). This method do \b not update its time label to avoid useless incrementation
 686    * of counter. So the caller is expected to call TimeLabel::declareAsNew on \a this at the end of the push session.
 687    *
 688    *  \param [in] valsBg - an array of values to push at the end of \c this.
 689    *  \param [in] valsEnd - specifies the end of the array \a valsBg, so that
 690    *              the last value of \a valsBg is \a valsEnd[ -1 ].
 691    * \throw If \a this has already been allocated with number of components different from one.
 692    * \sa DataArrayDouble::pushBackSilent
 693    */
 694   template<class T>
 695   void DataArrayTemplate<T>::pushBackValsSilent(const T *valsBg, const T *valsEnd)
 696   {
 697     int nbCompo(getNumberOfComponents());
 698     if(nbCompo==1)
 699       _mem.insertAtTheEnd(valsBg,valsEnd);
 700     else if(nbCompo==0)
 701       {
 702         _info_on_compo.resize(1);
 703         _mem.insertAtTheEnd(valsBg,valsEnd);
 704       }
 705     else
 706       {
 707         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::pushBackValsSilent : not available for DataArrayDouble with number of components different than 1 !";
 708         throw INTERP_KERNEL::Exception(oss.str().c_str());
 709       }
 710   }
 711
 712   /*!
 713    * This method returns silently ( without updating time label in \a this ) the last value, if any and suppress it.
 714    * \throw If \a this is already empty.
 715    * \throw If \a this has number of components different from one.
 716    */
 717   template<class T>
 718   T DataArrayTemplate<T>::popBackSilent()
 719   {
 720     if(getNumberOfComponents()==1)
 721       return _mem.popBack();
 722     else
 723       {
 724         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::popBackSilent : not available for DataArrayDouble with number of components different than 1 !";
 725         throw INTERP_KERNEL::Exception(oss.str().c_str());
 726       }
 727   }
 728
 729   /*!
 730    * Allocates the raw data in memory. If exactly same memory as needed already
 731    * allocated, it is not re-allocated.
 732    *  \param [in] nbOfTuple - number of tuples of data to allocate.
 733    *  \param [in] nbOfCompo - number of components of data to allocate.
 734    *  \throw If \a nbOfTuple < 0 or \a nbOfCompo < 0.
 735    */
 736   template<class T>
 737   void DataArrayTemplate<T>::allocIfNecessary(int nbOfTuple, int nbOfCompo)
 738   {
 739     if(isAllocated())
 740       {
 741         if(nbOfTuple!=getNumberOfTuples() || nbOfCompo!=getNumberOfComponents())
 742           alloc(nbOfTuple,nbOfCompo);
 743       }
 744     else
 745       alloc(nbOfTuple,nbOfCompo);
 746   }
 747
 748   /*!
 749    * Checks the number of tuples.
 750    *  \return bool - \a true if getNumberOfTuples() == 0, \a false else.
 751    *  \throw If \a this is not allocated.
 752    */
 753   template<class T>
 754   bool DataArrayTemplate<T>::empty() const
 755   {
 756     checkAllocated();
 757     return getNumberOfTuples()==0;
 758   }
 759
 760   /*!
 761    * Copies all the data from another DataArrayDouble. For more info see
 762    * \ref MEDCouplingArrayBasicsCopyDeepAssign.
 763    *  \param [in] other - another instance of DataArrayDouble to copy data from.
 764    *  \throw If the \a other is not allocated.
 765    */
 766   template<class T>
 767   void DataArrayTemplate<T>::deepCopyFrom(const DataArrayTemplate<T>& other)
 768   {
 769     other.checkAllocated();
 770     int nbOfTuples(other.getNumberOfTuples()),nbOfComp(other.getNumberOfComponents());
 771     allocIfNecessary(nbOfTuples,nbOfComp);
 772     std::size_t nbOfElems((std::size_t)nbOfTuples*nbOfComp);
 773     T *pt(getPointer());
 774     const T *ptI(other.begin());
 775     for(std::size_t i=0;i<nbOfElems;i++)
 776       pt[i]=ptI[i];
 777     copyStringInfoFrom(other);
 778   }
 779
 780   /*!
 781    * Reverse the array values.
 782    *  \throw If \a this->getNumberOfComponents() < 1.
 783    *  \throw If \a this is not allocated.
 784    */
 785   template<class T>
 786   void DataArrayTemplate<T>::reverse()
 787   {
 788     checkAllocated();
 789     _mem.reverse(getNumberOfComponents());
 790     declareAsNew();
 791   }
 792
 793   /*!
 794    * Assign \a val to all values in \a this array. To know more on filling arrays see
 795    * \ref MEDCouplingArrayFill.
 796    *  \param [in] val - the value to fill with.
 797    *  \throw If \a this is not allocated.
 798    */
 799   template<class T>
 800   void DataArrayTemplate<T>::fillWithValue(T val)
 801   {
 802     checkAllocated();
 803     _mem.fillWithValue(val);
 804     declareAsNew();
 805   }
 806
 807   /*!
 808    * Changes number of tuples in the array. If the new number of tuples is smaller
 809    * than the current number the array is truncated, otherwise the array is extended.
 810    *  \param [in] nbOfTuples - new number of tuples.
 811    *  \throw If \a this is not allocated.
 812    *  \throw If \a nbOfTuples is negative.
 813    */
 814   template<class T>
 815   void DataArrayTemplate<T>::reAlloc(int nbOfTuples)
 816   {
 817     if(nbOfTuples<0)
 818       {
 819         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::reAlloc : input new number of tuples should be >=0 !";
 820         throw INTERP_KERNEL::Exception(oss.str().c_str());
 821       }
 822     checkAllocated();
 823     _mem.reAlloc(getNumberOfComponents()*(std::size_t)nbOfTuples);
 824     declareAsNew();
 825   }
 826
 827   /*!
 828    * Permutes values of \a this array as required by \a old2New array. The values are
 829    * permuted so that \c new[ \a old2New[ i ]] = \c old[ i ]. Number of tuples remains
 830    * the same as in \c this one.
 831    * If a permutation reduction is needed, subArray() or selectByTupleId() should be used.
 832    * For more info on renumbering see \ref numbering.
 833    *  \param [in] old2New - C array of length equal to \a this->getNumberOfTuples()
 834    *     giving a new position for i-th old value.
 835    */
 836   template<class T>
 837   void DataArrayTemplate<T>::renumberInPlace(const int *old2New)
 838   {
 839     checkAllocated();
 840     int nbTuples(getNumberOfTuples()),nbOfCompo(getNumberOfComponents());
 841     T *tmp(new T[nbTuples*nbOfCompo]);
 842     const T *iptr(begin());
 843     for(int i=0;i<nbTuples;i++)
 844       {
 845         int v=old2New[i];
 846         if(v>=0 && v<nbTuples)
 847           std::copy(iptr+nbOfCompo*i,iptr+nbOfCompo*(i+1),tmp+nbOfCompo*v);
 848         else
 849           {
 850             std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::renumberInPlace : At place #" << i << " value is " << v << " ! Should be in [0," << nbTuples << ") !";
 851             throw INTERP_KERNEL::Exception(oss.str().c_str());
 852           }
 853       }
 854     std::copy(tmp,tmp+nbTuples*nbOfCompo,getPointer());
 855     delete [] tmp;
 856     declareAsNew();
 857   }
 858
 859
 860   /*!
 861    * Permutes values of \a this array as required by \a new2Old array. The values are
 862    * permuted so that \c new[ i ] = \c old[ \a new2Old[ i ]]. Number of tuples remains
 863    * the same as in \c this one.
 864    * For more info on renumbering see \ref numbering.
 865    *  \param [in] new2Old - C array of length equal to \a this->getNumberOfTuples()
 866    *     giving a previous position of i-th new value.
 867    *  \return DataArrayDouble * - the new instance of DataArrayDouble that the caller
 868    *          is to delete using decrRef() as it is no more needed.
 869    */
 870   template<class T>
 871   void DataArrayTemplate<T>::renumberInPlaceR(const int *new2Old)
 872   {
 873     checkAllocated();
 874     int nbTuples(getNumberOfTuples()),nbOfCompo(getNumberOfComponents());
 875     T *tmp(new T[nbTuples*nbOfCompo]);
 876     const T *iptr(begin());
 877     for(int i=0;i<nbTuples;i++)
 878       {
 879         int v=new2Old[i];
 880         if(v>=0 && v<nbTuples)
 881           std::copy(iptr+nbOfCompo*v,iptr+nbOfCompo*(v+1),tmp+nbOfCompo*i);
 882         else
 883           {
 884             std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::renumberInPlaceR : At place #" << i << " value is " << v << " ! Should be in [0," << nbTuples << ") !";
 885             throw INTERP_KERNEL::Exception(oss.str().c_str());
 886           }
 887       }
 888     std::copy(tmp,tmp+nbTuples*nbOfCompo,getPointer());
 889     delete [] tmp;
 890     declareAsNew();
 891   }
 892
 893   /*!
 894    * Sorts values of the array.
 895    *  \param [in] asc - \a true means ascending order, \a false, descending.
 896    *  \throw If \a this is not allocated.
 897    *  \throw If \a this->getNumberOfComponents() != 1.
 898    */
 899   template<class T>
 900   void DataArrayTemplate<T>::sort(bool asc)
 901   {
 902     checkAllocated();
 903     if(getNumberOfComponents()!=1)
 904       {
 905         std::ostringstream oss; oss << Traits<T>::ArrayTypeName << "::sort : only supported with 'this' array with ONE component !";
 906         throw INTERP_KERNEL::Exception(oss.str().c_str());
 907       }
 908     _mem.sort(asc);
 909     declareAsNew();
 910   }
 911
 912   /*!
 913    * Returns a copy of \a this array with values permuted as required by \a old2New array.
 914    * The values are permuted so that  \c new[ \a old2New[ i ]] = \c old[ i ].
 915    * Number of tuples in the result array remains the same as in \c this one.
 916    * If a permutation reduction is needed, renumberAndReduce() should be used.
 917    * For more info on renumbering see \ref numbering.
 918    *  \param [in] old2New - C array of length equal to \a this->getNumberOfTuples()
 919    *          giving a new position for i-th old value.
 920    *  \return DataArrayDouble * - the new instance of DataArrayDouble that the caller
 921    *          is to delete using decrRef() as it is no more needed.
 922    *  \throw If \a this is not allocated.
 923    */
 924   template<class T>
 925   typename Traits<T>::ArrayType *DataArrayTemplate<T>::renumber(const int *old2New) const
 926   {
 927     checkAllocated();
 928     int nbTuples(getNumberOfTuples()),nbOfCompo(getNumberOfComponents());
 929     MCAuto<DataArray> ret0(buildNewEmptyInstance());
 930     MCAuto< typename Traits<T>::ArrayType > ret(DynamicCastSafe<DataArray,typename Traits<T>::ArrayType>(ret0));
 931     ret->alloc(nbTuples,nbOfCompo);
 932     ret->copyStringInfoFrom(*this);
 933     const T *iptr(begin());
 934     T *optr(ret->getPointer());
 935     for(int i=0;i<nbTuples;i++)
 936       std::copy(iptr+nbOfCompo*i,iptr+nbOfCompo*(i+1),optr+nbOfCompo*old2New[i]);
 937     ret->copyStringInfoFrom(*this);
 938     return ret.retn();
 939   }
 940 }
 941
 942 #endif