]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/MEDCoupling/MEDCouplingCMesh.cxx
Salome HOME
Little refactoring of progeny mechanism to avoid if.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingCMesh.cxx
1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay (CEA/DEN)
20
21 #include "MEDCouplingCMesh.hxx"
22 #include "MEDCouplingMemArray.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24
25 #include <functional>
26 #include <algorithm>
27 #include <sstream>
28 #include <numeric>
29
30 using namespace ParaMEDMEM;
31
32 MEDCouplingCMesh::MEDCouplingCMesh():_x_array(0),_y_array(0),_z_array(0)
33 {
34 }
35
36 MEDCouplingCMesh::MEDCouplingCMesh(const MEDCouplingCMesh& other, bool deepCopy):MEDCouplingStructuredMesh(other,deepCopy)
37 {
38   if(deepCopy)
39     {
40       if(other._x_array)
41         _x_array=other._x_array->deepCpy();
42       else
43         _x_array=0;
44       if(other._y_array)
45         _y_array=other._y_array->deepCpy();
46       else
47         _y_array=0;
48       if(other._z_array)
49         _z_array=other._z_array->deepCpy();
50       else
51         _z_array=0;
52     }
53   else
54     {
55       _x_array=other._x_array;
56       if(_x_array)
57         _x_array->incrRef();
58       _y_array=other._y_array;
59       if(_y_array)
60         _y_array->incrRef();
61       _z_array=other._z_array;
62       if(_z_array)
63         _z_array->incrRef();
64     }
65 }
66
67 MEDCouplingCMesh::~MEDCouplingCMesh()
68 {
69   if(_x_array)
70     _x_array->decrRef();
71   if(_y_array)
72     _y_array->decrRef();
73   if(_z_array)
74     _z_array->decrRef();
75 }
76
77 MEDCouplingCMesh *MEDCouplingCMesh::New()
78 {
79   return new MEDCouplingCMesh;
80 }
81
82 MEDCouplingCMesh *MEDCouplingCMesh::New(const std::string& meshName)
83 {
84   MEDCouplingCMesh *ret=new MEDCouplingCMesh;
85   ret->setName(meshName);
86   return ret;
87 }
88
89 MEDCouplingMesh *MEDCouplingCMesh::deepCpy() const
90 {
91   return clone(true);
92 }
93
94 MEDCouplingCMesh *MEDCouplingCMesh::clone(bool recDeepCpy) const
95 {
96   return new MEDCouplingCMesh(*this,recDeepCpy);
97 }
98
99 void MEDCouplingCMesh::updateTime() const
100 {
101   if(_x_array)
102     updateTimeWith(*_x_array);
103   if(_y_array)
104     updateTimeWith(*_y_array);
105   if(_z_array)
106     updateTimeWith(*_z_array);
107 }
108
109 std::size_t MEDCouplingCMesh::getHeapMemorySizeWithoutChildren() const
110 {
111   return MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren();
112 }
113
114 std::vector<const BigMemoryObject *> MEDCouplingCMesh::getDirectChildrenWithNull() const
115 {
116   std::vector<const BigMemoryObject *> ret;
117   ret.push_back(_x_array);
118   ret.push_back(_y_array);
119   ret.push_back(_z_array);
120   return ret;
121 }
122
123 /*!
124  * This method copyies all tiny strings from other (name and components name).
125  * @throw if other and this have not same mesh type.
126  */
127 void MEDCouplingCMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
128 {
129   MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
130   const MEDCouplingCMesh *otherC(dynamic_cast<const MEDCouplingCMesh *>(other));
131   if(!otherC)
132     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::copyTinyStringsFrom : meshes have not same type !");
133   if(_x_array && otherC->_x_array)
134     _x_array->copyStringInfoFrom(*otherC->_x_array);
135   if(_y_array && otherC->_y_array)
136     _y_array->copyStringInfoFrom(*otherC->_y_array);
137   if(_z_array && otherC->_z_array)
138     _z_array->copyStringInfoFrom(*otherC->_z_array);
139 }
140
141 bool MEDCouplingCMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
142 {
143   if(!other)
144     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::isEqualIfNotWhy : input other pointer is null !");
145   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
146   if(!otherC)
147     {
148       reason="mesh given in input is not castable in MEDCouplingCMesh !";
149       return false;
150     }
151   if(!MEDCouplingStructuredMesh::isEqualIfNotWhy(other,prec,reason))
152     return false;
153   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
154   const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
155   std::ostringstream oss; oss.precision(15);
156   for(int i=0;i<3;i++)
157     {
158       if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
159         {
160           oss << "Only one CMesh between the two this and other has its coordinates of rank" << i << " defined !";
161           reason=oss.str();
162           return false;
163         }
164       if(thisArr[i])
165         if(!thisArr[i]->isEqualIfNotWhy(*otherArr[i],prec,reason))
166           {
167             oss << "Coordinates DataArrayDouble of rank #" << i << " differ :";
168             reason.insert(0,oss.str());
169             return false;
170           }
171     }
172   return true;
173 }
174
175 bool MEDCouplingCMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
176 {
177   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
178   if(!otherC)
179     return false;
180   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
181   const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
182   for(int i=0;i<3;i++)
183     {
184       if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
185         return false;
186       if(thisArr[i])
187         if(!thisArr[i]->isEqualWithoutConsideringStr(*otherArr[i],prec))
188           return false;
189     }
190   return true;
191 }
192
193 void MEDCouplingCMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
194                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const
195 {
196   if(!isEqualWithoutConsideringStr(other,prec))
197     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalWith : Meshes are not the same !");
198 }
199
200 /*!
201  * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingCMesh instance too).
202  * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingCMesh, \a this and \a other are the same !
203  */
204 void MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
205                                                        DataArrayInt *&cellCor) const
206 {
207   if(!isEqualWithoutConsideringStr(other,prec))
208     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith : Meshes are not the same !");
209 }
210
211 void MEDCouplingCMesh::checkCoherency() const
212 {
213   const char msg0[]="Invalid ";
214   const char msg1[]=" array ! Must contain more than 1 element.";
215   const char msg2[]=" array ! Must be with only one component.";
216   getSpaceDimension();// here to check that no holes in arrays !
217   if(_x_array)
218     {
219       if(_x_array->getNbOfElems()<2)
220         {
221           std::ostringstream os; os << msg0 << 'X' << msg1;
222           throw INTERP_KERNEL::Exception(os.str().c_str());
223         }
224       if(_x_array->getNumberOfComponents()!=1)
225         {
226           std::ostringstream os; os << msg0 << 'X' << msg2;
227           throw INTERP_KERNEL::Exception(os.str().c_str());
228         }
229     }
230   if(_y_array)
231     {
232       if(_y_array->getNbOfElems()<2)
233         {
234           std::ostringstream os; os << msg0 << 'Y' << msg1;
235           throw INTERP_KERNEL::Exception(os.str().c_str());
236         }
237       if(_y_array->getNumberOfComponents()!=1)
238         {
239           std::ostringstream os; os << msg0 << 'Y' << msg2;
240           throw INTERP_KERNEL::Exception(os.str().c_str());
241         }
242     }
243   if(_z_array)
244     {
245       if(_z_array->getNbOfElems()<2)
246         {
247           std::ostringstream os; os << msg0 << 'Z' << msg1;
248           throw INTERP_KERNEL::Exception(os.str().c_str());
249         }
250       if(_z_array->getNumberOfComponents()!=1)
251         {
252           std::ostringstream os; os << msg0 << 'Z' << msg2;
253           throw INTERP_KERNEL::Exception(os.str().c_str());
254         }
255     }
256 }
257
258 void MEDCouplingCMesh::checkCoherency1(double eps) const
259 {
260   checkCoherency();
261   if(_x_array)
262     _x_array->checkMonotonic(true, eps);
263   if(_y_array)
264     _y_array->checkMonotonic(true, eps);
265   if(_z_array)
266     _z_array->checkMonotonic(true, eps);
267 }
268
269 void MEDCouplingCMesh::checkCoherency2(double eps) const
270 {
271   checkCoherency1(eps);
272 }
273
274 void MEDCouplingCMesh::getNodeGridStructure(int *res) const
275 {
276   std::vector<int> ret(getNodeGridStructure());
277   std::copy(ret.begin(),ret.end(),res);
278 }
279
280 std::vector<int> MEDCouplingCMesh::getNodeGridStructure() const
281 {
282   static const char MSG[]="MEDCouplingCMesh::getNodeGridStructure : mesh is invalid ! null vectors (X, Y or Z) must be put contiguously at the end !";
283   std::vector<int> ret;
284   bool isOK(true);
285   if(_x_array)
286     {
287       if(!_x_array->isAllocated() || _x_array->getNumberOfComponents()!=1)
288         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getNodeGridStructure : X array exits but it is not allocated or with nb of components equal to one !");
289       ret.push_back(_x_array->getNumberOfTuples());
290     }
291   else
292     isOK=false;
293   if(_y_array)
294     {
295       if(!_y_array->isAllocated() || _y_array->getNumberOfComponents()!=1)
296         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getNodeGridStructure : Y array exits but it is not allocated or with nb of components equal to one !");
297       if(!isOK)
298         throw INTERP_KERNEL::Exception(MSG);
299       ret.push_back(_y_array->getNumberOfTuples());
300     }
301   else
302     isOK=false;
303   if(_z_array)
304     {
305       if(!_z_array->isAllocated() || _z_array->getNumberOfComponents()!=1)
306         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getNodeGridStructure : Z array exits but it is not allocated or with nb of components equal to one !");
307       if(!isOK)
308         throw INTERP_KERNEL::Exception(MSG);
309       ret.push_back(_z_array->getNumberOfTuples());
310     }
311   return ret;
312 }
313
314 MEDCouplingStructuredMesh *MEDCouplingCMesh::buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const
315 {
316   checkCoherency();
317   int dim(getSpaceDimension());
318   if(dim!=(int)cellPart.size())
319     {
320       std::ostringstream oss; oss << "MEDCouplingCMesh::buildStructuredSubPart : the space dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
321       throw INTERP_KERNEL::Exception(oss.str().c_str());
322     }
323   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCMesh> ret(dynamic_cast<MEDCouplingCMesh *>(deepCpy()));
324   for(int i=0;i<dim;i++)
325     {
326       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp(ret->getCoordsAt(i)->selectByTupleId2(cellPart[i].first,cellPart[i].second+1,1));
327       ret->setCoordsAt(i,tmp);
328     }
329   return ret.retn();
330 }
331
332 /*!
333  * Return the space dimension of \a this. It only considers the arrays along X, Y and Z to deduce that.
334  * This method throws exceptions if the not null arrays defining this are not contiguously at the end. For example X!=0,Y==0,Z!=0 will throw.
335  */
336 int MEDCouplingCMesh::getSpaceDimension() const
337 {
338   return (int)getNodeGridStructure().size();
339 }
340
341 void MEDCouplingCMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
342 {
343   int tmp[3];
344   int spaceDim=getSpaceDimension();
345   getSplitNodeValues(tmp);
346   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
347   int tmp2[3];
348   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
349   for(int j=0;j<spaceDim;j++)
350     if(tabs[j])
351       coo.push_back(tabs[j]->getConstPointer()[tmp2[j]]);
352 }
353
354 std::string MEDCouplingCMesh::simpleRepr() const
355 {
356   std::ostringstream ret;
357   ret << "Cartesian mesh with name : \"" << getName() << "\"\n";
358   ret << "Description of mesh : \"" << getDescription() << "\"\n";
359   int tmpp1,tmpp2;
360   double tt=getTime(tmpp1,tmpp2);
361   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
362   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
363   ret << "Space dimension : " << getSpaceDimension() << "\n\nArrays :\n________\n\n";
364   if(_x_array)
365     {
366       ret << "X Array :\n";
367       _x_array->reprZipWithoutNameStream(ret);
368     }
369   if(_y_array)
370     {
371       ret << "Y Array :\n";
372       _y_array->reprZipWithoutNameStream(ret);
373     }
374   if(_z_array)
375     {
376       ret << "Z Array :\n";
377       _z_array->reprZipWithoutNameStream(ret);
378     }
379   return ret.str();
380 }
381
382 std::string MEDCouplingCMesh::advancedRepr() const
383 {
384   return simpleRepr();
385 }
386
387 /*!
388  * Returns a DataArrayDouble holding positions of nodes along a given axis.
389  * For more info on Cartesian meshes, see \ref MEDCouplingCMeshPage.
390  *  \param [in] i - an index of axis, a value from [0,1,2].
391  *  \return const DataArrayDouble * - a pointer to the data array of node coordinates
392  *         referred by \a this mesh.
393  *  \throw If \a i is not one of [0,1,2].
394  *
395  *  \if ENABLE_EXAMPLES
396  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
397  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
398  *  \endif
399  */
400 const DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) const
401 {
402   switch(i)
403   {
404     case 0:
405       return _x_array;
406     case 1:
407       return _y_array;
408     case 2:
409       return _z_array;
410     default:
411       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
412   }
413 }
414
415 /*!
416  * Returns a DataArrayDouble holding positions of nodes along a given axis.
417  * For more info on Cartesian meshes, see \ref MEDCouplingCMeshPage.
418  *  \param [in] i - an index of axis, a value from [0,1,2].
419  *  \return const DataArrayDouble * - a pointer to the data array of node coordinates
420  *         referred by \a this mesh.
421  *  \throw If \a i is not one of [0,1,2].
422  *
423  *  \if ENABLE_EXAMPLES
424  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
425  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
426  *  \endif
427  */
428 DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i)
429 {
430   switch(i)
431   {
432     case 0:
433       return _x_array;
434     case 1:
435       return _y_array;
436     case 2:
437       return _z_array;
438     default:
439       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
440   }
441 }
442
443 /*!
444  * Sets node coordinates along a given axis. For more info on Cartesian meshes, see 
445  * \ref MEDCouplingCMeshPage.
446  *  \param [in] i - an index of axis, a value in range [0,1,2].
447  *  \param [in] arr - DataArrayDouble holding positions of nodes along the i-th
448  *         axis. It must be an array of one component.
449  *  \throw If \a arr->getNumberOfComponents() != 1.
450  *  \throw If \a i is not one of [0,1,2].
451  *
452  *  \if ENABLE_EXAMPLES
453  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
454  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
455  *  \endif
456  */
457 void MEDCouplingCMesh::setCoordsAt(int i, const DataArrayDouble *arr)
458 {
459   if(arr)
460     arr->checkNbOfComps(1,"MEDCouplingCMesh::setCoordsAt");
461   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
462   if(i<0 || i>2)
463     throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
464   if(arr!=*(thisArr[i]))
465     {
466       if(*(thisArr[i]))
467         (*(thisArr[i]))->decrRef();
468       (*(thisArr[i]))=const_cast<DataArrayDouble *>(arr);
469       if(*(thisArr[i]))
470         (*(thisArr[i]))->incrRef();
471       declareAsNew();
472     }
473 }
474
475 /*!
476  * Sets node coordinates along some of the tree axes. This method updates all the
477  * three node coordinates arrays at once. For more info on Cartesian meshes, see 
478  * \ref MEDCouplingCMeshPage.
479  *  \param [in] coordsX - DataArrayDouble holding positions of nodes along the X
480  *         axis. It must be an array of one component or \c NULL.
481  *  \param [in] coordsY - DataArrayDouble holding positions of nodes along the Y
482  *         axis. It must be an array of one component or \c NULL.
483  *  \param [in] coordsZ - DataArrayDouble holding positions of nodes along the Z
484  *         axis. It must be an array of one component or \c NULL.
485  *  \throw If \a coords*->getNumberOfComponents() != 1.
486  *
487  *  \if ENABLE_EXAMPLES
488  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
489  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
490  *  \endif
491  */
492 void MEDCouplingCMesh::setCoords(const DataArrayDouble *coordsX, const DataArrayDouble *coordsY, const DataArrayDouble *coordsZ)
493 {
494   if(coordsX)
495     coordsX->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsX");
496   if(coordsY)
497     coordsY->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsY");
498   if(coordsZ)
499     coordsZ->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsZ");
500   if(_x_array)
501     _x_array->decrRef();
502   _x_array=const_cast<DataArrayDouble *>(coordsX);
503   if(_x_array)
504     _x_array->incrRef();
505   if(_y_array)
506     _y_array->decrRef();
507   _y_array=const_cast<DataArrayDouble *>(coordsY);
508   if(_y_array)
509     _y_array->incrRef();
510   if(_z_array)
511     _z_array->decrRef();
512   _z_array=const_cast<DataArrayDouble *>(coordsZ);
513   if(_z_array)
514     _z_array->incrRef();
515   declareAsNew();
516 }
517
518 void MEDCouplingCMesh::getBoundingBox(double *bbox) const
519 {
520   int dim=getSpaceDimension();
521   int j=0;
522   for (int idim=0; idim<dim; idim++)
523     {
524       const DataArrayDouble *c=getCoordsAt(idim);
525       if(c)
526         {
527           const double *coords=c->getConstPointer();
528           int nb=c->getNbOfElems();
529           bbox[2*j]=coords[0];
530           bbox[2*j+1]=coords[nb-1];
531           j++;
532         }
533     }
534 }
535
536 /*!
537  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
538  * mesh.<br>
539  * For 1D cells, the returned field contains lengths.<br>
540  * For 2D cells, the returned field contains areas.<br>
541  * For 3D cells, the returned field contains volumes.
542  *  \param [in] isAbs - a not used parameter.
543  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
544  *         and one time . The caller is to delete this field using decrRef() as it is no
545  *         more needed.
546  */
547 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
548 {
549   std::string name="MeasureOfMesh_";
550   name+=getName();
551   int nbelem=getNumberOfCells();
552   MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
553   field->setName(name);
554   DataArrayDouble* array=DataArrayDouble::New();
555   array->alloc(nbelem,1);
556   double *area_vol=array->getPointer();
557   field->setArray(array) ;
558   array->decrRef();
559   field->setMesh(const_cast<MEDCouplingCMesh *>(this));
560   field->synchronizeTimeWithMesh();
561   int tmp[3];
562   getSplitCellValues(tmp);
563   int dim=getSpaceDimension();
564   const double **thisArr=new const double *[dim];
565   const DataArrayDouble *thisArr2[3]={_x_array,_y_array,_z_array};
566   for(int i=0;i<dim;i++)
567     thisArr[i]=thisArr2[i]->getConstPointer();
568   for(int icell=0;icell<nbelem;icell++)
569     {
570       int tmp2[3];
571       GetPosFromId(icell,dim,tmp,tmp2);
572       area_vol[icell]=1.;
573       for(int i=0;i<dim;i++)
574         area_vol[icell]*=thisArr[i][tmp2[i]+1]-thisArr[i][tmp2[i]];
575     }
576   delete [] thisArr;
577   return field;
578 }
579
580 /*!
581  * not implemented yet !
582  */
583 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
584 {
585   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getMeasureFieldOnNode : not implemented yet !");
586   //return 0;
587 }
588
589 int MEDCouplingCMesh::getCellContainingPoint(const double *pos, double eps) const
590 {
591   int dim=getSpaceDimension();
592   int ret=0;
593   int coeff=1;
594   for(int i=0;i<dim;i++)
595     {
596       const double *d=getCoordsAt(i)->getConstPointer();
597       int nbOfNodes=getCoordsAt(i)->getNbOfElems();
598       double ref=pos[i];
599       const double *w=std::find_if(d,d+nbOfNodes,std::bind2nd(std::greater_equal<double>(),ref));
600       int w2=(int)std::distance(d,w);
601       if(w2<nbOfNodes)
602         {
603           if(w2==0)
604             {
605               if(ref>d[0]-eps)
606                 w2=1;
607               else
608                 return -1;
609             }
610           ret+=coeff*(w2-1);
611           coeff*=nbOfNodes-1;
612         }
613       else
614         return -1;
615     }
616   return ret;
617 }
618
619 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
620 {
621   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to untructured mesh to apply it !");
622 }
623
624 /*!
625  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
626  * component of the \a vector to all node coordinates of a corresponding axis.
627  *  \param [in] vector - the translation vector whose size must be not less than \a
628  *         this->getSpaceDimension().
629  */
630 void MEDCouplingCMesh::translate(const double *vector)
631 {
632   if(_x_array)
633     std::transform(_x_array->getConstPointer(),_x_array->getConstPointer()+_x_array->getNbOfElems(),
634                    _x_array->getPointer(),std::bind2nd(std::plus<double>(),vector[0]));
635   if(_y_array)
636     std::transform(_y_array->getConstPointer(),_y_array->getConstPointer()+_y_array->getNbOfElems(),
637                    _y_array->getPointer(),std::bind2nd(std::plus<double>(),vector[1]));
638   if(_z_array)
639     std::transform(_z_array->getConstPointer(),_z_array->getConstPointer()+_z_array->getNbOfElems(),
640                    _z_array->getPointer(),std::bind2nd(std::plus<double>(),vector[2]));
641 }
642
643 /*!
644  * Applies scaling transformation to all nodes of \a this mesh.
645  *  \param [in] point - coordinates of a scaling center. This array is to be of
646  *         size \a this->getSpaceDimension() at least.
647  *  \param [in] factor - a scale factor.
648  */
649 void MEDCouplingCMesh::scale(const double *point, double factor)
650 {
651   for(int i=0;i<3;i++)
652     {
653       DataArrayDouble *c=getCoordsAt(i);
654       if(c)
655         {
656           double *coords=c->getPointer();
657           int lgth=c->getNbOfElems();
658           std::transform(coords,coords+lgth,coords,std::bind2nd(std::minus<double>(),point[i]));
659           std::transform(coords,coords+lgth,coords,std::bind2nd(std::multiplies<double>(),factor));
660           std::transform(coords,coords+lgth,coords,std::bind2nd(std::plus<double>(),point[i]));
661           c->declareAsNew();
662         }
663     }
664   updateTime();
665 }
666
667 MEDCouplingMesh *MEDCouplingCMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
668 {
669   //not implemented yet !
670   return 0;
671 }
672
673 /*!
674  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
675  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
676  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
677  *          components. The caller is to delete this array using decrRef() as it is
678  *          no more needed.
679  */
680 DataArrayDouble *MEDCouplingCMesh::getCoordinatesAndOwner() const
681 {
682   DataArrayDouble *ret=DataArrayDouble::New();
683   int spaceDim=getSpaceDimension();
684   int nbNodes=getNumberOfNodes();
685   ret->alloc(nbNodes,spaceDim);
686   double *pt=ret->getPointer();
687   int tmp[3];
688   getSplitNodeValues(tmp);
689   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
690   const double *tabsPtr[3];
691   for(int j=0;j<spaceDim;j++)
692     {
693       tabsPtr[j]=tabs[j]->getConstPointer();
694       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
695     }
696   int tmp2[3];
697   for(int i=0;i<nbNodes;i++)
698     {
699       GetPosFromId(i,spaceDim,tmp,tmp2);
700       for(int j=0;j<spaceDim;j++)
701         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
702     }
703   return ret;
704 }
705
706 /*!
707  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
708  * computed by averaging coordinates of cell nodes.
709  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
710  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
711  *          components. The caller is to delete this array using decrRef() as it is
712  *          no more needed.
713  */
714 DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const
715 {
716   DataArrayDouble *ret=DataArrayDouble::New();
717   int spaceDim=getSpaceDimension();
718   int nbCells=getNumberOfCells();
719   ret->alloc(nbCells,spaceDim);
720   double *pt=ret->getPointer();
721   int tmp[3];
722   getSplitCellValues(tmp);
723   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
724   std::vector<double> tabsPtr[3];
725   for(int j=0;j<spaceDim;j++)
726     {
727       int sz=tabs[j]->getNbOfElems()-1;
728       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
729       const double *srcPtr=tabs[j]->getConstPointer();
730       tabsPtr[j].insert(tabsPtr[j].end(),srcPtr,srcPtr+sz);
731       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),srcPtr+1,tabsPtr[j].begin(),std::plus<double>());
732       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),tabsPtr[j].begin(),std::bind2nd(std::multiplies<double>(),0.5));
733     }
734   int tmp2[3];
735   for(int i=0;i<nbCells;i++)
736     {
737       GetPosFromId(i,spaceDim,tmp,tmp2);
738       for(int j=0;j<spaceDim;j++)
739         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
740     }
741   return ret;
742 }
743
744 DataArrayDouble *MEDCouplingCMesh::computeIsoBarycenterOfNodesPerCell() const
745 {
746   return MEDCouplingCMesh::getBarycenterAndOwner();
747 }
748
749 void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check)
750 {
751   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !");
752 }
753
754 void MEDCouplingCMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
755 {
756   int it,order;
757   double time=getTime(it,order);
758   tinyInfo.clear();
759   tinyInfoD.clear();
760   littleStrings.clear();
761   littleStrings.push_back(getName());
762   littleStrings.push_back(getDescription());
763   littleStrings.push_back(getTimeUnit());
764   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
765   for(int i=0;i<3;i++)
766     {
767       int val=-1;
768       std::string st;
769       if(thisArr[i])
770         {
771           val=thisArr[i]->getNumberOfTuples();
772           st=thisArr[i]->getInfoOnComponent(0);
773         }
774       tinyInfo.push_back(val);
775       littleStrings.push_back(st);
776     }
777   tinyInfo.push_back(it);
778   tinyInfo.push_back(order);
779   tinyInfoD.push_back(time);
780 }
781
782 void MEDCouplingCMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
783 {
784   a1->alloc(0,1);
785   int sum=0;
786   for(int i=0;i<3;i++)
787     if(tinyInfo[i]!=-1)
788       sum+=tinyInfo[i];
789   a2->alloc(sum,1);
790 }
791
792 void MEDCouplingCMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
793 {
794   a1=DataArrayInt::New();
795   a1->alloc(0,1);
796   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
797   int sz=0;
798   for(int i=0;i<3;i++)
799     {
800       if(thisArr[i])
801         sz+=thisArr[i]->getNumberOfTuples();
802     }
803   a2=DataArrayDouble::New();
804   a2->alloc(sz,1);
805   double *a2Ptr=a2->getPointer();
806   for(int i=0;i<3;i++)
807     if(thisArr[i])
808       a2Ptr=std::copy(thisArr[i]->getConstPointer(),thisArr[i]->getConstPointer()+thisArr[i]->getNumberOfTuples(),a2Ptr);
809 }
810
811 void MEDCouplingCMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
812                                        const std::vector<std::string>& littleStrings)
813 {
814   setName(littleStrings[0]);
815   setDescription(littleStrings[1]);
816   setTimeUnit(littleStrings[2]);
817   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
818   const double *data=a2->getConstPointer();
819   for(int i=0;i<3;i++)
820     {
821       if(tinyInfo[i]!=-1)
822         {
823           (*(thisArr[i]))=DataArrayDouble::New();
824           (*(thisArr[i]))->alloc(tinyInfo[i],1);
825           (*(thisArr[i]))->setInfoOnComponent(0,littleStrings[i+3]);
826           std::copy(data,data+tinyInfo[i],(*(thisArr[i]))->getPointer());
827           data+=tinyInfo[i];
828         }
829     }
830   setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
831 }
832
833 void MEDCouplingCMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
834 {
835   std::ostringstream extent;
836   DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
837   for(int i=0;i<3;i++)
838     {
839       if(thisArr[i])
840         { extent << "0 " <<  thisArr[i]->getNumberOfTuples()-1 << " "; }
841       else
842         { extent << "0 0 "; }
843     }
844   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
845   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
846   ofs << "      <PointData>\n" << pointData << std::endl;
847   ofs << "      </PointData>\n";
848   ofs << "      <CellData>\n" << cellData << std::endl;
849   ofs << "      </CellData>\n";
850   ofs << "      <Coordinates>\n";
851   for(int i=0;i<3;i++)
852     {
853       if(thisArr[i])
854         thisArr[i]->writeVTK(ofs,8,"Array",byteData);
855       else
856         {
857           MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::New(); coo->alloc(1,1);
858           coo->setIJ(0,0,0.);
859           coo->writeVTK(ofs,8,"Array",byteData);
860         }
861     }
862   ofs << "      </Coordinates>\n";
863   ofs << "    </Piece>\n";
864   ofs << "  </" << getVTKDataSetType() << ">\n";
865 }
866
867 void MEDCouplingCMesh::reprQuickOverview(std::ostream& stream) const
868 {
869   stream << "MEDCouplingCMesh C++ instance at " << this << ". Name : \"" << getName() << "\".";
870   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
871   std::ostringstream stream2[3];
872   bool isDef[3];
873   int nbOfCells=1,nbOfNodes=1;
874   for(int i=0;i<3;i++)
875     {
876       isDef[i]=thisArr[i]!=0;
877       if(isDef[i])
878         {    
879           char tmp='X'+i;
880           stream2[i] << tmp << " positions array ";
881           if(!thisArr[i]->isAllocated())
882             stream2[i] << "set but not allocated.";
883           else
884             {
885               int nbCompo=thisArr[i]->getNumberOfComponents();
886               if(nbCompo==1)
887                 {
888                   int nbTuples=thisArr[i]->getNumberOfTuples();
889                   if(nbTuples<1)
890                     { stream2[i] << "set and allocated - WARNING number of elements < 1 !"; nbOfCells=-1; nbOfNodes=-1; }
891                   else
892                     {
893                       stream2[i] << "(length=" << nbTuples << ")" << ": ";
894                       thisArr[i]->reprQuickOverviewData(stream2[i],200);
895                       if(nbOfCells!=-1)
896                         { nbOfNodes*=nbTuples; nbOfCells*=nbTuples-1; }
897                     }
898                 }
899               else
900                 { stream2[i] << "set and allocated - WARNING number of components != 1 !"; nbOfCells=-1; nbOfNodes=-1; }
901             }
902         }
903     }
904   if(!isDef[0] && !isDef[1] && !isDef[2])
905     { stream << " No arrays set !"; return; }
906   if(nbOfCells>=0)
907     { stream << std::endl << "Number of cells : " << nbOfCells << ". Number of nodes : " << nbOfNodes << "."; }
908   for(int i=0;i<3;i++)
909     {
910       if(isDef[i])
911         stream << std::endl << stream2[i].str();
912     }
913 }
914
915 std::string MEDCouplingCMesh::getVTKFileExtension() const
916 {
917   return std::string("vtr");
918 }
919
920 std::string MEDCouplingCMesh::getVTKDataSetType() const
921 {
922   return std::string("RectilinearGrid");
923 }