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