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