src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx

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