src/MEDCoupling/MEDCouplingPointSet.cxx

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