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