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