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