src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx

   1 // Copyright (C) 2007-2014  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
  20
  21 #include "MEDCouplingCartesianAMRMesh.hxx"
  22 #include "MEDCouplingFieldDouble.hxx"
  23 #include "MEDCoupling1GTUMesh.hxx"
  24 #include "MEDCouplingIMesh.hxx"
  25 #include "MEDCouplingUMesh.hxx"
  26
  27 #include <limits>
  28 #include <sstream>
  29 #include <numeric>
  30
  31 using namespace ParaMEDMEM;
  32
  33 /// @cond INTERNAL
  34
  35 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithOverlap() const
  36 {
  37   return _mesh->getNumberOfCellsRecursiveWithOverlap();
  38 }
  39
  40 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithoutOverlap() const
  41 {
  42   return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
  43 }
  44
  45 int MEDCouplingCartesianAMRPatchGen::getMaxNumberOfLevelsRelativeToThis() const
  46 {
  47   return _mesh->getMaxNumberOfLevelsRelativeToThis();
  48 }
  49
  50 MEDCouplingCartesianAMRPatchGen::MEDCouplingCartesianAMRPatchGen(MEDCouplingCartesianAMRMeshGen *mesh)
  51 {
  52   if(!mesh)
  53     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen constructor : input mesh is NULL !");
  54   _mesh=mesh; _mesh->incrRef();
  55 }
  56
  57 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe() const
  58 {
  59   const MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
  60   if(!mesh)
  61     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe const : the mesh is NULL !");
  62   return mesh;
  63 }
  64
  65 MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe()
  66 {
  67   MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
  68     if(!mesh)
  69       throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe : the mesh is NULL !");
  70     return mesh;
  71 }
  72
  73 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatchGen::getDirectChildren() const
  74 {
  75   std::vector<const BigMemoryObject *> ret;
  76   if((const MEDCouplingCartesianAMRMeshGen *)_mesh)
  77     ret.push_back((const MEDCouplingCartesianAMRMeshGen *)_mesh);
  78   return ret;
  79 }
  80
  81 /*!
  82  * \param [in] mesh not null pointer of refined mesh replacing the cell range of \a father defined by the bottom left and top right just after.
  83  * \param [in] bottomLeftTopRight a vector equal to the space dimension of \a mesh that specifies for each dimension, the included cell start of the range for the first element of the pair,
  84  *                                a the end cell (\b excluded) of the range for the second element of the pair.
  85  */
  86 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMeshGen *mesh, const std::vector< std::pair<int,int> >& bottomLeftTopRight):MEDCouplingCartesianAMRPatchGen(mesh),_bl_tr(bottomLeftTopRight)
  87 {
  88   int dim((int)bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
  89   if(dim!=dimExp)
  90     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
  91 }
  92
  93 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
  94 {
  95   return getMeshSafe()->addPatch(bottomLeftTopRight,factors);
  96 }
  97
  98 int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
  99 {
 100   return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
 101 }
 102
 103 /*!
 104  * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
 105  * the must be >= 0. \b WARNING this method only works if \a this and \a other share the same father (no check of this will be done !).
 106  * Call isInMyNeighborhoodExt to deal with 2 patches not sharing directly the same father.
 107  *
 108  * \param [in] other - The other patch
 109  * \param [in] ghostLev - The size of the neighborhood zone.
 110  *
 111  * \throw if \a this or \a other are invalid (end before start).
 112  * \throw if \a ghostLev is \b not >= 0 .
 113  * \throw if \a this and \a other have not the same space dimension.
 114  *
 115  * \sa isInMyNeighborhoodExt
 116  */
 117 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
 118 {
 119   if(ghostLev<0)
 120     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the size of the neighborhood must be >= 0 !");
 121   if(!other)
 122     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the input patch is NULL !");
 123   const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
 124   const std::vector< std::pair<int,int> >& otherp(other->getBLTRRange());
 125   return IsInMyNeighborhood(ghostLev,thisp,otherp);
 126 }
 127
 128 /*!
 129  * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
 130  * the must be >= 0. This method works even if \a this and \a other does not share the same father.
 131  *
 132  * \param [in] other - The other patch
 133  * \param [in] ghostLev - The size of the neighborhood zone.
 134  *
 135  * \throw if \a this or \a other are invalid (end before start).
 136  * \throw if \a ghostLev is \b not >= 0 .
 137  * \throw if \a this and \a other have not the same space dimension.
 138  * \throw if there is not common ancestor of \a this and \a other.
 139  *
 140  * \sa isInMyNeighborhood
 141  */
 142 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
 143 {
 144   if(ghostLev<0)
 145     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the size of the neighborhood must be >= 0 !");
 146   if(!other)
 147     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhoodExt : the input patch is NULL !");
 148   int lev;
 149   const MEDCouplingCartesianAMRMeshGen *com(FindCommonAncestor(this,other,lev));//check that factors are OK
 150   if(lev==0)
 151     return isInMyNeighborhood(other,ghostLev);
 152   std::vector<int> offset(ComputeOffsetFromTwoToOne(com,lev,this,other));
 153   const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
 154   std::vector< std::pair<int,int> > otherp(other->getBLTRRange());
 155   std::size_t sz(offset.size());
 156   for(std::size_t i=0;i<sz;i++)
 157     {
 158       otherp[i].first+=offset[i];
 159       otherp[i].second+=offset[i];
 160     }
 161   return IsInMyNeighborhood(ghostLev,thisp,otherp);
 162 }
 163
 164 bool MEDCouplingCartesianAMRPatch::IsInMyNeighborhood(int ghostLev, const std::vector< std::pair<int,int> >& p1, const std::vector< std::pair<int,int> >& p2)
 165 {
 166   std::size_t thispsize(p1.size());
 167   if(thispsize!=p2.size())
 168     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the dimensions must be the same !");
 169   for(std::size_t i=0;i<thispsize;i++)
 170     {
 171       const std::pair<int,int>& thispp(p1[i]);
 172       const std::pair<int,int>& otherpp(p2[i]);
 173       if(thispp.second<thispp.first)
 174         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
 175       if(otherpp.second<otherpp.first)
 176         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
 177       if(otherpp.first==thispp.second+ghostLev-1)
 178         continue;
 179       if(otherpp.second+ghostLev-1==thispp.first)
 180         continue;
 181       int start(std::max(thispp.first,otherpp.first)),end(std::min(thispp.second,otherpp.second));
 182       if(end<start)
 183         return false;
 184     }
 185   return true;
 186 }
 187
 188 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
 189 {
 190   std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
 191   ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
 192   return ret;
 193 }
 194
 195 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatch::FindCommonAncestor(const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2, int& lev)
 196 {
 197   const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
 198   lev=0;
 199   while(f1!=f2 || f1==0 || f2==0)
 200     {
 201       f1=f1->getFather(); f2=f2->getFather();
 202       if(f1->getFactors()!=f2->getFactors())
 203         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : factors differ !");
 204       lev++;
 205     }
 206   if(f1!=f2)
 207     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::FindCommonAncestor : no common ancestor between p1 and p2 !");
 208   return f1;
 209 }
 210
 211 std::vector<int> MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne(const MEDCouplingCartesianAMRMeshGen *comAncestor, int lev, const MEDCouplingCartesianAMRPatch *p1, const MEDCouplingCartesianAMRPatch *p2)
 212 {
 213   if(lev<=0)
 214     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : this method is useful only for lev > 0 !");
 215   int dim(p1->getMesh()->getSpaceDimension());
 216   if(p2->getMesh()->getSpaceDimension()!=dim)
 217     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::ComputeOffsetFromTwoToOne : dimension must be the same !");
 218   std::vector< int > ret(dim,0);
 219
 220   for(int i=0;i<lev;i++)
 221     {
 222       const MEDCouplingCartesianAMRMeshGen *f1(p1->_mesh),*f2(p2->_mesh);
 223       const MEDCouplingCartesianAMRPatch *p1h(0),*p2h(0);
 224       for(int j=0;j<lev-i;j++)
 225         {
 226           const MEDCouplingCartesianAMRMeshGen *f1tmp(f1->getFather()),*f2tmp(f2->getFather());
 227           int pid1(f1tmp->getPatchIdFromChildMesh(f1)),pid2(f2tmp->getPatchIdFromChildMesh(f2));
 228           p1h=f1tmp->getPatch(pid1); p2h=f2tmp->getPatch(pid2);
 229           f1=f1tmp; f2=f2tmp;
 230         }
 231       std::vector< std::pair<int,int> > p2c(p2h->getBLTRRange());
 232       for(int k=0;k<dim;k++)
 233         {
 234           p2c[k].first+=ret[k];
 235           p2c[k].second+=ret[k];
 236         }
 237       for(int k=0;k<dim;k++)
 238         {
 239           ret[k]=p2c[k].first-p1h->getBLTRRange()[k].first;
 240           ret[k]*=f1->getFactors()[k];
 241         }
 242     }
 243   return ret;
 244 }
 245
 246 MEDCouplingCartesianAMRPatchGF::MEDCouplingCartesianAMRPatchGF(MEDCouplingCartesianAMRMesh *mesh):MEDCouplingCartesianAMRPatchGen(mesh)
 247 {
 248 }
 249
 250 std::size_t MEDCouplingCartesianAMRPatchGF::getHeapMemorySizeWithoutChildren() const
 251 {
 252   return sizeof(MEDCouplingCartesianAMRPatchGF);
 253 }
 254
 255 MEDCouplingDataForGodFather::MEDCouplingDataForGodFather(MEDCouplingCartesianAMRMesh *gf):_gf(gf),_tlc(gf)
 256 {
 257   if(!_gf)
 258     throw INTERP_KERNEL::Exception("MEDCouplingDataForGodFather constructor : A data has to be attached to a AMR Mesh instance !");
 259   _gf->setData(this);
 260 }
 261
 262 void MEDCouplingDataForGodFather::checkGodFatherFrozen() const
 263 {
 264   _tlc.checkConst();
 265 }
 266
 267 bool MEDCouplingDataForGodFather::changeGodFather(MEDCouplingCartesianAMRMesh *gf)
 268 {
 269   bool ret(_tlc.keepTrackOfNewTL(gf));
 270   if(ret)
 271     {
 272       _gf=gf;
 273     }
 274   return ret;
 275 }
 276
 277 /// @endcond
 278
 279 int MEDCouplingCartesianAMRMeshGen::getSpaceDimension() const
 280 {
 281   return _mesh->getSpaceDimension();
 282 }
 283
 284 void MEDCouplingCartesianAMRMeshGen::setFactors(const std::vector<int>& newFactors)
 285 {
 286   if(getSpaceDimension()!=(int)newFactors.size())
 287     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : size of input factors is not equal to the space dimension !");
 288   if(_factors.empty())
 289     {
 290       _factors=newFactors;
 291       return ;
 292     }
 293   if(_factors==newFactors)
 294     return ;
 295   if(!_patches.empty())
 296     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : modification of factors is not allowed when presence of patches !");
 297   _factors=newFactors;
 298   declareAsNew();
 299 }
 300
 301 int MEDCouplingCartesianAMRMeshGen::getMaxNumberOfLevelsRelativeToThis() const
 302 {
 303   int ret(1);
 304   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 305     ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
 306   return ret;
 307 }
 308
 309 /*!
 310  * This method returns the number of cells of \a this with the help of the MEDCouplingIMesh instance representing \a this.
 311  * The patches in \a this are ignored here.
 312  * \sa getNumberOfCellsAtCurrentLevelGhost, getNumberOfCellsRecursiveWithOverlap
 313  */
 314 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevel() const
 315 {
 316   return _mesh->getNumberOfCells();
 317 }
 318
 319 /*!
 320  * This method returns the number of cells of \a this with the help of the MEDCouplingIMesh instance representing \a this enlarged by \a ghostLev size
 321  * to take into account of the ghost cells for future computation.
 322  * The patches in \a this are ignored here.
 323  *
 324  * \sa getNumberOfCellsAtCurrentLevel
 325  */
 326 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevelGhost(int ghostLev) const
 327 {
 328   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> tmp(_mesh->buildWithGhost(ghostLev));
 329   return tmp->getNumberOfCells();
 330 }
 331
 332 /*!
 333  * This method returns the number of cells including the current level but \b also \b including recursively all cells of other levels
 334  * starting from this. The set of cells which size is returned here are generally overlapping each other.
 335  */
 336 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithOverlap() const
 337 {
 338   int ret(_mesh->getNumberOfCells());
 339   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 340     {
 341       ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
 342     }
 343   return ret;
 344 }
 345
 346 /*!
 347  * This method returns the max number of cells covering all the space without overlapping.
 348  * It returns the number of cells of the mesh with the highest resolution.
 349  * The returned value is equal to the number of cells of mesh returned by buildUnstructured.
 350  *
 351  * \sa buildUnstructured
 352  */
 353 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithoutOverlap() const
 354 {
 355   int ret(_mesh->getNumberOfCells());
 356   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 357     {
 358       ret-=(*it)->getNumberOfOverlapedCellsForFather();
 359       ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
 360     }
 361   return ret;
 362 }
 363
 364 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getFather() const
 365 {
 366   return _father;
 367 }
 368
 369 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getGodFather() const
 370 {
 371   if(_father==0)
 372     return this;
 373   else
 374     return _father->getGodFather();
 375 }
 376
 377 /*!
 378  * This method returns the level of \a this. 0 for god father. -1 for children of god father ...
 379  */
 380 int MEDCouplingCartesianAMRMeshGen::getAbsoluteLevel() const
 381 {
 382   if(_father==0)
 383     return 0;
 384   else
 385     return _father->getAbsoluteLevel()-1;
 386 }
 387
 388 /*!
 389  * This method returns grids relative to god father to specified level \a absoluteLev.
 390  *
 391  * \return std::vector<MEDCouplingCartesianAMRPatchGen *> - objects in vector are to be managed (decrRef) by the caller.
 392  */
 393 std::vector<MEDCouplingCartesianAMRPatchGen *> MEDCouplingCartesianAMRMeshGen::retrieveGridsAt(int absoluteLev) const
 394 {
 395   if(absoluteLev<0)
 396     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::retrieveGridsAt : absolute level must be >=0 !");
 397   if(_father)
 398     return getGodFather()->retrieveGridsAt(absoluteLev);
 399   //
 400   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> > rets;
 401   retrieveGridsAtInternal(absoluteLev,rets);
 402   std::vector< MEDCouplingCartesianAMRPatchGen * > ret(rets.size());
 403   for(std::size_t i=0;i<rets.size();i++)
 404     {
 405       ret[i]=rets[i].retn();
 406     }
 407   return ret;
 408 }
 409
 410 void MEDCouplingCartesianAMRMeshGen::detachFromFather()
 411 {
 412   _father=0;
 413   declareAsNew();
 414 }
 415
 416 /*!
 417  * \param [in] bottomLeftTopRight a vector equal to the space dimension of \a mesh that specifies for each dimension, the included cell start of the range for the first element of the pair,
 418  *                                a the end cell (\b excluded) of the range for the second element of the pair.
 419  * \param [in] factors The factor of refinement per axis (different from 0).
 420  */
 421 void MEDCouplingCartesianAMRMeshGen::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
 422 {
 423   checkFactorsAndIfNotSetAssign(factors);
 424   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
 425   mesh->refineWithFactor(factors);
 426   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMeshSub> zeMesh(new MEDCouplingCartesianAMRMeshSub(this,mesh));
 427   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight));
 428   _patches.push_back(elt);
 429   declareAsNew();
 430 }
 431
 432 /// @cond INTERNAL
 433
 434 class InternalPatch : public RefCountObjectOnly
 435 {
 436 public:
 437   InternalPatch():_nb_of_true(0) { }
 438   int getDimension() const { return (int)_part.size(); }
 439   double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
 440   int getNumberOfCells() const { return (int)_crit.size(); }
 441   void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
 442   std::vector<bool>& getCriterion() { return _crit; }
 443   const std::vector<bool>& getConstCriterion() const { return _crit; }
 444   void setPart(const std::vector< std::pair<int,int> >& part) { _part=part; }
 445   std::vector< std::pair<int,int> >& getPart() { return _part; }
 446   const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
 447   bool presenceOfTrue() const { return _nb_of_true>0; }
 448   std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
 449   std::vector< std::vector<int> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
 450   double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
 451   void zipToFitOnCriterion();
 452   void updateNumberOfTrue() const;
 453   MEDCouplingAutoRefCountObjectPtr<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
 454   MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
 455 protected:
 456   ~InternalPatch() { }
 457 private:
 458   mutable int _nb_of_true;
 459   std::vector<bool> _crit;
 460   //! _part is global
 461   std::vector< std::pair<int,int> > _part;
 462 };
 463
 464 void InternalPatch::zipToFitOnCriterion()
 465 {
 466   std::vector<int> cgs(computeCGS());
 467   std::vector<bool> newCrit;
 468   std::vector< std::pair<int,int> > newPart,newPart2;
 469   int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart));
 470   MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
 471   if(newNbOfTrue!=_nb_of_true)
 472     throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
 473   _crit=newCrit; _part=newPart2;
 474 }
 475
 476 void InternalPatch::updateNumberOfTrue() const
 477 {
 478   _nb_of_true=(int)std::count(_crit.begin(),_crit.end(),true);
 479 }
 480
 481 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const
 482 {
 483   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
 484   std::vector<int> cgs(computeCGS());
 485   std::vector< std::pair<int,int> > newPart;
 486   MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
 487   MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
 488   ret->setPart(partInGlobal);
 489   ret->updateNumberOfTrue();
 490   return ret;
 491 }
 492
 493 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::deepCpy() const
 494 {
 495   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
 496   (*ret)=*this;
 497   return ret;
 498 }
 499
 500 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
 501 {
 502   cutFound=false; cutPlace=-1;
 503   std::vector<double> ratio(rangeOfAxisId-1);
 504   for(int id=0;id<rangeOfAxisId-1;id++)
 505     {
 506       double efficiency[2];
 507       for(int h=0;h<2;h++)
 508         {
 509           std::vector< std::pair<int,int> > rectH(patchToBeSplit->getConstPart());
 510           if(h==0)
 511             rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id;
 512           else
 513             rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id;
 514           MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(patchToBeSplit->deepCpy());
 515           p->zipToFitOnCriterion();
 516           //anouar rectH ?
 517           efficiency[h]=p->getEfficiencyPerAxis(axisId);
 518         }
 519       ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]);
 520     }
 521   int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
 522   int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
 523   std::vector<double> ratio_inner(dimRatioInner);
 524   double minRatio(1.e10);
 525   for(int i=0; i<dimRatioInner; i++)
 526     {
 527       if(ratio[minCellDirection-1+i]<minRatio)
 528         {
 529           minRatio=ratio[minCellDirection-1+i];
 530           indexMin=i+minCellDirection;
 531         }
 532     }
 533   cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
 534 }
 535
 536 void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace)
 537 {
 538   cutPlace=-1; cutFound=false;
 539   int minCellDirection(bso.getMinCellDirection());
 540   const int dim(patchToBeSplit->getDimension());
 541   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
 542   for(int id=0;id<dim;id++)
 543     {
 544       const std::vector<int>& signature(signatures[id]);
 545       std::vector<int> hole;
 546       std::vector<double> distance;
 547       int len((int)signature.size());
 548       for(int i=0;i<len;i++)
 549         if(signature[i]==0)
 550           if(len>= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1)
 551             hole.push_back(i);
 552       if(!hole.empty())
 553         {
 554           double center(((double)len/2.));
 555           for(std::size_t i=0;i<hole.size();i++)
 556             distance.push_back(fabs(hole[i]+1.-center));
 557
 558           std::size_t posDistanceMin(std::distance(distance.begin(),std::min_element(distance.begin(),distance.end())));
 559           cutFound=true;
 560           axisId=id;
 561           cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
 562           return ;
 563         }
 564     }
 565 }
 566
 567 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
 568 {
 569   cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
 570
 571   const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
 572   int sign,minCellDirection(bso.getMinCellDirection());
 573   const int dim(patchToBeSplit->getDimension());
 574
 575   std::vector<int> zeroCrossDims(dim,-1);
 576   std::vector<int> zeroCrossVals(dim,-1);
 577   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
 578   for (int id=0;id<dim;id++)
 579     {
 580       const std::vector<int>& signature(signatures[id]);
 581
 582       std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
 583       std::vector<double> distance ;
 584
 585       for (unsigned int i=1;i<signature.size()-1;i++)
 586         derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
 587
 588       // Gradient absolute value
 589       for ( unsigned int i=1;i<derivate_second_order.size();i++)
 590         gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
 591       if(derivate_second_order.empty())
 592         continue;
 593       for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
 594         {
 595           if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
 596             sign = -1 ;
 597           if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
 598             sign = 1 ;
 599           if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
 600             sign = 0 ;
 601           if ( sign==0 || sign==-1 )
 602             if ( i >= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 )
 603               {
 604                 zero_cross.push_back(i) ;
 605                 edge.push_back(gradient_absolute[i]) ;
 606               }
 607           signe_change.push_back(sign) ;
 608         }
 609       if ( zero_cross.size() > 0 )
 610         {
 611           int max_cross=*max_element(edge.begin(),edge.end()) ;
 612           for (unsigned int i=0;i<edge.size();i++)
 613             if (edge[i]==max_cross)
 614               max_cross_list.push_back(zero_cross[i]+1) ;
 615
 616           double center((signature.size()/2.0));
 617           for (unsigned int i=0;i<max_cross_list.size();i++)
 618             distance.push_back(fabs(max_cross_list[i]+1-center));
 619
 620           float distance_min=*min_element(distance.begin(),distance.end()) ;
 621           int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
 622           int best_place = max_cross_list[pos_distance_min] + part[id].first ;
 623           if ( max_cross >=0 )
 624             {
 625               zeroCrossDims[id] = best_place ;
 626               zeroCrossVals[id] = max_cross ;
 627             }
 628         }
 629       derivate_second_order.clear() ;
 630       gradient_absolute.clear() ;
 631       signe_change.clear() ;
 632       zero_cross.clear() ;
 633       edge.clear() ;
 634       max_cross_list.clear() ;
 635       distance.clear() ;
 636     }
 637
 638   if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1  )
 639     {
 640       int max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
 641
 642       if (zeroCrossVals[0]==max_cross_dims &&  zeroCrossVals[1]==max_cross_dims )
 643         {
 644           int nl_left(part[0].second-part[0].first);
 645           int nc_left(part[1].second-part[1].first);
 646           if ( nl_left >=  nc_left )
 647             max_cross_dims = 0 ;
 648           else
 649             max_cross_dims = 1 ;
 650         }
 651       else
 652         max_cross_dims=std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin();
 653       cutFound=true;
 654       cutPlace=zeroCrossDims[max_cross_dims];
 655       axisId=max_cross_dims ;
 656     }
 657 }
 658
 659 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
 660 {
 661   cutFound=false;
 662   if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd())
 663     {
 664       if(rangeOfAxisId>=2*bso.getMinCellDirection())
 665         {
 666           cutFound=true;
 667           cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
 668         }
 669     }
 670   else
 671     {
 672       if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells())
 673         {
 674           DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace);
 675         }
 676     }
 677 }
 678
 679 MEDCouplingAutoRefCountObjectPtr<InternalPatch> DealWithNoCut(const InternalPatch *patch)
 680 {
 681   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(const_cast<InternalPatch *>(patch));
 682   ret->incrRef();
 683   return ret;
 684 }
 685
 686 void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
 687 {
 688   MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
 689   std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
 690   std::vector< std::pair<int,int> > leftRect(rect),rightRect(rect);
 691   leftRect[axisId].second=cutPlace+1;
 692   rightRect[axisId].first=cutPlace+1;
 693   leftPart=patchToBeSplit->extractPart(leftRect);
 694   rightPart=patchToBeSplit->extractPart(rightRect);
 695   leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion();
 696   listOfPatches.push_back(leftPart);
 697   listOfPatches.push_back(rightPart);
 698 }
 699
 700 /// @endcond
 701
 702 /*!
 703  * This method creates patches in \a this (by destroying the patches if any). This method uses \a criterion array as a field on cells on this level.
 704  */
 705 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector<bool>& criterion, const std::vector<int>& factors)
 706 {
 707   int nbCells(getNumberOfCellsAtCurrentLevel());
 708   if(nbCells!=(int)criterion.size())
 709     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the number of tuples of criterion array must be equal to the number of cells at the current level !");
 710   _patches.clear();
 711   std::vector<int> cgs(_mesh->getCellGridStructure());
 712   std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
 713   //
 714   MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
 715   p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,criterion,p->getCriterion(),p->getPart()));
 716   if(p->presenceOfTrue())
 717     listOfPatches.push_back(p);
 718   while(!listOfPatches.empty())
 719     {
 720       std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
 721       for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
 722         {
 723           //
 724           int axisId,rangeOfAxisId,cutPlace;
 725           bool cutFound;
 726           MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
 727           if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
 728             { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
 729           FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
 730           if(cutFound)
 731             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
 732           FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
 733           if(cutFound)
 734             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
 735           TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
 736           if(cutFound)
 737             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
 738           listOfPatchesOK.push_back(DealWithNoCut(*it));
 739         }
 740       listOfPatches=listOfPatchesTmp;
 741     }
 742   for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
 743     addPatch((*it)->getConstPart(),factors);
 744   declareAsNew();
 745 }
 746
 747 /*!
 748  * This method creates patches in \a this (by destroying the patches if any). This method uses \a criterion array as a field on cells on this level.
 749  */
 750 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
 751 {
 752   if(!criterion || !criterion->isAllocated())
 753     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
 754   std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
 755   createPatchesFromCriterion(bso,crit,factors);
 756   declareAsNew();
 757 }
 758
 759 void MEDCouplingCartesianAMRMeshGen::removeAllPatches()
 760 {
 761   _patches.clear();
 762   declareAsNew();
 763 }
 764
 765 void MEDCouplingCartesianAMRMeshGen::removePatch(int patchId)
 766 {
 767   checkPatchId(patchId);
 768   int sz((int)_patches.size()),j(0);
 769   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
 770   for(int i=0;i<sz;i++)
 771     if(i!=patchId)
 772       patches[j++]=_patches[i];
 773   (const_cast<MEDCouplingCartesianAMRMeshGen *>(_patches[patchId]->getMesh()))->detachFromFather();
 774   _patches=patches;
 775   declareAsNew();
 776 }
 777
 778 int MEDCouplingCartesianAMRMeshGen::getNumberOfPatches() const
 779 {
 780   return (int)_patches.size();
 781 }
 782
 783 int MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh(const MEDCouplingCartesianAMRMeshGen *mesh) const
 784 {
 785   int ret(0);
 786   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ret++)
 787     {
 788       if((*it)->getMesh()==mesh)
 789         return ret;
 790     }
 791   throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::getPatchIdFromChildMesh : no such a mesh in my direct progeny !");
 792 }
 793
 794 std::vector< const MEDCouplingCartesianAMRPatch *> MEDCouplingCartesianAMRMeshGen::getPatches() const
 795 {
 796   std::size_t sz(_patches.size());
 797   std::vector< const MEDCouplingCartesianAMRPatch *> ret(sz);
 798   for(std::size_t i=0;i<sz;i++)
 799     ret[i]=_patches[i];
 800   return ret;
 801 }
 802
 803 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatch(int patchId) const
 804 {
 805   checkPatchId(patchId);
 806   return _patches[patchId];
 807 }
 808
 809 /*!
 810  * This method states if patch2 (with id \a patchId2) is in the neighborhood of patch1 (with id \a patchId1).
 811  * The neighborhood size is defined by \a ghostLev in the reference of \a this ( \b not in the reference of patches !).
 812  */
 813 bool MEDCouplingCartesianAMRMeshGen::isPatchInNeighborhoodOf(int patchId1, int patchId2, int ghostLev) const
 814 {
 815   const MEDCouplingCartesianAMRPatch *p1(getPatch(patchId1)),*p2(getPatch(patchId2));
 816   return p1->isInMyNeighborhood(p2,GetGhostLevelInFineRef(ghostLev,_factors));
 817 }
 818
 819 /*!
 820  * This method creates a new cell field array on given \a patchId patch in \a this starting from a coarse cell field on \a this \a cellFieldOnThis.
 821  * This method can be seen as a fast projection from the cell field \a cellFieldOnThis on \c this->getImageMesh() to a refined part of \a this
 822  * defined by the patch with id \a patchId.
 823  *
 824  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 825  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 826  * \return DataArrayDouble * - The array of the cell field on the requested patch
 827  *
 828  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
 829  * \throw if \a cellFieldOnThis is NULL or not allocated
 830  * \sa fillCellFieldOnPatch, MEDCouplingIMesh::SpreadCoarseToFine
 831  */
 832 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::createCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis) const
 833 {
 834   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
 835     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
 836   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 837   const MEDCouplingIMesh *fine(patch->getMesh()->getImageMesh());
 838   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(fine->getNumberOfCells(),cellFieldOnThis->getNumberOfComponents());
 839   ret->copyStringInfoFrom(*cellFieldOnThis);
 840   MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),ret,patch->getBLTRRange(),getFactors());
 841   return ret.retn();
 842 }
 843
 844 /*!
 845  * This method is equivalent to MEDCouplingCartesianAMRMesh::createCellFieldOnPatch except that here instead of creating a new instance
 846  * it fills the value into the \a cellFieldOnPatch data.
 847  *
 848  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 849  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 850  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
 851  *
 852  * \sa createCellFieldOnPatch, fillCellFieldComingFromPatch
 853  */
 854 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch) const
 855 {
 856   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
 857     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
 858   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 859   MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors());
 860 }
 861
 862 /*!
 863  * This method is the generalization of fillCellFieldOnPatch method. This method only projects coarse to fine without considering the
 864  * potential neighbor patches covered by the ghost cells of patch with id \a patchId.
 865  *
 866  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 867  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 868  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
 869  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
 870  *
 871  * \sa fillCellFieldOnPatch, fillCellFieldOnPatchGhostAdv
 872  */
 873 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhost(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev) const
 874 {
 875   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
 876     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatchGhost : the input cell field array is NULL or not allocated !");
 877   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 878   MEDCouplingIMesh::SpreadCoarseToFineGhost(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
 879 }
 880
 881 /*!
 882  * This method is equivalent to  fillCellFieldOnPatchGhost except that here \b ONLY \b the \b ghost \b zone will be updated
 883  * in \a cellFieldOnPatch.
 884  *
 885  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 886  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 887  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled \b only \b in \b the \b ghost \b zone.
 888  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
 889  */
 890 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyOnGhostZone(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev) const
 891 {
 892   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
 893     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyOnGhostZone : the input cell field array is NULL or not allocated !");
 894   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 895   MEDCouplingIMesh::SpreadCoarseToFineGhostZone(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
 896 }
 897
 898 /*!
 899  * This method is a refinement of fillCellFieldOnPatchGhost. fillCellFieldOnPatchGhost is first called.
 900  * Then for all other patches than those pointed by \a patchId that overlap the ghost zone of the patch impact the ghost zone adequately.
 901  *
 902  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 903  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 904  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
 905  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
 906  * \param [in] arrsOnPatches - \b WARNING arrsOnPatches[patchId] is \b NOT \b const. All others are const.
 907  *
 908  * \sa fillCellFieldOnPatchOnlyGhostAdv
 909  */
 910 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhostAdv(int patchId, const DataArrayDouble *cellFieldOnThis, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches) const
 911 {
 912   int nbp(getNumberOfPatches());
 913   if(nbp!=(int)arrsOnPatches.size())
 914     {
 915       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
 916       throw INTERP_KERNEL::Exception(oss.str().c_str());
 917     }
 918   DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
 919   // first, do as usual
 920   fillCellFieldOnPatchGhost(patchId,cellFieldOnThis,theFieldToFill,ghostLev);
 921   fillCellFieldOnPatchOnlyGhostAdv(patchId,ghostLev,arrsOnPatches);
 922 }
 923
 924 /*!
 925  * This method updates the patch with id \a patchId considering the only the all the patches in \a this to fill ghost zone.
 926  * So \b warning, the DataArrayDouble instance \a arrsOnPatches[patchId] is non const.
 927  *
 928  * \sa getPatchIdsInTheNeighborhoodOf
 929  */
 930 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchOnlyGhostAdv(int patchId, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches) const
 931 {
 932   int nbp(getNumberOfPatches()),dim(getSpaceDimension());
 933   if(nbp!=(int)arrsOnPatches.size())
 934     {
 935       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchOnlyGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
 936       throw INTERP_KERNEL::Exception(oss.str().c_str());
 937     }
 938   DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
 939   const MEDCouplingCartesianAMRPatch *refP(getPatch(patchId));
 940   const std::vector< std::pair<int,int> >& refBLTR(refP->getBLTRRange());//[(1,4),(2,4)]
 941   std::vector<int> dimsCoarse(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(refBLTR));//[3,2]
 942   std::transform(dimsCoarse.begin(),dimsCoarse.end(),_factors.begin(),dimsCoarse.begin(),std::multiplies<int>());//[12,8]
 943   std::transform(dimsCoarse.begin(),dimsCoarse.end(),dimsCoarse.begin(),std::bind2nd(std::plus<int>(),2*ghostLev));//[14,10]
 944   std::vector< std::pair<int,int> > rangeCoarse(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsCoarse));//[(0,14),(0,10)]
 945   std::vector<int> fakeFactors(dim,1),ids(getPatchIdsInTheNeighborhoodOf(patchId,ghostLev));
 946   //
 947   for(std::vector<int>::const_iterator it=ids.begin();it!=ids.end();it++)
 948     {
 949       const MEDCouplingCartesianAMRPatch *otherP(getPatch(*it));
 950       const std::vector< std::pair<int,int> >& otherBLTR(otherP->getBLTRRange());//[(4,5),(3,4)]
 951       std::vector< std::pair<int,int> > tmp0,tmp1,tmp2;
 952       MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(refBLTR,otherBLTR,tmp0,false);//tmp0=[(3,4),(1,2)]
 953       ApplyFactorsOnCompactFrmt(tmp0,_factors);//tmp0=[(12,16),(4,8)]
 954       ApplyGhostOnCompactFrmt(tmp0,ghostLev);//tmp0=[(13,17),(5,9)]
 955       std::vector< std::pair<int,int> > interstRange(MEDCouplingStructuredMesh::IntersectRanges(tmp0,rangeCoarse));//interstRange=[(13,14),(5,9)]
 956       MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(otherBLTR,refBLTR,tmp1,false);//tmp1=[(-3,0),(-1,1)]
 957       ApplyFactorsOnCompactFrmt(tmp1,_factors);//tmp1=[(-12,-4),(-4,0)]
 958       MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(tmp1,interstRange,tmp2,false);//tmp2=[(1,2),(1,5)]
 959       //
 960       std::vector< std::pair<int,int> > dimsFine(otherBLTR);
 961       ApplyFactorsOnCompactFrmt(dimsFine,_factors);
 962       ApplyAllGhostOnCompactFrmt(dimsFine,ghostLev);
 963       //
 964       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ghostVals(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(dimsFine),arrsOnPatches[*it],tmp2));
 965       MEDCouplingIMesh::CondenseFineToCoarse(dimsCoarse,ghostVals,interstRange,fakeFactors,theFieldToFill);
 966     }
 967 }
 968
 969 /*!
 970  * This method updates \a cellFieldOnThis part of values coming from the cell field \a cellFieldOnPatch lying on patch having id \a patchId.
 971  *
 972  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 973  * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
 974  * \param [in,out] cellFieldOnThis The array of the cell field on \a this to be updated only on the part concerning the patch with id \a patchId.
 975  *
 976  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
 977  * \throw if \a cellFieldOnPatch is NULL or not allocated
 978  * \sa createCellFieldOnPatch, MEDCouplingIMesh::CondenseFineToCoarse,fillCellFieldComingFromPatchGhost
 979  */
 980 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatch(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis) const
 981 {
 982   if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
 983       throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch : the input cell field array is NULL or not allocated !");
 984   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 985   MEDCouplingIMesh::CondenseFineToCoarse(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis);
 986 }
 987
 988 /*!
 989  * This method is the extension of MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch managing the ghost cells. If this
 990  * method is called with \a ghostLev equal to 0 it behaves exactly as MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch.
 991  *
 992  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 993  * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
 994  * \param [in,out] cellFieldOnThis The array of the cell field on \a this to be updated only on the part concerning the patch with id \a patchId.
 995  * \param [in] ghostLev The size of ghost zone (must be >= 0 !)
 996  *
 997  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
 998  * \throw if \a cellFieldOnPatch is NULL or not allocated
 999  * \sa fillCellFieldComingFromPatch
1000  */
1001 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatchGhost(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, int ghostLev) const
1002 {
1003   if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
1004     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatchGhost : the input cell field array is NULL or not allocated !");
1005   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
1006   MEDCouplingIMesh::CondenseFineToCoarseGhost(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis,ghostLev);
1007 }
1008
1009 /*!
1010  * This method finds all patches (located by their ids) that are in the neighborhood of patch with id \a patchId. The neighborhood size is
1011  * defined by ghostLev.
1012  *
1013  * \param [in] patchId - the id of the considered patch.
1014  * \param [in] ghostLev - the size of the neighborhood.
1015  * \return DataArrayInt * - the newly allocated array containing the list of patches in the neighborhood of the considered patch. This array is to be deallocated by the caller.
1016  */
1017 DataArrayInt *MEDCouplingCartesianAMRMeshGen::findPatchesInTheNeighborhoodOf(int patchId, int ghostLev) const
1018 {
1019   int nbp(getNumberOfPatches());
1020   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
1021   for(int i=0;i<nbp;i++)
1022     {
1023       if(i!=patchId)
1024         if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1025           ret->pushBackSilent(i);
1026     }
1027   return ret.retn();
1028 }
1029
1030 MEDCouplingUMesh *MEDCouplingCartesianAMRMeshGen::buildUnstructured() const
1031 {
1032   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1033   std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1034   std::vector<int> cgs(_mesh->getCellGridStructure());
1035   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
1036   std::size_t ii(0);
1037   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1038     {
1039       MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1040       msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
1041     }
1042   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
1043   msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1044   std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
1045   for(std::size_t i=0;i<msSafe.size();i++)
1046     ms[i]=msSafe[i];
1047   return MEDCouplingUMesh::MergeUMeshes(ms);
1048 }
1049
1050 /*!
1051  * This method returns a mesh containing as cells that there is patches at the current level.
1052  * The patches are seen like 'boxes' that is too say the refinement will not appear here.
1053  *
1054  * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
1055  */
1056 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshFromPatchEnvelop() const
1057 {
1058   std::vector<const MEDCoupling1SGTUMesh *> cells;
1059   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > cellsSafe;
1060   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1061     {
1062       const MEDCouplingCartesianAMRPatch *patch(*it);
1063       if(patch)
1064         {
1065           MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
1066           MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
1067           cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
1068         }
1069     }
1070   return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
1071 }
1072
1073 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshOfDirectChildrenOnly() const
1074 {
1075   std::vector<const MEDCoupling1SGTUMesh *> patches;
1076   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > patchesSafe;
1077   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1078       {
1079         const MEDCouplingCartesianAMRPatch *patch(*it);
1080         if(patch)
1081           {
1082             MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> patchMesh(patch->getMesh()->getImageMesh()->build1SGTUnstructured());
1083             patchesSafe.push_back(patchMesh); patches.push_back(patchMesh);
1084           }
1085       }
1086     return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(patches);
1087 }
1088
1089 /*!
1090  * This method works same as buildUnstructured except that arrays are given in input to build a field on cell in output.
1091  * \return MEDCouplingFieldDouble * - a newly created instance the caller has reponsability to deal with.
1092  * \sa buildUnstructured
1093  */
1094 MEDCouplingFieldDouble *MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost(int ghostSz, const std::vector<const DataArrayDouble *>& recurseArrs) const
1095 {
1096   if(recurseArrs.empty())
1097     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::buildCellFieldOnRecurseWithoutOverlapWithoutGhost : array is empty ! Should never happen !");
1098   //
1099   std::vector<bool> bs(_mesh->getNumberOfCells(),false);
1100   std::vector<int> cgs(_mesh->getCellGridStructure());
1101   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> > msSafe(_patches.size()+1);
1102   std::size_t ii(0);
1103   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
1104     {
1105       MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
1106       std::vector<const DataArrayDouble *> tmpArrs(extractSubTreeFromGlobalFlatten((*it)->getMesh(),recurseArrs));
1107       msSafe[ii+1]=(*it)->getMesh()->buildCellFieldOnRecurseWithoutOverlapWithoutGhost(ghostSz,tmpArrs);
1108     }
1109   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
1110   //
1111   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_CELLS));
1112   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2(extractGhostFrom(ghostSz,recurseArrs[0]));
1113   arr2=arr2->selectByTupleIdSafe(eltsOff->begin(),eltsOff->end());
1114   ret->setArray(arr2);
1115   ret->setName(arr2->getName());
1116   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
1117   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> mesh(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
1118   ret->setMesh(mesh);
1119   msSafe[0]=ret;
1120   //
1121   std::vector< const MEDCouplingFieldDouble * > ms(msSafe.size());
1122   for(std::size_t i=0;i<msSafe.size();i++)
1123     ms[i]=msSafe[i];
1124   //
1125   return MEDCouplingFieldDouble::MergeFields(ms);
1126 }
1127
1128 /*!
1129  * This method extracts from \arr arr the part inside \a arr by cutting the \a ghostSz external part.
1130  * \arr is expected to be an array having a number of tuples equal to \c getImageMesh()->buildWithGhost(ghostSz).
1131  */
1132 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::extractGhostFrom(int ghostSz, const DataArrayDouble *arr) const
1133 {
1134   std::vector<int> st(_mesh->getCellGridStructure());
1135   std::vector< std::pair<int,int> > p(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(st));
1136   std::transform(st.begin(),st.end(),st.begin(),std::bind2nd(std::plus<int>(),2*ghostSz));
1137   ApplyGhostOnCompactFrmt(p,ghostSz);
1138   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(st,arr,p));
1139   return ret.retn();
1140 }
1141
1142 /*!
1143  * This method returns all the patches in \a this not equal to \a patchId that are in neighborhood of patch with id \a patchId.
1144  *
1145  * \sa fillCellFieldOnPatchOnlyGhostAdv
1146  */
1147 std::vector<int> MEDCouplingCartesianAMRMeshGen::getPatchIdsInTheNeighborhoodOf(int patchId, int ghostLev) const
1148 {
1149   std::vector<int> ret;
1150   int nbp(getNumberOfPatches());
1151   //
1152   for(int i=0;i<nbp;i++)
1153     {
1154       if(i!=patchId)
1155         if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
1156           ret.push_back(i);
1157     }
1158   return ret;
1159 }
1160
1161 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1162                                                                const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
1163 {
1164   _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1165 }
1166
1167 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):_father(father)
1168 {
1169   if(!_father)
1170     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : empty father !");
1171   if(!mesh)
1172     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
1173   mesh->checkCoherency();
1174   _mesh=mesh; _mesh->incrRef();
1175 }
1176
1177 void MEDCouplingCartesianAMRMeshGen::checkPatchId(int patchId) const
1178 {
1179   int sz(getNumberOfPatches());
1180   if(patchId<0 || patchId>=sz)
1181     {
1182       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMeshGen::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
1183       throw INTERP_KERNEL::Exception(oss.str().c_str());
1184     }
1185 }
1186
1187 void MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign(const std::vector<int>& factors)
1188 {
1189   if(getSpaceDimension()!=(int)factors.size())
1190     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign : invalid size of factors ! size must be equal to the spaceDimension !");
1191   if(_factors.empty())
1192     {
1193       _factors=factors;
1194     }
1195   else
1196     {
1197       if(_factors!=factors)
1198         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : the factors ");
1199     }
1200 }
1201
1202 void MEDCouplingCartesianAMRMeshGen::retrieveGridsAtInternal(int lev, std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> >& grids) const
1203 {
1204   if(lev==0)
1205     {
1206       const MEDCouplingCartesianAMRMesh *thisc(dynamic_cast<const MEDCouplingCartesianAMRMesh *>(this));//tony
1207       MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGF> elt(new MEDCouplingCartesianAMRPatchGF(const_cast<MEDCouplingCartesianAMRMesh *>(thisc)));
1208       grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatchGF,MEDCouplingCartesianAMRPatchGen>(elt));
1209     }
1210   else if(lev==1)
1211     {
1212       for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1213         {
1214           const MEDCouplingCartesianAMRPatch *pt(*it);
1215           if(pt)
1216             {
1217               MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> tmp1(*it);
1218               grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatch,MEDCouplingCartesianAMRPatchGen>(tmp1));
1219             }
1220         }
1221     }
1222   else
1223     {
1224       for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1225         {
1226           const MEDCouplingCartesianAMRPatch *pt(*it);
1227           if(pt)
1228             pt->getMesh()->retrieveGridsAtInternal(lev-1,grids);
1229         }
1230     }
1231 }
1232
1233 /*!
1234  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in refined reference.
1235  * \param [in] factors - the factors per axis.
1236  */
1237 void MEDCouplingCartesianAMRMeshGen::ApplyFactorsOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, const std::vector<int>& factors)
1238 {
1239   std::size_t sz(factors.size());
1240   if(sz!=partBeforeFact.size())
1241     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::ApplyFactorsOnCompactFrmt : size of input vectors must be the same !");
1242   for(std::size_t i=0;i<sz;i++)
1243     {
1244       partBeforeFact[i].first*=factors[i];
1245       partBeforeFact[i].second*=factors[i];
1246     }
1247 }
1248
1249 /*!
1250  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
1251  * \param [in] ghostSize - the ghost size of zone for all axis.
1252  */
1253 void MEDCouplingCartesianAMRMeshGen::ApplyGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
1254 {
1255   if(ghostSize<0)
1256     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::ApplyGhostOnCompactFrmt : ghost size must be >= 0 !");
1257   std::size_t sz(partBeforeFact.size());
1258   for(std::size_t i=0;i<sz;i++)
1259     {
1260       partBeforeFact[i].first+=ghostSize;
1261       partBeforeFact[i].second+=ghostSize;
1262     }
1263 }
1264
1265 /*!
1266  * This method is different than ApplyGhostOnCompactFrmt. The \a partBeforeFact parameter is enlarger contrary to ApplyGhostOnCompactFrmt.
1267  *
1268  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
1269  * \param [in] ghostSize - the ghost size of zone for all axis.
1270  */
1271 void MEDCouplingCartesianAMRMeshGen::ApplyAllGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
1272 {
1273   if(ghostSize<0)
1274     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::ApplyAllGhostOnCompactFrmt : ghost size must be >= 0 !");
1275   std::size_t sz(partBeforeFact.size());
1276   for(std::size_t i=0;i<sz;i++)
1277     {
1278       partBeforeFact[i].first-=ghostSize;
1279       partBeforeFact[i].second+=ghostSize;
1280     }
1281 }
1282
1283 int MEDCouplingCartesianAMRMeshGen::GetGhostLevelInFineRef(int ghostLev, const std::vector<int>& factors)
1284 {
1285   if(ghostLev<0)
1286     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : the ghost size must be >=0 !");
1287   if(factors.empty())
1288     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::GetGhostLevelInFineRef : no factors defined !");
1289   int ghostLevInPatchRef;
1290   if(ghostLev==0)
1291     ghostLevInPatchRef=0;
1292   else
1293     {
1294       ghostLevInPatchRef=(ghostLev-1)/factors[0]+1;
1295       for(std::size_t i=0;i<factors.size();i++)
1296         ghostLevInPatchRef=std::max(ghostLevInPatchRef,(ghostLev-1)/factors[i]+1);
1297     }
1298   return ghostLevInPatchRef;
1299 }
1300
1301 /*!
1302  * This method returns a sub set of \a all. The subset is defined by the \a head in the tree defined by \a this.
1303  * Elements in \a all are expected to be sorted from god father to most refined structure.
1304  */
1305 std::vector<const DataArrayDouble *> MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten(const MEDCouplingCartesianAMRMeshGen *head, const std::vector<const DataArrayDouble *>& all) const
1306 {
1307   int maxLev(getMaxNumberOfLevelsRelativeToThis());
1308   std::vector<const DataArrayDouble *> ret;
1309   std::vector<const MEDCouplingCartesianAMRMeshGen *> meshes(1,this);
1310   std::size_t kk(0);
1311   for(int i=0;i<maxLev;i++)
1312     {
1313       std::vector<const MEDCouplingCartesianAMRMeshGen *> meshesTmp;
1314       for(std::vector<const MEDCouplingCartesianAMRMeshGen *>::const_iterator it=meshes.begin();it!=meshes.end();it++)
1315         {
1316           if((*it)==head || head->isObjectInTheProgeny(*it))
1317             ret.push_back(all[kk]);
1318           kk++;
1319           std::vector< const MEDCouplingCartesianAMRPatch *> ps((*it)->getPatches());
1320           for(std::vector< const MEDCouplingCartesianAMRPatch *>::const_iterator it0=ps.begin();it0!=ps.end();it0++)
1321             {
1322               const MEDCouplingCartesianAMRMeshGen *mesh((*it0)->getMesh());
1323               meshesTmp.push_back(mesh);
1324             }
1325         }
1326       meshes=meshesTmp;
1327     }
1328   if(kk!=all.size())
1329     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::extractSubTreeFromGlobalFlatten : the size of input vector is not compatible with number of leaves in this !");
1330   return ret;
1331 }
1332
1333 std::size_t MEDCouplingCartesianAMRMeshGen::getHeapMemorySizeWithoutChildren() const
1334 {
1335   return sizeof(MEDCouplingCartesianAMRMeshGen);
1336 }
1337
1338 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMeshGen::getDirectChildren() const
1339 {
1340   std::vector<const BigMemoryObject *> ret;
1341   if((const MEDCouplingIMesh *)_mesh)
1342     ret.push_back((const MEDCouplingIMesh *)_mesh);
1343   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1344     {
1345       if((const MEDCouplingCartesianAMRPatch*)*it)
1346         ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
1347     }
1348   return ret;
1349 }
1350
1351 void MEDCouplingCartesianAMRMeshGen::updateTime() const
1352 {
1353   if((const MEDCouplingIMesh *)_mesh)
1354     updateTimeWith(*_mesh);
1355   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1356     {
1357       const MEDCouplingCartesianAMRPatch *elt(*it);
1358       if(!elt)
1359         continue;
1360       const MEDCouplingCartesianAMRMeshGen *mesh(elt->getMesh());
1361       if(mesh)
1362         updateTimeWith(*mesh);
1363     }
1364 }
1365
1366 MEDCouplingCartesianAMRMeshSub::MEDCouplingCartesianAMRMeshSub(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):MEDCouplingCartesianAMRMeshGen(father,mesh)
1367 {
1368 }
1369
1370 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1371                                                               const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
1372 {
1373   return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1374 }
1375
1376 void MEDCouplingCartesianAMRMesh::setData(MEDCouplingDataForGodFather *data)
1377 {
1378   MEDCouplingDataForGodFather *myData(_data);
1379   if(myData==data)
1380     return ;
1381   if(myData)
1382     myData->changeGodFather(0);
1383   _data=data;
1384   if(data)
1385     data->incrRef();
1386 }
1387
1388 void MEDCouplingCartesianAMRMesh::allocData(int ghostLev) const
1389 {
1390   checkData();
1391   _data->alloc(ghostLev);
1392 }
1393
1394 void MEDCouplingCartesianAMRMesh::deallocData() const
1395 {
1396   checkData();
1397   _data->dealloc();
1398 }
1399
1400 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1401                                                          const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):MEDCouplingCartesianAMRMeshGen(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop)
1402 {
1403 }
1404
1405 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
1406 {
1407   std::vector<const BigMemoryObject *> ret(MEDCouplingCartesianAMRMeshGen::getDirectChildren());
1408   const MEDCouplingDataForGodFather *pt(_data);
1409   if(pt)
1410     ret.push_back(pt);
1411   return ret;
1412 }
1413
1414 void MEDCouplingCartesianAMRMesh::checkData() const
1415 {
1416   const MEDCouplingDataForGodFather *data(_data);
1417   if(!data)
1418     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkData : no data set !");
1419 }
1420
1421 MEDCouplingCartesianAMRMesh::~MEDCouplingCartesianAMRMesh()
1422 {
1423   MEDCouplingDataForGodFather *data(_data);
1424   if(data)
1425     data->changeGodFather(0);
1426 }