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 "MEDCoupling1GTUMesh.hxx"
  23 #include "MEDCouplingIMesh.hxx"
  24 #include "MEDCouplingUMesh.hxx"
  25
  26 #include <limits>
  27 #include <sstream>
  28 #include <numeric>
  29
  30 using namespace ParaMEDMEM;
  31
  32 /// @cond INTERNAL
  33
  34 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithOverlap() const
  35 {
  36   return _mesh->getNumberOfCellsRecursiveWithOverlap();
  37 }
  38
  39 int MEDCouplingCartesianAMRPatchGen::getNumberOfCellsRecursiveWithoutOverlap() const
  40 {
  41   return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
  42 }
  43
  44 int MEDCouplingCartesianAMRPatchGen::getMaxNumberOfLevelsRelativeToThis() const
  45 {
  46   return _mesh->getMaxNumberOfLevelsRelativeToThis();
  47 }
  48
  49 MEDCouplingCartesianAMRPatchGen::MEDCouplingCartesianAMRPatchGen(MEDCouplingCartesianAMRMeshGen *mesh)
  50 {
  51   if(!mesh)
  52     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen constructor : input mesh is NULL !");
  53   _mesh=mesh; _mesh->incrRef();
  54 }
  55
  56 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe() const
  57 {
  58   const MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
  59   if(!mesh)
  60     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe const : the mesh is NULL !");
  61   return mesh;
  62 }
  63
  64 MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRPatchGen::getMeshSafe()
  65 {
  66   MEDCouplingCartesianAMRMeshGen *mesh(_mesh);
  67     if(!mesh)
  68       throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatchGen::getMeshSafe : the mesh is NULL !");
  69     return mesh;
  70 }
  71
  72 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatchGen::getDirectChildren() const
  73 {
  74   std::vector<const BigMemoryObject *> ret;
  75   if((const MEDCouplingCartesianAMRMeshGen *)_mesh)
  76     ret.push_back((const MEDCouplingCartesianAMRMeshGen *)_mesh);
  77   return ret;
  78 }
  79
  80 /*!
  81  * \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.
  82  * \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,
  83  *                                a the end cell (\b excluded) of the range for the second element of the pair.
  84  */
  85 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMeshGen *mesh, const std::vector< std::pair<int,int> >& bottomLeftTopRight):MEDCouplingCartesianAMRPatchGen(mesh),_bl_tr(bottomLeftTopRight)
  86 {
  87   int dim((int)bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
  88   if(dim!=dimExp)
  89     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
  90 }
  91
  92 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
  93 {
  94   return getMeshSafe()->addPatch(bottomLeftTopRight,factors);
  95 }
  96
  97 int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
  98 {
  99   return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
 100 }
 101
 102 /*!
 103  * This method states if \a other patch is in the neighborhood of \a this. The neighborhood zone is defined by \a ghostLev parameter
 104  * the must be >= 0.
 105  *
 106  * \param [in] other - The other patch
 107  * \param [in] ghostLev - The size of the neighborhood zone.
 108  *
 109  * \throw if \a this or \a other are invalid (end before start).
 110  * \throw if \a ghostLev is \b not >= 0 .
 111  * \throw if \a this and \a other have not the same space dimension.
 112  */
 113 bool MEDCouplingCartesianAMRPatch::isInMyNeighborhood(const MEDCouplingCartesianAMRPatch *other, int ghostLev) const
 114 {
 115   if(ghostLev<0)
 116     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the size of the neighborhood must be >= 0 !");
 117   if(!other)
 118     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the input patch is NULL !");
 119   const std::vector< std::pair<int,int> >& thisp(getBLTRRange());
 120   const std::vector< std::pair<int,int> >& otherp(other->getBLTRRange());
 121   std::size_t thispsize(thisp.size());
 122   if(thispsize!=otherp.size())
 123     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : the dimensions must be the same !");
 124   for(std::size_t i=0;i<thispsize;i++)
 125     {
 126       const std::pair<int,int>& thispp(thisp[i]);
 127       const std::pair<int,int>& otherpp(otherp[i]);
 128       if(thispp.second<thispp.first)
 129         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
 130       if(otherpp.second<otherpp.first)
 131         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch::isInMyNeighborhood : this patch is invalid !");
 132       if(otherpp.first==thispp.second+ghostLev-1)
 133         continue;
 134       if(otherpp.second+ghostLev-1==thispp.first)
 135         continue;
 136       int start(std::max(thispp.first,otherpp.first)),end(std::min(thispp.second,otherpp.second));
 137       if(end<start)
 138         return false;
 139     }
 140   return true;
 141 }
 142
 143 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
 144 {
 145   std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
 146   ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
 147   return ret;
 148 }
 149
 150 MEDCouplingCartesianAMRPatchGF::MEDCouplingCartesianAMRPatchGF(MEDCouplingCartesianAMRMesh *mesh):MEDCouplingCartesianAMRPatchGen(mesh)
 151 {
 152 }
 153
 154 std::size_t MEDCouplingCartesianAMRPatchGF::getHeapMemorySizeWithoutChildren() const
 155 {
 156   return sizeof(MEDCouplingCartesianAMRPatchGF);
 157 }
 158
 159 /// @endcond
 160
 161 int MEDCouplingCartesianAMRMeshGen::getSpaceDimension() const
 162 {
 163   return _mesh->getSpaceDimension();
 164 }
 165
 166 void MEDCouplingCartesianAMRMeshGen::setFactors(const std::vector<int>& newFactors)
 167 {
 168   if(getSpaceDimension()!=(int)newFactors.size())
 169     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : size of input factors is not equal to the space dimension !");
 170   if(_factors.empty())
 171     {
 172       _factors=newFactors;
 173       return ;
 174     }
 175   if(_factors==newFactors)
 176     return ;
 177   if(!_patches.empty())
 178     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::setFactors : modification of factors is not allowed when presence of patches !");
 179   _factors=newFactors;
 180 }
 181
 182 int MEDCouplingCartesianAMRMeshGen::getMaxNumberOfLevelsRelativeToThis() const
 183 {
 184   int ret(1);
 185   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 186     ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
 187   return ret;
 188 }
 189
 190 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsAtCurrentLevel() const
 191 {
 192   return _mesh->getNumberOfCells();
 193 }
 194
 195 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithOverlap() const
 196 {
 197   int ret(_mesh->getNumberOfCells());
 198   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 199     {
 200       ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
 201     }
 202   return ret;
 203 }
 204
 205 int MEDCouplingCartesianAMRMeshGen::getNumberOfCellsRecursiveWithoutOverlap() const
 206 {
 207   int ret(_mesh->getNumberOfCells());
 208   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 209     {
 210       ret-=(*it)->getNumberOfOverlapedCellsForFather();
 211       ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
 212     }
 213   return ret;
 214 }
 215
 216 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getFather() const
 217 {
 218   return _father;
 219 }
 220
 221 const MEDCouplingCartesianAMRMeshGen *MEDCouplingCartesianAMRMeshGen::getGodFather() const
 222 {
 223   if(_father==0)
 224     return this;
 225   else
 226     return _father->getGodFather();
 227 }
 228
 229 /*!
 230  * This method returns the level of \a this. 0 for god father. -1 for children of god father ...
 231  */
 232 int MEDCouplingCartesianAMRMeshGen::getAbsoluteLevel() const
 233 {
 234   if(_father==0)
 235     return 0;
 236   else
 237     return _father->getAbsoluteLevel()-1;
 238 }
 239
 240 /*!
 241  * This method returns grids relative to god father to specified level \a absoluteLev.
 242  *
 243  * \return std::vector<MEDCouplingCartesianAMRPatchGen *> - objects in vector are to be managed (decrRef) by the caller.
 244  */
 245 std::vector<MEDCouplingCartesianAMRPatchGen *> MEDCouplingCartesianAMRMeshGen::retrieveGridsAt(int absoluteLev) const
 246 {
 247   if(absoluteLev<0)
 248     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::retrieveGridsAt : absolute level must be >=0 !");
 249   if(_father)
 250     return getGodFather()->retrieveGridsAt(absoluteLev);
 251   //
 252   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> > rets;
 253   retrieveGridsAtInternal(absoluteLev,rets);
 254   std::vector< MEDCouplingCartesianAMRPatchGen * > ret(rets.size());
 255   for(std::size_t i=0;i<rets.size();i++)
 256     {
 257       ret[i]=rets[i].retn();
 258     }
 259   return ret;
 260 }
 261
 262 void MEDCouplingCartesianAMRMeshGen::detachFromFather()
 263 {
 264   _father=0;
 265 }
 266
 267 /*!
 268  * \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,
 269  *                                a the end cell (\b excluded) of the range for the second element of the pair.
 270  * \param [in] factors The factor of refinement per axis (different from 0).
 271  */
 272 void MEDCouplingCartesianAMRMeshGen::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
 273 {
 274   checkFactorsAndIfNotSetAssign(factors);
 275   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
 276   mesh->refineWithFactor(factors);
 277   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMeshSub> zeMesh(new MEDCouplingCartesianAMRMeshSub(this,mesh));
 278   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight));
 279   _patches.push_back(elt);
 280 }
 281
 282 /// @cond INTERNAL
 283
 284 class InternalPatch : public RefCountObjectOnly
 285 {
 286 public:
 287   InternalPatch():_nb_of_true(0) { }
 288   int getDimension() const { return (int)_part.size(); }
 289   double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
 290   int getNumberOfCells() const { return (int)_crit.size(); }
 291   void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
 292   std::vector<bool>& getCriterion() { return _crit; }
 293   const std::vector<bool>& getConstCriterion() const { return _crit; }
 294   void setPart(const std::vector< std::pair<int,int> >& part) { _part=part; }
 295   std::vector< std::pair<int,int> >& getPart() { return _part; }
 296   const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
 297   bool presenceOfTrue() const { return _nb_of_true>0; }
 298   std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
 299   std::vector< std::vector<int> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
 300   double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
 301   void zipToFitOnCriterion();
 302   void updateNumberOfTrue() const;
 303   MEDCouplingAutoRefCountObjectPtr<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
 304   MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
 305 protected:
 306   ~InternalPatch() { }
 307 private:
 308   mutable int _nb_of_true;
 309   std::vector<bool> _crit;
 310   //! _part is global
 311   std::vector< std::pair<int,int> > _part;
 312 };
 313
 314 void InternalPatch::zipToFitOnCriterion()
 315 {
 316   std::vector<int> cgs(computeCGS());
 317   std::vector<bool> newCrit;
 318   std::vector< std::pair<int,int> > newPart,newPart2;
 319   int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart));
 320   MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
 321   if(newNbOfTrue!=_nb_of_true)
 322     throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
 323   _crit=newCrit; _part=newPart2;
 324 }
 325
 326 void InternalPatch::updateNumberOfTrue() const
 327 {
 328   _nb_of_true=(int)std::count(_crit.begin(),_crit.end(),true);
 329 }
 330
 331 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const
 332 {
 333   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
 334   std::vector<int> cgs(computeCGS());
 335   std::vector< std::pair<int,int> > newPart;
 336   MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
 337   MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
 338   ret->setPart(partInGlobal);
 339   ret->updateNumberOfTrue();
 340   return ret;
 341 }
 342
 343 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::deepCpy() const
 344 {
 345   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
 346   (*ret)=*this;
 347   return ret;
 348 }
 349
 350 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
 351 {
 352   cutFound=false; cutPlace=-1;
 353   std::vector<double> ratio(rangeOfAxisId-1);
 354   for(int id=0;id<rangeOfAxisId-1;id++)
 355     {
 356       double efficiency[2];
 357       for(int h=0;h<2;h++)
 358         {
 359           std::vector< std::pair<int,int> > rectH(patchToBeSplit->getConstPart());
 360           if(h==0)
 361             rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id;
 362           else
 363             rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id;
 364           MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(patchToBeSplit->deepCpy());
 365           p->zipToFitOnCriterion();
 366           //anouar rectH ?
 367           efficiency[h]=p->getEfficiencyPerAxis(axisId);
 368         }
 369       ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]);
 370     }
 371   int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
 372   int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
 373   std::vector<double> ratio_inner(dimRatioInner);
 374   double minRatio(1.e10);
 375   for(int i=0; i<dimRatioInner; i++)
 376     {
 377       if(ratio[minCellDirection-1+i]<minRatio)
 378         {
 379           minRatio=ratio[minCellDirection-1+i];
 380           indexMin=i+minCellDirection;
 381         }
 382     }
 383   cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
 384 }
 385
 386 void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int& axisId, bool& cutFound, int& cutPlace)
 387 {
 388   cutPlace=-1; cutFound=false;
 389   int minCellDirection(bso.getMinCellDirection());
 390   const int dim(patchToBeSplit->getDimension());
 391   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
 392   for(int id=0;id<dim;id++)
 393     {
 394       const std::vector<int>& signature(signatures[id]);
 395       std::vector<int> hole;
 396       std::vector<double> distance;
 397       int len((int)signature.size());
 398       for(int i=0;i<len;i++)
 399         if(signature[i]==0)
 400           if(len>= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1)
 401             hole.push_back(i);
 402       if(!hole.empty())
 403         {
 404           double center(((double)len/2.));
 405           for(std::size_t i=0;i<hole.size();i++)
 406             distance.push_back(fabs(hole[i]+1.-center));
 407
 408           std::size_t posDistanceMin(std::distance(distance.begin(),std::min_element(distance.begin(),distance.end())));
 409           cutFound=true;
 410           axisId=id;
 411           cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
 412           return ;
 413         }
 414     }
 415 }
 416
 417 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
 418 {
 419   cutFound=false; cutPlace=-1;// do not set axisId before to be sure that cutFound was set to true
 420
 421   const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
 422   int sign,minCellDirection(bso.getMinCellDirection());
 423   const int dim(patchToBeSplit->getDimension());
 424
 425   std::vector<int> zeroCrossDims(dim,-1);
 426   std::vector<int> zeroCrossVals(dim,-1);
 427   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
 428   for (int id=0;id<dim;id++)
 429     {
 430       const std::vector<int>& signature(signatures[id]);
 431
 432       std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
 433       std::vector<double> distance ;
 434
 435       for (unsigned int i=1;i<signature.size()-1;i++)
 436         derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
 437
 438       // Gradient absolute value
 439       for ( unsigned int i=1;i<derivate_second_order.size();i++)
 440         gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
 441       if(derivate_second_order.empty())
 442         continue;
 443       for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
 444         {
 445           if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
 446             sign = -1 ;
 447           if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
 448             sign = 1 ;
 449           if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
 450             sign = 0 ;
 451           if ( sign==0 || sign==-1 )
 452             if ( i >= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 )
 453               {
 454                 zero_cross.push_back(i) ;
 455                 edge.push_back(gradient_absolute[i]) ;
 456               }
 457           signe_change.push_back(sign) ;
 458         }
 459       if ( zero_cross.size() > 0 )
 460         {
 461           int max_cross=*max_element(edge.begin(),edge.end()) ;
 462           for (unsigned int i=0;i<edge.size();i++)
 463             if (edge[i]==max_cross)
 464               max_cross_list.push_back(zero_cross[i]+1) ;
 465
 466           double center((signature.size()/2.0));
 467           for (unsigned int i=0;i<max_cross_list.size();i++)
 468             distance.push_back(fabs(max_cross_list[i]+1-center));
 469
 470           float distance_min=*min_element(distance.begin(),distance.end()) ;
 471           int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
 472           int best_place = max_cross_list[pos_distance_min] + part[id].first ;
 473           if ( max_cross >=0 )
 474             {
 475               zeroCrossDims[id] = best_place ;
 476               zeroCrossVals[id] = max_cross ;
 477             }
 478         }
 479       derivate_second_order.clear() ;
 480       gradient_absolute.clear() ;
 481       signe_change.clear() ;
 482       zero_cross.clear() ;
 483       edge.clear() ;
 484       max_cross_list.clear() ;
 485       distance.clear() ;
 486     }
 487
 488   if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1  )
 489     {
 490       int max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
 491
 492       if (zeroCrossVals[0]==max_cross_dims &&  zeroCrossVals[1]==max_cross_dims )
 493         {
 494           int nl_left(part[0].second-part[0].first);
 495           int nc_left(part[1].second-part[1].first);
 496           if ( nl_left >=  nc_left )
 497             max_cross_dims = 0 ;
 498           else
 499             max_cross_dims = 1 ;
 500         }
 501       else
 502         max_cross_dims=std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin();
 503       cutFound=true;
 504       cutPlace=zeroCrossDims[max_cross_dims];
 505       axisId=max_cross_dims ;
 506     }
 507 }
 508
 509 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
 510 {
 511   cutFound=false;
 512   if(patchToBeSplit->getEfficiency()<=bso.getEffeciencySnd())
 513     {
 514       if(rangeOfAxisId>=2*bso.getMinCellDirection())
 515         {
 516           cutFound=true;
 517           cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
 518         }
 519     }
 520   else
 521     {
 522       if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells())
 523         {
 524           DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace);
 525         }
 526     }
 527 }
 528
 529 MEDCouplingAutoRefCountObjectPtr<InternalPatch> DealWithNoCut(const InternalPatch *patch)
 530 {
 531   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(const_cast<InternalPatch *>(patch));
 532   ret->incrRef();
 533   return ret;
 534 }
 535
 536 void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
 537 {
 538   MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
 539   std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
 540   std::vector< std::pair<int,int> > leftRect(rect),rightRect(rect);
 541   leftRect[axisId].second=cutPlace+1;
 542   rightRect[axisId].first=cutPlace+1;
 543   leftPart=patchToBeSplit->extractPart(leftRect);
 544   rightPart=patchToBeSplit->extractPart(rightRect);
 545   leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion();
 546   listOfPatches.push_back(leftPart);
 547   listOfPatches.push_back(rightPart);
 548 }
 549
 550 /// @endcond
 551
 552 /*!
 553  * 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.
 554  */
 555 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const std::vector<bool>& criterion, const std::vector<int>& factors)
 556 {
 557   int nbCells(getNumberOfCellsAtCurrentLevel());
 558   if(nbCells!=(int)criterion.size())
 559     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 !");
 560   _patches.clear();
 561   std::vector<int> cgs(_mesh->getCellGridStructure());
 562   std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
 563   //
 564   MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
 565   p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,criterion,p->getCriterion(),p->getPart()));
 566   if(p->presenceOfTrue())
 567     listOfPatches.push_back(p);
 568   while(!listOfPatches.empty())
 569     {
 570       std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
 571       for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
 572         {
 573           //
 574           int axisId,rangeOfAxisId,cutPlace;
 575           bool cutFound;
 576           MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
 577           if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
 578             { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
 579           FindHole(bso,*it,axisId,cutFound,cutPlace);//axisId overwritten here if FindHole equal to true !
 580           if(cutFound)
 581             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
 582           FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
 583           if(cutFound)
 584             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
 585           TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
 586           if(cutFound)
 587             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
 588           listOfPatchesOK.push_back(DealWithNoCut(*it));
 589         }
 590       listOfPatches=listOfPatchesTmp;
 591     }
 592   for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
 593     addPatch((*it)->getConstPart(),factors);
 594 }
 595
 596 /*!
 597  * 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.
 598  */
 599 void MEDCouplingCartesianAMRMeshGen::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
 600 {
 601   if(!criterion || !criterion->isAllocated())
 602     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
 603   std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
 604   createPatchesFromCriterion(bso,crit,factors);
 605 }
 606
 607 void MEDCouplingCartesianAMRMeshGen::removeAllPatches()
 608 {
 609   _patches.clear();
 610   declareAsNew();
 611 }
 612
 613 void MEDCouplingCartesianAMRMeshGen::removePatch(int patchId)
 614 {
 615   checkPatchId(patchId);
 616   int sz((int)_patches.size()),j(0);
 617   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
 618   for(int i=0;i<sz;i++)
 619     if(i!=patchId)
 620       patches[j++]=_patches[i];
 621   (const_cast<MEDCouplingCartesianAMRMeshGen *>(_patches[patchId]->getMesh()))->detachFromFather();
 622   _patches=patches;
 623   declareAsNew();
 624 }
 625
 626 int MEDCouplingCartesianAMRMeshGen::getNumberOfPatches() const
 627 {
 628   return (int)_patches.size();
 629 }
 630
 631 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMeshGen::getPatch(int patchId) const
 632 {
 633   checkPatchId(patchId);
 634   return _patches[patchId];
 635 }
 636
 637 /*!
 638  * This method states if patch2 (with id \a patchId2) is in the neighborhood of patch1 (with id \a patchId1).
 639  * The neighborhood size is defined by \a ghostLev in the reference of \a this ( \b not in the reference of patches !).
 640  */
 641 bool MEDCouplingCartesianAMRMeshGen::isPatchInNeighborhoodOf(int patchId1, int patchId2, int ghostLev) const
 642 {
 643   if(ghostLev<0)
 644     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::isPatchInNeighborhoodOf : the ghost size must be >=0 !");
 645   const MEDCouplingCartesianAMRPatch *p1(getPatch(patchId1)),*p2(getPatch(patchId2));
 646   if(_factors.empty())
 647     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::isPatchInNeighborhoodOf : no factors defined !");
 648   int ghostLevInPatchRef;
 649   if(ghostLev==0)
 650     ghostLevInPatchRef=0;
 651   else
 652     {
 653       ghostLevInPatchRef=(ghostLev-1)/_factors[0]+1;
 654       for(std::size_t i=0;i<_factors.size();i++)
 655         ghostLevInPatchRef=std::max(ghostLevInPatchRef,(ghostLev-1)/_factors[i]+1);
 656     }
 657   return p1->isInMyNeighborhood(p2,ghostLevInPatchRef);
 658 }
 659
 660 /*!
 661  * 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.
 662  * 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
 663  * defined by the patch with id \a patchId.
 664  *
 665  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 666  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 667  * \return DataArrayDouble * - The array of the cell field on the requested patch
 668  *
 669  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
 670  * \throw if \a cellFieldOnThis is NULL or not allocated
 671  * \sa fillCellFieldOnPatch, MEDCouplingIMesh::SpreadCoarseToFine
 672  */
 673 DataArrayDouble *MEDCouplingCartesianAMRMeshGen::createCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis) const
 674 {
 675   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
 676     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
 677   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 678   const MEDCouplingIMesh *fine(patch->getMesh()->getImageMesh());
 679   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(fine->getNumberOfCells(),cellFieldOnThis->getNumberOfComponents());
 680   ret->copyStringInfoFrom(*cellFieldOnThis);
 681   MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),ret,patch->getBLTRRange(),getFactors());
 682   return ret.retn();
 683 }
 684
 685 /*!
 686  * This method is equivalent to MEDCouplingCartesianAMRMesh::createCellFieldOnPatch except that here instead of creating a new instance
 687  * it fills the value into the \a cellFieldOnPatch data.
 688  *
 689  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 690  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 691  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
 692  *
 693  * \sa createCellFieldOnPatch, fillCellFieldComingFromPatch
 694  */
 695 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatch(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch) const
 696 {
 697   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
 698     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatch : the input cell field array is NULL or not allocated !");
 699   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 700   MEDCouplingIMesh::SpreadCoarseToFine(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors());
 701 }
 702
 703 /*!
 704  * This method is the generalization of fillCellFieldOnPatch method. This method only projects coarse to fine without considering the
 705  * potential neighbor patches covered by the ghost cells of patch with id \a patchId.
 706  *
 707  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 708  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 709  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
 710  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
 711  *
 712  * \sa fillCellFieldOnPatch, fillCellFieldOnPatchGhostAdv
 713  */
 714 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhost(int patchId, const DataArrayDouble *cellFieldOnThis, DataArrayDouble *cellFieldOnPatch, int ghostLev) const
 715 {
 716   if(!cellFieldOnThis || !cellFieldOnThis->isAllocated())
 717     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createCellFieldOnPatchGhost : the input cell field array is NULL or not allocated !");
 718   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 719   MEDCouplingIMesh::SpreadCoarseToFineGhost(cellFieldOnThis,_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),ghostLev);
 720 }
 721
 722 /*!
 723  * This method is a refinement of fillCellFieldOnPatchGhost. fillCellFieldOnPatchGhost is first called.
 724  * Then for all other patches than those pointed by \a patchId that overlap the ghost zone of the patch impact the ghost zone adequately.
 725  *
 726  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 727  * \param [in] cellFieldOnThis - The array of the cell field on \c this->getImageMesh() to be projected to patch having id \a patchId.
 728  * \param [in,out] cellFieldOnPatch - The array of the cell field on the requested patch to be filled.
 729  * \param [in] ghostLev - The size of the ghost zone (must be >=0 !)
 730  * \param [in] arrsOnPatches - \b WARNING arrsOnPatches[patchId] is \b NOT \b const. All others are const.
 731  */
 732 void MEDCouplingCartesianAMRMeshGen::fillCellFieldOnPatchGhostAdv(int patchId, const DataArrayDouble *cellFieldOnThis, int ghostLev, const std::vector<const DataArrayDouble *>& arrsOnPatches) const
 733 {
 734   int nbp(getNumberOfPatches()),dim(getSpaceDimension());
 735   if(nbp!=(int)arrsOnPatches.size())
 736     {
 737       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::fillCellFieldOnPatchGhostAdv : there are " << nbp << " patches in this and " << arrsOnPatches.size() << " arrays in the last parameter !";
 738       throw INTERP_KERNEL::Exception(oss.str().c_str());
 739     }
 740   DataArrayDouble *theFieldToFill(const_cast<DataArrayDouble *>(arrsOnPatches[patchId]));
 741   // first, do as usual
 742   fillCellFieldOnPatchGhost(patchId,cellFieldOnThis,theFieldToFill,ghostLev);
 743   // all reference patch stuff
 744   const MEDCouplingCartesianAMRPatch *refP(getPatch(patchId));
 745   const std::vector< std::pair<int,int> >& refBLTR(refP->getBLTRRange());//[(1,4),(2,4)]
 746   std::vector<int> dimsCoarse(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(refBLTR));//[3,2]
 747   std::transform(dimsCoarse.begin(),dimsCoarse.end(),_factors.begin(),dimsCoarse.begin(),std::multiplies<int>());//[12,8]
 748   std::transform(dimsCoarse.begin(),dimsCoarse.end(),dimsCoarse.begin(),std::bind2nd(std::plus<int>(),2*ghostLev));//[14,10]
 749   std::vector< std::pair<int,int> > rangeCoarse(MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(dimsCoarse));//[(0,14),(0,10)]
 750   std::vector<int> fakeFactors(dim,1);
 751   //
 752   for(int i=0;i<nbp;i++)
 753     {
 754       if(i!=patchId)
 755         if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
 756           {
 757             const MEDCouplingCartesianAMRPatch *otherP(getPatch(i));
 758             const std::vector< std::pair<int,int> >& otherBLTR(otherP->getBLTRRange());//[(4,5),(3,4)]
 759             std::vector< std::pair<int,int> > tmp0,tmp1,tmp2;
 760             MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(refBLTR,otherBLTR,tmp0,false);//tmp0=[(3,4),(1,2)]
 761             ApplyFactorsOnCompactFrmt(tmp0,_factors);//tmp0=[(12,16),(4,8)]
 762             ApplyGhostOnCompactFrmt(tmp0,ghostLev);//tmp0=[(13,17),(5,9)]
 763             std::vector< std::pair<int,int> > interstRange(MEDCouplingStructuredMesh::IntersectRanges(tmp0,rangeCoarse));//interstRange=[(13,14),(5,9)]
 764             MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(otherBLTR,refBLTR,tmp1,false);//tmp1=[(-3,0),(-1,1)]
 765             ApplyFactorsOnCompactFrmt(tmp1,_factors);//tmp1=[(-12,-4),(-4,0)]
 766             MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(tmp1,interstRange,tmp2,false);//tmp2=[(1,2),(1,5)]
 767             //
 768             std::vector< std::pair<int,int> > dimsFine(otherBLTR);
 769             ApplyFactorsOnCompactFrmt(dimsFine,_factors);
 770             ApplyAllGhostOnCompactFrmt(dimsFine,ghostLev);
 771             //
 772             MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ghostVals(MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(dimsFine),arrsOnPatches[i],tmp2));
 773             MEDCouplingIMesh::CondenseFineToCoarse(dimsCoarse,ghostVals,interstRange,fakeFactors,theFieldToFill);
 774           }
 775     }
 776 }
 777
 778 /*!
 779  * This method updates \a cellFieldOnThis part of values coming from the cell field \a cellFieldOnPatch lying on patch having id \a patchId.
 780  *
 781  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 782  * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
 783  * \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.
 784  *
 785  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
 786  * \throw if \a cellFieldOnPatch is NULL or not allocated
 787  * \sa createCellFieldOnPatch, MEDCouplingIMesh::CondenseFineToCoarse,fillCellFieldComingFromPatchGhost
 788  */
 789 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatch(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis) const
 790 {
 791   if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
 792       throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch : the input cell field array is NULL or not allocated !");
 793   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 794   MEDCouplingIMesh::CondenseFineToCoarse(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis);
 795 }
 796
 797 /*!
 798  * This method is the extension of MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch managing the ghost cells. If this
 799  * method is called with \a ghostLev equal to 0 it behaves exactly as MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatch.
 800  *
 801  * \param [in] patchId - The id of the patch \a cellFieldOnThis has to be put on.
 802  * \param [in] cellFieldOnPatch - The array of the cell field on patch with id \a patchId.
 803  * \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.
 804  * \param [in] ghostLev The size of ghost zone (must be >= 0 !)
 805  *
 806  * \throw if \a patchId is not in [ 0 , \c this->getNumberOfPatches() )
 807  * \throw if \a cellFieldOnPatch is NULL or not allocated
 808  * \sa fillCellFieldComingFromPatch
 809  */
 810 void MEDCouplingCartesianAMRMeshGen::fillCellFieldComingFromPatchGhost(int patchId, const DataArrayDouble *cellFieldOnPatch, DataArrayDouble *cellFieldOnThis, int ghostLev) const
 811 {
 812   if(!cellFieldOnPatch || !cellFieldOnPatch->isAllocated())
 813     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::fillCellFieldComingFromPatchGhost : the input cell field array is NULL or not allocated !");
 814   const MEDCouplingCartesianAMRPatch *patch(getPatch(patchId));
 815   MEDCouplingIMesh::CondenseFineToCoarseGhost(_mesh->getCellGridStructure(),cellFieldOnPatch,patch->getBLTRRange(),getFactors(),cellFieldOnThis,ghostLev);
 816 }
 817
 818 /*!
 819  * This method finds all patches (located by their ids) that are in the neighborhood of patch with id \a patchId. The neighborhood size is
 820  * defined by ghostLev.
 821  *
 822  * \param [in] patchId - the id of the considered patch.
 823  * \param [in] ghostLev - the size of the neighborhood.
 824  * \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.
 825  */
 826 DataArrayInt *MEDCouplingCartesianAMRMeshGen::findPatchesInTheNeighborhoodOf(int patchId, int ghostLev) const
 827 {
 828   int nbp(getNumberOfPatches());
 829   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
 830   for(int i=0;i<nbp;i++)
 831     {
 832       if(i!=patchId)
 833         if(isPatchInNeighborhoodOf(i,patchId,ghostLev))
 834           ret->pushBackSilent(i);
 835     }
 836   return ret.retn();
 837 }
 838
 839 MEDCouplingUMesh *MEDCouplingCartesianAMRMeshGen::buildUnstructured() const
 840 {
 841   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
 842   std::vector<bool> bs(_mesh->getNumberOfCells(),false);
 843   std::vector<int> cgs(_mesh->getCellGridStructure());
 844   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
 845   std::size_t ii(0);
 846   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
 847     {
 848       MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
 849       msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
 850     }
 851   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
 852   msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
 853   std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
 854   for(std::size_t i=0;i<msSafe.size();i++)
 855     ms[i]=msSafe[i];
 856   return MEDCouplingUMesh::MergeUMeshes(ms);
 857 }
 858
 859 /*!
 860  * This method returns a mesh containing as cells that there is patches at the current level.
 861  * The patches are seen like 'boxes' that is too say the refinement will not appear here.
 862  *
 863  * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
 864  */
 865 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshFromPatchEnvelop() const
 866 {
 867   std::vector<const MEDCoupling1SGTUMesh *> cells;
 868   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > cellsSafe;
 869   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 870     {
 871       const MEDCouplingCartesianAMRPatch *patch(*it);
 872       if(patch)
 873         {
 874           MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
 875           MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
 876           cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
 877         }
 878     }
 879   return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
 880 }
 881
 882 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMeshGen::buildMeshOfDirectChildrenOnly() const
 883 {
 884   std::vector<const MEDCoupling1SGTUMesh *> patches;
 885   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > patchesSafe;
 886   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 887       {
 888         const MEDCouplingCartesianAMRPatch *patch(*it);
 889         if(patch)
 890           {
 891             MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> patchMesh(patch->getMesh()->getImageMesh()->build1SGTUnstructured());
 892             patchesSafe.push_back(patchMesh); patches.push_back(patchMesh);
 893           }
 894       }
 895     return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(patches);
 896 }
 897
 898 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
 899                                                                const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
 900 {
 901   _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
 902 }
 903
 904 MEDCouplingCartesianAMRMeshGen::MEDCouplingCartesianAMRMeshGen(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):_father(father)
 905 {
 906   if(!_father)
 907     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : empty father !");
 908   if(!mesh)
 909     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
 910   mesh->checkCoherency();
 911   _mesh=mesh; _mesh->incrRef();
 912 }
 913
 914 void MEDCouplingCartesianAMRMeshGen::checkPatchId(int patchId) const
 915 {
 916   int sz(getNumberOfPatches());
 917   if(patchId<0 || patchId>=sz)
 918     {
 919       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMeshGen::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
 920       throw INTERP_KERNEL::Exception(oss.str().c_str());
 921     }
 922 }
 923
 924 void MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign(const std::vector<int>& factors)
 925 {
 926   if(getSpaceDimension()!=(int)factors.size())
 927     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::checkFactorsAndIfNotSetAssign : invalid size of factors ! size must be equal to the spaceDimension !");
 928   if(_factors.empty())
 929     {
 930       _factors=factors;
 931     }
 932   else
 933     {
 934       if(_factors!=factors)
 935         throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::checkFactorsAndIfNotSetAssign : the factors ");
 936     }
 937 }
 938
 939 void MEDCouplingCartesianAMRMeshGen::retrieveGridsAtInternal(int lev, std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGen> >& grids) const
 940 {
 941   if(lev==0)
 942     {
 943       const MEDCouplingCartesianAMRMesh *thisc(dynamic_cast<const MEDCouplingCartesianAMRMesh *>(this));//tony
 944       MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatchGF> elt(new MEDCouplingCartesianAMRPatchGF(const_cast<MEDCouplingCartesianAMRMesh *>(thisc)));
 945       grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatchGF,MEDCouplingCartesianAMRPatchGen>(elt));
 946     }
 947   else if(lev==1)
 948     {
 949       for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 950         {
 951           const MEDCouplingCartesianAMRPatch *pt(*it);
 952           if(pt)
 953             {
 954               MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> tmp1(*it);
 955               grids.push_back(DynamicCastSafe<MEDCouplingCartesianAMRPatch,MEDCouplingCartesianAMRPatchGen>(tmp1));
 956             }
 957         }
 958     }
 959   else
 960     {
 961       for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 962         {
 963           const MEDCouplingCartesianAMRPatch *pt(*it);
 964           if(pt)
 965             pt->getMesh()->retrieveGridsAtInternal(lev-1,grids);
 966         }
 967     }
 968 }
 969
 970 /*!
 971  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in refined reference.
 972  * \param [in] factors - the factors per axis.
 973  */
 974 void MEDCouplingCartesianAMRMeshGen::ApplyFactorsOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, const std::vector<int>& factors)
 975 {
 976   std::size_t sz(factors.size());
 977   if(sz!=partBeforeFact.size())
 978     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::ApplyFactorsOnCompactFrmt : size of input vectors must be the same !");
 979   for(std::size_t i=0;i<sz;i++)
 980     {
 981       partBeforeFact[i].first*=factors[i];
 982       partBeforeFact[i].second*=factors[i];
 983     }
 984 }
 985
 986 /*!
 987  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
 988  * \param [in] ghostSize - the ghost size of zone for all axis.
 989  */
 990 void MEDCouplingCartesianAMRMeshGen::ApplyGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
 991 {
 992   if(ghostSize<0)
 993     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::ApplyGhostOnCompactFrmt : ghost size must be >= 0 !");
 994   std::size_t sz(partBeforeFact.size());
 995   for(std::size_t i=0;i<sz;i++)
 996     {
 997       partBeforeFact[i].first+=ghostSize;
 998       partBeforeFact[i].second+=ghostSize;
 999     }
1000 }
1001
1002 /*!
1003  * This method is different than ApplyGhostOnCompactFrmt
1004  *
1005  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
1006  * \param [in] ghostSize - the ghost size of zone for all axis.
1007  */
1008 void MEDCouplingCartesianAMRMeshGen::ApplyAllGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
1009 {
1010   if(ghostSize<0)
1011     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMeshGen::ApplyAllGhostOnCompactFrmt : ghost size must be >= 0 !");
1012   std::size_t sz(partBeforeFact.size());
1013   for(std::size_t i=0;i<sz;i++)
1014     {
1015       partBeforeFact[i].first-=ghostSize;
1016       partBeforeFact[i].second+=ghostSize;
1017     }
1018 }
1019
1020 std::size_t MEDCouplingCartesianAMRMeshGen::getHeapMemorySizeWithoutChildren() const
1021 {
1022   return sizeof(MEDCouplingCartesianAMRMeshGen);
1023 }
1024
1025 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMeshGen::getDirectChildren() const
1026 {
1027   std::vector<const BigMemoryObject *> ret;
1028   if((const MEDCouplingIMesh *)_mesh)
1029     ret.push_back((const MEDCouplingIMesh *)_mesh);
1030   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1031     {
1032       if((const MEDCouplingCartesianAMRPatch*)*it)
1033         ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
1034     }
1035   return ret;
1036 }
1037
1038 void MEDCouplingCartesianAMRMeshGen::updateTime() const
1039 {
1040   if((const MEDCouplingIMesh *)_mesh)
1041     updateTimeWith(*_mesh);
1042   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
1043     {
1044       const MEDCouplingCartesianAMRPatch *elt(*it);
1045       if(!elt)
1046         continue;
1047       const MEDCouplingCartesianAMRMeshGen *mesh(elt->getMesh());
1048       if(mesh)
1049         updateTimeWith(*mesh);
1050     }
1051 }
1052
1053 MEDCouplingCartesianAMRMeshSub::MEDCouplingCartesianAMRMeshSub(MEDCouplingCartesianAMRMeshGen *father, MEDCouplingIMesh *mesh):MEDCouplingCartesianAMRMeshGen(father,mesh)
1054 {
1055 }
1056
1057 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1058                                                               const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
1059 {
1060   return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
1061 }
1062
1063 void MEDCouplingCartesianAMRMesh::setData(MEDCouplingDataForGodFather *data)
1064 {
1065   _data=data;
1066 }
1067
1068 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
1069                                                          const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):MEDCouplingCartesianAMRMeshGen(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop)
1070 {
1071 }
1072
1073 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
1074 {
1075   std::vector<const BigMemoryObject *> ret(MEDCouplingCartesianAMRMeshGen::getDirectChildren());
1076   const MEDCouplingDataForGodFather *pt(_data);
1077   if(pt)
1078     ret.push_back(pt);
1079   return ret;
1080 }