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 "MEDCouplingIMesh.hxx"
  23 #include "MEDCouplingUMesh.hxx"
  24
  25 #include <limits>
  26 #include <sstream>
  27
  28 using namespace ParaMEDMEM;
  29
  30 /// @cond INTERNAL
  31
  32 /*!
  33  * \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.
  34  * \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,
  35  *                                a the end cell (\b excluded) of the range for the second element of the pair.
  36  * \param [in] factors The refinement per axis relative to the father of \a this.
  37  */
  38 MEDCouplingCartesianAMRPatch::MEDCouplingCartesianAMRPatch(MEDCouplingCartesianAMRMesh *mesh, const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
  39 {
  40   if(!mesh)
  41     throw INTERP_KERNEL::Exception("EDCouplingCartesianAMRPatch constructor : input mesh is NULL !");
  42   _mesh=mesh; _mesh->incrRef();
  43   int dim((int)bottomLeftTopRight.size()),dimExp(_mesh->getSpaceDimension());
  44   if(dim!=dimExp)
  45     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input bottomLeft/topRight size mismatches !");
  46   _bl_tr=bottomLeftTopRight;
  47   if((int)factors.size()!=dimExp)
  48     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRPatch constructor : space dimension of father and input factors per axis size mismatches !");
  49   _factors=factors;
  50 }
  51
  52 int MEDCouplingCartesianAMRPatch::getNumberOfCellsRecursiveWithOverlap() const
  53 {
  54   return _mesh->getNumberOfCellsRecursiveWithOverlap();
  55 }
  56
  57 int MEDCouplingCartesianAMRPatch::getNumberOfCellsRecursiveWithoutOverlap() const
  58 {
  59   return _mesh->getNumberOfCellsRecursiveWithoutOverlap();
  60 }
  61
  62 int MEDCouplingCartesianAMRPatch::getMaxNumberOfLevelsRelativeToThis() const
  63 {
  64   return _mesh->getMaxNumberOfLevelsRelativeToThis();
  65 }
  66
  67 void MEDCouplingCartesianAMRPatch::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
  68 {
  69   return _mesh->addPatch(bottomLeftTopRight,factors);
  70 }
  71
  72 int MEDCouplingCartesianAMRPatch::getNumberOfOverlapedCellsForFather() const
  73 {
  74   return MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(_bl_tr);
  75 }
  76
  77 std::size_t MEDCouplingCartesianAMRPatch::getHeapMemorySizeWithoutChildren() const
  78 {
  79   std::size_t ret(sizeof(MEDCouplingCartesianAMRPatch));
  80   ret+=_bl_tr.capacity()*sizeof(std::pair<int,int>);
  81   return ret;
  82 }
  83
  84 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRPatch::getDirectChildren() const
  85 {
  86   std::vector<const BigMemoryObject *> ret;
  87   if((const MEDCouplingCartesianAMRMesh *)_mesh)
  88     ret.push_back((const MEDCouplingCartesianAMRMesh *)_mesh);
  89   return ret;
  90 }
  91
  92 /// @endcond
  93
  94
  95 MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
  96                                                               const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
  97 {
  98   return new MEDCouplingCartesianAMRMesh(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
  99 }
 100
 101 int MEDCouplingCartesianAMRMesh::getSpaceDimension() const
 102 {
 103   return _mesh->getSpaceDimension();
 104 }
 105
 106 int MEDCouplingCartesianAMRMesh::getMaxNumberOfLevelsRelativeToThis() const
 107 {
 108   int ret(1);
 109   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 110     ret=std::max(ret,(*it)->getMaxNumberOfLevelsRelativeToThis()+1);
 111   return ret;
 112 }
 113
 114 int MEDCouplingCartesianAMRMesh::getNumberOfCellsAtCurrentLevel() const
 115 {
 116   return _mesh->getNumberOfCells();
 117 }
 118
 119 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithOverlap() const
 120 {
 121   int ret(_mesh->getNumberOfCells());
 122   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 123     {
 124       ret+=(*it)->getNumberOfCellsRecursiveWithOverlap();
 125     }
 126   return ret;
 127 }
 128
 129 int MEDCouplingCartesianAMRMesh::getNumberOfCellsRecursiveWithoutOverlap() const
 130 {
 131   int ret(_mesh->getNumberOfCells());
 132   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 133     {
 134       ret-=(*it)->getNumberOfOverlapedCellsForFather();
 135       ret+=(*it)->getNumberOfCellsRecursiveWithoutOverlap();
 136     }
 137   return ret;
 138 }
 139
 140 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getFather() const
 141 {
 142   return _father;
 143 }
 144
 145 const MEDCouplingCartesianAMRMesh *MEDCouplingCartesianAMRMesh::getGodFather() const
 146 {
 147   if(_father==0)
 148     return this;
 149   else
 150     return _father->getGodFather();
 151 }
 152
 153 void MEDCouplingCartesianAMRMesh::detachFromFather()
 154 {
 155   _father=0;
 156 }
 157
 158 /*!
 159  * \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,
 160  *                                a the end cell (\b excluded) of the range for the second element of the pair.
 161  * \param [in] factors The != 0 factor of refinement per axis.
 162  */
 163 void MEDCouplingCartesianAMRMesh::addPatch(const std::vector< std::pair<int,int> >& bottomLeftTopRight, const std::vector<int>& factors)
 164 {
 165   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> mesh(static_cast<MEDCouplingIMesh *>(_mesh->buildStructuredSubPart(bottomLeftTopRight)));
 166   mesh->refineWithFactor(factors);
 167   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRMesh> zeMesh(new MEDCouplingCartesianAMRMesh(this,mesh));
 168   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> elt(new MEDCouplingCartesianAMRPatch(zeMesh,bottomLeftTopRight,factors));
 169   _patches.push_back(elt);
 170 }
 171
 172 /// @cond INTERNAL
 173
 174 class InternalPatch : public RefCountObjectOnly
 175 {
 176 public:
 177   InternalPatch():_nb_of_true(0) { }
 178   double getEfficiency() const { return (double)_nb_of_true/(double)_crit.size(); }
 179   int getNumberOfCells() const { return (int)_crit.size(); }
 180   void setNumberOfTrue(int nboft) { _nb_of_true=nboft; }
 181   std::vector<bool>& getCriterion() { return _crit; }
 182   std::vector< std::pair<int,int> >& getPart() { return _part; }
 183   const std::vector< std::pair<int,int> >& getConstPart() const { return _part; }
 184   bool presenceOfTrue() const { return _nb_of_true>0; }
 185 protected:
 186   ~InternalPatch() { }
 187 private:
 188   int _nb_of_true;
 189   std::vector<bool> _crit;
 190   std::vector< std::pair<int,int> > _part;
 191 };
 192
 193 #if 0
 194 void dissectBigPatch (const Mesh& mesh, const Field& fieldFlag, const unsigned int minCellDirection,
 195                       const unsigned int big_dims, const int dissect_direction, int cut[3] ) const
 196 {
 197   int cut_found = 0 ;
 198   int cut_place = -1 ;
 199   float * ratio = NULL ;
 200   float * ratio_inner = NULL ;
 201
 202   ratio = new float [big_dims-1];
 203   for(unsigned int id=0; id<big_dims-1; id++)
 204     {
 205       float efficiency[2] ;
 206       for(int h=0; h<2; h++)
 207         {
 208           int rect_h[4] ;
 209           copy(getIndexCorners(),getIndexCorners()+4,rect_h) ;
 210           if (h == 0 )
 211             rect_h[dissect_direction+2] = _indexCorners[dissect_direction]+id ;
 212           else if ( h == 1)
 213             rect_h[dissect_direction] =  _indexCorners[dissect_direction]+id+1;
 214
 215           Patch patch_h(rect_h);
 216           patch_h.computeMesh(mesh);
 217           patch_h.computeFieldFlag(fieldFlag);
 218
 219           int nb_cells_h ;
 220           if ( dissect_direction == 0 )
 221             nb_cells_h = patch_h.getNx() ;
 222           else
 223             nb_cells_h = patch_h.getNy() ;
 224
 225           int nb_cells_flag_h = patch_h.getNumberOfCellsFlags();
 226           efficiency[h] = float (nb_cells_flag_h) / float(nb_cells_h) ;
 227         }
 228       ratio[id] = max(efficiency[0],efficiency[1])/
 229           min(efficiency[0],efficiency[1]) ;
 230     }
 231
 232   int dim_ratio_inner = big_dims-1-2*(minCellDirection-1) ;
 233   ratio_inner = new float [dim_ratio_inner];
 234   float min_ratio = 1.E10 ;
 235   int index_min = -1 ;
 236   for(int i=0; i<dim_ratio_inner; i++)
 237     {
 238       if ( ratio[minCellDirection-1+i] < min_ratio )
 239         {
 240           min_ratio = ratio[minCellDirection-1+i] ;
 241           index_min = i+minCellDirection ;
 242         }
 243     }
 244   cut_found = 1 ;
 245   cut_place = index_min + _indexCorners[dissect_direction] - 1 ;
 246   cut[0] = cut_found ;
 247   cut[1] = cut_place ;
 248   cut[2] = dissect_direction ;
 249   delete [] ratio ;
 250   delete [] ratio_inner ;
 251 }
 252 #endif
 253 /// @endcond
 254
 255 /*!
 256  * This method creates patches in \a this (by destroying the patches if any). This method uses \a criterion
 257  */
 258 void MEDCouplingCartesianAMRMesh::createPatchesFromCriterion(const INTERP_KERNEL::BoxSplittingOptions& bso, const DataArrayByte *criterion, const std::vector<int>& factors)
 259 {
 260   if(!criterion || !criterion->isAllocated())
 261     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh::createPatchesFromCriterion : the criterion DataArrayByte instance must be allocated and not NULL !");
 262   int nbCells(getNumberOfCellsAtCurrentLevel());
 263   if(nbCells!=criterion->getNumberOfTuples())
 264     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 !");
 265   _patches.clear();
 266   //
 267   std::vector<int> cgs(_mesh->getCellGridStructure());
 268   std::vector<bool> crit(criterion->toVectorOfBool());//check that criterion has one component.
 269   std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatches,listOfPatchesOK;
 270   //
 271   MEDCouplingAutoRefCountObjectPtr<InternalPatch> p(new InternalPatch);
 272   p->setNumberOfTrue(MEDCouplingStructuredMesh::FindMinimalPartOf(cgs,crit,p->getCriterion(),p->getPart()));
 273   if(p->presenceOfTrue())
 274     listOfPatches.push_back(p);
 275   while(!listOfPatches.empty())
 276     {
 277       std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> > listOfPatchesTmp;
 278       for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::iterator it=listOfPatches.begin();it!=listOfPatches.end();it++)
 279         {
 280           //
 281           if((*it)->getEfficiency()>=bso.getEffeciency())
 282             {
 283               if((*it)->getNumberOfCells()>=bso.getMaxCells())
 284                 {
 285
 286                 }
 287             }
 288         }
 289       listOfPatches=listOfPatchesTmp;
 290     }
 291   for(std::vector< MEDCouplingAutoRefCountObjectPtr<InternalPatch> >::const_iterator it=listOfPatchesOK.begin();it!=listOfPatchesOK.end();it++)
 292     {
 293       addPatch((*it)->getConstPart(),factors);
 294     }
 295 }
 296
 297 void MEDCouplingCartesianAMRMesh::removePatch(int patchId)
 298 {
 299   checkPatchId(patchId);
 300   int sz((int)_patches.size()),j(0);
 301   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> > patches(sz-1);
 302   for(int i=0;i<sz;i++)
 303     if(i!=patchId)
 304       patches[j++]=_patches[i];
 305   (const_cast<MEDCouplingCartesianAMRMesh *>(_patches[patchId]->getMesh()))->detachFromFather();
 306   _patches=patches;
 307   declareAsNew();
 308 }
 309
 310 int MEDCouplingCartesianAMRMesh::getNumberOfPatches() const
 311 {
 312   return (int)_patches.size();
 313 }
 314
 315 const MEDCouplingCartesianAMRPatch *MEDCouplingCartesianAMRMesh::getPatch(int patchId) const
 316 {
 317   checkPatchId(patchId);
 318   return _patches[patchId];
 319 }
 320
 321 MEDCouplingUMesh *MEDCouplingCartesianAMRMesh::buildUnstructured() const
 322 {
 323   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> part(_mesh->buildUnstructured());
 324   std::vector<bool> bs(_mesh->getNumberOfCells(),false);
 325   std::vector<int> cgs(_mesh->getCellGridStructure());
 326   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > msSafe(_patches.size()+1);
 327   std::size_t ii(0);
 328   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++,ii++)
 329     {
 330       MEDCouplingStructuredMesh::SwitchOnIdsFrom(cgs,(*it)->getBLTRRange(),bs);
 331       msSafe[ii+1]=(*it)->getMesh()->buildUnstructured();
 332     }
 333   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsOff(DataArrayInt::BuildListOfSwitchedOff(bs));
 334   msSafe[0]=static_cast<MEDCouplingUMesh *>(part->buildPartOfMySelf(eltsOff->begin(),eltsOff->end(),false));
 335   std::vector< const MEDCouplingUMesh * > ms(msSafe.size());
 336   for(std::size_t i=0;i<msSafe.size();i++)
 337     ms[i]=msSafe[i];
 338   return MEDCouplingUMesh::MergeUMeshes(ms);
 339 }
 340
 341 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
 342                                                          const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop):_father(0)
 343 {
 344   _mesh=MEDCouplingIMesh::New(meshName,spaceDim,nodeStrctStart,nodeStrctStop,originStart,originStop,dxyzStart,dxyzStop);
 345 }
 346
 347 MEDCouplingCartesianAMRMesh::MEDCouplingCartesianAMRMesh(MEDCouplingCartesianAMRMesh *father, MEDCouplingIMesh *mesh):_father(father)
 348 {
 349   if(!_father)
 350     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : empty father !");
 351   if(!mesh)
 352     throw INTERP_KERNEL::Exception("MEDCouplingCartesianAMRMesh(MEDCouplingIMesh *mesh) constructor : The input mesh is null !");
 353   mesh->checkCoherency();
 354   _mesh=mesh; _mesh->incrRef();
 355 }
 356
 357 void MEDCouplingCartesianAMRMesh::checkPatchId(int patchId) const
 358 {
 359   int sz(getNumberOfPatches());
 360   if(patchId<0 || patchId>=sz)
 361     {
 362       std::ostringstream oss; oss << "MEDCouplingCartesianAMRMesh::checkPatchId : invalid patchId (" << patchId << ") ! Must be in [0," << sz << ") !";
 363       throw INTERP_KERNEL::Exception(oss.str().c_str());
 364     }
 365 }
 366
 367 std::size_t MEDCouplingCartesianAMRMesh::getHeapMemorySizeWithoutChildren() const
 368 {
 369   return sizeof(MEDCouplingCartesianAMRMesh);
 370 }
 371
 372 std::vector<const BigMemoryObject *> MEDCouplingCartesianAMRMesh::getDirectChildren() const
 373 {
 374   std::vector<const BigMemoryObject *> ret;
 375   if((const MEDCouplingIMesh *)_mesh)
 376     ret.push_back((const MEDCouplingIMesh *)_mesh);
 377   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 378     {
 379       if((const MEDCouplingCartesianAMRPatch*)*it)
 380         ret.push_back((const MEDCouplingCartesianAMRPatch*)*it);
 381     }
 382   return ret;
 383 }
 384
 385 void MEDCouplingCartesianAMRMesh::updateTime() const
 386 {
 387   if((const MEDCouplingIMesh *)_mesh)
 388     updateTimeWith(*_mesh);
 389   for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingCartesianAMRPatch> >::const_iterator it=_patches.begin();it!=_patches.end();it++)
 390     {
 391       const MEDCouplingCartesianAMRPatch *elt(*it);
 392       if(!elt)
 393         continue;
 394       const MEDCouplingCartesianAMRMesh *mesh(elt->getMesh());
 395       if(mesh)
 396         updateTimeWith(*mesh);
 397     }
 398 }