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 MEDCoupling1GTUMesh::MEDCoupling1GTUMesh(const char *name, const INTERP_KERNEL::CellModel& cm):_cm(&cm)
  30 {
  31   setName(name);
  32 }
  33
  34 MEDCoupling1GTUMesh::MEDCoupling1GTUMesh(const MEDCoupling1GTUMesh& other, bool recDeepCpy):MEDCouplingPointSet(other,recDeepCpy),_cm(other._cm)
  35 {
  36 }
  37
  38 MEDCoupling1GTUMesh *MEDCoupling1GTUMesh::New(const char *name, INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
  39 {
  40   if(type==INTERP_KERNEL::NORM_ERROR)
  41     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::New : NORM_ERROR is not a valid type to be used as base geometric type for a mesh !");
  42   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
  43   if(!cm.isDynamic())
  44     return MEDCoupling1SGTUMesh::New(name,type);
  45   else
  46     return MEDCoupling1DGTUMesh::New(name,type);
  47 }
  48
  49 const INTERP_KERNEL::CellModel& MEDCoupling1GTUMesh::getCellModel() const throw(INTERP_KERNEL::Exception)
  50 {
  51   return *_cm;
  52 }
  53
  54 INTERP_KERNEL::NormalizedCellType MEDCoupling1GTUMesh::getCellModelEnum() const throw(INTERP_KERNEL::Exception)
  55 {
  56   return _cm->getEnum();
  57 }
  58
  59 int MEDCoupling1GTUMesh::getMeshDimension() const
  60 {
  61   return (int)_cm->getDimension();
  62 }
  63
  64 /*!
  65  * This method returns a newly allocated array containing cell ids (ascendingly sorted) whose geometric type are equal to type.
  66  * This method does not throw exception if geometric type \a type is not in \a this.
  67  * This method throws an INTERP_KERNEL::Exception if meshdimension of \b this is not equal to those of \b type.
  68  * The coordinates array is not considered here.
  69  *
  70  * \param [in] type the geometric type
  71  * \return cell ids in this having geometric type \a type.
  72  */
  73 DataArrayInt *MEDCoupling1GTUMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
  74 {
  75   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
  76   if(type==getCellModelEnum())
  77     ret->alloc(getNumberOfCells(),1);
  78   else
  79     ret->alloc(0,1);
  80   ret->iota();
  81   return ret.retn();
  82 }
  83
  84 /*!
  85  * Returns nb of cells having the geometric type \a type. No throw if no cells in \a this has the geometric type \a type.
  86  */
  87 int MEDCoupling1GTUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
  88 {
  89   return type==getCellModelEnum()?getNumberOfCells():0;
  90 }
  91
  92 /*!
  93  * Returns a type of a cell by its id.
  94  *  \param [in] cellId - the id of the cell of interest.
  95  *  \return INTERP_KERNEL::NormalizedCellType - enumeration item describing the cell type.
  96  *  \throw If \a cellId is invalid. Valid range is [0, \a this->getNumberOfCells() ).
  97  */
  98 INTERP_KERNEL::NormalizedCellType MEDCoupling1GTUMesh::getTypeOfCell(int cellId) const
  99 {
 100   if(cellId>=0 && cellId<getNumberOfCells())
 101     return getCellModelEnum();
 102   std::ostringstream oss; oss << "MEDCoupling1GTUMesh::getTypeOfCell : Requesting type of cell #" << cellId << " but it should be in [0," << getNumberOfCells() << ") !";
 103   throw INTERP_KERNEL::Exception(oss.str().c_str());
 104 }
 105
 106 /*!
 107  * Returns a set of all cell types available in \a this mesh.
 108  * \return std::set<INTERP_KERNEL::NormalizedCellType> - the set of cell types.
 109  * \warning this method does not throw any exception even if \a this is not defined.
 110  */
 111 std::set<INTERP_KERNEL::NormalizedCellType> MEDCoupling1GTUMesh::getAllGeoTypes() const
 112 {
 113   std::set<INTERP_KERNEL::NormalizedCellType> ret;
 114   ret.insert(getCellModelEnum());
 115   return ret;
 116 }
 117
 118 /*!
 119  * This method expects that \a this is sorted by types. If not an exception will be thrown.
 120  * This method returns in the same format as code (see MEDCouplingUMesh::checkTypeConsistencyAndContig or MEDCouplingUMesh::splitProfilePerType) how
 121  * \a this is composed in cell types.
 122  * The returned array is of size 3*n where n is the number of different types present in \a this.
 123  * For every k in [0,n] ret[3*k+2]==0 because it has no sense here.
 124  * This parameter is kept only for compatibility with other methode listed above.
 125  */
 126 std::vector<int> MEDCoupling1GTUMesh::getDistributionOfTypes() const throw(INTERP_KERNEL::Exception)
 127 {
 128   std::vector<int> ret(3);
 129   ret[0]=(int)getCellModelEnum(); ret[1]=getNumberOfCells(); ret[2]=0;
 130   return ret;
 131 }
 132
 133 /*!
 134  * 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.
 135  * 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.
 136  * This method has 1 input \a profile and 3 outputs \a code \a idsInPflPerType and \a idsPerType.
 137  *
 138  * \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.
 139  * \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,
 140  *              \a idsInPflPerType[i] stores the tuple ids in \a profile that correspond to the geometric type code[3*i+0]
 141  * \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.
 142  *              This vector can be empty in case of all geometric type cells are fully covered in ascending in the given input \a profile.
 143  *
 144  * \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.
 145  *
 146  * \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
 147  *
 148  *  \b Example1: <br>
 149  *          - Before \a this has 3 cells \a profile contains [0,1,2]
 150  *          - After \a code contains [NORM_...,nbCells,-1], \a idsInPflPerType [[0,1,2]] and \a idsPerType is empty <br>
 151  *
 152  *  \b Example2: <br>
 153  *          - Before \a this has 3 cells \a profile contains [1,2]
 154  *          - After \a code contains [NORM_...,nbCells,0], \a idsInPflPerType [[0,1]] and \a idsPerType is [[1,2]] <br>
 155
 156  */
 157 void MEDCoupling1GTUMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 158 {
 159   if(!profile)
 160     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::splitProfilePerType : input profile is NULL !");
 161   if(profile->getNumberOfComponents()!=1)
 162     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::splitProfilePerType : input profile should have exactly one component !");
 163   int nbTuples=profile->getNumberOfTuples();
 164   int nbOfCells=getNumberOfCells();
 165   code.resize(3); idsInPflPerType.resize(1);
 166   code[0]=(int)getCellModelEnum(); code[1]=nbTuples;
 167   idsInPflPerType.resize(1);
 168   if(profile->isIdentity() && nbTuples==nbOfCells)
 169     {
 170       code[2]=-1;
 171       idsInPflPerType[0]=const_cast<DataArrayInt *>(profile); idsInPflPerType[0]->incrRef();
 172       idsPerType.clear();
 173       return ;
 174     }
 175   code[2]=0;
 176   profile->checkAllIdsInRange(0,nbOfCells);
 177   idsPerType.resize(1);
 178   idsPerType[0]=const_cast<DataArrayInt *>(profile); idsPerType[0]->incrRef();
 179   idsInPflPerType[0]=DataArrayInt::Range(0,nbTuples,1);
 180 }
 181
 182 /*!
 183  * This method tries to minimize at most the number of deep copy.
 184  * So if \a idsPerType is not empty it can be returned directly (without copy, but with ref count incremented) in return.
 185  *
 186  * \sa MEDCouplingUMesh::checkTypeConsistencyAndContig
 187  */
 188 DataArrayInt *MEDCoupling1GTUMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 189 {
 190   int nbOfCells=getNumberOfCells();
 191   if(code.size()!=3)
 192     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : invalid input code should be exactly of size 3 !");
 193   if(code[0]!=(int)getCellModelEnum())
 194     {
 195       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() << ") !";
 196       throw INTERP_KERNEL::Exception(oss.str().c_str());
 197     }
 198   if(code[2]==-1)
 199     {
 200       if(code[1]==nbOfCells)
 201         return 0;
 202       else
 203         {
 204           std::ostringstream oss; oss << "MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : mismatch between the number of cells in this (" << nbOfCells << ") and the number of non profile (" << code[1] << ") !";
 205           throw INTERP_KERNEL::Exception(oss.str().c_str());
 206         }
 207     }
 208   if(code[2]!=0)
 209     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : single geo type mesh ! 0 or -1 is expected at pos #2 of input code !");
 210   if(idsPerType.size()!=1)
 211     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : input code points to DataArrayInt #0 whereas the size of idsPerType is not equal to 1 !");
 212   const DataArrayInt *pfl=idsPerType[0];
 213   if(!pfl)
 214     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : the input code points to a NULL DataArrayInt at rank 0 !");
 215   if(pfl->getNumberOfComponents()!=1)
 216     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::checkTypeConsistencyAndContig : input profile should have exactly one component !");
 217   pfl->checkAllIdsInRange(0,nbOfCells);
 218   pfl->incrRef();
 219   return const_cast<DataArrayInt *>(pfl);
 220 }
 221
 222 void MEDCoupling1GTUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
 223 {
 224   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 225   m->writeVTKLL(ofs,cellData,pointData);
 226 }
 227
 228 std::string MEDCoupling1GTUMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
 229 {
 230   return std::string("UnstructuredGrid");
 231 }
 232
 233 std::size_t MEDCoupling1GTUMesh::getHeapMemorySize() const
 234 {
 235   return MEDCouplingPointSet::getHeapMemorySize();
 236 }
 237
 238 bool MEDCoupling1GTUMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
 239 {
 240   if(!MEDCouplingPointSet::isEqualIfNotWhy(other,prec,reason))
 241     return false;
 242   if(!other)
 243     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::isEqualIfNotWhy : input other pointer is null !");
 244   const MEDCoupling1GTUMesh *otherC=dynamic_cast<const MEDCoupling1GTUMesh *>(other);
 245   if(!otherC)
 246     {
 247       reason="mesh given in input is not castable in MEDCouplingSGTUMesh !";
 248       return false;
 249     }
 250   if(_cm!=otherC->_cm)
 251     {
 252       reason="mismatch in geometric type !";
 253       return false;
 254     }
 255   return true;
 256 }
 257
 258 bool MEDCoupling1GTUMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 259 {
 260   if(!MEDCouplingPointSet::isEqualWithoutConsideringStr(other,prec))
 261     return false;
 262   if(!other)
 263     throw INTERP_KERNEL::Exception("MEDCoupling1GTUMesh::isEqualWithoutConsideringStr : input other pointer is null !");
 264   const MEDCoupling1GTUMesh *otherC=dynamic_cast<const MEDCoupling1GTUMesh *>(other);
 265   if(!otherC)
 266     return false;
 267   if(_cm!=otherC->_cm)
 268     return false;
 269   return true;
 270 }
 271
 272 void MEDCoupling1GTUMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
 273 {
 274   MEDCouplingPointSet::checkCoherency();
 275 }
 276
 277 DataArrayDouble *MEDCoupling1GTUMesh::getBarycenterAndOwner() const
 278 {
 279   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 280   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=m->getBarycenterAndOwner();
 281   return ret.retn();
 282 }
 283
 284 MEDCouplingFieldDouble *MEDCoupling1GTUMesh::getMeasureField(bool isAbs) const
 285 {
 286   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 287   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=m->getMeasureField(isAbs);
 288   ret->setMesh(this);
 289   return ret.retn();
 290 }
 291
 292 MEDCouplingFieldDouble *MEDCoupling1GTUMesh::getMeasureFieldOnNode(bool isAbs) const
 293 {
 294   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 295   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=m->getMeasureFieldOnNode(isAbs);
 296   ret->setMesh(this);
 297   return ret.retn();
 298 }
 299
 300 /*!
 301  * to improve perf !
 302  */
 303 int MEDCoupling1GTUMesh::getCellContainingPoint(const double *pos, double eps) const
 304 {
 305   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 306   return m->getCellContainingPoint(pos,eps);
 307 }
 308
 309 MEDCouplingFieldDouble *MEDCoupling1GTUMesh::buildOrthogonalField() const
 310 {
 311   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 312   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=m->buildOrthogonalField();
 313   ret->setMesh(this);
 314   return ret.retn();
 315 }
 316
 317 DataArrayInt *MEDCoupling1GTUMesh::getCellsInBoundingBox(const double *bbox, double eps) const
 318 {
 319   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 320   return m->getCellsInBoundingBox(bbox,eps);
 321 }
 322
 323 DataArrayInt *MEDCoupling1GTUMesh::getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps)
 324 {
 325   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 326   return m->getCellsInBoundingBox(bbox,eps);
 327 }
 328
 329 MEDCouplingPointSet *MEDCoupling1GTUMesh::buildFacePartOfMySelfNode(const int *start, const int *end, bool fullyIn) const
 330 {
 331   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 332   return m->buildFacePartOfMySelfNode(start,end,fullyIn);
 333 }
 334
 335 DataArrayInt *MEDCoupling1GTUMesh::findBoundaryNodes() const
 336 {
 337   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 338   return m->findBoundaryNodes();
 339 }
 340
 341 MEDCouplingPointSet *MEDCoupling1GTUMesh::buildBoundaryMesh(bool keepCoords) const
 342 {
 343   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 344   return m->buildBoundaryMesh(keepCoords);
 345 }
 346
 347 void MEDCoupling1GTUMesh::findCommonCells(int compType, int startCellId, DataArrayInt *& commonCellsArr, DataArrayInt *& commonCellsIArr) const throw(INTERP_KERNEL::Exception)
 348 {
 349   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
 350   m->findCommonCells(compType,startCellId,commonCellsArr,commonCellsIArr);
 351 }
 352
 353 //==
 354
 355 MEDCoupling1SGTUMesh::MEDCoupling1SGTUMesh(const MEDCoupling1SGTUMesh& other, bool recDeepCpy):MEDCoupling1GTUMesh(other,recDeepCpy),_conn(other._conn)
 356 {
 357   if(recDeepCpy)
 358     {
 359       const DataArrayInt *c(other._conn);
 360       if(c)
 361         _conn=c->deepCpy();
 362     }
 363 }
 364
 365 MEDCoupling1SGTUMesh::MEDCoupling1SGTUMesh(const char *name, const INTERP_KERNEL::CellModel& cm):MEDCoupling1GTUMesh(name,cm)
 366 {
 367 }
 368
 369 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::New(const char *name, INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 370 {
 371   if(type==INTERP_KERNEL::NORM_ERROR)
 372     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::New : NORM_ERROR is not a valid type to be used as base geometric type for a mesh !");
 373   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 374   if(cm.isDynamic())
 375     {
 376       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::New : the input geometric type " << cm.getRepr() << " is dynamic ! Only static types are allowed here !";
 377       throw INTERP_KERNEL::Exception(oss.str().c_str());
 378     }
 379   return new MEDCoupling1SGTUMesh(name,cm);
 380 }
 381
 382 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::clone(bool recDeepCpy) const
 383 {
 384   return new MEDCoupling1SGTUMesh(*this,recDeepCpy);
 385 }
 386
 387 void MEDCoupling1SGTUMesh::shallowCopyConnectivityFrom(const MEDCouplingPointSet *other) throw(INTERP_KERNEL::Exception)
 388 {
 389   if(!other)
 390     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::shallowCopyConnectivityFrom : input pointer is null !");
 391   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
 392   if(!otherC)
 393     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::shallowCopyConnectivityFrom : input pointer is not an MEDCoupling1SGTUMesh instance !");
 394   setNodalConnectivity(otherC->getNodalConnectivity());
 395 }
 396
 397 void MEDCoupling1SGTUMesh::updateTime() const
 398 {
 399   MEDCoupling1GTUMesh::updateTime();
 400   const DataArrayInt *c(_conn);
 401   if(c)
 402     updateTimeWith(*c);
 403 }
 404
 405 std::size_t MEDCoupling1SGTUMesh::getHeapMemorySize() const
 406 {
 407   std::size_t ret=0;
 408   const DataArrayInt *c(_conn);
 409   if(c)
 410     ret+=c->getHeapMemorySize();
 411   return MEDCoupling1GTUMesh::getHeapMemorySize()+ret;
 412 }
 413
 414 MEDCouplingMesh *MEDCoupling1SGTUMesh::deepCpy() const
 415 {
 416   return clone(true);
 417 }
 418
 419 bool MEDCoupling1SGTUMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
 420 {
 421   if(!other)
 422     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::isEqualIfNotWhy : input other pointer is null !");
 423   std::ostringstream oss; oss.precision(15);
 424   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
 425   if(!otherC)
 426     {
 427       reason="mesh given in input is not castable in MEDCoupling1SGTUMesh !";
 428       return false;
 429     }
 430   if(!MEDCoupling1GTUMesh::isEqualIfNotWhy(other,prec,reason))
 431     return false;
 432   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
 433   if(c1==c2)
 434     return true;
 435   if(!c1 || !c2)
 436     {
 437       reason="in connectivity of single static geometric type exactly one among this and other is null !";
 438       return false;
 439     }
 440   if(!c1->isEqualIfNotWhy(*c2,reason))
 441     {
 442       reason.insert(0,"Nodal connectivity DataArrayInt differ : ");
 443       return false;
 444     }
 445   return true;
 446 }
 447
 448 bool MEDCoupling1SGTUMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 449 {
 450   if(!other)
 451     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::isEqualWithoutConsideringStr : input other pointer is null !");
 452   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
 453   if(!otherC)
 454     return false;
 455   if(!MEDCoupling1GTUMesh::isEqualWithoutConsideringStr(other,prec))
 456     return false;
 457   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
 458   if(c1==c2)
 459     return true;
 460   if(!c1 || !c2)
 461     return false;
 462   if(!c1->isEqualWithoutConsideringStr(*c2))
 463     return false;
 464   return true;
 465 }
 466
 467 void MEDCoupling1SGTUMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
 468 {
 469   MEDCouplingPointSet::checkCoherency();
 470   const DataArrayInt *c1(_conn);
 471   if(c1)
 472     {
 473       if(c1->getNumberOfComponents()!=1)
 474         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to be with number of components set to one !");
 475       if(c1->getInfoOnComponent(0)!="")
 476         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to have no info on its single component !");
 477       c1->checkAllocated();
 478     }
 479   else
 480     throw INTERP_KERNEL::Exception("Nodal connectivity array not defined !");
 481 }
 482
 483 void MEDCoupling1SGTUMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
 484 {
 485   checkCoherency();
 486   const DataArrayInt *c1(_conn);
 487   int nbOfTuples=c1->getNumberOfTuples();
 488   int nbOfNodesPerCell=(int)_cm->getNumberOfNodes();
 489   if(nbOfTuples%nbOfNodesPerCell!=0)
 490     {
 491       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 !";
 492       throw INTERP_KERNEL::Exception(oss.str().c_str());
 493     }
 494   int nbOfNodes=getNumberOfNodes();
 495   int nbOfCells=nbOfTuples/nbOfNodesPerCell;
 496   const int *w(c1->begin());
 497   for(int i=0;i<nbOfCells;i++)
 498     for(int j=0;j<nbOfNodesPerCell;j++,w++)
 499       {
 500         if(*w<0 || *w>=nbOfNodes)
 501           {
 502             std::ostringstream oss; oss << "At node #" << j << " of  cell #" << i << ", is equal to " << *w << " must be in [0," << nbOfNodes << ") !";
 503             throw INTERP_KERNEL::Exception(oss.str().c_str());
 504           }
 505       }
 506 }
 507
 508 void MEDCoupling1SGTUMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
 509 {
 510   checkCoherency1(eps);
 511 }
 512
 513 int MEDCoupling1SGTUMesh::getNumberOfCells() const
 514 {
 515   int nbOfTuples=getNodalConnectivityLength();
 516   int nbOfNodesPerCell=getNumberOfNodesPerCell();
 517   if(nbOfTuples%nbOfNodesPerCell!=0)
 518     {
 519       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 !";
 520       throw INTERP_KERNEL::Exception(oss.str().c_str());
 521     }
 522   return nbOfTuples/nbOfNodesPerCell;
 523 }
 524
 525 int MEDCoupling1SGTUMesh::getNumberOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
 526 {
 527   checkNonDynamicGeoType();
 528   return (int)_cm->getNumberOfNodes();
 529 }
 530
 531 int MEDCoupling1SGTUMesh::getNodalConnectivityLength() const throw(INTERP_KERNEL::Exception)
 532 {
 533   const DataArrayInt *c1(_conn);
 534   if(!c1)
 535     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::getNodalConnectivityLength : no connectivity set !");
 536   if(c1->getNumberOfComponents()!=1)
 537     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::getNodalConnectivityLength : Nodal connectivity array set must have exactly one component !");
 538   if(!c1->isAllocated())
 539     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::getNodalConnectivityLength : Nodal connectivity array must be allocated !");
 540   return c1->getNumberOfTuples();
 541 }
 542
 543 DataArrayInt *MEDCoupling1SGTUMesh::computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
 544 {
 545   checkNonDynamicGeoType();
 546   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 547   ret->alloc(getNumberOfCells(),1);
 548   ret->fillWithValue((int)_cm->getNumberOfNodes());
 549   return ret.retn();
 550 }
 551
 552 DataArrayInt *MEDCoupling1SGTUMesh::computeNbOfFacesPerCell() const throw(INTERP_KERNEL::Exception)
 553 {
 554   checkNonDynamicGeoType();
 555   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 556   ret->alloc(getNumberOfCells(),1);
 557   ret->fillWithValue((int)_cm->getNumberOfSons());
 558   return ret.retn();
 559 }
 560
 561 void MEDCoupling1SGTUMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
 562 {
 563   int sz=getNumberOfNodesPerCell();
 564   conn.resize(sz);
 565   if(cellId>=0 && cellId<getNumberOfCells())
 566     std::copy(_conn->begin()+cellId*sz,_conn->begin()+(cellId+1)*sz,conn.begin());
 567   else
 568     {
 569       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getNodeIdsOfCell : request for cellId #" << cellId << " must be in [0," << getNumberOfCells() << ") !";
 570       throw INTERP_KERNEL::Exception(oss.str().c_str());
 571     }
 572 }
 573
 574 void MEDCoupling1SGTUMesh::checkNonDynamicGeoType() const throw(INTERP_KERNEL::Exception)
 575 {
 576   if(_cm->isDynamic())
 577     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkNonDynamicGeoType : internal error ! the internal geo type is dynamic ! should be static !");
 578 }
 579
 580 std::string MEDCoupling1SGTUMesh::simpleRepr() const
 581 {
 582   static const char msg0[]="No coordinates specified !";
 583   std::ostringstream ret;
 584   ret << "Single static geometic type (" << _cm->getRepr() << ") unstructured mesh with name : \"" << getName() << "\"\n";
 585   ret << "Description of mesh : \"" << getDescription() << "\"\n";
 586   int tmpp1,tmpp2;
 587   double tt=getTime(tmpp1,tmpp2);
 588   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
 589   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
 590   ret << "Mesh dimension : " << getMeshDimension() << "\nSpace dimension : ";
 591   if(_coords!=0)
 592     {
 593       const int spaceDim=getSpaceDimension();
 594       ret << spaceDim << "\nInfo attached on space dimension : ";
 595       for(int i=0;i<spaceDim;i++)
 596         ret << "\"" << _coords->getInfoOnComponent(i) << "\" ";
 597       ret << "\n";
 598     }
 599   else
 600     ret << msg0 << "\n";
 601   ret << "Number of nodes : ";
 602   if(_coords!=0)
 603     ret << getNumberOfNodes() << "\n";
 604   else
 605     ret << msg0 << "\n";
 606   ret << "Number of cells : ";
 607   if((const DataArrayInt *)_conn)
 608     {
 609       if(_conn->isAllocated())
 610         {
 611           if(_conn->getNumberOfComponents()==1)
 612             ret << getNumberOfCells() << "\n";
 613           else
 614             ret << "Nodal connectivity array specified and allocated but with not exactly one component !" << "\n";
 615         }
 616       else
 617         ret << "Nodal connectivity array specified but not allocated !" << "\n";
 618     }
 619   else
 620     ret << "No connectivity specified !" << "\n";
 621   ret << "Cell type : " << _cm->getRepr() << "\n";
 622   return ret.str();
 623 }
 624
 625 std::string MEDCoupling1SGTUMesh::advancedRepr() const
 626 {
 627   std::ostringstream ret;
 628   ret << simpleRepr();
 629   ret << "\nCoordinates array : \n___________________\n\n";
 630   if(_coords)
 631     _coords->reprWithoutNameStream(ret);
 632   else
 633     ret << "No array set !\n";
 634   ret << "\n\nConnectivity array : \n____________________\n\n";
 635   //
 636   if((const DataArrayInt *)_conn)
 637     {
 638       if(_conn->isAllocated())
 639         {
 640           if(_conn->getNumberOfComponents()==1)
 641             {
 642              int nbOfCells=getNumberOfCells();
 643              int sz=getNumberOfNodesPerCell();
 644              const int *connPtr=_conn->begin();
 645              for(int i=0;i<nbOfCells;i++,connPtr+=sz)
 646                {
 647                  ret << "Cell #" << i << " : ";
 648                  std::copy(connPtr,connPtr+sz,std::ostream_iterator<int>(ret," "));
 649                  ret << "\n";
 650                }
 651             }
 652           else
 653             ret << "Nodal connectivity array specified and allocated but with not exactly one component !" << "\n";
 654         }
 655       else
 656         ret << "Nodal connectivity array specified but not allocated !" << "\n";
 657     }
 658   else
 659     ret << "No connectivity specified !" << "\n";
 660   return ret.str();
 661 }
 662
 663 DataArrayDouble *MEDCoupling1SGTUMesh::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
 664 {
 665   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
 666   int spaceDim=getSpaceDimension();
 667   int nbOfCells=getNumberOfCells();//checkCoherency()
 668   int nbOfNodes=getNumberOfNodes();
 669   ret->alloc(nbOfCells,spaceDim);
 670   double *ptToFill=ret->getPointer();
 671   const double *coor=_coords->begin();
 672   const int *nodal=_conn->begin();
 673   int sz=getNumberOfNodesPerCell();
 674   double coeff=1./(double)sz;
 675   for(int i=0;i<nbOfCells;i++,ptToFill+=spaceDim)
 676     {
 677       std::fill(ptToFill,ptToFill+spaceDim,0.);
 678       for(int j=0;j<sz;j++,nodal++)
 679         if(*nodal>=0 && *nodal<nbOfNodes)
 680           std::transform(coor+spaceDim*nodal[0],coor+spaceDim*(nodal[0]+1),ptToFill,ptToFill,std::plus<double>());
 681         else
 682           {
 683             std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of nodeId #" << *nodal << " should be in [0," <<   nbOfNodes << ") !";
 684             throw INTERP_KERNEL::Exception(oss.str().c_str());
 685           }
 686       std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),coeff));
 687     }
 688   return ret.retn();
 689 }
 690
 691 void MEDCoupling1SGTUMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 692 {
 693   int nbCells=getNumberOfCells();
 694   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> o2n=DataArrayInt::New();
 695   o2n->useArray(old2NewBg,false,C_DEALLOC,nbCells,1);
 696   if(check)
 697     o2n=o2n->checkAndPreparePermutation();
 698   //
 699   const int *conn=_conn->begin();
 700   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> n2o=o2n->invertArrayO2N2N2O(nbCells);
 701   const int *n2oPtr=n2o->begin();
 702   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New();
 703   newConn->alloc(_conn->getNumberOfTuples(),1);
 704   newConn->copyStringInfoFrom(*_conn);
 705   int sz=getNumberOfNodesPerCell();
 706   //
 707   int *newC=newConn->getPointer();
 708   for(int i=0;i<nbCells;i++,newC+=sz)
 709     {
 710       int pos=n2oPtr[i];
 711       std::copy(conn+pos*sz,conn+(pos+1)*sz,newC);
 712     }
 713   _conn=newConn;
 714 }
 715
 716 /*!
 717  * Keeps from \a this only cells which constituing point id are in the ids specified by [\a begin,\a end).
 718  * The resulting cell ids are stored at the end of the 'cellIdsKept' parameter.
 719  * Parameter \a fullyIn specifies if a cell that has part of its nodes in ids array is kept or not.
 720  * If \a fullyIn is true only cells whose ids are \b fully contained in [\a begin,\a end) tab will be kept.
 721  *
 722  * \param [in] begin input start of array of node ids.
 723  * \param [in] end input end of array of node ids.
 724  * \param [in] fullyIn input that specifies if all node ids must be in [\a begin,\a end) array to consider cell to be in.
 725  * \param [in,out] cellIdsKeptArr array where all candidate cell ids are put at the end.
 726  */
 727 void MEDCoupling1SGTUMesh::fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const
 728 {
 729   int nbOfCells=getNumberOfCells();
 730   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIdsKept=DataArrayInt::New(); cellIdsKept->alloc(0,1);
 731   int tmp=-1;
 732   int sz=_conn->getMaxValue(tmp); sz=std::max(sz,0)+1;
 733   std::vector<bool> fastFinder(sz,false);
 734   for(const int *work=begin;work!=end;work++)
 735     if(*work>=0 && *work<sz)
 736       fastFinder[*work]=true;
 737   const int *conn=_conn->begin();
 738   int nbNodesPerCell=getNumberOfNodesPerCell();
 739   for(int i=0;i<nbOfCells;i++,conn+=nbNodesPerCell)
 740     {
 741       int ref=0,nbOfHit=0;
 742       for(int j=0;j<nbNodesPerCell;j++)
 743         if(conn[j]>=0)
 744           {
 745             ref++;
 746             if(fastFinder[conn[j]])
 747               nbOfHit++;
 748           }
 749       if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
 750         cellIdsKept->pushBackSilent(i);
 751     }
 752   cellIdsKeptArr=cellIdsKept.retn();
 753 }
 754
 755 MEDCouplingMesh *MEDCoupling1SGTUMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 756 {
 757   if(other->getType()!=SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED)
 758     throw INTERP_KERNEL::Exception("Merge of umesh only available with umesh single static geo type each other !");
 759   const MEDCoupling1SGTUMesh *otherC=static_cast<const MEDCoupling1SGTUMesh *>(other);
 760   return Merge1SGTUMeshes(this,otherC);
 761 }
 762
 763 MEDCouplingUMesh *MEDCoupling1SGTUMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
 764 {
 765   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(getName(),getMeshDimension());
 766   ret->setCoords(getCoords());
 767   const int *nodalConn=_conn->begin();
 768   int nbCells=getNumberOfCells();
 769   int nbNodesPerCell=getNumberOfNodesPerCell();
 770   int geoType=(int)getCellModelEnum();
 771   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c=DataArrayInt::New(); c->alloc(nbCells*(nbNodesPerCell+1),1);
 772   int *cPtr=c->getPointer();
 773   for(int i=0;i<nbCells;i++,nodalConn+=nbNodesPerCell)
 774     {
 775       *cPtr++=geoType;
 776       cPtr=std::copy(nodalConn,nodalConn+nbNodesPerCell,cPtr);
 777     }
 778   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cI=DataArrayInt::Range(0,(nbCells+1)*(nbNodesPerCell+1),nbNodesPerCell+1);
 779   ret->setConnectivity(c,cI,true);
 780   return ret.retn();
 781 }
 782
 783 DataArrayInt *MEDCoupling1SGTUMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
 784 {
 785   switch(policy)
 786     {
 787     case 0:
 788       return simplexizePol0();
 789     case 1:
 790       return simplexizePol1();
 791     case (int) INTERP_KERNEL::PLANAR_FACE_5:
 792       return simplexizePlanarFace5();
 793     case (int) INTERP_KERNEL::PLANAR_FACE_6:
 794       return simplexizePlanarFace6();
 795     default:
 796       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)");
 797     }
 798 }
 799
 800 /// @cond INTERNAL
 801
 802 struct MEDCouplingAccVisit
 803 {
 804   MEDCouplingAccVisit():_new_nb_of_nodes(0) { }
 805   int operator()(int val) { if(val!=-1) return _new_nb_of_nodes++; else return -1; }
 806   int _new_nb_of_nodes;
 807 };
 808
 809 /// @endcond
 810
 811 /*!
 812  * Finds nodes not used in any cell and returns an array giving a new id to every node
 813  * by excluding the unused nodes, for which the array holds -1. The result array is
 814  * a mapping in "Old to New" mode.
 815  *  \param [out] nbrOfNodesInUse - number of node ids present in the nodal connectivity.
 816  *  \return DataArrayInt * - a new instance of DataArrayInt. Its length is \a
 817  *          this->getNumberOfNodes(). It holds for each node of \a this mesh either -1
 818  *          if the node is unused or a new id else. The caller is to delete this
 819  *          array using decrRef() as it is no more needed.
 820  *  \throw If the coordinates array is not set.
 821  *  \throw If the nodal connectivity of cells is not defined.
 822  *  \throw If the nodal connectivity includes an invalid id.
 823  */
 824 DataArrayInt *MEDCoupling1SGTUMesh::getNodeIdsInUse(int& nbrOfNodesInUse) const throw(INTERP_KERNEL::Exception)
 825 {
 826   nbrOfNodesInUse=-1;
 827   int nbOfNodes=getNumberOfNodes();
 828   int nbOfCells=getNumberOfCells();
 829   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 830   ret->alloc(nbOfNodes,1);
 831   int *traducer=ret->getPointer();
 832   std::fill(traducer,traducer+nbOfNodes,-1);
 833   const int *conn=_conn->begin();
 834   int nbNodesPerCell=getNumberOfNodesPerCell();
 835   for(int i=0;i<nbOfCells;i++)
 836     for(int j=0;j<nbNodesPerCell;j++,conn++)
 837       if(*conn>=0 && *conn<nbOfNodes)
 838         traducer[*conn]=1;
 839       else
 840         {
 841           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getNodeIdsInUse : In cell #" << i  << " presence of node id " <<  conn[j] << " not in [0," << nbOfNodes << ") !";
 842           throw INTERP_KERNEL::Exception(oss.str().c_str());
 843         }
 844   nbrOfNodesInUse=(int)std::count(traducer,traducer+nbOfNodes,1);
 845   std::transform(traducer,traducer+nbOfNodes,traducer,MEDCouplingAccVisit());
 846   return ret.retn();
 847 }
 848
 849 /*!
 850  * Changes ids of nodes within the nodal connectivity arrays according to a permutation
 851  * array in "Old to New" mode. The node coordinates array is \b not changed by this method.
 852  * This method is a generalization of shiftNodeNumbersInConn().
 853  *  \warning This method performs no check of validity of new ids. **Use it with care !**
 854  *  \param [in] newNodeNumbersO2N - a permutation array, of length \a
 855  *         this->getNumberOfNodes(), in "Old to New" mode.
 856  *         See \ref MEDCouplingArrayRenumbering for more info on renumbering modes.
 857  *  \throw If the nodal connectivity of cells is not defined.
 858  */
 859 void MEDCoupling1SGTUMesh::renumberNodesInConn(const int *newNodeNumbersO2N)
 860 {
 861   getNumberOfCells();//only to check that all is well defined.
 862   _conn->transformWithIndArr(newNodeNumbersO2N,newNodeNumbersO2N+getNumberOfNodes());
 863   updateTime();
 864 }
 865
 866 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshes(const MEDCoupling1SGTUMesh *mesh1, const MEDCoupling1SGTUMesh *mesh2) throw(INTERP_KERNEL::Exception)
 867 {
 868   std::vector<const MEDCoupling1SGTUMesh *> tmp(2);
 869   tmp[0]=const_cast<MEDCoupling1SGTUMesh *>(mesh1); tmp[1]=const_cast<MEDCoupling1SGTUMesh *>(mesh2);
 870   return Merge1SGTUMeshes(tmp);
 871 }
 872
 873 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshes(std::vector<const MEDCoupling1SGTUMesh *>& a) throw(INTERP_KERNEL::Exception)
 874 {
 875   std::size_t sz=a.size();
 876   if(sz==0)
 877     return Merge1SGTUMeshesLL(a);
 878   for(std::size_t ii=0;ii<sz;ii++)
 879     if(!a[ii])
 880       {
 881         std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::Merge1SGTUMeshes : item #" << ii << " in input array of size "<< sz << " is empty !";
 882         throw INTERP_KERNEL::Exception(oss.str().c_str());
 883       }
 884   const INTERP_KERNEL::CellModel *cm=&(a[0]->getCellModel());
 885   for(std::size_t ii=0;ii<sz;ii++)
 886     if(&(a[ii]->getCellModel())!=cm)
 887       throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshes : all items must have the same geo type !");
 888   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> > bb(sz);
 889   std::vector< const MEDCoupling1SGTUMesh * > aa(sz);
 890   int spaceDim=-3;
 891   for(std::size_t i=0;i<sz && spaceDim==-3;i++)
 892     {
 893       const MEDCoupling1SGTUMesh *cur=a[i];
 894       const DataArrayDouble *coo=cur->getCoords();
 895       if(coo)
 896         spaceDim=coo->getNumberOfComponents();
 897     }
 898   if(spaceDim==-3)
 899     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshes : no spaceDim specified ! unable to perform merge !");
 900   for(std::size_t i=0;i<sz;i++)
 901     {
 902       bb[i]=a[i]->buildSetInstanceFromThis(spaceDim);
 903       aa[i]=bb[i];
 904     }
 905   return Merge1SGTUMeshesLL(aa);
 906 }
 907
 908 /*!
 909  * \throw If presence of a null instance in the input vector \a a.
 910  * \throw If a is empty
 911  */
 912 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords(std::vector<const MEDCoupling1SGTUMesh *>& a) throw(INTERP_KERNEL::Exception)
 913 {
 914   if(a.empty())
 915     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : input array must be NON EMPTY !");
 916   std::vector<const MEDCoupling1SGTUMesh *>::const_iterator it=a.begin();
 917   if(!(*it))
 918     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : null instance in the first element of input vector !");
 919   std::vector<const DataArrayInt *> ncs(a.size());
 920   int nbOfCells=(*it)->getNumberOfCells();
 921   const DataArrayDouble *coords=(*it)->getCoords();
 922   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
 923   int nbNodesPerCell=(*it)->getNumberOfNodesPerCell();
 924   ncs[0]=(*it)->getNodalConnectivity();
 925   it++;
 926   for(int i=1;it!=a.end();i++,it++)
 927     {
 928       if(!(*it))
 929         throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : presence of a null instance in the input vector !");
 930       if(cm!=&((*it)->getCellModel()))
 931         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1SGTUMeshes impossible !");
 932       (*it)->getNumberOfCells();//to check that all is OK
 933       ncs[i]=(*it)->getNodalConnectivity();
 934       if(coords!=(*it)->getCoords())
 935         throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords : not lying on same coords !");
 936     }
 937   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh("merge",*cm));
 938   ret->setCoords(coords);
 939   ret->_conn=DataArrayInt::Aggregate(ncs);
 940   return ret.retn();
 941 }
 942
 943 /*!
 944  * 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)
 945  */
 946 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::Merge1SGTUMeshesLL(std::vector<const MEDCoupling1SGTUMesh *>& a) throw(INTERP_KERNEL::Exception)
 947 {
 948   if(a.empty())
 949     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::Merge1SGTUMeshes : input array must be NON EMPTY !");
 950   std::vector<const MEDCoupling1SGTUMesh *>::const_iterator it=a.begin();
 951   int nbOfCells=(*it)->getNumberOfCells();
 952   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
 953   int nbNodesPerCell=(*it)->getNumberOfNodesPerCell();
 954   it++;
 955   for(;it!=a.end();it++)
 956     {
 957       if(cm!=&((*it)->getCellModel()))
 958         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1SGTUMeshes impossible !");
 959       nbOfCells+=(*it)->getNumberOfCells();
 960     }
 961   std::vector<const MEDCouplingPointSet *> aps(a.size());
 962   std::copy(a.begin(),a.end(),aps.begin());
 963   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> pts=MergeNodesArray(aps);
 964   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh("merge",*cm));
 965   ret->setCoords(pts);
 966   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c=DataArrayInt::New();
 967   c->alloc(nbOfCells*nbNodesPerCell,1);
 968   int *cPtr=c->getPointer();
 969   int offset=0;
 970   for(it=a.begin();it!=a.end();it++)
 971     {
 972       int curConnLgth=(*it)->getNodalConnectivityLength();
 973       const int *curC=(*it)->_conn->begin();
 974       cPtr=std::transform(curC,curC+curConnLgth,cPtr,std::bind2nd(std::plus<int>(),offset));
 975       offset+=(*it)->getNumberOfNodes();
 976     }
 977   //
 978   ret->setNodalConnectivity(c);
 979   return ret.retn();
 980 }
 981
 982 MEDCouplingPointSet *MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoords(const int *begin, const int *end) const
 983 {
 984   int ncell=getNumberOfCells();
 985   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh(getName(),*_cm));
 986   ret->setCoords(_coords);
 987   std::size_t nbOfElemsRet=std::distance(begin,end);
 988   const int *inConn=_conn->getConstPointer();
 989   int sz=getNumberOfNodesPerCell();
 990   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connRet=DataArrayInt::New(); connRet->alloc((int)nbOfElemsRet*sz,1);
 991   int *connPtr=connRet->getPointer();
 992   for(const int *work=begin;work!=end;work++,connPtr+=sz)
 993     {
 994       if(*work>=0 && *work<ncell)
 995         std::copy(inConn+(work[0])*sz,inConn+(work[0]+1)*sz,connPtr);
 996       else
 997         {
 998           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoords : On pos #" << std::distance(begin,work) << " input cell id =" << *work << " should be in [0," << ncell << ") !";
 999           throw INTERP_KERNEL::Exception(oss.str().c_str());
1000         }
1001     }
1002   ret->_conn=connRet;
1003   ret->copyTinyInfoFrom(this);
1004   return ret.retn();
1005 }
1006
1007 MEDCouplingPointSet *MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoords2(int start, int end, int step) const
1008 {
1009   int ncell=getNumberOfCells();
1010   int nbOfElemsRet=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoords2 : ");
1011   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh(getName(),*_cm));
1012   ret->setCoords(_coords);
1013   const int *inConn=_conn->getConstPointer();
1014   int sz=getNumberOfNodesPerCell();
1015   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connRet=DataArrayInt::New(); connRet->alloc((int)nbOfElemsRet*sz,1);
1016   int *connPtr=connRet->getPointer();
1017   int curId=start;
1018   for(int i=0;i<nbOfElemsRet;i++,connPtr+=sz,curId+=step)
1019     {
1020       if(curId>=0 && curId<ncell)
1021         std::copy(inConn+curId*sz,inConn+(curId+1)*sz,connPtr);
1022       else
1023         {
1024           std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::buildPartOfMySelfKeepCoords2 : On pos #" << i << " input cell id =" << curId  << " should be in [0," << ncell << ") !";
1025           throw INTERP_KERNEL::Exception(oss.str().c_str());
1026         }
1027     }
1028   ret->_conn=connRet;
1029   ret->copyTinyInfoFrom(this);
1030   return ret.retn();
1031 }
1032
1033 MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::buildSetInstanceFromThis(int spaceDim) const throw(INTERP_KERNEL::Exception)
1034 {
1035   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(new MEDCoupling1SGTUMesh(getName(),*_cm));
1036   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1;
1037   const DataArrayInt *nodalConn(_conn);
1038   if(!nodalConn)
1039     {
1040       tmp1=DataArrayInt::New(); tmp1->alloc(0,1);
1041     }
1042   else
1043     tmp1=_conn;
1044   ret->_conn=tmp1;
1045   if(!_coords)
1046     {
1047       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=DataArrayDouble::New(); coords->alloc(0,spaceDim);
1048       ret->setCoords(coords);
1049     }
1050   else
1051     ret->setCoords(_coords);
1052   return ret.retn();
1053 }
1054
1055 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePol0() throw(INTERP_KERNEL::Exception)
1056 {
1057   int nbOfCells=getNumberOfCells();
1058   if(getCellModelEnum()!=INTERP_KERNEL::NORM_QUAD4)
1059     return DataArrayInt::Range(0,nbOfCells,1);
1060   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(2*3*nbOfCells,1);
1061   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(2*nbOfCells,1);
1062   const int *c(_conn->begin());
1063   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1064   for(int i=0;i<nbOfCells;i++,c+=4,newConnPtr+=6,retPtr+=2)
1065     {
1066       newConnPtr[0]=c[0]; newConnPtr[1]=c[1]; newConnPtr[2]=c[2];
1067       newConnPtr[3]=c[0]; newConnPtr[4]=c[2]; newConnPtr[5]=c[3];
1068       retPtr[0]=i; retPtr[1]=i;
1069     }
1070   _conn=newConn;
1071   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TRI3);
1072   updateTime();
1073   return ret.retn();
1074 }
1075
1076 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePol1() throw(INTERP_KERNEL::Exception)
1077 {
1078   int nbOfCells=getNumberOfCells();
1079   if(getCellModelEnum()!=INTERP_KERNEL::NORM_QUAD4)
1080     return DataArrayInt::Range(0,nbOfCells,1);
1081   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(2*3*nbOfCells,1);
1082   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(2*nbOfCells,1);
1083   const int *c(_conn->begin());
1084   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1085   for(int i=0;i<nbOfCells;i++,c+=4,newConnPtr+=6,retPtr+=2)
1086     {
1087       newConnPtr[0]=c[0]; newConnPtr[1]=c[1]; newConnPtr[2]=c[3];
1088       newConnPtr[3]=c[1]; newConnPtr[4]=c[2]; newConnPtr[5]=c[3];
1089       retPtr[0]=i; retPtr[1]=i;
1090     }
1091   _conn=newConn;
1092   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TRI3);
1093   updateTime();
1094   return ret.retn();
1095 }
1096
1097 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePlanarFace5() throw(INTERP_KERNEL::Exception)
1098 {
1099   int nbOfCells=getNumberOfCells();
1100   if(getCellModelEnum()!=INTERP_KERNEL::NORM_HEXA8)
1101     return DataArrayInt::Range(0,nbOfCells,1);
1102   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(5*4*nbOfCells,1);
1103   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(5*nbOfCells,1);
1104   const int *c(_conn->begin());
1105   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1106   for(int i=0;i<nbOfCells;i++,c+=8,newConnPtr+=20,retPtr+=5)
1107     {
1108       for(int j=0;j<20;j++)
1109         newConnPtr[j]=c[INTERP_KERNEL::SPLIT_NODES_5_WO[j]];
1110       retPtr[0]=i; retPtr[1]=i; retPtr[2]=i; retPtr[3]=i; retPtr[4]=i;
1111     }
1112   _conn=newConn;
1113   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TETRA4);
1114   updateTime();
1115   return ret.retn();
1116 }
1117
1118 DataArrayInt *MEDCoupling1SGTUMesh::simplexizePlanarFace6() throw(INTERP_KERNEL::Exception)
1119 {
1120   int nbOfCells=getNumberOfCells();
1121   if(getCellModelEnum()!=INTERP_KERNEL::NORM_HEXA8)
1122     return DataArrayInt::Range(0,nbOfCells,1);
1123   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(6*4*nbOfCells,1);
1124   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(6*nbOfCells,1);
1125   const int *c(_conn->begin());
1126   int *retPtr(ret->getPointer()),*newConnPtr(newConn->getPointer());
1127   for(int i=0;i<nbOfCells;i++,c+=8,newConnPtr+=24,retPtr+=6)
1128     {
1129       for(int j=0;j<24;j++)
1130         newConnPtr[j]=c[INTERP_KERNEL::SPLIT_NODES_6_WO[j]];
1131       retPtr[0]=i; retPtr[1]=i; retPtr[2]=i; retPtr[3]=i; retPtr[4]=i; retPtr[5]=i;
1132     }
1133   _conn=newConn;
1134   _cm=&INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_TETRA4);
1135   updateTime();
1136   return ret.retn();
1137 }
1138
1139 void MEDCoupling1SGTUMesh::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
1140 {
1141   stream << "MEDCoupling1SGTUMesh C++ instance at " << this << ". Type=" << _cm->getRepr() << ". Name : \"" << getName() << "\".";
1142   stream << " Mesh dimension : " << getMeshDimension() << ".";
1143   if(!_coords)
1144     { stream << " No coordinates set !"; return ; }
1145   if(!_coords->isAllocated())
1146     { stream << " Coordinates set but not allocated !"; return ; }
1147   stream << " Space dimension : " << _coords->getNumberOfComponents() << "." << std::endl;
1148   stream << "Number of nodes : " << _coords->getNumberOfTuples() << ".";
1149   if(!(const DataArrayInt *)_conn)
1150     { stream << std::endl << "Nodal connectivity NOT set !"; return ; }
1151   if(_conn->isAllocated())
1152     {
1153       if(_conn->getNumberOfComponents()==1)
1154         stream << std::endl << "Number of cells : " << getNumberOfCells() << ".";
1155     }
1156 }
1157
1158 void MEDCoupling1SGTUMesh::checkFullyDefined() const throw(INTERP_KERNEL::Exception)
1159 {
1160   if(!((const DataArrayInt *)_conn) || !((const DataArrayDouble *)_coords))
1161     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFullyDefined : part of this is not fully defined.");
1162 }
1163
1164 /*!
1165  * First step of unserialization process.
1166  */
1167 bool MEDCoupling1SGTUMesh::isEmptyMesh(const std::vector<int>& tinyInfo) const
1168 {
1169   throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::isEmptyMesh : not implemented yet !");
1170 }
1171
1172 /*!
1173  * Checks if \a this and \a other meshes are geometrically equivalent with high
1174  * probability, else an exception is thrown. The meshes are considered equivalent if
1175  * (1) meshes contain the same number of nodes and the same number of elements of the
1176  * same types (2) three cells of the two meshes (first, last and middle) are based
1177  * on coincident nodes (with a specified precision).
1178  *  \param [in] other - the mesh to compare with.
1179  *  \param [in] prec - the precision used to compare nodes of the two meshes.
1180  *  \throw If the two meshes do not match.
1181  */
1182 void MEDCoupling1SGTUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
1183 {
1184   MEDCouplingPointSet::checkFastEquivalWith(other,prec);
1185   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
1186   if(!otherC)
1187     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : Two meshes are not not unstructured with single static geometric type !");
1188   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
1189   if(c1==c2)
1190     return;
1191   if(!c1 || !c2)
1192     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : presence of nodal connectivity only in one of the 2 meshes !");
1193   if((c1->isAllocated() && !c2->isAllocated()) || (!c1->isAllocated() && c2->isAllocated()))
1194     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : in nodal connectivity, only one is allocated !");
1195   if(c1->getNumberOfComponents()!=1 || c1->getNumberOfComponents()!=1)
1196     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : in nodal connectivity, must have 1 and only 1 component !");
1197   if(c1->getHashCode()!=c2->getHashCode())
1198     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::checkFastEquivalWith : nodal connectivity differs");
1199 }
1200
1201 MEDCouplingPointSet *MEDCoupling1SGTUMesh::mergeMyselfWithOnSameCoords(const MEDCouplingPointSet *other) const
1202 {
1203   if(!other)
1204     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::mergeMyselfWithOnSameCoords : input other is null !");
1205   const MEDCoupling1SGTUMesh *otherC=dynamic_cast<const MEDCoupling1SGTUMesh *>(other);
1206   if(!otherC)
1207     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::mergeMyselfWithOnSameCoords : the input other mesh is not of type single statuc geo type unstructured !");
1208   std::vector<const MEDCoupling1SGTUMesh *> ms(2);
1209   ms[0]=this;
1210   ms[1]=otherC;
1211   return Merge1SGTUMeshesOnSameCoords(ms);
1212 }
1213
1214 void MEDCoupling1SGTUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const throw(INTERP_KERNEL::Exception)
1215 {
1216   checkFullyDefined();
1217   int nbOfNodes=getNumberOfNodes();
1218   int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int));
1219   revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1);
1220   std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
1221   const int *conn=_conn->begin();
1222   int nbOfCells=getNumberOfCells();
1223   int nbOfEltsInRevNodal=0;
1224   int nbOfNodesPerCell=getNumberOfNodesPerCell();
1225   for(int eltId=0;eltId<nbOfCells;eltId++)
1226     {
1227       for(int j=0;j<nbOfNodesPerCell;j++,conn++)
1228         {
1229           if(conn[0]>=0 && conn[0]<nbOfNodes)
1230             {
1231               nbOfEltsInRevNodal++;
1232               revNodalIndxPtr[conn[0]+1]++;
1233             }
1234           else
1235             {
1236               std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getReverseNodalConnectivity : At cell #" << eltId << " presence of nodeId #" << conn[0] << " should be in [0," << nbOfNodes << ") !";
1237               throw INTERP_KERNEL::Exception(oss.str().c_str());
1238             }
1239         }
1240     }
1241   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<int>());
1242   conn=_conn->begin();
1243   int *revNodalPtr=(int *)malloc((nbOfEltsInRevNodal)*sizeof(int));
1244   revNodal->useArray(revNodalPtr,true,C_DEALLOC,nbOfEltsInRevNodal,1);
1245   std::fill(revNodalPtr,revNodalPtr+nbOfEltsInRevNodal,-1);
1246   for(int eltId=0;eltId<nbOfCells;eltId++)
1247     {
1248       for(int j=0;j<nbOfNodesPerCell;j++,conn++)
1249         {
1250           *std::find_if(revNodalPtr+revNodalIndxPtr[*conn],revNodalPtr+revNodalIndxPtr[*conn+1],std::bind2nd(std::equal_to<int>(),-1))=eltId;
1251         }
1252     }
1253 }
1254
1255 /*!
1256  * Use \a nodalConn array as nodal connectivity of \a this. The input \a nodalConn pointer can be null.
1257  */
1258 void MEDCoupling1SGTUMesh::setNodalConnectivity(DataArrayInt *nodalConn) throw(INTERP_KERNEL::Exception)
1259 {
1260   if(nodalConn)
1261     nodalConn->incrRef();
1262   _conn=nodalConn;
1263   declareAsNew();
1264 }
1265
1266 /*!
1267  * \return DataArrayInt * - the internal reference to the nodal connectivity. The caller is not reponsible to deallocate it.
1268  */
1269 DataArrayInt *MEDCoupling1SGTUMesh::getNodalConnectivity() const throw(INTERP_KERNEL::Exception)
1270 {
1271   const DataArrayInt *ret(_conn);
1272   return const_cast<DataArrayInt *>(ret);
1273 }
1274
1275 /*!
1276  * Allocates memory to store an estimation of the given number of cells. Closer is the estimation to the number of cells effectively inserted,
1277  * less will be the needs to realloc. If the number of cells to be inserted is not known simply put 0 to this parameter.
1278  * If a nodal connectivity previouly existed before the call of this method, it will be reset.
1279  *
1280  *  \param [in] nbOfCells - estimation of the number of cell \a this mesh will contain.
1281  */
1282 void MEDCoupling1SGTUMesh::allocateCells(int nbOfCells) throw(INTERP_KERNEL::Exception)
1283 {
1284   if(nbOfCells<0)
1285     throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::allocateCells : the input number of cells should be >= 0 !");
1286   _conn=DataArrayInt::New();
1287   _conn->reserve(getNumberOfNodesPerCell()*nbOfCells);
1288   declareAsNew();
1289 }
1290
1291 /*!
1292  * Appends at the end of \a this a cell having nodal connectivity array defined in [ \a nodalConnOfCellBg, \a nodalConnOfCellEnd ).
1293  *
1294  * \param [in] nodalConnOfCellBg - the begin (included) of nodal connectivity of the cell to add.
1295  * \param [in] nodalConnOfCellEnd - the end (excluded) of nodal connectivity of the cell to add.
1296  * \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
1297  *        attached to \a this.
1298  * \thow If the nodal connectivity array in \a this is null (call MEDCoupling1SGTUMesh::allocateCells before).
1299  */
1300 void MEDCoupling1SGTUMesh::insertNextCell(const int *nodalConnOfCellBg, const int *nodalConnOfCellEnd) throw(INTERP_KERNEL::Exception)
1301 {
1302   int sz=(int)std::distance(nodalConnOfCellBg,nodalConnOfCellEnd);
1303   int ref=getNumberOfNodesPerCell();
1304   if(sz==ref)
1305     {
1306       DataArrayInt *c(_conn);
1307       if(c)
1308         c->pushBackValsSilent(nodalConnOfCellBg,nodalConnOfCellEnd);
1309       else
1310         throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::insertNextCell : nodal connectivity array is null ! Call MEDCoupling1SGTUMesh::allocateCells before !");
1311     }
1312   else
1313     {
1314       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::insertNextCell : input nodal size (" << sz << ") does not match number of nodes per cell of this (";
1315       oss << ref << ") !";
1316       throw INTERP_KERNEL::Exception(oss.str().c_str());
1317     }
1318 }
1319
1320 //== find static tony
1321
1322 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::New(const char *name, INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
1323 {
1324   if(type==INTERP_KERNEL::NORM_ERROR)
1325     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::New : NORM_ERROR is not a valid type to be used as base geometric type for a mesh !");
1326   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1327   if(!cm.isDynamic())
1328     {
1329       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::New : the input geometric type " << cm.getRepr() << " is static ! Only dynamic types are allowed here !";
1330       throw INTERP_KERNEL::Exception(oss.str().c_str());
1331     }
1332   return new MEDCoupling1DGTUMesh(name,cm);
1333 }
1334
1335 MEDCoupling1DGTUMesh::MEDCoupling1DGTUMesh(const char *name, const INTERP_KERNEL::CellModel& cm):MEDCoupling1GTUMesh(name,cm)
1336 {
1337 }
1338
1339 MEDCoupling1DGTUMesh::MEDCoupling1DGTUMesh(const MEDCoupling1DGTUMesh& other, bool recDeepCpy):MEDCoupling1GTUMesh(other,recDeepCpy),_conn(other._conn)
1340 {
1341   if(recDeepCpy)
1342     {
1343       const DataArrayInt *c(other._conn);
1344       if(c)
1345         _conn=c->deepCpy();
1346       c=other._conn_indx;
1347       if(c)
1348         _conn_indx=c->deepCpy();
1349     }
1350 }
1351
1352 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::clone(bool recDeepCpy) const
1353 {
1354   return new MEDCoupling1DGTUMesh(*this,recDeepCpy);
1355 }
1356
1357 void MEDCoupling1DGTUMesh::updateTime() const
1358 {
1359   MEDCoupling1GTUMesh::updateTime();
1360   const DataArrayInt *c(_conn);
1361   if(c)
1362     updateTimeWith(*c);
1363   c=_conn_indx;
1364   if(c)
1365     updateTimeWith(*c);
1366 }
1367
1368 std::size_t MEDCoupling1DGTUMesh::getHeapMemorySize() const
1369 {
1370   std::size_t ret=0;
1371   const DataArrayInt *c(_conn);
1372   if(c)
1373     ret+=c->getHeapMemorySize();
1374   c=_conn_indx;
1375   if(c)
1376     ret+=c->getHeapMemorySize();
1377   return MEDCoupling1GTUMesh::getHeapMemorySize()+ret;
1378 }
1379
1380 MEDCouplingMesh *MEDCoupling1DGTUMesh::deepCpy() const
1381 {
1382   return clone(true);
1383 }
1384
1385 bool MEDCoupling1DGTUMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
1386 {
1387   if(!other)
1388     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::isEqualIfNotWhy : input other pointer is null !");
1389   std::ostringstream oss; oss.precision(15);
1390   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
1391   if(!otherC)
1392     {
1393       reason="mesh given in input is not castable in MEDCoupling1DGTUMesh !";
1394       return false;
1395     }
1396   if(!MEDCoupling1GTUMesh::isEqualIfNotWhy(other,prec,reason))
1397     return false;
1398   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
1399   if(c1==c2)
1400     return true;
1401   if(!c1 || !c2)
1402     {
1403       reason="in connectivity of single dynamic geometric type exactly one among this and other is null !";
1404       return false;
1405     }
1406   if(!c1->isEqualIfNotWhy(*c2,reason))
1407     {
1408       reason.insert(0,"Nodal connectivity DataArrayInt differs : ");
1409       return false;
1410     }
1411   c1=_conn_indx; c2=otherC->_conn_indx;
1412   if(c1==c2)
1413     return true;
1414   if(!c1 || !c2)
1415     {
1416       reason="in connectivity index of single dynamic geometric type exactly one among this and other is null !";
1417       return false;
1418     }
1419   if(!c1->isEqualIfNotWhy(*c2,reason))
1420     {
1421       reason.insert(0,"Nodal connectivity index DataArrayInt differs : ");
1422       return false;
1423     }
1424   return true;
1425 }
1426
1427 bool MEDCoupling1DGTUMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
1428 {
1429   if(!other)
1430     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::isEqualWithoutConsideringStr : input other pointer is null !");
1431   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
1432   if(!otherC)
1433     return false;
1434   if(!MEDCoupling1GTUMesh::isEqualWithoutConsideringStr(other,prec))
1435     return false;
1436   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
1437   if(c1==c2)
1438     return true;
1439   if(!c1 || !c2)
1440     return false;
1441   if(!c1->isEqualWithoutConsideringStr(*c2))
1442     return false;
1443   return true;
1444   c1=_conn_indx; c2=otherC->_conn_indx;
1445   if(c1==c2)
1446     return true;
1447   if(!c1 || !c2)
1448     return false;
1449   if(!c1->isEqualWithoutConsideringStr(*c2))
1450     return false;
1451   return true;
1452 }
1453
1454 /*!
1455  * Checks if \a this and \a other meshes are geometrically equivalent with high
1456  * probability, else an exception is thrown. The meshes are considered equivalent if
1457  * (1) meshes contain the same number of nodes and the same number of elements of the
1458  * same types (2) three cells of the two meshes (first, last and middle) are based
1459  * on coincident nodes (with a specified precision).
1460  *  \param [in] other - the mesh to compare with.
1461  *  \param [in] prec - the precision used to compare nodes of the two meshes.
1462  *  \throw If the two meshes do not match.
1463  */
1464 void MEDCoupling1DGTUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
1465 {
1466   MEDCouplingPointSet::checkFastEquivalWith(other,prec);
1467   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
1468   if(!otherC)
1469     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : Two meshes are not not unstructured with single static geometric type !");
1470   const DataArrayInt *c1(_conn),*c2(otherC->_conn);
1471   if(c1!=c2)
1472     {
1473       if(!c1 || !c2)
1474         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : presence of nodal connectivity only in one of the 2 meshes !");
1475       if((c1->isAllocated() && !c2->isAllocated()) || (!c1->isAllocated() && c2->isAllocated()))
1476         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity, only one is allocated !");
1477       if(c1->getNumberOfComponents()!=1 || c1->getNumberOfComponents()!=1)
1478         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity, must have 1 and only 1 component !");
1479       if(c1->getHashCode()!=c2->getHashCode())
1480         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : nodal connectivity differs");
1481     }
1482   c1=_conn_indx; c2=otherC->_conn_indx;
1483   if(c1!=c2)
1484     {
1485       if(!c1 || !c2)
1486         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : presence of nodal connectivity index only in one of the 2 meshes !");
1487       if((c1->isAllocated() && !c2->isAllocated()) || (!c1->isAllocated() && c2->isAllocated()))
1488         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity index, only one is allocated !");
1489       if(c1->getNumberOfComponents()!=1 || c1->getNumberOfComponents()!=1)
1490         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : in nodal connectivity index, must have 1 and only 1 component !");
1491       if(c1->getHashCode()!=c2->getHashCode())
1492         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFastEquivalWith : nodal connectivity index differs");
1493     }
1494 }
1495
1496 /*!
1497  * 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.
1498  * In addition you are sure that the length of nodal connectivity index array is bigger than or equal to one.
1499  * In addition you are also sure that length of nodal connectivity is coherent with the content of the last value in the index array.
1500  */
1501 void MEDCoupling1DGTUMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
1502 {
1503   MEDCouplingPointSet::checkCoherency();
1504   const DataArrayInt *c1(_conn);
1505   if(c1)
1506     {
1507       if(c1->getNumberOfComponents()!=1)
1508         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to be with number of components set to one !");
1509       if(c1->getInfoOnComponent(0)!="")
1510         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to have no info on its single component !");
1511       c1->checkAllocated();
1512     }
1513   else
1514     throw INTERP_KERNEL::Exception("Nodal connectivity array not defined !");
1515   //
1516   int sz2=_conn->getNumberOfTuples();
1517   c1=_conn_indx;
1518   if(c1)
1519     {
1520       if(c1->getNumberOfComponents()!=1)
1521         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to be with number of components set to one !");
1522       c1->checkAllocated();
1523       if(c1->getNumberOfTuples()<1)
1524         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to have a a size of 1 at least !");
1525       if(c1->getInfoOnComponent(0)!="")
1526         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to have no info on its single component !");
1527       int f=c1->front(),ll=c1->back();
1528       if(f<0 || f>=sz2)
1529         {
1530           std::ostringstream oss; oss << "Nodal connectivity index array first value (" << f << ") is expected to be exactly in [0," << sz2 << ") !";
1531           throw INTERP_KERNEL::Exception(oss.str().c_str());
1532         }
1533       if(ll<0 || ll>sz2)
1534         {
1535           std::ostringstream oss; oss << "Nodal connectivity index array last value (" << ll << ") is expected to be exactly in [0," << sz2 << "] !";
1536           throw INTERP_KERNEL::Exception(oss.str().c_str());
1537         }
1538       if(f>ll)
1539         {
1540           std::ostringstream oss; oss << "Nodal connectivity index array looks very bad (not increasing monotonic) because front (" << f << ") is greater that back (" << ll << ") !";
1541           throw INTERP_KERNEL::Exception(oss.str().c_str());
1542         }
1543     }
1544   else
1545     throw INTERP_KERNEL::Exception("Nodal connectivity index array not defined !");
1546   int szOfC1Exp=_conn_indx->back();
1547   if(sz2<szOfC1Exp)
1548     {
1549       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::checkCoherency : The expected length of nodal connectivity array regarding index is " << szOfC1Exp << " but the actual size of it is " << c1->getNumberOfTuples() << " !";
1550       throw INTERP_KERNEL::Exception(oss.str().c_str());
1551     }
1552 }
1553
1554 void MEDCoupling1DGTUMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
1555 {
1556   checkCoherency();
1557   const DataArrayInt *c1(_conn),*c2(_conn_indx);
1558   if(!c2->isMonotonic(true))
1559     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkCoherency1 : the nodal connectivity index is expected to be increasing monotinic !");
1560   //
1561   int nbOfTuples=c1->getNumberOfTuples();
1562   int nbOfNodes=getNumberOfNodes();
1563   const int *w(c1->begin());
1564   for(int i=0;i<nbOfTuples;i++,w++)
1565     {
1566       if(*w==-1) continue;
1567       if(*w<0 || *w>=nbOfNodes)
1568         {
1569           std::ostringstream oss; oss << "At pos #" << i << " of nodal connectivity array references to node id #" << *w << " must be in [0," << nbOfNodes << ") !";
1570           throw INTERP_KERNEL::Exception(oss.str().c_str());
1571         }
1572     }
1573 }
1574
1575 void MEDCoupling1DGTUMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
1576 {
1577   checkCoherency1(eps);
1578 }
1579
1580 int MEDCoupling1DGTUMesh::getNumberOfCells() const
1581 {
1582   checkCoherency();//do not remove
1583   return _conn_indx->getNumberOfTuples()-1;
1584 }
1585
1586 /*!
1587  * This method returns a newly allocated array containing this->getNumberOfCells() tuples and 1 component.
1588  * For each cell in \b this the number of nodes constituting cell is computed.
1589  * For each polyhedron cell, the sum of the number of nodes of each face constituting polyhedron cell is returned.
1590  * So for pohyhedrons some nodes can be counted several times in the returned result.
1591  *
1592  * \return a newly allocated array
1593  */
1594 DataArrayInt *MEDCoupling1DGTUMesh::computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
1595 {
1596   checkCoherency();
1597   _conn_indx->checkMonotonic(true);
1598   if(getCellModelEnum()!=INTERP_KERNEL::NORM_POLYHED)
1599     return _conn_indx->deltaShiftIndex();
1600   // for polyhedrons
1601   int nbOfCells=_conn_indx->getNumberOfTuples()-1;
1602   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1603   ret->alloc(nbOfCells,1);
1604   int *retPtr=ret->getPointer();
1605   const int *ci=_conn_indx->begin(),*c=_conn->begin();
1606   for(int i=0;i<nbOfCells;i++,retPtr++,ci++)
1607     *retPtr=ci[1]-ci[0]-std::count(c+ci[0],c+ci[1],-1);
1608   return ret.retn();
1609 }
1610
1611 /*!
1612  * This method returns a newly allocated array containing this->getNumberOfCells() tuples and 1 component.
1613  * For each cell in \b this the number of faces constituting (entity of dimension this->getMeshDimension()-1) cell is computed.
1614  *
1615  * \return a newly allocated array
1616  */
1617 DataArrayInt *MEDCoupling1DGTUMesh::computeNbOfFacesPerCell() const throw(INTERP_KERNEL::Exception)
1618 {
1619   checkCoherency();
1620   _conn_indx->checkMonotonic(true);
1621   if(getCellModelEnum()!=INTERP_KERNEL::NORM_POLYHED && getCellModelEnum()!=INTERP_KERNEL::NORM_QPOLYG)
1622     return _conn_indx->deltaShiftIndex();
1623   if(getCellModelEnum()==INTERP_KERNEL::NORM_QPOLYG)
1624     {
1625       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=_conn_indx->deltaShiftIndex();
1626       ret->applyDivideBy(2);
1627       return ret.retn();
1628     }
1629   // for polyhedrons
1630   int nbOfCells=_conn_indx->getNumberOfTuples()-1;
1631   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1632   ret->alloc(nbOfCells,1);
1633   int *retPtr=ret->getPointer();
1634   const int *ci=_conn_indx->begin(),*c=_conn->begin();
1635   for(int i=0;i<nbOfCells;i++,retPtr++,ci++)
1636     *retPtr=std::count(c+ci[0],c+ci[1],-1)+1;
1637   return ret.retn();
1638 }
1639
1640 void MEDCoupling1DGTUMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
1641 {
1642   int nbOfCells=getNumberOfCells();//performs checks
1643   if(cellId>=0 && cellId<nbOfCells)
1644     {
1645       int strt=_conn_indx->getIJ(cellId,0),stp=_conn_indx->getIJ(cellId+1,0);
1646       int nbOfNodes=stp-strt;
1647       if(nbOfNodes<0)
1648         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::getNodeIdsOfCell : the index array is invalid ! Should be increasing monotonic !");
1649       conn.resize(nbOfNodes);
1650       std::copy(_conn->begin()+strt,_conn->begin()+stp,conn.begin());
1651     }
1652   else
1653     {
1654       std::ostringstream oss; oss << "MEDCoupling1SGTUMesh::getNodeIdsOfCell : request for cellId #" << cellId << " must be in [0," << nbOfCells << ") !";
1655       throw INTERP_KERNEL::Exception(oss.str().c_str());
1656     }
1657 }
1658
1659 std::string MEDCoupling1DGTUMesh::simpleRepr() const
1660 {
1661   static const char msg0[]="No coordinates specified !";
1662   std::ostringstream ret;
1663   ret << "Single dynamic geometic type (" << _cm->getRepr() << ") unstructured mesh with name : \"" << getName() << "\"\n";
1664   ret << "Description of mesh : \"" << getDescription() << "\"\n";
1665   int tmpp1,tmpp2;
1666   double tt=getTime(tmpp1,tmpp2);
1667   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
1668   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
1669   ret << "Mesh dimension : " << getMeshDimension() << "\nSpace dimension : ";
1670   if(_coords!=0)
1671     {
1672       const int spaceDim=getSpaceDimension();
1673       ret << spaceDim << "\nInfo attached on space dimension : ";
1674       for(int i=0;i<spaceDim;i++)
1675         ret << "\"" << _coords->getInfoOnComponent(i) << "\" ";
1676       ret << "\n";
1677     }
1678   else
1679     ret << msg0 << "\n";
1680   ret << "Number of nodes : ";
1681   if(_coords!=0)
1682     ret << getNumberOfNodes() << "\n";
1683   else
1684     ret << msg0 << "\n";
1685   ret << "Number of cells : ";
1686   bool isOK=true;
1687   try { checkCoherency(); } catch(INTERP_KERNEL::Exception& e)
1688     {
1689       ret << "Nodal connectivity arrays are not set or badly set !\n";
1690       isOK=false;
1691     }
1692   if(isOK)
1693     ret << getNumberOfCells() << "\n";
1694   ret << "Cell type : " << _cm->getRepr() << "\n";
1695   return ret.str();
1696 }
1697
1698 std::string MEDCoupling1DGTUMesh::advancedRepr() const
1699 {
1700   std::ostringstream ret;
1701   ret << simpleRepr();
1702   ret << "\nCoordinates array : \n___________________\n\n";
1703   if(_coords)
1704     _coords->reprWithoutNameStream(ret);
1705   else
1706     ret << "No array set !\n";
1707   ret << "\n\nNodal Connectivity : \n____________________\n\n";
1708   //
1709   bool isOK=true;
1710   try { checkCoherency1(); } catch(INTERP_KERNEL::Exception& e)
1711     {
1712       ret << "Nodal connectivity arrays are not set or badly set !\n";
1713       isOK=false;
1714     }
1715   if(!isOK)
1716     return ret.str();
1717   int nbOfCells=getNumberOfCells();
1718   const int *ci=_conn_indx->begin(),*c=_conn->begin();
1719   for(int i=0;i<nbOfCells;i++,ci++)
1720     {
1721       ret << "Cell #" << i << " : ";
1722       std::copy(c+ci[0],c+ci[1],std::ostream_iterator<int>(ret," "));
1723       ret << "\n";
1724     }
1725   return ret.str();
1726 }
1727
1728 DataArrayDouble *MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
1729 {
1730   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1731   int spaceDim=getSpaceDimension();
1732   int nbOfCells=getNumberOfCells();//checkCoherency()
1733   int nbOfNodes=getNumberOfNodes();
1734   ret->alloc(nbOfCells,spaceDim);
1735   double *ptToFill=ret->getPointer();
1736   const double *coor=_coords->begin();
1737   const int *nodal=_conn->begin(),*nodali=_conn_indx->begin();
1738   nodal+=nodali[0];
1739   if(getCellModelEnum()!=INTERP_KERNEL::NORM_POLYHED)
1740     {
1741       for(int i=0;i<nbOfCells;i++,ptToFill+=spaceDim,nodali++)
1742         {
1743           std::fill(ptToFill,ptToFill+spaceDim,0.);
1744           if(nodali[0]>=nodali[1])// >= to avoid division by 0.
1745             {
1746               for(int j=nodali[0];j<nodali[1];j++,nodal++)
1747                 {
1748                   if(*nodal>=0 && *nodal<nbOfNodes)
1749                     std::transform(coor+spaceDim*nodal[0],coor+spaceDim*(nodal[0]+1),ptToFill,ptToFill,std::plus<double>());
1750                   else
1751                     {
1752                       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of nodeId #" << *nodal << " should be in [0," <<   nbOfNodes << ") !";
1753                       throw INTERP_KERNEL::Exception(oss.str().c_str());
1754                     }
1755                   std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(nodali[1]-nodali[0])));
1756                 }
1757             }
1758           else
1759             {
1760               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell : at cell #" << i << " the nodal index array is invalid !";
1761               throw INTERP_KERNEL::Exception(oss.str().c_str());
1762             }
1763         }
1764     }
1765   else
1766     {
1767       for(int i=0;i<nbOfCells;i++,ptToFill+=spaceDim,nodali++)
1768         {
1769           std::fill(ptToFill,ptToFill+spaceDim,0.);
1770           if(nodali[0]>=nodali[1])// >= to avoid division by 0.
1771             {
1772               int nbOfNod=0;
1773               for(int j=nodali[0];j<nodali[1];j++,nodal++)
1774                 {
1775                   if(*nodal==-1) continue;
1776                   if(*nodal>=0 && *nodal<nbOfNodes)
1777                     {
1778                       std::transform(coor+spaceDim*nodal[0],coor+spaceDim*(nodal[0]+1),ptToFill,ptToFill,std::plus<double>());
1779                       nbOfNod++;
1780                     }
1781                   else
1782                     {
1783                       std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell (polyhedron) : on cell #" << i << " presence of nodeId #" << *nodal << " should be in [0," <<   nbOfNodes << ") !";
1784                       throw INTERP_KERNEL::Exception(oss.str().c_str());
1785                     }
1786                 }
1787               if(nbOfNod!=0)
1788                 std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./nbOfNod));
1789               else
1790                 {
1791                   std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell (polyhedron) : no nodes in cell #" << i << " !";
1792                   throw INTERP_KERNEL::Exception(oss.str().c_str());
1793                 }
1794             }
1795           else
1796             {
1797               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::computeIsoBarycenterOfNodesPerCell (polyhedron)  : at cell #" << i << " the nodal index array is invalid !";
1798               throw INTERP_KERNEL::Exception(oss.str().c_str());
1799             }
1800         }
1801     }
1802   return ret.retn();
1803 }
1804
1805 void MEDCoupling1DGTUMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1806 {
1807   int nbCells=getNumberOfCells();
1808   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> o2n=DataArrayInt::New();
1809   o2n->useArray(old2NewBg,false,C_DEALLOC,nbCells,1);
1810   if(check)
1811     o2n=o2n->checkAndPreparePermutation();
1812   //
1813   const int *conn=_conn->begin(),*conni=_conn_indx->begin();
1814   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> n2o=o2n->invertArrayO2N2N2O(nbCells);
1815   const int *n2oPtr=n2o->begin();
1816   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New();
1817   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConnI=DataArrayInt::New();
1818   newConn->alloc(_conn->getNumberOfTuples(),1); newConnI->alloc(nbCells,1);
1819   newConn->copyStringInfoFrom(*_conn); newConnI->copyStringInfoFrom(*_conn_indx);
1820   //
1821   int *newC=newConn->getPointer(),*newCI=newConnI->getPointer();
1822   for(int i=0;i<nbCells;i++)
1823     {
1824       int sz=conni[i+1]-conni[i];
1825       if(sz>=0)
1826         newCI[n2oPtr[i]]=sz;
1827       else
1828         {
1829           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::renumberCells : the index nodal array is invalid for cell #" << i << " !";
1830           throw INTERP_KERNEL::Exception(oss.str().c_str());
1831         }
1832     }
1833   newConnI->computeOffsets2(); newCI=newConnI->getPointer();
1834   //
1835   for(int i=0;i<nbCells;i++,conni++)
1836     {
1837       int sz=conni[1]-conni[0];
1838       int newp=n2oPtr[i];
1839       std::copy(conn+conni[0],conn+conni[1],newC+newCI[newp]);
1840     }
1841   _conn=newConn;
1842 }
1843
1844 MEDCouplingMesh *MEDCoupling1DGTUMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
1845 {
1846   if(other->getType()!=SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED)
1847     throw INTERP_KERNEL::Exception("Merge of umesh only available with umesh single dynamic geo type each other !");
1848   const MEDCoupling1DGTUMesh *otherC=static_cast<const MEDCoupling1DGTUMesh *>(other);
1849   return Merge1DGTUMeshes(this,otherC);
1850 }
1851
1852 MEDCouplingUMesh *MEDCoupling1DGTUMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
1853 {
1854   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(getName(),getMeshDimension());
1855   ret->setCoords(getCoords());
1856   const int *nodalConn=_conn->begin(),*nodalConnI=_conn_indx->begin();
1857   int nbCells=getNumberOfCells();//checkCoherency
1858   int geoType=(int)getCellModelEnum();
1859   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c=DataArrayInt::New(); c->alloc(nbCells+_conn->getNumberOfTuples(),1);
1860   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cI=DataArrayInt::New(); cI->alloc(nbCells+1);
1861   int *cPtr=c->getPointer(),*ciPtr=cI->getPointer();
1862   ciPtr[0]=0;
1863   for(int i=0;i<nbCells;i++,ciPtr++)
1864     {
1865       int sz=nodalConnI[i+1]-nodalConnI[i];
1866       if(sz>=0)
1867         {
1868           *cPtr++=geoType;
1869           cPtr=std::copy(nodalConn+nodalConnI[i],nodalConn+nodalConnI[i+1],cPtr);
1870           ciPtr[1]=ciPtr[0]+sz+1;
1871         }
1872       else
1873         {
1874           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::buildUnstructured : Invalid for nodal index for cell #" << i << " !";
1875           throw INTERP_KERNEL::Exception(oss.str().c_str());
1876         }
1877     }
1878   ret->setConnectivity(c,cI,true);
1879   return ret.retn();
1880 }
1881
1882 /*!
1883  * Do nothing for the moment, because there is no policy that allows to split polygons, polyhedrons ... into simplexes
1884  */
1885 DataArrayInt *MEDCoupling1DGTUMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
1886 {
1887   int nbOfCells=getNumberOfCells();
1888   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1889   ret->alloc(nbOfCells,1);
1890   ret->iota(0);
1891   return ret.retn();
1892 }
1893
1894 void MEDCoupling1DGTUMesh::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
1895 {
1896   stream << "MEDCoupling1DGTUMesh C++ instance at " << this << ". Type=" << _cm->getRepr() << ". Name : \"" << getName() << "\".";
1897   stream << " Mesh dimension : " << getMeshDimension() << ".";
1898   if(!_coords)
1899     { stream << " No coordinates set !"; return ; }
1900   if(!_coords->isAllocated())
1901     { stream << " Coordinates set but not allocated !"; return ; }
1902   stream << " Space dimension : " << _coords->getNumberOfComponents() << "." << std::endl;
1903   stream << "Number of nodes : " << _coords->getNumberOfTuples() << ".";
1904   bool isOK=true;
1905   try { checkCoherency(); } catch(INTERP_KERNEL::Exception& e)
1906     {
1907       stream << std::endl << "Nodal connectivity NOT set properly !\n";
1908       isOK=false;
1909     }
1910   if(isOK)
1911     stream << std::endl << "Number of cells : " << getNumberOfCells() << ".";
1912 }
1913
1914 void MEDCoupling1DGTUMesh::shallowCopyConnectivityFrom(const MEDCouplingPointSet *other) throw(INTERP_KERNEL::Exception)
1915 {
1916   if(!other)
1917     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::shallowCopyConnectivityFrom : input pointer is null !");
1918   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
1919   if(!otherC)
1920     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::shallowCopyConnectivityFrom : input pointer is not an MEDCoupling1DGTUMesh instance !");
1921   setNodalConnectivity(otherC->getNodalConnectivity(),otherC->getNodalConnectivityIndex());
1922 }
1923
1924 MEDCouplingPointSet *MEDCoupling1DGTUMesh::mergeMyselfWithOnSameCoords(const MEDCouplingPointSet *other) const
1925 {
1926   if(!other)
1927     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::mergeMyselfWithOnSameCoords : input other is null !");
1928   const MEDCoupling1DGTUMesh *otherC=dynamic_cast<const MEDCoupling1DGTUMesh *>(other);
1929   if(!otherC)
1930     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::mergeMyselfWithOnSameCoords : the input other mesh is not of type single statuc geo type unstructured !");
1931   std::vector<const MEDCoupling1DGTUMesh *> ms(2);
1932   ms[0]=this;
1933   ms[1]=otherC;
1934   return Merge1DGTUMeshesOnSameCoords(ms);
1935 }
1936
1937 MEDCouplingPointSet *MEDCoupling1DGTUMesh::buildPartOfMySelfKeepCoords(const int *begin, const int *end) const
1938 {
1939   checkCoherency();
1940   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
1941   ret->setCoords(_coords);
1942   DataArrayInt *c=0,*ci=0;
1943   MEDCouplingUMesh::ExtractFromIndexedArrays(begin,end,_conn,_conn_indx,c,ci);
1944   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cSafe(c),ciSafe(ci);
1945   ret->setNodalConnectivity(c,ci);
1946   return ret.retn();
1947 }
1948
1949 MEDCouplingPointSet *MEDCoupling1DGTUMesh::buildPartOfMySelfKeepCoords2(int start, int end, int step) const
1950 {
1951   checkCoherency();
1952   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
1953   ret->setCoords(_coords);
1954   DataArrayInt *c=0,*ci=0;
1955   MEDCouplingUMesh::ExtractFromIndexedArrays2(start,end,step,_conn,_conn_indx,c,ci);
1956   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cSafe(c),ciSafe(ci);
1957   ret->setNodalConnectivity(c,ci);
1958   return ret.retn();
1959 }
1960
1961 void MEDCoupling1DGTUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const throw(INTERP_KERNEL::Exception)
1962 {
1963   checkFullyDefined();
1964   int nbOfNodes=getNumberOfNodes();
1965   int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int));
1966   revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1);
1967   std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
1968   const int *conn=_conn->begin(),*conni=_conn_indx->begin();
1969   int nbOfCells=getNumberOfCells();
1970   int nbOfEltsInRevNodal=0;
1971   for(int eltId=0;eltId<nbOfCells;eltId++)
1972     {
1973       int nbOfNodesPerCell=conni[eltId+1]-conni[eltId];
1974       if(nbOfNodesPerCell>=0)
1975         {
1976           for(int j=0;j<nbOfNodesPerCell;j++)
1977             {
1978               int nodeId=conn[conni[eltId]+j];
1979               if(nodeId==-1) continue;
1980               if(nodeId>=0 && nodeId<nbOfNodes)
1981                 {
1982                   nbOfEltsInRevNodal++;
1983                   revNodalIndxPtr[nodeId+1]++;
1984                 }
1985               else
1986                 {
1987                   std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getReverseNodalConnectivity : At cell #" << eltId << " presence of nodeId #" << conn[0] << " should be in [0," << nbOfNodes << ") !";
1988                   throw INTERP_KERNEL::Exception(oss.str().c_str());
1989                 }
1990             }
1991         }
1992       else
1993         {
1994           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getReverseNodalConnectivity : At cell #" << eltId << "nodal connectivity is invalid !";
1995           throw INTERP_KERNEL::Exception(oss.str().c_str());
1996         }
1997     }
1998   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<int>());
1999   conn=_conn->begin();
2000   int *revNodalPtr=(int *)malloc((nbOfEltsInRevNodal)*sizeof(int));
2001   revNodal->useArray(revNodalPtr,true,C_DEALLOC,nbOfEltsInRevNodal,1);
2002   std::fill(revNodalPtr,revNodalPtr+nbOfEltsInRevNodal,-1);
2003   for(int eltId=0;eltId<nbOfCells;eltId++)
2004     {
2005       int nbOfNodesPerCell=conni[eltId+1]-conni[eltId];
2006       for(int j=0;j<nbOfNodesPerCell;j++)
2007         {
2008           int nodeId=conn[conni[eltId]+j];
2009           if(nodeId!=-1)
2010             *std::find_if(revNodalPtr+revNodalIndxPtr[nodeId],revNodalPtr+revNodalIndxPtr[nodeId+1],std::bind2nd(std::equal_to<int>(),-1))=eltId;
2011         }
2012     }
2013 }
2014
2015 void MEDCoupling1DGTUMesh::checkFullyDefined() const throw(INTERP_KERNEL::Exception)
2016 {
2017   if(!((const DataArrayInt *)_conn) || !((const DataArrayInt *)_conn_indx) || !((const DataArrayDouble *)_coords))
2018     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::checkFullyDefined : part of this is not fully defined.");
2019 }
2020
2021 bool MEDCoupling1DGTUMesh::isEmptyMesh(const std::vector<int>& tinyInfo) const
2022 {
2023   throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::isEmptyMesh : not implemented yet !");
2024 }
2025
2026 /*!
2027  * Finds nodes not used in any cell and returns an array giving a new id to every node
2028  * by excluding the unused nodes, for which the array holds -1. The result array is
2029  * a mapping in "Old to New" mode.
2030  *  \param [out] nbrOfNodesInUse - number of node ids present in the nodal connectivity.
2031  *  \return DataArrayInt * - a new instance of DataArrayInt. Its length is \a
2032  *          this->getNumberOfNodes(). It holds for each node of \a this mesh either -1
2033  *          if the node is unused or a new id else. The caller is to delete this
2034  *          array using decrRef() as it is no more needed.
2035  *  \throw If the coordinates array is not set.
2036  *  \throw If the nodal connectivity of cells is not defined.
2037  *  \throw If the nodal connectivity includes an invalid id.
2038  */
2039 DataArrayInt *MEDCoupling1DGTUMesh::getNodeIdsInUse(int& nbrOfNodesInUse) const throw(INTERP_KERNEL::Exception)
2040 {
2041   nbrOfNodesInUse=-1;
2042   int nbOfNodes=getNumberOfNodes();
2043   int nbOfCells=getNumberOfCells();//checkCoherency
2044   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
2045   ret->alloc(nbOfNodes,1);
2046   int *traducer=ret->getPointer();
2047   std::fill(traducer,traducer+nbOfNodes,-1);
2048   const int *conn=_conn->begin(),*conni(_conn_indx->begin());
2049   for(int i=0;i<nbOfCells;i++,conni++)
2050     {
2051       int nbNodesPerCell=conni[1]-conni[0];
2052       for(int j=0;j<nbNodesPerCell;j++)
2053         {
2054           int nodeId=conn[conni[0]+j];
2055           if(nodeId==-1) continue;
2056           if(nodeId>=0 && nodeId<nbOfNodes)
2057             traducer[*conn]=1;
2058           else
2059             {
2060               std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::getNodeIdsInUse : In cell #" << i  << " presence of node id " <<  nodeId << " not in [0," << nbOfNodes << ") !";
2061               throw INTERP_KERNEL::Exception(oss.str().c_str());
2062             }
2063         }
2064     }
2065   nbrOfNodesInUse=(int)std::count(traducer,traducer+nbOfNodes,1);
2066   std::transform(traducer,traducer+nbOfNodes,traducer,MEDCouplingAccVisit());
2067   return ret.retn();
2068 }
2069
2070 /*!
2071  * Changes ids of nodes within the nodal connectivity arrays according to a permutation
2072  * array in "Old to New" mode. The node coordinates array is \b not changed by this method.
2073  * This method is a generalization of shiftNodeNumbersInConn().
2074  *  \warning This method performs no check of validity of new ids. **Use it with care !**
2075  *  \param [in] newNodeNumbersO2N - a permutation array, of length \a
2076  *         this->getNumberOfNodes(), in "Old to New" mode.
2077  *         See \ref MEDCouplingArrayRenumbering for more info on renumbering modes.
2078  *  \throw If the nodal connectivity of cells is not defined.
2079  */
2080 void MEDCoupling1DGTUMesh::renumberNodesInConn(const int *newNodeNumbersO2N)
2081 {
2082   getNumberOfCells();//only to check that all is well defined.
2083   //
2084   int nbElemsIn=getNumberOfNodes();
2085   int nbOfTuples=_conn->getNumberOfTuples();
2086   int *pt=_conn->getPointer();
2087   for(int i=0;i<nbOfTuples;i++,pt++)
2088     {
2089       if(*pt==-1) continue;
2090       if(*pt>=0 && *pt<nbElemsIn)
2091         *pt=newNodeNumbersO2N[*pt];
2092       else
2093         {
2094           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::renumberNodesInConn : error on tuple #" << i << " value is " << *pt << " and indirectionnal array as a size equal to " << nbElemsIn;
2095           throw INTERP_KERNEL::Exception(oss.str().c_str());
2096         }
2097     }
2098   _conn->declareAsNew();
2099   //
2100   updateTime();
2101 }
2102
2103 /*!
2104  * Keeps from \a this only cells which constituing point id are in the ids specified by [\a begin,\a end).
2105  * The resulting cell ids are stored at the end of the 'cellIdsKept' parameter.
2106  * Parameter \a fullyIn specifies if a cell that has part of its nodes in ids array is kept or not.
2107  * If \a fullyIn is true only cells whose ids are \b fully contained in [\a begin,\a end) tab will be kept.
2108  *
2109  * \param [in] begin input start of array of node ids.
2110  * \param [in] end input end of array of node ids.
2111  * \param [in] fullyIn input that specifies if all node ids must be in [\a begin,\a end) array to consider cell to be in.
2112  * \param [in,out] cellIdsKeptArr array where all candidate cell ids are put at the end.
2113  */
2114 void MEDCoupling1DGTUMesh::fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const
2115 {
2116   int nbOfCells=getNumberOfCells();
2117   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIdsKept=DataArrayInt::New(); cellIdsKept->alloc(0,1);
2118   int tmp=-1;
2119   int sz=_conn->getMaxValue(tmp); sz=std::max(sz,0)+1;
2120   std::vector<bool> fastFinder(sz,false);
2121   for(const int *work=begin;work!=end;work++)
2122     if(*work>=0 && *work<sz)
2123       fastFinder[*work]=true;
2124   const int *conn=_conn->begin(),*conni=_conn_indx->begin();
2125   for(int i=0;i<nbOfCells;i++,conni++)
2126     {
2127       int ref=0,nbOfHit=0;
2128       int nbNodesPerCell=conni[1]-conni[0];
2129       if(nbNodesPerCell>=0)
2130         {
2131           for(int j=0;j<nbNodesPerCell;j++)
2132             {
2133               int nodeId=conn[conni[0]+j];
2134               if(nodeId>=0)
2135                 {
2136                   ref++;
2137                   if(fastFinder[nodeId])
2138                     nbOfHit++;
2139                 }
2140             }
2141         }
2142       else
2143         {
2144           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::fillCellIdsToKeepFromNodeIds : invalid index array for cell #" << i << " !";
2145           throw INTERP_KERNEL::Exception(oss.str().c_str());
2146         }
2147       if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
2148         cellIdsKept->pushBackSilent(i);
2149     }
2150   cellIdsKeptArr=cellIdsKept.retn();
2151 }
2152
2153 void MEDCoupling1DGTUMesh::allocateCells(int nbOfCells) throw(INTERP_KERNEL::Exception)
2154 {
2155   if(nbOfCells<0)
2156     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::allocateCells : the input number of cells should be >= 0 !");
2157   _conn=DataArrayInt::New();
2158   _conn->reserve(nbOfCells*3);
2159   _conn_indx=DataArrayInt::New();
2160   _conn_indx->reserve(nbOfCells+1); _conn_indx->pushBackSilent(0);
2161   declareAsNew();
2162 }
2163
2164 /*!
2165  * Appends at the end of \a this a cell having nodal connectivity array defined in [ \a nodalConnOfCellBg, \a nodalConnOfCellEnd ).
2166  *
2167  * \param [in] nodalConnOfCellBg - the begin (included) of nodal connectivity of the cell to add.
2168  * \param [in] nodalConnOfCellEnd - the end (excluded) of nodal connectivity of the cell to add.
2169  * \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
2170  *        attached to \a this.
2171  * \thow If the nodal connectivity array in \a this is null (call MEDCoupling1SGTUMesh::allocateCells before).
2172  */
2173 void MEDCoupling1DGTUMesh::insertNextCell(const int *nodalConnOfCellBg, const int *nodalConnOfCellEnd) throw(INTERP_KERNEL::Exception)
2174 {
2175   int sz=(int)std::distance(nodalConnOfCellBg,nodalConnOfCellEnd);
2176   DataArrayInt *c(_conn),*c2(_conn_indx);
2177   if(c && c2)
2178     {
2179       int pos=c2->back();
2180       if(pos==c->getNumberOfTuples())
2181         {
2182           c->pushBackValsSilent(nodalConnOfCellBg,nodalConnOfCellEnd);
2183           c2->pushBackSilent(pos+sz);
2184         }
2185       else
2186         {
2187           std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::insertNextCell : The nodal index array (end=" << pos << ") mismatches with nodal array (length=" << c->getNumberOfTuples() << ") !";
2188           throw INTERP_KERNEL::Exception(oss.str().c_str());
2189         }
2190     }
2191   else
2192     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::insertNextCell : nodal connectivity array is null ! Call MEDCoupling1DGTUMesh::allocateCells before !");
2193 }
2194
2195 void MEDCoupling1DGTUMesh::setNodalConnectivity(DataArrayInt *nodalConn, DataArrayInt *nodalConnIndex) throw(INTERP_KERNEL::Exception)
2196 {
2197   if(nodalConn)
2198     nodalConn->incrRef();
2199   _conn=nodalConn;
2200   if(nodalConnIndex)
2201     nodalConnIndex->incrRef();
2202   _conn_indx=nodalConnIndex;
2203   declareAsNew();
2204 }
2205
2206 /*!
2207  * \return DataArrayInt * - the internal reference to the nodal connectivity. The caller is not reponsible to deallocate it.
2208  */
2209 DataArrayInt *MEDCoupling1DGTUMesh::getNodalConnectivity() const throw(INTERP_KERNEL::Exception)
2210 {
2211   const DataArrayInt *ret(_conn);
2212   return const_cast<DataArrayInt *>(ret);
2213 }
2214
2215 /*!
2216  * \return DataArrayInt * - the internal reference to the nodal connectivity index. The caller is not reponsible to deallocate it.
2217  */
2218 DataArrayInt *MEDCoupling1DGTUMesh::getNodalConnectivityIndex() const throw(INTERP_KERNEL::Exception)
2219 {
2220   const DataArrayInt *ret(_conn_indx);
2221   return const_cast<DataArrayInt *>(ret);
2222 }
2223
2224 /*!
2225  * See the definition of the nodal connectivity pack \ref MEDCoupling1DGTUMesh::isPacked "here".
2226  * This method tries to build a new instance geometrically equivalent to \a this, by limiting at most the number of new object (nodal connectivity).
2227  * 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.
2228  *
2229  * 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.
2230  *
2231  * \param [out] isShallowCpyOfNodalConnn - tells if the returned instance share the same pair of nodal connectivity arrays (true) or if nodal
2232  *              connectivity arrays are different (false)
2233  * \return a new object to be managed by the caller.
2234  *
2235  * \sa MEDCoupling1DGTUMesh::retrievePackedNodalConnectivity, MEDCoupling1DGTUMesh::isPacked
2236  */
2237 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::copyWithNodalConnectivityPacked(bool& isShallowCpyOfNodalConnn) const throw(INTERP_KERNEL::Exception)
2238 {
2239   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
2240   DataArrayInt *nc=0,*nci=0;
2241   isShallowCpyOfNodalConnn=retrievePackedNodalConnectivity(nc,nci);
2242   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ncs(nc),ncis(nci);
2243   ret->_conn=ncs; ret->_conn_indx=ncis;
2244   ret->setCoords(getCoords());
2245   return ret.retn();
2246 }
2247
2248 /*!
2249  * This method allows to compute, if needed, the packed nodal connectivity pair.
2250  * Indeed, it is possible to store in \a this a nodal connectivity array bigger than ranges convered by nodal connectivity index array.
2251  * 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.
2252  *
2253  * 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)
2254  * true will be returned and respectively \a this->_conn and \a this->_conn_indx (with ref counter incremented). This is the classical case.
2255  *
2256  * 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
2257  * will be returned.
2258  *
2259  * This method return 3 elements.
2260  * \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
2261  *                          this pointer can be seen as a new object, that is to managed by the caller.
2262  * \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
2263  *                              this pointer can be seen as a new object, that is to managed by the caller.
2264  * \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
2265  * output parameters are newly created objects.
2266  *
2267  * \throw if \a this does not pass MEDCoupling1DGTUMesh::checkCoherency test
2268  */
2269 bool MEDCoupling1DGTUMesh::retrievePackedNodalConnectivity(DataArrayInt *&nodalConn, DataArrayInt *&nodalConnIndx) const throw(INTERP_KERNEL::Exception)
2270 {
2271   if(isPacked())//performs the checkCoherency
2272     {
2273       const DataArrayInt *c0(_conn),*c1(_conn_indx);
2274       nodalConn=const_cast<DataArrayInt *>(c0); nodalConnIndx=const_cast<DataArrayInt *>(c1);
2275       nodalConn->incrRef(); nodalConnIndx->incrRef();
2276       return true;
2277     }
2278   int bg=_conn_indx->front(),end=_conn_indx->back();
2279   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nc(_conn->selectByTupleId2(bg,end,1));
2280   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nci(_conn_indx->deepCpy());
2281   nci->applyLin(1,-bg);
2282   nodalConn=nc.retn(); nodalConnIndx=nci.retn();
2283   return false;
2284 }
2285
2286 /*
2287  * 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)
2288  * true will be returned and respectively \a this->_conn and \a this->_conn_indx (with ref counter incremented). This is the classical case.
2289  * If nodal connectivity index points to a subpart of nodal connectivity index false will be returned.
2290  * \return bool - true if \a this looks packed, false is not.
2291  *
2292  * \throw if \a this does not pass MEDCoupling1DGTUMesh::checkCoherency test
2293  */
2294 bool MEDCoupling1DGTUMesh::isPacked() const throw(INTERP_KERNEL::Exception)
2295 {
2296   checkCoherency();
2297   return _conn_indx->front()==0 && _conn_indx->back()==_conn->getNumberOfTuples();
2298 }
2299
2300 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshes(const MEDCoupling1DGTUMesh *mesh1, const MEDCoupling1DGTUMesh *mesh2) throw(INTERP_KERNEL::Exception)
2301 {
2302   std::vector<const MEDCoupling1DGTUMesh *> tmp(2);
2303   tmp[0]=const_cast<MEDCoupling1DGTUMesh *>(mesh1); tmp[1]=const_cast<MEDCoupling1DGTUMesh *>(mesh2);
2304   return Merge1DGTUMeshes(tmp);
2305 }
2306
2307 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshes(std::vector<const MEDCoupling1DGTUMesh *>& a) throw(INTERP_KERNEL::Exception)
2308 {
2309   std::size_t sz=a.size();
2310   if(sz==0)
2311     return Merge1DGTUMeshesLL(a);
2312   for(std::size_t ii=0;ii<sz;ii++)
2313     if(!a[ii])
2314       {
2315         std::ostringstream oss; oss << "MEDCoupling1DGTUMesh::Merge1DGTUMeshes : item #" << ii << " in input array of size "<< sz << " is empty !";
2316         throw INTERP_KERNEL::Exception(oss.str().c_str());
2317       }
2318   const INTERP_KERNEL::CellModel *cm=&(a[0]->getCellModel());
2319   for(std::size_t ii=0;ii<sz;ii++)
2320     if(&(a[ii]->getCellModel())!=cm)
2321       throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshes : all items must have the same geo type !");
2322   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> > bb(sz);
2323   std::vector< const MEDCoupling1DGTUMesh * > aa(sz);
2324   int spaceDim=-3;
2325   for(std::size_t i=0;i<sz && spaceDim==-3;i++)
2326     {
2327       const MEDCoupling1DGTUMesh *cur=a[i];
2328       const DataArrayDouble *coo=cur->getCoords();
2329       if(coo)
2330         spaceDim=coo->getNumberOfComponents();
2331     }
2332   if(spaceDim==-3)
2333     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshes : no spaceDim specified ! unable to perform merge !");
2334   for(std::size_t i=0;i<sz;i++)
2335     {
2336       bb[i]=a[i]->buildSetInstanceFromThis(spaceDim);
2337       aa[i]=bb[i];
2338     }
2339   return Merge1DGTUMeshesLL(aa);
2340 }
2341
2342 /*!
2343  * \throw If presence of a null instance in the input vector \a a.
2344  * \throw If a is empty
2345  */
2346 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords(std::vector<const MEDCoupling1DGTUMesh *>& a) throw(INTERP_KERNEL::Exception)
2347 {
2348   if(a.empty())
2349     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : input array must be NON EMPTY !");
2350   std::vector<const MEDCoupling1DGTUMesh *>::const_iterator it=a.begin();
2351   if(!(*it))
2352     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : null instance in the first element of input vector !");
2353   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> > objs(a.size());
2354   std::vector<const DataArrayInt *> ncs(a.size()),ncis(a.size());
2355   int nbOfCells=(*it)->getNumberOfCells();
2356   const DataArrayDouble *coords=(*it)->getCoords();
2357   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
2358   bool tmp;
2359   objs[0]=(*it)->copyWithNodalConnectivityPacked(tmp);
2360   ncs[0]=objs[0]->getNodalConnectivity(); ncis[0]=objs[0]->getNodalConnectivityIndex();
2361   it++;
2362   for(int i=1;it!=a.end();i++,it++)
2363     {
2364       if(!(*it))
2365         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : presence of null instance !");
2366       if(cm!=&((*it)->getCellModel()))
2367         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1DGTUMeshes impossible !");
2368       (*it)->getNumberOfCells();//to check that all is OK
2369       objs[i]=(*it)->copyWithNodalConnectivityPacked(tmp);
2370       ncs[i]=objs[i]->getNodalConnectivity(); ncis[i]=objs[i]->getNodalConnectivityIndex();
2371       if(coords!=(*it)->getCoords())
2372         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshesOnSameCoords : not lying on same coords !");
2373     }
2374   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh("merge",*cm));
2375   ret->setCoords(coords);
2376   ret->_conn=DataArrayInt::Aggregate(ncs);
2377   ret->_conn_indx=DataArrayInt::AggregateIndexes(ncis);
2378   return ret.retn();
2379 }
2380
2381 /*!
2382  * 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)
2383  */
2384 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::Merge1DGTUMeshesLL(std::vector<const MEDCoupling1DGTUMesh *>& a) throw(INTERP_KERNEL::Exception)
2385 {
2386   if(a.empty())
2387     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::Merge1DGTUMeshes : input array must be NON EMPTY !");
2388   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> > objs(a.size());
2389   std::vector<const DataArrayInt *> ncs(a.size()),ncis(a.size());
2390   std::vector<const MEDCoupling1DGTUMesh *>::const_iterator it=a.begin();
2391   std::vector<int> nbNodesPerElt(a.size());
2392   int nbOfCells=(*it)->getNumberOfCells();
2393   bool tmp;
2394   objs[0]=(*it)->copyWithNodalConnectivityPacked(tmp);
2395   ncs[0]=objs[0]->getNodalConnectivity(); ncis[0]=objs[0]->getNodalConnectivityIndex();
2396   nbNodesPerElt[0]=0;
2397   int prevNbOfNodes=(*it)->getNumberOfNodes();
2398   const INTERP_KERNEL::CellModel *cm=&((*it)->getCellModel());
2399   it++;
2400   for(int i=1;it!=a.end();i++,it++)
2401     {
2402       if(cm!=&((*it)->getCellModel()))
2403         throw INTERP_KERNEL::Exception("Geometric types mismatches, Merge1DGTUMeshes impossible !");
2404       objs[i]=(*it)->copyWithNodalConnectivityPacked(tmp);
2405       ncs[i]=objs[i]->getNodalConnectivity(); ncis[i]=objs[i]->getNodalConnectivityIndex();
2406       nbOfCells+=(*it)->getNumberOfCells();
2407       nbNodesPerElt[i]=nbNodesPerElt[i-1]+prevNbOfNodes;
2408       prevNbOfNodes=(*it)->getNumberOfNodes();
2409     }
2410   std::vector<const MEDCouplingPointSet *> aps(a.size());
2411   std::copy(a.begin(),a.end(),aps.begin());
2412   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> pts=MergeNodesArray(aps);
2413   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh("merge",*cm));
2414   ret->setCoords(pts);
2415   ret->_conn=AggregateNodalConnAndShiftNodeIds(ncs,nbNodesPerElt);
2416   ret->_conn_indx=DataArrayInt::AggregateIndexes(ncis);
2417   return ret.retn();
2418 }
2419
2420 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::buildSetInstanceFromThis(int spaceDim) const throw(INTERP_KERNEL::Exception)
2421 {
2422   MEDCouplingAutoRefCountObjectPtr<MEDCoupling1DGTUMesh> ret(new MEDCoupling1DGTUMesh(getName(),*_cm));
2423   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1,tmp2;
2424   const DataArrayInt *nodalConn(_conn),*nodalConnI(_conn_indx);
2425   if(!nodalConn)
2426     {
2427       tmp1=DataArrayInt::New(); tmp1->alloc(0,1);
2428     }
2429   else
2430     tmp1=_conn;
2431   ret->_conn=tmp1;
2432   //
2433   if(!nodalConnI)
2434     {
2435       tmp2=DataArrayInt::New(); tmp2->alloc(1,1); tmp2->setIJ(0,0,0);
2436     }
2437   else
2438     tmp2=_conn_indx;
2439   ret->_conn_indx=tmp2;
2440   //
2441   if(!_coords)
2442     {
2443       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=DataArrayDouble::New(); coords->alloc(0,spaceDim);
2444       ret->setCoords(coords);
2445     }
2446   else
2447     ret->setCoords(_coords);
2448   return ret.retn();
2449 }
2450
2451 /*!
2452  * This method performs an aggregation of \a nodalConns (as DataArrayInt::Aggregate does) but in addition of that a shift is applied on the
2453  * values contained in \a nodalConns using corresponding offset specified in input \a offsetInNodeIdsPerElt.
2454  * But it also manage the values -1, that have a semantic in MEDCoupling1DGTUMesh class (separator for polyhedron).
2455  *
2456  * \param [in] nodalConns - a list of nodal connectivity arrays same size than \a offsetInNodeIdsPerElt.
2457  * \param [in] offsetInNodeIdsPerElt - a list of offsets to apply.
2458  * \return DataArrayInt * - A new object (to be managed by the caller) that is the result of the aggregation.
2459  * \throw If \a nodalConns or \a offsetInNodeIdsPerElt are empty.
2460  * \throw If \a nodalConns and \a offsetInNodeIdsPerElt have not the same size.
2461  * \throw If presence of null pointer in \a nodalConns.
2462  * \throw If presence of not allocated or array with not exactly one component in \a nodalConns.
2463  */
2464 DataArrayInt *MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds(const std::vector<const DataArrayInt *>& nodalConns, const std::vector<int>& offsetInNodeIdsPerElt) throw(INTERP_KERNEL::Exception)
2465 {
2466   std::size_t sz1(nodalConns.size()),sz2(offsetInNodeIdsPerElt.size());
2467   if(sz1!=sz2)
2468     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : input vectors do not have the same size !");
2469   if(sz1==0)
2470     throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : empty vectors in input !");
2471   int nbOfTuples=0;
2472   for(std::vector<const DataArrayInt *>::const_iterator it=nodalConns.begin();it!=nodalConns.end();it++)
2473     {
2474       if(!(*it))
2475         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : presence of null pointer in input vector !");
2476       if(!(*it)->isAllocated())
2477         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : presence of non allocated array in input vector !");
2478       if((*it)->getNumberOfComponents()!=1)
2479         throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::AggregateNodalConnAndShiftNodeIds : presence of array with not exactly one component !");
2480       nbOfTuples+=(*it)->getNumberOfTuples();
2481     }
2482   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(nbOfTuples,1);
2483   int *pt=ret->getPointer();
2484   int i=0;
2485   for(std::vector<const DataArrayInt *>::const_iterator it=nodalConns.begin();it!=nodalConns.end();it++,i++)
2486     {
2487       int curNbt=(*it)->getNumberOfTuples();
2488       const int *inPt=(*it)->begin();
2489       int offset=offsetInNodeIdsPerElt[i];
2490       for(int j=0;j<curNbt;j++,pt++)
2491         {
2492           if(inPt[j]!=-1)
2493             *pt=inPt[j]+offset;
2494           else
2495             *pt=-1;
2496         }
2497     }
2498   return ret.retn();
2499 }