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