Salome HOME
Minor doc: isIota()
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingPointSet.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 (CEA/DEN)
20
21 #include "MEDCouplingPointSet.hxx"
22 #include "MCAuto.hxx"
23 #include "MEDCoupling1GTUMesh.hxx"
24 #include "MEDCouplingUMesh.hxx"
25 #include "MEDCouplingMemArray.hxx"
26 #include "PlanarIntersector.txx"
27 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
28 #include "InterpKernelGeo2DNode.hxx"
29 #include "DirectedBoundingBox.hxx"
30 #include "InterpKernelAutoPtr.hxx"
31
32 #include <cmath>
33 #include <limits>
34 #include <numeric>
35
36 using namespace MEDCoupling;
37
38 MEDCouplingPointSet::MEDCouplingPointSet():_coords(0)
39 {
40 }
41
42 MEDCouplingPointSet::MEDCouplingPointSet(const MEDCouplingPointSet& other, bool deepCopy):MEDCouplingMesh(other),_coords(0)
43 {
44   if(other._coords)
45     _coords=other._coords->performCopyOrIncrRef(deepCopy);
46 }
47
48 MEDCouplingPointSet::~MEDCouplingPointSet()
49 {
50   if(_coords)
51     _coords->decrRef();
52 }
53
54 int MEDCouplingPointSet::getNumberOfNodes() const
55 {
56   if(_coords)
57     return _coords->getNumberOfTuples();
58   else
59     throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !");
60 }
61
62 int MEDCouplingPointSet::getSpaceDimension() const
63 {
64   if(_coords)
65     return _coords->getNumberOfComponents();
66   else
67     throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !");
68 }
69
70 void MEDCouplingPointSet::updateTime() const
71 {
72   if(_coords)
73     {
74       updateTimeWith(*_coords);
75     }
76 }
77
78 std::size_t MEDCouplingPointSet::getHeapMemorySizeWithoutChildren() const
79 {
80   return MEDCouplingMesh::getHeapMemorySizeWithoutChildren();
81 }
82
83 std::vector<const BigMemoryObject *> MEDCouplingPointSet::getDirectChildrenWithNull() const
84 {
85   std::vector<const BigMemoryObject *> ret;
86   ret.push_back(_coords);
87   return ret;
88 }
89
90 void MEDCouplingPointSet::setCoords(const DataArrayDouble *coords)
91 {
92   if( coords != _coords )
93     {
94       if (_coords)
95         _coords->decrRef();
96       _coords=const_cast<DataArrayDouble *>(coords);
97       if(_coords)
98         _coords->incrRef();
99       declareAsNew();
100     }
101 }
102
103 /*!
104  * Returns a pointer to the array of point coordinates held by \a this.
105  *  \return DataArrayDouble * - the pointer to the array of point coordinates. The
106  *          caller is to delete this array using decrRef() as it is no more needed.
107  */
108 DataArrayDouble *MEDCouplingPointSet::getCoordinatesAndOwner() const
109 {
110   if(_coords)
111     _coords->incrRef();
112   return _coords;
113 }
114
115 /*!
116  * Copies string attributes from an \a other mesh. The copied strings are
117  * - mesh name
118  * - mesh description
119  * - time units
120  * - textual data of the coordinates array (name and components info)
121  *
122  *  \param [in] other - the mesh to copy string attributes from.
123  */
124 void MEDCouplingPointSet::copyTinyStringsFrom(const MEDCouplingMesh *other)
125 {
126   MEDCouplingMesh::copyTinyStringsFrom(other);
127   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
128   if(!otherC)
129     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::copyTinyStringsFrom : meshes have not same type !");
130   if(_coords && otherC->_coords)
131     _coords->copyStringInfoFrom(*otherC->_coords);
132 }
133
134 bool MEDCouplingPointSet::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
135 {
136   if(!other)
137     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::isEqualIfNotWhy : null mesh instance in input !");
138   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
139   if(!otherC)
140     {
141       reason="mesh given in input is not castable in MEDCouplingPointSet !";
142       return false;
143     }
144   if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
145     return false;
146   if(!areCoordsEqualIfNotWhy(*otherC,prec,reason))
147     return false;
148   return true;
149 }
150
151 /*!
152  * Checks equality of point coordinates with coordinates of an \a other mesh.
153  *        None textual data is considered.
154  *  \param [in] other - the mesh to compare coordinates with \a this one.
155  *  \param [in] prec - precision value to compare coordinates.
156  *  \return bool - \a true if coordinates of points are equal, \a false else.
157  */
158 bool MEDCouplingPointSet::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
159 {
160   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
161   if(!otherC)
162     return false;
163   if(!areCoordsEqualWithoutConsideringStr(*otherC,prec))
164     return false;
165   return true;
166 }
167
168 bool MEDCouplingPointSet::areCoordsEqualIfNotWhy(const MEDCouplingPointSet& other, double prec, std::string& reason) const
169 {
170   if(_coords==0 && other._coords==0)
171     return true;
172   if(_coords==0 || other._coords==0)
173     {
174       reason="Only one PointSet between the two this and other has coordinate defined !";
175       return false;
176     }
177   if(_coords==other._coords)
178     return true;
179   bool ret=_coords->isEqualIfNotWhy(*other._coords,prec,reason);
180   if(!ret)
181     reason.insert(0,"Coordinates DataArray do not match : ");
182   return ret;
183 }
184
185 /*!
186  * Checks equality of point coordinates with \a other point coordinates.
187  *        Textual data (name and components info) \b is compared as well.
188  *  \param [in] other - the point coordinates to compare with \a this one.
189  *  \param [in] prec - precision value to compare coordinates.
190  *  \return bool - \a true if coordinates of points are equal, \a false else.
191  */
192 bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const
193 {
194   std::string tmp;
195   return areCoordsEqualIfNotWhy(other,prec,tmp);
196 }
197
198 /*!
199  * Checks equality of point coordinates with \a other point coordinates.
200  *        None textual data is considered.
201  *  \param [in] other - the point coordinates to compare with \a this one.
202  *  \param [in] prec - precision value to compare coordinates.
203  *  \return bool - \a true if coordinates of points are equal, \a false else.
204  */
205 bool MEDCouplingPointSet::areCoordsEqualWithoutConsideringStr(const MEDCouplingPointSet& other, double prec) const
206 {
207   if(_coords==0 && other._coords==0)
208     return true;
209   if(_coords==0 || other._coords==0)
210     return false;
211   if(_coords==other._coords)
212     return true;
213   return _coords->isEqualWithoutConsideringStr(*other._coords,prec);
214 }
215
216 /*!
217  * Returns coordinates of \a nodeId-th node.
218  *  \param [in] nodeId - the ID of the node of interest.
219  *  \param [in, out] coo - the array filled with coordinates of the \a nodeId-th
220  *         node. This array is not cleared before filling in, the coordinates are
221  *         appended to its end.
222  *  \throw If the coordinates array is not set.
223  *  \throw If \a nodeId is not a valid index for the coordinates array.
224  *
225  *  \if ENABLE_EXAMPLES
226  *  \ref cpp_mcpointset_getcoordinatesofnode "Here is a C++ example".<br>
227  *  \ref  py_mcpointset_getcoordinatesofnode "Here is a Python example".
228  *  \endif
229  */
230 void MEDCouplingPointSet::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
231 {
232   if(!_coords)
233     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCoordinatesOfNode : no coordinates array set !");
234   int nbNodes=getNumberOfNodes();
235   if(nodeId>=0 && nodeId<nbNodes)
236     {
237       const double *cooPtr=_coords->getConstPointer();
238       int spaceDim=getSpaceDimension();
239       coo.insert(coo.end(),cooPtr+spaceDim*nodeId,cooPtr+spaceDim*(nodeId+1));
240     }
241   else
242     {
243       std::ostringstream oss; oss << "MEDCouplingPointSet::getCoordinatesOfNode : request of nodeId \"" << nodeId << "\" but it should be in [0,"<< nbNodes << ") !";
244       throw INTERP_KERNEL::Exception(oss.str().c_str());
245     }
246 }
247
248 /*!
249  * Finds nodes equal within \a precision and returns an array describing the 
250  * permutation to remove duplicated nodes.
251  *  \param [in] precision - minimal absolute distance between two nodes at which they are
252  *              considered not coincident.
253  *  \param [in] limitNodeId - limit node id. If all nodes within a group of coincident
254  *              nodes have id strictly lower than \a limitTupleId then they are not
255  *              returned. Put -1 to this parameter to have all nodes returned.
256  *  \param [out] areNodesMerged - is set to \a true if any coincident nodes found.
257  *  \param [out] newNbOfNodes - returns number of unique nodes.
258  *  \return DataArrayInt * - the permutation array in "Old to New" mode. For more 
259  *          info on "Old to New" mode see \ref numbering. The caller
260  *          is to delete this array using decrRef() as it is no more needed.
261  *  \throw If the coordinates array is not set.
262  */
263 DataArrayInt *MEDCouplingPointSet::buildPermArrayForMergeNode(double precision, int limitNodeId, bool& areNodesMerged, int& newNbOfNodes) const
264 {
265   DataArrayInt *comm,*commI;
266   findCommonNodes(precision,limitNodeId,comm,commI);
267   int oldNbOfNodes=getNumberOfNodes();
268   MCAuto<DataArrayInt> ret=buildNewNumberingFromCommonNodesFormat(comm,commI,newNbOfNodes);
269   areNodesMerged=(oldNbOfNodes!=newNbOfNodes);
270   comm->decrRef();
271   commI->decrRef();
272   return ret.retn();
273 }
274
275 /*!
276  * Finds nodes coincident within \a prec tolerance.
277  * Ids of coincident nodes are stored in output arrays in the \ref numbering-indirect format.
278  *  \param [in] prec - minimal absolute distance (using infinite norm) between two nodes at which they are
279  *              considered not coincident.
280  *  \param [in] limitNodeId - limit node id. If all nodes within a group of coincident
281  *              nodes have id strictly lower than \a limitTupleId then they are not
282  *              returned. Put -1 to this parameter to have all nodes treated.
283  *  \param [out] comm - the array holding ids of coincident nodes.
284  *               \a comm->getNumberOfComponents() == 1. 
285  *               \a comm->getNumberOfTuples() == \a commIndex->back(). The caller
286  *               is to delete this array using decrRef() as it is no more needed.
287  *  \param [out] commIndex - the array dividing all ids stored in \a comm into
288  *               groups of (ids of) coincident nodes (\ref numbering-indirect). Its every value is a tuple
289  *               index where a next group of nodes begins. For example the second
290  *               group of nodes in \a comm is described by following range of indices:
291  *               [ \a commIndex[1], \a commIndex[2] ). \a commIndex->getNumberOfTuples()-1
292  *               gives the number of groups of coincident nodes. The caller
293  *               is to delete this array using decrRef() as it is no more needed.
294  *  \throw If the coordinates array is not set.
295  *
296  *  \if ENABLE_EXAMPLES
297  *  \ref cpp_mcpointset_findcommonnodes "Here is a C++ example".<br>
298  *  \ref  py_mcpointset_findcommonnodes "Here is a Python example".
299  *  \endif
300  */
301 void MEDCouplingPointSet::findCommonNodes(double prec, int limitNodeId, DataArrayInt *&comm, DataArrayInt *&commIndex) const
302 {
303   if(!_coords)
304     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findCommonNodes : no coords specified !");
305   _coords->findCommonTuples(prec,limitNodeId,comm,commIndex);
306 }
307
308 /*!
309  * Finds nodes located at distances lower that \a eps from a given point.
310  *  \param [in] pos - pointer to coordinates of the point.  This array is expected to
311  *         be of length \a this->getSpaceDimension() at least, else the
312  *         behavior is not warranted.
313  *  \param [in] eps - the lowest distance between a point and a node (using infinite norm) at which the node is
314  *          not returned by this method.
315  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of nodes
316  *          close to the point. The caller is to delete this
317  *          array using decrRef() as it is no more needed.
318  *  \throw If the coordinates array is not set.
319  *
320  *  \if ENABLE_EXAMPLES
321  *  \ref cpp_mcpointset_getnodeidsnearpoint "Here is a C++ example".<br>
322  *  \ref  py_mcpointset_getnodeidsnearpoint "Here is a Python example".
323  *  \endif
324  */
325 DataArrayInt *MEDCouplingPointSet::getNodeIdsNearPoint(const double *pos, double eps) const
326 {
327   DataArrayInt *c=0,*cI=0;
328   getNodeIdsNearPoints(pos,1,eps,c,cI);
329   MCAuto<DataArrayInt> cITmp(cI);
330   return c;
331 }
332
333 /*!
334  * Finds nodes located at distances lower that \a eps from given points.
335  *  \param [in] pos - pointer to coordinates of the points. This array is expected to
336  *         be of length \a nbOfPoints * \a this->getSpaceDimension() at least, else the
337  *         behavior is not warranted.
338  *  \param [in] nbOfPoints - number of points whose coordinates are given by \a pos
339  *         parameter. 
340  *  \param [in] eps - the lowest distance between (using infinite norm) a point and a node at which the node is
341  *         not returned by this method.
342  *  \param [out] c - array (\ref numbering-indirect) returning ids of nodes located closer than \a eps to the
343  *         given points. The caller
344  *         is to delete this array using decrRef() as it is no more needed.
345  *  \param [out] cI - for each i-th given point, the array specifies tuples of \a c
346  *         holding ids of nodes close to the i-th point (\ref numbering-indirect). <br>The i-th value of \a cI is an
347  *         index of tuple of \a c holding id of a first (if any) node close to the
348  *         i-th given point. Difference between the i-th and (i+1)-th value of \a cI
349  *         (i.e. \a cI[ i+1 ] - \a cI[ i ]) defines number of nodes close to the i-th
350  *         point (that can be zero!). For example, the group of nodes close to the
351  *         second point is described by following range of indices [ \a cI[1], \a cI[2] ).
352  *         The caller is to delete this array using decrRef() as it is no more needed.
353  *  \throw If the coordinates array is not set.
354  *
355  *  \if ENABLE_EXAMPLES
356  *  \ref cpp_mcpointset_getnodeidsnearpoints "Here is a C++ example".<br>
357  *  \ref  py_mcpointset_getnodeidsnearpoints "Here is a Python example".
358  *  \endif
359  */
360 void MEDCouplingPointSet::getNodeIdsNearPoints(const double *pos, int nbOfPoints, double eps, DataArrayInt *& c, DataArrayInt *& cI) const
361 {
362   if(!_coords)
363     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : no coordiantes set !");
364   int spaceDim=getSpaceDimension();
365   MCAuto<DataArrayDouble> points=DataArrayDouble::New();
366   points->useArray(pos,false,CPP_DEALLOC,nbOfPoints,spaceDim);
367   _coords->computeTupleIdsNearTuples(points,eps,c,cI);
368 }
369
370 /*!
371  * @param comm in param in the same format than one returned by findCommonNodes method (\ref numbering-indirect).
372  * @param commI in param in the same format than one returned by findCommonNodes method (\ref numbering-indirect).
373  * @return the old to new correspondance array.
374  */
375 DataArrayInt *MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat(const DataArrayInt *comm, const DataArrayInt *commIndex,
376                                                                           int& newNbOfNodes) const
377 {
378   if(!_coords)
379     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat : no coords specified !");
380   return DataArrayInt::ConvertIndexArrayToO2N(getNumberOfNodes(),comm->begin(),commIndex->begin(),commIndex->end(),newNbOfNodes);
381 }
382
383 /*!
384  * Permutes and possibly removes nodes as specified by \a newNodeNumbers array.
385  * If \a newNodeNumbers[ i ] < 0 then the i-th node is removed, 
386  * else \a newNodeNumbers[ i ] is a new id of the i-th node. The nodal connectivity
387  * array is modified accordingly.
388  *  \param [in] newNodeNumbers - a permutation array, of length \a
389  *         this->getNumberOfNodes(), in "Old to New" mode. 
390  *         See \ref numbering for more info on renumbering modes.
391  *  \param [in] newNbOfNodes - number of nodes remaining after renumbering.
392  *  \throw If the coordinates array is not set.
393  *  \throw If the nodal connectivity of cells is not defined.
394  *
395  *  \if ENABLE_EXAMPLES
396  *  \ref cpp_mcumesh_renumberNodes "Here is a C++ example".<br>
397  *  \ref  py_mcumesh_renumberNodes "Here is a Python example".
398  *  \endif
399  */
400 void MEDCouplingPointSet::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
401 {
402   if(!_coords)
403     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::renumberNodes : no coords specified !");
404   MCAuto<DataArrayDouble> newCoords=_coords->renumberAndReduce(newNodeNumbers,newNbOfNodes);
405   renumberNodesInConn(newNodeNumbers);
406   setCoords(newCoords);//let it here not before renumberNodesInConn because old number of nodes is sometimes used...
407 }
408
409 /*!
410  * Permutes and possibly removes nodes as specified by \a newNodeNumbers array.
411  * If \a newNodeNumbers[ i ] < 0 then the i-th node is removed, 
412  * else \a newNodeNumbers[ i ] is a new id of the i-th node. The nodal connectivity
413  * array is modified accordingly. In contrast to renumberNodes(), location
414  * of merged nodes (whose new ids coincide) is changed to be at their barycenter.
415  *  \param [in] newNodeNumbers - a permutation array, of length \a
416  *         this->getNumberOfNodes(), in "Old to New" mode. 
417  *         See \ref numbering for more info on renumbering modes.
418  *  \param [in] newNbOfNodes - number of nodes remaining after renumbering, which is
419  *         actually one more than the maximal id in \a newNodeNumbers.
420  *  \throw If the coordinates array is not set.
421  *  \throw If the nodal connectivity of cells is not defined.
422  *
423  *  \if ENABLE_EXAMPLES
424  *  \ref cpp_mcumesh_renumberNodes "Here is a C++ example".<br>
425  *  \ref  py_mcumesh_renumberNodes "Here is a Python example".
426  *  \endif
427  */
428 void MEDCouplingPointSet::renumberNodesCenter(const int *newNodeNumbers, int newNbOfNodes)
429 {
430   DataArrayDouble *newCoords=DataArrayDouble::New();
431   std::vector<int> div(newNbOfNodes);
432   int spaceDim=getSpaceDimension();
433   newCoords->alloc(newNbOfNodes,spaceDim);
434   newCoords->copyStringInfoFrom(*_coords);
435   newCoords->fillWithZero();
436   int oldNbOfNodes=getNumberOfNodes();
437   double *ptToFill=newCoords->getPointer();
438   const double *oldCoordsPtr=_coords->getConstPointer();
439   for(int i=0;i<oldNbOfNodes;i++)
440     {
441       std::transform(oldCoordsPtr+i*spaceDim,oldCoordsPtr+(i+1)*spaceDim,ptToFill+newNodeNumbers[i]*spaceDim,
442                      ptToFill+newNodeNumbers[i]*spaceDim,std::plus<double>());
443       div[newNodeNumbers[i]]++;
444     }
445   for(int i=0;i<newNbOfNodes;i++)
446     ptToFill=std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)div[i]));
447   setCoords(newCoords);
448   newCoords->decrRef();
449   renumberNodesInConn(newNodeNumbers);
450 }
451
452 /*!
453  * Computes the minimum box bounding all nodes. The edges of the box are parallel to
454  * the Cartesian coordinate axes. The bounding box is described by coordinates of its
455  * two extremum points with minimal and maximal coordinates.
456  *  \param [out] bbox - array filled with coordinates of extremum points in "no
457  *         interlace" mode, i.e. xMin, xMax, yMin, yMax, zMin, zMax (if in 3D). This
458  *         array, of length 2 * \a this->getSpaceDimension() at least, is to be
459  *         pre-allocated by the caller.
460  *  \throw If the coordinates array is not set.
461  *
462  *  \if ENABLE_EXAMPLES
463  *  \ref cpp_mcpointset_getBoundingBox "Here is a C++ example".<br>
464  *  \ref  py_mcpointset_getBoundingBox "Here is a Python example".
465  *  \endif
466  */
467 void MEDCouplingPointSet::getBoundingBox(double *bbox) const
468 {
469   if(!_coords)
470     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getBoundingBox : Coordinates not set !");
471   _coords->getMinMaxPerComponent(bbox);
472 }
473
474 /*!
475  * Removes "free" nodes, i.e. nodes not used to define any element.
476  *  \throw If the coordinates array is not set.
477  *  \throw If the elements are not defined.
478  */
479 void MEDCouplingPointSet::zipCoords()
480 {
481   checkFullyDefined();
482   DataArrayInt *traducer=zipCoordsTraducer();
483   traducer->decrRef();
484 }
485
486 /*! \cond HIDDEN_ITEMS */
487 struct MEDCouplingCompAbs
488 {
489   bool operator()(double x, double y) { return std::abs(x)<std::abs(y);}
490 };
491 /*! \endcond */
492
493 /*!
494  * Returns the carateristic dimension of \a this point set, that is a maximal
495  * absolute values of node coordinates.
496  *  \throw If the coordinates array is not set.
497  */
498 double MEDCouplingPointSet::getCaracteristicDimension() const
499 {
500   if(!_coords)
501     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCaracteristicDimension : Coordinates not set !");
502   const double *coords=_coords->getConstPointer();
503   int nbOfValues=_coords->getNbOfElems();
504   return std::abs(*std::max_element(coords,coords+nbOfValues,MEDCouplingCompAbs()));
505 }
506
507 /*!
508  * This method recenter coordinates of nodes in \b this in order to be centered at the origin to benefit about the advantages of the precision to be around the box
509  * around origin of 'radius' 1.
510  *
511  * \warning this method is non const and alterates coordinates in \b this without modifying.
512  * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
513  *
514  */
515 void MEDCouplingPointSet::recenterForMaxPrecision(double eps)
516 {
517   if(!_coords)
518     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::recenterForMaxPrecision : Coordinates not set !");
519   _coords->recenterForMaxPrecision(eps);
520   updateTime();
521 }
522
523 /*!
524  * Rotates \a this set of nodes by \a angle around either an axis (in 3D) or a point
525  * (in 2D). 
526  *  \param [in] center - coordinates either of an origin of rotation axis (in 3D) or
527  *         of center of rotation (in 2D). This array is to be of size \a
528  *         this->getSpaceDimension() at least.
529  *  \param [in] vector - 3 components of a vector defining direction of the rotation
530  *         axis in 3D. In 2D this parameter is not used.
531  *  \param [in] angle - the rotation angle in radians.
532  *  \throw If the coordinates array is not set.
533  *  \throw If \a this->getSpaceDimension() != 2 && \a this->getSpaceDimension() != 3.
534  *  \throw If \a center == NULL
535  *  \throw If \a vector == NULL && \a this->getSpaceDimension() == 3.
536  *  \throw If Magnitude of \a vector is zero.
537  *
538  *  \if ENABLE_EXAMPLES
539  *  \ref cpp_mcpointset_rotate "Here is a C++ example".<br>
540  *  \ref  py_mcpointset_rotate "Here is a Python example".
541  *  \endif
542  */
543 void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
544 {
545   int spaceDim=getSpaceDimension();
546   if(spaceDim==3)
547     rotate3D(center,vector,angle);
548   else if(spaceDim==2)
549     rotate2D(center,angle);
550   else
551     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
552   _coords->declareAsNew();
553   updateTime();
554 }
555
556 /*!
557  * Translates \a this set of nodes. 
558  *  \param [in] vector - components of a translation vector. This array is to be of
559  *         size \a this->getSpaceDimension() at least. 
560  *  \throw If the coordinates array is not set.
561  *  \throw If \a vector == NULL.
562  *
563  *  \if ENABLE_EXAMPLES
564  *  \ref cpp_mcpointset_translate "Here is a C++ example".<br>
565  *  \ref  py_mcpointset_translate "Here is a Python example".
566  *  \endif
567  */
568 void MEDCouplingPointSet::translate(const double *vector)
569 {
570   if(!vector)
571     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::translate : NULL input vector !");
572   if(!_coords)
573     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::translate : no coordinates set !");
574   double *coords=_coords->getPointer();
575   int nbNodes=getNumberOfNodes();
576   int dim=getSpaceDimension();
577   for(int i=0; i<nbNodes; i++)
578     for(int idim=0; idim<dim;idim++)
579       coords[i*dim+idim]+=vector[idim];
580   _coords->declareAsNew();
581   updateTime();
582 }
583
584
585 /*!
586  * Applies scaling transformation to \a this set of nodes. 
587  *  \param [in] point - coordinates of a scaling center. This array is to be of
588  *         size \a this->getSpaceDimension() at least.
589  *  \param [in] factor - a scale factor.
590  *  \throw If the coordinates array is not set.
591  *  \throw If \a point == NULL.
592  *
593  *  \if ENABLE_EXAMPLES
594  *  \ref cpp_mcpointset_scale "Here is a C++ example".<br>
595  *  \ref  py_mcpointset_scale "Here is a Python example".
596  *  \endif
597  */
598 void MEDCouplingPointSet::scale(const double *point, double factor)
599 {
600   if(!point)
601     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::scale : NULL input point !");
602   if(!_coords)
603     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::scale : no coordinates set !");
604   double *coords=_coords->getPointer();
605   int nbNodes=getNumberOfNodes();
606   int dim=getSpaceDimension();
607   for(int i=0;i<nbNodes;i++)
608     {
609       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
610       std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
611       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
612     }
613   _coords->declareAsNew();
614   updateTime();
615 }
616
617 /*!
618  * Converts \a this set of points to an other dimension by changing number of
619  * components of point coordinates. If the dimension increases, added components
620  * are filled with \a dftValue. If the dimension decreases, last components are lost.
621  * If the new dimension is same as \a this->getSpaceDimension(), nothing is done.
622  *  \param [in] newSpaceDim - the new space dimension.
623  *  \param [in] dftValue - the value to assign to added components of point coordinates
624  *         (if the dimension increases).
625  *  \throw If the coordinates array is not set.
626  *  \throw If \a newSpaceDim < 1.
627  */
628 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue)
629 {
630   if(getCoords()==0)
631     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
632   if(newSpaceDim<1)
633     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
634   int oldSpaceDim=getSpaceDimension();
635   if(newSpaceDim==oldSpaceDim)
636     return ;
637   DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
638   setCoords(newCoords);
639   newCoords->decrRef();
640   updateTime();
641 }
642
643 /*!
644  * Substitutes \a this->_coords with \a other._coords provided that coordinates of
645  * the two point sets match with a specified precision, else an exception is thrown.
646  *  \param [in] other - the other point set whose coordinates array will be used by
647  *         \a this point set in case of their equality.
648  *  \param [in] epsilon - the precision used to compare coordinates.
649  *  \throw If the coordinates array of \a this is not set.
650  *  \throw If the coordinates array of \a other is not set.
651  *  \throw If the coordinates of \a this and \a other do not match.
652  */
653 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon)
654 {
655   if(_coords==other._coords)
656     return ;
657   if(!_coords)
658     throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
659   if(!other._coords)
660     throw INTERP_KERNEL::Exception("Other instance has no coords whereas current has !");
661   if(!_coords->isEqualWithoutConsideringStr(*other._coords,epsilon))
662     throw INTERP_KERNEL::Exception("Coords are not the same !");
663   setCoords(other._coords);
664 }
665
666 /*!
667  * This method duplicates the nodes whose ids are in [\b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd) and put the result of their duplication at the end
668  * of existing node ids.
669  * 
670  * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
671  * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
672  */
673 void MEDCouplingPointSet::duplicateNodesInCoords(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd)
674 {
675   if(!_coords)
676     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::duplicateNodesInCoords : no coords set !");
677   MCAuto<DataArrayDouble> newCoords=_coords->selectByTupleIdSafe(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd);
678   MCAuto<DataArrayDouble> newCoords2=DataArrayDouble::Aggregate(_coords,newCoords);
679   setCoords(newCoords2);
680 }
681
682 /*!
683  * Finds nodes located at distance lower that \a eps from a specified plane.
684  *  \param [in] pt - 3 components of a point defining location of the plane.
685  *  \param [in] vec - 3 components of a normal vector to the plane. Vector magnitude
686  *         must be greater than 10*\a eps.
687  *  \param [in] eps - maximal distance of a node from the plane at which the node is
688  *         considered to lie on the plane.
689  *  \param [in,out] nodes - a vector returning ids of found nodes. This vector is not
690  *         cleared before filling in.
691  *  \throw If the coordinates array is not set.
692  *  \throw If \a pt == NULL.
693  *  \throw If \a vec == NULL.
694  *  \throw If the magnitude of \a vec is zero.
695  *  \throw If \a this->getSpaceDimension() != 3.
696  */
697 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const
698 {
699   if(getSpaceDimension()!=3)
700     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
701   if(!pt)
702     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : NULL point pointer specified !");
703   if(!vec)
704     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : NULL vector pointer specified !");
705   int nbOfNodes=getNumberOfNodes();
706   double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
707   double deno=sqrt(a*a+b*b+c*c);
708   if(deno<std::numeric_limits<double>::min())
709     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : vector pointer specified has norm equal to 0. !");
710   const double *work=_coords->getConstPointer();
711   for(int i=0;i<nbOfNodes;i++)
712     {
713       if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
714         nodes.push_back(i);
715       work+=3;
716     }
717 }
718
719 /*!
720  * Finds nodes located at distance lower that \a eps from a specified line in 2D and 3D.
721  *  \param [in] pt - components of coordinates of an initial point of the line. This
722  *         array is to be of size \a this->getSpaceDimension() at least.
723  *  \param [in] vec - components of a vector defining the line direction. This array
724  *         is to be of size \a this->getSpaceDimension() at least. Vector magnitude 
725  *         must be greater than 10*\a eps.
726  *  \param [in] eps - maximal distance of a node from the line at which the node is
727  *         considered to lie on the line.
728  *  \param [in,out] nodes - a vector returning ids of found nodes. This vector is not
729  *         cleared before filling in.
730  *  \throw If the coordinates array is not set.
731  *  \throw If \a pt == NULL.
732  *  \throw If \a vec == NULL.
733  *  \throw If the magnitude of \a vec is zero.
734  *  \throw If ( \a this->getSpaceDimension() != 3 && \a this->getSpaceDimension() != 2 ).
735  */
736 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const
737 {
738   int spaceDim=getSpaceDimension();
739   if(spaceDim!=2 && spaceDim!=3)
740     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
741   if(!pt)
742     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : NULL point pointer specified !");
743   if(!vec)
744     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : NULL vector pointer specified !");
745   int nbOfNodes=getNumberOfNodes();
746   double den=0.;
747   for(int i=0;i<spaceDim;i++)
748     den+=vec[i]*vec[i];
749   double deno=sqrt(den);
750   if(deno<10.*eps)
751     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
752   INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
753   for(int i=0;i<spaceDim;i++)
754     vecn[i]=vec[i]/deno;
755   const double *work=_coords->getConstPointer();
756   if(spaceDim==2)
757     {
758       for(int i=0;i<nbOfNodes;i++)
759         {
760           if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
761             nodes.push_back(i);
762           work+=2;
763         }
764     }
765   else
766     {
767       for(int i=0;i<nbOfNodes;i++)
768         {
769           double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
770           double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
771           double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
772           if(std::sqrt(a*a+b*b+c*c)<eps)
773             nodes.push_back(i);
774           work+=3;
775         }
776     }
777 }
778
779 /*!
780  * Returns a new array of node coordinates by concatenating node coordinates of two
781  * given point sets, so that (1) the number of nodes in the result array is a sum of the
782  * number of nodes of given point sets and (2) the number of component in the result array
783  * is same as that of each of given point sets. Info on components is copied from the first
784  * of the given point set. Space dimension of the given point sets must be the same. 
785  *  \param [in] m1 - a point set whose coordinates will be included in the result array.
786  *  \param [in] m2 - another point set whose coordinates will be included in the
787  *         result array. 
788  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
789  *          The caller is to delete this result array using decrRef() as it is no more
790  *          needed.
791  *  \throw If both \a m1 and \a m2 are NULL.
792  *  \throw If \a m1->getSpaceDimension() != \a m2->getSpaceDimension().
793  */
794 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2)
795 {
796   int spaceDim=m1->getSpaceDimension();
797   if(spaceDim!=m2->getSpaceDimension())
798     throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
799   return DataArrayDouble::Aggregate(m1->getCoords(),m2->getCoords());
800 }
801
802 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms)
803 {
804   if(ms.empty())
805     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : input array must be NON EMPTY !");
806   std::vector<const MEDCouplingPointSet *>::const_iterator it=ms.begin();
807   std::vector<const DataArrayDouble *> coo(ms.size());
808   int spaceDim=(*it)->getSpaceDimension();
809   coo[0]=(*it++)->getCoords();
810   if(!coo[0]->isAllocated())
811     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : first element in coordinates is not allocated !");
812   for(int i=1;it!=ms.end();it++,i++)
813     {
814       const DataArrayDouble *tmp=(*it)->getCoords();
815       if(tmp)
816         {
817           if(tmp->isAllocated())
818             {
819               if((*it)->getSpaceDimension()==spaceDim)
820                 coo[i]=tmp;
821               else
822                 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : Mismatch in SpaceDim !");
823             }
824           else
825             throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : Presence of a non allocated array !");
826         }
827       else
828         throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : Empty coords detected !");
829     }
830   return DataArrayDouble::Aggregate(coo);
831 }
832
833 /*!
834  * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
835  * This method is used during unserialization process.
836  */
837 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
838 {
839   switch(type)
840     {
841     case UNSTRUCTURED:
842       return MEDCouplingUMesh::New();
843     case SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED:
844       return MEDCoupling1SGTUMesh::New();
845     case SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED:
846       return MEDCoupling1DGTUMesh::New();
847     default:
848       throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
849     }
850 }
851
852 /*!
853  * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
854  */
855 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
856 {
857   int it,order;
858   double time=getTime(it,order);
859   if(_coords)
860     {
861       int spaceDim=getSpaceDimension();
862       littleStrings.resize(spaceDim+4);
863       littleStrings[0]=getName();
864       littleStrings[1]=getDescription();
865       littleStrings[2]=_coords->getName();
866       littleStrings[3]=getTimeUnit();
867       for(int i=0;i<spaceDim;i++)
868         littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
869       tinyInfo.clear();
870       tinyInfo.push_back(getType());
871       tinyInfo.push_back(spaceDim);
872       tinyInfo.push_back(getNumberOfNodes());
873       tinyInfo.push_back(it);
874       tinyInfo.push_back(order);
875       tinyInfoD.push_back(time);
876     }
877   else
878     {
879       littleStrings.resize(3);
880       littleStrings[0]=getName();
881       littleStrings[1]=getDescription();
882       littleStrings[2]=getTimeUnit();
883       tinyInfo.clear();
884       tinyInfo.push_back(getType());
885       tinyInfo.push_back(-1);
886       tinyInfo.push_back(-1);
887       tinyInfo.push_back(it);
888       tinyInfo.push_back(order);
889       tinyInfoD.push_back(time);
890     }
891 }
892
893 /*!
894  * Third and final step of serialization process.
895  */
896 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
897 {
898   if(_coords)
899     {
900       a2=const_cast<DataArrayDouble *>(getCoords());
901       a2->incrRef();
902     }
903   else
904     a2=0;
905 }
906
907 /*!
908  * Second step of serialization process.
909  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
910  */
911 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
912 {
913   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
914     {
915       a2->alloc(tinyInfo[2],tinyInfo[1]);
916       littleStrings.resize(tinyInfo[1]+4);
917     }
918   else
919     {
920       littleStrings.resize(3);
921     }
922 }
923
924 /*!
925  * Second and final unserialization process.
926  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
927  */
928 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
929 {
930   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
931     {
932       setCoords(a2);
933       setName(littleStrings[0]);
934       setDescription(littleStrings[1]);
935       a2->setName(littleStrings[2]);
936       setTimeUnit(littleStrings[3]);
937       for(int i=0;i<tinyInfo[1];i++)
938         getCoords()->setInfoOnComponent(i,littleStrings[i+4]);
939       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
940     }
941   else
942     {
943       setName(littleStrings[0]);
944       setDescription(littleStrings[1]);
945       setTimeUnit(littleStrings[2]);
946       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
947     }
948 }
949
950 void MEDCouplingPointSet::checkConsistencyLight() const
951 {
952   if(!_coords)
953     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkConsistencyLight : no coordinates set !");
954 }
955
956 /*!
957  * Intersect Bounding Box given 2 Bounding Boxes.
958  */
959 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
960 {
961   double* bbtemp = new double[2*dim];
962   double deltamax=0.0;
963
964   for (int i=0; i< dim; i++)
965     {
966       double delta = bb1[2*i+1]-bb1[2*i];
967       if ( delta > deltamax )
968         {
969           deltamax = delta ;
970         }
971     }
972   for (int i=0; i<dim; i++)
973     {
974       bbtemp[i*2]=bb1[i*2]-deltamax*eps;
975       bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
976     }
977   
978   for (int idim=0; idim < dim; idim++)
979     {
980       bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
981         && (bb2[idim*2]<bbtemp[idim*2+1]) ;
982       if (!intersects)
983         {
984           delete [] bbtemp;
985           return false; 
986         }
987     }
988   delete [] bbtemp;
989   return true;
990 }
991
992 /*!
993  * Intersect 2 given Bounding Boxes.
994  */
995 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
996 {
997   double* bbtemp(new double[2*dim]);
998   double deltamax(0.0);
999
1000   for (int i=0; i< dim; i++)
1001     {
1002       double delta = bb2[2*i+1]-bb2[2*i];
1003       if ( delta > deltamax )
1004         {
1005           deltamax = delta ;
1006         }
1007     }
1008   for (int i=0; i<dim; i++)
1009     {
1010       bbtemp[i*2]=bb2[i*2]-deltamax*eps;
1011       bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
1012     }
1013   
1014   bool intersects(!bb1.isDisjointWith(bbtemp));
1015   delete [] bbtemp;
1016   return intersects;
1017 }
1018
1019 /*!
1020  * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
1021  */
1022 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
1023 {
1024   double *coords(_coords->getPointer());
1025   int nbNodes(getNumberOfNodes());
1026   DataArrayDouble::Rotate3DAlg(center,vect,angle,nbNodes,coords,coords);
1027 }
1028
1029 /*!
1030  * This method allows to give for each cell in \a trgMesh, how much it interacts with cells of \a srcMesh.
1031  * The returned array can be seen as a weighted array on the target cells of \a trgMesh input parameter.
1032  *
1033  * \param [in] srcMesh - source mesh
1034  * \param [in] trgMesh - target mesh
1035  * \param [in] eps - precision of the detection
1036  * \return DataArrayInt * - An array that gives for each cell of \a trgMesh, how many cells in \a srcMesh (regarding the precision of detection \a eps) can interacts.
1037  * 
1038  * \throw If \a srcMesh and \a trgMesh have not the same space dimension.
1039  */
1040 DataArrayInt *MEDCouplingPointSet::ComputeNbOfInteractionsWithSrcCells(const MEDCouplingPointSet *srcMesh, const MEDCouplingPointSet *trgMesh, double eps)
1041 {
1042   if(!srcMesh || !trgMesh)
1043     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::ComputeNbOfInteractionsWithSrcCells : the input meshes must be not NULL !");
1044   MCAuto<DataArrayDouble> sbbox(srcMesh->getBoundingBoxForBBTree()),tbbox(trgMesh->getBoundingBoxForBBTree());
1045   return tbbox->computeNbOfInteractionsWith(sbbox,eps);
1046 }
1047
1048 /*!
1049  * Creates a new MEDCouplingMesh containing a part of cells of \a this mesh. The new
1050  * mesh shares a coordinates array with \a this one. The cells to include to the
1051  * result mesh are specified by an array of cell ids.
1052  *  \param [in] start - an array of cell ids to include to the result mesh.
1053  *  \param [in] end - specifies the end of the array \a start, so that
1054  *              the last value of \a start is \a end[ -1 ].
1055  *  \return MEDCouplingMesh * - a new instance of MEDCouplingMesh. The caller is to
1056  *         delete this mesh using decrRef() as it is no more needed. 
1057  */
1058 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
1059 {
1060   return buildPartOfMySelf(start,end,true);
1061 }
1062
1063 /*!
1064  * Creates a new MEDCouplingMesh containing a part of cells of \a this mesh. The
1065  * cells to include to the result mesh are specified by an array of cell ids. 
1066  * <br> This method additionally returns a renumbering map in "Old to New" mode
1067  * which allows the caller to know the mapping between nodes in \a this and the result mesh.
1068  *  \param [in] start - an array of cell ids to include to the result mesh.
1069  *  \param [in] end - specifies the end of the array \a start, so that
1070  *              the last value of \a start is \a end[ -1 ].
1071  *  \param [out] arr - a new DataArrayInt that is the "Old to New" renumbering
1072  *         map. The caller is to delete this array using decrRef() as it is no more needed.
1073  *  \return MEDCouplingMesh * - a new instance of MEDCouplingMesh. The caller is to
1074  *         delete this mesh using decrRef() as it is no more needed. 
1075  */
1076 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
1077 {
1078   MCAuto<MEDCouplingPointSet> ret=buildPartOfMySelf(start,end,true);
1079   arr=ret->zipCoordsTraducer();
1080   return ret.retn();
1081 }
1082
1083 /*!
1084  * This method specialized the MEDCouplingMesh::buildPartRange.
1085  * This method is equivalent to MEDCouplingMesh::buildPart method except that here the cell ids are specified using slice
1086  * \a beginCellIds \a endCellIds and \a stepCellIds.
1087  * \b WARNING , there is a big difference compared to MEDCouplingMesh::buildPart method.
1088  * If the input range is equal all cells in \a this, \a this is returned !
1089  *
1090  * \return a new ref to be managed by the caller. Warning this ref can be equal to \a this if input slice is exactly equal to the whole cells in the same order.
1091  *
1092  * \sa MEDCouplingUMesh::buildPartOfMySelfSlice
1093  */
1094 MEDCouplingMesh *MEDCouplingPointSet::buildPartRange(int beginCellIds, int endCellIds, int stepCellIds) const
1095 {
1096   if(beginCellIds==0 && endCellIds==getNumberOfCells() && stepCellIds==1)
1097     {
1098       MEDCouplingMesh *ret(const_cast<MEDCouplingPointSet *>(this));
1099       ret->incrRef();
1100       return ret;
1101     }
1102   else
1103     {
1104       return buildPartOfMySelfSlice(beginCellIds,endCellIds,stepCellIds,true);
1105     }
1106 }
1107
1108 /*!
1109  * This method specialized the MEDCouplingMesh::buildPartRangeAndReduceNodes
1110  *
1111  * \param [out] beginOut valid only if \a arr not NULL !
1112  * \param [out] endOut valid only if \a arr not NULL !
1113  * \param [out] stepOut valid only if \a arr not NULL !
1114  * \param [out] arr correspondance old to new in node ids.
1115  * 
1116  * \sa MEDCouplingUMesh::buildPartOfMySelfSlice
1117  */
1118 MEDCouplingMesh *MEDCouplingPointSet::buildPartRangeAndReduceNodes(int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt*& arr) const
1119 {
1120   MCAuto<MEDCouplingPointSet> ret(buildPartOfMySelfSlice(beginCellIds,endCellIds,stepCellIds,true));
1121   arr=ret->zipCoordsTraducer();
1122   return ret.retn();
1123 }
1124
1125 /*!
1126  * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
1127  */
1128 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
1129 {
1130   double *coords(_coords->getPointer());
1131   int nbNodes(getNumberOfNodes());
1132   DataArrayDouble::Rotate2DAlg(center,angle,nbNodes,coords,coords);
1133 }
1134
1135 /// @cond INTERNAL
1136
1137 class DummyClsMCPS
1138 {
1139 public:
1140   static const int MY_SPACEDIM=3;
1141   static const int MY_MESHDIM=2;
1142   typedef int MyConnType;
1143   static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
1144 };
1145
1146 /// @endcond
1147
1148 /*!
1149  * res should be an empty vector before calling this method.
1150  * This method returns all the node coordinates included in _coords which ids are in [startConn;endConn) and put it into 'res' vector.
1151  * If spaceDim==3 a projection will be done for each nodes on the middle plane containing these all nodes in [startConn;endConn).
1152  * And after each projected nodes are moved to Oxy plane in order to consider these nodes as 2D nodes.
1153  */
1154 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
1155 {
1156   const double *coords(_coords->getConstPointer());
1157   int spaceDim(getSpaceDimension());
1158   for(const int *it=startConn;it!=endConn;it++)
1159     res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
1160   if(spaceDim==2)
1161     return ;
1162   if(spaceDim==3)
1163     {
1164       std::vector<double> cpy(res);
1165       int nbNodes=(int)std::distance(startConn,endConn);
1166       INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::Projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0./*max distance*/,-1./*min dot*/,0.,true);
1167       res.resize(2*nbNodes);
1168       for(int i=0;i<nbNodes;i++)
1169         {
1170           res[2*i]=cpy[3*i];
1171           res[2*i+1]=cpy[3*i+1];
1172         }
1173       return ;
1174     }
1175   throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
1176 }
1177
1178 /*!
1179  * low level method that checks that the 2D cell is not a butterfly cell.
1180  */
1181 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
1182 {
1183   std::size_t nbOfNodes(res.size()/2);
1184   std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
1185   for(std::size_t i=0;i<nbOfNodes;i++)
1186     {
1187       INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
1188       nodes[i]=tmp;
1189     }
1190   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1191   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1192   INTERP_KERNEL::QuadraticPolygon *pol=0;
1193   if(isQuad)
1194     pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
1195   else
1196     pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
1197   bool ret(pol->isButterflyAbs());
1198   delete pol;
1199   return ret;
1200 }
1201
1202 /*!
1203  * This method compares 2 cells coming from two unstructured meshes : \a this and \a other.
1204  * This method compares 2 cells having the same id 'cellId' in \a this and \a other.
1205  */
1206 bool MEDCouplingPointSet::areCellsFrom2MeshEqual(const MEDCouplingPointSet *other, int cellId, double prec) const
1207 {
1208   if(getTypeOfCell(cellId)!=other->getTypeOfCell(cellId))
1209     return false;
1210   std::vector<int> c1,c2;
1211   getNodeIdsOfCell(cellId,c1);
1212   other->getNodeIdsOfCell(cellId,c2);
1213   std::size_t sz(c1.size());
1214   if(sz!=c2.size())
1215     return false;
1216   for(std::size_t i=0;i<sz;i++)
1217     {
1218       std::vector<double> n1,n2;
1219       getCoordinatesOfNode(c1[0],n1);
1220       other->getCoordinatesOfNode(c2[0],n2);
1221       std::transform(n1.begin(),n1.end(),n2.begin(),n1.begin(),std::minus<double>());
1222       std::transform(n1.begin(),n1.end(),n1.begin(),std::ptr_fun<double,double>(fabs));
1223       if(*std::max_element(n1.begin(),n1.end())>prec)
1224         return false;
1225     }
1226   return true;
1227 }
1228
1229 /*!
1230  * Substitutes node coordinates array of \a this mesh with that of \a other mesh
1231  * (i.e. \a this->_coords with \a other._coords) provided that coordinates of the two
1232  * meshes match with a specified precision, else an exception is thrown and \a this
1233  * remains unchanged. In case of success the nodal connectivity of \a this mesh
1234  * is permuted according to new order of nodes.
1235  * Contrary to tryToShareSameCoords() this method makes a deeper analysis of
1236  * coordinates (and so more expensive) than simple equality.
1237  *  \param [in] other - the other mesh whose node coordinates array will be used by
1238  *         \a this mesh in case of their equality.
1239  *  \param [in] epsilon - the precision used to compare coordinates (using infinite norm).
1240  *  \throw If the coordinates array of \a this is not set.
1241  *  \throw If the coordinates array of \a other is not set.
1242  *  \throw If the coordinates of \a this and \a other do not match.
1243  */
1244 void MEDCouplingPointSet::tryToShareSameCoordsPermute(const MEDCouplingPointSet& other, double epsilon)
1245 {
1246   const DataArrayDouble *coords=other.getCoords();
1247   if(!coords)
1248     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute : No coords specified in other !");
1249   if(!_coords)
1250     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute : No coords specified in this whereas there is any in other !");
1251   int otherNbOfNodes=other.getNumberOfNodes();
1252   MCAuto<DataArrayDouble> newCoords=MergeNodesArray(&other,this);
1253   _coords->incrRef();
1254   MCAuto<DataArrayDouble> oldCoords=_coords;
1255   setCoords(newCoords);
1256   bool areNodesMerged;
1257   int newNbOfNodes;
1258   MCAuto<DataArrayInt> da=buildPermArrayForMergeNode(epsilon,otherNbOfNodes,areNodesMerged,newNbOfNodes);
1259   if(!areNodesMerged)
1260     {
1261       setCoords(oldCoords);
1262       throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute fails : no nodes are mergeable with specified given epsilon !");
1263     }
1264   int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+otherNbOfNodes);
1265   const int *pt=std::find_if(da->getConstPointer()+otherNbOfNodes,da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
1266   if(pt!=da->getConstPointer()+da->getNbOfElems())
1267     {
1268       setCoords(oldCoords);
1269       throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute fails : some nodes in this are not in other !");
1270     }
1271   setCoords(oldCoords);
1272   renumberNodesInConn(da->getConstPointer()+otherNbOfNodes);
1273   setCoords(coords);
1274 }
1275
1276 MEDCouplingPointSet *MEDCouplingPointSet::buildPartOfMySelf(const int *begin, const int *end, bool keepCoords) const
1277 {
1278   MCAuto<MEDCouplingPointSet> ret=buildPartOfMySelfKeepCoords(begin,end);
1279   if(!keepCoords)
1280     ret->zipCoords();
1281   return ret.retn();
1282 }
1283
1284 MEDCouplingPointSet *MEDCouplingPointSet::buildPartOfMySelfSlice(int start, int end, int step, bool keepCoords) const
1285 {
1286   MCAuto<MEDCouplingPointSet> ret=buildPartOfMySelfKeepCoordsSlice(start,end,step);
1287   if(!keepCoords)
1288     ret->zipCoords();
1289   return ret.retn();
1290 }
1291
1292 /*!
1293  Creates a new MEDCouplingUMesh containing some cells of \a this mesh. The cells to
1294  copy are selected basing on specified node ids and the value of \a fullyIn
1295  parameter. If \a fullyIn ==\c true, a cell is copied if its all nodes are in the 
1296  array \a begin of node ids. If \a fullyIn ==\c false, a cell is copied if any its
1297  node is in the array of node ids. The created mesh shares the node coordinates array
1298  with \a this mesh.
1299  *  \param [in] begin - the array of node ids.
1300  *  \param [in] end - a pointer to the (last+1)-th element of \a begin.
1301  *  \param [in] fullyIn - if \c true, then cells whose all nodes are in the
1302  *         array \a begin are copied, else cells whose any node is in the
1303  *         array \a begin are copied.
1304  *  \return MEDCouplingPointSet * - new instance of MEDCouplingUMesh. The caller is
1305  *         to delete this mesh using decrRef() as it is no more needed. 
1306  *  \throw If the coordinates array is not set.
1307  *  \throw If the nodal connectivity of cells is not defined.
1308  *  \throw If any node id in \a begin is not valid.
1309  *
1310  *  \if ENABLE_EXAMPLES
1311  *  \ref cpp_mcumesh_buildPartOfMySelfNode "Here is a C++ example".<br>
1312  *  \ref  py_mcumesh_buildPartOfMySelfNode "Here is a Python example".
1313  *  \endif
1314  */
1315 MEDCouplingPointSet *MEDCouplingPointSet::buildPartOfMySelfNode(const int *begin, const int *end, bool fullyIn) const
1316 {
1317   DataArrayInt *cellIdsKept=0;
1318   fillCellIdsToKeepFromNodeIds(begin,end,fullyIn,cellIdsKept);
1319   MCAuto<DataArrayInt> cellIdsKept2(cellIdsKept);
1320   return buildPartOfMySelf(cellIdsKept->begin(),cellIdsKept->end(),true);
1321 }
1322
1323 /*!
1324  * Removes duplicates of cells from \a this mesh and returns an array mapping between
1325  * new and old cell ids in "Old to New" mode. Nothing is changed in \a this mesh if no
1326  * equal cells found.
1327  *  \warning Cells of the result mesh are \b not sorted by geometric type, hence,
1328  *           to write this mesh to the MED file, its cells must be sorted using
1329  *           sortCellsInMEDFileFrmt().
1330  *  \param [in] compType - specifies a cell comparison technique. Meaning of its
1331  *          valid values [0,1,2] is as follows.
1332  *   - 0 : "exact". Two cells are considered equal \c iff they have exactly same nodal
1333  *         connectivity and type. This is the strongest policy.
1334  *   - 1 : "permuted same orientation". Two cells are considered equal \c iff they
1335  *         are based on same nodes and have the same type and orientation.
1336  *   - 2 : "nodal". Two cells are considered equal \c iff they
1337  *         are based on same nodes and have the same type. This is the weakest
1338  *         policy, it can be used by users not sensitive to cell orientation.
1339  *  \param [in] startCellId - specifies the cell id at which search for equal cells
1340  *         starts. By default it is 0, which means that all cells in \a this will be
1341  *         scanned. 
1342  *  \return DataArrayInt - a new instance of DataArrayInt, of length \a
1343  *           this->getNumberOfCells() before call of this method. The caller is to
1344  *           delete this array using decrRef() as it is no more needed. 
1345  *  \throw If the coordinates array is not set.
1346  *  \throw If the nodal connectivity of cells is not defined.
1347  *  \throw If the nodal connectivity includes an invalid id.
1348  *
1349  *  \if ENABLE_EXAMPLES
1350  *  \ref cpp_mcumesh_zipConnectivityTraducer "Here is a C++ example".<br>
1351  *  \ref  py_mcumesh_zipConnectivityTraducer "Here is a Python example".
1352  *  \endif
1353  */
1354 DataArrayInt *MEDCouplingPointSet::zipConnectivityTraducer(int compType, int startCellId)
1355 {
1356   DataArrayInt *commonCells=0,*commonCellsI=0;
1357   findCommonCells(compType,startCellId,commonCells,commonCellsI);
1358   MCAuto<DataArrayInt> commonCellsTmp(commonCells),commonCellsITmp(commonCellsI);
1359   int newNbOfCells=-1;
1360   MCAuto<DataArrayInt> ret=DataArrayInt::ConvertIndexArrayToO2N(getNumberOfCells(),commonCells->begin(),commonCellsI->begin(),
1361                                                                                                           commonCellsI->end(),newNbOfCells);
1362   MCAuto<DataArrayInt> ret2=ret->invertArrayO2N2N2O(newNbOfCells);
1363   MCAuto<MEDCouplingPointSet> self=buildPartOfMySelf(ret2->begin(),ret2->end(),true);
1364   shallowCopyConnectivityFrom(self);
1365   return ret.retn();
1366 }
1367
1368 /*!
1369  * This const method states if the nodal connectivity of this fetches all nodes in \a this.
1370  * In other words, this method looks is there are no orphan nodes in \a this.
1371  * \sa zipCoordsTraducer, getNodeIdsInUse, computeFetchedNodeIds.
1372  */
1373 bool MEDCouplingPointSet::areAllNodesFetched() const
1374 {
1375   checkFullyDefined();
1376   int nbNodes(getNumberOfNodes());
1377   std::vector<bool> fetchedNodes(nbNodes,false);
1378   computeNodeIdsAlg(fetchedNodes);
1379   return std::find(fetchedNodes.begin(),fetchedNodes.end(),false)==fetchedNodes.end();
1380 }
1381
1382 /*!
1383  * Checks if \a this and \a other meshes are geometrically equivalent, else an
1384  * exception is thrown. The meshes are
1385  * considered equivalent if (1) \a this mesh contains the same nodes as the \a other
1386  * mesh (with a specified precision) and (2) \a this mesh contains the same cells as
1387  * the \a other mesh (with use of a specified cell comparison technique). The mapping 
1388  * from \a other to \a this for nodes and cells is returned via out parameters.
1389  *  \param [in] other - the mesh to compare with.
1390  *  \param [in] cellCompPol - id [0-2] of cell comparison method. See meaning of
1391  *         each method in description of MEDCouplingPointSet::zipConnectivityTraducer().
1392  *  \param [in] prec - the precision used to compare nodes of the two meshes.
1393  *  \param [out] cellCor - a cell permutation array in "Old to New" mode. The caller is
1394  *         to delete this array using decrRef() as it is no more needed.
1395  *  \param [out] nodeCor - a node permutation array in "Old to New" mode. The caller is
1396  *         to delete this array using decrRef() as it is no more needed.
1397  *  \throw If the two meshes do not match.
1398  *
1399  *  \if ENABLE_EXAMPLES
1400  *  \ref cpp_mcumesh_checkDeepEquivalWith "Here is a C++ example".<br>
1401  *  \ref  py_mcumesh_checkDeepEquivalWith "Here is a Python example".
1402  *  \endif
1403  */
1404 void MEDCouplingPointSet::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
1405                                                DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const
1406 {
1407   if(!other)
1408     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalWith : input is null !");
1409   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
1410   if(!otherC)
1411     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalWith : other is not a PointSet mesh !");
1412   MCAuto<MEDCouplingPointSet> m=dynamic_cast<MEDCouplingPointSet *>(mergeMyselfWith(otherC));
1413   bool areNodesMerged;
1414   int newNbOfNodes;
1415   int oldNbOfNodes=getNumberOfNodes();
1416   MCAuto<DataArrayInt> da=m->buildPermArrayForMergeNode(prec,oldNbOfNodes,areNodesMerged,newNbOfNodes);
1417   //mergeNodes
1418   if(!areNodesMerged && oldNbOfNodes != 0)
1419     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Nodes are incompatible ! ");
1420   const int *pt=std::find_if(da->getConstPointer()+oldNbOfNodes,da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),oldNbOfNodes-1));
1421   if(pt!=da->getConstPointer()+da->getNbOfElems())
1422     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some nodes in other are not in this !");
1423   m->renumberNodes(da->getConstPointer(),newNbOfNodes);
1424   //
1425   MCAuto<DataArrayInt> nodeCor2=da->subArray(oldNbOfNodes);
1426   da=m->mergeNodes(prec,areNodesMerged,newNbOfNodes);
1427   //
1428   da=m->zipConnectivityTraducer(cellCompPol);
1429   int nbCells=getNumberOfCells();
1430   if (nbCells != other->getNumberOfCells())
1431     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some cells in other are not in this !");
1432   int dan(da->getNumberOfTuples());
1433   if (dan)
1434     {
1435       MCAuto<DataArrayInt> da1(DataArrayInt::New()),da2(DataArrayInt::New());
1436       da1->alloc(dan/2,1); da2->alloc(dan/2,1);
1437       std::copy(da->getConstPointer(), da->getConstPointer()+dan/2, da1->getPointer());
1438       std::copy(da->getConstPointer()+dan/2, da->getConstPointer()+dan, da2->getPointer());
1439       da1->sort(); da2->sort();
1440       if (!da1->isEqualWithoutConsideringStr(*da2))
1441         throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some cells in other are not in this !");
1442     }
1443   MCAuto<DataArrayInt> cellCor2=da->selectByTupleIdSafeSlice(nbCells,da->getNbOfElems(),1);
1444   nodeCor=nodeCor2->isIota(nodeCor2->getNumberOfTuples())?0:nodeCor2.retn();
1445   cellCor=cellCor2->isIota(cellCor2->getNumberOfTuples())?0:cellCor2.retn();
1446 }
1447
1448 /*!
1449  * Checks if \a this and \a other meshes are geometrically equivalent, else an
1450  * exception is thrown. The meshes are considered equivalent if (1) they share one
1451  * node coordinates array and (2) they contain the same cells (with use of a specified
1452  * cell comparison technique). The mapping from cells of the \a other to ones of \a this 
1453  * is returned via an out parameter.
1454  *  \param [in] other - the mesh to compare with.
1455  *  \param [in] cellCompPol - id [0-2] of cell comparison method. See the meaning of
1456  *         each method in description of MEDCouplingPointSet::zipConnectivityTraducer().
1457  *  \param [in] prec - a not used parameter.
1458  *  \param [out] cellCor - the permutation array in "Old to New" mode. The caller is
1459  *         to delete this array using decrRef() as it is no more needed.
1460  *  \throw If the two meshes do not match.
1461  *
1462  * \if ENABLE_EXAMPLES
1463  * \ref cpp_mcumesh_checkDeepEquivalWith "Here is a C++ example".<br>
1464  * \ref  py_mcumesh_checkDeepEquivalWith "Here is a Python example".
1465  * \endif
1466  */
1467 void MEDCouplingPointSet::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
1468                                                        DataArrayInt *&cellCor) const
1469 {
1470   if(!other)
1471     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalOnSameNodesWith : input is null !");
1472   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
1473   if(!otherC)
1474     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalOnSameNodesWith : other is not a PointSet mesh !");
1475   if(_coords!=otherC->_coords)
1476     throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : meshes do not share the same coordinates ! Use tryToShareSameCoordinates or call checkDeepEquivalWith !");
1477   MCAuto<MEDCouplingPointSet> m=mergeMyselfWithOnSameCoords(otherC);
1478   MCAuto<DataArrayInt> da=m->zipConnectivityTraducer(cellCompPol);
1479   int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+getNumberOfCells());
1480   const int *pt=std::find_if(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
1481   if(pt!=da->getConstPointer()+da->getNbOfElems())
1482     {
1483       throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : some cells in other are not in this !");
1484     }
1485   MCAuto<DataArrayInt> cellCor2=da->selectByTupleIdSafeSlice(getNumberOfCells(),da->getNbOfElems(),1);
1486   cellCor=cellCor2->isIota(cellCor2->getNumberOfTuples())?0:cellCor2.retn();
1487 }
1488
1489 void MEDCouplingPointSet::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const
1490 {
1491   MEDCouplingMesh::checkFastEquivalWith(other,prec);
1492   //other not null checked by the line before
1493   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
1494   if(!otherC)
1495     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkFastEquivalWith : fails because other is not a pointset mesh !");
1496   int nbOfCells=getNumberOfCells();
1497   if(nbOfCells<1)
1498     return ;
1499   bool status=true;
1500   status&=areCellsFrom2MeshEqual(otherC,0,prec);
1501   status&=areCellsFrom2MeshEqual(otherC,nbOfCells/2,prec);
1502   status&=areCellsFrom2MeshEqual(otherC,nbOfCells-1,prec);
1503   if(!status)
1504     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Two meshes are not equal because on 3 test cells some difference have been detected !");
1505 }
1506
1507 /*!
1508  * Finds cells whose all or some nodes are in a given array of node ids.
1509  *  \param [in] begin - the array of node ids.
1510  *  \param [in] end - a pointer to the (last+1)-th element of \a begin.
1511  *  \param [in] fullyIn - if \c true, then cells whose all nodes are in the
1512  *         array \a begin are returned only, else cells whose any node is in the
1513  *         array \a begin are returned.
1514  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of found
1515  *         cells. The caller is to delete this array using decrRef() as it is no more
1516  *         needed. 
1517  *  \throw If the coordinates array is not set.
1518  *  \throw If the nodal connectivity of cells is not defined.
1519  *  \throw If any cell id in \a begin is not valid.
1520  *
1521  * \sa MEDCouplingPointSet::getCellIdsFullyIncludedInNodeIds
1522  *
1523  *  \if ENABLE_EXAMPLES
1524  *  \ref cpp_mcumesh_getCellIdsLyingOnNodes "Here is a C++ example".<br>
1525  *  \ref  py_mcumesh_getCellIdsLyingOnNodes "Here is a Python example".
1526  *  \endif
1527  */
1528 DataArrayInt *MEDCouplingPointSet::getCellIdsLyingOnNodes(const int *begin, const int *end, bool fullyIn) const
1529 {
1530   DataArrayInt *cellIdsKept=0;
1531   fillCellIdsToKeepFromNodeIds(begin,end,fullyIn,cellIdsKept);
1532   cellIdsKept->setName(getName());
1533   return cellIdsKept;
1534 }
1535
1536 /*!
1537  * Finds cells whose all nodes are in a given array of node ids.
1538  * This method is a specialization of MEDCouplingPointSet::getCellIdsLyingOnNodes (true
1539  * as last input argument).
1540  *  \param [in] partBg - the array of node ids.
1541  *  \param [in] partEnd - a pointer to a (last+1)-th element of \a partBg.
1542  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of found
1543  *          cells. The caller is to delete this array using decrRef() as it is no
1544  *          more needed.
1545  *  \throw If the coordinates array is not set.
1546  *  \throw If the nodal connectivity of cells is not defined.
1547  *  \throw If any cell id in \a partBg is not valid.
1548  * 
1549  * \sa MEDCouplingPointSet::getCellIdsLyingOnNodes
1550  *
1551  *  \if ENABLE_EXAMPLES
1552  *  \ref cpp_mcumesh_getCellIdsFullyIncludedInNodeIds "Here is a C++ example".<br>
1553  *  \ref  py_mcumesh_getCellIdsFullyIncludedInNodeIds "Here is a Python example".
1554  *  \endif
1555  */
1556 DataArrayInt *MEDCouplingPointSet::getCellIdsFullyIncludedInNodeIds(const int *partBg, const int *partEnd) const
1557 {
1558   return getCellIdsLyingOnNodes(partBg,partEnd,true);
1559 }
1560
1561 /*!
1562  * Removes unused nodes (the node coordinates array is shorten) and returns an array
1563  * mapping between new and old node ids in "Old to New" mode. -1 values in the returned
1564  * array mean that the corresponding old node is no more used. 
1565  *  \return DataArrayInt * - a new instance of DataArrayInt of length \a
1566  *           this->getNumberOfNodes() before call of this method. The caller is to
1567  *           delete this array using decrRef() as it is no more needed. 
1568  *  \throw If the coordinates array is not set.
1569  *  \throw If the nodal connectivity of cells is not defined.
1570  *  \throw If the nodal connectivity includes an invalid id.
1571  *
1572  *  \if ENABLE_EXAMPLES
1573  *  \ref cpp_mcumesh_zipCoordsTraducer "Here is a C++ example".<br>
1574  *  \ref  py_mcumesh_zipCoordsTraducer "Here is a Python example".
1575  *  \endif
1576  */
1577 DataArrayInt *MEDCouplingPointSet::zipCoordsTraducer()
1578 {
1579   int newNbOfNodes=-1;
1580   MCAuto<DataArrayInt> traducer=getNodeIdsInUse(newNbOfNodes);
1581   renumberNodes(traducer->getConstPointer(),newNbOfNodes);
1582   return traducer.retn();
1583 }
1584
1585 /*!
1586  * Merges nodes equal within \a precision and returns an array describing the 
1587  * permutation used to remove duplicate nodes.
1588  *  \param [in] precision - minimal absolute distance between two nodes at which they are
1589  *              considered not coincident.
1590  *  \param [out] areNodesMerged - is set to \c true if any coincident nodes removed.
1591  *  \param [out] newNbOfNodes - number of nodes remaining after the removal.
1592  *  \return DataArrayInt * - the permutation array in "Old to New" mode. For more 
1593  *          info on "Old to New" mode see \ref numbering. The caller
1594  *          is to delete this array using decrRef() as it is no more needed.
1595  *  \throw If the coordinates array is not set.
1596  *  \throw If the nodal connectivity of cells is not defined.
1597  *
1598  *  \if ENABLE_EXAMPLES
1599  *  \ref cpp_mcumesh_mergeNodes "Here is a C++ example".<br>
1600  *  \ref  py_mcumesh_mergeNodes "Here is a Python example".
1601  *  \endif
1602  */
1603 DataArrayInt *MEDCouplingPointSet::mergeNodes(double precision, bool& areNodesMerged, int& newNbOfNodes)
1604 {
1605   MCAuto<DataArrayInt> ret=buildPermArrayForMergeNode(precision,-1,areNodesMerged,newNbOfNodes);
1606   if(areNodesMerged)
1607     renumberNodes(ret->begin(),newNbOfNodes);
1608   return ret.retn();
1609 }
1610
1611 /*!
1612  * Merges nodes equal within \a precision and returns an array describing the 
1613  * permutation used to remove duplicate nodes. In contrast to mergeNodes(), location
1614  *  of merged nodes is changed to be at their barycenter.
1615  *  \param [in] precision - minimal absolute distance between two nodes at which they are
1616  *              considered not coincident.
1617  *  \param [out] areNodesMerged - is set to \c true if any coincident nodes removed.
1618  *  \param [out] newNbOfNodes - number of nodes remaining after the removal.
1619  *  \return DataArrayInt * - the permutation array in "Old to New" mode. For more 
1620  *          info on "Old to New" mode see \ref numbering. The caller
1621  *          is to delete this array using decrRef() as it is no more needed.
1622  *  \throw If the coordinates array is not set.
1623  *  \throw If the nodal connectivity of cells is not defined.
1624  *
1625  *  \if ENABLE_EXAMPLES
1626  *  \ref cpp_mcumesh_mergeNodes "Here is a C++ example".<br>
1627  *  \ref  py_mcumesh_mergeNodes "Here is a Python example".
1628  *  \endif
1629  */
1630 DataArrayInt *MEDCouplingPointSet::mergeNodesCenter(double precision, bool& areNodesMerged, int& newNbOfNodes)
1631 {
1632   DataArrayInt *ret=buildPermArrayForMergeNode(precision,-1,areNodesMerged,newNbOfNodes);
1633   if(areNodesMerged)
1634     renumberNodesCenter(ret->getConstPointer(),newNbOfNodes);
1635   return ret;
1636 }