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