src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx

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