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 /*!
  35  * \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.
  36  * \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,
  37  *                                a the end cell (\b excluded) of the range for the second element of the pair.
  38  * \param [in] factors The refinement per axis relative to the father of \a this.
  39  */
  40 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMesh *mesh, const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
  41 {
  42   if(!mesh)
  43     throw INTERP_KERNEL::Exception("EDCouplingCartesianAMRPatch constructor : input mesh is NULL !");
  44   _mesh=mesh; _mesh->incrRef();
  45   int dim((int)bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
  46   if(dim!=dimExp)
  47     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
  48   _bl_tr=bottomLeftTopRight;
  49   if((int)factors.size()!=dimExp)
  50     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input factors per axis size mismatches !");
  51   _factors=factors;
  52 }
  53
  54 int MEDCouplingCartesianAMRPatch::getNumberOfCellsRecursiveWithOverlap() const
  55 {
  56   return _mesh->getNumberOfCellsRecursiveWithOverlap();
  57 }
  58
  59 int MEDCouplingCartesianAMRPatch::getNumberOfCellsRecursiveWithoutOverlap() const
  60 {
  61   return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
  62 }
  63
  64 int MEDCouplingCartesianAMRPatch::getMaxNumberOfLevelsRelativeToThis() const
  65 {
  66   return _mesh->getMaxNumberOfLevelsRelativeToThis();
  67 }
  68
  69 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
  70 {
  71   return _mesh->addPatch(bottomLeftTopRight,factors);
  72 }
  73
  74 int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
  75 {
  76   return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
  77 }
  78
  79 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
  80 {
  81   std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
  82   ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
  83   return ret;
  84 }
  85
  86 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatch::getDirectChildren() const
  87 {
  88   std::vector<const BigMemoryObject *> ret;
  89   if((const MEDCouplingCartesianAMRMesh *)_mesh)
  90     ret.push_back((const MEDCouplingCartesianAMRMesh *)_mesh);
  91   return ret;
  92 }
  93
  94 /// @endcond
  95
  96
  97 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
  98                                                               const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
  99 {
 100   return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
 101 }
 102
 103 int MEDCouplingCartesianAMRMesh::getSpaceDimension() const
 104 {
 105   return _mesh->getSpaceDimension();
 106 }
 107
 108 int MEDCouplingCartesianAMRMesh::getMaxNumberOfLevelsRelativeToThis() const
 109 {
 110   int ret(1);
 111   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 112     ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
 113   return ret;
 114 }
 115
 116 int MEDCouplingCartesianAMRMesh::getNumberOfCellsAtCurrentLevel() const
 117 {
 118   return _mesh->getNumberOfCells();
 119 }
 120
 121 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithOverlap() const
 122 {
 123   int ret(_mesh->getNumberOfCells());
 124   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 125     {
 126       ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
 127     }
 128   return ret;
 129 }
 130
 131 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithoutOverlap() const
 132 {
 133   int ret(_mesh->getNumberOfCells());
 134   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 135     {
 136       ret-=(*it)->getNumberOfOverlapedCellsForFather();
 137       ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
 138     }
 139   return ret;
 140 }
 141
 142 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getFather() const
 143 {
 144   return _father;
 145 }
 146
 147 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getGodFather() const
 148 {
 149   if(_father==0)
 150     return this;
 151   else
 152     return _father->getGodFather();
 153 }
 154
 155 void MEDCouplingCartesianAMRMesh::detachFromFather()
 156 {
 157   _father=0;
 158 }
 159
 160 /*!
 161  * \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,
 162  *                                a the end cell (\b excluded) of the range for the second element of the pair.
 163  * \param [in] factors The != 0 factor of refinement per axis.
 164  */
 165 void MEDCouplingCartesianAMRMesh::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
 166 {
 167   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
 168   mesh->refineWithFactor(factors);
 169   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMesh> zeMesh(new MEDCouplingCartesianAMRMesh(this,mesh));
 170   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight,factors));
 171   _patches.push_back(elt);
 172 }
 173
 174 /// @cond INTERNAL
 175
 176 class InternalPatch : public RefCountObjectOnly
 177 {
 178 public:
 179   InternalPatch():_nb_of_true(0) { }
 180   double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
 181   int getNumberOfCells() const { return (int)_crit.size(); }
 182   void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
 183   std::vector<bool>& getCriterion() { return _crit; }
 184   const std::vector<bool>& getConstCriterion() const { return _crit; }
 185   void setPart(const std::vector< std::pair<int,int> >& part) { _part=part; }
 186   std::vector< std::pair<int,int> >& getPart() { return _part; }
 187   const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
 188   bool presenceOfTrue() const { return _nb_of_true>0; }
 189   std::vector<int> computeCGS() const { return MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(_part); }
 190   std::vector< std::vector<int> > computeSignature() const { return MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(computeCGS(),getConstCriterion()); }
 191   double getEfficiencyPerAxis(int axisId) const { return (double)_nb_of_true/((double)(_part[axisId].second-_part[axisId].first)); }
 192   void zipToFitOnCriterion();
 193   void updateNumberOfTrue() const;
 194   MEDCouplingAutoRefCountObjectPtr<InternalPatch> extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const;
 195   MEDCouplingAutoRefCountObjectPtr<InternalPatch> deepCpy() const;
 196 protected:
 197   ~InternalPatch() { }
 198 private:
 199   mutable int _nb_of_true;
 200   std::vector<bool> _crit;
 201   //! _part is global
 202   std::vector< std::pair<int,int> > _part;
 203 };
 204
 205 void InternalPatch::zipToFitOnCriterion()
 206 {
 207   std::vector<int> cgs(computeCGS());
 208   std::vector<bool> newCrit;
 209   std::vector< std::pair<int,int> > newPart,newPart2;
 210   int newNbOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,_crit,newCrit,newPart));
 211   MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(_part,newPart,newPart2);
 212   if(newNbOfTrue!=_nb_of_true)
 213     throw INTERP_KERNEL::Exception("InternalPatch::zipToFitOnCrit : internal error !");
 214   _crit=newCrit; _part=newPart2;
 215 }
 216
 217 void InternalPatch::updateNumberOfTrue() const
 218 {
 219   _nb_of_true=(int)std::count(_crit.begin(),_crit.end(),true);
 220 }
 221
 222 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::extractPart(const std::vector< std::pair<int,int> >&partInGlobal) const
 223 {
 224   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
 225   std::vector<int> cgs(computeCGS());
 226   std::vector< std::pair<int,int> > newPart;
 227   MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(_part,partInGlobal,newPart);
 228   MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(cgs,_crit,newPart,ret->getCriterion());
 229   ret->setPart(partInGlobal);
 230   ret->updateNumberOfTrue();
 231   return ret;
 232 }
 233
 234 MEDCouplingAutoRefCountObjectPtr<InternalPatch> InternalPatch::deepCpy() const
 235 {
 236   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(new InternalPatch);
 237   (*ret)=*this;
 238   return ret;
 239 }
 240
 241 void DissectBigPatch(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
 242 {
 243   cutFound=false; cutPlace=-1;
 244   std::vector<double> ratio(rangeOfAxisId-1);
 245   for(int id=0;id<rangeOfAxisId-1;id++)
 246     {
 247       double efficiency[2];
 248       for(int h=0;h<2;h++)
 249         {
 250           std::vector< std::pair<int,int> > rectH(patchToBeSplit->getConstPart());
 251           if(h==0)
 252             rectH[axisId].second=patchToBeSplit->getConstPart()[axisId].first+id;
 253           else
 254             rectH[axisId].first=patchToBeSplit->getConstPart()[axisId].first+id;
 255           MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(patchToBeSplit->deepCpy());
 256           p->zipToFitOnCriterion();
 257           //anouar rectH ?
 258           efficiency[h]=p->getEfficiencyPerAxis(axisId);
 259         }
 260       ratio[id]=std::max(efficiency[0],efficiency[1])/std::min(efficiency[0],efficiency[1]);
 261     }
 262   int minCellDirection(bso.getMinCellDirection()),indexMin(-1);
 263   int dimRatioInner(rangeOfAxisId-1-2*(minCellDirection-1));
 264   std::vector<double > ratio_inner(dimRatioInner);
 265   double minRatio(1.e10);
 266   for(int i=0; i<dimRatioInner; i++)
 267     {
 268       if(ratio[minCellDirection-1+i]<minRatio)
 269         {
 270           minRatio=ratio[minCellDirection-1+i];
 271           indexMin=i+minCellDirection;
 272         }
 273     }
 274   cutFound=true; cutPlace=indexMin+patchToBeSplit->getConstPart()[axisId].first-1;
 275 }
 276
 277 void FindHole(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, const int axisId, bool& cutFound, int& cutPlace)
 278 {
 279   cutPlace=-1; cutFound=false;
 280   int minCellDirection(bso.getMinCellDirection());
 281
 282   int sortedDims[2];
 283   sortedDims[0]=axisId==0?1:0;
 284   sortedDims[1]=axisId==0?0:1;
 285   static const int dim(2);
 286   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
 287   for(int id=0;id<dim;id++)
 288     {
 289       const std::vector<int>& signature(signatures[sortedDims[id]]);
 290       std::vector<int> hole;
 291       std::vector<double> distance ;
 292       int len((int)signature.size());
 293       for(int i=0;i<len;i++)
 294         if(signature[i]==0)
 295           if(len>= 2*minCellDirection && i >= minCellDirection-1 && i <= len-minCellDirection-1)
 296             hole.push_back(i);
 297       if (hole.size()>0)
 298         {
 299           double center(((double)len/2.)+0.5);
 300           for(std::size_t i=0;i<hole.size();i++)
 301             distance.push_back(fabs(hole[i]+1+0.5-center));//anouar ! why 0.5 ?
 302
 303           double distanceMin=*std::min_element(distance.begin(),distance.end());
 304           int posDistanceMin=std::find(distance.begin(),distance.end(),distanceMin)-distance.begin()-1;
 305
 306           cutFound=true;
 307           cutPlace=hole[posDistanceMin]+patchToBeSplit->getConstPart()[axisId].first+1;
 308           return ;
 309         }
 310     }
 311 }
 312
 313 void FindInflection(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, bool& cutFound, int& cutPlace, int& axisId)
 314 {
 315   cutFound=false; cutPlace=-1; axisId=-1;
 316
 317   const std::vector< std::pair<int,int> >& part(patchToBeSplit->getConstPart());
 318   int sign,minCellDirection(bso.getMinCellDirection());
 319   static const int dim = 2 ;
 320
 321   std::vector<int> zeroCrossDims(2,-1);
 322   std::vector<int> zeroCrossVals(2,-1);
 323   std::vector< std::vector<int> > signatures(patchToBeSplit->computeSignature());
 324   for (int id=0;id<dim;id++)
 325     {
 326       const std::vector<int>& signature(signatures[id]);
 327
 328       std::vector<int> derivate_second_order,gradient_absolute,signe_change,zero_cross,edge,max_cross_list ;
 329       std::vector<double> distance ;
 330
 331       for (unsigned int i=1;i<signature.size()-1;i++)
 332         derivate_second_order.push_back(signature[i-1]-2*signature[i]+signature[i+1]) ;
 333
 334       // Gradient absolute value
 335       for ( unsigned int i=1;i<derivate_second_order.size();i++)
 336         gradient_absolute.push_back(fabs(derivate_second_order[i]-derivate_second_order[i-1])) ;
 337       for (unsigned int i=0;i<derivate_second_order.size()-1;i++)
 338         {
 339           if (derivate_second_order[i]*derivate_second_order[i+1] < 0 )
 340             sign = -1 ;
 341           if (derivate_second_order[i]*derivate_second_order[i+1] > 0 )
 342             sign = 1 ;
 343           if (derivate_second_order[i]*derivate_second_order[i+1] == 0 )
 344             sign = 0 ;
 345           if ( sign==0 || sign==-1 )
 346             if ( i >= (unsigned int)minCellDirection-2 && i <= signature.size()-minCellDirection-2 )
 347               {
 348                 zero_cross.push_back(i) ;
 349                 edge.push_back(gradient_absolute[i]) ;
 350               }
 351           signe_change.push_back(sign) ;
 352         }
 353       if ( zero_cross.size() > 0 )
 354         {
 355           int max_cross=*max_element(edge.begin(),edge.end()) ;
 356           for (unsigned int i=0;i<edge.size();i++)
 357             if (edge[i]==max_cross)
 358               max_cross_list.push_back(zero_cross[i]+1) ;
 359
 360           float center = (signature.size()/2.0)+0.5;
 361           for (unsigned int i=0;i<max_cross_list.size();i++)
 362             distance.push_back(fabs(max_cross_list[i]+1+0.5-center));
 363
 364           float distance_min=*min_element(distance.begin(),distance.end()) ;
 365           int pos_distance_min=find(distance.begin(),distance.end(),distance_min)-distance.begin();
 366           int best_place = max_cross_list[pos_distance_min] + part[id].first ;
 367           if ( max_cross >=0 )
 368             {
 369               zeroCrossDims[id] = best_place ;
 370               zeroCrossVals[id] = max_cross ;
 371             }
 372         }
 373       derivate_second_order.clear() ;
 374       gradient_absolute.clear() ;
 375       signe_change.clear() ;
 376       zero_cross.clear() ;
 377       edge.clear() ;
 378       max_cross_list.clear() ;
 379       distance.clear() ;
 380     }
 381
 382   if ( zeroCrossDims[0]!=-1 || zeroCrossDims[1]!=-1  )
 383     {
 384       int max_cross_dims = *max_element(zeroCrossVals.begin(),zeroCrossVals.end()) ;
 385
 386       if (zeroCrossVals[0]==max_cross_dims &&  zeroCrossVals[1]==max_cross_dims )
 387         {
 388           int nl_left(part[0].second-part[0].first);
 389           int nc_left(part[1].second-part[1].first);
 390           if ( nl_left >=  nc_left )
 391             max_cross_dims = 0 ;
 392           else
 393             max_cross_dims = 1 ;
 394         }
 395       else
 396         max_cross_dims=std::find(zeroCrossVals.begin(),zeroCrossVals.end(),max_cross_dims)-zeroCrossVals.begin();
 397       cutFound=true;
 398       cutPlace=zeroCrossDims[max_cross_dims];
 399       axisId=max_cross_dims ;
 400     }
 401 }
 402
 403 void TryAction4(const INTERP_KERNEL::BoxSplittingOptions& bso, const InternalPatch *patchToBeSplit, int axisId, int rangeOfAxisId, bool& cutFound, int& cutPlace)
 404 {
 405   cutFound=false;
 406   if(patchToBeSplit->getEfficiency() <= bso.getEffeciencySnd())
 407     {
 408       if(rangeOfAxisId>=2*bso.getMinCellDirection())
 409         {
 410           cutFound=true;
 411           cutPlace=rangeOfAxisId/2+patchToBeSplit->getConstPart()[axisId].first-1;
 412         }
 413     }
 414   else
 415     {
 416       if(patchToBeSplit->getNumberOfCells()>bso.getMaxCells())
 417         {
 418           DissectBigPatch(bso,patchToBeSplit,axisId,rangeOfAxisId,cutFound,cutPlace);
 419         }
 420     }
 421 }
 422
 423 MEDCouplingAutoRefCountObjectPtr<InternalPatch> DealWithNoCut(const InternalPatch *patch)
 424 {
 425   MEDCouplingAutoRefCountObjectPtr<InternalPatch> ret(const_cast<InternalPatch *>(patch));
 426   ret->incrRef();
 427   return ret;
 428 }
 429
 430 void DealWithCut(const InternalPatch *patchToBeSplit, int axisId, int cutPlace, std::vector<MEDCouplingAutoRefCountObjectPtr<InternalPatch> >& listOfPatches)
 431 {
 432   MEDCouplingAutoRefCountObjectPtr<InternalPatch> leftPart,rightPart;
 433   std::vector< std::pair<int,int> > rect(patchToBeSplit->getConstPart());
 434   std::vector< std::pair<int,int> > leftRect(rect),rightRect(rect);
 435   leftRect[axisId].second=cutPlace+1;
 436   rightRect[axisId].first=cutPlace+1;
 437   leftPart=patchToBeSplit->extractPart(leftRect);
 438   rightPart=patchToBeSplit->extractPart(rightRect);
 439   leftPart->zipToFitOnCriterion(); rightPart->zipToFitOnCriterion();
 440   listOfPatches.push_back(leftPart);
 441   listOfPatches.push_back(rightPart);
 442 }
 443
 444 /// @endcond
 445
 446 /*!
 447  * 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.
 448  */
 449 void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
 450 {
 451   if(!criterion || !criterion->isAllocated())
 452     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
 453   int nbCells(getNumberOfCellsAtCurrentLevel());
 454   if(nbCells!=criterion->getNumberOfTuples())
 455     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 !");
 456   _patches.clear();
 457   //
 458   std::vector<int> cgs(_mesh->getCellGridStructure());
 459   std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
 460   std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
 461   //
 462   MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
 463   p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,crit,p->getCriterion(),p->getPart()));
 464   if(p->presenceOfTrue())
 465     listOfPatches.push_back(p);
 466   while(!listOfPatches.empty())
 467     {
 468       std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
 469       for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
 470         {
 471           //
 472           int axisId,rangeOfAxisId;
 473           bool cutFound;
 474           int cutPlace;
 475           MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt((*it)->getConstPart(),axisId,rangeOfAxisId);
 476           if((*it)->getEfficiency()>=bso.getEffeciency() && (*it)->getNumberOfCells()<bso.getMaxCells())
 477             { listOfPatchesOK.push_back(DealWithNoCut(*it)); continue; }//action 1
 478           FindHole(bso,*it,axisId,cutFound,cutPlace);
 479           if(cutFound)
 480             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 2
 481           FindInflection(bso,*it,cutFound,cutPlace,axisId);//axisId overwritten here if cutFound equal to true !
 482           if(cutFound)
 483             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 3
 484           TryAction4(bso,*it,axisId,rangeOfAxisId,cutFound,cutPlace);
 485           if(cutFound)
 486             { DealWithCut(*it,axisId,cutPlace,listOfPatchesTmp); continue; }//action 4
 487           listOfPatchesOK.push_back(DealWithNoCut(*it));
 488         }
 489       listOfPatches=listOfPatchesTmp;
 490     }
 491   for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
 492     {
 493       addPatch((*it)->getConstPart(),factors);
 494     }
 495 }
 496
 497 void MEDCouplingCartesianAMRMesh::removePatch(int patchId)
 498 {
 499   checkPatchId(patchId);
 500   int sz((int)_patches.size()),j(0);
 501   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
 502   for(int i=0;i<sz;i++)
 503     if(i!=patchId)
 504       patches[j++]=_patches[i];
 505   (const_cast<MEDCouplingCartesianAMRMesh *>(_patches[patchId]->getMesh()))->detachFromFather();
 506   _patches=patches;
 507   declareAsNew();
 508 }
 509
 510 int MEDCouplingCartesianAMRMesh::getNumberOfPatches() const
 511 {
 512   return (int)_patches.size();
 513 }
 514
 515 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMesh::getPatch(int patchId) const
 516 {
 517   checkPatchId(patchId);
 518   return _patches[patchId];
 519 }
 520
 521 MEDCouplingUMesh *MEDCouplingCartesianAMRMesh::buildUnstructured() const
 522 {
 523   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
 524   std::vector<bool> bs(_mesh->getNumberOfCells(),false);
 525   std::vector<int> cgs(_mesh->getCellGridStructure());
 526   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
 527   std::size_t ii(0);
 528   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
 529     {
 530       MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
 531       msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
 532     }
 533   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
 534   msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
 535   std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
 536   for(std::size_t i=0;i<msSafe.size();i++)
 537     ms[i]=msSafe[i];
 538   return MEDCouplingUMesh::MergeUMeshes(ms);
 539 }
 540
 541 /*!
 542  * This method returns a mesh containing as cells that there is patches at the current level.
 543  * The patches are seen like 'boxes' that is too say the refinement will not appear here.
 544  *
 545  * \return MEDCoupling1SGTUMesh * - A new object to be managed by the caller containing as cells as there are patches in \a this.
 546  */
 547 MEDCoupling1SGTUMesh *MEDCouplingCartesianAMRMesh::buildMeshFromPatchEnvelop() const
 548 {
 549   std::vector<const MEDCoupling1SGTUMesh *> cells;
 550   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > cellsSafe;
 551   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 552     {
 553       const MEDCouplingCartesianAMRPatch *patch(*it);
 554       if(patch)
 555         {
 556           MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> cell(patch->getMesh()->getImageMesh()->asSingleCell());
 557           MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> cell1SGT(cell->build1SGTUnstructured());
 558           cellsSafe.push_back(cell1SGT); cells.push_back(cell1SGT);
 559         }
 560     }
 561   return MEDCoupling1SGTUMesh::Merge1SGTUMeshes(cells);
 562 }
 563
 564 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
 565                                                          const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
 566 {
 567   _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
 568 }
 569
 570 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(MEDCouplingCartesianAMRMesh *father, MEDCouplingIMesh *mesh):_father(father)
 571 {
 572   if(!_father)
 573     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : empty father !");
 574   if(!mesh)
 575     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
 576   mesh->checkCoherency();
 577   _mesh=mesh; _mesh->incrRef();
 578 }
 579
 580 void MEDCouplingCartesianAMRMesh::checkPatchId(int patchId) const
 581 {
 582   int sz(getNumberOfPatches());
 583   if(patchId<0 || patchId>=sz)
 584     {
 585       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
 586       throw INTERP_KERNEL::Exception(oss.str().c_str());
 587     }
 588 }
 589
 590 std::size_t MEDCouplingCartesianAMRMesh::getHeapMemorySizeWithoutChildren() const
 591 {
 592   return sizeof(MEDCouplingCartesianAMRMesh);
 593 }
 594
 595 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
 596 {
 597   std::vector<const BigMemoryObject *> ret;
 598   if((const MEDCouplingIMesh *)_mesh)
 599     ret.push_back((const MEDCouplingIMesh *)_mesh);
 600   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 601     {
 602       if((const MEDCouplingCartesianAMRPatch*)*it)
 603         ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
 604     }
 605   return ret;
 606 }
 607
 608 void MEDCouplingCartesianAMRMesh::updateTime() const
 609 {
 610   if((const MEDCouplingIMesh *)_mesh)
 611     updateTimeWith(*_mesh);
 612   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 613     {
 614       const MEDCouplingCartesianAMRPatch *elt(*it);
 615       if(!elt)
 616         continue;
 617       const MEDCouplingCartesianAMRMesh *mesh(elt->getMesh());
 618       if(mesh)
 619         updateTimeWith(*mesh);
 620     }
 621 }