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