src/MEDCoupling/MEDCouplingStructuredMesh.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 (CEA/DEN)
  20
  21 #include "MEDCouplingStructuredMesh.hxx"
  22 #include "MEDCouplingFieldDouble.hxx"
  23 #include "MEDCouplingMemArray.hxx"
  24 #include "MEDCoupling1GTUMesh.hxx"
  25 #include "MEDCouplingUMesh.hxx"
  26 #include "MEDCouplingIMesh.hxx"//tony to throw when optimization will be performed in AssignPartOfFieldOfDoubleUsing
  27
  28 #include <numeric>
  29
  30 using namespace ParaMEDMEM;
  31
  32 MEDCouplingStructuredMesh::MEDCouplingStructuredMesh()
  33 {
  34 }
  35
  36 MEDCouplingStructuredMesh::MEDCouplingStructuredMesh(const MEDCouplingStructuredMesh& other, bool deepCopy):MEDCouplingMesh(other)
  37 {
  38 }
  39
  40 MEDCouplingStructuredMesh::~MEDCouplingStructuredMesh()
  41 {
  42 }
  43
  44 std::size_t MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren() const
  45 {
  46   return MEDCouplingMesh::getHeapMemorySizeWithoutChildren();
  47 }
  48
  49 void MEDCouplingStructuredMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
  50 {
  51   MEDCouplingMesh::copyTinyStringsFrom(other);
  52 }
  53
  54 bool MEDCouplingStructuredMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
  55 {
  56   return MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason);
  57 }
  58
  59 INTERP_KERNEL::NormalizedCellType MEDCouplingStructuredMesh::getTypeOfCell(int cellId) const
  60 {
  61   return GetGeoTypeGivenMeshDimension(getMeshDimension());
  62 }
  63
  64 INTERP_KERNEL::NormalizedCellType MEDCouplingStructuredMesh::GetGeoTypeGivenMeshDimension(int meshDim)
  65 {
  66   switch(meshDim)
  67   {
  68     case 3:
  69       return INTERP_KERNEL::NORM_HEXA8;
  70     case 2:
  71       return INTERP_KERNEL::NORM_QUAD4;
  72     case 1:
  73       return INTERP_KERNEL::NORM_SEG2;
  74     case 0:
  75       return INTERP_KERNEL::NORM_POINT1;
  76     default:
  77       throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingStructuredMesh::GetGeoTypeGivenMeshDimension !");
  78   }
  79 }
  80
  81 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingStructuredMesh::getAllGeoTypes() const
  82 {
  83   std::set<INTERP_KERNEL::NormalizedCellType> ret2;
  84   ret2.insert(getTypeOfCell(0));
  85   return ret2;
  86 }
  87
  88 int MEDCouplingStructuredMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
  89 {
  90   int ret=getNumberOfCells();
  91   if(type==getTypeOfCell(0))
  92     return ret;
  93   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(getTypeOfCell(0));
  94   std::ostringstream oss; oss << "MEDCouplingStructuredMesh::getNumberOfCellsWithType : no specified type ! Type available is " << cm.getRepr() << " !";
  95   throw INTERP_KERNEL::Exception(oss.str().c_str());
  96 }
  97
  98 DataArrayInt *MEDCouplingStructuredMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
  99 {
 100   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 101   if(getTypeOfCell(0)==type)
 102     {
 103       ret->alloc(getNumberOfCells(),1);
 104       ret->iota(0);
 105     }
 106   else
 107     ret->alloc(0,1);
 108   return ret.retn();
 109 }
 110
 111 DataArrayInt *MEDCouplingStructuredMesh::computeNbOfNodesPerCell() const
 112 {
 113   int nbCells=getNumberOfCells();
 114   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 115   ret->alloc(nbCells,1);
 116   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(getTypeOfCell(0));
 117   ret->fillWithValue((int)cm.getNumberOfNodes());
 118   return ret.retn();
 119 }
 120
 121 DataArrayInt *MEDCouplingStructuredMesh::computeNbOfFacesPerCell() const
 122 {
 123   int nbCells=getNumberOfCells();
 124   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 125   ret->alloc(nbCells,1);
 126   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(getTypeOfCell(0));
 127   ret->fillWithValue((int)cm.getNumberOfSons());
 128   return ret.retn();
 129 }
 130
 131 /*!
 132  * This method computes effective number of nodes per cell. That is to say nodes appearing several times in nodal connectivity of a cell,
 133  * will be counted only once here whereas it will be counted several times in MEDCouplingMesh::computeNbOfNodesPerCell method.
 134  * Here for structured mesh it returns exactly as MEDCouplingStructuredMesh::computeNbOfNodesPerCell does.
 135  *
 136  * \return DataArrayInt * - new object to be deallocated by the caller.
 137  */
 138 DataArrayInt *MEDCouplingStructuredMesh::computeEffectiveNbOfNodesPerCell() const
 139 {
 140   return computeNbOfNodesPerCell();
 141 }
 142
 143 void MEDCouplingStructuredMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
 144 {
 145   int meshDim=getMeshDimension();
 146   int tmpCell[3],tmpNode[3];
 147   getSplitCellValues(tmpCell);
 148   getSplitNodeValues(tmpNode);
 149   int tmp2[3];
 150   GetPosFromId(cellId,meshDim,tmpCell,tmp2);
 151   switch(meshDim)
 152   {
 153     case 1:
 154       conn.push_back(tmp2[0]); conn.push_back(tmp2[0]+1);
 155       break;
 156     case 2:
 157       conn.push_back(tmp2[1]*tmpNode[1]+tmp2[0]); conn.push_back(tmp2[1]*tmpNode[1]+tmp2[0]+1);
 158       conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+1); conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]);
 159       break;
 160     case 3:
 161       conn.push_back(tmp2[1]*tmpNode[1]+tmp2[0]+tmp2[2]*tmpNode[2]); conn.push_back(tmp2[1]*tmpNode[1]+tmp2[0]+1+tmp2[2]*tmpNode[2]);
 162       conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+1+tmp2[2]*tmpNode[2]); conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+tmp2[2]*tmpNode[2]);
 163       conn.push_back(tmp2[1]*tmpNode[1]+tmp2[0]+(tmp2[2]+1)*tmpNode[2]); conn.push_back(tmp2[1]*tmpNode[1]+tmp2[0]+1+(tmp2[2]+1)*tmpNode[2]);
 164       conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+1+(tmp2[2]+1)*tmpNode[2]); conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+(tmp2[2]+1)*tmpNode[2]);
 165       break;
 166     default:
 167       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::getNodeIdsOfCell : big problem spacedim must be in 1,2 or 3 !");
 168   };
 169 }
 170
 171 /*!
 172  * This method returns the mesh dimension of \a this. It can be different from space dimension in case of a not null dimension contains only one node.
 173  */
 174 int MEDCouplingStructuredMesh::getMeshDimension() const
 175 {
 176   std::vector<int> ngs(getNodeGridStructure());
 177   int ret(0),pos(0);
 178   for(std::vector<int>::const_iterator it=ngs.begin();it!=ngs.end();it++,pos++)
 179     {
 180       if(*it<=0)
 181         {
 182           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::getMeshDimension : At pos #" << pos << " number of nodes is " << *it << " ! Must be > 0 !";
 183           throw INTERP_KERNEL::Exception(oss.str().c_str());
 184         }
 185       if(*it>1)
 186         ret++;
 187     }
 188   return ret;
 189 }
 190
 191 /*!
 192  * This method returns the space dimension by only considering the node grid structure.
 193  * For cartesian mesh the returned value is equal to those returned by getSpaceDimension.
 194  * But for curvelinear is could be different !
 195  */
 196 int MEDCouplingStructuredMesh::getSpaceDimensionOnNodeStruct() const
 197 {
 198   std::vector<int> nodeStr(getNodeGridStructure());
 199   int spd1(0),pos(0);
 200   for(std::vector<int>::const_iterator it=nodeStr.begin();it!=nodeStr.end();it++,pos++)
 201     {
 202       int elt(*it);
 203       if(elt<=0)
 204         {
 205           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::getSpaceDimensionOnNodeStruct : At pos #" << pos << " value of node grid structure is " << *it << " ! must be >=1 !";
 206           throw INTERP_KERNEL::Exception(oss.str().c_str());
 207         }
 208       spd1++;
 209     }
 210   return spd1;
 211 }
 212
 213 void MEDCouplingStructuredMesh::getSplitCellValues(int *res) const
 214 {
 215   std::vector<int> strct(getCellGridStructure());
 216   std::vector<int> ret(MEDCouplingStructuredMesh::GetSplitVectFromStruct(strct));
 217   std::copy(ret.begin(),ret.end(),res);
 218 }
 219
 220 void MEDCouplingStructuredMesh::getSplitNodeValues(int *res) const
 221 {
 222   std::vector<int> strct(getNodeGridStructure());
 223   std::vector<int> ret(MEDCouplingStructuredMesh::GetSplitVectFromStruct(strct));
 224   std::copy(ret.begin(),ret.end(),res);
 225 }
 226
 227 /*!
 228  * This method returns the number of cells of unstructured sub level mesh, without building it.
 229  */
 230 int MEDCouplingStructuredMesh::getNumberOfCellsOfSubLevelMesh() const
 231 {
 232   std::vector<int> cgs(getCellGridStructure());
 233   return GetNumberOfCellsOfSubLevelMesh(cgs,getMeshDimension());
 234 }
 235
 236 /*!
 237  * See MEDCouplingUMesh::getDistributionOfTypes for more information
 238  */
 239 std::vector<int> MEDCouplingStructuredMesh::getDistributionOfTypes() const
 240 {
 241   //only one type of cell
 242   std::vector<int> ret(3);
 243   ret[0]=getTypeOfCell(0);
 244   ret[1]=getNumberOfCells();
 245   ret[2]=-1; //ret[3*k+2]==-1 because it has no sense here
 246   return ret;
 247 }
 248
 249 /*!
 250  * This method tries to minimize at most the number of deep copy.
 251  * So if \a idsPerType is not empty it can be returned directly (without copy, but with ref count incremented) in return.
 252  *
 253  * See MEDCouplingUMesh::checkTypeConsistencyAndContig for more information
 254  */
 255 DataArrayInt *MEDCouplingStructuredMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const
 256 {
 257   int nbOfCells=getNumberOfCells();
 258   if(code.size()!=3)
 259     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::checkTypeConsistencyAndContig : invalid input code should be exactly of size 3 !");
 260   if(code[0]!=(int)getTypeOfCell(0))
 261     {
 262       std::ostringstream oss; oss << "MEDCouplingStructuredMesh::checkTypeConsistencyAndContig : Mismatch of geometric type ! Asking for " << code[0] << " whereas the geometric type is \a this is " << getTypeOfCell(0) << " !";
 263       throw INTERP_KERNEL::Exception(oss.str().c_str());
 264     }
 265   if(code[2]==-1)
 266     {
 267       if(code[1]==nbOfCells)
 268         return 0;
 269       else
 270         {
 271           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::checkTypeConsistencyAndContig : mismatch between the number of cells in this (" << nbOfCells << ") and the number of non profile (" << code[1] << ") !";
 272           throw INTERP_KERNEL::Exception(oss.str().c_str());
 273         }
 274     }
 275   if(code[2]!=0)
 276     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::checkTypeConsistencyAndContig : single geo type mesh ! 0 or -1 is expected at pos #2 of input code !");
 277   if(idsPerType.size()!=1)
 278     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::checkTypeConsistencyAndContig : input code points to DataArrayInt #0 whereas the size of idsPerType is not equal to 1 !");
 279   const DataArrayInt *pfl=idsPerType[0];
 280   if(!pfl)
 281     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::checkTypeConsistencyAndContig : the input code points to a NULL DataArrayInt at rank 0 !");
 282   if(pfl->getNumberOfComponents()!=1)
 283     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::checkTypeConsistencyAndContig : input profile should have exactly one component !");
 284   pfl->checkAllIdsInRange(0,nbOfCells);
 285   pfl->incrRef();
 286   return const_cast<DataArrayInt *>(pfl);
 287 }
 288
 289 /*!
 290  * This method is the opposite of MEDCouplingUMesh::checkTypeConsistencyAndContig method. Given a list of cells in \a profile it returns a list of sub-profiles sorted by geo type.
 291  * The result is put in the array \a idsPerType. In the returned parameter \a code, foreach i \a code[3*i+2] refers (if different from -1) to a location into the \a idsPerType.
 292  * This method has 1 input \a profile and 3 outputs \a code \a idsInPflPerType and \a idsPerType.
 293  *
 294  * \param [out] code is a vector of size 3*n where n is the number of different geometric type in \a this \b reduced to the profile \a profile. \a code has exactly the same semantic than in MEDCouplingUMesh::checkTypeConsistencyAndContig method.
 295  * \param [out] idsInPflPerType is a vector of size of different geometric type in the subpart defined by \a profile of \a this ( equal to \a code.size()/3). For each i,
 296  *              \a idsInPflPerType[i] stores the tuple ids in \a profile that correspond to the geometric type code[3*i+0]
 297  * \param [out] idsPerType is a vector of size of different sub profiles needed to be defined to represent the profile \a profile for a given geometric type.
 298  *              This vector can be empty in case of all geometric type cells are fully covered in ascending in the given input \a profile.
 299  *
 300  * \warning for performance reasons no deep copy will be performed, if \a profile can been used as this in output parameters \a idsInPflPerType and \a idsPerType.
 301  *
 302  * \throw if \a profile has not exactly one component. It throws too, if \a profile contains some values not in [0,getNumberOfCells()) or if \a this is not fully defined
 303  *
 304  *  \b Example1: <br>
 305  *          - Before \a this has 3 cells \a profile contains [0,1,2]
 306  *          - After \a code contains [NORM_...,nbCells,-1], \a idsInPflPerType [[0,1,2]] and \a idsPerType is empty <br>
 307  *
 308  *  \b Example2: <br>
 309  *          - Before \a this has 3 cells \a profile contains [1,2]
 310  *          - After \a code contains [NORM_...,nbCells,0], \a idsInPflPerType [[0,1]] and \a idsPerType is [[1,2]] <br>
 311
 312  */
 313 void MEDCouplingStructuredMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const
 314 {
 315   if(!profile || !profile->isAllocated())
 316     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::splitProfilePerType : input profile is NULL or not allocated !");
 317   if(profile->getNumberOfComponents()!=1)
 318     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::splitProfilePerType : input profile should have exactly one component !");
 319   int nbTuples=profile->getNumberOfTuples();
 320   int nbOfCells=getNumberOfCells();
 321   code.resize(3); idsInPflPerType.resize(1);
 322   code[0]=(int)getTypeOfCell(0); code[1]=nbOfCells;
 323   idsInPflPerType.resize(1);
 324   if(profile->isIdentity() && nbTuples==nbOfCells)
 325     {
 326       code[2]=-1;
 327       idsInPflPerType[0]=0;
 328       idsPerType.clear();
 329       return ;
 330     }
 331   code[1]=profile->getNumberOfTuples();
 332   code[2]=0;
 333   profile->checkAllIdsInRange(0,nbOfCells);
 334   idsPerType.resize(1);
 335   idsPerType[0]=profile->deepCpy();
 336   idsInPflPerType[0]=DataArrayInt::Range(0,nbTuples,1);
 337 }
 338
 339 /*!
 340  * Creates a new unstructured mesh (MEDCoupling1SGTUMesh) from \a this structured one.
 341  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingUMesh. The caller is to
 342  * delete this array using decrRef() as it is no more needed.
 343  *  \throw If \a this->getMeshDimension() is not among [1,2,3].
 344  */
 345 MEDCoupling1SGTUMesh *MEDCouplingStructuredMesh::build1SGTUnstructured() const
 346 {
 347   int meshDim(getMeshDimension()),spaceDim(getSpaceDimensionOnNodeStruct());
 348   if((meshDim<0 || meshDim>3) || (spaceDim<0 || spaceDim>3))
 349     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::build1SGTUnstructured : meshdim and spacedim must be in [1,2,3] !");
 350   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords(getCoordinatesAndOwner());
 351   int ns[3];
 352   getNodeGridStructure(ns);
 353   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(Build1GTNodalConnectivity(ns,ns+spaceDim));
 354   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName(),GetGeoTypeGivenMeshDimension(meshDim)));
 355   ret->setNodalConnectivity(conn); ret->setCoords(coords);
 356   try
 357     { ret->copyTinyInfoFrom(this); }
 358   catch(INTERP_KERNEL::Exception&) { }
 359   return ret.retn();
 360 }
 361
 362 /*!
 363  * This method returns the unstructured mesh (having single geometric type) of the sub level mesh of \a this.
 364  * This method is equivalent to computing MEDCouplingUMesh::buildDescendingConnectivity on the unstructurized \a this mesh.
 365  *
 366  * The caller is to delete the returned mesh using decrRef() as it is no more needed.
 367  */
 368 MEDCoupling1SGTUMesh *MEDCouplingStructuredMesh::build1SGTSubLevelMesh() const
 369 {
 370   int meshDim(getMeshDimension());
 371   if(meshDim<1 || meshDim>3)
 372     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::build1SGTSubLevelMesh : meshdim must be in [2,3] !");
 373   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords(getCoordinatesAndOwner());
 374   int ns[3];
 375   getNodeGridStructure(ns);
 376   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(Build1GTNodalConnectivityOfSubLevelMesh(ns,ns+meshDim));
 377   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName(),GetGeoTypeGivenMeshDimension(meshDim-1)));
 378   ret->setNodalConnectivity(conn); ret->setCoords(coords);
 379   return ret.retn();
 380 }
 381
 382 /*!
 383  * Creates a new unstructured mesh (MEDCouplingUMesh) from \a this structured one.
 384  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingUMesh. The caller is to
 385  * delete this array using decrRef() as it is no more needed.
 386  *  \throw If \a this->getMeshDimension() is not among [1,2,3].
 387  */
 388 MEDCouplingUMesh *MEDCouplingStructuredMesh::buildUnstructured() const
 389 {
 390   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret0(build1SGTUnstructured());
 391   return ret0->buildUnstructured();
 392 }
 393
 394 /*!
 395  * Creates a new MEDCouplingUMesh containing a part of cells of \a this mesh.
 396  * The cells to include to the
 397  * result mesh are specified by an array of cell ids.
 398  *  \param [in] start - an array of cell ids to include to the result mesh.
 399  *  \param [in] end - specifies the end of the array \a start, so that
 400  *              the last value of \a start is \a end[ -1 ].
 401  *  \return MEDCouplingMesh * - a new instance of MEDCouplingUMesh. The caller is to
 402  *         delete this mesh using decrRef() as it is no more needed.
 403  */
 404 MEDCouplingMesh *MEDCouplingStructuredMesh::buildPart(const int *start, const int *end) const
 405 {
 406   MEDCouplingUMesh *um=buildUnstructured();
 407   MEDCouplingMesh *ret=um->buildPart(start,end);
 408   um->decrRef();
 409   return ret;
 410 }
 411
 412 MEDCouplingMesh *MEDCouplingStructuredMesh::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
 413 {
 414   std::vector<int> cgs(getCellGridStructure());
 415   std::vector< std::pair<int,int> > cellPartFormat,nodePartFormat;
 416   if(IsPartStructured(start,end,cgs,cellPartFormat))
 417     {
 418       MEDCouplingAutoRefCountObjectPtr<MEDCouplingStructuredMesh> ret(buildStructuredSubPart(cellPartFormat));
 419       nodePartFormat=cellPartFormat;
 420       for(std::vector< std::pair<int,int> >::iterator it=nodePartFormat.begin();it!=nodePartFormat.end();it++)
 421         (*it).second++;
 422       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1(BuildExplicitIdsFrom(getNodeGridStructure(),nodePartFormat));
 423       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp2(DataArrayInt::New()); tmp2->alloc(getNumberOfNodes(),1);
 424       tmp2->fillWithValue(-1);
 425       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp3(DataArrayInt::New()); tmp3->alloc(tmp1->getNumberOfTuples(),1); tmp3->iota(0);
 426       tmp2->setPartOfValues3(tmp3,tmp1->begin(),tmp1->end(),0,1,1);
 427       arr=tmp2.retn();
 428       return ret.retn();
 429     }
 430   else
 431     {
 432       MEDCouplingUMesh *um=buildUnstructured();
 433       MEDCouplingMesh *ret=um->buildPartAndReduceNodes(start,end,arr);
 434       um->decrRef();
 435       return ret;
 436     }
 437 }
 438
 439 DataArrayInt *MEDCouplingStructuredMesh::simplexize(int policy)
 440 {
 441   throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::simplexize : not available for Cartesian mesh !");
 442 }
 443
 444 /*!
 445  * Returns a new MEDCouplingFieldDouble holding normal vectors to cells of \a this
 446  * 2D mesh. The computed vectors have 3 components and are normalized.
 447  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on
 448  *          cells and one time. The caller is to delete this field using decrRef() as
 449  *          it is no more needed.
 450  *  \throw If \a this->getMeshDimension() != 2.
 451  */
 452 MEDCouplingFieldDouble *MEDCouplingStructuredMesh::buildOrthogonalField() const
 453 {
 454   if(getMeshDimension()!=2)
 455     throw INTERP_KERNEL::Exception("Expected a MEDCouplingStructuredMesh with meshDim == 2 !");
 456   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 457   DataArrayDouble *array=DataArrayDouble::New();
 458   int nbOfCells=getNumberOfCells();
 459   array->alloc(nbOfCells,3);
 460   double *vals=array->getPointer();
 461   for(int i=0;i<nbOfCells;i++)
 462     { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
 463   ret->setArray(array);
 464   array->decrRef();
 465   ret->setMesh(this);
 466   return ret;
 467 }
 468
 469 void MEDCouplingStructuredMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const
 470 {
 471   std::vector<int> ngs(getNodeGridStructure());
 472   int dim(getSpaceDimension());
 473   switch(dim)
 474   {
 475     case 1:
 476       return GetReverseNodalConnectivity1(ngs,revNodal,revNodalIndx);
 477     case 2:
 478       return GetReverseNodalConnectivity2(ngs,revNodal,revNodalIndx);
 479     case 3:
 480       return GetReverseNodalConnectivity3(ngs,revNodal,revNodalIndx);
 481     default:
 482       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::getReverseNodalConnectivity : only dimensions 1, 2 and 3 are supported !");
 483   }
 484 }
 485
 486 void MEDCouplingStructuredMesh::GetReverseNodalConnectivity1(const std::vector<int>& ngs, DataArrayInt *revNodal, DataArrayInt *revNodalIndx)
 487 {
 488   int nbNodes(ngs[0]);
 489   revNodalIndx->alloc(nbNodes+1,1);
 490   if(nbNodes==0)
 491     { revNodal->alloc(0,1); revNodalIndx->setIJ(0,0,0); return ; }
 492   if(nbNodes==1)
 493     { revNodal->alloc(1,1); revNodal->setIJ(0,0,0); revNodalIndx->setIJ(0,0,0); revNodalIndx->setIJ(1,0,1); return ; }
 494   revNodal->alloc(2*(nbNodes-1),1);
 495   int *rn(revNodal->getPointer()),*rni(revNodalIndx->getPointer());
 496   *rni++=0; *rni=1; *rn++=0;
 497   for(int i=1;i<nbNodes-1;i++,rni++)
 498     {
 499       rn[0]=i-1; rn[1]=i;
 500       rni[1]=rni[0]+2;
 501       rn+=2;
 502     }
 503   rn[0]=nbNodes-2; rni[1]=rni[0]+1;
 504 }
 505
 506 void MEDCouplingStructuredMesh::GetReverseNodalConnectivity2(const std::vector<int>& ngs, DataArrayInt *revNodal, DataArrayInt *revNodalIndx)
 507 {
 508   int nbNodesX(ngs[0]),nbNodesY(ngs[1]);
 509   int nbNodes(nbNodesX*nbNodesY);
 510   if(nbNodesX==0 || nbNodesY==0)
 511     { revNodal->alloc(0,1); revNodalIndx->setIJ(0,0,0); return ; }
 512   if(nbNodesX==1 || nbNodesY==1)
 513     { std::vector<int> ngs2(1); ngs2[0]=std::max(nbNodesX,nbNodesY); return GetReverseNodalConnectivity1(ngs2,revNodal,revNodalIndx); }
 514   revNodalIndx->alloc(nbNodes+1,1);
 515   int nbCellsX(nbNodesX-1),nbCellsY(nbNodesY-1);
 516   revNodal->alloc(4*(nbNodesX-2)*(nbNodesY-2)+2*2*(nbNodesX-2)+2*2*(nbNodesY-2)+4,1);
 517   int *rn(revNodal->getPointer()),*rni(revNodalIndx->getPointer());
 518   *rni++=0; *rni=1; *rn++=0;
 519   for(int i=1;i<nbNodesX-1;i++,rni++,rn+=2)
 520     {
 521       rn[0]=i-1; rn[1]=i;
 522       rni[1]=rni[0]+2;
 523     }
 524   rni[1]=rni[0]+1; *rn++=nbCellsX-1;
 525   rni++;
 526   for(int j=1;j<nbNodesY-1;j++)
 527     {
 528       int off(nbCellsX*(j-1)),off2(nbCellsX*j);
 529       rni[1]=rni[0]+2; rn[0]=off; rn[1]=off2;
 530       rni++; rn+=2;
 531       for(int i=1;i<nbNodesX-1;i++,rni++,rn+=4)
 532         {
 533           rn[0]=i-1+off; rn[1]=i+off; rn[2]=i-1+off2; rn[3]=i+off2;
 534           rni[1]=rni[0]+4;
 535         }
 536       rni[1]=rni[0]+2; rn[0]=off+nbCellsX-1; rn[1]=off2+nbCellsX-1;
 537       rni++; rn+=2;
 538     }
 539   int off3(nbCellsX*(nbCellsY-1));
 540   rni[1]=rni[0]+1;
 541   rni++; *rn++=off3;
 542   for(int i=1;i<nbNodesX-1;i++,rni++,rn+=2)
 543     {
 544       rn[0]=i-1+off3; rn[1]=i+off3;
 545       rni[1]=rni[0]+2;
 546     }
 547   rni[1]=rni[0]+1; rn[0]=nbCellsX*nbCellsY-1;
 548 }
 549
 550 void MEDCouplingStructuredMesh::GetReverseNodalConnectivity3(const std::vector<int>& ngs, DataArrayInt *revNodal, DataArrayInt *revNodalIndx)
 551 {
 552   int nbNodesX(ngs[0]),nbNodesY(ngs[1]),nbNodesZ(ngs[2]);
 553   int nbNodes(nbNodesX*nbNodesY*nbNodesZ);
 554   if(nbNodesX==0 || nbNodesY==0 || nbNodesZ==0)
 555     { revNodal->alloc(0,1); revNodalIndx->setIJ(0,0,0); return ; }
 556   if(nbNodesX==1 || nbNodesY==1 || nbNodesZ==1)
 557     {
 558       std::vector<int> ngs2(2);
 559       int pos(0);
 560       bool pass(false);
 561       for(int i=0;i<3;i++)
 562         {
 563           if(pass)
 564             { ngs2[pos++]=ngs[i]; }
 565           else
 566             {
 567               pass=ngs[i]==1;
 568               if(!pass)
 569                 { ngs2[pos++]=ngs[i]; }
 570             }
 571         }
 572       return GetReverseNodalConnectivity2(ngs2,revNodal,revNodalIndx);
 573     }
 574   revNodalIndx->alloc(nbNodes+1,1);
 575   int nbCellsX(nbNodesX-1),nbCellsY(nbNodesY-1),nbCellsZ(nbNodesZ-1);
 576   revNodal->alloc(8*(nbNodesX-2)*(nbNodesY-2)*(nbNodesZ-2)+4*(2*(nbNodesX-2)*(nbNodesY-2)+2*(nbNodesX-2)*(nbNodesZ-2)+2*(nbNodesY-2)*(nbNodesZ-2))+2*4*(nbNodesX-2)+2*4*(nbNodesY-2)+2*4*(nbNodesZ-2)+8,1);
 577   int *rn(revNodal->getPointer()),*rni(revNodalIndx->getPointer());
 578   *rni=0;
 579   for(int k=0;k<nbNodesZ;k++)
 580     {
 581       bool factZ(k!=0 && k!=nbNodesZ-1);
 582       int offZ0((k-1)*nbCellsX*nbCellsY),offZ1(k*nbCellsX*nbCellsY);
 583       for(int j=0;j<nbNodesY;j++)
 584         {
 585           bool factYZ(factZ && (j!=0 && j!=nbNodesY-1));
 586           int off00((j-1)*nbCellsX+offZ0),off01(j*nbCellsX+offZ0),off10((j-1)*nbCellsX+offZ1),off11(j*nbCellsX+offZ1);
 587           for(int i=0;i<nbNodesX;i++,rni++)
 588             {
 589               int fact(factYZ && (i!=0 && i!=nbNodesX-1));
 590               if(fact)
 591                 {//most of points fall in this part of code
 592                   rn[0]=off00+i-1; rn[1]=off00+i; rn[2]=off01+i-1; rn[3]=off01+i;
 593                   rn[4]=off10+i-1; rn[5]=off10+i; rn[6]=off11+i-1; rn[7]=off11+i;
 594                   rni[1]=rni[0]+8;
 595                   rn+=8;
 596                 }
 597               else
 598                 {
 599                   int *rnRef(rn);
 600                   if(k>=1 && j>=1 && i>=1)
 601                     *rn++=off00+i-1;
 602                   if(k>=1 && j>=1 && i<nbCellsX)
 603                     *rn++=off00+i;
 604                   if(k>=1 && j<nbCellsY && i>=1)
 605                     *rn++=off01+i-1;
 606                   if(k>=1 && j<nbCellsY && i<nbCellsX)
 607                     *rn++=off01+i;
 608                   //
 609                   if(k<nbCellsZ && j>=1 && i>=1)
 610                     *rn++=off10+i-1;
 611                   if(k<nbCellsZ && j>=1 && i<nbCellsX)
 612                     *rn++=off10+i;
 613                   if(k<nbCellsZ && j<nbCellsY && i>=1)
 614                     *rn++=off11+i-1;
 615                   if(k<nbCellsZ && j<nbCellsY && i<nbCellsX)
 616                     *rn++=off11+i;
 617                   rni[1]=rni[0]+(int)(std::distance(rnRef,rn));
 618                 }
 619             }
 620         }
 621     }
 622 }
 623
 624 /*!
 625  * \return DataArrayInt * - newly allocated instance of nodal connectivity compatible for MEDCoupling1SGTMesh instance
 626  */
 627 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity(const int *nodeStBg, const int *nodeStEnd)
 628 {
 629   int zippedNodeSt[3];
 630   int dim(ZipNodeStructure(nodeStBg,nodeStEnd,zippedNodeSt));
 631   switch(dim)
 632   {
 633     case 0:
 634       {
 635         MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
 636         conn->alloc(1,1); conn->setIJ(0,0,0);
 637         return conn.retn();
 638       }
 639     case 1:
 640       return Build1GTNodalConnectivity1D(zippedNodeSt);
 641     case 2:
 642       return Build1GTNodalConnectivity2D(zippedNodeSt);
 643     case 3:
 644       return Build1GTNodalConnectivity3D(zippedNodeSt);
 645     default:
 646       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::Build1GTNodalConnectivity : only dimension in [0,1,2,3] supported !");
 647   }
 648 }
 649
 650 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivityOfSubLevelMesh(const int *nodeStBg, const int *nodeStEnd)
 651 {
 652   std::size_t dim(std::distance(nodeStBg,nodeStEnd));
 653   switch(dim)
 654   {
 655     case 3:
 656       return Build1GTNodalConnectivityOfSubLevelMesh3D(nodeStBg);
 657     case 2:
 658       return Build1GTNodalConnectivityOfSubLevelMesh2D(nodeStBg);
 659     default:
 660       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::Build1GTNodalConnectivityOfSubLevelMesh: only dimension in [2,3] supported !");
 661   }
 662 }
 663
 664 /*!
 665  * This method returns the list of ids sorted ascendingly of entities that are in the corner in ghost zone.
 666  * The ids are returned in a newly created DataArrayInt having a single component.
 667  *
 668  * \param [in] st - The structure \b without ghost cells.
 669  * \param [in] ghostLev - The size of the ghost zone (>=0)
 670  * \return DataArrayInt * - The DataArray containing all the ids the caller is to deallocate.
 671  */
 672 DataArrayInt *MEDCouplingStructuredMesh::ComputeCornersGhost(const std::vector<int>& st, int ghostLev)
 673 {
 674   if(ghostLev<0)
 675     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeCornersGhost : ghost lev must be >= 0 !");
 676   std::size_t dim(st.size());
 677   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New());
 678   switch(dim)
 679   {
 680     case 1:
 681       {
 682         ret->alloc(2*ghostLev,1);
 683         int *ptr(ret->getPointer());
 684         for(int i=0;i<ghostLev;i++,ptr++)
 685           *ptr=i;
 686         int offset(st[0]);
 687         if(offset<0)
 688           throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeCornersGhost : element in 1D structure must be >= 0 !");
 689         for(int i=0;i<ghostLev;i++,ptr++)
 690           *ptr=offset+ghostLev+i;
 691         break;
 692       }
 693     case 2:
 694       {
 695         int offsetX(st[0]),offsetY(st[1]);
 696         if(offsetX<0 || offsetY<0)
 697           throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeCornersGhost : elements in 2D structure must be >= 0 !");
 698         ret->alloc(4*ghostLev,1);
 699         int *ptr(ret->getPointer());
 700         for(int i=0;i<ghostLev;i++)
 701           {
 702             *ptr++=i*(2*ghostLev+offsetX+1);
 703             *ptr++=offsetX+2*ghostLev-1+i*(2*ghostLev+offsetX-1);
 704           }
 705         for(int i=0;i<ghostLev;i++)
 706           {
 707             *ptr++=(2*ghostLev+offsetX)*(offsetY+ghostLev)+ghostLev-1+i*(2*ghostLev+offsetX-1);
 708             *ptr++=(2*ghostLev+offsetX)*(offsetY+ghostLev)+offsetX+ghostLev+i*(2*ghostLev+offsetX+1);
 709           }
 710         break;
 711       }
 712     case 3:
 713       {
 714         int offsetX(st[0]),offsetY(st[1]),offsetZ(st[2]);
 715         if(offsetX<0 || offsetY<0 || offsetZ<0)
 716           throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeCornersGhost : elements in 3D structure must be >= 0 !");
 717         ret->alloc(8*ghostLev,1);
 718         int *ptr(ret->getPointer());
 719         int zeOffsetZ((offsetX+2*ghostLev)*(offsetY+2*ghostLev));
 720         for(int i=0;i<ghostLev;i++)
 721           {
 722             *ptr++=i*(2*ghostLev+offsetX+1)+i*zeOffsetZ;
 723             *ptr++=offsetX+2*ghostLev-1+i*(2*ghostLev+offsetX-1)+i*zeOffsetZ;
 724             *ptr++=(2*ghostLev+offsetX)*(offsetY+ghostLev)+ghostLev-1+(ghostLev-i-1)*(2*ghostLev+offsetX-1)+i*zeOffsetZ;
 725             *ptr++=(2*ghostLev+offsetX)*(offsetY+ghostLev)+offsetX+ghostLev+(ghostLev-i-1)*(2*ghostLev+offsetX+1)+i*zeOffsetZ;
 726           }
 727         int j(0),zeOffsetZ2(zeOffsetZ*(offsetZ+ghostLev));
 728         for(int i=ghostLev-1;i>=0;i--,j++)
 729           {
 730             *ptr++=i*(2*ghostLev+offsetX+1)+j*zeOffsetZ+zeOffsetZ2;
 731             *ptr++=offsetX+2*ghostLev-1+i*(2*ghostLev+offsetX-1)+j*zeOffsetZ+zeOffsetZ2;
 732             *ptr++=(2*ghostLev+offsetX)*(offsetY+ghostLev)+ghostLev-1+(ghostLev-i-1)*(2*ghostLev+offsetX-1)+j*zeOffsetZ+zeOffsetZ2;
 733             *ptr++=(2*ghostLev+offsetX)*(offsetY+ghostLev)+offsetX+ghostLev+(ghostLev-i-1)*(2*ghostLev+offsetX+1)+j*zeOffsetZ+zeOffsetZ2;
 734           }
 735         break;
 736       }
 737     default:
 738       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeCornersGhost : Only dimensions 1, 2 and 3 are supported actually !");
 739   }
 740   return ret.retn();
 741 }
 742
 743 /*!
 744  * This method retrieves the number of entities (it can be cells or nodes) given a range in compact standard format
 745  * used in methods like BuildExplicitIdsFrom,IsPartStructured.
 746  *
 747  * \sa BuildExplicitIdsFrom,IsPartStructured
 748  */
 749 int MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(const std::vector< std::pair<int,int> >& partCompactFormat)
 750 {
 751   int ret(1);
 752   std::size_t ii(0);
 753   for(std::vector< std::pair<int,int> >::const_iterator it=partCompactFormat.begin();it!=partCompactFormat.end();it++,ii++)
 754     {
 755       int a((*it).first),b((*it).second);
 756       if(a<0 || b<0 || b-a<0)
 757         {
 758           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt : invalid input at dimension " << ii << " !";
 759           throw INTERP_KERNEL::Exception(oss.str().c_str());
 760         }
 761       ret*=(b-a);
 762     }
 763   return ret;
 764 }
 765
 766 int MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(const std::vector<int>& st)
 767 {
 768   int ret(1);
 769   bool isFetched(false);
 770   for(std::size_t i=0;i<st.size();i++)
 771     {
 772       if(st[i]<0)
 773         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure : presence of a negative value in structure !");
 774       ret*=st[i];
 775       isFetched=true;
 776     }
 777   return isFetched?ret:0;
 778 }
 779
 780 void MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt(const std::vector< std::pair<int,int> >& partCompactFormat, int& axisId, int& sizeOfRange)
 781 {
 782     int dim((int)partCompactFormat.size());
 783     int ret(-1);
 784     for(int i=0;i<dim;i++)
 785       {
 786         int curDelta(partCompactFormat[i].second-partCompactFormat[i].first);
 787         if(curDelta<0)
 788           {
 789             std::ostringstream oss; oss << "MEDCouplingStructuredMesh::FindTheWidestAxisOfGivenRangeInCompactFrmt : at axis #" << i << " the range is invalid (first value < second value) !";
 790             throw INTERP_KERNEL::Exception(oss.str().c_str());
 791           }
 792         if(curDelta>ret)
 793           {
 794             axisId=i; sizeOfRange=curDelta;
 795             ret=curDelta;
 796           }
 797       }
 798 }
 799
 800 /*!
 801  * This method is \b NOT wrapped in python because it has no sense in python (for performance reasons).
 802  * This method starts from a structured mesh with structure \a st on which a boolean field \a crit is set.
 803  * This method find for such minimalist information of mesh and field which is the part of the mesh, given by the range per axis in output parameter
 804  * \a partCompactFormat that contains all the True in \a crit. The returned vector of boolean is the field reduced to that part.
 805  * So the number of True is equal in \a st and in returned vector of boolean.
 806  *
 807  * \param [in] minPatchLgth - minimum length that the patch may have for all directions.
 808  * \param [in] st - The structure per axis of the structured mesh considered.
 809  * \param [in] crit - The field of boolean (for performance reasons) lying on the mesh defined by \a st.
 810  * \param [out] partCompactFormat - The minimal part of \a st containing all the true of \a crit.
 811  * \param [out] reducedCrit - The reduction of \a criterion on \a partCompactFormat.
 812  * \return - The number of True in \a st (that is equal to those in \a reducedCrit)
 813  */
 814 int MEDCouplingStructuredMesh::FindMinimalPartOf(int minPatchLgth, const std::vector<int>& st, const std::vector<bool>& crit, std::vector<bool>& reducedCrit, std::vector< std::pair<int,int> >& partCompactFormat)
 815 {
 816   if(minPatchLgth<0)
 817     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : the input minPatchLgth has to be >=0 !");
 818   if((int)crit.size()!=DeduceNumberOfGivenStructure(st))
 819     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : size of vector of boolean is invalid regarding the declared structure !");
 820   int ret(-1);
 821   switch((int)st.size())
 822   {
 823     case 1:
 824       {
 825         ret=FindMinimalPartOf1D(st,crit,partCompactFormat);
 826         break;
 827       }
 828     case 2:
 829       {
 830         ret=FindMinimalPartOf2D(st,crit,partCompactFormat);
 831         break;
 832       }
 833     case 3:
 834       {
 835         ret=FindMinimalPartOf3D(st,crit,partCompactFormat);
 836         break;
 837       }
 838     default:
 839       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : only dimension 1, 2 and 3 are supported actually !");
 840   }
 841   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(partCompactFormat));
 842   int i(0);
 843   for(std::vector< std::pair<int,int> >::iterator it=partCompactFormat.begin();it!=partCompactFormat.end();it++,i++)
 844     {
 845       if(st[i]<minPatchLgth)
 846         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf : the input patch is tinier than the min length constraint !");
 847       int start((*it).first),stop((*it).second),middle((start+stop)/2);
 848       if(stop-start<minPatchLgth)
 849         {
 850           (*it).first=middle-minPatchLgth/2;
 851           (*it).second=middle+minPatchLgth-minPatchLgth/2;
 852           if((*it).first<0)
 853             {
 854               (*it).second+=-(*it).first;
 855               (*it).first=0;
 856             }
 857           if((*it).second>st[i])
 858             {
 859               (*it).first-=(*it).second-st[i];
 860               (*it).second=st[i];
 861             }
 862         }
 863     }
 864   ExtractFieldOfBoolFrom(st,crit,partCompactFormat,reducedCrit);
 865   return ret;
 866 }
 867
 868 /*!
 869  * This method is \b NOT wrapped in python.
 870  * This method considers \a crit input parameter as a matrix having dimensions specified by \a st. This method returns for each axis
 871  * the signature, that is to say the number of elems equal to true in \a crit along this axis.
 872  */
 873 std::vector< std::vector<int> > MEDCouplingStructuredMesh::ComputeSignaturePerAxisOf(const std::vector<int>& st, const std::vector<bool>& crit)
 874 {
 875   int dim((int)st.size());
 876   std::vector< std::vector<int> > ret(dim);
 877   switch(dim)
 878   {
 879     case 1:
 880       {
 881         int nx(st[0]);
 882         ret[0].resize(nx);
 883         std::vector<int>& retX(ret[0]);
 884         for(int i=0;i<nx;i++)
 885           retX[i]=crit[i]?1:0;
 886         break;
 887       }
 888     case 2:
 889       {
 890         int nx(st[0]),ny(st[1]);
 891         ret[0].resize(nx); ret[1].resize(ny);
 892         std::vector<int>& retX(ret[0]);
 893         for(int i=0;i<nx;i++)
 894           {
 895             int cnt(0);
 896             for(int j=0;j<ny;j++)
 897               if(crit[j*nx+i])
 898                 cnt++;
 899             retX[i]=cnt;
 900           }
 901         std::vector<int>& retY(ret[1]);
 902         for(int j=0;j<ny;j++)
 903           {
 904             int cnt(0);
 905             for(int i=0;i<nx;i++)
 906               if(crit[j*nx+i])
 907                 cnt++;
 908             retY[j]=cnt;
 909           }
 910         break;
 911       }
 912     case 3:
 913       {
 914         int nx(st[0]),ny(st[1]),nz(st[2]);
 915         ret[0].resize(nx); ret[1].resize(ny); ret[2].resize(nz);
 916         std::vector<int>& retX(ret[0]);
 917         for(int i=0;i<nx;i++)
 918           {
 919             int cnt(0);
 920             for(int k=0;k<nz;k++)
 921               {
 922                 int offz(k*nx*ny+i);
 923                 for(int j=0;j<ny;j++)
 924                   if(crit[offz+j*nx])
 925                     cnt++;
 926               }
 927             retX[i]=cnt;
 928           }
 929         std::vector<int>& retY(ret[1]);
 930         for(int j=0;j<ny;j++)
 931           {
 932             int cnt(0),offy(j*nx);
 933             for(int k=0;k<nz;k++)
 934               {
 935                 int offz(k*nx*ny+offy);
 936                 for(int i=0;i<nx;i++)
 937                   if(crit[offz+i])
 938                     cnt++;
 939               }
 940             retY[j]=cnt;
 941           }
 942         std::vector<int>& retZ(ret[2]);
 943         for(int k=0;k<nz;k++)
 944           {
 945             int cnt(0),offz(k*nx*ny);
 946             for(int j=0;j<ny;j++)
 947               {
 948                 int offy(offz+j*nx);
 949                 for(int i=0;i<nx;i++)
 950                   if(crit[offy+i])
 951                     cnt++;
 952               }
 953             retZ[k]=cnt;
 954           }
 955         break;
 956       }
 957     default:
 958        throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeSignatureOf : only dimensions 1, 2 and 3 are supported !");
 959   }
 960   return ret;
 961 }
 962
 963 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity1D(const int *nodeStBg)
 964 {
 965   int nbOfCells(*nodeStBg-1);
 966   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
 967   conn->alloc(2*nbOfCells,1);
 968   int *cp=conn->getPointer();
 969   for(int i=0;i<nbOfCells;i++)
 970     {
 971       cp[2*i+0]=i;
 972       cp[2*i+1]=i+1;
 973     }
 974   return conn.retn();
 975 }
 976
 977 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity2D(const int *nodeStBg)
 978 {
 979   int n1=nodeStBg[0]-1;
 980   int n2=nodeStBg[1]-1;
 981   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
 982   conn->alloc(4*n1*n2,1);
 983   int *cp=conn->getPointer();
 984   int pos=0;
 985   for(int j=0;j<n2;j++)
 986     for(int i=0;i<n1;i++,pos++)
 987       {
 988         cp[4*pos+0]=i+1+j*(n1+1);
 989         cp[4*pos+1]=i+j*(n1+1);
 990         cp[4*pos+2]=i+(j+1)*(n1+1);
 991         cp[4*pos+3]=i+1+(j+1)*(n1+1);
 992       }
 993   return conn.retn();
 994 }
 995
 996 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity3D(const int *nodeStBg)
 997 {
 998   int n1=nodeStBg[0]-1;
 999   int n2=nodeStBg[1]-1;
1000   int n3=nodeStBg[2]-1;
1001   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
1002   conn->alloc(8*n1*n2*n3,1);
1003   int *cp=conn->getPointer();
1004   int pos=0;
1005   for(int k=0;k<n3;k++)
1006     for(int j=0;j<n2;j++)
1007       for(int i=0;i<n1;i++,pos++)
1008         {
1009           int tmp=(n1+1)*(n2+1);
1010           cp[8*pos+0]=i+1+j*(n1+1)+k*tmp;
1011           cp[8*pos+1]=i+j*(n1+1)+k*tmp;
1012           cp[8*pos+2]=i+(j+1)*(n1+1)+k*tmp;
1013           cp[8*pos+3]=i+1+(j+1)*(n1+1)+k*tmp;
1014           cp[8*pos+4]=i+1+j*(n1+1)+(k+1)*tmp;
1015           cp[8*pos+5]=i+j*(n1+1)+(k+1)*tmp;
1016           cp[8*pos+6]=i+(j+1)*(n1+1)+(k+1)*tmp;
1017           cp[8*pos+7]=i+1+(j+1)*(n1+1)+(k+1)*tmp;
1018         }
1019   return conn.retn();
1020 }
1021
1022 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivityOfSubLevelMesh3D(const int *nodeStBg)
1023 {
1024   std::vector<int> ngs(3);
1025   int n0(nodeStBg[0]-1),n1(nodeStBg[1]-1),n2(nodeStBg[2]-1); ngs[0]=n0; ngs[1]=n1; ngs[2]=n2;
1026   int off0(nodeStBg[0]),off1(nodeStBg[0]*nodeStBg[1]);
1027   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
1028   conn->alloc(4*GetNumberOfCellsOfSubLevelMesh(ngs,3));
1029   int *cp(conn->getPointer());
1030   //X
1031   for(int i=0;i<nodeStBg[0];i++)
1032     for(int j=0;j<n1;j++)
1033       for(int k=0;k<n2;k++,cp+=4)
1034         { cp[0]=k*off1+j*off0+i; cp[1]=(k+1)*off1+j*off0+i; cp[2]=(k+1)*off1+(j+1)*off0+i; cp[3]=k*off1+(j+1)*off0+i; }
1035   //Y
1036   for(int j=0;j<nodeStBg[1];j++)
1037     for(int i=0;i<n0;i++)
1038       for(int k=0;k<n2;k++,cp+=4)
1039         { cp[0]=k*off1+j*off0+i; cp[1]=(k+1)*off1+j*off0+i; cp[2]=(k+1)*off1+j*off0+(i+1); cp[3]=k*off1+j*off0+(i+1); }
1040   //Z
1041   for(int k=0;k<nodeStBg[2];k++)
1042     for(int i=0;i<n0;i++)
1043       for(int j=0;j<n1;j++,cp+=4)
1044         { cp[0]=k*off1+j*off0+i; cp[1]=k*off1+j*off0+(i+1); cp[2]=k*off1+(j+1)*off0+(i+1); cp[3]=k*off1+(j+1)*off0+i; }
1045   return conn.retn();
1046 }
1047
1048 /*!
1049  * \sa MEDCouplingStructuredMesh::FindMinimalPartOf
1050  */
1051 int MEDCouplingStructuredMesh::FindMinimalPartOf1D(const std::vector<int>& st, const std::vector<bool>& crit, std::vector< std::pair<int,int> >& partCompactFormat)
1052 {
1053   if(st.size()!=1)
1054     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf1D : the input size of st must be equal to 1 !");
1055   int nxMin(std::numeric_limits<int>::max()),nxMax(-std::numeric_limits<int>::max());
1056   int nx(st[0]),ret(0);
1057   for(int i=0;i<nx;i++)
1058     {
1059       if(crit[i])
1060         {
1061           nxMin=std::min(nxMin,i); nxMax=std::max(nxMax,i);
1062           ret++;
1063         }
1064     }
1065   if(ret==0)
1066     {
1067       std::size_t sz(st.size());
1068       partCompactFormat.resize(sz);
1069       for(std::size_t i=0;i<sz;i++)
1070         {
1071           partCompactFormat[i].first=st[i]/2;
1072           partCompactFormat[i].second=st[i]/2;
1073         }
1074       return ret;
1075     }
1076   partCompactFormat.resize(1);
1077   partCompactFormat[0].first=nxMin; partCompactFormat[0].second=nxMax+1;
1078   return ret;
1079 }
1080
1081 /*!
1082  * \sa MEDCouplingStructuredMesh::FindMinimalPartOf
1083  */
1084 int MEDCouplingStructuredMesh::FindMinimalPartOf2D(const std::vector<int>& st, const std::vector<bool>& crit, std::vector< std::pair<int,int> >& partCompactFormat)
1085 {
1086   if(st.size()!=2)
1087     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf2D : the input size of st must be equal to 2 !");
1088   int nxMin(std::numeric_limits<int>::max()),nxMax(-std::numeric_limits<int>::max()),nyMin(std::numeric_limits<int>::max()),nyMax(-std::numeric_limits<int>::max());
1089   int it(0),nx(st[0]),ny(st[1]);
1090   int ret(0);
1091   for(int i=0;i<ny;i++)
1092     for(int j=0;j<nx;j++,it++)
1093       {
1094         if(crit[it])
1095           {
1096             nxMin=std::min(nxMin,j); nxMax=std::max(nxMax,j);
1097             nyMin=std::min(nyMin,i); nyMax=std::max(nyMax,i);
1098             ret++;
1099           }
1100       }
1101   if(ret==0)
1102     {
1103       std::size_t sz(st.size());
1104       partCompactFormat.resize(sz);
1105       for(std::size_t i=0;i<sz;i++)
1106         {
1107           partCompactFormat[i].first=st[i]/2;
1108           partCompactFormat[i].second=st[i]/2;
1109         }
1110       return ret;
1111     }
1112   partCompactFormat.resize(2);
1113   partCompactFormat[0].first=nxMin; partCompactFormat[0].second=nxMax+1;
1114   partCompactFormat[1].first=nyMin; partCompactFormat[1].second=nyMax+1;
1115   return ret;
1116 }
1117
1118 /*!
1119  * \sa MEDCouplingStructuredMesh::FindMinimalPartOf
1120  */
1121 int MEDCouplingStructuredMesh::FindMinimalPartOf3D(const std::vector<int>& st, const std::vector<bool>& crit, std::vector< std::pair<int,int> >& partCompactFormat)
1122 {
1123   if(st.size()!=3)
1124     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindMinimalPartOf3D : the input size of st must be equal to 3 !");
1125   int nxMin(std::numeric_limits<int>::max()),nxMax(-std::numeric_limits<int>::max()),nyMin(std::numeric_limits<int>::max()),nyMax(-std::numeric_limits<int>::max()),nzMin(std::numeric_limits<int>::max()),nzMax(-std::numeric_limits<int>::max());
1126   int it(0),nx(st[0]),ny(st[1]),nz(st[2]);
1127   int ret(0);
1128   for(int i=0;i<nz;i++)
1129     for(int j=0;j<ny;j++)
1130       for(int k=0;k<nx;k++,it++)
1131         {
1132           if(crit[it])
1133             {
1134               nxMin=std::min(nxMin,k); nxMax=std::max(nxMax,k);
1135               nyMin=std::min(nyMin,j); nyMax=std::max(nyMax,j);
1136               nzMin=std::min(nzMin,i); nzMax=std::max(nzMax,i);
1137               ret++;
1138             }
1139         }
1140   if(ret==0)
1141     {
1142       std::size_t sz(st.size());
1143       partCompactFormat.resize(sz);
1144       for(std::size_t i=0;i<sz;i++)
1145         {
1146           partCompactFormat[i].first=st[i]/2;
1147           partCompactFormat[i].second=st[i]/2;
1148         }
1149       return ret;
1150     }
1151   partCompactFormat.resize(3);
1152   partCompactFormat[0].first=nxMin; partCompactFormat[0].second=nxMax+1;
1153   partCompactFormat[1].first=nyMin; partCompactFormat[1].second=nyMax+1;
1154   partCompactFormat[2].first=nzMin; partCompactFormat[2].second=nzMax+1;
1155   return ret;
1156 }
1157
1158 /*!
1159  * This method computes given the nodal structure defined by [ \a nodeStBg , \a nodeStEnd ) the zipped form.
1160  * std::distance( \a nodeStBg, \a nodeStEnd ) is equal to the space dimension. The returned value is equal to
1161  * the meshDimension (or the zipped spaceDimension).
1162  *
1163  * \param [out] zipNodeSt - The zipped node strucutre
1164  * \return int - the
1165  */
1166 int MEDCouplingStructuredMesh::ZipNodeStructure(const int *nodeStBg, const int *nodeStEnd, int zipNodeSt[3])
1167 {
1168   int spaceDim((int)std::distance(nodeStBg,nodeStEnd));
1169   if(spaceDim>3 || spaceDim<1)
1170     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ZipNodeStructure : spaceDim must in [1,2,3] !");
1171   zipNodeSt[0]=0; zipNodeSt[1]=0; zipNodeSt[2]=0;
1172   int zippedI(0);
1173   for(int i=0;i<spaceDim;i++)
1174     {
1175       int elt(nodeStBg[i]);
1176       if(elt<1)
1177         {
1178           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::ZipNodeStructure : the input nodal structure at pos#" << i << "(" << nodeStBg[i] << ") is invalid !";
1179           throw INTERP_KERNEL::Exception(oss.str().c_str());
1180         }
1181       if(elt>=2)
1182         zipNodeSt[zippedI++]=elt;
1183     }
1184   return zippedI;
1185 }
1186
1187 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivityOfSubLevelMesh2D(const int *nodeStBg)
1188 {
1189   std::vector<int> ngs(2);
1190   int n0(nodeStBg[0]-1),n1(nodeStBg[1]-1); ngs[0]=n0; ngs[1]=n1;
1191   int off0(nodeStBg[0]);
1192   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
1193   conn->alloc(2*GetNumberOfCellsOfSubLevelMesh(ngs,2));
1194   int *cp(conn->getPointer());
1195   //X
1196   for(int i=0;i<nodeStBg[0];i++)
1197     for(int j=0;j<n1;j++,cp+=2)
1198       { cp[0]=j*off0+i; cp[1]=(j+1)*off0+i; }
1199   //Y
1200   for(int j=0;j<nodeStBg[1];j++)
1201     for(int i=0;i<n0;i++,cp+=2)
1202       { cp[0]=j*off0+i; cp[1]=j*off0+(i+1); }
1203   return conn.retn();
1204 }
1205
1206 /*!
1207  * Returns a cell id by its (i,j,k) index. The cell is located between the i-th and
1208  * ( i + 1 )-th nodes along X axis etc.
1209  *  \param [in] i - a index of node coordinates array along X axis.
1210  *  \param [in] j - a index of node coordinates array along Y axis.
1211  *  \param [in] k - a index of node coordinates array along Z axis.
1212  *  \return int - a cell id in \a this mesh.
1213  */
1214 int MEDCouplingStructuredMesh::getCellIdFromPos(int i, int j, int k) const
1215 {
1216   int tmp[3]={i,j,k};
1217   int tmp2[3];
1218   int meshDim(getMeshDimension());
1219   getSplitCellValues(tmp2);
1220   std::transform(tmp,tmp+meshDim,tmp2,tmp,std::multiplies<int>());
1221   return std::accumulate(tmp,tmp+meshDim,0);
1222 }
1223
1224 /*!
1225  * Returns a node id by its (i,j,k) index.
1226  *  \param [in] i - a index of node coordinates array along X axis.
1227  *  \param [in] j - a index of node coordinates array along Y axis.
1228  *  \param [in] k - a index of node coordinates array along Z axis.
1229  *  \return int - a node id in \a this mesh.
1230  */
1231 int MEDCouplingStructuredMesh::getNodeIdFromPos(int i, int j, int k) const
1232 {
1233   int tmp[3]={i,j,k};
1234   int tmp2[3];
1235   int spaceDim(getSpaceDimension());
1236   getSplitNodeValues(tmp2);
1237   std::transform(tmp,tmp+spaceDim,tmp2,tmp,std::multiplies<int>());
1238   return std::accumulate(tmp,tmp+spaceDim,0);
1239 }
1240
1241
1242 int MEDCouplingStructuredMesh::getNumberOfCells() const
1243 {
1244   std::vector<int> ngs(getNodeGridStructure());
1245   int ret(1);
1246   bool isCatched(false);
1247   std::size_t ii(0);
1248   for(std::vector<int>::const_iterator it=ngs.begin();it!=ngs.end();it++,ii++)
1249     {
1250       int elt(*it);
1251       if(elt<=0)
1252         {
1253           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::getNumberOfCells : at pos #" << ii << " the number of nodes in nodeStructure is " << *it << " ! Must be > 0 !";
1254           throw INTERP_KERNEL::Exception(oss.str().c_str());
1255         }
1256       if(elt>1)
1257         {
1258           ret*=elt-1;
1259           isCatched=true;
1260         }
1261     }
1262   return isCatched?ret:0;
1263 }
1264
1265 int MEDCouplingStructuredMesh::getNumberOfNodes() const
1266 {
1267   std::vector<int> ngs(getNodeGridStructure());
1268   int ret(1);
1269   for(std::vector<int>::const_iterator it=ngs.begin();it!=ngs.end();it++)
1270     ret*=*it;
1271   return ret;
1272 }
1273
1274 void MEDCouplingStructuredMesh::GetPosFromId(int nodeId, int meshDim, const int *split, int *res)
1275 {
1276   int work=nodeId;
1277   for(int i=meshDim-1;i>=0;i--)
1278     {
1279       int pos=work/split[i];
1280       work=work%split[i];
1281       res[i]=pos;
1282     }
1283 }
1284
1285 std::vector<int> MEDCouplingStructuredMesh::getCellGridStructure() const
1286 {
1287   std::vector<int> ret(getNodeGridStructure());
1288   std::transform(ret.begin(),ret.end(),ret.begin(),std::bind2nd(std::plus<int>(),-1));
1289   return ret;
1290 }
1291
1292 /*!
1293  * This method returns the squareness of \a this (quadrature). \a this is expected to be with a mesh dimension equal to 2 or 3.
1294  */
1295 double MEDCouplingStructuredMesh::computeSquareness() const
1296 {
1297   std::vector<int> cgs(getCellGridStructure());
1298   if(cgs.empty())
1299     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::computeSquareness : empty mesh !");
1300   std::size_t dim(cgs.size());
1301   if(dim==1)
1302     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::computeSquareness : A segment cannot be square !");
1303   if(dim<4)
1304     {
1305       int minAx(cgs[0]),maxAx(cgs[0]);
1306       for(std::size_t i=1;i<dim;i++)
1307         {
1308           minAx=std::min(minAx,cgs[i]);
1309           maxAx=std::max(maxAx,cgs[i]);
1310         }
1311       return (double)minAx/(double)maxAx;
1312     }
1313   throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::computeSquareness : only dimension 2 and 3 supported !");
1314 }
1315
1316 /*!
1317  * Given a struct \a strct it returns a split vector [1,strct[0],strct[0]*strct[1]...]
1318  * This decomposition allows to quickly find i,j,k given a global id.
1319  */
1320 std::vector<int> MEDCouplingStructuredMesh::GetSplitVectFromStruct(const std::vector<int>& strct)
1321 {
1322   int spaceDim((int)strct.size());
1323   std::vector<int> res(spaceDim);
1324   for(int l=0;l<spaceDim;l++)
1325     {
1326       int val=1;
1327       for(int p=0;p<spaceDim-l-1;p++)
1328         val*=strct[p];
1329       res[spaceDim-l-1]=val;
1330     }
1331   return res;
1332 }
1333
1334 /*!
1335  * This method states if given part ids [ \a startIds, \a stopIds) and a structure \a st returns if it can be considered as a structured dataset.
1336  * If true is returned \a partCompactFormat will contain the information to build the corresponding part.
1337  *
1338  * \sa MEDCouplingStructuredMesh::BuildExplicitIdsFrom, MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt
1339  */
1340 bool MEDCouplingStructuredMesh::IsPartStructured(const int *startIds, const int *stopIds, const std::vector<int>& st, std::vector< std::pair<int,int> >& partCompactFormat)
1341 {
1342   int dim((int)st.size());
1343   partCompactFormat.resize(dim);
1344   if(dim<1 || dim>3)
1345     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::isPartStructured : input structure must be of dimension in [1,2,3] !");
1346   std::vector<int> tmp2(dim),tmp(dim),tmp3(dim),tmp4(dim); tmp2[0]=1;
1347   for(int i=1;i<dim;i++)
1348     tmp2[i]=tmp2[i-1]*st[i-1];
1349   std::size_t sz(std::distance(startIds,stopIds));
1350   if(sz==0)
1351     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::IsPartStructured : empty input !");
1352   GetPosFromId(*startIds,dim,&tmp2[0],&tmp[0]);
1353   partCompactFormat.resize(dim);
1354   for(int i=0;i<dim;i++)
1355     partCompactFormat[i].first=tmp[i];
1356   if(tmp[dim-1]<0 || tmp[dim-1]>=st[dim-1])
1357     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::IsPartStructured : first id in input is not in valid range !");
1358   if(sz==1)
1359     {
1360       for(int i=0;i<dim;i++)
1361         partCompactFormat[i].second=tmp[i]+1;
1362       return true;
1363     }
1364   GetPosFromId(startIds[sz-1],dim,&tmp2[0],&tmp3[0]);
1365   int szExp(1);
1366   for(int i=0;i<dim;i++)
1367     {
1368       if(tmp3[i]<0 || tmp3[i]>=st[i])
1369         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::IsPartStructured : last id in input is not in valid range !");
1370       partCompactFormat[i].second=tmp3[i]+1;
1371       tmp4[i]=partCompactFormat[i].second-partCompactFormat[i].first;
1372       if(tmp4[i]<=0)
1373         return false;
1374       szExp*=tmp4[i];
1375     }
1376   if(szExp!=(int)sz)
1377     return false;
1378   const int *w(startIds);
1379   switch(dim)
1380   {
1381     case 3:
1382       {
1383         for(int i=0;i<tmp4[2];i++)
1384           {
1385             int a=tmp2[2]*(partCompactFormat[2].first+i);
1386             for(int j=0;j<tmp4[1];j++)
1387               {
1388                 int b=tmp2[1]*(partCompactFormat[1].first+j);
1389                 for(int k=0;k<tmp4[0];k++,w++)
1390                   {
1391                     if(partCompactFormat[0].first+k+b+a!=*w)
1392                       return false;
1393                   }
1394               }
1395           }
1396         return true;
1397       }
1398     case 2:
1399       {
1400         for(int j=0;j<tmp4[1];j++)
1401           {
1402             int b=tmp2[1]*(partCompactFormat[1].first+j);
1403             for(int k=0;k<tmp4[0];k++,w++)
1404               {
1405                 if(partCompactFormat[0].first+k+b!=*w)
1406                   return false;
1407               }
1408           }
1409         return true;
1410       }
1411     case 1:
1412       {
1413         for(int k=0;k<tmp4[0];k++,w++)
1414           {
1415             if(partCompactFormat[0].first+k!=*w)
1416               return false;
1417           }
1418         return true;
1419       }
1420     default:
1421       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::IsPartStructured : internal error !");
1422   }
1423 }
1424
1425 /*!
1426  * This method takes in input a compact format [[Xmax,Xmin),[Ymin,Ymax)] and returns the corresponding dimensions for each axis that is to say
1427  * [Xmax-Xmin,Ymax-Ymin].
1428  *
1429  * \throw if an axis range is so that max<min
1430  * \sa GetCompactFrmtFromDimensions
1431  */
1432 std::vector<int> MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(const std::vector< std::pair<int,int> >& partCompactFormat)
1433 {
1434   std::vector<int> ret(partCompactFormat.size());
1435   for(std::size_t i=0;i<partCompactFormat.size();i++)
1436     {
1437       if(partCompactFormat[i].first>partCompactFormat[i].second)
1438         {
1439           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt : For axis #" << i << " end is before start !";
1440           throw INTERP_KERNEL::Exception(oss.str().c_str());
1441         }
1442       ret[i]=partCompactFormat[i].second-partCompactFormat[i].first;
1443     }
1444   return ret;
1445 }
1446
1447 /*!
1448  * This method takes in input a vector giving the number of entity per axis and returns for each axis a range starting from [0,0...]
1449  *
1450  * \throw if there is an axis in \a dims that is < 0.
1451  * \sa GetDimensionsFromCompactFrmt, ChangeReferenceFromGlobalOfCompactFrmt, ChangeReferenceToGlobalOfCompactFrmt
1452  */
1453 std::vector< std::pair<int,int> > MEDCouplingStructuredMesh::GetCompactFrmtFromDimensions(const std::vector<int>& dims)
1454 {
1455   std::size_t sz(dims.size());
1456   std::vector< std::pair<int,int> > ret(sz);
1457   for(std::size_t i=0;i<sz;i++)
1458     {
1459       if(dims[i]<0)
1460         {
1461           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt : For axis #" << i << " dimension < 0 !";
1462           throw INTERP_KERNEL::Exception(oss.str().c_str());
1463         }
1464       ret[i].first=0;
1465       ret[i].second=dims[i];
1466     }
1467   return ret;
1468 }
1469
1470 /*!
1471  * This method returns the intersection zone of two ranges (in compact format) \a r1 and \a r2.
1472  * This method will throw exception if on one axis the intersection is empty.
1473  *
1474  * \sa AreRangesIntersect
1475  */
1476 std::vector< std::pair<int,int> > MEDCouplingStructuredMesh::IntersectRanges(const std::vector< std::pair<int,int> >& r1, const std::vector< std::pair<int,int> >& r2)
1477 {
1478   std::size_t sz(r1.size());
1479   if(sz!=r2.size())
1480     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::IntersectRanges : the two ranges must have the same dimension !");
1481   std::vector< std::pair<int,int> > ret(sz);
1482   for(std::size_t i=0;i<sz;i++)
1483     {
1484       if(r1[i].first>r1[i].second)
1485         {
1486           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::IntersectRanges : On axis " << i << " of range r1, end is before start !";
1487           throw INTERP_KERNEL::Exception(oss.str().c_str());
1488         }
1489       if(r2[i].first>r2[i].second)
1490         {
1491           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::IntersectRanges : On axis " << i << " of range r2, end is before start !";
1492           throw INTERP_KERNEL::Exception(oss.str().c_str());
1493         }
1494       ret[i].first=std::max(r1[i].first,r2[i].first);
1495       ret[i].second=std::min(r1[i].second,r2[i].second);
1496       if(ret[i].first>ret[i].second)
1497         {
1498           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::IntersectRanges : On axis " << i << " the intersection of r1 and r2 is empty !";
1499           throw INTERP_KERNEL::Exception(oss.str().c_str());
1500         }
1501     }
1502   return ret;
1503 }
1504
1505 /*!
1506  * This method states if \a r1 and \a r2 do overlap of not. If yes you can call IntersectRanges to know the intersection area.
1507  *
1508  * \sa IntersectRanges
1509  */
1510 bool MEDCouplingStructuredMesh::AreRangesIntersect(const std::vector< std::pair<int,int> >& r1, const std::vector< std::pair<int,int> >& r2)
1511 {
1512   std::size_t sz(r1.size());
1513   if(sz!=r2.size())
1514     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::AreRangesIntersect : the two ranges must have the same dimension !");
1515   for(std::size_t i=0;i<sz;i++)
1516     {
1517       if(r1[i].first>r1[i].second)
1518         {
1519           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::AreRangesIntersect : On axis " << i << " of range r1, end is before start !";
1520           throw INTERP_KERNEL::Exception(oss.str().c_str());
1521         }
1522       if(r2[i].first>r2[i].second)
1523         {
1524           std::ostringstream oss; oss << "MEDCouplingStructuredMesh::AreRangesIntersect : On axis " << i << " of range r2, end is before start !";
1525           throw INTERP_KERNEL::Exception(oss.str().c_str());
1526         }
1527       if(r1[i].second<=r2[i].first)
1528         return false;
1529       if(r1[i].first>=r2[i].second)
1530         return false;
1531     }
1532   return true;
1533 }
1534
1535 /*!
1536  * This method is close to BuildExplicitIdsFrom except that instead of returning a DataArrayInt instance containing explicit ids it
1537  * enable elems in the vector of booleans (for performance reasons). As it is method for performance, this method is \b not
1538  * available in python.
1539  *
1540  * \param [in] st The entity structure.
1541  * \param [in] partCompactFormat The compact subpart to be enabled.
1542  * \param [in,out] vectToSwitchOn Vector which fetched items are enabled.
1543  *
1544  * \sa MEDCouplingStructuredMesh::BuildExplicitIdsFrom, ExtractFieldOfBoolFrom
1545  */
1546 void MEDCouplingStructuredMesh::SwitchOnIdsFrom(const std::vector<int>& st, const std::vector< std::pair<int,int> >& partCompactFormat, std::vector<bool>& vectToSwitchOn)
1547 {
1548   if(st.size()!=partCompactFormat.size())
1549     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::SwitchOnIdsFrom : input arrays must have the same size !");
1550   if((int)vectToSwitchOn.size()!=DeduceNumberOfGivenStructure(st))
1551     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::SwitchOnIdsFrom : invalid size of input vector of boolean regarding the structure !");
1552   std::vector<int> dims(GetDimensionsFromCompactFrmt(partCompactFormat));
1553   switch(st.size())
1554   {
1555     case 3:
1556       {
1557         for(int i=0;i<dims[2];i++)
1558           {
1559             int a=(partCompactFormat[2].first+i)*st[0]*st[1];
1560             for(int j=0;j<dims[1];j++)
1561               {
1562                 int b=(partCompactFormat[1].first+j)*st[0];
1563                 for(int k=0;k<dims[0];k++)
1564                   vectToSwitchOn[partCompactFormat[0].first+k+b+a]=true;
1565               }
1566           }
1567         break;
1568       }
1569     case 2:
1570       {
1571         for(int j=0;j<dims[1];j++)
1572           {
1573             int b=(partCompactFormat[1].first+j)*st[0];
1574             for(int k=0;k<dims[0];k++)
1575               vectToSwitchOn[partCompactFormat[0].first+k+b]=true;
1576           }
1577         break;
1578       }
1579     case 1:
1580       {
1581         for(int k=0;k<dims[0];k++)
1582           vectToSwitchOn[partCompactFormat[0].first+k]=true;
1583         break;
1584       }
1585     default:
1586       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::SwitchOnIdsFrom : Dimension supported are 1,2 or 3 !");
1587   }
1588 }
1589
1590 /*!
1591  * Obviously this method is \b NOT wrapped in python.
1592  * This method is close to SwitchOnIdsFrom except that here, a sub field \a fieldOut is built starting from the input field \a fieldOfBool having the structure \a st.
1593  * The extraction is defined by \a partCompactFormat.
1594  *
1595  * \param [in] st The entity structure.
1596  * \param [in] fieldOfBool field of booleans having the size equal to \c MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(st).
1597  * \param [in] partCompactFormat The compact subpart to be enabled.
1598  * \param [out] fieldOut the result of the extraction.
1599  *
1600  * \sa MEDCouplingStructuredMesh::BuildExplicitIdsFrom, SwitchOnIdsFrom, ExtractFieldOfDoubleFrom
1601  */
1602 void MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom(const std::vector<int>& st, const std::vector<bool>& fieldOfBool, const std::vector< std::pair<int,int> >& partCompactFormat, std::vector<bool>& fieldOut)
1603 {
1604   if(st.size()!=partCompactFormat.size())
1605     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom : input arrays must have the same size !");
1606   if((int)fieldOfBool.size()!=DeduceNumberOfGivenStructure(st))
1607     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom : invalid size of input field of boolean regarding the structure !");
1608   std::vector<int> dims(GetDimensionsFromCompactFrmt(partCompactFormat));
1609   int nbOfTuplesOfOutField(DeduceNumberOfGivenStructure(dims));
1610   fieldOut.resize(nbOfTuplesOfOutField);
1611   int it(0);
1612   switch(st.size())
1613   {
1614     case 3:
1615       {
1616         for(int i=0;i<dims[2];i++)
1617           {
1618             int a=(partCompactFormat[2].first+i)*st[0]*st[1];
1619             for(int j=0;j<dims[1];j++)
1620               {
1621                 int b=(partCompactFormat[1].first+j)*st[0];
1622                 for(int k=0;k<dims[0];k++)
1623                   fieldOut[it++]=fieldOfBool[partCompactFormat[0].first+k+b+a];
1624               }
1625           }
1626         break;
1627       }
1628     case 2:
1629       {
1630         for(int j=0;j<dims[1];j++)
1631           {
1632             int b=(partCompactFormat[1].first+j)*st[0];
1633             for(int k=0;k<dims[0];k++)
1634               fieldOut[it++]=fieldOfBool[partCompactFormat[0].first+k+b];
1635           }
1636         break;
1637       }
1638     case 1:
1639       {
1640         for(int k=0;k<dims[0];k++)
1641           fieldOut[it++]=fieldOfBool[partCompactFormat[0].first+k];
1642         break;
1643       }
1644     default:
1645       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfBoolFrom : Dimension supported are 1,2 or 3 !");
1646   }
1647 }
1648
1649 /*!
1650  * This method is close to SwitchOnIdsFrom except that here, a sub field \a fieldOut is built starting from the input field \a fieldOfDbl having the structure \a st.
1651  * The extraction is defined by \a partCompactFormat.
1652  *
1653  * \param [in] st The entity structure.
1654  * \param [in] fieldOfDbl field of doubles having a number of tuples equal to \c MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(st).
1655  * \param [in] partCompactFormat The compact subpart to be enabled.
1656  * \return DataArrayDouble * -the result of the extraction.
1657  *
1658  * \sa MEDCouplingStructuredMesh::BuildExplicitIdsFrom, SwitchOnIdsFrom, ExtractFieldOfBoolFrom
1659  */
1660 DataArrayDouble *MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(const std::vector<int>& st, const DataArrayDouble *fieldOfDbl, const std::vector< std::pair<int,int> >& partCompactFormat)
1661 {
1662   if(!fieldOfDbl || !fieldOfDbl->isAllocated())
1663     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom : input array of double is NULL or not allocated!");
1664   if(st.size()!=partCompactFormat.size())
1665     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom : input arrays must have the same size !");
1666   if(fieldOfDbl->getNumberOfTuples()!=DeduceNumberOfGivenStructure(st))
1667     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom : invalid size of input array of double regarding the structure !");
1668   std::vector<int> dims(GetDimensionsFromCompactFrmt(partCompactFormat));
1669   int nbOfTuplesOfOutField(DeduceNumberOfGivenStructure(dims)),nbComp(fieldOfDbl->getNumberOfComponents());
1670   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuplesOfOutField,nbComp);
1671   ret->copyStringInfoFrom(*fieldOfDbl);
1672   double *ptRet(ret->getPointer());
1673   const double *fieldOfDblPtr(fieldOfDbl->begin());
1674   switch(st.size())
1675   {
1676     case 3:
1677       {
1678         for(int i=0;i<dims[2];i++)
1679           {
1680             int a=(partCompactFormat[2].first+i)*st[0]*st[1];
1681             for(int j=0;j<dims[1];j++)
1682               {
1683                 int b=(partCompactFormat[1].first+j)*st[0];
1684                 for(int k=0;k<dims[0];k++)
1685                   ptRet=std::copy(fieldOfDblPtr+(partCompactFormat[0].first+k+b+a)*nbComp,fieldOfDblPtr+(partCompactFormat[0].first+k+b+a+1)*nbComp,ptRet);
1686               }
1687           }
1688         break;
1689       }
1690     case 2:
1691       {
1692         for(int j=0;j<dims[1];j++)
1693           {
1694             int b=(partCompactFormat[1].first+j)*st[0];
1695             for(int k=0;k<dims[0];k++)
1696               ptRet=std::copy(fieldOfDblPtr+(partCompactFormat[0].first+k+b)*nbComp,fieldOfDblPtr+(partCompactFormat[0].first+k+b+1)*nbComp,ptRet);
1697           }
1698         break;
1699       }
1700     case 1:
1701       {
1702         for(int k=0;k<dims[0];k++)
1703           ptRet=std::copy(fieldOfDblPtr+(partCompactFormat[0].first+k)*nbComp,fieldOfDblPtr+(partCompactFormat[0].first+k+1)*nbComp,ptRet);
1704         break;
1705       }
1706     default:
1707       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom : Dimension supported are 1,2 or 3 !");
1708   }
1709   return ret.retn();
1710 }
1711
1712 /*!
1713  * This method assign a part of values in \a fieldOfDbl using entirely values of \b other.
1714  *
1715  * \param [in] st - the structure of \a fieldOfDbl.
1716  * \param [in,out] fieldOfDbl - the array that will be partially filled using \a other.
1717  * \param [in] partCompactFormat - the specification of the part.
1718  * \param [in] other - the array that will be used to fill \a fieldOfDbl.
1719  */
1720 void MEDCouplingStructuredMesh::AssignPartOfFieldOfDoubleUsing(const std::vector<int>& st, DataArrayDouble *fieldOfDbl, const std::vector< std::pair<int,int> >& partCompactFormat, const DataArrayDouble *other)
1721 {//to be optimized
1722   std::vector<int> facts(st.size(),1.);
1723   MEDCouplingIMesh::CondenseFineToCoarse(st,other,partCompactFormat,facts,fieldOfDbl);
1724 }
1725
1726 /*!
1727  * This method changes the reference of a part of structured mesh \a partOfBigInAbs define in absolute reference to a new reference \a bigInAbs.
1728  * So this method only performs a translation by doing \a partOfBigRelativeToBig = \a partOfBigInAbs - \a bigInAbs
1729  * This method also checks (if \a check=true) that \a partOfBigInAbs is included in \a bigInAbs.
1730  * This method is useful to extract a part from a field lying on a big mesh.
1731  *
1732  * \sa ChangeReferenceToGlobalOfCompactFrmt, BuildExplicitIdsFrom, SwitchOnIdsFrom, ExtractFieldOfBoolFrom, ExtractFieldOfDoubleFrom
1733  */
1734 void MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt(const std::vector< std::pair<int,int> >& bigInAbs, const std::vector< std::pair<int,int> >& partOfBigInAbs, std::vector< std::pair<int,int> >& partOfBigRelativeToBig, bool check)
1735 {
1736   std::size_t dim(bigInAbs.size());
1737   if(dim!=partOfBigInAbs.size())
1738     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt : The size of parts (dimension) must be the same !");
1739   partOfBigRelativeToBig.resize(dim);
1740   for(std::size_t i=0;i<dim;i++)
1741     {
1742       if(check)
1743         {
1744           if(bigInAbs[i].first>bigInAbs[i].second)
1745             {
1746               std::ostringstream oss; oss << "MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt : Error at axis #" << i << " the input big part invalid, end before start !";
1747               throw INTERP_KERNEL::Exception(oss.str().c_str());
1748             }
1749           if(partOfBigInAbs[i].first<bigInAbs[i].first || partOfBigInAbs[i].first>=bigInAbs[i].second)
1750             {
1751               std::ostringstream oss; oss << "MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt : Error at axis #" << i << " the part is not included in the big one (start) !";
1752               throw INTERP_KERNEL::Exception(oss.str().c_str());
1753             }
1754         }
1755       partOfBigRelativeToBig[i].first=partOfBigInAbs[i].first-bigInAbs[i].first;
1756       if(check)
1757         {
1758           if(partOfBigInAbs[i].second<partOfBigInAbs[i].first || partOfBigInAbs[i].second>bigInAbs[i].second)
1759             {
1760               std::ostringstream oss; oss << "MEDCouplingStructuredMesh::ChangeReferenceFromGlobalOfCompactFrmt : Error at axis #" << i << " the part is not included in the big one (end) !";
1761               throw INTERP_KERNEL::Exception(oss.str().c_str());
1762             }
1763         }
1764       partOfBigRelativeToBig[i].second=partOfBigInAbs[i].second-bigInAbs[i].first;
1765     }
1766 }
1767
1768 /*
1769  * This method is performs the opposite reference modification than explained in ChangeReferenceFromGlobalOfCompactFrmt.
1770  *
1771  * \sa ChangeReferenceFromGlobalOfCompactFrmt
1772  */
1773 void MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt(const std::vector< std::pair<int,int> >& bigInAbs, const std::vector< std::pair<int,int> >& partOfBigRelativeToBig, std::vector< std::pair<int,int> >& partOfBigInAbs, bool check)
1774 {
1775   std::size_t dim(bigInAbs.size());
1776   if(dim!=partOfBigRelativeToBig.size())
1777     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt : The size of parts (dimension) must be the same !");
1778   partOfBigInAbs.resize(dim);
1779   for(std::size_t i=0;i<dim;i++)
1780     {
1781       if(check)
1782         {
1783           if(bigInAbs[i].first>bigInAbs[i].second)
1784             {
1785               std::ostringstream oss; oss << "MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt : Error at axis #" << i << " the input big part invalid, end before start !";
1786               throw INTERP_KERNEL::Exception(oss.str().c_str());
1787             }
1788           if(partOfBigRelativeToBig[i].first<0 || partOfBigRelativeToBig[i].first>=bigInAbs[i].second-bigInAbs[i].first)
1789             {
1790               std::ostringstream oss; oss << "MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt : Error at axis #" << i << " the start of part is not in the big one !";
1791               throw INTERP_KERNEL::Exception(oss.str().c_str());
1792             }
1793         }
1794       partOfBigInAbs[i].first=partOfBigRelativeToBig[i].first+bigInAbs[i].first;
1795       if(check)
1796         {
1797           if(partOfBigRelativeToBig[i].second<partOfBigRelativeToBig[i].first || partOfBigRelativeToBig[i].second>bigInAbs[i].second-bigInAbs[i].first)
1798             {
1799               std::ostringstream oss; oss << "MEDCouplingStructuredMesh::ChangeReferenceToGlobalOfCompactFrmt : Error at axis #" << i << " the end of part is not in the big one !";
1800               throw INTERP_KERNEL::Exception(oss.str().c_str());
1801             }
1802         }
1803       partOfBigInAbs[i].second=partOfBigRelativeToBig[i].second+bigInAbs[i].first;
1804     }
1805 }
1806
1807 /*!
1808  * This method performs a translation (defined by \a translation) of \a part and returns the result of translated part.
1809  *
1810  * \sa FindTranslationFrom
1811  */
1812 std::vector< std::pair<int,int> > MEDCouplingStructuredMesh::TranslateCompactFrmt(const std::vector< std::pair<int,int> >& part, const std::vector<int>& translation)
1813 {
1814   std::size_t sz(part.size());
1815   if(translation.size()!=sz)
1816     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::TranslateCompactFrmt : the sizes are not equal !");
1817   std::vector< std::pair<int,int> > ret(sz);
1818   for(std::size_t i=0;i<sz;i++)
1819     {
1820       ret[i].first=part[i].first+translation[i];
1821       ret[i].second=part[i].second+translation[i];
1822     }
1823   return ret;
1824 }
1825
1826 /*!
1827  * \sa TranslateCompactFrmt
1828  */
1829 std::vector<int> MEDCouplingStructuredMesh::FindTranslationFrom(const std::vector< std::pair<int,int> >& startingFrom, const std::vector< std::pair<int,int> >& goingTo)
1830 {
1831   std::size_t sz(startingFrom.size());
1832   if(goingTo.size()!=sz)
1833     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::FindTranslationFrom : the sizes are not equal !");
1834   std::vector< int > ret(sz);
1835   for(std::size_t i=0;i<sz;i++)
1836     {
1837       ret[i]=goingTo[i].first-startingFrom[i].first;
1838     }
1839   return ret;
1840 }
1841
1842 /*!
1843  * This method builds the explicit entity array from the structure in \a st and the range in \a partCompactFormat.
1844  * If the range contains invalid values regarding sructure an exception will be thrown.
1845  *
1846  * \return DataArrayInt * - a new object.
1847  * \sa MEDCouplingStructuredMesh::IsPartStructured, MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt, SwitchOnIdsFrom, ExtractFieldOfBoolFrom, ExtractFieldOfDoubleFrom, MultiplyPartOf
1848  */
1849 DataArrayInt *MEDCouplingStructuredMesh::BuildExplicitIdsFrom(const std::vector<int>& st, const std::vector< std::pair<int,int> >& partCompactFormat)
1850 {
1851   if(st.size()!=partCompactFormat.size())
1852     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::BuildExplicitIdsFrom : input arrays must have the same size !");
1853   int nbOfItems(1);
1854   std::vector<int> dims(st.size());
1855   for(std::size_t i=0;i<st.size();i++)
1856     {
1857       if(partCompactFormat[i].first<0 || partCompactFormat[i].first>st[i])
1858         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::BuildExplicitIdsFrom : invalid input range 1 !");
1859       if(partCompactFormat[i].second<0 || partCompactFormat[i].second>st[i])
1860         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::BuildExplicitIdsFrom : invalid input range 2 !");
1861       if(partCompactFormat[i].second<partCompactFormat[i].first)
1862         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::BuildExplicitIdsFrom : invalid input range 3 !");
1863       dims[i]=partCompactFormat[i].second-partCompactFormat[i].first;
1864       nbOfItems*=dims[i];
1865     }
1866   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New());
1867   ret->alloc(nbOfItems,1);
1868   int *pt(ret->getPointer());
1869   switch(st.size())
1870   {
1871     case 3:
1872       {
1873         for(int i=0;i<dims[2];i++)
1874           {
1875             int a=(partCompactFormat[2].first+i)*st[0]*st[1];
1876             for(int j=0;j<dims[1];j++)
1877               {
1878                 int b=(partCompactFormat[1].first+j)*st[0];
1879                 for(int k=0;k<dims[0];k++,pt++)
1880                   *pt=partCompactFormat[0].first+k+b+a;
1881               }
1882           }
1883         break;
1884       }
1885     case 2:
1886       {
1887         for(int j=0;j<dims[1];j++)
1888           {
1889             int b=(partCompactFormat[1].first+j)*st[0];
1890             for(int k=0;k<dims[0];k++,pt++)
1891               *pt=partCompactFormat[0].first+k+b;
1892           }
1893         break;
1894       }
1895     case 1:
1896       {
1897         for(int k=0;k<dims[0];k++,pt++)
1898           *pt=partCompactFormat[0].first+k;
1899         break;
1900       }
1901     default:
1902       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::BuildExplicitIdsFrom : Dimension supported are 1,2 or 3 !");
1903   }
1904   return ret.retn();
1905 }
1906
1907 /*!
1908  * This method multiplies by \a factor values in tuples located by \a part in \a da.
1909  *
1910  * \param [in] st - the structure of grid ( \b without considering ghost cells).
1911  * \param [in] part - the part in the structure ( \b without considering ghost cells) contained in grid whose structure is defined by \a st.
1912  * \param [in] factor - the factor, the tuples in \a da will be multiply by.
1913  * \param [in,out] da - The DataArray in wich only tuples specified by \a part will be modified.
1914  *
1915  * \sa BuildExplicitIdsFrom
1916  */
1917 void MEDCouplingStructuredMesh::MultiplyPartOf(const std::vector<int>& st, const std::vector< std::pair<int,int> >& part, double factor, DataArrayDouble *da)
1918 {
1919   if(!da || !da->isAllocated())
1920     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::MultiplyPartOf : DataArrayDouble instance must be not NULL and allocated !");
1921   if(st.size()!=part.size())
1922     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::MultiplyPartOf : input arrays must have the same size !");
1923   std::vector<int> dims(st.size());
1924   for(std::size_t i=0;i<st.size();i++)
1925     {
1926       if(part[i].first<0 || part[i].first>st[i])
1927         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::MultiplyPartOf : invalid input range 1 !");
1928       if(part[i].second<0 || part[i].second>st[i])
1929         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::MultiplyPartOf : invalid input range 2 !");
1930       if(part[i].second<part[i].first)
1931         throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::MultiplyPartOf : invalid input range 3 !");
1932       dims[i]=part[i].second-part[i].first;
1933     }
1934   int nbOfTuplesExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(st)),nbCompo(da->getNumberOfComponents());
1935   if(da->getNumberOfTuples()!=nbOfTuplesExp)
1936     {
1937       std::ostringstream oss; oss << "MEDCouplingStructuredMesh::MultiplyPartOf : invalid nb of tuples ! Expected " << nbOfTuplesExp << " having " << da->getNumberOfTuples() << " !";
1938       throw INTERP_KERNEL::Exception(oss.str().c_str());
1939     }
1940   double *pt(da->getPointer());
1941   switch(st.size())
1942   {
1943     case 3:
1944       {
1945         for(int i=0;i<dims[2];i++)
1946           {
1947             int a=(part[2].first+i)*st[0]*st[1];
1948             for(int j=0;j<dims[1];j++)
1949               {
1950                 int b=(part[1].first+j)*st[0];
1951                 for(int k=0;k<dims[0];k++)
1952                   {
1953                     int offset(part[0].first+k+b+a);
1954                     std::transform(pt+nbCompo*offset,pt+nbCompo*(offset+1),pt+nbCompo*offset,std::bind2nd(std::multiplies<double>(),factor));
1955                   }
1956               }
1957           }
1958         break;
1959       }
1960     case 2:
1961       {
1962         for(int j=0;j<dims[1];j++)
1963           {
1964             int b=(part[1].first+j)*st[0];
1965             for(int k=0;k<dims[0];k++)
1966               {
1967                 int offset(part[0].first+k+b);
1968                 std::transform(pt+nbCompo*offset,pt+nbCompo*(offset+1),pt+nbCompo*offset,std::bind2nd(std::multiplies<double>(),factor));
1969               }
1970           }
1971         break;
1972       }
1973     case 1:
1974       {
1975         for(int k=0;k<dims[0];k++)
1976           {
1977             int offset(part[0].first+k);
1978             std::transform(pt+nbCompo*offset,pt+nbCompo*(offset+1),pt+nbCompo*offset,std::bind2nd(std::multiplies<double>(),factor));
1979           }
1980         break;
1981       }
1982     default:
1983       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::MultiplyPartOf : Dimension supported are 1,2 or 3 !");
1984   }
1985 }
1986
1987 /*!
1988  * This method multiplies by \a factor values in tuples located by \a part in \a da.
1989  *
1990  * \param [in] st - the structure of grid ( \b without considering ghost cells).
1991  * \param [in] part - the part in the structure ( \b without considering ghost cells) contained in grid whose structure is defined by \a st.
1992  * \param [in] ghostSize - \a ghostSize must be >= 0.
1993  * \param [in] factor - the factor, the tuples in \a da will be multiply by.
1994  * \param [in,out] da - The DataArray in wich only tuples specified by \a part will be modified.
1995  *
1996  * \sa MultiplyPartOf, PutInGhostFormat
1997  */
1998 void MEDCouplingStructuredMesh::MultiplyPartOfByGhost(const std::vector<int>& st, const std::vector< std::pair<int,int> >& part, int ghostSize, double factor, DataArrayDouble *da)
1999 {
2000   std::vector<int> stWG;
2001   std::vector< std::pair<int,int> > partWG;
2002   PutInGhostFormat(ghostSize,st,part,stWG,partWG);
2003   MultiplyPartOf(stWG,partWG,factor,da);
2004 }
2005
2006 /*!
2007  * This method multiplies by \a factor values in tuples located by \a part in \a da.
2008  *
2009  * \param [in] st - the structure of grid ( \b without considering ghost cells).
2010  * \param [in] part - the part in the structure ( \b without considering ghost cells) contained in grid whose structure is defined by \a st.
2011  * \param [in] ghostSize - \a ghostSize must be >= 0.
2012  * \param [out] stWithGhost - the structure considering ghost cells.
2013  * \param [out] partWithGhost - the part considering the ghost cells.
2014  *
2015  * \sa MultiplyPartOf, PutInGhostFormat
2016  */
2017 void MEDCouplingStructuredMesh::PutInGhostFormat(int ghostSize, const std::vector<int>& st, const std::vector< std::pair<int,int> >& part, std::vector<int>& stWithGhost, std::vector< std::pair<int,int> >&partWithGhost)
2018 {
2019   if(ghostSize<0)
2020     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::PutInGhostFormat : ghost size must be >= 0 !");
2021   std::size_t dim(part.size());
2022   if(st.size()!=dim)
2023     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::PutInGhostFormat : the dimension of input vectors must be the same !");
2024   for(std::size_t i=0;i<dim;i++)
2025     if(part[i].first<0 || part[i].first>part[i].second || part[i].second>st[i])
2026       throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::PutInGhostFormat : the specified part is invalid ! The begin must be >= 0 and <= end ! The end must be <= to the size at considered dimension !");
2027   stWithGhost.resize(st.size());
2028   std::transform(st.begin(),st.end(),stWithGhost.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
2029   partWithGhost=part;
2030   ApplyGhostOnCompactFrmt(partWithGhost,ghostSize);
2031 }
2032
2033 /*!
2034  * \param [in,out] partBeforeFact - the part of a image mesh in compact format that will be put in ghost reference.
2035  * \param [in] ghostSize - the ghost size of zone for all axis.
2036  */
2037 void MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt(std::vector< std::pair<int,int> >& partBeforeFact, int ghostSize)
2038 {
2039   if(ghostSize<0)
2040     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ApplyGhostOnCompactFrmt : ghost size must be >= 0 !");
2041   std::size_t sz(partBeforeFact.size());
2042   for(std::size_t i=0;i<sz;i++)
2043     {
2044       partBeforeFact[i].first+=ghostSize;
2045       partBeforeFact[i].second+=ghostSize;
2046     }
2047 }
2048
2049 int MEDCouplingStructuredMesh::GetNumberOfCellsOfSubLevelMesh(const std::vector<int>& cgs, int mdim)
2050 {
2051   int ret(0);
2052   for(int i=0;i<mdim;i++)
2053     {
2054       int locRet(1);
2055       for(int j=0;j<mdim;j++)
2056         if(j!=i)
2057           locRet*=cgs[j];
2058         else
2059           locRet*=cgs[j]+1;
2060       ret+=locRet;
2061     }
2062   return ret;
2063 }