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