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   renumberNodesInConn(newNodeNumbers);
 393   setCoords(newCoords);//let it here not before renumberNodesInConn because old number of nodes is sometimes used...
 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   if(!coo[0]->isAllocated())
 786     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : first element in coordinates is not allocated !");
 787   for(int i=1;it!=ms.end();it++,i++)
 788     {
 789       const DataArrayDouble *tmp=(*it)->getCoords();
 790       if(tmp)
 791         {
 792           if(tmp->isAllocated())
 793             {
 794               if((*it)->getSpaceDimension()==spaceDim)
 795                 coo[i]=tmp;
 796               else
 797                 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : Mismatch in SpaceDim !");
 798             }
 799           else
 800             throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : Presence of a non allocated array !");
 801         }
 802       else
 803         throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : Empty coords detected !");
 804     }
 805   return DataArrayDouble::Aggregate(coo);
 806 }
 807
 808 /*!
 809  * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
 810  * This method is used during unserialization process.
 811  */
 812 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
 813 {
 814   switch(type)
 815     {
 816     case UNSTRUCTURED:
 817       return MEDCouplingUMesh::New();
 818     case UNSTRUCTURED_DESC:
 819       return MEDCouplingUMeshDesc::New();
 820     default:
 821       throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
 822     }
 823 }
 824
 825 /*!
 826  * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
 827  */
 828 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 829 {
 830   int it,order;
 831   double time=getTime(it,order);
 832   if(_coords)
 833     {
 834       int spaceDim=getSpaceDimension();
 835       littleStrings.resize(spaceDim+4);
 836       littleStrings[0]=getName();
 837       littleStrings[1]=getDescription();
 838       littleStrings[2]=_coords->getName();
 839       littleStrings[3]=getTimeUnit();
 840       for(int i=0;i<spaceDim;i++)
 841         littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
 842       tinyInfo.clear();
 843       tinyInfo.push_back(getType());
 844       tinyInfo.push_back(spaceDim);
 845       tinyInfo.push_back(getNumberOfNodes());
 846       tinyInfo.push_back(it);
 847       tinyInfo.push_back(order);
 848       tinyInfoD.push_back(time);
 849     }
 850   else
 851     {
 852       littleStrings.resize(3);
 853       littleStrings[0]=getName();
 854       littleStrings[1]=getDescription();
 855       littleStrings[2]=getTimeUnit();
 856       tinyInfo.clear();
 857       tinyInfo.push_back(getType());
 858       tinyInfo.push_back(-1);
 859       tinyInfo.push_back(-1);
 860       tinyInfo.push_back(it);
 861       tinyInfo.push_back(order);
 862       tinyInfoD.push_back(time);
 863     }
 864 }
 865
 866 /*!
 867  * Third and final step of serialization process.
 868  */
 869 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 870 {
 871   if(_coords)
 872     {
 873       a2=const_cast<DataArrayDouble *>(getCoords());
 874       a2->incrRef();
 875     }
 876   else
 877     a2=0;
 878 }
 879
 880 /*!
 881  * Second step of serialization process.
 882  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 883  */
 884 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 885 {
 886   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 887     {
 888       a2->alloc(tinyInfo[2],tinyInfo[1]);
 889       littleStrings.resize(tinyInfo[1]+4);
 890     }
 891   else
 892     {
 893       littleStrings.resize(3);
 894     }
 895 }
 896
 897 /*!
 898  * Second and final unserialization process.
 899  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 900  */
 901 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
 902 {
 903   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 904     {
 905       setCoords(a2);
 906       setName(littleStrings[0].c_str());
 907       setDescription(littleStrings[1].c_str());
 908       a2->setName(littleStrings[2].c_str());
 909       setTimeUnit(littleStrings[3].c_str());
 910       for(int i=0;i<tinyInfo[1];i++)
 911         getCoords()->setInfoOnComponent(i,littleStrings[i+4].c_str());
 912       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 913     }
 914   else
 915     {
 916       setName(littleStrings[0].c_str());
 917       setDescription(littleStrings[1].c_str());
 918       setTimeUnit(littleStrings[2].c_str());
 919       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 920     }
 921 }
 922
 923 void MEDCouplingPointSet::checkCoherency() const throw(INTERP_KERNEL::Exception)
 924 {
 925   if(!_coords)
 926     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkCoherency : no coordinates set !");
 927 }
 928
 929 /*!
 930  * Intersect Bounding Box given 2 Bounding Boxes.
 931  */
 932 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
 933 {
 934   double* bbtemp = new double[2*dim];
 935   double deltamax=0.0;
 936
 937   for (int i=0; i< dim; i++)
 938     {
 939       double delta = bb1[2*i+1]-bb1[2*i];
 940       if ( delta > deltamax )
 941         {
 942           deltamax = delta ;
 943         }
 944     }
 945   for (int i=0; i<dim; i++)
 946     {
 947       bbtemp[i*2]=bb1[i*2]-deltamax*eps;
 948       bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
 949     }
 950
 951   for (int idim=0; idim < dim; idim++)
 952     {
 953       bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
 954         && (bb2[idim*2]<bbtemp[idim*2+1]) ;
 955       if (!intersects)
 956         {
 957           delete [] bbtemp;
 958           return false;
 959         }
 960     }
 961   delete [] bbtemp;
 962   return true;
 963 }
 964
 965 /*!
 966  * Intersect 2 given Bounding Boxes.
 967  */
 968 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
 969 {
 970   double* bbtemp = new double[2*dim];
 971   double deltamax=0.0;
 972
 973   for (int i=0; i< dim; i++)
 974     {
 975       double delta = bb2[2*i+1]-bb2[2*i];
 976       if ( delta > deltamax )
 977         {
 978           deltamax = delta ;
 979         }
 980     }
 981   for (int i=0; i<dim; i++)
 982     {
 983       bbtemp[i*2]=bb2[i*2]-deltamax*eps;
 984       bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
 985     }
 986
 987   bool intersects = !bb1.isDisjointWith( bbtemp );
 988   delete [] bbtemp;
 989   return intersects;
 990 }
 991
 992 /*!
 993  * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
 994  */
 995 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
 996 {
 997   double *coords=_coords->getPointer();
 998   int nbNodes=getNumberOfNodes();
 999   Rotate3DAlg(center,vect,angle,nbNodes,coords);
1000 }
1001
1002 /*!
1003  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
1004  * around an axe ('center','vect') and with angle 'angle'.
1005  */
1006 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
1007 {
1008   if(!center || !vect)
1009     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : null vector in input !");
1010   double sina=sin(angle);
1011   double cosa=cos(angle);
1012   double vectorNorm[3];
1013   double matrix[9];
1014   double matrixTmp[9];
1015   double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
1016   if(norm<std::numeric_limits<double>::min())
1017     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : magnitude of input vector is too close of 0. !");
1018   std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
1019   //rotation matrix computation
1020   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;
1021   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
1022   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
1023   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
1024   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
1025   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
1026   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
1027   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
1028   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
1029   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
1030   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
1031   //rotation matrix computed.
1032   double tmp[3];
1033   for(int i=0; i<nbNodes; i++)
1034     {
1035       std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
1036       coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
1037       coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
1038       coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
1039     }
1040 }
1041
1042 /*!
1043  * Creates a new MEDCouplingMesh containing a part of cells of \a this mesh. The new
1044  * mesh shares a coordinates array with \a this one. The cells to include to the
1045  * result mesh are specified by an array of cell ids.
1046  *  \param [in] start - an array of cell ids to include to the result mesh.
1047  *  \param [in] end - specifies the end of the array \a start, so that
1048  *              the last value of \a start is \a end[ -1 ].
1049  *  \return MEDCouplingMesh * - a new instance of MEDCouplingMesh. The caller is to
1050  *         delete this mesh using decrRef() as it is no more needed.
1051  */
1052 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
1053 {
1054   return buildPartOfMySelf(start,end,true);
1055 }
1056
1057 /*!
1058  * Creates a new MEDCouplingMesh containing a part of cells of \a this mesh. The
1059  * cells to include to the result mesh are specified by an array of cell ids.
1060  * <br> This method additionally returns a renumbering map in "Old to New" mode
1061  * which allows the caller to know the mapping between nodes in \a this and the result mesh.
1062  *  \param [in] start - an array of cell ids to include to the result mesh.
1063  *  \param [in] end - specifies the end of the array \a start, so that
1064  *              the last value of \a start is \a end[ -1 ].
1065  *  \param [out] arr - a new DataArrayInt that is the "Old to New" renumbering
1066  *         map. The caller is to delete this array using decrRef() as it is no more needed.
1067  *  \return MEDCouplingMesh * - a new instance of MEDCouplingMesh. The caller is to
1068  *         delete this mesh using decrRef() as it is no more needed.
1069  */
1070 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
1071 {
1072   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> ret=buildPartOfMySelf(start,end,true);
1073   arr=ret->zipCoordsTraducer();
1074   return ret.retn();
1075 }
1076
1077 /*!
1078  * This method specialized the MEDCouplingMesh::buildPartRange
1079  *
1080  * \sa MEDCouplingUMesh::buildPartOfMySelf2
1081  */
1082 MEDCouplingMesh *MEDCouplingPointSet::buildPartRange(int beginCellIds, int endCellIds, int stepCellIds) const throw(INTERP_KERNEL::Exception)
1083 {
1084   return buildPartOfMySelf2(beginCellIds,endCellIds,stepCellIds,true);
1085 }
1086
1087 /*!
1088  * This method specialized the MEDCouplingMesh::buildPartRangeAndReduceNodes
1089  *
1090  * \param [out] beginOut valid only if \a arr not NULL !
1091  * \param [out] endOut valid only if \a arr not NULL !
1092  * \param [out] stepOut valid only if \a arr not NULL !
1093  * \param [out] arr correspondance old to new in node ids.
1094  *
1095  * \sa MEDCouplingUMesh::buildPartOfMySelf2
1096  */
1097 MEDCouplingMesh *MEDCouplingPointSet::buildPartRangeAndReduceNodes(int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt*& arr) const throw(INTERP_KERNEL::Exception)
1098 {
1099   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> ret=buildPartOfMySelf2(beginCellIds,endCellIds,stepCellIds,true);
1100   arr=ret->zipCoordsTraducer();
1101   return ret.retn();
1102 }
1103
1104 /*!
1105  * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
1106  */
1107 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
1108 {
1109   double *coords=_coords->getPointer();
1110   int nbNodes=getNumberOfNodes();
1111   Rotate2DAlg(center,angle,nbNodes,coords);
1112 }
1113
1114 /*!
1115  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
1116  * around the center point 'center' and with angle 'angle'.
1117  */
1118 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
1119 {
1120   double cosa=cos(angle);
1121   double sina=sin(angle);
1122   double matrix[4];
1123   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
1124   double tmp[2];
1125   for(int i=0; i<nbNodes; i++)
1126     {
1127       std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
1128       coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
1129       coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
1130     }
1131 }
1132
1133 /// @cond INTERNAL
1134
1135 class DummyClsMCPS
1136 {
1137 public:
1138   static const int MY_SPACEDIM=3;
1139   static const int MY_MESHDIM=2;
1140   typedef int MyConnType;
1141   static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
1142 };
1143
1144 /// @endcond
1145
1146 /*!
1147  * res should be an empty vector before calling this method.
1148  * This method returns all the node coordinates included in _coords which ids are in [startConn;endConn) and put it into 'res' vector.
1149  * If spaceDim==3 a projection will be done for each nodes on the middle plane containing these all nodes in [startConn;endConn).
1150  * And after each projected nodes are moved to Oxy plane in order to consider these nodes as 2D nodes.
1151  */
1152 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
1153 {
1154   const double *coords=_coords->getConstPointer();
1155   int spaceDim=getSpaceDimension();
1156   for(const int *it=startConn;it!=endConn;it++)
1157     res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
1158   if(spaceDim==2)
1159     return ;
1160   if(spaceDim==3)
1161     {
1162       std::vector<double> cpy(res);
1163       int nbNodes=(int)std::distance(startConn,endConn);
1164       INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0.,0.,true);
1165       res.resize(2*nbNodes);
1166       for(int i=0;i<nbNodes;i++)
1167         {
1168           res[2*i]=cpy[3*i];
1169           res[2*i+1]=cpy[3*i+1];
1170         }
1171       return ;
1172     }
1173   throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
1174 }
1175
1176 /*!
1177  * low level method that checks that the 2D cell is not a butterfly cell.
1178  */
1179 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
1180 {
1181   std::size_t nbOfNodes=res.size()/2;
1182   std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
1183   for(std::size_t i=0;i<nbOfNodes;i++)
1184     {
1185       INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
1186       nodes[i]=tmp;
1187     }
1188   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1189   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1190   INTERP_KERNEL::QuadraticPolygon *pol=0;
1191   if(isQuad)
1192     pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
1193   else
1194     pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
1195   bool ret=pol->isButterflyAbs();
1196   delete pol;
1197   return ret;
1198 }
1199
1200 /*!
1201  * This method compares 2 cells coming from two unstructured meshes : \a this and \a other.
1202  * This method compares 2 cells having the same id 'cellId' in \a this and \a other.
1203  */
1204 bool MEDCouplingPointSet::areCellsFrom2MeshEqual(const MEDCouplingPointSet *other, int cellId, double prec) const
1205 {
1206   if(getTypeOfCell(cellId)!=other->getTypeOfCell(cellId))
1207     return false;
1208   std::vector<int> c1,c2;
1209   getNodeIdsOfCell(cellId,c1);
1210   other->getNodeIdsOfCell(cellId,c2);
1211   std::size_t sz=c1.size();
1212   if(sz!=c2.size())
1213     return false;
1214   for(std::size_t i=0;i<sz;i++)
1215     {
1216       std::vector<double> n1,n2;
1217       getCoordinatesOfNode(c1[0],n1);
1218       other->getCoordinatesOfNode(c2[0],n2);
1219       std::transform(n1.begin(),n1.end(),n2.begin(),n1.begin(),std::minus<double>());
1220       std::transform(n1.begin(),n1.end(),n1.begin(),std::ptr_fun<double,double>(fabs));
1221       if(*std::max_element(n1.begin(),n1.end())>prec)
1222         return false;
1223     }
1224   return true;
1225 }
1226
1227 /*!
1228  * Substitutes node coordinates array of \a this mesh with that of \a other mesh
1229  * (i.e. \a this->_coords with \a other._coords) provided that coordinates of the two
1230  * meshes match with a specified precision, else an exception is thrown and \a this
1231  * remains unchanged. In case of success the nodal connectivity of \a this mesh
1232  * is permuted according to new order of nodes.
1233  * Contrary to tryToShareSameCoords() this method makes a deeper analysis of
1234  * coordinates (and so more expensive) than simple equality.
1235  *  \param [in] other - the other mesh whose node coordinates array will be used by
1236  *         \a this mesh in case of their equality.
1237  *  \param [in] epsilon - the precision used to compare coordinates (using infinite norm).
1238  *  \throw If the coordinates array of \a this is not set.
1239  *  \throw If the coordinates array of \a other is not set.
1240  *  \throw If the coordinates of \a this and \a other do not match.
1241  */
1242 void MEDCouplingPointSet::tryToShareSameCoordsPermute(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
1243 {
1244   const DataArrayDouble *coords=other.getCoords();
1245   if(!coords)
1246     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute : No coords specified in other !");
1247   if(!_coords)
1248     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute : No coords specified in this whereas there is any in other !");
1249   int otherNbOfNodes=other.getNumberOfNodes();
1250   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=MergeNodesArray(&other,this);
1251   _coords->incrRef();
1252   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> oldCoords=_coords;
1253   setCoords(newCoords);
1254   bool areNodesMerged;
1255   int newNbOfNodes;
1256   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=buildPermArrayForMergeNode(epsilon,otherNbOfNodes,areNodesMerged,newNbOfNodes);
1257   if(!areNodesMerged)
1258     {
1259       setCoords(oldCoords);
1260       throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute fails : no nodes are mergeable with specified given epsilon !");
1261     }
1262   int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+otherNbOfNodes);
1263   const int *pt=std::find_if(da->getConstPointer()+otherNbOfNodes,da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
1264   if(pt!=da->getConstPointer()+da->getNbOfElems())
1265     {
1266       setCoords(oldCoords);
1267       throw INTERP_KERNEL::Exception("MEDCouplingPointSet::tryToShareSameCoordsPermute fails : some nodes in this are not in other !");
1268     }
1269   setCoords(oldCoords);
1270   renumberNodesInConn(da->getConstPointer()+otherNbOfNodes);
1271   setCoords(coords);
1272 }
1273
1274 MEDCouplingPointSet *MEDCouplingPointSet::buildPartOfMySelf(const int *begin, const int *end, bool keepCoords) const
1275 {
1276   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> ret=buildPartOfMySelfKeepCoords(begin,end);
1277   if(!keepCoords)
1278     ret->zipCoords();
1279   return ret.retn();
1280 }
1281
1282 MEDCouplingPointSet *MEDCouplingPointSet::buildPartOfMySelf2(int start, int end, int step, bool keepCoords) const throw(INTERP_KERNEL::Exception)
1283 {
1284   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> ret=buildPartOfMySelfKeepCoords2(start,end,step);
1285   if(!keepCoords)
1286     ret->zipCoords();
1287   return ret.retn();
1288 }
1289
1290 /*!
1291  Creates a new MEDCouplingUMesh containing some cells of \a this mesh. The cells to
1292  copy are selected basing on specified node ids and the value of \a fullyIn
1293  parameter. If \a fullyIn ==\c true, a cell is copied if its all nodes are in the
1294  array \a begin of node ids. If \a fullyIn ==\c false, a cell is copied if any its
1295  node is in the array of node ids. The created mesh shares the node coordinates array
1296  with \a this mesh.
1297  *  \param [in] begin - the array of node ids.
1298  *  \param [in] end - a pointer to the (last+1)-th element of \a begin.
1299  *  \param [in] fullyIn - if \c true, then cells whose all nodes are in the
1300  *         array \a begin are copied, else cells whose any node is in the
1301  *         array \a begin are copied.
1302  *  \return MEDCouplingPointSet * - new instance of MEDCouplingUMesh. The caller is
1303  *         to delete this mesh using decrRef() as it is no more needed.
1304  *  \throw If the coordinates array is not set.
1305  *  \throw If the nodal connectivity of cells is not defined.
1306  *  \throw If any node id in \a begin is not valid.
1307  *
1308  *  \ref cpp_mcumesh_buildPartOfMySelfNode "Here is a C++ example".<br>
1309  *  \ref  py_mcumesh_buildPartOfMySelfNode "Here is a Python example".
1310  */
1311 MEDCouplingPointSet *MEDCouplingPointSet::buildPartOfMySelfNode(const int *begin, const int *end, bool fullyIn) const
1312 {
1313   DataArrayInt *cellIdsKept=0;
1314   fillCellIdsToKeepFromNodeIds(begin,end,fullyIn,cellIdsKept);
1315   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIdsKept2(cellIdsKept);
1316   return buildPartOfMySelf(cellIdsKept->begin(),cellIdsKept->end(),true);
1317 }
1318
1319 /*!
1320  * Removes duplicates of cells from \a this mesh and returns an array mapping between
1321  * new and old cell ids in "Old to New" mode. Nothing is changed in \a this mesh if no
1322  * equal cells found.
1323  *  \warning Cells of the result mesh are \b not sorted by geometric type, hence,
1324  *           to write this mesh to the MED file, its cells must be sorted using
1325  *           sortCellsInMEDFileFrmt().
1326  *  \param [in] compType - specifies a cell comparison technique. Meaning of its
1327  *          valid values [0,1,2] is as follows.
1328  *   - 0 : "exact". Two cells are considered equal \c iff they have exactly same nodal
1329  *         connectivity and type. This is the strongest policy.
1330  *   - 1 : "permuted same orientation". Two cells are considered equal \c iff they
1331  *         are based on same nodes and have the same type and orientation.
1332  *   - 2 : "nodal". Two cells are considered equal \c iff they
1333  *         are based on same nodes and have the same type. This is the weakest
1334  *         policy, it can be used by users not sensitive to cell orientation.
1335  *  \param [in] startCellId - specifies the cell id at which search for equal cells
1336  *         starts. By default it is 0, which means that all cells in \a this will be
1337  *         scanned.
1338  *  \return DataArrayInt - a new instance of DataArrayInt, of length \a
1339  *           this->getNumberOfCells() before call of this method. The caller is to
1340  *           delete this array using decrRef() as it is no more needed.
1341  *  \throw If the coordinates array is not set.
1342  *  \throw If the nodal connectivity of cells is not defined.
1343  *  \throw If the nodal connectivity includes an invalid id.
1344  *
1345  *  \ref cpp_mcumesh_zipConnectivityTraducer "Here is a C++ example".<br>
1346  *  \ref  py_mcumesh_zipConnectivityTraducer "Here is a Python example".
1347  */
1348 DataArrayInt *MEDCouplingPointSet::zipConnectivityTraducer(int compType, int startCellId) throw(INTERP_KERNEL::Exception)
1349 {
1350   DataArrayInt *commonCells=0,*commonCellsI=0;
1351   findCommonCells(compType,startCellId,commonCells,commonCellsI);
1352   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> commonCellsTmp(commonCells),commonCellsITmp(commonCellsI);
1353   int newNbOfCells=-1;
1354   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::BuildOld2NewArrayFromSurjectiveFormat2(getNumberOfCells(),commonCells->begin(),commonCellsI->begin(),
1355                                                                                                           commonCellsI->end(),newNbOfCells);
1356   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2=ret->invertArrayO2N2N2O(newNbOfCells);
1357   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> self=buildPartOfMySelf(ret2->begin(),ret2->end(),true);
1358   shallowCopyConnectivityFrom(self);
1359   return ret.retn();
1360 }
1361
1362 /*!
1363  * Checks if \a this and \a other meshes are geometrically equivalent, else an
1364  * exception is thrown. The meshes are
1365  * considered equivalent if (1) \a this mesh contains the same nodes as the \a other
1366  * mesh (with a specified precision) and (2) \a this mesh contains the same cells as
1367  * the \a other mesh (with use of a specified cell comparison technique). The mapping
1368  * from \a other to \a this for nodes and cells is returned via out parameters.
1369  *  \param [in] other - the mesh to compare with.
1370  *  \param [in] cellCompPol - id [0-2] of cell comparison method. See meaning of
1371  *         each method in description of MEDCouplingPointSet::zipConnectivityTraducer().
1372  *  \param [in] prec - the precision used to compare nodes of the two meshes.
1373  *  \param [out] cellCor - a cell permutation array in "Old to New" mode. The caller is
1374  *         to delete this array using decrRef() as it is no more needed.
1375  *  \param [out] nodeCor - a node permutation array in "Old to New" mode. The caller is
1376  *         to delete this array using decrRef() as it is no more needed.
1377  *  \throw If the two meshes do not match.
1378  *
1379  *  \ref cpp_mcumesh_checkDeepEquivalWith "Here is a C++ example".<br>
1380  *  \ref  py_mcumesh_checkDeepEquivalWith "Here is a Python example".
1381  */
1382 void MEDCouplingPointSet::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
1383                                                DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
1384 {
1385   if(!other)
1386     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalWith : input is null !");
1387   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
1388   if(!otherC)
1389     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalWith : other is not a PointSet mesh !");
1390   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> m=dynamic_cast<MEDCouplingPointSet *>(mergeMyselfWith(otherC));
1391   bool areNodesMerged;
1392   int newNbOfNodes;
1393   int oldNbOfNodes=getNumberOfNodes();
1394   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=m->buildPermArrayForMergeNode(prec,oldNbOfNodes,areNodesMerged,newNbOfNodes);
1395   //mergeNodes
1396   if(!areNodesMerged)
1397     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Nodes are incompatible ! ");
1398   const int *pt=std::find_if(da->getConstPointer()+oldNbOfNodes,da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),oldNbOfNodes-1));
1399   if(pt!=da->getConstPointer()+da->getNbOfElems())
1400     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some nodes in other are not in this !");
1401   m->renumberNodes(da->getConstPointer(),newNbOfNodes);
1402   //
1403   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nodeCor2=da->substr(oldNbOfNodes);
1404   da=m->mergeNodes(prec,areNodesMerged,newNbOfNodes);
1405   //
1406   da=m->zipConnectivityTraducer(cellCompPol);
1407   int nbCells=getNumberOfCells();
1408   int maxId=-1;
1409   if(nbCells!=0)
1410     maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+nbCells);
1411   pt=std::find_if(da->getConstPointer()+nbCells,da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
1412   if(pt!=da->getConstPointer()+da->getNbOfElems())
1413     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some cells in other are not in this !");
1414   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellCor2=da->selectByTupleId2(nbCells,da->getNbOfElems(),1);
1415   nodeCor=nodeCor2->isIdentity()?0:nodeCor2.retn();
1416   cellCor=cellCor2->isIdentity()?0:cellCor2.retn();
1417 }
1418
1419 /*!
1420  * Checks if \a this and \a other meshes are geometrically equivalent, else an
1421  * exception is thrown. The meshes are considered equivalent if (1) they share one
1422  * node coordinates array and (2) they contain the same cells (with use of a specified
1423  * cell comparison technique). The mapping from cells of the \a other to ones of \a this
1424  * is returned via an out parameter.
1425  *  \param [in] other - the mesh to compare with.
1426  *  \param [in] cellCompPol - id [0-2] of cell comparison method. See the meaning of
1427  *         each method in description of MEDCouplingPointSet::zipConnectivityTraducer().
1428  *  \param [in] prec - a not used parameter.
1429  *  \param [out] cellCor - the permutation array in "Old to New" mode. The caller is
1430  *         to delete this array using decrRef() as it is no more needed.
1431  *  \throw If the two meshes do not match.
1432  *
1433  * \ref cpp_mcumesh_checkDeepEquivalWith "Here is a C++ example".<br>
1434  * \ref  py_mcumesh_checkDeepEquivalWith "Here is a Python example".
1435  */
1436 void MEDCouplingPointSet::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
1437                                                        DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception)
1438 {
1439   if(!other)
1440     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalOnSameNodesWith : input is null !");
1441   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
1442   if(!otherC)
1443     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkDeepEquivalOnSameNodesWith : other is not a PointSet mesh !");
1444   if(_coords!=otherC->_coords)
1445     throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : meshes do not share the same coordinates ! Use tryToShareSameCoordinates or call checkDeepEquivalWith !");
1446   MEDCouplingAutoRefCountObjectPtr<MEDCouplingPointSet> m=mergeMyselfWithOnSameCoords(otherC);
1447   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=m->zipConnectivityTraducer(cellCompPol);
1448   int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+getNumberOfCells());
1449   const int *pt=std::find_if(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
1450   if(pt!=da->getConstPointer()+da->getNbOfElems())
1451     {
1452       throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : some cells in other are not in this !");
1453     }
1454   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellCor2=da->selectByTupleId2(getNumberOfCells(),da->getNbOfElems(),1);
1455   cellCor=cellCor2->isIdentity()?0:cellCor2.retn();
1456 }
1457
1458 void MEDCouplingPointSet::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
1459 {
1460   MEDCouplingMesh::checkFastEquivalWith(other,prec);
1461   //other not null checked by the line before
1462   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
1463   if(!otherC)
1464     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkFastEquivalWith : fails because other is not a pointset mesh !");
1465   int nbOfCells=getNumberOfCells();
1466   if(nbOfCells<1)
1467     return ;
1468   bool status=true;
1469   status&=areCellsFrom2MeshEqual(otherC,0,prec);
1470   status&=areCellsFrom2MeshEqual(otherC,nbOfCells/2,prec);
1471   status&=areCellsFrom2MeshEqual(otherC,nbOfCells-1,prec);
1472   if(!status)
1473     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Two meshes are not equal because on 3 test cells some difference have been detected !");
1474 }
1475
1476 /*!
1477  * Finds cells whose all or some nodes are in a given array of node ids.
1478  *  \param [in] begin - the array of node ids.
1479  *  \param [in] end - a pointer to the (last+1)-th element of \a begin.
1480  *  \param [in] fullyIn - if \c true, then cells whose all nodes are in the
1481  *         array \a begin are returned only, else cells whose any node is in the
1482  *         array \a begin are returned.
1483  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of found
1484  *         cells. The caller is to delete this array using decrRef() as it is no more
1485  *         needed.
1486  *  \throw If the coordinates array is not set.
1487  *  \throw If the nodal connectivity of cells is not defined.
1488  *  \throw If any cell id in \a begin is not valid.
1489  *
1490  * \sa MEDCouplingPointSet::getCellIdsFullyIncludedInNodeIds
1491  *
1492  *  \ref cpp_mcumesh_getCellIdsLyingOnNodes "Here is a C++ example".<br>
1493  *  \ref  py_mcumesh_getCellIdsLyingOnNodes "Here is a Python example".
1494  */
1495 DataArrayInt *MEDCouplingPointSet::getCellIdsLyingOnNodes(const int *begin, const int *end, bool fullyIn) const
1496 {
1497   DataArrayInt *cellIdsKept=0;
1498   fillCellIdsToKeepFromNodeIds(begin,end,fullyIn,cellIdsKept);
1499   cellIdsKept->setName(getName());
1500   return cellIdsKept;
1501 }
1502
1503 /*!
1504  * Finds cells whose all nodes are in a given array of node ids.
1505  * This method is a specialization of MEDCouplingPointSet::getCellIdsLyingOnNodes (true
1506  * as last input argument).
1507  *  \param [in] partBg - the array of node ids.
1508  *  \param [in] partEnd - a pointer to a (last+1)-th element of \a partBg.
1509  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of found
1510  *          cells. The caller is to delete this array using decrRef() as it is no
1511  *          more needed.
1512  *  \throw If the coordinates array is not set.
1513  *  \throw If the nodal connectivity of cells is not defined.
1514  *  \throw If any cell id in \a partBg is not valid.
1515  *
1516  * \sa MEDCouplingPointSet::getCellIdsLyingOnNodes
1517  *
1518  *  \ref cpp_mcumesh_getCellIdsFullyIncludedInNodeIds "Here is a C++ example".<br>
1519  *  \ref  py_mcumesh_getCellIdsFullyIncludedInNodeIds "Here is a Python example".
1520  */
1521 DataArrayInt *MEDCouplingPointSet::getCellIdsFullyIncludedInNodeIds(const int *partBg, const int *partEnd) const
1522 {
1523   return getCellIdsLyingOnNodes(partBg,partEnd,true);
1524 }
1525
1526 /*!
1527  * Removes unused nodes (the node coordinates array is shorten) and returns an array
1528  * mapping between new and old node ids in "Old to New" mode. -1 values in the returned
1529  * array mean that the corresponding old node is no more used.
1530  *  \return DataArrayInt * - a new instance of DataArrayInt of length \a
1531  *           this->getNumberOfNodes() before call of this method. The caller is to
1532  *           delete this array using decrRef() as it is no more needed.
1533  *  \throw If the coordinates array is not set.
1534  *  \throw If the nodal connectivity of cells is not defined.
1535  *  \throw If the nodal connectivity includes an invalid id.
1536  *
1537  *  \ref cpp_mcumesh_zipCoordsTraducer "Here is a C++ example".<br>
1538  *  \ref  py_mcumesh_zipCoordsTraducer "Here is a Python example".
1539  */
1540 DataArrayInt *MEDCouplingPointSet::zipCoordsTraducer() throw(INTERP_KERNEL::Exception)
1541 {
1542   int newNbOfNodes=-1;
1543   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> traducer=getNodeIdsInUse(newNbOfNodes);
1544   renumberNodes(traducer->getConstPointer(),newNbOfNodes);
1545   return traducer.retn();
1546 }
1547
1548 /*!
1549  * Merges nodes equal within \a precision and returns an array describing the
1550  * permutation used to remove duplicate nodes.
1551  *  \param [in] precision - minimal absolute distance between two nodes at which they are
1552  *              considered not coincident.
1553  *  \param [out] areNodesMerged - is set to \c true if any coincident nodes removed.
1554  *  \param [out] newNbOfNodes - number of nodes remaining after the removal.
1555  *  \return DataArrayInt * - the permutation array in "Old to New" mode. For more
1556  *          info on "Old to New" mode see \ref MEDCouplingArrayRenumbering. The caller
1557  *          is to delete this array using decrRef() as it is no more needed.
1558  *  \throw If the coordinates array is not set.
1559  *  \throw If the nodal connectivity of cells is not defined.
1560  *
1561  *  \ref cpp_mcumesh_mergeNodes "Here is a C++ example".<br>
1562  *  \ref  py_mcumesh_mergeNodes "Here is a Python example".
1563  */
1564 DataArrayInt *MEDCouplingPointSet::mergeNodes(double precision, bool& areNodesMerged, int& newNbOfNodes)
1565 {
1566   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=buildPermArrayForMergeNode(precision,-1,areNodesMerged,newNbOfNodes);
1567   if(areNodesMerged)
1568     renumberNodes(ret->begin(),newNbOfNodes);
1569   return ret.retn();
1570 }
1571
1572 /*!
1573  * Merges nodes equal within \a precision and returns an array describing the
1574  * permutation used to remove duplicate nodes. In contrast to mergeNodes(), location
1575  *  of merged nodes is changed to be at their barycenter.
1576  *  \param [in] precision - minimal absolute distance between two nodes at which they are
1577  *              considered not coincident.
1578  *  \param [out] areNodesMerged - is set to \c true if any coincident nodes removed.
1579  *  \param [out] newNbOfNodes - number of nodes remaining after the removal.
1580  *  \return DataArrayInt * - the permutation array in "Old to New" mode. For more
1581  *          info on "Old to New" mode see \ref MEDCouplingArrayRenumbering. The caller
1582  *          is to delete this array using decrRef() as it is no more needed.
1583  *  \throw If the coordinates array is not set.
1584  *  \throw If the nodal connectivity of cells is not defined.
1585  *
1586  *  \ref cpp_mcumesh_mergeNodes "Here is a C++ example".<br>
1587  *  \ref  py_mcumesh_mergeNodes "Here is a Python example".
1588  */
1589 DataArrayInt *MEDCouplingPointSet::mergeNodes2(double precision, bool& areNodesMerged, int& newNbOfNodes)
1590 {
1591   DataArrayInt *ret=buildPermArrayForMergeNode(precision,-1,areNodesMerged,newNbOfNodes);
1592   if(areNodesMerged)
1593     renumberNodes2(ret->getConstPointer(),newNbOfNodes);
1594   return ret;
1595 }