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