src/MEDCoupling/MEDCouplingCartesianAMRMesh.cxx

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