src/MEDCoupling/MEDCoupling1GTUMesh.cxx

   1 // Copyright (C) 2007-2016  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 "MEDCoupling1GTUMesh.hxx"
  22 #include "MEDCouplingUMesh.hxx"
  23 #include "MEDCouplingFieldDouble.hxx"
  24 #include "MEDCouplingCMesh.hxx"
  25
  26 #include "SplitterTetra.hxx"
  27 #include "DiameterCalculator.hxx"
  28 #include "InterpKernelAutoPtr.hxx"
  29
  30 using namespace MEDCoupling;
  31
  32 const int MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS[6]={0,1,2,4,3,5};
  33
  34 MEDCoupling1GTUMesh::MEDCoupling1GTUMesh():_cm(0)
  35 {
  36 }
  37
  38 MEDCoupling1GTUMesh::MEDCoupling1GTUMesh(const std::string& name, const INTERP_KERNEL::CellModel& cm):_cm(&cm)
  39 {
  40   setName(name);
  41 }
  42
  43 MEDCoupling1GTUMesh::MEDCoupling1GTUMesh(const MEDCoupling1GTUMesh& other, bool recDeepCpy):MEDCouplingPointSet(other,recDeepCpy),_cm(other._cm)
  44 {
  45 }
  46
  47 MEDCoupling1GTUMesh *MEDCoupling1GTUMesh::New(const std::string& name, INTERP_KERNEL::NormalizedCellType type)
  48 {
  49   if(type==INTERP_KERNEL::NORM_ERROR)
  50     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::New : NORM_ERROR is not a valid type to be used as base geometric type for a mesh !");
  51   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
  52   if(!cm.isDynamic())
  53     return MEDCoupling1SGTUMesh::New(name,type);
  54   else
  55     return MEDCoupling1DGTUMesh::New(name,type);
  56 }
  57
  58 MEDCoupling1GTUMesh *MEDCoupling1GTUMesh::New(const MEDCouplingUMesh *m)
  59 {
  60   if(!m)
  61     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::New : input mesh is null !");
  62   std::set<INTERP_KERNEL::NormalizedCellType> gts(m->getAllGeoTypes());
  63   if(gts.size()!=1)
  64     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::New : input mesh must have exactly one geometric type !");
  65   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(*gts.begin());
  66   if(!cm.isDynamic())
  67     return MEDCoupling1SGTUMesh::New(m);
  68   else
  69     return MEDCoupling1DGTUMesh::New(m);
  70 }
  71
  72 const INTERP_KERNEL::CellModel& MEDCoupling1GTUMesh::getCellModel() const
  73 {
  74   return *_cm;
  75 }
  76
  77 INTERP_KERNEL::NormalizedCellType MEDCoupling1GTUMesh::getCellModelEnum() const
  78 {
  79   return _cm->getEnum();
  80 }
  81
  82 int MEDCoupling1GTUMesh::getMeshDimension() const
  83 {
  84   return (int)_cm->getDimension();
  85 }
  86
  87 /*!
  88  * This method returns a newly allocated array containing cell ids (ascendingly sorted) whose geometric type are equal to type.
  89  * This method does not throw exception if geometric type \a type is not in \a this.
  90  * This method throws an INTERP_KERNEL::Exception if meshdimension of \b this is not equal to those of \b type.
  91  * The coordinates array is not considered here.
  92  *
  93  * \param [in] type the geometric type
  94  * \return cell ids in this having geometric type \a type.
  95  */
  96 DataArrayInt *MEDCoupling1GTUMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
  97 {
  98   MCAuto<DataArrayInt> ret=DataArrayInt::New();
  99   if(type==getCellModelEnum())
 100     ret->alloc(getNumberOfCells(),1);
 101   else
 102     ret->alloc(0,1);
 103   ret->iota();
 104   return ret.retn();
 105 }
 106
 107 /*!
 108  * Returns nb of cells having the geometric type \a type. No throw if no cells in \a this has the geometric type \a type.
 109  */
 110 int MEDCoupling1GTUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
 111 {
 112   return type==getCellModelEnum()?getNumberOfCells():0;
 113 }
 114
 115 /*!
 116  * Returns a type of a cell by its id.
 117  *  \param [in] cellId - the id of the cell of interest.
 118  *  \return INTERP_KERNEL::NormalizedCellType - enumeration item describing the cell type.
 119  *  \throw If \a cellId is invalid. Valid range is [0, \a this->getNumberOfCells() ).
 120  */
 121 INTERP_KERNEL::NormalizedCellType MEDCoupling1GTUMesh::getTypeOfCell(int cellId) const
 122 {
 123   if(cellId>=0 && cellId<getNumberOfCells())
 124     return getCellModelEnum();
 125   std::ostringstream oss; oss << "MEDCoupling1GTUMesh::getTypeOfCell : Requesting type of cell #" << cellId << " but it should be in [0," << getNumberOfCells() << ") !";
 126   throw INTERP_KERNEL::Exception(oss.str().c_str());
 127 }
 128
 129 /*!
 130  * Returns a set of all cell types available in \a this mesh.
 131  * \return std::set<INTERP_KERNEL::NormalizedCellType> - the set of cell types.
 132  * \warning this method does not throw any exception even if \a this is not defined.
 133  */
 134 std::set<INTERP_KERNEL::NormalizedCellType> MEDCoupling1GTUMesh::getAllGeoTypes() const
 135 {
 136   std::set<INTERP_KERNEL::NormalizedCellType> ret;
 137   ret.insert(getCellModelEnum());
 138   return ret;
 139 }
 140
 141 /*!
 142  * This method expects that \a this is sorted by types. If not an exception will be thrown.
 143  * This method returns in the same format as code (see MEDCouplingUMesh::checkTypeConsistencyAndContig or MEDCouplingUMesh::splitProfilePerType) how
 144  * \a this is composed in cell types.
 145  * The returned array is of size 3*n where n is the number of different types present in \a this.
 146  * For every k in [0,n] ret[3*k+2]==-1 because it has no sense here.
 147  * This parameter is kept only for compatibility with other methode listed above.
 148  */
 149 std::vector<int> MEDCoupling1GTUMesh::getDistributionOfTypes() const
 150 {
 151   std::vector<int> ret(3);
 152   ret[0]=(int)getCellModelEnum(); ret[1]=getNumberOfCells(); ret[2]=-1;
 153   return ret;
 154 }
 155
 156 /*!
 157  * 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.
 158  * 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.
 159  * This method has 1 input \a profile and 3 outputs \a code \a idsInPflPerType and \a idsPerType.
 160  *
 161  * \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.
 162  * \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,
 163  *              \a idsInPflPerType[i] stores the tuple ids in \a profile that correspond to the geometric type code[3*i+0]
 164  * \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.
 165  *              This vector can be empty in case of all geometric type cells are fully covered in ascending in the given input \a profile.
 166  *
 167  * \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.
 168  *
 169  * \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
 170  *
 171  *  \b Example1: <br>
 172  *          - Before \a this has 3 cells \a profile contains [0,1,2]
 173  *          - After \a code contains [NORM_...,nbCells,-1], \a idsInPflPerType [[0,1,2]] and \a idsPerType is empty <br>
 174  *
 175  *  \b Example2: <br>
 176  *          - Before \a this has 3 cells \a profile contains [1,2]
 177  *          - After \a code contains [NORM_...,nbCells,0], \a idsInPflPerType [[0,1]] and \a idsPerType is [[1,2]] <br>
 178
 179  */
 180 void MEDCoupling1GTUMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const
 181 {
 182   if(!profile)
 183     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::splitProfilePerType : input profile is NULL !");
 184   if(profile->getNumberOfComponents()!=1)
 185     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::splitProfilePerType : input profile should have exactly one component !");
 186   int nbTuples(profile->getNumberOfTuples()),nbOfCells(getNumberOfCells());
 187   code.resize(3); idsInPflPerType.resize(1);
 188   code[0]=(int)getCellModelEnum(); code[1]=nbTuples;
 189   idsInPflPerType.resize(1);
 190   if(profile->isIota(nbOfCells))
 191     {
 192       code[2]=-1;
 193       idsInPflPerType[0]=const_cast<DataArrayInt *>(profile); idsInPflPerType[0]->incrRef();
 194       idsPerType.clear();
 195       return ;
 196     }
 197   code[2]=0;
 198   profile->checkAllIdsInRange(0,nbOfCells);
 199   idsPerType.resize(1);
 200   idsPerType[0]=const_cast<DataArrayInt *>(profile); idsPerType[0]->incrRef();
 201   idsInPflPerType[0]=DataArrayInt::Range(0,nbTuples,1);
 202 }
 203
 204 /*!
 205  * This method tries to minimize at most the number of deep copy.
 206  * So if \a idsPerType is not empty it can be returned directly (without copy, but with ref count incremented) in return.
 207  *
 208  * \sa MEDCouplingUMesh::checkTypeConsistencyAndContig
 209  */
 210 DataArrayInt *MEDCoupling1GTUMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const
 211 {
 212   int nbOfCells=getNumberOfCells();
 213   if(code.size()!=3)
 214     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : invalid input code should be exactly of size 3 !");
 215   if(code[0]!=(int)getCellModelEnum())
 216     {
 217       std::ostringstream oss; oss << "MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : Mismatch of geometric type ! Asking for " << code[0] << " whereas the geometric type is \a this is " << getCellModelEnum() << " (" << _cm->getRepr() << ") !";
 218       throw INTERP_KERNEL::Exception(oss.str().c_str());
 219     }
 220   if(code[2]==-1)
 221     {
 222       if(code[1]==nbOfCells)
 223         return 0;
 224       else
 225         {
 226           std::ostringstream oss; oss << "MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : mismatch between the number of cells in this (" << nbOfCells << ") and the number of non profile (" << code[1] << ") !";
 227           throw INTERP_KERNEL::Exception(oss.str().c_str());
 228         }
 229     }
 230   if(code[2]!=0)
 231     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : single geo type mesh ! 0 or -1 is expected at pos #2 of input code !");
 232   if(idsPerType.size()!=1)
 233     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : input code points to DataArrayInt #0 whereas the size of idsPerType is not equal to 1 !");
 234   const DataArrayInt *pfl=idsPerType[0];
 235   if(!pfl)
 236     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : the input code points to a NULL DataArrayInt at rank 0 !");
 237   if(pfl->getNumberOfComponents()!=1)
 238     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : input profile should have exactly one component !");
 239   pfl->checkAllIdsInRange(0,nbOfCells);
 240   pfl->incrRef();
 241   return const_cast<DataArrayInt *>(pfl);
 242 }
 243
 244 void MEDCoupling1GTUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
 245 {
 246   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 247   m->writeVTKLL(ofs,cellData,pointData,byteData);
 248 }
 249
 250 std::string MEDCoupling1GTUMesh::getVTKDataSetType() const
 251 {
 252   return std::string("UnstructuredGrid");
 253 }
 254
 255 std::string MEDCoupling1GTUMesh::getVTKFileExtension() const
 256 {
 257   return std::string("vtu");
 258 }
 259
 260 std::size_t MEDCoupling1GTUMesh::getHeapMemorySizeWithoutChildren() const
 261 {
 262   return MEDCouplingPointSet::getHeapMemorySizeWithoutChildren();
 263 }
 264
 265 bool MEDCoupling1GTUMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
 266 {
 267   if(!MEDCouplingPointSet::isEqualIfNotWhy(other,prec,reason))
 268     return false;
 269   if(!other)
 270     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::isEqualIfNotWhy : input other pointer is null !");
 271   const MEDCoupling1GTUMesh *otherC=dynamic_cast<const MEDCoupling1GTUMesh *>(other);
 272   if(!otherC)
 273     {
 274       reason="mesh given in input is not castable in MEDCouplingSGTUMesh !";
 275       return false;
 276     }
 277   if(_cm!=otherC->_cm)
 278     {
 279       reason="mismatch in geometric type !";
 280       return false;
 281     }
 282   return true;
 283 }
 284
 285 bool MEDCoupling1GTUMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 286 {
 287   if(!MEDCouplingPointSet::isEqualWithoutConsideringStr(other,prec))
 288     return false;
 289   if(!other)
 290     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::isEqualWithoutConsideringStr : input other pointer is null !");
 291   const MEDCoupling1GTUMesh *otherC=dynamic_cast<const MEDCoupling1GTUMesh *>(other);
 292   if(!otherC)
 293     return false;
 294   if(_cm!=otherC->_cm)
 295     return false;
 296   return true;
 297 }
 298
 299 void MEDCoupling1GTUMesh::checkConsistencyLight() const
 300 {
 301   MEDCouplingPointSet::checkConsistencyLight();
 302 }
 303
 304 DataArrayDouble *MEDCoupling1GTUMesh::computeCellCenterOfMass() const
 305 {
 306   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 307   MCAuto<DataArrayDouble> ret=m->computeCellCenterOfMass();
 308   return ret.retn();
 309 }
 310
 311 MEDCouplingFieldDouble *MEDCoupling1GTUMesh::getMeasureField(bool isAbs) const
 312 {
 313   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 314   MCAuto<MEDCouplingFieldDouble> ret=m->getMeasureField(isAbs);
 315   ret->setMesh(this);
 316   return ret.retn();
 317 }
 318
 319 MEDCouplingFieldDouble *MEDCoupling1GTUMesh::getMeasureFieldOnNode(bool isAbs) const
 320 {
 321   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 322   MCAuto<MEDCouplingFieldDouble> ret=m->getMeasureFieldOnNode(isAbs);
 323   ret->setMesh(this);
 324   return ret.retn();
 325 }
 326
 327 /*!
 328  * to improve perf !
 329  */
 330 int MEDCoupling1GTUMesh::getCellContainingPoint(const double *pos, double eps) const
 331 {
 332   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 333   return m->getCellContainingPoint(pos,eps);
 334 }
 335
 336 MEDCouplingFieldDouble *MEDCoupling1GTUMesh::buildOrthogonalField() const
 337 {
 338   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 339   MCAuto<MEDCouplingFieldDouble> ret=m->buildOrthogonalField();
 340   ret->setMesh(this);
 341   return ret.retn();
 342 }
 343
 344 DataArrayInt *MEDCoupling1GTUMesh::getCellsInBoundingBox(const double *bbox, double eps) const
 345 {
 346   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 347   return m->getCellsInBoundingBox(bbox,eps);
 348 }
 349
 350 DataArrayInt *MEDCoupling1GTUMesh::getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps)
 351 {
 352   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 353   return m->getCellsInBoundingBox(bbox,eps);
 354 }
 355
 356 MEDCouplingPointSet *MEDCoupling1GTUMesh::buildFacePartOfMySelfNode(const int *start, const int *end, bool fullyIn) const
 357 {
 358   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 359   return m->buildFacePartOfMySelfNode(start,end,fullyIn);
 360 }
 361
 362 DataArrayInt *MEDCoupling1GTUMesh::findBoundaryNodes() const
 363 {
 364   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 365   return m->findBoundaryNodes();
 366 }
 367
 368 MEDCouplingPointSet *MEDCoupling1GTUMesh::buildBoundaryMesh(bool keepCoords) const
 369 {
 370   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 371   return m->buildBoundaryMesh(keepCoords);
 372 }
 373
 374 void MEDCoupling1GTUMesh::findCommonCells(int compType, int startCellId, DataArrayInt *& commonCellsArr, DataArrayInt *& commonCellsIArr) const
 375 {
 376   MCAuto<MEDCouplingUMesh> m=buildUnstructured();
 377   m->findCommonCells(compType,startCellId,commonCellsArr,commonCellsIArr);
 378 }
 379
 380 int MEDCoupling1GTUMesh::getNodalConnectivityLength() const
 381 {
 382   const DataArrayInt *c1(getNodalConnectivity());
 383   if(!c1)
 384     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::getNodalConnectivityLength : no connectivity set !");
 385   if(c1->getNumberOfComponents()!=1)
 386     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::getNodalConnectivityLength : Nodal connectivity array set must have exactly one component !");
 387   if(!c1->isAllocated())
 388     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::getNodalConnectivityLength : Nodal connectivity array must be allocated !");
 389   return c1->getNumberOfTuples();
 390 }
 391
 392 /*!
 393  * This method aggregates all the meshes in \a parts to put them in a single unstructured mesh (those returned).
 394  * The order of cells is the returned instance is those in the order of instances in \a parts.
 395  *
 396  * \param [in] parts - all not null parts of single geo type meshes to be aggreagated having the same mesh dimension and same coordinates.
 397  * \return MEDCouplingUMesh * - new object to be dealt by the caller.
 398  *
 399  * \throw If one element is null in \a parts.
 400  * \throw If not all the parts do not have the same mesh dimension.
 401  * \throw If not all the parts do not share the same coordinates.
 402  * \throw If not all the parts have their connectivity set properly.
 403  * \throw If \a parts is empty.
 404  */
 405 MEDCouplingUMesh *MEDCoupling1GTUMesh::AggregateOnSameCoordsToUMesh(const std::vector< const MEDCoupling1GTUMesh *>& parts)
 406 {
 407   if(parts.empty())
 408     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::AggregateOnSameCoordsToUMesh : input parts vector is empty !");
 409   const MEDCoupling1GTUMesh *firstPart(parts[0]);
 410   if(!firstPart)
 411     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::AggregateOnSameCoordsToUMesh : the first instance in input parts is null !");
 412   const DataArrayDouble *coords(firstPart->getCoords());
 413   int meshDim(firstPart->getMeshDimension());
 414   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(firstPart->getName(),meshDim)); ret->setDescription(firstPart->getDescription());
 415   ret->setCoords(coords);
 416   int nbOfCells(0),connSize(0);
 417   for(std::vector< const MEDCoupling1GTUMesh *>::const_iterator it=parts.begin();it!=parts.end();it++)
 418     {
 419       if(!(*it))
 420         throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::AggregateOnSameCoordsToUMesh : presence of null pointer in input vector !");
 421       if((*it)->getMeshDimension()!=meshDim)
 422         throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::AggregateOnSameCoordsToUMesh : all the instances in input vector must have same mesh dimension !");
 423       if((*it)->getCoords()!=coords)
 424         throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::AggregateOnSameCoordsToUMesh : all the instances must share the same coordinates pointer !");
 425       nbOfCells+=(*it)->getNumberOfCells();
 426       connSize+=(*it)->getNodalConnectivityLength();
 427     }
 428   MCAuto<DataArrayInt> conn(DataArrayInt::New()),connI(DataArrayInt::New());
 429   connI->alloc(nbOfCells+1,1); conn->alloc(connSize+nbOfCells,1);
 430   int *c(conn->getPointer()),*ci(connI->getPointer()); *ci=0;
 431   for(std::vector< const MEDCoupling1GTUMesh *>::const_iterator it=parts.begin();it!=parts.end();it++)
 432     {
 433       int curNbCells((*it)->getNumberOfCells());
 434       int geoType((int)(*it)->getCellModelEnum());
 435       const int *cinPtr((*it)->getNodalConnectivity()->begin());
 436       const MEDCoupling1SGTUMesh *ps(dynamic_cast<const MEDCoupling1SGTUMesh *>(*it));
 437       const MEDCoupling1DGTUMesh *pd(dynamic_cast<const MEDCoupling1DGTUMesh *>(*it));
 438       if(ps && !pd)
 439         {
 440           int nNodesPerCell(ps->getNumberOfNodesPerCell());
 441           for(int i=0;i<curNbCells;i++,ci++,cinPtr+=nNodesPerCell)
 442             {
 443               *c++=geoType;
 444               c=std::copy(cinPtr,cinPtr+nNodesPerCell,c);
 445               ci[1]=ci[0]+nNodesPerCell+1;
 446             }
 447         }
 448       else if(!ps && pd)
 449         {
 450           const int *ciinPtr(pd->getNodalConnectivityIndex()->begin());
 451           for(int i=0;i<curNbCells;i++,ci++,ciinPtr++)
 452             {
 453               *c++=geoType;
 454               c=std::copy(cinPtr+ciinPtr[0],cinPtr+ciinPtr[1],c);
 455               ci[1]=ci[0]+ciinPtr[1]-ciinPtr[0]+1;
 456             }
 457         }
 458       else
 459         throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::AggregateOnSameCoordsToUMesh : presence of instance which type is not in [MEDCoupling1SGTUMesh,MEDCoupling1DGTUMesh] !");
 460     }
 461   ret->setConnectivity(conn,connI,true);
 462   return ret.retn();
 463 }
 464
 465 //==
 466
 467 MEDCoupling1SGTUMesh::MEDCoupling1SGTUMesh(const MEDCoupling1SGTUMesh& other, bool recDeepCpy):MEDCoupling1GTUMesh(other,recDeepCpy),_conn(other._conn)
 468 {
 469   if(recDeepCpy)
 470     {
 471       const DataArrayInt *c(other._conn);
 472       if(c)
 473         _conn=c->deepCopy();
 474     }
 475 }
 476
 477 MEDCoupling1SGTUMesh::MEDCoupling1SGTUMesh(const std::string& name, const INTERP_KERNEL::CellModel& cm):MEDCoupling1GTUMesh(name,cm)
 478 {
 479 }
 480
 481 MEDCoupling1SGTUMesh::MEDCoupling1SGTUMesh()
 482 {
 483 }
 484
 485 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::New()
 486 {
 487   return new MEDCoupling1SGTUMesh;
 488 }
 489
 490 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::New(const std::string& name, INTERP_KERNEL::NormalizedCellType type)
 491 {
 492   if(type==INTERP_KERNEL::NORM_ERROR)
 493     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::New : NORM_ERROR is not a valid type to be used as base geometric type for a mesh !");
 494   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 495   if(cm.isDynamic())
 496     {
 497       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::New : the input geometric type " << cm.getRepr() << " is dynamic ! Only static types are allowed here !";
 498       throw INTERP_KERNEL::Exception(oss.str().c_str());
 499     }
 500   return new MEDCoupling1SGTUMesh(name,cm);
 501 }
 502
 503 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::New(const MEDCouplingUMesh *m)
 504 {
 505   if(!m)
 506     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::New : input mesh is null !");
 507   std::set<INTERP_KERNEL::NormalizedCellType> gts(m->getAllGeoTypes());
 508   if(gts.size()!=1)
 509     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::New : input mesh must have exactly one geometric type !");
 510   int geoType((int)*gts.begin());
 511   MCAuto<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(m->getName(),*gts.begin()));
 512   ret->setCoords(m->getCoords()); ret->setDescription(m->getDescription());
 513   int nbCells(m->getNumberOfCells());
 514   int nbOfNodesPerCell(ret->getNumberOfNodesPerCell());
 515   MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc(nbCells*nbOfNodesPerCell,1);
 516   int *c(conn->getPointer());
 517   const int *cin(m->getNodalConnectivity()->begin()),*ciin(m->getNodalConnectivityIndex()->begin());
 518   for(int i=0;i<nbCells;i++,ciin++)
 519     {
 520       if(cin[ciin[0]]==geoType)
 521         {
 522           if(ciin[1]-ciin[0]==nbOfNodesPerCell+1)
 523             c=std::copy(cin+ciin[0]+1,cin+ciin[1],c);
 524           else
 525             {
 526               std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::New(const MEDCouplingUMesh *m) : something is wrong in the input mesh at cell #" << i << " ! The size of cell is not those expected (" << nbOfNodesPerCell << ") !";
 527               throw INTERP_KERNEL::Exception(oss.str().c_str());
 528             }
 529         }
 530       else
 531         {
 532           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::New(const MEDCouplingUMesh *m) : something is wrong in the input mesh at cell #" << i << " ! The geometric type is not those expected !";
 533           throw INTERP_KERNEL::Exception(oss.str().c_str());
 534         }
 535     }
 536   ret->setNodalConnectivity(conn);
 537   try
 538   { ret->copyTinyInfoFrom(m); }
 539   catch(INTERP_KERNEL::Exception&) { }
 540   return ret.retn();
 541 }
 542
 543 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::clone(bool recDeepCpy) const
 544 {
 545   return new MEDCoupling1SGTUMesh(*this,recDeepCpy);
 546 }
 547
 548 /*!
 549  * This method behaves mostly like MEDCoupling1SGTUMesh::deepCopy method, except that only nodal connectivity arrays are deeply copied.
 550  * The coordinates are shared between \a this and the returned instance.
 551  *
 552  * \return MEDCoupling1SGTUMesh * - A new object instance holding the copy of \a this (deep for connectivity, shallow for coordiantes)
 553  * \sa MEDCoupling1SGTUMesh::deepCopy
 554  */
 555 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::deepCopyConnectivityOnly() const
 556 {
 557   checkConsistencyLight();
 558   MCAuto<MEDCoupling1SGTUMesh> ret(clone(false));
 559   MCAuto<DataArrayInt> c(_conn->deepCopy());
 560   ret->setNodalConnectivity(c);
 561   return ret.retn();
 562 }
 563
 564 void MEDCoupling1SGTUMesh::shallowCopyConnectivityFrom(const MEDCouplingPointSet *other)
 565 {
 566   if(!other)
 567     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::shallowCopyConnectivityFrom : input pointer is null !");
 568   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
 569   if(!otherC)
 570     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::shallowCopyConnectivityFrom : input pointer is not an MEDCoupling1SGTUMesh instance !");
 571   setNodalConnectivity(otherC->getNodalConnectivity());
 572 }
 573
 574 void MEDCoupling1SGTUMesh::updateTime() const
 575 {
 576   MEDCoupling1GTUMesh::updateTime();
 577   const DataArrayInt *c(_conn);
 578   if(c)
 579     updateTimeWith(*c);
 580 }
 581
 582 std::size_t MEDCoupling1SGTUMesh::getHeapMemorySizeWithoutChildren() const
 583 {
 584   return MEDCoupling1GTUMesh::getHeapMemorySizeWithoutChildren();
 585 }
 586
 587 std::vector<const BigMemoryObject *> MEDCoupling1SGTUMesh::getDirectChildrenWithNull() const
 588 {
 589   std::vector<const BigMemoryObject *> ret(MEDCoupling1GTUMesh::getDirectChildrenWithNull());
 590   ret.push_back((const DataArrayInt *)_conn);
 591   return ret;
 592 }
 593
 594 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::deepCopy() const
 595 {
 596   return clone(true);
 597 }
 598
 599 bool MEDCoupling1SGTUMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
 600 {
 601   if(!other)
 602     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::isEqualIfNotWhy : input other pointer is null !");
 603   std::ostringstream oss; oss.precision(15);
 604   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
 605   if(!otherC)
 606     {
 607       reason="mesh given in input is not castable in MEDCoupling1SGTUMesh !";
 608       return false;
 609     }
 610   if(!MEDCoupling1GTUMesh::isEqualIfNotWhy(other,prec,reason))
 611     return false;
 612   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
 613   if(c1==c2)
 614     return true;
 615   if(!c1 || !c2)
 616     {
 617       reason="in connectivity of single static geometric type exactly one among this and other is null !";
 618       return false;
 619     }
 620   if(!c1->isEqualIfNotWhy(*c2,reason))
 621     {
 622       reason.insert(0,"Nodal connectivity DataArrayInt differ : ");
 623       return false;
 624     }
 625   return true;
 626 }
 627
 628 bool MEDCoupling1SGTUMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 629 {
 630   if(!other)
 631     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::isEqualWithoutConsideringStr : input other pointer is null !");
 632   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
 633   if(!otherC)
 634     return false;
 635   if(!MEDCoupling1GTUMesh::isEqualWithoutConsideringStr(other,prec))
 636     return false;
 637   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
 638   if(c1==c2)
 639     return true;
 640   if(!c1 || !c2)
 641     return false;
 642   if(!c1->isEqualWithoutConsideringStr(*c2))
 643     return false;
 644   return true;
 645 }
 646
 647 void MEDCoupling1SGTUMesh::checkConsistencyOfConnectivity() const
 648 {
 649   const DataArrayInt *c1(_conn);
 650   if(c1)
 651     {
 652       if(c1->getNumberOfComponents()!=1)
 653         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to be with number of components set to one !");
 654       if(c1->getInfoOnComponent(0)!="")
 655         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to have no info on its single component !");
 656       c1->checkAllocated();
 657     }
 658   else
 659     throw INTERP_KERNEL::Exception("Nodal connectivity array not defined !");
 660 }
 661
 662 void MEDCoupling1SGTUMesh::checkConsistencyLight() const
 663 {
 664   MEDCouplingPointSet::checkConsistencyLight();
 665   checkConsistencyOfConnectivity();
 666 }
 667
 668 void MEDCoupling1SGTUMesh::checkConsistency(double eps) const
 669 {
 670   checkConsistencyLight();
 671   const DataArrayInt *c1(_conn);
 672   int nbOfTuples=c1->getNumberOfTuples();
 673   int nbOfNodesPerCell=(int)_cm->getNumberOfNodes();
 674   if(nbOfTuples%nbOfNodesPerCell!=0)
 675     {
 676       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::checkConsistency : the nb of tuples in conn is " << nbOfTuples << " and number of nodes per cell is " << nbOfNodesPerCell << ". But " << nbOfTuples << "%" << nbOfNodesPerCell << " !=0 !";
 677       throw INTERP_KERNEL::Exception(oss.str().c_str());
 678     }
 679   int nbOfNodes=getNumberOfNodes();
 680   int nbOfCells=nbOfTuples/nbOfNodesPerCell;
 681   const int *w(c1->begin());
 682   for(int i=0;i<nbOfCells;i++)
 683     for(int j=0;j<nbOfNodesPerCell;j++,w++)
 684       {
 685         if(*w<0 || *w>=nbOfNodes)
 686           {
 687             std::ostringstream oss; oss << "At node #" << j << " of  cell #" << i << ", is equal to " << *w << " must be in [0," << nbOfNodes << ") !";
 688             throw INTERP_KERNEL::Exception(oss.str().c_str());
 689           }
 690       }
 691 }
 692
 693 int MEDCoupling1SGTUMesh::getNumberOfCells() const
 694 {
 695   int nbOfTuples=getNodalConnectivityLength();
 696   int nbOfNodesPerCell=getNumberOfNodesPerCell();
 697   if(nbOfTuples%nbOfNodesPerCell!=0)
 698     {
 699       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh:getNumberOfCells: : the nb of tuples in conn is " << nbOfTuples << " and number of nodes per cell is " << nbOfNodesPerCell << ". But " << nbOfTuples << "%" << nbOfNodesPerCell << " !=0 !";
 700       throw INTERP_KERNEL::Exception(oss.str().c_str());
 701     }
 702   return nbOfTuples/nbOfNodesPerCell;
 703 }
 704
 705 int MEDCoupling1SGTUMesh::getNumberOfNodesInCell(int cellId) const
 706 {
 707   return getNumberOfNodesPerCell();
 708 }
 709
 710 int MEDCoupling1SGTUMesh::getNumberOfNodesPerCell() const
 711 {
 712   checkNonDynamicGeoType();
 713   return (int)_cm->getNumberOfNodes();
 714 }
 715
 716 DataArrayInt *MEDCoupling1SGTUMesh::computeNbOfNodesPerCell() const
 717 {
 718   checkNonDynamicGeoType();
 719   MCAuto<DataArrayInt> ret=DataArrayInt::New();
 720   ret->alloc(getNumberOfCells(),1);
 721   ret->fillWithValue((int)_cm->getNumberOfNodes());
 722   return ret.retn();
 723 }
 724
 725 DataArrayInt *MEDCoupling1SGTUMesh::computeNbOfFacesPerCell() const
 726 {
 727   checkNonDynamicGeoType();
 728   MCAuto<DataArrayInt> ret=DataArrayInt::New();
 729   ret->alloc(getNumberOfCells(),1);
 730   ret->fillWithValue((int)_cm->getNumberOfSons());
 731   return ret.retn();
 732 }
 733
 734 DataArrayInt *MEDCoupling1SGTUMesh::computeEffectiveNbOfNodesPerCell() const
 735 {
 736   checkNonDynamicGeoType();
 737   MCAuto<DataArrayInt> ret=DataArrayInt::New();
 738   int nbCells(getNumberOfCells());
 739   ret->alloc(nbCells,1);
 740   int *retPtr(ret->getPointer());
 741   int nbNodesPerCell(getNumberOfNodesPerCell());
 742   const int *conn(_conn->begin());
 743   for(int i=0;i<nbCells;i++,conn+=nbNodesPerCell,retPtr++)
 744     {
 745       std::set<int> s(conn,conn+nbNodesPerCell);
 746       *retPtr=(int)s.size();
 747     }
 748   return ret.retn();
 749 }
 750
 751 void MEDCoupling1SGTUMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
 752 {
 753   int sz=getNumberOfNodesPerCell();
 754   conn.resize(sz);
 755   if(cellId>=0 && cellId<getNumberOfCells())
 756     std::copy(_conn->begin()+cellId*sz,_conn->begin()+(cellId+1)*sz,conn.begin());
 757   else
 758     {
 759       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getNodeIdsOfCell : request for cellId #" << cellId << " must be in [0," << getNumberOfCells() << ") !";
 760       throw INTERP_KERNEL::Exception(oss.str().c_str());
 761     }
 762 }
 763
 764 void MEDCoupling1SGTUMesh::checkNonDynamicGeoType() const
 765 {
 766   if(_cm->isDynamic())
 767     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkNonDynamicGeoType : internal error ! the internal geo type is dynamic ! should be static !");
 768 }
 769
 770 std::string MEDCoupling1SGTUMesh::simpleRepr() const
 771 {
 772   static const char msg0[]="No coordinates specified !";
 773   std::ostringstream ret;
 774   ret << "Single static geometic type (" << _cm->getRepr() << ") unstructured mesh with name : \"" << getName() << "\"\n";
 775   ret << "Description of mesh : \"" << getDescription() << "\"\n";
 776   int tmpp1,tmpp2;
 777   double tt=getTime(tmpp1,tmpp2);
 778   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
 779   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
 780   ret << "Mesh dimension : " << getMeshDimension() << "\nSpace dimension : ";
 781   if(_coords!=0)
 782     {
 783       const int spaceDim=getSpaceDimension();
 784       ret << spaceDim << "\nInfo attached on space dimension : ";
 785       for(int i=0;i<spaceDim;i++)
 786         ret << "\"" << _coords->getInfoOnComponent(i) << "\" ";
 787       ret << "\n";
 788     }
 789   else
 790     ret << msg0 << "\n";
 791   ret << "Number of nodes : ";
 792   if(_coords!=0)
 793     ret << getNumberOfNodes() << "\n";
 794   else
 795     ret << msg0 << "\n";
 796   ret << "Number of cells : ";
 797   if((const DataArrayInt *)_conn)
 798     {
 799       if(_conn->isAllocated())
 800         {
 801           if(_conn->getNumberOfComponents()==1)
 802             ret << getNumberOfCells() << "\n";
 803           else
 804             ret << "Nodal connectivity array specified and allocated but with not exactly one component !" << "\n";
 805         }
 806       else
 807         ret << "Nodal connectivity array specified but not allocated !" << "\n";
 808     }
 809   else
 810     ret << "No connectivity specified !" << "\n";
 811   ret << "Cell type : " << _cm->getRepr() << "\n";
 812   return ret.str();
 813 }
 814
 815 std::string MEDCoupling1SGTUMesh::advancedRepr() const
 816 {
 817   std::ostringstream ret;
 818   ret << simpleRepr();
 819   ret << "\nCoordinates array : \n___________________\n\n";
 820   if(_coords)
 821     _coords->reprWithoutNameStream(ret);
 822   else
 823     ret << "No array set !\n";
 824   ret << "\n\nConnectivity array : \n____________________\n\n";
 825   //
 826   if((const DataArrayInt *)_conn)
 827     {
 828       if(_conn->isAllocated())
 829         {
 830           if(_conn->getNumberOfComponents()==1)
 831             {
 832               int nbOfCells=getNumberOfCells();
 833               int sz=getNumberOfNodesPerCell();
 834               const int *connPtr=_conn->begin();
 835               for(int i=0;i<nbOfCells;i++,connPtr+=sz)
 836                 {
 837                   ret << "Cell #" << i << " : ";
 838                   std::copy(connPtr,connPtr+sz,std::ostream_iterator<int>(ret," "));
 839                   ret << "\n";
 840                 }
 841             }
 842           else
 843             ret << "Nodal connectivity array specified and allocated but with not exactly one component !" << "\n";
 844         }
 845       else
 846         ret << "Nodal connectivity array specified but not allocated !" << "\n";
 847     }
 848   else
 849     ret << "No connectivity specified !" << "\n";
 850   return ret.str();
 851 }
 852
 853 DataArrayDouble *MEDCoupling1SGTUMesh::computeIsoBarycenterOfNodesPerCell() const
 854 {
 855   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
 856   int spaceDim=getSpaceDimension();
 857   int nbOfCells=getNumberOfCells();//checkConsistencyLight()
 858   int nbOfNodes=getNumberOfNodes();
 859   ret->alloc(nbOfCells,spaceDim);
 860   double *ptToFill=ret->getPointer();
 861   const double *coor=_coords->begin();
 862   const int *nodal=_conn->begin();
 863   int sz=getNumberOfNodesPerCell();
 864   double coeff=1./(double)sz;
 865   for(int i=0;i<nbOfCells;i++,ptToFill+=spaceDim)
 866     {
 867       std::fill(ptToFill,ptToFill+spaceDim,0.);
 868       for(int j=0;j<sz;j++,nodal++)
 869         if(*nodal>=0 && *nodal<nbOfNodes)
 870           std::transform(coor+spaceDim*nodal[0],coor+spaceDim*(nodal[0]+1),ptToFill,ptToFill,std::plus<double>());
 871         else
 872           {
 873             std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of nodeId #" << *nodal << " should be in [0," <<   nbOfNodes << ") !";
 874             throw INTERP_KERNEL::Exception(oss.str().c_str());
 875           }
 876       std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),coeff));
 877     }
 878   return ret.retn();
 879 }
 880
 881 void MEDCoupling1SGTUMesh::renumberCells(const int *old2NewBg, bool check)
 882 {
 883   int nbCells=getNumberOfCells();
 884   MCAuto<DataArrayInt> o2n=DataArrayInt::New();
 885   o2n->useArray(old2NewBg,false,C_DEALLOC,nbCells,1);
 886   if(check)
 887     o2n=o2n->checkAndPreparePermutation();
 888   //
 889   const int *conn=_conn->begin();
 890   MCAuto<DataArrayInt> n2o=o2n->invertArrayO2N2N2O(nbCells);
 891   const int *n2oPtr=n2o->begin();
 892   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
 893   newConn->alloc(_conn->getNumberOfTuples(),1);
 894   newConn->copyStringInfoFrom(*_conn);
 895   int sz=getNumberOfNodesPerCell();
 896   //
 897   int *newC=newConn->getPointer();
 898   for(int i=0;i<nbCells;i++,newC+=sz)
 899     {
 900       int pos=n2oPtr[i];
 901       std::copy(conn+pos*sz,conn+(pos+1)*sz,newC);
 902     }
 903   _conn=newConn;
 904 }
 905
 906 /*!
 907  * Keeps from \a this only cells which constituing point id are in the ids specified by [\a begin,\a end).
 908  * The resulting cell ids are stored at the end of the 'cellIdsKept' parameter.
 909  * Parameter \a fullyIn specifies if a cell that has part of its nodes in ids array is kept or not.
 910  * If \a fullyIn is true only cells whose ids are \b fully contained in [\a begin,\a end) tab will be kept.
 911  *
 912  * \param [in] begin input start of array of node ids.
 913  * \param [in] end input end of array of node ids.
 914  * \param [in] fullyIn input that specifies if all node ids must be in [\a begin,\a end) array to consider cell to be in.
 915  * \param [in,out] cellIdsKeptArr array where all candidate cell ids are put at the end.
 916  */
 917 void MEDCoupling1SGTUMesh::fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const
 918 {
 919   int nbOfCells=getNumberOfCells();
 920   MCAuto<DataArrayInt> cellIdsKept=DataArrayInt::New(); cellIdsKept->alloc(0,1);
 921   int tmp=-1;
 922   int sz=_conn->getMaxValue(tmp); sz=std::max(sz,0)+1;
 923   std::vector<bool> fastFinder(sz,false);
 924   for(const int *work=begin;work!=end;work++)
 925     if(*work>=0 && *work<sz)
 926       fastFinder[*work]=true;
 927   const int *conn=_conn->begin();
 928   int nbNodesPerCell=getNumberOfNodesPerCell();
 929   for(int i=0;i<nbOfCells;i++,conn+=nbNodesPerCell)
 930     {
 931       int ref=0,nbOfHit=0;
 932       for(int j=0;j<nbNodesPerCell;j++)
 933         if(conn[j]>=0)
 934           {
 935             ref++;
 936             if(fastFinder[conn[j]])
 937               nbOfHit++;
 938           }
 939       if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
 940         cellIdsKept->pushBackSilent(i);
 941     }
 942   cellIdsKeptArr=cellIdsKept.retn();
 943 }
 944
 945 MEDCouplingMesh *MEDCoupling1SGTUMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 946 {
 947   if(other->getType()!=SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED)
 948     throw INTERP_KERNEL::Exception("Merge of umesh only available with umesh single static geo type each other !");
 949   const MEDCoupling1SGTUMesh *otherC=static_cast<const MEDCoupling1SGTUMesh *>(other);
 950   return Merge1SGTUMeshes(this,otherC);
 951 }
 952
 953 MEDCouplingUMesh *MEDCoupling1SGTUMesh::buildUnstructured() const
 954 {
 955   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(getName(),getMeshDimension());
 956   ret->setCoords(getCoords());
 957   const int *nodalConn=_conn->begin();
 958   int nbCells=getNumberOfCells();
 959   int nbNodesPerCell=getNumberOfNodesPerCell();
 960   int geoType=(int)getCellModelEnum();
 961   MCAuto<DataArrayInt> c=DataArrayInt::New(); c->alloc(nbCells*(nbNodesPerCell+1),1);
 962   int *cPtr=c->getPointer();
 963   for(int i=0;i<nbCells;i++,nodalConn+=nbNodesPerCell)
 964     {
 965       *cPtr++=geoType;
 966       cPtr=std::copy(nodalConn,nodalConn+nbNodesPerCell,cPtr);
 967     }
 968   MCAuto<DataArrayInt> cI=DataArrayInt::Range(0,(nbCells+1)*(nbNodesPerCell+1),nbNodesPerCell+1);
 969   ret->setConnectivity(c,cI,true);
 970   try
 971   { ret->copyTinyInfoFrom(this); }
 972   catch(INTERP_KERNEL::Exception&) { }
 973   return ret.retn();
 974 }
 975
 976 DataArrayInt *MEDCoupling1SGTUMesh::simplexize(int policy)
 977 {
 978   switch(policy)
 979   {
 980     case 0:
 981       return simplexizePol0();
 982     case 1:
 983       return simplexizePol1();
 984     case (int) INTERP_KERNEL::PLANAR_FACE_5:
 985         return simplexizePlanarFace5();
 986     case (int) INTERP_KERNEL::PLANAR_FACE_6:
 987         return simplexizePlanarFace6();
 988     default:
 989       throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::simplexize : unrecognized policy ! Must be :\n  - 0 or 1 (only available for meshdim=2) \n  - PLANAR_FACE_5, PLANAR_FACE_6  (only for meshdim=3)");
 990   }
 991 }
 992
 993 /// @cond INTERNAL
 994
 995 struct MEDCouplingAccVisit
 996 {
 997   MEDCouplingAccVisit():_new_nb_of_nodes(0) { }
 998   int operator()(int val) { if(val!=-1) return _new_nb_of_nodes++; else return -1; }
 999   int _new_nb_of_nodes;
1000 };
1001
1002 /// @endcond
1003
1004 /*!
1005  * This method returns all node ids used in \b this. The data array returned has to be dealt by the caller.
1006  * The returned node ids are sortes ascendingly. This method is closed to MEDCoupling1SGTUMesh::getNodeIdsInUse except
1007  * the format of returned DataArrayInt instance.
1008  *
1009  * \return a newly allocated DataArrayInt sorted ascendingly of fetched node ids.
1010  * \sa MEDCoupling1SGTUMesh::getNodeIdsInUse, areAllNodesFetched
1011  */
1012 DataArrayInt *MEDCoupling1SGTUMesh::computeFetchedNodeIds() const
1013 {
1014   checkConsistencyOfConnectivity();
1015   int nbNodes(getNumberOfNodes());
1016   std::vector<bool> fetchedNodes(nbNodes,false);
1017   computeNodeIdsAlg(fetchedNodes);
1018   int sz((int)std::count(fetchedNodes.begin(),fetchedNodes.end(),true));
1019   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(sz,1);
1020   int *retPtr(ret->getPointer());
1021   for(int i=0;i<nbNodes;i++)
1022     if(fetchedNodes[i])
1023       *retPtr++=i;
1024   return ret.retn();
1025 }
1026
1027 /*!
1028  * Finds nodes not used in any cell and returns an array giving a new id to every node
1029  * by excluding the unused nodes, for which the array holds -1. The result array is
1030  * a mapping in "Old to New" mode.
1031  *  \param [out] nbrOfNodesInUse - number of node ids present in the nodal connectivity.
1032  *  \return DataArrayInt * - a new instance of DataArrayInt. Its length is \a
1033  *          this->getNumberOfNodes(). It holds for each node of \a this mesh either -1
1034  *          if the node is unused or a new id else. The caller is to delete this
1035  *          array using decrRef() as it is no more needed.
1036  *  \throw If the coordinates array is not set.
1037  *  \throw If the nodal connectivity of cells is not defined.
1038  *  \throw If the nodal connectivity includes an invalid id.
1039  *  \sa MEDCoupling1SGTUMesh::computeFetchedNodeIds, areAllNodesFetched
1040  */
1041 DataArrayInt *MEDCoupling1SGTUMesh::getNodeIdsInUse(int& nbrOfNodesInUse) const
1042 {
1043   nbrOfNodesInUse=-1;
1044   int nbOfNodes=getNumberOfNodes();
1045   int nbOfCells=getNumberOfCells();
1046   MCAuto<DataArrayInt> ret(DataArrayInt::New());
1047   ret->alloc(nbOfNodes,1);
1048   int *traducer=ret->getPointer();
1049   std::fill(traducer,traducer+nbOfNodes,-1);
1050   const int *conn=_conn->begin();
1051   int nbNodesPerCell=getNumberOfNodesPerCell();
1052   for(int i=0;i<nbOfCells;i++)
1053     for(int j=0;j<nbNodesPerCell;j++,conn++)
1054       if(*conn>=0 && *conn<nbOfNodes)
1055         traducer[*conn]=1;
1056       else
1057         {
1058           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getNodeIdsInUse : In cell #" << i  << " presence of node id " <<  conn[j] << " not in [0," << nbOfNodes << ") !";
1059           throw INTERP_KERNEL::Exception(oss.str().c_str());
1060         }
1061   nbrOfNodesInUse=(int)std::count(traducer,traducer+nbOfNodes,1);
1062   std::transform(traducer,traducer+nbOfNodes,traducer,MEDCouplingAccVisit());
1063   return ret.retn();
1064 }
1065
1066 /*!
1067  * This method renumbers only nodal connectivity in \a this. The renumbering is only an offset applied. So this method is a specialization of
1068  * \a renumberNodesInConn. \b WARNING, this method does not check that the resulting node ids in the nodal connectivity is in a valid range !
1069  *
1070  * \param [in] offset - specifies the offset to be applied on each element of connectivity.
1071  *
1072  * \sa renumberNodesInConn
1073  */
1074 void MEDCoupling1SGTUMesh::renumberNodesWithOffsetInConn(int offset)
1075 {
1076   getNumberOfCells();//only to check that all is well defined.
1077   _conn->applyLin(1,offset);
1078   updateTime();
1079 }
1080
1081 /*!
1082  *  Same than renumberNodesInConn(const int *) except that here the format of old-to-new traducer is using map instead
1083  *  of array. This method is dedicated for renumbering from a big set of nodes the a tiny set of nodes which is the case during extraction
1084  *  of a big mesh.
1085  */
1086 void MEDCoupling1SGTUMesh::renumberNodesInConn(const INTERP_KERNEL::HashMap<int,int>& newNodeNumbersO2N)
1087 {
1088   getNumberOfCells();//only to check that all is well defined.
1089   int *begPtr(_conn->getPointer());
1090   int nbElt(_conn->getNumberOfTuples());
1091   int *endPtr(begPtr+nbElt);
1092   for(int *it=begPtr;it!=endPtr;it++)
1093     {
1094       INTERP_KERNEL::HashMap<int,int>::const_iterator it2(newNodeNumbersO2N.find(*it));
1095       if(it2!=newNodeNumbersO2N.end())
1096         {
1097           *it=(*it2).second;
1098         }
1099       else
1100         {
1101           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::renumberNodesInConn : At pos #" << std::distance(begPtr,it) << " of nodal connectivity value is " << *it << ". Not in map !";
1102           throw INTERP_KERNEL::Exception(oss.str().c_str());
1103         }
1104     }
1105   updateTime();
1106 }
1107
1108 /*!
1109  * Changes ids of nodes within the nodal connectivity arrays according to a permutation
1110  * array in "Old to New" mode. The node coordinates array is \b not changed by this method.
1111  * This method is a generalization of shiftNodeNumbersInConn().
1112  *  \warning This method performs no check of validity of new ids. **Use it with care !**
1113  *  \param [in] newNodeNumbersO2N - a permutation array, of length \a
1114  *         this->getNumberOfNodes(), in "Old to New" mode.
1115  *         See \ref numbering for more info on renumbering modes.
1116  *  \throw If the nodal connectivity of cells is not defined.
1117  */
1118 void MEDCoupling1SGTUMesh::renumberNodesInConn(const int *newNodeNumbersO2N)
1119 {
1120   getNumberOfCells();//only to check that all is well defined.
1121   _conn->transformWithIndArr(newNodeNumbersO2N,newNodeNumbersO2N+getNumberOfNodes());
1122   updateTime();
1123 }
1124
1125 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshes(const MEDCoupling1SGTUMesh *mesh1, const MEDCoupling1SGTUMesh *mesh2)
1126 {
1127   std::vector<const MEDCoupling1SGTUMesh *> tmp(2);
1128   tmp[0]=const_cast<MEDCoupling1SGTUMesh *>(mesh1); tmp[1]=const_cast<MEDCoupling1SGTUMesh *>(mesh2);
1129   return Merge1SGTUMeshes(tmp);
1130 }
1131
1132 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshes(std::vector<const MEDCoupling1SGTUMesh *>& a)
1133 {
1134   std::size_t sz=a.size();
1135   if(sz==0)
1136     return Merge1SGTUMeshesLL(a);
1137   for(std::size_t ii=0;ii<sz;ii++)
1138     if(!a[ii])
1139       {
1140         std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::Merge1SGTUMeshes : item #" << ii << " in input array of size "<< sz << " is empty !";
1141         throw INTERP_KERNEL::Exception(oss.str().c_str());
1142       }
1143   const INTERP_KERNEL::CellModel *cm=&(a[0]->getCellModel());
1144   for(std::size_t ii=0;ii<sz;ii++)
1145     if(&(a[ii]->getCellModel())!=cm)
1146       throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshes : all items must have the same geo type !");
1147   std::vector< MCAuto<MEDCoupling1SGTUMesh> > bb(sz);
1148   std::vector< const MEDCoupling1SGTUMesh * > aa(sz);
1149   int spaceDim=-3;
1150   for(std::size_t i=0;i<sz && spaceDim==-3;i++)
1151     {
1152       const MEDCoupling1SGTUMesh *cur=a[i];
1153       const DataArrayDouble *coo=cur->getCoords();
1154       if(coo)
1155         spaceDim=coo->getNumberOfComponents();
1156     }
1157   if(spaceDim==-3)
1158     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshes : no spaceDim specified ! unable to perform merge !");
1159   for(std::size_t i=0;i<sz;i++)
1160     {
1161       bb[i]=a[i]->buildSetInstanceFromThis(spaceDim);
1162       aa[i]=bb[i];
1163     }
1164   return Merge1SGTUMeshesLL(aa);
1165 }
1166
1167 /*!
1168  * \throw If presence of a null instance in the input vector \a a.
1169  * \throw If a is empty
1170  */
1171 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords(std::vector<const MEDCoupling1SGTUMesh *>& a)
1172 {
1173   if(a.empty())
1174     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : input array must be NON EMPTY !");
1175   std::vector<const MEDCoupling1SGTUMesh *>::const_iterator it=a.begin();
1176   if(!(*it))
1177     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : null instance in the first element of input vector !");
1178   std::vector<const DataArrayInt *> ncs(a.size());
1179   (*it)->getNumberOfCells();//to check that all is OK
1180   const DataArrayDouble *coords=(*it)->getCoords();
1181   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
1182   ncs[0]=(*it)->getNodalConnectivity();
1183   it++;
1184   for(int i=1;it!=a.end();i++,it++)
1185     {
1186       if(!(*it))
1187         throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : presence of a null instance in the input vector !");
1188       if(cm!=&((*it)->getCellModel()))
1189         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1SGTUMeshes impossible !");
1190       (*it)->getNumberOfCells();//to check that all is OK
1191       ncs[i]=(*it)->getNodalConnectivity();
1192       if(coords!=(*it)->getCoords())
1193         throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : not lying on same coords !");
1194     }
1195   MCAuto<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh("merge",*cm));
1196   ret->setCoords(coords);
1197   ret->_conn=DataArrayInt::Aggregate(ncs);
1198   return ret.retn();
1199 }
1200
1201 /*!
1202  * Assume that all instances in \a a are non null. If null it leads to a crash. That's why this method is assigned to be low level (LL)
1203  */
1204 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshesLL(std::vector<const MEDCoupling1SGTUMesh *>& a)
1205 {
1206   if(a.empty())
1207     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshes : input array must be NON EMPTY !");
1208   std::vector<const MEDCoupling1SGTUMesh *>::const_iterator it=a.begin();
1209   int nbOfCells=(*it)->getNumberOfCells();
1210   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
1211   int nbNodesPerCell=(*it)->getNumberOfNodesPerCell();
1212   it++;
1213   for(;it!=a.end();it++)
1214     {
1215       if(cm!=&((*it)->getCellModel()))
1216         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1SGTUMeshes impossible !");
1217       nbOfCells+=(*it)->getNumberOfCells();
1218     }
1219   std::vector<const MEDCouplingPointSet *> aps(a.size());
1220   std::copy(a.begin(),a.end(),aps.begin());
1221   MCAuto<DataArrayDouble> pts=MergeNodesArray(aps);
1222   MCAuto<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh("merge",*cm));
1223   ret->setCoords(pts);
1224   MCAuto<DataArrayInt> c=DataArrayInt::New();
1225   c->alloc(nbOfCells*nbNodesPerCell,1);
1226   int *cPtr=c->getPointer();
1227   int offset=0;
1228   for(it=a.begin();it!=a.end();it++)
1229     {
1230       int curConnLgth=(*it)->getNodalConnectivityLength();
1231       const int *curC=(*it)->_conn->begin();
1232       cPtr=std::transform(curC,curC+curConnLgth,cPtr,std::bind2nd(std::plus<int>(),offset));
1233       offset+=(*it)->getNumberOfNodes();
1234     }
1235   //
1236   ret->setNodalConnectivity(c);
1237   return ret.retn();
1238 }
1239
1240 MEDCouplingPointSet *MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoords(const int *begin, const int *end) const
1241 {
1242   int ncell=getNumberOfCells();
1243   MCAuto<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh(getName(),*_cm));
1244   ret->setCoords(_coords);
1245   std::size_t nbOfElemsRet=std::distance(begin,end);
1246   const int *inConn=_conn->getConstPointer();
1247   int sz=getNumberOfNodesPerCell();
1248   MCAuto<DataArrayInt> connRet=DataArrayInt::New(); connRet->alloc((int)nbOfElemsRet*sz,1);
1249   int *connPtr=connRet->getPointer();
1250   for(const int *work=begin;work!=end;work++,connPtr+=sz)
1251     {
1252       if(*work>=0 && *work<ncell)
1253         std::copy(inConn+(work[0])*sz,inConn+(work[0]+1)*sz,connPtr);
1254       else
1255         {
1256           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoords : On pos #" << std::distance(begin,work) << " input cell id =" << *work << " should be in [0," << ncell << ") !";
1257           throw INTERP_KERNEL::Exception(oss.str().c_str());
1258         }
1259     }
1260   ret->_conn=connRet;
1261   ret->copyTinyInfoFrom(this);
1262   return ret.retn();
1263 }
1264
1265 MEDCouplingPointSet *MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoordsSlice(int start, int end, int step) const
1266 {
1267   int ncell=getNumberOfCells();
1268   int nbOfElemsRet=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoordsSlice : ");
1269   MCAuto<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh(getName(),*_cm));
1270   ret->setCoords(_coords);
1271   const int *inConn=_conn->getConstPointer();
1272   int sz=getNumberOfNodesPerCell();
1273   MCAuto<DataArrayInt> connRet=DataArrayInt::New(); connRet->alloc((int)nbOfElemsRet*sz,1);
1274   int *connPtr=connRet->getPointer();
1275   int curId=start;
1276   for(int i=0;i<nbOfElemsRet;i++,connPtr+=sz,curId+=step)
1277     {
1278       if(curId>=0 && curId<ncell)
1279         std::copy(inConn+curId*sz,inConn+(curId+1)*sz,connPtr);
1280       else
1281         {
1282           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoordsSlice : On pos #" << i << " input cell id =" << curId  << " should be in [0," << ncell << ") !";
1283           throw INTERP_KERNEL::Exception(oss.str().c_str());
1284         }
1285     }
1286   ret->_conn=connRet;
1287   ret->copyTinyInfoFrom(this);
1288   return ret.retn();
1289 }
1290
1291 void MEDCoupling1SGTUMesh::computeNodeIdsAlg(std::vector<bool>& nodeIdsInUse) const
1292 {
1293   int sz((int)nodeIdsInUse.size());
1294   for(const int *conn=_conn->begin();conn!=_conn->end();conn++)
1295     {
1296       if(*conn>=0 && *conn<sz)
1297        nodeIdsInUse[*conn]=true;
1298       else
1299         {
1300           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::computeFetchedNodeIds : At pos #" << std::distance(_conn->begin(),conn) << " value is " << *conn << " must be in [0," << sz << ") !";
1301           throw INTERP_KERNEL::Exception(oss.str().c_str());
1302         }
1303     }
1304 }
1305
1306 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::buildSetInstanceFromThis(int spaceDim) const
1307 {
1308   MCAuto<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh(getName(),*_cm));
1309   MCAuto<DataArrayInt> tmp1;
1310   const DataArrayInt *nodalConn(_conn);
1311   if(!nodalConn)
1312     {
1313       tmp1=DataArrayInt::New(); tmp1->alloc(0,1);
1314     }
1315   else
1316     tmp1=_conn;
1317   ret->_conn=tmp1;
1318   if(!_coords)
1319     {
1320       MCAuto<DataArrayDouble> coords=DataArrayDouble::New(); coords->alloc(0,spaceDim);
1321       ret->setCoords(coords);
1322     }
1323   else
1324     ret->setCoords(_coords);
1325   return ret.retn();
1326 }
1327
1328 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePol0()
1329 {
1330   int nbOfCells=getNumberOfCells();
1331   if(getCellModelEnum()!=INTERP_KERNEL::NORM_QUAD4)
1332     return DataArrayInt::Range(0,nbOfCells,1);
1333   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(2*3*nbOfCells,1);
1334   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(2*nbOfCells,1);
1335   const int *c(_conn->begin());
1336   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1337   for(int i=0;i<nbOfCells;i++,c+=4,newConnPtr+=6,retPtr+=2)
1338     {
1339       newConnPtr[0]=c[0]; newConnPtr[1]=c[1]; newConnPtr[2]=c[2];
1340       newConnPtr[3]=c[0]; newConnPtr[4]=c[2]; newConnPtr[5]=c[3];
1341       retPtr[0]=i; retPtr[1]=i;
1342     }
1343   _conn=newConn;
1344   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TRI3);
1345   updateTime();
1346   return ret.retn();
1347 }
1348
1349 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePol1()
1350 {
1351   int nbOfCells=getNumberOfCells();
1352   if(getCellModelEnum()!=INTERP_KERNEL::NORM_QUAD4)
1353     return DataArrayInt::Range(0,nbOfCells,1);
1354   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(2*3*nbOfCells,1);
1355   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(2*nbOfCells,1);
1356   const int *c(_conn->begin());
1357   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1358   for(int i=0;i<nbOfCells;i++,c+=4,newConnPtr+=6,retPtr+=2)
1359     {
1360       newConnPtr[0]=c[0]; newConnPtr[1]=c[1]; newConnPtr[2]=c[3];
1361       newConnPtr[3]=c[1]; newConnPtr[4]=c[2]; newConnPtr[5]=c[3];
1362       retPtr[0]=i; retPtr[1]=i;
1363     }
1364   _conn=newConn;
1365   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TRI3);
1366   updateTime();
1367   return ret.retn();
1368 }
1369
1370 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePlanarFace5()
1371 {
1372   int nbOfCells=getNumberOfCells();
1373   if(getCellModelEnum()!=INTERP_KERNEL::NORM_HEXA8)
1374     return DataArrayInt::Range(0,nbOfCells,1);
1375   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(5*4*nbOfCells,1);
1376   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(5*nbOfCells,1);
1377   const int *c(_conn->begin());
1378   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1379   for(int i=0;i<nbOfCells;i++,c+=8,newConnPtr+=20,retPtr+=5)
1380     {
1381       for(int j=0;j<20;j++)
1382         newConnPtr[j]=c[INTERP_KERNEL::SPLIT_NODES_5_WO[j]];
1383       retPtr[0]=i; retPtr[1]=i; retPtr[2]=i; retPtr[3]=i; retPtr[4]=i;
1384     }
1385   _conn=newConn;
1386   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TETRA4);
1387   updateTime();
1388   return ret.retn();
1389 }
1390
1391 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePlanarFace6()
1392 {
1393   int nbOfCells=getNumberOfCells();
1394   if(getCellModelEnum()!=INTERP_KERNEL::NORM_HEXA8)
1395     return DataArrayInt::Range(0,nbOfCells,1);
1396   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(6*4*nbOfCells,1);
1397   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(6*nbOfCells,1);
1398   const int *c(_conn->begin());
1399   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1400   for(int i=0;i<nbOfCells;i++,c+=8,newConnPtr+=24,retPtr+=6)
1401     {
1402       for(int j=0;j<24;j++)
1403         newConnPtr[j]=c[INTERP_KERNEL::SPLIT_NODES_6_WO[j]];
1404       retPtr[0]=i; retPtr[1]=i; retPtr[2]=i; retPtr[3]=i; retPtr[4]=i; retPtr[5]=i;
1405     }
1406   _conn=newConn;
1407   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TETRA4);
1408   updateTime();
1409   return ret.retn();
1410 }
1411
1412 void MEDCoupling1SGTUMesh::reprQuickOverview(std::ostream& stream) const
1413 {
1414   stream << "MEDCoupling1SGTUMesh C++ instance at " << this << ". Type=" << _cm->getRepr() << ". Name : \"" << getName() << "\".";
1415   stream << " Mesh dimension : " << getMeshDimension() << ".";
1416   if(!_coords)
1417     { stream << " No coordinates set !"; return ; }
1418   if(!_coords->isAllocated())
1419     { stream << " Coordinates set but not allocated !"; return ; }
1420   stream << " Space dimension : " << _coords->getNumberOfComponents() << "." << std::endl;
1421   stream << "Number of nodes : " << _coords->getNumberOfTuples() << ".";
1422   if(!(const DataArrayInt *)_conn)
1423     { stream << std::endl << "Nodal connectivity NOT set !"; return ; }
1424   if(_conn->isAllocated())
1425     {
1426       if(_conn->getNumberOfComponents()==1)
1427         stream << std::endl << "Number of cells : " << getNumberOfCells() << ".";
1428     }
1429 }
1430
1431 void MEDCoupling1SGTUMesh::checkFullyDefined() const
1432 {
1433   if(!((const DataArrayInt *)_conn) || !((const DataArrayDouble *)_coords))
1434     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFullyDefined : part of this is not fully defined.");
1435 }
1436
1437 /*!
1438  * First step of unserialization process.
1439  */
1440 bool MEDCoupling1SGTUMesh::isEmptyMesh(const std::vector<int>& tinyInfo) const
1441 {
1442   throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::isEmptyMesh : not implemented yet !");
1443 }
1444
1445 void MEDCoupling1SGTUMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
1446 {
1447   int it,order;
1448   double time=getTime(it,order);
1449   tinyInfo.clear(); tinyInfoD.clear(); littleStrings.clear();
1450   //
1451   littleStrings.push_back(getName());
1452   littleStrings.push_back(getDescription());
1453   littleStrings.push_back(getTimeUnit());
1454   //
1455   std::vector<std::string> littleStrings2,littleStrings3;
1456   if((const DataArrayDouble *)_coords)
1457     _coords->getTinySerializationStrInformation(littleStrings2);
1458   if((const DataArrayInt *)_conn)
1459     _conn->getTinySerializationStrInformation(littleStrings3);
1460   int sz0((int)littleStrings2.size()),sz1((int)littleStrings3.size());
1461   littleStrings.insert(littleStrings.end(),littleStrings2.begin(),littleStrings2.end());
1462   littleStrings.insert(littleStrings.end(),littleStrings3.begin(),littleStrings3.end());
1463   //
1464   tinyInfo.push_back(getCellModelEnum());
1465   tinyInfo.push_back(it);
1466   tinyInfo.push_back(order);
1467   std::vector<int> tinyInfo2,tinyInfo3;
1468   if((const DataArrayDouble *)_coords)
1469     _coords->getTinySerializationIntInformation(tinyInfo2);
1470   if((const DataArrayInt *)_conn)
1471     _conn->getTinySerializationIntInformation(tinyInfo3);
1472   int sz2((int)tinyInfo2.size()),sz3((int)tinyInfo3.size());
1473   tinyInfo.push_back(sz0); tinyInfo.push_back(sz1); tinyInfo.push_back(sz2); tinyInfo.push_back(sz3);
1474   tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
1475   tinyInfo.insert(tinyInfo.end(),tinyInfo3.begin(),tinyInfo3.end());
1476   //
1477   tinyInfoD.push_back(time);
1478 }
1479
1480 void MEDCoupling1SGTUMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
1481 {
1482   std::vector<int> tinyInfo2(tinyInfo.begin()+7,tinyInfo.begin()+7+tinyInfo[5]);
1483   std::vector<int> tinyInfo1(tinyInfo.begin()+7+tinyInfo[5],tinyInfo.begin()+7+tinyInfo[5]+tinyInfo[6]);
1484   a1->resizeForUnserialization(tinyInfo1);
1485   a2->resizeForUnserialization(tinyInfo2);
1486 }
1487
1488 void MEDCoupling1SGTUMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
1489 {
1490   int sz(0);
1491   if((const DataArrayInt *)_conn)
1492     if(_conn->isAllocated())
1493       sz=_conn->getNbOfElems();
1494   a1=DataArrayInt::New();
1495   a1->alloc(sz,1);
1496   if(sz!=0 && (const DataArrayInt *)_conn)
1497     std::copy(_conn->begin(),_conn->end(),a1->getPointer());
1498   sz=0;
1499   if((const DataArrayDouble *)_coords)
1500     if(_coords->isAllocated())
1501       sz=_coords->getNbOfElems();
1502   a2=DataArrayDouble::New();
1503   a2->alloc(sz,1);
1504   if(sz!=0 && (const DataArrayDouble *)_coords)
1505     std::copy(_coords->begin(),_coords->end(),a2->getPointer());
1506 }
1507
1508 void MEDCoupling1SGTUMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
1509                                            const std::vector<std::string>& littleStrings)
1510 {
1511   INTERP_KERNEL::NormalizedCellType gt((INTERP_KERNEL::NormalizedCellType)tinyInfo[0]);
1512   _cm=&INTERP_KERNEL::CellModel::GetCellModel(gt);
1513   setName(littleStrings[0]);
1514   setDescription(littleStrings[1]);
1515   setTimeUnit(littleStrings[2]);
1516   setTime(tinyInfoD[0],tinyInfo[1],tinyInfo[2]);
1517   int sz0(tinyInfo[3]),sz1(tinyInfo[4]),sz2(tinyInfo[5]),sz3(tinyInfo[6]);
1518   //
1519   _coords=DataArrayDouble::New();
1520   std::vector<int> tinyInfo2(tinyInfo.begin()+7,tinyInfo.begin()+7+sz2);
1521   _coords->resizeForUnserialization(tinyInfo2);
1522   std::copy(a2->begin(),a2->end(),_coords->getPointer());
1523   _conn=DataArrayInt::New();
1524   std::vector<int> tinyInfo3(tinyInfo.begin()+7+sz2,tinyInfo.begin()+7+sz2+sz3);
1525   _conn->resizeForUnserialization(tinyInfo3);
1526   std::copy(a1->begin(),a1->end(),_conn->getPointer());
1527   std::vector<std::string> littleStrings2(littleStrings.begin()+3,littleStrings.begin()+3+sz0);
1528   _coords->finishUnserialization(tinyInfo2,littleStrings2);
1529   std::vector<std::string> littleStrings3(littleStrings.begin()+3+sz0,littleStrings.begin()+3+sz0+sz1);
1530   _conn->finishUnserialization(tinyInfo3,littleStrings3);
1531 }
1532
1533 /*!
1534  * Checks if \a this and \a other meshes are geometrically equivalent with high
1535  * probability, else an exception is thrown. The meshes are considered equivalent if
1536  * (1) meshes contain the same number of nodes and the same number of elements of the
1537  * same types (2) three cells of the two meshes (first, last and middle) are based
1538  * on coincident nodes (with a specified precision).
1539  *  \param [in] other - the mesh to compare with.
1540  *  \param [in] prec - the precision used to compare nodes of the two meshes.
1541  *  \throw If the two meshes do not match.
1542  */
1543 void MEDCoupling1SGTUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const
1544 {
1545   MEDCouplingPointSet::checkFastEquivalWith(other,prec);
1546   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
1547   if(!otherC)
1548     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : Two meshes are not unstructured with single static geometric type !");
1549   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
1550   if(c1==c2)
1551     return;
1552   if(!c1 || !c2)
1553     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : presence of nodal connectivity only in one of the 2 meshes !");
1554   if((c1->isAllocated() && !c2->isAllocated()) || (!c1->isAllocated() && c2->isAllocated()))
1555     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : in nodal connectivity, only one is allocated !");
1556   if(c1->getNumberOfComponents()!=1 || c1->getNumberOfComponents()!=1)
1557     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : in nodal connectivity, must have 1 and only 1 component !");
1558   if(c1->getHashCode()!=c2->getHashCode())
1559     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : nodal connectivity differs");
1560 }
1561
1562 MEDCouplingPointSet *MEDCoupling1SGTUMesh::mergeMyselfWithOnSameCoords(const MEDCouplingPointSet *other) const
1563 {
1564   if(!other)
1565     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::mergeMyselfWithOnSameCoords : input other is null !");
1566   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
1567   if(!otherC)
1568     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::mergeMyselfWithOnSameCoords : the input other mesh is not of type single statuc geo type unstructured !");
1569   std::vector<const MEDCoupling1SGTUMesh *> ms(2);
1570   ms[0]=this;
1571   ms[1]=otherC;
1572   return Merge1SGTUMeshesOnSameCoords(ms);
1573 }
1574
1575 void MEDCoupling1SGTUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const
1576 {
1577   checkFullyDefined();
1578   int nbOfNodes=getNumberOfNodes();
1579   int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int));
1580   revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1);
1581   std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
1582   const int *conn=_conn->begin();
1583   int nbOfCells=getNumberOfCells();
1584   int nbOfEltsInRevNodal=0;
1585   int nbOfNodesPerCell=getNumberOfNodesPerCell();
1586   for(int eltId=0;eltId<nbOfCells;eltId++)
1587     {
1588       for(int j=0;j<nbOfNodesPerCell;j++,conn++)
1589         {
1590           if(conn[0]>=0 && conn[0]<nbOfNodes)
1591             {
1592               nbOfEltsInRevNodal++;
1593               revNodalIndxPtr[conn[0]+1]++;
1594             }
1595           else
1596             {
1597               std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getReverseNodalConnectivity : At cell #" << eltId << " presence of nodeId #" << conn[0] << " should be in [0," << nbOfNodes << ") !";
1598               throw INTERP_KERNEL::Exception(oss.str().c_str());
1599             }
1600         }
1601     }
1602   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<int>());
1603   conn=_conn->begin();
1604   int *revNodalPtr=(int *)malloc((nbOfEltsInRevNodal)*sizeof(int));
1605   revNodal->useArray(revNodalPtr,true,C_DEALLOC,nbOfEltsInRevNodal,1);
1606   std::fill(revNodalPtr,revNodalPtr+nbOfEltsInRevNodal,-1);
1607   for(int eltId=0;eltId<nbOfCells;eltId++)
1608     {
1609       for(int j=0;j<nbOfNodesPerCell;j++,conn++)
1610         {
1611           *std::find_if(revNodalPtr+revNodalIndxPtr[*conn],revNodalPtr+revNodalIndxPtr[*conn+1],std::bind2nd(std::equal_to<int>(),-1))=eltId;
1612         }
1613     }
1614 }
1615
1616 /*!
1617  * Use \a nodalConn array as nodal connectivity of \a this. The input \a nodalConn pointer can be null.
1618  */
1619 void MEDCoupling1SGTUMesh::setNodalConnectivity(DataArrayInt *nodalConn)
1620 {
1621   if(nodalConn)
1622     nodalConn->incrRef();
1623   _conn=nodalConn;
1624   declareAsNew();
1625 }
1626
1627 /*!
1628  * \return DataArrayInt * - the internal reference to the nodal connectivity. The caller is not reponsible to deallocate it.
1629  */
1630 DataArrayInt *MEDCoupling1SGTUMesh::getNodalConnectivity() const
1631 {
1632   const DataArrayInt *ret(_conn);
1633   return const_cast<DataArrayInt *>(ret);
1634 }
1635
1636 /*!
1637  * Allocates memory to store an estimation of the given number of cells. Closer is the estimation to the number of cells effectively inserted,
1638  * less will be the needs to realloc. If the number of cells to be inserted is not known simply put 0 to this parameter.
1639  * If a nodal connectivity previouly existed before the call of this method, it will be reset.
1640  *
1641  *  \param [in] nbOfCells - estimation of the number of cell \a this mesh will contain.
1642  */
1643 void MEDCoupling1SGTUMesh::allocateCells(int nbOfCells)
1644 {
1645   if(nbOfCells<0)
1646     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::allocateCells : the input number of cells should be >= 0 !");
1647   _conn=DataArrayInt::New();
1648   _conn->reserve(getNumberOfNodesPerCell()*nbOfCells);
1649   declareAsNew();
1650 }
1651
1652 /*!
1653  * Appends at the end of \a this a cell having nodal connectivity array defined in [ \a nodalConnOfCellBg, \a nodalConnOfCellEnd ).
1654  *
1655  * \param [in] nodalConnOfCellBg - the begin (included) of nodal connectivity of the cell to add.
1656  * \param [in] nodalConnOfCellEnd - the end (excluded) of nodal connectivity of the cell to add.
1657  * \throw If the length of the input nodal connectivity array of the cell to add is not equal to number of nodes per cell relative to the unique geometric type
1658  *        attached to \a this.
1659  * \thow If the nodal connectivity array in \a this is null (call MEDCoupling1SGTUMesh::allocateCells before).
1660  */
1661 void MEDCoupling1SGTUMesh::insertNextCell(const int *nodalConnOfCellBg, const int *nodalConnOfCellEnd)
1662 {
1663   int sz=(int)std::distance(nodalConnOfCellBg,nodalConnOfCellEnd);
1664   int ref=getNumberOfNodesPerCell();
1665   if(sz==ref)
1666     {
1667       DataArrayInt *c(_conn);
1668       if(c)
1669         c->pushBackValsSilent(nodalConnOfCellBg,nodalConnOfCellEnd);
1670       else
1671         throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::insertNextCell : nodal connectivity array is null ! Call MEDCoupling1SGTUMesh::allocateCells before !");
1672     }
1673   else
1674     {
1675       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::insertNextCell : input nodal size (" << sz << ") does not match number of nodes per cell of this (";
1676       oss << ref << ") !";
1677       throw INTERP_KERNEL::Exception(oss.str().c_str());
1678     }
1679 }
1680
1681 /*!
1682  * This method builds the dual mesh of \a this and returns it.
1683  *
1684  * \return MEDCoupling1SGTUMesh * - newly object created to be managed by the caller.
1685  * \throw If \a this is not a mesh containing only simplex cells.
1686  * \throw If \a this is not correctly allocated (coordinates and connectivities have to be correctly set !).
1687  * \throw If at least one node in \a this is orphan (without any simplex cell lying on it !)
1688  */
1689 MEDCoupling1GTUMesh *MEDCoupling1SGTUMesh::computeDualMesh() const
1690 {
1691   const INTERP_KERNEL::CellModel& cm(getCellModel());
1692   if(!cm.isSimplex())
1693     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::computeDualMesh : this mesh is not a simplex mesh ! Please invoke simplexize of tetrahedrize on this before calling this method !");
1694   switch(getMeshDimension())
1695   {
1696     case 3:
1697       return computeDualMesh3D();
1698     case 2:
1699       return computeDualMesh2D();
1700     default:
1701       throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::computeDualMesh : meshdimension must be in [2,3] !");
1702   }
1703 }
1704
1705 /*!
1706  * This method explode each NORM_HEXA8 cells in \a this into 6 NORM_QUAD4 cells and put the result into the MEDCoupling1SGTUMesh returned instance.
1707  *
1708  * \return MEDCoupling1SGTUMesh * - a newly allocated instances (to be managed by the caller) storing the result of the explosion.
1709  * \throw If \a this is not a mesh containing only NORM_HEXA8 cells.
1710  * \throw If \a this is not properly allocated.
1711  */
1712 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::explodeEachHexa8To6Quad4() const
1713 {
1714   const INTERP_KERNEL::CellModel& cm(getCellModel());
1715   if(cm.getEnum()!=INTERP_KERNEL::NORM_HEXA8)
1716     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::explodeEachHexa8To6Quad4 : this method can be applied only on HEXA8 mesh !");
1717   int nbHexa8(getNumberOfCells());
1718   const int *inConnPtr(getNodalConnectivity()->begin());
1719   MCAuto<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName(),INTERP_KERNEL::NORM_QUAD4));
1720   MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc(nbHexa8*6*4,1);
1721   int *cPtr(c->getPointer());
1722   for(int i=0;i<nbHexa8;i++,inConnPtr+=8)
1723     {
1724       for(int j=0;j<6;j++,cPtr+=4)
1725         cm.fillSonCellNodalConnectivity(j,inConnPtr,cPtr);
1726     }
1727   ret->setCoords(getCoords());
1728   ret->setNodalConnectivity(c);
1729   return ret.retn();
1730 }
1731
1732 /*!
1733  * This method starts from an unstructured mesh that hides in reality a cartesian mesh.
1734  * If it is not the case, an exception will be thrown.
1735  * This method returns three objects : The cartesian mesh geometrically equivalent to \a this (within a precision of \a eps) and a permutation of cells
1736  * and a permutation of nodes.
1737  *
1738  * - this[cellPerm[i]]=ret[i]
1739  *
1740  * \param [out] cellPerm the permutation array of size \c this->getNumberOfCells()
1741  * \param [out] nodePerm the permutation array of size \c this->getNumberOfNodes()
1742  * \return MEDCouplingCMesh * - a newly allocated mesh that is the result of the structurization of \a this.
1743  */
1744 MEDCouplingCMesh *MEDCoupling1SGTUMesh::structurizeMe(DataArrayInt *& cellPerm, DataArrayInt *& nodePerm, double eps) const
1745 {
1746   checkConsistencyLight();
1747   int spaceDim(getSpaceDimension()),meshDim(getMeshDimension()),nbNodes(getNumberOfNodes());
1748   if(MEDCouplingStructuredMesh::GetGeoTypeGivenMeshDimension(meshDim)!=getCellModelEnum())
1749     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::structurizeMe : the unique geo type in this is not compatible with the geometric type regarding mesh dimension !");
1750   MCAuto<MEDCouplingCMesh> cm(MEDCouplingCMesh::New());
1751   for(int i=0;i<spaceDim;i++)
1752     {
1753       std::vector<int> tmp(1,i);
1754       MCAuto<DataArrayDouble> elt(static_cast<DataArrayDouble*>(getCoords()->keepSelectedComponents(tmp)));
1755       elt=elt->getDifferentValues(eps);
1756       elt->sort(true);
1757       cm->setCoordsAt(i,elt);
1758     }
1759   if(nbNodes!=cm->getNumberOfNodes())
1760     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::structurizeMe : considering the number of nodes after split per components in space this can't be a cartesian mesh ! Maybe your epsilon parameter is invalid ?");
1761   try
1762   { cm->copyTinyInfoFrom(this); }
1763   catch(INTERP_KERNEL::Exception&) { }
1764   MCAuto<MEDCouplingUMesh> um(cm->buildUnstructured()),self(buildUnstructured());
1765   self->checkGeoEquivalWith(um,12,eps,cellPerm,nodePerm);
1766   return cm.retn();
1767 }
1768
1769 /// @cond INTERNAL
1770
1771 bool UpdateHexa8Cell(int validAxis, int neighId, const int *validConnQuad4NeighSide, int *allFacesNodalConn, int *myNeighbours)
1772 {
1773   static const int TAB[48]={
1774     0,1,2,3,4,5,6,7,//0
1775     4,7,6,5,0,3,2,1,//1
1776     0,3,7,4,1,2,6,5,//2
1777     4,0,3,7,5,1,2,6,//3
1778     5,1,0,4,6,2,3,7,//4
1779     3,7,4,0,2,6,5,1 //5
1780   };
1781   static const int TAB2[6]={0,0,3,3,3,3};
1782   if(myNeighbours[validAxis]==neighId && allFacesNodalConn[4*validAxis+0]==validConnQuad4NeighSide[TAB2[validAxis]])
1783     return true;
1784   int oldAxis((int)std::distance(myNeighbours,std::find(myNeighbours,myNeighbours+6,neighId)));
1785   std::size_t pos(std::distance(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS,std::find(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS,MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS+6,oldAxis)));
1786   std::size_t pos0(pos/2),pos1(pos%2);
1787   int oldAxisOpp(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS[2*pos0+(pos1+1)%2]);
1788   int oldConn[8],myConn2[8]={-1,-1,-1,-1,-1,-1,-1,-1},myConn[8],edgeConn[2],allFacesTmp[24],neighTmp[6];
1789   oldConn[0]=allFacesNodalConn[0]; oldConn[1]=allFacesNodalConn[1]; oldConn[2]=allFacesNodalConn[2]; oldConn[3]=allFacesNodalConn[3];
1790   oldConn[4]=allFacesNodalConn[4]; oldConn[5]=allFacesNodalConn[7]; oldConn[6]=allFacesNodalConn[6]; oldConn[7]=allFacesNodalConn[5];
1791   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_HEXA8));
1792   for(int i=0;i<4;i++)
1793     myConn2[i]=validConnQuad4NeighSide[(4-i+TAB2[validAxis])%4];
1794   for(int i=0;i<4;i++)
1795     {
1796       int nodeId(myConn2[i]);//the node id for which the opposite one will be found
1797       bool found(false);
1798       INTERP_KERNEL::NormalizedCellType typeOfSon;
1799       for(int j=0;j<12 && !found;j++)
1800         {
1801           cm.fillSonEdgesNodalConnectivity3D(j,oldConn,-1,edgeConn,typeOfSon);
1802           if(edgeConn[0]==nodeId || edgeConn[1]==nodeId)
1803             {
1804               if(std::find(allFacesNodalConn+4*oldAxisOpp,allFacesNodalConn+4*oldAxisOpp+4,edgeConn[0]==nodeId?edgeConn[1]:edgeConn[0])!=allFacesNodalConn+4*oldAxisOpp+4)
1805                 {
1806                   myConn2[i+4]=edgeConn[0]==nodeId?edgeConn[1]:edgeConn[0];
1807                   found=true;
1808                 }
1809             }
1810         }
1811       if(!found)
1812         throw INTERP_KERNEL::Exception("UpdateHexa8Cell : Internal Error !");
1813     }
1814   const int *myTab(TAB+8*validAxis);
1815   for(int i=0;i<8;i++)
1816     myConn[i]=myConn2[myTab[i]];
1817   for(int i=0;i<6;i++)
1818     {
1819       cm.fillSonCellNodalConnectivity(i,myConn,allFacesTmp+4*i);
1820       std::set<int> s(allFacesTmp+4*i,allFacesTmp+4*i+4);
1821       bool found(false);
1822       for(int j=0;j<6 && !found;j++)
1823         {
1824           std::set<int> s1(allFacesNodalConn+4*j,allFacesNodalConn+4*j+4);
1825           if(s==s1)
1826             {
1827               neighTmp[i]=myNeighbours[j];
1828               found=true;
1829             }
1830         }
1831       if(!found)
1832         throw INTERP_KERNEL::Exception("UpdateHexa8Cell : Internal Error #2 !");
1833     }
1834   std::copy(allFacesTmp,allFacesTmp+24,allFacesNodalConn);
1835   std::copy(neighTmp,neighTmp+6,myNeighbours);
1836   return false;
1837 }
1838
1839 /// @endcond
1840
1841 /*!
1842  * This method expects the \a this contains NORM_HEXA8 cells only. This method will sort each cells in \a this so that their numbering was
1843  * homogeneous. If it succeeds the result of MEDCouplingUMesh::tetrahedrize will return a conform mesh.
1844  *
1845  * \return DataArrayInt * - a newly allocated array (to be managed by the caller) containing renumbered cell ids.
1846  *
1847  * \throw If \a this is not a mesh containing only NORM_HEXA8 cells.
1848  * \throw If \a this is not properly allocated.
1849  * \sa MEDCouplingUMesh::tetrahedrize, MEDCouplingUMesh::simplexize.
1850  */
1851 DataArrayInt *MEDCoupling1SGTUMesh::sortHexa8EachOther()
1852 {
1853   MCAuto<MEDCoupling1SGTUMesh> quads(explodeEachHexa8To6Quad4());//checks that only hexa8
1854   int nbHexa8(getNumberOfCells()),*cQuads(quads->getNodalConnectivity()->getPointer());
1855   MCAuto<DataArrayInt> neighOfQuads(DataArrayInt::New()); neighOfQuads->alloc(nbHexa8*6,1); neighOfQuads->fillWithValue(-1);
1856   int *ptNeigh(neighOfQuads->getPointer());
1857   {//neighOfQuads tells for each face of each Quad8 which cell (if!=-1) is connected to this face.
1858     MCAuto<MEDCouplingUMesh> quadsTmp(quads->buildUnstructured());
1859     MCAuto<DataArrayInt> ccSafe,cciSafe;
1860     DataArrayInt *cc(0),*cci(0);
1861     quadsTmp->findCommonCells(3,0,cc,cci);
1862     ccSafe=cc; cciSafe=cci;
1863     const int *ccPtr(ccSafe->begin()),nbOfPair(cci->getNumberOfTuples()-1);
1864     for(int i=0;i<nbOfPair;i++)
1865       { ptNeigh[ccPtr[2*i+0]]=ccPtr[2*i+1]/6; ptNeigh[ccPtr[2*i+1]]=ccPtr[2*i+0]/6; }
1866   }
1867   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
1868   std::vector<bool> fetched(nbHexa8,false);
1869   std::vector<bool>::iterator it(std::find(fetched.begin(),fetched.end(),false));
1870   while(it!=fetched.end())//it will turns as time as number of connected zones
1871     {
1872       int cellId((int)std::distance(fetched.begin(),it));//it is the seed of the connected zone.
1873       std::set<int> s; s.insert(cellId);//s contains already organized.
1874       while(!s.empty())
1875         {
1876           std::set<int> sNext;
1877           for(std::set<int>::const_iterator it0=s.begin();it0!=s.end();it0++)
1878             {
1879               fetched[*it0]=true;
1880               int *myNeighb(ptNeigh+6*(*it0));
1881               for(int i=0;i<6;i++)
1882                 {
1883                   if(myNeighb[i]!=-1 && !fetched[myNeighb[i]])
1884                     {
1885                       std::size_t pos(std::distance(HEXA8_FACE_PAIRS,std::find(HEXA8_FACE_PAIRS,HEXA8_FACE_PAIRS+6,i)));
1886                       std::size_t pos0(pos/2),pos1(pos%2);
1887                       if(!UpdateHexa8Cell(HEXA8_FACE_PAIRS[2*pos0+(pos1+1)%2],*it0,cQuads+6*4*(*it0)+4*i,cQuads+6*4*myNeighb[i],ptNeigh+6*myNeighb[i]))
1888                         ret->pushBackSilent(myNeighb[i]);
1889                       fetched[myNeighb[i]]=true;
1890                       sNext.insert(myNeighb[i]);
1891                     }
1892                 }
1893             }
1894           s=sNext;
1895         }
1896       it=std::find(fetched.begin(),fetched.end(),false);
1897     }
1898   if(!ret->empty())
1899     {
1900       int *conn(getNodalConnectivity()->getPointer());
1901       for(const int *pt=ret->begin();pt!=ret->end();pt++)
1902         {
1903           int cellId(*pt);
1904           conn[8*cellId+0]=cQuads[24*cellId+0]; conn[8*cellId+1]=cQuads[24*cellId+1]; conn[8*cellId+2]=cQuads[24*cellId+2]; conn[8*cellId+3]=cQuads[24*cellId+3];
1905           conn[8*cellId+4]=cQuads[24*cellId+4]; conn[8*cellId+5]=cQuads[24*cellId+7]; conn[8*cellId+6]=cQuads[24*cellId+6]; conn[8*cellId+7]=cQuads[24*cellId+5];
1906         }
1907       declareAsNew();
1908     }
1909   return ret.retn();
1910 }
1911
1912 MEDCoupling1DGTUMesh *MEDCoupling1SGTUMesh::computeDualMesh3D() const
1913 {
1914   static const int DUAL_TETRA_0[36]={
1915     4,1,0, 6,0,3, 7,3,1,
1916     4,0,1, 5,2,0, 8,1,2,
1917     6,3,0, 5,0,2, 9,2,3,
1918     7,1,3, 9,3,2, 8,2,1
1919   };
1920   static const int DUAL_TETRA_1[36]={
1921     8,4,10, 11,5,8, 10,7,11,
1922     9,4,8, 8,5,12, 12,6,9,
1923     10,4,9, 9,6,13, 13,7,10,
1924     12,5,11, 13,6,12, 11,7,13
1925   };
1926   static const int FACEID_NOT_SH_NODE[4]={2,3,1,0};
1927   if(getCellModelEnum()!=INTERP_KERNEL::NORM_TETRA4)
1928     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::computeDualMesh3D : only TETRA4 supported !");
1929   checkFullyDefined();
1930   MCAuto<MEDCouplingUMesh> thisu(buildUnstructured());
1931   MCAuto<DataArrayInt> revNodArr(DataArrayInt::New()),revNodIArr(DataArrayInt::New());
1932   thisu->getReverseNodalConnectivity(revNodArr,revNodIArr);
1933   const int *revNod(revNodArr->begin()),*revNodI(revNodIArr->begin()),*nodal(_conn->begin());
1934   MCAuto<DataArrayInt> d1Arr(DataArrayInt::New()),di1Arr(DataArrayInt::New()),rd1Arr(DataArrayInt::New()),rdi1Arr(DataArrayInt::New());
1935   MCAuto<MEDCouplingUMesh> edges(thisu->explode3DMeshTo1D(d1Arr,di1Arr,rd1Arr,rdi1Arr));
1936   const int *d1(d1Arr->begin());
1937   MCAuto<DataArrayInt> d2Arr(DataArrayInt::New()),di2Arr(DataArrayInt::New()),rd2Arr(DataArrayInt::New()),rdi2Arr(DataArrayInt::New());
1938   MCAuto<MEDCouplingUMesh> faces(thisu->buildDescendingConnectivity(d2Arr,di2Arr,rd2Arr,rdi2Arr));  thisu=0;
1939   const int *d2(d2Arr->begin()),*rdi2(rdi2Arr->begin());
1940   MCAuto<DataArrayDouble> edgesBaryArr(edges->computeCellCenterOfMass()),facesBaryArr(faces->computeCellCenterOfMass()),baryArr(computeCellCenterOfMass());
1941   const int nbOfNodes(getNumberOfNodes()),offset0(nbOfNodes+faces->getNumberOfCells()),offset1(offset0+edges->getNumberOfCells());
1942   edges=0; faces=0;
1943   std::vector<const DataArrayDouble *> v(4); v[0]=getCoords(); v[1]=facesBaryArr; v[2]=edgesBaryArr; v[3]=baryArr;
1944   MCAuto<DataArrayDouble> zeArr(DataArrayDouble::Aggregate(v)); baryArr=0; edgesBaryArr=0; facesBaryArr=0;
1945   std::string name("DualOf_"); name+=getName();
1946   MCAuto<MEDCoupling1DGTUMesh> ret(MEDCoupling1DGTUMesh::New(name,INTERP_KERNEL::NORM_POLYHED)); ret->setCoords(zeArr);
1947   MCAuto<DataArrayInt> cArr(DataArrayInt::New()),ciArr(DataArrayInt::New()); ciArr->alloc(nbOfNodes+1,1); ciArr->setIJ(0,0,0); cArr->alloc(0,1);
1948   for(int i=0;i<nbOfNodes;i++,revNodI++)
1949     {
1950       int nbOfCellsSharingNode(revNodI[1]-revNodI[0]);
1951       if(nbOfCellsSharingNode==0)
1952         {
1953           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::computeDualMesh3D : Node #" << i << " is orphan !";
1954           throw INTERP_KERNEL::Exception(oss.str().c_str());
1955         }
1956       for(int j=0;j<nbOfCellsSharingNode;j++)
1957         {
1958           int curCellId(revNod[revNodI[0]+j]);
1959           const int *connOfCurCell(nodal+4*curCellId);
1960           std::size_t nodePosInCurCell(std::distance(connOfCurCell,std::find(connOfCurCell,connOfCurCell+4,i)));
1961           if(j!=0) cArr->pushBackSilent(-1);
1962           int tmp[14];
1963           //
1964           tmp[0]=d1[6*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+0]-4]+offset0; tmp[1]=d2[4*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+1]]+nbOfNodes;
1965           tmp[2]=curCellId+offset1; tmp[3]=d2[4*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+2]]+nbOfNodes;
1966           tmp[4]=-1;
1967           tmp[5]=d1[6*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+3]-4]+offset0; tmp[6]=d2[4*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+4]]+nbOfNodes;
1968           tmp[7]=curCellId+offset1; tmp[8]=d2[4*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+5]]+nbOfNodes;
1969           tmp[9]=-1;
1970           tmp[10]=d1[6*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+6]-4]+offset0; tmp[11]=d2[4*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+7]]+nbOfNodes;
1971           tmp[12]=curCellId+offset1; tmp[13]=d2[4*curCellId+DUAL_TETRA_0[nodePosInCurCell*9+8]]+nbOfNodes;
1972           cArr->insertAtTheEnd(tmp,tmp+14);
1973           int kk(0);
1974           for(int k=0;k<4;k++)
1975             {
1976               if(FACEID_NOT_SH_NODE[nodePosInCurCell]!=k)
1977                 {
1978                   const int *faceId(d2+4*curCellId+k);
1979                   if(rdi2[*faceId+1]-rdi2[*faceId]==1)
1980                     {
1981                       int tmp2[5]; tmp2[0]=-1; tmp2[1]=i;
1982                       tmp2[2]=d1[6*curCellId+DUAL_TETRA_1[9*nodePosInCurCell+3*kk+0]-8]+offset0;
1983                       tmp2[3]=d2[4*curCellId+DUAL_TETRA_1[9*nodePosInCurCell+3*kk+1]-4]+nbOfNodes;
1984                       tmp2[4]=d1[6*curCellId+DUAL_TETRA_1[9*nodePosInCurCell+3*kk+2]-8]+offset0;
1985                       cArr->insertAtTheEnd(tmp2,tmp2+5);
1986                     }
1987                   kk++;
1988                 }
1989             }
1990         }
1991       ciArr->setIJ(i+1,0,cArr->getNumberOfTuples());
1992     }
1993   ret->setNodalConnectivity(cArr,ciArr);
1994   return ret.retn();
1995 }
1996
1997 MEDCoupling1DGTUMesh *MEDCoupling1SGTUMesh::computeDualMesh2D() const
1998 {
1999   static const int DUAL_TRI_0[6]={0,2, 1,0, 2,1};
2000   static const int DUAL_TRI_1[6]={-3,+5, +3,-4, +4,-5};
2001   static const int FACEID_NOT_SH_NODE[3]={1,2,0};
2002   if(getCellModelEnum()!=INTERP_KERNEL::NORM_TRI3)
2003     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::computeDualMesh2D : only TRI3 supported !");
2004   checkFullyDefined();
2005   MCAuto<MEDCouplingUMesh> thisu(buildUnstructured());
2006   MCAuto<DataArrayInt> revNodArr(DataArrayInt::New()),revNodIArr(DataArrayInt::New());
2007   thisu->getReverseNodalConnectivity(revNodArr,revNodIArr);
2008   const int *revNod(revNodArr->begin()),*revNodI(revNodIArr->begin()),*nodal(_conn->begin());
2009   MCAuto<DataArrayInt> d2Arr(DataArrayInt::New()),di2Arr(DataArrayInt::New()),rd2Arr(DataArrayInt::New()),rdi2Arr(DataArrayInt::New());
2010   MCAuto<MEDCouplingUMesh> edges(thisu->buildDescendingConnectivity(d2Arr,di2Arr,rd2Arr,rdi2Arr));  thisu=0;
2011   const int *d2(d2Arr->begin()),*rdi2(rdi2Arr->begin());
2012   MCAuto<DataArrayDouble> edgesBaryArr(edges->computeCellCenterOfMass()),baryArr(computeCellCenterOfMass());
2013   const int nbOfNodes(getNumberOfNodes()),offset0(nbOfNodes+edges->getNumberOfCells());
2014   edges=0;
2015   std::vector<const DataArrayDouble *> v(3); v[0]=getCoords(); v[1]=edgesBaryArr; v[2]=baryArr;
2016   MCAuto<DataArrayDouble> zeArr(DataArrayDouble::Aggregate(v)); baryArr=0; edgesBaryArr=0;
2017   std::string name("DualOf_"); name+=getName();
2018   MCAuto<MEDCoupling1DGTUMesh> ret(MEDCoupling1DGTUMesh::New(name,INTERP_KERNEL::NORM_POLYGON)); ret->setCoords(zeArr);
2019   MCAuto<DataArrayInt> cArr(DataArrayInt::New()),ciArr(DataArrayInt::New()); ciArr->alloc(nbOfNodes+1,1); ciArr->setIJ(0,0,0); cArr->alloc(0,1);
2020   for(int i=0;i<nbOfNodes;i++,revNodI++)
2021     {
2022       int nbOfCellsSharingNode(revNodI[1]-revNodI[0]);
2023       if(nbOfCellsSharingNode==0)
2024         {
2025           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::computeDualMesh2D : Node #" << i << " is orphan !";
2026           throw INTERP_KERNEL::Exception(oss.str().c_str());
2027         }
2028       std::vector< std::vector<int> > polyg;
2029       for(int j=0;j<nbOfCellsSharingNode;j++)
2030         {
2031           int curCellId(revNod[revNodI[0]+j]);
2032           const int *connOfCurCell(nodal+3*curCellId);
2033           std::size_t nodePosInCurCell(std::distance(connOfCurCell,std::find(connOfCurCell,connOfCurCell+4,i)));
2034           std::vector<int> locV(3);
2035           locV[0]=d2[3*curCellId+DUAL_TRI_0[2*nodePosInCurCell+0]]+nbOfNodes; locV[1]=curCellId+offset0; locV[2]=d2[3*curCellId+DUAL_TRI_0[2*nodePosInCurCell+1]]+nbOfNodes;
2036           polyg.push_back(locV);
2037           int kk(0);
2038           for(int k=0;k<3;k++)
2039             {
2040               if(FACEID_NOT_SH_NODE[nodePosInCurCell]!=k)
2041                 {
2042                   const int *edgeId(d2+3*curCellId+k);
2043                   if(rdi2[*edgeId+1]-rdi2[*edgeId]==1)
2044                     {
2045                       std::vector<int> locV2(2);
2046                       int zeLocEdgeIdRel(DUAL_TRI_1[2*nodePosInCurCell+kk]);
2047                       if(zeLocEdgeIdRel>0)
2048                         {  locV2[0]=d2[3*curCellId+zeLocEdgeIdRel-3]+nbOfNodes;  locV2[1]=i; }
2049                       else
2050                         {  locV2[0]=i; locV2[1]=d2[3*curCellId-zeLocEdgeIdRel-3]+nbOfNodes; }
2051                       polyg.push_back(locV2);
2052                     }
2053                   kk++;
2054                 }
2055             }
2056         }
2057       std::vector<int> zePolyg(MEDCoupling1DGTUMesh::BuildAPolygonFromParts(polyg));
2058       cArr->insertAtTheEnd(zePolyg.begin(),zePolyg.end());
2059       ciArr->setIJ(i+1,0,cArr->getNumberOfTuples());
2060     }
2061   ret->setNodalConnectivity(cArr,ciArr);
2062   return ret.retn();
2063 }
2064
2065 /*!
2066  * This method aggregate the bbox of each cell and put it into bbox
2067  *
2068  * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
2069  *                         For all other cases this input parameter is ignored.
2070  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
2071  *
2072  * \throw If \a this is not fully set (coordinates and connectivity).
2073  * \throw If a cell in \a this has no valid nodeId.
2074  */
2075 DataArrayDouble *MEDCoupling1SGTUMesh::getBoundingBoxForBBTree(double arcDetEps) const
2076 {
2077   int spaceDim(getSpaceDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()),nbOfNodesPerCell(getNumberOfNodesPerCell());
2078   MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
2079   double *bbox(ret->getPointer());
2080   for(int i=0;i<nbOfCells*spaceDim;i++)
2081     {
2082       bbox[2*i]=std::numeric_limits<double>::max();
2083       bbox[2*i+1]=-std::numeric_limits<double>::max();
2084     }
2085   const double *coordsPtr(_coords->getConstPointer());
2086   const int *conn(_conn->getConstPointer());
2087   for(int i=0;i<nbOfCells;i++)
2088     {
2089       int kk(0);
2090       for(int j=0;j<nbOfNodesPerCell;j++,conn++)
2091         {
2092           int nodeId(*conn);
2093           if(nodeId>=0 && nodeId<nbOfNodes)
2094             {
2095               for(int k=0;k<spaceDim;k++)
2096                 {
2097                   bbox[2*spaceDim*i+2*k]=std::min(bbox[2*spaceDim*i+2*k],coordsPtr[spaceDim*nodeId+k]);
2098                   bbox[2*spaceDim*i+2*k+1]=std::max(bbox[2*spaceDim*i+2*k+1],coordsPtr[spaceDim*nodeId+k]);
2099                 }
2100               kk++;
2101             }
2102         }
2103       if(kk==0)
2104         {
2105           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getBoundingBoxForBBTree : cell #" << i << " contains no valid nodeId !";
2106           throw INTERP_KERNEL::Exception(oss.str().c_str());
2107         }
2108     }
2109   return ret.retn();
2110 }
2111
2112 /*!
2113  * Returns the cell field giving for each cell in \a this its diameter. Diameter means the max length of all possible SEG2 in the cell.
2114  *
2115  * \return a new instance of field containing the result. The returned instance has to be deallocated by the caller.
2116  */
2117 MEDCouplingFieldDouble *MEDCoupling1SGTUMesh::computeDiameterField() const
2118 {
2119   checkFullyDefined();
2120   MCAuto<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME));
2121   int nbCells(getNumberOfCells());
2122   MCAuto<DataArrayDouble> arr(DataArrayDouble::New());
2123   arr->alloc(nbCells,1);
2124   INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::DiameterCalculator> dc(_cm->buildInstanceOfDiameterCalulator(getSpaceDimension()));
2125   dc->computeFor1SGTUMeshFrmt(nbCells,_conn->begin(),getCoords()->begin(),arr->getPointer());
2126   ret->setMesh(this);
2127   ret->setArray(arr);
2128   ret->setName("Diameter");
2129   return ret.retn();
2130 }
2131
2132 //==
2133
2134 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::New()
2135 {
2136   return new MEDCoupling1DGTUMesh;
2137 }
2138
2139 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::New(const std::string& name, INTERP_KERNEL::NormalizedCellType type)
2140 {
2141   if(type==INTERP_KERNEL::NORM_ERROR)
2142     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::New : NORM_ERROR is not a valid type to be used as base geometric type for a mesh !");
2143   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
2144   if(!cm.isDynamic())
2145     {
2146       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::New : the input geometric type " << cm.getRepr() << " is static ! Only dynamic types are allowed here !";
2147       throw INTERP_KERNEL::Exception(oss.str().c_str());
2148     }
2149   return new MEDCoupling1DGTUMesh(name,cm);
2150 }
2151
2152 MEDCoupling1DGTUMesh::MEDCoupling1DGTUMesh()
2153 {
2154 }
2155
2156 MEDCoupling1DGTUMesh::MEDCoupling1DGTUMesh(const std::string& name, const INTERP_KERNEL::CellModel& cm):MEDCoupling1GTUMesh(name,cm)
2157 {
2158 }
2159
2160 MEDCoupling1DGTUMesh::MEDCoupling1DGTUMesh(const MEDCoupling1DGTUMesh& other, bool recDeepCpy):MEDCoupling1GTUMesh(other,recDeepCpy),_conn_indx(other._conn_indx),_conn(other._conn)
2161 {
2162   if(recDeepCpy)
2163     {
2164       const DataArrayInt *c(other._conn);
2165       if(c)
2166         _conn=c->deepCopy();
2167       c=other._conn_indx;
2168       if(c)
2169         _conn_indx=c->deepCopy();
2170     }
2171 }
2172
2173 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::clone(bool recDeepCpy) const
2174 {
2175   return new MEDCoupling1DGTUMesh(*this,recDeepCpy);
2176 }
2177
2178 /*!
2179  * This method behaves mostly like MEDCoupling1DGTUMesh::deepCopy method, except that only nodal connectivity arrays are deeply copied.
2180  * The coordinates are shared between \a this and the returned instance.
2181  *
2182  * \return MEDCoupling1DGTUMesh * - A new object instance holding the copy of \a this (deep for connectivity, shallow for coordiantes)
2183  * \sa MEDCoupling1DGTUMesh::deepCopy
2184  */
2185 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::deepCopyConnectivityOnly() const
2186 {
2187   checkConsistencyLight();
2188   MCAuto<MEDCoupling1DGTUMesh> ret(clone(false));
2189   MCAuto<DataArrayInt> c(_conn->deepCopy()),ci(_conn_indx->deepCopy());
2190   ret->setNodalConnectivity(c,ci);
2191   return ret.retn();
2192 }
2193
2194 void MEDCoupling1DGTUMesh::updateTime() const
2195 {
2196   MEDCoupling1GTUMesh::updateTime();
2197   const DataArrayInt *c(_conn);
2198   if(c)
2199     updateTimeWith(*c);
2200   c=_conn_indx;
2201   if(c)
2202     updateTimeWith(*c);
2203 }
2204
2205 std::size_t MEDCoupling1DGTUMesh::getHeapMemorySizeWithoutChildren() const
2206 {
2207   return MEDCoupling1GTUMesh::getHeapMemorySizeWithoutChildren();
2208 }
2209
2210 std::vector<const BigMemoryObject *> MEDCoupling1DGTUMesh::getDirectChildrenWithNull() const
2211 {
2212   std::vector<const BigMemoryObject *> ret(MEDCoupling1GTUMesh::getDirectChildrenWithNull());
2213   ret.push_back((const DataArrayInt *)_conn);
2214   ret.push_back((const DataArrayInt *)_conn_indx);
2215   return ret;
2216 }
2217
2218 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::deepCopy() const
2219 {
2220   return clone(true);
2221 }
2222
2223 bool MEDCoupling1DGTUMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
2224 {
2225   if(!other)
2226     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::isEqualIfNotWhy : input other pointer is null !");
2227   std::ostringstream oss; oss.precision(15);
2228   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
2229   if(!otherC)
2230     {
2231       reason="mesh given in input is not castable in MEDCoupling1DGTUMesh !";
2232       return false;
2233     }
2234   if(!MEDCoupling1GTUMesh::isEqualIfNotWhy(other,prec,reason))
2235     return false;
2236   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
2237   if(c1==c2)
2238     return true;
2239   if(!c1 || !c2)
2240     {
2241       reason="in connectivity of single dynamic geometric type exactly one among this and other is null !";
2242       return false;
2243     }
2244   if(!c1->isEqualIfNotWhy(*c2,reason))
2245     {
2246       reason.insert(0,"Nodal connectivity DataArrayInt differs : ");
2247       return false;
2248     }
2249   c1=_conn_indx; c2=otherC->_conn_indx;
2250   if(c1==c2)
2251     return true;
2252   if(!c1 || !c2)
2253     {
2254       reason="in connectivity index of single dynamic geometric type exactly one among this and other is null !";
2255       return false;
2256     }
2257   if(!c1->isEqualIfNotWhy(*c2,reason))
2258     {
2259       reason.insert(0,"Nodal connectivity index DataArrayInt differs : ");
2260       return false;
2261     }
2262   return true;
2263 }
2264
2265 bool MEDCoupling1DGTUMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
2266 {
2267   if(!other)
2268     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::isEqualWithoutConsideringStr : input other pointer is null !");
2269   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
2270   if(!otherC)
2271     return false;
2272   if(!MEDCoupling1GTUMesh::isEqualWithoutConsideringStr(other,prec))
2273     return false;
2274   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
2275   if(c1==c2)
2276     return true;
2277   if(!c1 || !c2)
2278     return false;
2279   if(!c1->isEqualWithoutConsideringStr(*c2))
2280     return false;
2281   return true;
2282   c1=_conn_indx; c2=otherC->_conn_indx;
2283   if(c1==c2)
2284     return true;
2285   if(!c1 || !c2)
2286     return false;
2287   if(!c1->isEqualWithoutConsideringStr(*c2))
2288     return false;
2289   return true;
2290 }
2291
2292 /*!
2293  * Checks if \a this and \a other meshes are geometrically equivalent with high
2294  * probability, else an exception is thrown. The meshes are considered equivalent if
2295  * (1) meshes contain the same number of nodes and the same number of elements of the
2296  * same types (2) three cells of the two meshes (first, last and middle) are based
2297  * on coincident nodes (with a specified precision).
2298  *  \param [in] other - the mesh to compare with.
2299  *  \param [in] prec - the precision used to compare nodes of the two meshes.
2300  *  \throw If the two meshes do not match.
2301  */
2302 void MEDCoupling1DGTUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const
2303 {
2304   MEDCouplingPointSet::checkFastEquivalWith(other,prec);
2305   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
2306   if(!otherC)
2307     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : Two meshes are not unstructured with single dynamic geometric type !");
2308   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
2309   if(c1!=c2)
2310     {
2311       if(!c1 || !c2)
2312         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : presence of nodal connectivity only in one of the 2 meshes !");
2313       if((c1->isAllocated() && !c2->isAllocated()) || (!c1->isAllocated() && c2->isAllocated()))
2314         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity, only one is allocated !");
2315       if(c1->getNumberOfComponents()!=1 || c1->getNumberOfComponents()!=1)
2316         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity, must have 1 and only 1 component !");
2317       if(c1->getHashCode()!=c2->getHashCode())
2318         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : nodal connectivity differs");
2319     }
2320   c1=_conn_indx; c2=otherC->_conn_indx;
2321   if(c1!=c2)
2322     {
2323       if(!c1 || !c2)
2324         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : presence of nodal connectivity index only in one of the 2 meshes !");
2325       if((c1->isAllocated() && !c2->isAllocated()) || (!c1->isAllocated() && c2->isAllocated()))
2326         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity index, only one is allocated !");
2327       if(c1->getNumberOfComponents()!=1 || c1->getNumberOfComponents()!=1)
2328         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity index, must have 1 and only 1 component !");
2329       if(c1->getHashCode()!=c2->getHashCode())
2330         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : nodal connectivity index differs");
2331     }
2332 }
2333
2334 void MEDCoupling1DGTUMesh::checkConsistencyOfConnectivity() const
2335 {
2336   const DataArrayInt *c1(_conn);
2337   if(c1)
2338     {
2339       if(c1->getNumberOfComponents()!=1)
2340         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to be with number of components set to one !");
2341       if(c1->getInfoOnComponent(0)!="")
2342         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to have no info on its single component !");
2343       c1->checkAllocated();
2344     }
2345   else
2346     throw INTERP_KERNEL::Exception("Nodal connectivity array not defined !");
2347   //
2348   int sz2=_conn->getNumberOfTuples();
2349   c1=_conn_indx;
2350   if(c1)
2351     {
2352       if(c1->getNumberOfComponents()!=1)
2353         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to be with number of components set to one !");
2354       c1->checkAllocated();
2355       if(c1->getNumberOfTuples()<1)
2356         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to have a a size of 1 at least !");
2357       if(c1->getInfoOnComponent(0)!="")
2358         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to have no info on its single component !");
2359       int f=c1->front(),ll=c1->back();
2360       if(f<0 || f>=sz2)
2361         {
2362           std::ostringstream oss; oss << "Nodal connectivity index array first value (" << f << ") is expected to be exactly in [0," << sz2 << ") !";
2363           throw INTERP_KERNEL::Exception(oss.str().c_str());
2364         }
2365       if(ll<0 || ll>sz2)
2366         {
2367           std::ostringstream oss; oss << "Nodal connectivity index array last value (" << ll << ") is expected to be exactly in [0," << sz2 << "] !";
2368           throw INTERP_KERNEL::Exception(oss.str().c_str());
2369         }
2370       if(f>ll)
2371         {
2372           std::ostringstream oss; oss << "Nodal connectivity index array looks very bad (not increasing monotonic) because front (" << f << ") is greater that back (" << ll << ") !";
2373           throw INTERP_KERNEL::Exception(oss.str().c_str());
2374         }
2375     }
2376   else
2377     throw INTERP_KERNEL::Exception("Nodal connectivity index array not defined !");
2378   int szOfC1Exp=_conn_indx->back();
2379   if(sz2<szOfC1Exp)
2380     {
2381       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::checkConsistencyOfConnectivity : The expected length of nodal connectivity array regarding index is " << szOfC1Exp << " but the actual size of it is " << c1->getNumberOfTuples() << " !";
2382       throw INTERP_KERNEL::Exception(oss.str().c_str());
2383     }
2384 }
2385
2386 /*!
2387  * If \a this pass this method, you are sure that connectivity arrays are not null, with exactly one component, no name, no component name, allocated.
2388  * In addition you are sure that the length of nodal connectivity index array is bigger than or equal to one.
2389  * In addition you are also sure that length of nodal connectivity is coherent with the content of the last value in the index array.
2390  */
2391 void MEDCoupling1DGTUMesh::checkConsistencyLight() const
2392 {
2393   MEDCouplingPointSet::checkConsistencyLight();
2394   checkConsistencyOfConnectivity();
2395 }
2396
2397 void MEDCoupling1DGTUMesh::checkConsistency(double eps) const
2398 {
2399   checkConsistencyLight();
2400   const DataArrayInt *c1(_conn),*c2(_conn_indx);
2401   if(!c2->isMonotonic(true))
2402     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkConsistency : the nodal connectivity index is expected to be increasing monotinic !");
2403   //
2404   int nbOfTuples=c1->getNumberOfTuples();
2405   int nbOfNodes=getNumberOfNodes();
2406   const int *w(c1->begin());
2407   for(int i=0;i<nbOfTuples;i++,w++)
2408     {
2409       if(*w==-1) continue;
2410       if(*w<0 || *w>=nbOfNodes)
2411         {
2412           std::ostringstream oss; oss << "At pos #" << i << " of nodal connectivity array references to node id #" << *w << " must be in [0," << nbOfNodes << ") !";
2413           throw INTERP_KERNEL::Exception(oss.str().c_str());
2414         }
2415     }
2416 }
2417
2418 int MEDCoupling1DGTUMesh::getNumberOfCells() const
2419 {
2420   checkConsistencyOfConnectivity();//do not remove
2421   return _conn_indx->getNumberOfTuples()-1;
2422 }
2423
2424 /*!
2425  * This method returns a newly allocated array containing this->getNumberOfCells() tuples and 1 component.
2426  * For each cell in \b this the number of nodes constituting cell is computed.
2427  * For each polyhedron cell, the sum of the number of nodes of each face constituting polyhedron cell is returned.
2428  * So for pohyhedrons some nodes can be counted several times in the returned result.
2429  *
2430  * \return a newly allocated array
2431  */
2432 DataArrayInt *MEDCoupling1DGTUMesh::computeNbOfNodesPerCell() const
2433 {
2434   checkConsistencyLight();
2435   _conn_indx->checkMonotonic(true);
2436   if(getCellModelEnum()!=INTERP_KERNEL::NORM_POLYHED)
2437     return _conn_indx->deltaShiftIndex();
2438   // for polyhedrons
2439   int nbOfCells=_conn_indx->getNumberOfTuples()-1;
2440   MCAuto<DataArrayInt> ret=DataArrayInt::New();
2441   ret->alloc(nbOfCells,1);
2442   int *retPtr=ret->getPointer();
2443   const int *ci=_conn_indx->begin(),*c=_conn->begin();
2444   for(int i=0;i<nbOfCells;i++,retPtr++,ci++)
2445     *retPtr=ci[1]-ci[0]-std::count(c+ci[0],c+ci[1],-1);
2446   return ret.retn();
2447 }
2448
2449 /*!
2450  * This method returns a newly allocated array containing this->getNumberOfCells() tuples and 1 component.
2451  * For each cell in \b this the number of faces constituting (entity of dimension this->getMeshDimension()-1) cell is computed.
2452  *
2453  * \return a newly allocated array
2454  */
2455 DataArrayInt *MEDCoupling1DGTUMesh::computeNbOfFacesPerCell() const
2456 {
2457   checkConsistencyLight();
2458   _conn_indx->checkMonotonic(true);
2459   if(getCellModelEnum()!=INTERP_KERNEL::NORM_POLYHED && getCellModelEnum()!=INTERP_KERNEL::NORM_QPOLYG)
2460     return _conn_indx->deltaShiftIndex();
2461   if(getCellModelEnum()==INTERP_KERNEL::NORM_QPOLYG)
2462     {
2463       MCAuto<DataArrayInt> ret=_conn_indx->deltaShiftIndex();
2464       ret->applyDivideBy(2);
2465       return ret.retn();
2466     }
2467   // for polyhedrons
2468   int nbOfCells=_conn_indx->getNumberOfTuples()-1;
2469   MCAuto<DataArrayInt> ret=DataArrayInt::New();
2470   ret->alloc(nbOfCells,1);
2471   int *retPtr=ret->getPointer();
2472   const int *ci=_conn_indx->begin(),*c=_conn->begin();
2473   for(int i=0;i<nbOfCells;i++,retPtr++,ci++)
2474     *retPtr=std::count(c+ci[0],c+ci[1],-1)+1;
2475   return ret.retn();
2476 }
2477
2478 /*!
2479  * This method computes effective number of nodes per cell. That is to say nodes appearing several times in nodal connectivity of a cell,
2480  * will be counted only once here whereas it will be counted several times in MEDCoupling1DGTUMesh::computeNbOfNodesPerCell method.
2481  *
2482  * \return DataArrayInt * - new object to be deallocated by the caller.
2483  * \sa MEDCoupling1DGTUMesh::computeNbOfNodesPerCell
2484  */
2485 DataArrayInt *MEDCoupling1DGTUMesh::computeEffectiveNbOfNodesPerCell() const
2486 {
2487   checkConsistencyLight();
2488   _conn_indx->checkMonotonic(true);
2489   int nbOfCells(_conn_indx->getNumberOfTuples()-1);
2490   MCAuto<DataArrayInt> ret=DataArrayInt::New();
2491   ret->alloc(nbOfCells,1);
2492   int *retPtr(ret->getPointer());
2493   const int *ci(_conn_indx->begin()),*c(_conn->begin());
2494   if(getCellModelEnum()!=INTERP_KERNEL::NORM_POLYHED)
2495     {
2496       for(int i=0;i<nbOfCells;i++,retPtr++,ci++)
2497         {
2498           std::set<int> s(c+ci[0],c+ci[1]);
2499           *retPtr=(int)s.size();
2500         }
2501     }
2502   else
2503     {
2504       for(int i=0;i<nbOfCells;i++,retPtr++,ci++)
2505         {
2506           std::set<int> s(c+ci[0],c+ci[1]); s.erase(-1);
2507           *retPtr=(int)s.size();
2508         }
2509     }
2510   return ret.retn();
2511 }
2512
2513 void MEDCoupling1DGTUMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
2514 {
2515   int nbOfCells(getNumberOfCells());//performs checks
2516   if(cellId>=0 && cellId<nbOfCells)
2517     {
2518       int strt=_conn_indx->getIJ(cellId,0),stp=_conn_indx->getIJ(cellId+1,0);
2519       int nbOfNodes=stp-strt;
2520       if(nbOfNodes<0)
2521         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::getNodeIdsOfCell : the index array is invalid ! Should be increasing monotonic !");
2522       conn.resize(nbOfNodes);
2523       std::copy(_conn->begin()+strt,_conn->begin()+stp,conn.begin());
2524     }
2525   else
2526     {
2527       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getNodeIdsOfCell : request for cellId #" << cellId << " must be in [0," << nbOfCells << ") !";
2528       throw INTERP_KERNEL::Exception(oss.str().c_str());
2529     }
2530 }
2531
2532 int MEDCoupling1DGTUMesh::getNumberOfNodesInCell(int cellId) const
2533 {
2534   int nbOfCells(getNumberOfCells());//performs checks
2535   if(cellId>=0 && cellId<nbOfCells)
2536     {
2537       const int *conn(_conn->begin());
2538       int strt=_conn_indx->getIJ(cellId,0),stp=_conn_indx->getIJ(cellId+1,0);
2539       return stp-strt-std::count(conn+strt,conn+stp,-1);
2540     }
2541   else
2542     {
2543       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getNumberOfNodesInCell : request for cellId #" << cellId << " must be in [0," << nbOfCells << ") !";
2544       throw INTERP_KERNEL::Exception(oss.str().c_str());
2545     }
2546 }
2547
2548 std::string MEDCoupling1DGTUMesh::simpleRepr() const
2549 {
2550   static const char msg0[]="No coordinates specified !";
2551   std::ostringstream ret;
2552   ret << "Single dynamic geometic type (" << _cm->getRepr() << ") unstructured mesh with name : \"" << getName() << "\"\n";
2553   ret << "Description of mesh : \"" << getDescription() << "\"\n";
2554   int tmpp1,tmpp2;
2555   double tt=getTime(tmpp1,tmpp2);
2556   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
2557   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
2558   ret << "Mesh dimension : " << getMeshDimension() << "\nSpace dimension : ";
2559   if(_coords!=0)
2560     {
2561       const int spaceDim=getSpaceDimension();
2562       ret << spaceDim << "\nInfo attached on space dimension : ";
2563       for(int i=0;i<spaceDim;i++)
2564         ret << "\"" << _coords->getInfoOnComponent(i) << "\" ";
2565       ret << "\n";
2566     }
2567   else
2568     ret << msg0 << "\n";
2569   ret << "Number of nodes : ";
2570   if(_coords!=0)
2571     ret << getNumberOfNodes() << "\n";
2572   else
2573     ret << msg0 << "\n";
2574   ret << "Number of cells : ";
2575   bool isOK=true;
2576   try { checkConsistencyLight(); } catch(INTERP_KERNEL::Exception& /* e */)
2577   {
2578       ret << "Nodal connectivity arrays are not set or badly set !\n";
2579       isOK=false;
2580   }
2581   if(isOK)
2582     ret << getNumberOfCells() << "\n";
2583   ret << "Cell type : " << _cm->getRepr() << "\n";
2584   return ret.str();
2585 }
2586
2587 std::string MEDCoupling1DGTUMesh::advancedRepr() const
2588 {
2589   std::ostringstream ret;
2590   ret << simpleRepr();
2591   ret << "\nCoordinates array : \n___________________\n\n";
2592   if(_coords)
2593     _coords->reprWithoutNameStream(ret);
2594   else
2595     ret << "No array set !\n";
2596   ret << "\n\nNodal Connectivity : \n____________________\n\n";
2597   //
2598   bool isOK=true;
2599   try { checkConsistency(); } catch(INTERP_KERNEL::Exception& /* e */)
2600   {
2601       ret << "Nodal connectivity arrays are not set or badly set !\n";
2602       isOK=false;
2603   }
2604   if(!isOK)
2605     return ret.str();
2606   int nbOfCells=getNumberOfCells();
2607   const int *ci=_conn_indx->begin(),*c=_conn->begin();
2608   for(int i=0;i<nbOfCells;i++,ci++)
2609     {
2610       ret << "Cell #" << i << " : ";
2611       std::copy(c+ci[0],c+ci[1],std::ostream_iterator<int>(ret," "));
2612       ret << "\n";
2613     }
2614   return ret.str();
2615 }
2616
2617 DataArrayDouble *MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell() const
2618 {
2619   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2620   int spaceDim=getSpaceDimension();
2621   int nbOfCells=getNumberOfCells();//checkConsistencyLight()
2622   int nbOfNodes=getNumberOfNodes();
2623   ret->alloc(nbOfCells,spaceDim);
2624   double *ptToFill=ret->getPointer();
2625   const double *coor=_coords->begin();
2626   const int *nodal=_conn->begin(),*nodali=_conn_indx->begin();
2627   nodal+=nodali[0];
2628   if(getCellModelEnum()!=INTERP_KERNEL::NORM_POLYHED)
2629     {
2630       for(int i=0;i<nbOfCells;i++,ptToFill+=spaceDim,nodali++)
2631         {
2632           std::fill(ptToFill,ptToFill+spaceDim,0.);
2633           if(nodali[0]<nodali[1])// >= to avoid division by 0.
2634             {
2635               for(int j=nodali[0];j<nodali[1];j++,nodal++)
2636                 {
2637                   if(*nodal>=0 && *nodal<nbOfNodes)
2638                     std::transform(coor+spaceDim*nodal[0],coor+spaceDim*(nodal[0]+1),ptToFill,ptToFill,std::plus<double>());
2639                   else
2640                     {
2641                       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of nodeId #" << *nodal << " should be in [0," <<   nbOfNodes << ") !";
2642                       throw INTERP_KERNEL::Exception(oss.str().c_str());
2643                     }
2644                   std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(nodali[1]-nodali[0])));
2645                 }
2646             }
2647           else
2648             {
2649               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell : at cell #" << i << " the nodal index array is invalid !";
2650               throw INTERP_KERNEL::Exception(oss.str().c_str());
2651             }
2652         }
2653     }
2654   else
2655     {
2656       for(int i=0;i<nbOfCells;i++,ptToFill+=spaceDim,nodali++)
2657         {
2658           std::fill(ptToFill,ptToFill+spaceDim,0.);
2659           if(nodali[0]<nodali[1])// >= to avoid division by 0.
2660             {
2661               int nbOfNod=0;
2662               for(int j=nodali[0];j<nodali[1];j++,nodal++)
2663                 {
2664                   if(*nodal==-1) continue;
2665                   if(*nodal>=0 && *nodal<nbOfNodes)
2666                     {
2667                       std::transform(coor+spaceDim*nodal[0],coor+spaceDim*(nodal[0]+1),ptToFill,ptToFill,std::plus<double>());
2668                       nbOfNod++;
2669                     }
2670                   else
2671                     {
2672                       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell (polyhedron) : on cell #" << i << " presence of nodeId #" << *nodal << " should be in [0," <<   nbOfNodes << ") !";
2673                       throw INTERP_KERNEL::Exception(oss.str().c_str());
2674                     }
2675                 }
2676               if(nbOfNod!=0)
2677                 std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./nbOfNod));
2678               else
2679                 {
2680                   std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell (polyhedron) : no nodes in cell #" << i << " !";
2681                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2682                 }
2683             }
2684           else
2685             {
2686               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell (polyhedron)  : at cell #" << i << " the nodal index array is invalid !";
2687               throw INTERP_KERNEL::Exception(oss.str().c_str());
2688             }
2689         }
2690     }
2691   return ret.retn();
2692 }
2693
2694 void MEDCoupling1DGTUMesh::renumberCells(const int *old2NewBg, bool check)
2695 {
2696   int nbCells=getNumberOfCells();
2697   MCAuto<DataArrayInt> o2n=DataArrayInt::New();
2698   o2n->useArray(old2NewBg,false,C_DEALLOC,nbCells,1);
2699   if(check)
2700     o2n=o2n->checkAndPreparePermutation();
2701   //
2702   const int *o2nPtr=o2n->getPointer();
2703   const int *conn=_conn->begin(),*conni=_conn_indx->begin();
2704   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
2705   MCAuto<DataArrayInt> newConnI=DataArrayInt::New();
2706   newConn->alloc(_conn->getNumberOfTuples(),1); newConnI->alloc(nbCells,1);
2707   newConn->copyStringInfoFrom(*_conn); newConnI->copyStringInfoFrom(*_conn_indx);
2708   //
2709   int *newC=newConn->getPointer(),*newCI=newConnI->getPointer();
2710   for(int i=0;i<nbCells;i++)
2711     {
2712       int newPos=o2nPtr[i];
2713       int sz=conni[i+1]-conni[i];
2714       if(sz>=0)
2715         newCI[newPos]=sz;
2716       else
2717         {
2718           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::renumberCells : the index nodal array is invalid for cell #" << i << " !";
2719           throw INTERP_KERNEL::Exception(oss.str().c_str());
2720         }
2721     }
2722   newConnI->computeOffsetsFull(); newCI=newConnI->getPointer();
2723   //
2724   for(int i=0;i<nbCells;i++,conni++)
2725     {
2726       int newp=o2nPtr[i];
2727       std::copy(conn+conni[0],conn+conni[1],newC+newCI[newp]);
2728     }
2729   _conn=newConn;
2730   _conn_indx=newConnI;
2731 }
2732
2733 MEDCouplingMesh *MEDCoupling1DGTUMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
2734 {
2735   if(other->getType()!=SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED)
2736     throw INTERP_KERNEL::Exception("Merge of umesh only available with umesh single dynamic geo type each other !");
2737   const MEDCoupling1DGTUMesh *otherC=static_cast<const MEDCoupling1DGTUMesh *>(other);
2738   return Merge1DGTUMeshes(this,otherC);
2739 }
2740
2741 MEDCouplingUMesh *MEDCoupling1DGTUMesh::buildUnstructured() const
2742 {
2743   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(getName(),getMeshDimension());
2744   ret->setCoords(getCoords());
2745   const int *nodalConn=_conn->begin(),*nodalConnI=_conn_indx->begin();
2746   int nbCells=getNumberOfCells();//checkConsistencyLight
2747   int geoType=(int)getCellModelEnum();
2748   MCAuto<DataArrayInt> c=DataArrayInt::New(); c->alloc(nbCells+_conn->getNumberOfTuples(),1);
2749   MCAuto<DataArrayInt> cI=DataArrayInt::New(); cI->alloc(nbCells+1);
2750   int *cPtr=c->getPointer(),*ciPtr=cI->getPointer();
2751   ciPtr[0]=0;
2752   for(int i=0;i<nbCells;i++,ciPtr++)
2753     {
2754       int sz=nodalConnI[i+1]-nodalConnI[i];
2755       if(sz>=0)
2756         {
2757           *cPtr++=geoType;
2758           cPtr=std::copy(nodalConn+nodalConnI[i],nodalConn+nodalConnI[i+1],cPtr);
2759           ciPtr[1]=ciPtr[0]+sz+1;
2760         }
2761       else
2762         {
2763           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::buildUnstructured : Invalid for nodal index for cell #" << i << " !";
2764           throw INTERP_KERNEL::Exception(oss.str().c_str());
2765         }
2766     }
2767   ret->setConnectivity(c,cI,true);
2768   try
2769   { ret->copyTinyInfoFrom(this); }
2770   catch(INTERP_KERNEL::Exception&) { }
2771   return ret.retn();
2772 }
2773
2774 /*!
2775  * Do nothing for the moment, because there is no policy that allows to split polygons, polyhedrons ... into simplexes
2776  */
2777 DataArrayInt *MEDCoupling1DGTUMesh::simplexize(int policy)
2778 {
2779   int nbOfCells=getNumberOfCells();
2780   MCAuto<DataArrayInt> ret=DataArrayInt::New();
2781   ret->alloc(nbOfCells,1);
2782   ret->iota(0);
2783   return ret.retn();
2784 }
2785
2786 void MEDCoupling1DGTUMesh::reprQuickOverview(std::ostream& stream) const
2787 {
2788   stream << "MEDCoupling1DGTUMesh C++ instance at " << this << ". Type=" << _cm->getRepr() << ". Name : \"" << getName() << "\".";
2789   stream << " Mesh dimension : " << getMeshDimension() << ".";
2790   if(!_coords)
2791     { stream << " No coordinates set !"; return ; }
2792   if(!_coords->isAllocated())
2793     { stream << " Coordinates set but not allocated !"; return ; }
2794   stream << " Space dimension : " << _coords->getNumberOfComponents() << "." << std::endl;
2795   stream << "Number of nodes : " << _coords->getNumberOfTuples() << ".";
2796   bool isOK=true;
2797   try { checkConsistencyLight(); } catch(INTERP_KERNEL::Exception&  /* e */)
2798   {
2799       stream << std::endl << "Nodal connectivity NOT set properly !\n";
2800       isOK=false;
2801   }
2802   if(isOK)
2803     stream << std::endl << "Number of cells : " << getNumberOfCells() << ".";
2804 }
2805
2806 void MEDCoupling1DGTUMesh::shallowCopyConnectivityFrom(const MEDCouplingPointSet *other)
2807 {
2808   if(!other)
2809     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::shallowCopyConnectivityFrom : input pointer is null !");
2810   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
2811   if(!otherC)
2812     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::shallowCopyConnectivityFrom : input pointer is not an MEDCoupling1DGTUMesh instance !");
2813   setNodalConnectivity(otherC->getNodalConnectivity(),otherC->getNodalConnectivityIndex());
2814 }
2815
2816 MEDCouplingPointSet *MEDCoupling1DGTUMesh::mergeMyselfWithOnSameCoords(const MEDCouplingPointSet *other) const
2817 {
2818   if(!other)
2819     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::mergeMyselfWithOnSameCoords : input other is null !");
2820   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
2821   if(!otherC)
2822     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::mergeMyselfWithOnSameCoords : the input other mesh is not of type single statuc geo type unstructured !");
2823   std::vector<const MEDCoupling1DGTUMesh *> ms(2);
2824   ms[0]=this;
2825   ms[1]=otherC;
2826   return Merge1DGTUMeshesOnSameCoords(ms);
2827 }
2828
2829 MEDCouplingPointSet *MEDCoupling1DGTUMesh::buildPartOfMySelfKeepCoords(const int *begin, const int *end) const
2830 {
2831   checkConsistencyLight();
2832   MCAuto<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
2833   ret->setCoords(_coords);
2834   DataArrayInt *c=0,*ci=0;
2835   MEDCouplingUMesh::ExtractFromIndexedArrays(begin,end,_conn,_conn_indx,c,ci);
2836   MCAuto<DataArrayInt> cSafe(c),ciSafe(ci);
2837   ret->setNodalConnectivity(c,ci);
2838   return ret.retn();
2839 }
2840
2841 MEDCouplingPointSet *MEDCoupling1DGTUMesh::buildPartOfMySelfKeepCoordsSlice(int start, int end, int step) const
2842 {
2843   checkConsistencyLight();
2844   MCAuto<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
2845   ret->setCoords(_coords);
2846   DataArrayInt *c=0,*ci=0;
2847   MEDCouplingUMesh::ExtractFromIndexedArraysSlice(start,end,step,_conn,_conn_indx,c,ci);
2848   MCAuto<DataArrayInt> cSafe(c),ciSafe(ci);
2849   ret->setNodalConnectivity(c,ci);
2850   return ret.retn();
2851 }
2852
2853 void MEDCoupling1DGTUMesh::computeNodeIdsAlg(std::vector<bool>& nodeIdsInUse) const
2854 {
2855   checkConsistency();
2856   int sz((int)nodeIdsInUse.size());
2857   for(const int *conn=_conn->begin();conn!=_conn->end();conn++)
2858     {
2859       if(*conn>=0 && *conn<sz)
2860         nodeIdsInUse[*conn]=true;
2861       else
2862         {
2863           if(*conn!=-1)
2864             {
2865               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeNodeIdsAlg : At pos #" << std::distance(_conn->begin(),conn) << " value is " << *conn << " must be in [0," << sz << ") !";
2866               throw INTERP_KERNEL::Exception(oss.str().c_str());
2867             }
2868         }
2869     }
2870 }
2871
2872 void MEDCoupling1DGTUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const
2873 {
2874   checkFullyDefined();
2875   int nbOfNodes=getNumberOfNodes();
2876   int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int));
2877   revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1);
2878   std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
2879   const int *conn=_conn->begin(),*conni=_conn_indx->begin();
2880   int nbOfCells=getNumberOfCells();
2881   int nbOfEltsInRevNodal=0;
2882   for(int eltId=0;eltId<nbOfCells;eltId++)
2883     {
2884       int nbOfNodesPerCell=conni[eltId+1]-conni[eltId];
2885       if(nbOfNodesPerCell>=0)
2886         {
2887           for(int j=0;j<nbOfNodesPerCell;j++)
2888             {
2889               int nodeId=conn[conni[eltId]+j];
2890               if(nodeId==-1) continue;
2891               if(nodeId>=0 && nodeId<nbOfNodes)
2892                 {
2893                   nbOfEltsInRevNodal++;
2894                   revNodalIndxPtr[nodeId+1]++;
2895                 }
2896               else
2897                 {
2898                   std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getReverseNodalConnectivity : At cell #" << eltId << " presence of nodeId #" << conn[0] << " should be in [0," << nbOfNodes << ") !";
2899                   throw INTERP_KERNEL::Exception(oss.str().c_str());
2900                 }
2901             }
2902         }
2903       else
2904         {
2905           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getReverseNodalConnectivity : At cell #" << eltId << "nodal connectivity is invalid !";
2906           throw INTERP_KERNEL::Exception(oss.str().c_str());
2907         }
2908     }
2909   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<int>());
2910   conn=_conn->begin();
2911   int *revNodalPtr=(int *)malloc((nbOfEltsInRevNodal)*sizeof(int));
2912   revNodal->useArray(revNodalPtr,true,C_DEALLOC,nbOfEltsInRevNodal,1);
2913   std::fill(revNodalPtr,revNodalPtr+nbOfEltsInRevNodal,-1);
2914   for(int eltId=0;eltId<nbOfCells;eltId++)
2915     {
2916       int nbOfNodesPerCell=conni[eltId+1]-conni[eltId];
2917       for(int j=0;j<nbOfNodesPerCell;j++)
2918         {
2919           int nodeId=conn[conni[eltId]+j];
2920           if(nodeId!=-1)
2921             *std::find_if(revNodalPtr+revNodalIndxPtr[nodeId],revNodalPtr+revNodalIndxPtr[nodeId+1],std::bind2nd(std::equal_to<int>(),-1))=eltId;
2922         }
2923     }
2924 }
2925
2926 void MEDCoupling1DGTUMesh::checkFullyDefined() const
2927 {
2928   if(!((const DataArrayInt *)_conn) || !((const DataArrayInt *)_conn_indx) || !((const DataArrayDouble *)_coords))
2929     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFullyDefined : part of this is not fully defined.");
2930 }
2931
2932 bool MEDCoupling1DGTUMesh::isEmptyMesh(const std::vector<int>& tinyInfo) const
2933 {
2934   throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::isEmptyMesh : not implemented yet !");
2935 }
2936
2937 void MEDCoupling1DGTUMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
2938 {
2939   int it,order;
2940   double time=getTime(it,order);
2941   tinyInfo.clear(); tinyInfoD.clear(); littleStrings.clear();
2942   //
2943   littleStrings.push_back(getName());
2944   littleStrings.push_back(getDescription());
2945   littleStrings.push_back(getTimeUnit());
2946   //
2947   std::vector<std::string> littleStrings2,littleStrings3,littleStrings4;
2948   if((const DataArrayDouble *)_coords)
2949     _coords->getTinySerializationStrInformation(littleStrings2);
2950   if((const DataArrayInt *)_conn)
2951     _conn->getTinySerializationStrInformation(littleStrings3);
2952   if((const DataArrayInt *)_conn_indx)
2953     _conn_indx->getTinySerializationStrInformation(littleStrings4);
2954   int sz0((int)littleStrings2.size()),sz1((int)littleStrings3.size()),sz2((int)littleStrings4.size());
2955   littleStrings.insert(littleStrings.end(),littleStrings2.begin(),littleStrings2.end());
2956   littleStrings.insert(littleStrings.end(),littleStrings3.begin(),littleStrings3.end());
2957   littleStrings.insert(littleStrings.end(),littleStrings4.begin(),littleStrings4.end());
2958   //
2959   tinyInfo.push_back(getCellModelEnum());
2960   tinyInfo.push_back(it);
2961   tinyInfo.push_back(order);
2962   std::vector<int> tinyInfo2,tinyInfo3,tinyInfo4;
2963   if((const DataArrayDouble *)_coords)
2964     _coords->getTinySerializationIntInformation(tinyInfo2);
2965   if((const DataArrayInt *)_conn)
2966     _conn->getTinySerializationIntInformation(tinyInfo3);
2967   if((const DataArrayInt *)_conn_indx)
2968     _conn_indx->getTinySerializationIntInformation(tinyInfo4);
2969   int sz3((int)tinyInfo2.size()),sz4((int)tinyInfo3.size()),sz5((int)tinyInfo4.size());
2970   tinyInfo.push_back(sz0); tinyInfo.push_back(sz1); tinyInfo.push_back(sz2); tinyInfo.push_back(sz3); tinyInfo.push_back(sz4);  tinyInfo.push_back(sz5);
2971   tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
2972   tinyInfo.insert(tinyInfo.end(),tinyInfo3.begin(),tinyInfo3.end());
2973   tinyInfo.insert(tinyInfo.end(),tinyInfo4.begin(),tinyInfo4.end());
2974   //
2975   tinyInfoD.push_back(time);
2976 }
2977
2978 void MEDCoupling1DGTUMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
2979 {
2980   std::vector<int> tinyInfo2(tinyInfo.begin()+9,tinyInfo.begin()+9+tinyInfo[6]);
2981   std::vector<int> tinyInfo1(tinyInfo.begin()+9+tinyInfo[6],tinyInfo.begin()+9+tinyInfo[6]+tinyInfo[7]);
2982   std::vector<int> tinyInfo12(tinyInfo.begin()+9+tinyInfo[6]+tinyInfo[7],tinyInfo.begin()+9+tinyInfo[6]+tinyInfo[7]+tinyInfo[8]);
2983   MCAuto<DataArrayInt> p1(DataArrayInt::New()); p1->resizeForUnserialization(tinyInfo1);
2984   MCAuto<DataArrayInt> p2(DataArrayInt::New()); p2->resizeForUnserialization(tinyInfo12);
2985   std::vector<const DataArrayInt *> v(2); v[0]=p1; v[1]=p2;
2986   p2=DataArrayInt::Aggregate(v);
2987   a2->resizeForUnserialization(tinyInfo2);
2988   a1->alloc(p2->getNbOfElems(),1);
2989 }
2990
2991 void MEDCoupling1DGTUMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
2992 {
2993   int sz(0);
2994   if((const DataArrayInt *)_conn)
2995     if(_conn->isAllocated())
2996       sz=_conn->getNbOfElems();
2997   if((const DataArrayInt *)_conn_indx)
2998     if(_conn_indx->isAllocated())
2999       sz+=_conn_indx->getNbOfElems();
3000   a1=DataArrayInt::New();
3001   a1->alloc(sz,1);
3002   int *work(a1->getPointer());
3003   if(sz!=0 && (const DataArrayInt *)_conn)
3004     work=std::copy(_conn->begin(),_conn->end(),a1->getPointer());
3005   if(sz!=0 && (const DataArrayInt *)_conn_indx)
3006     std::copy(_conn_indx->begin(),_conn_indx->end(),work);
3007   sz=0;
3008   if((const DataArrayDouble *)_coords)
3009     if(_coords->isAllocated())
3010       sz=_coords->getNbOfElems();
3011   a2=DataArrayDouble::New();
3012   a2->alloc(sz,1);
3013   if(sz!=0 && (const DataArrayDouble *)_coords)
3014     std::copy(_coords->begin(),_coords->end(),a2->getPointer());
3015 }
3016
3017 void MEDCoupling1DGTUMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
3018                                            const std::vector<std::string>& littleStrings)
3019 {
3020   INTERP_KERNEL::NormalizedCellType gt((INTERP_KERNEL::NormalizedCellType)tinyInfo[0]);
3021   _cm=&INTERP_KERNEL::CellModel::GetCellModel(gt);
3022   setName(littleStrings[0]);
3023   setDescription(littleStrings[1]);
3024   setTimeUnit(littleStrings[2]);
3025   setTime(tinyInfoD[0],tinyInfo[1],tinyInfo[2]);
3026   int sz0(tinyInfo[3]),sz1(tinyInfo[4]),sz2(tinyInfo[5]),sz3(tinyInfo[6]),sz4(tinyInfo[7]),sz5(tinyInfo[8]);
3027   //
3028   _coords=DataArrayDouble::New();
3029   std::vector<int> tinyInfo2(tinyInfo.begin()+9,tinyInfo.begin()+9+sz3);
3030   _coords->resizeForUnserialization(tinyInfo2);
3031   std::copy(a2->begin(),a2->end(),_coords->getPointer());
3032   _conn=DataArrayInt::New();
3033   std::vector<int> tinyInfo3(tinyInfo.begin()+9+sz3,tinyInfo.begin()+9+sz3+sz4);
3034   _conn->resizeForUnserialization(tinyInfo3);
3035   std::copy(a1->begin(),a1->begin()+_conn->getNbOfElems(),_conn->getPointer());
3036   _conn_indx=DataArrayInt::New();
3037   std::vector<int> tinyInfo4(tinyInfo.begin()+9+sz3+sz4,tinyInfo.begin()+9+sz3+sz4+sz5);
3038   _conn_indx->resizeForUnserialization(tinyInfo4);
3039   std::copy(a1->begin()+_conn->getNbOfElems(),a1->end(),_conn_indx->getPointer());
3040   std::vector<std::string> littleStrings2(littleStrings.begin()+3,littleStrings.begin()+3+sz0);
3041   _coords->finishUnserialization(tinyInfo2,littleStrings2);
3042   std::vector<std::string> littleStrings3(littleStrings.begin()+3+sz0,littleStrings.begin()+3+sz0+sz1);
3043   _conn->finishUnserialization(tinyInfo3,littleStrings3);
3044   std::vector<std::string> littleStrings4(littleStrings.begin()+3+sz0+sz1,littleStrings.begin()+3+sz0+sz1+sz2);
3045   _conn_indx->finishUnserialization(tinyInfo4,littleStrings4);
3046 }
3047
3048 /*!
3049  * Finds nodes not used in any cell and returns an array giving a new id to every node
3050  * by excluding the unused nodes, for which the array holds -1. The result array is
3051  * a mapping in "Old to New" mode.
3052  *  \param [out] nbrOfNodesInUse - number of node ids present in the nodal connectivity.
3053  *  \return DataArrayInt * - a new instance of DataArrayInt. Its length is \a
3054  *          this->getNumberOfNodes(). It holds for each node of \a this mesh either -1
3055  *          if the node is unused or a new id else. The caller is to delete this
3056  *          array using decrRef() as it is no more needed.
3057  *  \throw If the coordinates array is not set.
3058  *  \throw If the nodal connectivity of cells is not defined.
3059  *  \throw If the nodal connectivity includes an invalid id.
3060  *  \sa MEDCoupling1DGTUMesh::getNodeIdsInUse, areAllNodesFetched
3061  */
3062 DataArrayInt *MEDCoupling1DGTUMesh::computeFetchedNodeIds() const
3063 {
3064   checkConsistency();
3065   int nbNodes(getNumberOfNodes());
3066   std::vector<bool> fetchedNodes(nbNodes,false);
3067   computeNodeIdsAlg(fetchedNodes);
3068   int sz((int)std::count(fetchedNodes.begin(),fetchedNodes.end(),true));
3069   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(sz,1);
3070   int *retPtr(ret->getPointer());
3071   for(int i=0;i<nbNodes;i++)
3072     if(fetchedNodes[i])
3073       *retPtr++=i;
3074   return ret.retn();
3075 }
3076
3077 /*!
3078  * Finds nodes not used in any cell and returns an array giving a new id to every node
3079  * by excluding the unused nodes, for which the array holds -1. The result array is
3080  * a mapping in "Old to New" mode.
3081  *  \param [out] nbrOfNodesInUse - number of node ids present in the nodal connectivity.
3082  *  \return DataArrayInt * - a new instance of DataArrayInt. Its length is \a
3083  *          this->getNumberOfNodes(). It holds for each node of \a this mesh either -1
3084  *          if the node is unused or a new id else. The caller is to delete this
3085  *          array using decrRef() as it is no more needed.
3086  *  \throw If the coordinates array is not set.
3087  *  \throw If the nodal connectivity of cells is not defined.
3088  *  \throw If the nodal connectivity includes an invalid id.
3089  *  \sa MEDCoupling1DGTUMesh::computeFetchedNodeIds, areAllNodesFetched
3090  */
3091 DataArrayInt *MEDCoupling1DGTUMesh::getNodeIdsInUse(int& nbrOfNodesInUse) const
3092 {
3093   nbrOfNodesInUse=-1;
3094   int nbOfNodes=getNumberOfNodes();
3095   int nbOfCells=getNumberOfCells();//checkConsistencyLight
3096   MCAuto<DataArrayInt> ret=DataArrayInt::New();
3097   ret->alloc(nbOfNodes,1);
3098   int *traducer=ret->getPointer();
3099   std::fill(traducer,traducer+nbOfNodes,-1);
3100   const int *conn=_conn->begin(),*conni(_conn_indx->begin());
3101   for(int i=0;i<nbOfCells;i++,conni++)
3102     {
3103       int nbNodesPerCell=conni[1]-conni[0];
3104       for(int j=0;j<nbNodesPerCell;j++)
3105         {
3106           int nodeId=conn[conni[0]+j];
3107           if(nodeId==-1) continue;
3108           if(nodeId>=0 && nodeId<nbOfNodes)
3109             traducer[nodeId]=1;
3110           else
3111             {
3112               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getNodeIdsInUse : In cell #" << i  << " presence of node id " <<  nodeId << " not in [0," << nbOfNodes << ") !";
3113               throw INTERP_KERNEL::Exception(oss.str().c_str());
3114             }
3115         }
3116     }
3117   nbrOfNodesInUse=(int)std::count(traducer,traducer+nbOfNodes,1);
3118   std::transform(traducer,traducer+nbOfNodes,traducer,MEDCouplingAccVisit());
3119   return ret.retn();
3120 }
3121
3122 /*!
3123  * This method renumbers only nodal connectivity in \a this. The renumbering is only an offset applied. So this method is a specialization of
3124  * \a renumberNodesInConn. \b WARNING, this method does not check that the resulting node ids in the nodal connectivity is in a valid range !
3125  *
3126  * \param [in] offset - specifies the offset to be applied on each element of connectivity.
3127  *
3128  * \sa renumberNodesInConn
3129  */
3130 void MEDCoupling1DGTUMesh::renumberNodesWithOffsetInConn(int offset)
3131 {
3132   getNumberOfCells();//only to check that all is well defined.
3133   //
3134   int nbOfTuples(_conn->getNumberOfTuples());
3135   int *pt(_conn->getPointer());
3136   for(int i=0;i<nbOfTuples;i++,pt++)
3137     {
3138       if(*pt==-1) continue;
3139       *pt+=offset;
3140     }
3141   //
3142   updateTime();
3143 }
3144
3145 /*!
3146  *  Same than renumberNodesInConn(const int *) except that here the format of old-to-new traducer is using map instead
3147  *  of array. This method is dedicated for renumbering from a big set of nodes the a tiny set of nodes which is the case during extraction
3148  *  of a big mesh.
3149  */
3150 void MEDCoupling1DGTUMesh::renumberNodesInConn(const INTERP_KERNEL::HashMap<int,int>& newNodeNumbersO2N)
3151 {
3152   getNumberOfCells();//only to check that all is well defined.
3153   //
3154   int nbElemsIn(getNumberOfNodes()),nbOfTuples(_conn->getNumberOfTuples());
3155   int *pt(_conn->getPointer());
3156   for(int i=0;i<nbOfTuples;i++,pt++)
3157     {
3158       if(*pt==-1) continue;
3159       if(*pt>=0 && *pt<nbElemsIn)
3160         {
3161           INTERP_KERNEL::HashMap<int,int>::const_iterator it(newNodeNumbersO2N.find(*pt));
3162           if(it!=newNodeNumbersO2N.end())
3163             *pt=(*it).second;
3164           else
3165             {
3166               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::renumberNodesInConn : At pos #" << i << " of connectivity, node id is " << *pt << ". Not in keys of input map !";
3167               throw INTERP_KERNEL::Exception(oss.str().c_str());
3168             }
3169         }
3170       else
3171         {
3172           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::renumberNodesInConn : error on tuple #" << i << " value is " << *pt << " and indirectionnal array as a size equal to " << nbElemsIn;
3173           throw INTERP_KERNEL::Exception(oss.str().c_str());
3174         }
3175     }
3176   //
3177   updateTime();
3178 }
3179
3180 /*!
3181  * Changes ids of nodes within the nodal connectivity arrays according to a permutation
3182  * array in "Old to New" mode. The node coordinates array is \b not changed by this method.
3183  * This method is a generalization of shiftNodeNumbersInConn().
3184  *  \warning This method performs no check of validity of new ids. **Use it with care !**
3185  *  \param [in] newNodeNumbersO2N - a permutation array, of length \a
3186  *         this->getNumberOfNodes(), in "Old to New" mode.
3187  *         See \ref numbering for more info on renumbering modes.
3188  *  \throw If the nodal connectivity of cells is not defined.
3189  */
3190 void MEDCoupling1DGTUMesh::renumberNodesInConn(const int *newNodeNumbersO2N)
3191 {
3192   getNumberOfCells();//only to check that all is well defined.
3193   //
3194   int nbElemsIn(getNumberOfNodes()),nbOfTuples(_conn->getNumberOfTuples());
3195   int *pt(_conn->getPointer());
3196   for(int i=0;i<nbOfTuples;i++,pt++)
3197     {
3198       if(*pt==-1) continue;
3199       if(*pt>=0 && *pt<nbElemsIn)
3200         *pt=newNodeNumbersO2N[*pt];
3201       else
3202         {
3203           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::renumberNodesInConn : error on tuple #" << i << " value is " << *pt << " and indirectionnal array as a size equal to " << nbElemsIn;
3204           throw INTERP_KERNEL::Exception(oss.str().c_str());
3205         }
3206     }
3207   //
3208   updateTime();
3209 }
3210
3211 /*!
3212  * Keeps from \a this only cells which constituing point id are in the ids specified by [\a begin,\a end).
3213  * The resulting cell ids are stored at the end of the 'cellIdsKept' parameter.
3214  * Parameter \a fullyIn specifies if a cell that has part of its nodes in ids array is kept or not.
3215  * If \a fullyIn is true only cells whose ids are \b fully contained in [\a begin,\a end) tab will be kept.
3216  *
3217  * \param [in] begin input start of array of node ids.
3218  * \param [in] end input end of array of node ids.
3219  * \param [in] fullyIn input that specifies if all node ids must be in [\a begin,\a end) array to consider cell to be in.
3220  * \param [in,out] cellIdsKeptArr array where all candidate cell ids are put at the end.
3221  */
3222 void MEDCoupling1DGTUMesh::fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const
3223 {
3224   int nbOfCells=getNumberOfCells();
3225   MCAuto<DataArrayInt> cellIdsKept=DataArrayInt::New(); cellIdsKept->alloc(0,1);
3226   int tmp=-1;
3227   int sz=_conn->getMaxValue(tmp); sz=std::max(sz,0)+1;
3228   std::vector<bool> fastFinder(sz,false);
3229   for(const int *work=begin;work!=end;work++)
3230     if(*work>=0 && *work<sz)
3231       fastFinder[*work]=true;
3232   const int *conn=_conn->begin(),*conni=_conn_indx->begin();
3233   for(int i=0;i<nbOfCells;i++,conni++)
3234     {
3235       int ref=0,nbOfHit=0;
3236       int nbNodesPerCell=conni[1]-conni[0];
3237       if(nbNodesPerCell>=0)
3238         {
3239           for(int j=0;j<nbNodesPerCell;j++)
3240             {
3241               int nodeId=conn[conni[0]+j];
3242               if(nodeId>=0)
3243                 {
3244                   ref++;
3245                   if(fastFinder[nodeId])
3246                     nbOfHit++;
3247                 }
3248             }
3249         }
3250       else
3251         {
3252           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::fillCellIdsToKeepFromNodeIds : invalid index array for cell #" << i << " !";
3253           throw INTERP_KERNEL::Exception(oss.str().c_str());
3254         }
3255       if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
3256         cellIdsKept->pushBackSilent(i);
3257     }
3258   cellIdsKeptArr=cellIdsKept.retn();
3259 }
3260
3261 void MEDCoupling1DGTUMesh::allocateCells(int nbOfCells)
3262 {
3263   if(nbOfCells<0)
3264     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::allocateCells : the input number of cells should be >= 0 !");
3265   _conn=DataArrayInt::New();
3266   _conn->reserve(nbOfCells*3);
3267   _conn_indx=DataArrayInt::New();
3268   _conn_indx->reserve(nbOfCells+1); _conn_indx->pushBackSilent(0);
3269   declareAsNew();
3270 }
3271
3272 /*!
3273  * Appends at the end of \a this a cell having nodal connectivity array defined in [ \a nodalConnOfCellBg, \a nodalConnOfCellEnd ).
3274  *
3275  * \param [in] nodalConnOfCellBg - the begin (included) of nodal connectivity of the cell to add.
3276  * \param [in] nodalConnOfCellEnd - the end (excluded) of nodal connectivity of the cell to add.
3277  * \throw If the length of the input nodal connectivity array of the cell to add is not equal to number of nodes per cell relative to the unique geometric type
3278  *        attached to \a this.
3279  * \thow If the nodal connectivity array in \a this is null (call MEDCoupling1SGTUMesh::allocateCells before).
3280  */
3281 void MEDCoupling1DGTUMesh::insertNextCell(const int *nodalConnOfCellBg, const int *nodalConnOfCellEnd)
3282 {
3283   int sz=(int)std::distance(nodalConnOfCellBg,nodalConnOfCellEnd);
3284   DataArrayInt *c(_conn),*c2(_conn_indx);
3285   if(c && c2)
3286     {
3287       int pos=c2->back();
3288       if(pos==c->getNumberOfTuples())
3289         {
3290           c->pushBackValsSilent(nodalConnOfCellBg,nodalConnOfCellEnd);
3291           c2->pushBackSilent(pos+sz);
3292         }
3293       else
3294         {
3295           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::insertNextCell : The nodal index array (end=" << pos << ") mismatches with nodal array (length=" << c->getNumberOfTuples() << ") !";
3296           throw INTERP_KERNEL::Exception(oss.str().c_str());
3297         }
3298     }
3299   else
3300     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::insertNextCell : nodal connectivity array is null ! Call MEDCoupling1DGTUMesh::allocateCells before !");
3301 }
3302
3303 void MEDCoupling1DGTUMesh::setNodalConnectivity(DataArrayInt *nodalConn, DataArrayInt *nodalConnIndex)
3304 {
3305   if(nodalConn)
3306     nodalConn->incrRef();
3307   _conn=nodalConn;
3308   if(nodalConnIndex)
3309     nodalConnIndex->incrRef();
3310   _conn_indx=nodalConnIndex;
3311   declareAsNew();
3312 }
3313
3314 /*!
3315  * \return DataArrayInt * - the internal reference to the nodal connectivity. The caller is not reponsible to deallocate it.
3316  */
3317 DataArrayInt *MEDCoupling1DGTUMesh::getNodalConnectivity() const
3318 {
3319   const DataArrayInt *ret(_conn);
3320   return const_cast<DataArrayInt *>(ret);
3321 }
3322
3323 /*!
3324  * \return DataArrayInt * - the internal reference to the nodal connectivity index. The caller is not reponsible to deallocate it.
3325  */
3326 DataArrayInt *MEDCoupling1DGTUMesh::getNodalConnectivityIndex() const
3327 {
3328   const DataArrayInt *ret(_conn_indx);
3329   return const_cast<DataArrayInt *>(ret);
3330 }
3331
3332 /*!
3333  * See the definition of the nodal connectivity pack \ref MEDCoupling1DGTUMesh::isPacked "here".
3334  * This method tries to build a new instance geometrically equivalent to \a this, by limiting at most the number of new object (nodal connectivity).
3335  * Geometrically the returned mesh is equal to \a this. So if \a this is already packed, the return value is a shallow copy of \a this.
3336  *
3337  * Whatever the status of pack of \a this, the coordinates array of the returned newly created instance is the same than those in \a this.
3338  *
3339  * \param [out] isShallowCpyOfNodalConnn - tells if the returned instance share the same pair of nodal connectivity arrays (true) or if nodal
3340  *              connectivity arrays are different (false)
3341  * \return a new object to be managed by the caller.
3342  *
3343  * \sa MEDCoupling1DGTUMesh::retrievePackedNodalConnectivity, MEDCoupling1DGTUMesh::isPacked
3344  */
3345 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::copyWithNodalConnectivityPacked(bool& isShallowCpyOfNodalConnn) const
3346 {
3347   MCAuto<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
3348   DataArrayInt *nc=0,*nci=0;
3349   isShallowCpyOfNodalConnn=retrievePackedNodalConnectivity(nc,nci);
3350   MCAuto<DataArrayInt> ncs(nc),ncis(nci);
3351   ret->_conn=ncs; ret->_conn_indx=ncis;
3352   ret->setCoords(getCoords());
3353   return ret.retn();
3354 }
3355
3356 /*!
3357  * This method allows to compute, if needed, the packed nodal connectivity pair.
3358  * Indeed, it is possible to store in \a this a nodal connectivity array bigger than ranges convered by nodal connectivity index array.
3359  * It is typically the case when nodalConnIndx starts with an id greater than 0, and finishes with id less than number of tuples in \c this->_conn.
3360  *
3361  * If \a this looks packed (the front of nodal connectivity index equal to 0 and back of connectivity index equal to number of tuple of nodal connectivity array)
3362  * true will be returned and respectively \a this->_conn and \a this->_conn_indx (with ref counter incremented). This is the classical case.
3363  *
3364  * If nodal connectivity index points to a subpart of nodal connectivity index the packed pair of arrays will be computed (new objects) and returned and false
3365  * will be returned.
3366  *
3367  * This method return 3 elements.
3368  * \param [out] nodalConn - a pointer that can be equal to \a this->_conn if true is returned (general case). Whatever the value of return parameter
3369  *                          this pointer can be seen as a new object, that is to managed by the caller.
3370  * \param [out] nodalConnIndx - a pointer that can be equal to \a this->_conn_indx if true is returned (general case). Whatever the value of return parameter
3371  *                              this pointer can be seen as a new object, that is to managed by the caller.
3372  * \return bool - an indication of the content of the 2 output parameters. If true, \a this looks packed (general case), if true, \a this is not packed then
3373  * output parameters are newly created objects.
3374  *
3375  * \throw if \a this does not pass MEDCoupling1DGTUMesh::checkConsistencyLight test
3376  */
3377 bool MEDCoupling1DGTUMesh::retrievePackedNodalConnectivity(DataArrayInt *&nodalConn, DataArrayInt *&nodalConnIndx) const
3378 {
3379   if(isPacked())//performs the checkConsistencyLight
3380     {
3381       const DataArrayInt *c0(_conn),*c1(_conn_indx);
3382       nodalConn=const_cast<DataArrayInt *>(c0); nodalConnIndx=const_cast<DataArrayInt *>(c1);
3383       nodalConn->incrRef(); nodalConnIndx->incrRef();
3384       return true;
3385     }
3386   int bg=_conn_indx->front(),end=_conn_indx->back();
3387   MCAuto<DataArrayInt> nc(_conn->selectByTupleIdSafeSlice(bg,end,1));
3388   MCAuto<DataArrayInt> nci(_conn_indx->deepCopy());
3389   nci->applyLin(1,-bg);
3390   nodalConn=nc.retn(); nodalConnIndx=nci.retn();
3391   return false;
3392 }
3393
3394 /*
3395  * If \a this looks packed (the front of nodal connectivity index equal to 0 and back of connectivity index equal to number of tuple of nodal connectivity array)
3396  * true will be returned and respectively \a this->_conn and \a this->_conn_indx (with ref counter incremented). This is the classical case.
3397  * If nodal connectivity index points to a subpart of nodal connectivity index false will be returned.
3398  * \return bool - true if \a this looks packed, false is not.
3399  *
3400  * \throw if \a this does not pass MEDCoupling1DGTUMesh::checkConsistencyLight test
3401  */
3402 bool MEDCoupling1DGTUMesh::isPacked() const
3403 {
3404   checkConsistencyLight();
3405   return _conn_indx->front()==0 && _conn_indx->back()==_conn->getNumberOfTuples();
3406 }
3407
3408 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshes(const MEDCoupling1DGTUMesh *mesh1, const MEDCoupling1DGTUMesh *mesh2)
3409 {
3410   std::vector<const MEDCoupling1DGTUMesh *> tmp(2);
3411   tmp[0]=const_cast<MEDCoupling1DGTUMesh *>(mesh1); tmp[1]=const_cast<MEDCoupling1DGTUMesh *>(mesh2);
3412   return Merge1DGTUMeshes(tmp);
3413 }
3414
3415 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshes(std::vector<const MEDCoupling1DGTUMesh *>& a)
3416 {
3417   std::size_t sz=a.size();
3418   if(sz==0)
3419     return Merge1DGTUMeshesLL(a);
3420   for(std::size_t ii=0;ii<sz;ii++)
3421     if(!a[ii])
3422       {
3423         std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::Merge1DGTUMeshes : item #" << ii << " in input array of size "<< sz << " is empty !";
3424         throw INTERP_KERNEL::Exception(oss.str().c_str());
3425       }
3426   const INTERP_KERNEL::CellModel *cm=&(a[0]->getCellModel());
3427   for(std::size_t ii=0;ii<sz;ii++)
3428     if(&(a[ii]->getCellModel())!=cm)
3429       throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshes : all items must have the same geo type !");
3430   std::vector< MCAuto<MEDCoupling1DGTUMesh> > bb(sz);
3431   std::vector< const MEDCoupling1DGTUMesh * > aa(sz);
3432   int spaceDim=-3;
3433   for(std::size_t i=0;i<sz && spaceDim==-3;i++)
3434     {
3435       const MEDCoupling1DGTUMesh *cur=a[i];
3436       const DataArrayDouble *coo=cur->getCoords();
3437       if(coo)
3438         spaceDim=coo->getNumberOfComponents();
3439     }
3440   if(spaceDim==-3)
3441     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshes : no spaceDim specified ! unable to perform merge !");
3442   for(std::size_t i=0;i<sz;i++)
3443     {
3444       bb[i]=a[i]->buildSetInstanceFromThis(spaceDim);
3445       aa[i]=bb[i];
3446     }
3447   return Merge1DGTUMeshesLL(aa);
3448 }
3449
3450 /*!
3451  * \throw If presence of a null instance in the input vector \a a.
3452  * \throw If a is empty
3453  */
3454 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords(std::vector<const MEDCoupling1DGTUMesh *>& a)
3455 {
3456   if(a.empty())
3457     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : input array must be NON EMPTY !");
3458   std::vector<const MEDCoupling1DGTUMesh *>::const_iterator it=a.begin();
3459   if(!(*it))
3460     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : null instance in the first element of input vector !");
3461   std::vector< MCAuto<MEDCoupling1DGTUMesh> > objs(a.size());
3462   std::vector<const DataArrayInt *> ncs(a.size()),ncis(a.size());
3463   (*it)->getNumberOfCells();//to check that all is OK
3464   const DataArrayDouble *coords=(*it)->getCoords();
3465   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
3466   bool tmp;
3467   objs[0]=(*it)->copyWithNodalConnectivityPacked(tmp);
3468   ncs[0]=objs[0]->getNodalConnectivity(); ncis[0]=objs[0]->getNodalConnectivityIndex();
3469   it++;
3470   for(int i=1;it!=a.end();i++,it++)
3471     {
3472       if(!(*it))
3473         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : presence of null instance !");
3474       if(cm!=&((*it)->getCellModel()))
3475         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1DGTUMeshes impossible !");
3476       (*it)->getNumberOfCells();//to check that all is OK
3477       objs[i]=(*it)->copyWithNodalConnectivityPacked(tmp);
3478       ncs[i]=objs[i]->getNodalConnectivity(); ncis[i]=objs[i]->getNodalConnectivityIndex();
3479       if(coords!=(*it)->getCoords())
3480         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : not lying on same coords !");
3481     }
3482   MCAuto<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh("merge",*cm));
3483   ret->setCoords(coords);
3484   ret->_conn=DataArrayInt::Aggregate(ncs);
3485   ret->_conn_indx=DataArrayInt::AggregateIndexes(ncis);
3486   return ret.retn();
3487 }
3488
3489 /*!
3490  * Assume that all instances in \a a are non null. If null it leads to a crash. That's why this method is assigned to be low level (LL)
3491  */
3492 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshesLL(std::vector<const MEDCoupling1DGTUMesh *>& a)
3493 {
3494   if(a.empty())
3495     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshes : input array must be NON EMPTY !");
3496   std::vector< MCAuto<MEDCoupling1DGTUMesh> > objs(a.size());
3497   std::vector<const DataArrayInt *> ncs(a.size()),ncis(a.size());
3498   std::vector<const MEDCoupling1DGTUMesh *>::const_iterator it=a.begin();
3499   std::vector<int> nbNodesPerElt(a.size());
3500   int nbOfCells=(*it)->getNumberOfCells();
3501   bool tmp;
3502   objs[0]=(*it)->copyWithNodalConnectivityPacked(tmp);
3503   ncs[0]=objs[0]->getNodalConnectivity(); ncis[0]=objs[0]->getNodalConnectivityIndex();
3504   nbNodesPerElt[0]=0;
3505   int prevNbOfNodes=(*it)->getNumberOfNodes();
3506   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
3507   it++;
3508   for(int i=1;it!=a.end();i++,it++)
3509     {
3510       if(cm!=&((*it)->getCellModel()))
3511         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1DGTUMeshes impossible !");
3512       objs[i]=(*it)->copyWithNodalConnectivityPacked(tmp);
3513       ncs[i]=objs[i]->getNodalConnectivity(); ncis[i]=objs[i]->getNodalConnectivityIndex();
3514       nbOfCells+=(*it)->getNumberOfCells();
3515       nbNodesPerElt[i]=nbNodesPerElt[i-1]+prevNbOfNodes;
3516       prevNbOfNodes=(*it)->getNumberOfNodes();
3517     }
3518   std::vector<const MEDCouplingPointSet *> aps(a.size());
3519   std::copy(a.begin(),a.end(),aps.begin());
3520   MCAuto<DataArrayDouble> pts=MergeNodesArray(aps);
3521   MCAuto<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh("merge",*cm));
3522   ret->setCoords(pts);
3523   ret->_conn=AggregateNodalConnAndShiftNodeIds(ncs,nbNodesPerElt);
3524   ret->_conn_indx=DataArrayInt::AggregateIndexes(ncis);
3525   return ret.retn();
3526 }
3527
3528 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::buildSetInstanceFromThis(int spaceDim) const
3529 {
3530   MCAuto<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
3531   MCAuto<DataArrayInt> tmp1,tmp2;
3532   const DataArrayInt *nodalConn(_conn),*nodalConnI(_conn_indx);
3533   if(!nodalConn)
3534     {
3535       tmp1=DataArrayInt::New(); tmp1->alloc(0,1);
3536     }
3537   else
3538     tmp1=_conn;
3539   ret->_conn=tmp1;
3540   //
3541   if(!nodalConnI)
3542     {
3543       tmp2=DataArrayInt::New(); tmp2->alloc(1,1); tmp2->setIJ(0,0,0);
3544     }
3545   else
3546     tmp2=_conn_indx;
3547   ret->_conn_indx=tmp2;
3548   //
3549   if(!_coords)
3550     {
3551       MCAuto<DataArrayDouble> coords=DataArrayDouble::New(); coords->alloc(0,spaceDim);
3552       ret->setCoords(coords);
3553     }
3554   else
3555     ret->setCoords(_coords);
3556   return ret.retn();
3557 }
3558
3559 /*!
3560  * This method aggregate the bbox of each cell and put it into bbox parameter.
3561  *
3562  * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
3563  *                         For all other cases this input parameter is ignored.
3564  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
3565  *
3566  * \throw If \a this is not fully set (coordinates and connectivity).
3567  * \throw If a cell in \a this has no valid nodeId.
3568  */
3569 DataArrayDouble *MEDCoupling1DGTUMesh::getBoundingBoxForBBTree(double arcDetEps) const
3570 {
3571   checkFullyDefined();
3572   int spaceDim(getSpaceDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
3573   MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
3574   double *bbox(ret->getPointer());
3575   for(int i=0;i<nbOfCells*spaceDim;i++)
3576     {
3577       bbox[2*i]=std::numeric_limits<double>::max();
3578       bbox[2*i+1]=-std::numeric_limits<double>::max();
3579     }
3580   const double *coordsPtr(_coords->getConstPointer());
3581   const int *conn(_conn->getConstPointer()),*connI(_conn_indx->getConstPointer());
3582   for(int i=0;i<nbOfCells;i++)
3583     {
3584       int offset=connI[i];
3585       int nbOfNodesForCell(connI[i+1]-offset),kk(0);
3586       for(int j=0;j<nbOfNodesForCell;j++)
3587         {
3588           int nodeId=conn[offset+j];
3589           if(nodeId>=0 && nodeId<nbOfNodes)
3590             {
3591               for(int k=0;k<spaceDim;k++)
3592                 {
3593                   bbox[2*spaceDim*i+2*k]=std::min(bbox[2*spaceDim*i+2*k],coordsPtr[spaceDim*nodeId+k]);
3594                   bbox[2*spaceDim*i+2*k+1]=std::max(bbox[2*spaceDim*i+2*k+1],coordsPtr[spaceDim*nodeId+k]);
3595                 }
3596               kk++;
3597             }
3598         }
3599       if(kk==0)
3600         {
3601           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getBoundingBoxForBBTree : cell #" << i << " contains no valid nodeId !";
3602           throw INTERP_KERNEL::Exception(oss.str().c_str());
3603         }
3604     }
3605   return ret.retn();
3606 }
3607
3608 /*!
3609  * Returns the cell field giving for each cell in \a this its diameter. Diameter means the max length of all possible SEG2 in the cell.
3610  *
3611  * \return a new instance of field containing the result. The returned instance has to be deallocated by the caller.
3612  */
3613 MEDCouplingFieldDouble *MEDCoupling1DGTUMesh::computeDiameterField() const
3614 {
3615   throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::computeDiameterField : not implemented yet for dynamic types !");
3616 }
3617
3618 std::vector<int> MEDCoupling1DGTUMesh::BuildAPolygonFromParts(const std::vector< std::vector<int> >& parts)
3619 {
3620   std::vector<int> ret;
3621   if(parts.empty())
3622     return ret;
3623   ret.insert(ret.end(),parts[0].begin(),parts[0].end());
3624   int ref(ret.back());
3625   std::size_t sz(parts.size()),nbh(1);
3626   std::vector<bool> b(sz,true); b[0]=false;
3627   while(nbh<sz)
3628     {
3629       std::size_t i(0);
3630       for(;i<sz;i++) if(b[i] && parts[i].front()==ref) { ret.insert(ret.end(),parts[i].begin()+1,parts[i].end()); nbh++; break; }
3631       if(i<sz)
3632         ref=ret.back();
3633       else
3634         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::BuildAPolygonFromParts : the input vector is not a part of a single polygon !");
3635     }
3636   if(ret.back()==ret.front())
3637     ret.pop_back();
3638   return ret;
3639 }
3640
3641 /*!
3642  * This method performs an aggregation of \a nodalConns (as DataArrayInt::Aggregate does) but in addition of that a shift is applied on the
3643  * values contained in \a nodalConns using corresponding offset specified in input \a offsetInNodeIdsPerElt.
3644  * But it also manage the values -1, that have a semantic in MEDCoupling1DGTUMesh class (separator for polyhedron).
3645  *
3646  * \param [in] nodalConns - a list of nodal connectivity arrays same size than \a offsetInNodeIdsPerElt.
3647  * \param [in] offsetInNodeIdsPerElt - a list of offsets to apply.
3648  * \return DataArrayInt * - A new object (to be managed by the caller) that is the result of the aggregation.
3649  * \throw If \a nodalConns or \a offsetInNodeIdsPerElt are empty.
3650  * \throw If \a nodalConns and \a offsetInNodeIdsPerElt have not the same size.
3651  * \throw If presence of null pointer in \a nodalConns.
3652  * \throw If presence of not allocated or array with not exactly one component in \a nodalConns.
3653  */
3654 DataArrayInt *MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds(const std::vector<const DataArrayInt *>& nodalConns, const std::vector<int>& offsetInNodeIdsPerElt)
3655 {
3656   std::size_t sz1(nodalConns.size()),sz2(offsetInNodeIdsPerElt.size());
3657   if(sz1!=sz2)
3658     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : input vectors do not have the same size !");
3659   if(sz1==0)
3660     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : empty vectors in input !");
3661   int nbOfTuples=0;
3662   for(std::vector<const DataArrayInt *>::const_iterator it=nodalConns.begin();it!=nodalConns.end();it++)
3663     {
3664       if(!(*it))
3665         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : presence of null pointer in input vector !");
3666       if(!(*it)->isAllocated())
3667         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : presence of non allocated array in input vector !");
3668       if((*it)->getNumberOfComponents()!=1)
3669         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : presence of array with not exactly one component !");
3670       nbOfTuples+=(*it)->getNumberOfTuples();
3671     }
3672   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(nbOfTuples,1);
3673   int *pt=ret->getPointer();
3674   int i=0;
3675   for(std::vector<const DataArrayInt *>::const_iterator it=nodalConns.begin();it!=nodalConns.end();it++,i++)
3676     {
3677       int curNbt=(*it)->getNumberOfTuples();
3678       const int *inPt=(*it)->begin();
3679       int offset=offsetInNodeIdsPerElt[i];
3680       for(int j=0;j<curNbt;j++,pt++)
3681         {
3682           if(inPt[j]!=-1)
3683             *pt=inPt[j]+offset;
3684           else
3685             *pt=-1;
3686         }
3687     }
3688   return ret.retn();
3689 }
3690
3691 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::New(const MEDCouplingUMesh *m)
3692 {
3693   if(!m)
3694     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::New : input mesh is null !");
3695   std::set<INTERP_KERNEL::NormalizedCellType> gts(m->getAllGeoTypes());
3696   if(gts.size()!=1)
3697     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::New : input mesh must have exactly one geometric type !");
3698   int geoType((int)*gts.begin());
3699   MCAuto<MEDCoupling1DGTUMesh> ret(MEDCoupling1DGTUMesh::New(m->getName(),*gts.begin()));
3700   ret->setCoords(m->getCoords()); ret->setDescription(m->getDescription());
3701   int nbCells(m->getNumberOfCells());
3702   MCAuto<DataArrayInt> conn(DataArrayInt::New()),connI(DataArrayInt::New());
3703   conn->alloc(m->getNodalConnectivityArrayLen()-nbCells,1); connI->alloc(nbCells+1,1);
3704   int *c(conn->getPointer()),*ci(connI->getPointer()); *ci=0;
3705   const int *cin(m->getNodalConnectivity()->begin()),*ciin(m->getNodalConnectivityIndex()->begin());
3706   for(int i=0;i<nbCells;i++,ciin++,ci++)
3707     {
3708       if(cin[ciin[0]]==geoType)
3709         {
3710           if(ciin[1]-ciin[0]>=1)
3711             {
3712               c=std::copy(cin+ciin[0]+1,cin+ciin[1],c);
3713               ci[1]=ci[0]+ciin[1]-ciin[0]-1;
3714             }
3715           else
3716             {
3717               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::New(const MEDCouplingUMesh *m) : something is wrong in the input mesh at cell #" << i << " ! The size of cell is not >=0 !";
3718               throw INTERP_KERNEL::Exception(oss.str().c_str());
3719             }
3720         }
3721       else
3722         {
3723           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::New(const MEDCouplingUMesh *m) : something is wrong in the input mesh at cell #" << i << " ! The geometric type is not those expected !";
3724           throw INTERP_KERNEL::Exception(oss.str().c_str());
3725         }
3726     }
3727   ret->setNodalConnectivity(conn,connI);
3728   return ret.retn();
3729 }