src/MEDCoupling/MEDCoupling1GTUMesh.cxx

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