src/MEDCoupling/MEDCouplingPointSet.cxx

   1 // Copyright (C) 2007-2012  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. Nodes with id strictly lower than \a
 248  *              limitTupleId are \b not considered. Put -1 to this parameter to have
 249  *              all nodes treated.
 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   DataArrayInt *ret=buildNewNumberingFromCommonNodesFormat(comm,commI,newNbOfNodes);
 263   areNodesMerged=(oldNbOfNodes!=newNbOfNodes);
 264   comm->decrRef();
 265   commI->decrRef();
 266   return ret;
 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 between two nodes at which they are
 274  *              considered not coincident.
 275  *  \param [in] limitNodeId - limit node id. Nodes with id strictly lower than \a
 276  *              limitTupleId are \b not considered. Put -1 to this parameter to have
 277  *              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 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 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  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
 374  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
 375  * This value is asked because often known by the caller of this method.
 376  * @param newNodeNumbers array specifying the new numbering in old2New convention..
 377  * @param newNbOfNodes the new number of nodes.
 378  */
 379 void MEDCouplingPointSet::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
 380 {
 381   if(!_coords)
 382     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::renumberNodes : no coords specified !");
 383   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->renumberAndReduce(newNodeNumbers,newNbOfNodes);
 384   setCoords(newCoords);
 385 }
 386
 387 /*
 388  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
 389  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
 390  * This value is asked because often known by the caller of this method.
 391  * Contrary to ParaMEDMEM::MEDCouplingPointSet::renumberNodes method for merged nodes the barycenter of them is computed here.
 392  *
 393  * @param newNodeNumbers array specifying the new numbering.
 394  * @param newNbOfNodes the new number of nodes.
 395  */
 396 void MEDCouplingPointSet::renumberNodes2(const int *newNodeNumbers, int newNbOfNodes)
 397 {
 398   DataArrayDouble *newCoords=DataArrayDouble::New();
 399   std::vector<int> div(newNbOfNodes);
 400   int spaceDim=getSpaceDimension();
 401   newCoords->alloc(newNbOfNodes,spaceDim);
 402   newCoords->copyStringInfoFrom(*_coords);
 403   newCoords->fillWithZero();
 404   int oldNbOfNodes=getNumberOfNodes();
 405   double *ptToFill=newCoords->getPointer();
 406   const double *oldCoordsPtr=_coords->getConstPointer();
 407   for(int i=0;i<oldNbOfNodes;i++)
 408     {
 409       std::transform(oldCoordsPtr+i*spaceDim,oldCoordsPtr+(i+1)*spaceDim,ptToFill+newNodeNumbers[i]*spaceDim,
 410                      ptToFill+newNodeNumbers[i]*spaceDim,std::plus<double>());
 411       div[newNodeNumbers[i]]++;
 412     }
 413   for(int i=0;i<newNbOfNodes;i++)
 414     ptToFill=std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)div[i]));
 415   setCoords(newCoords);
 416   newCoords->decrRef();
 417 }
 418
 419 /*!
 420  * Computes the minimum box bounding all nodes. The edges of the box are parallel to
 421  * the Cartesian coordinate axes. The bounding box is described by coordinates of its
 422  * two extremum points with minimal and maximal coordinates.
 423  *  \param [out] bbox - array filled with coordinates of extremum points in "no
 424  *         interlace" mode, i.e. xMin, xMax, yMin, yMax, zMin, zMax (if in 3D). This
 425  *         array, of length 2 * \a this->getSpaceDimension() at least, is to be
 426  *         pre-allocated by the caller.
 427  *  \throw If the coordinates array is not set.
 428  *
 429  *  \ref cpp_mcpointset_getBoundingBox "Here is a C++ example".<br>
 430  *  \ref  py_mcpointset_getBoundingBox "Here is a Python example".
 431  */
 432 void MEDCouplingPointSet::getBoundingBox(double *bbox) const throw(INTERP_KERNEL::Exception)
 433 {
 434   if(!_coords)
 435     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getBoundingBox : Coordinates not set !");
 436   _coords->getMinMaxPerComponent(bbox);
 437 }
 438
 439 /*!
 440  * Removes "free" nodes, i.e. nodes not used to define any element.
 441  *  \throw If the coordinates array is not set.
 442  *  \throw If the elements are not defined.
 443  */
 444 void MEDCouplingPointSet::zipCoords()
 445 {
 446   checkFullyDefined();
 447   DataArrayInt *traducer=zipCoordsTraducer();
 448   traducer->decrRef();
 449 }
 450
 451 struct MEDCouplingCompAbs
 452 {
 453   bool operator()(double x, double y) { return std::abs(x)<std::abs(y);}
 454 };
 455
 456 /*!
 457  * Returns the carateristic dimension of \a this point set, that is a maximal
 458  * absolute values of node coordinates.
 459  *  \throw If the coordinates array is not set.
 460  */
 461 double MEDCouplingPointSet::getCaracteristicDimension() const
 462 {
 463   if(!_coords)
 464     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCaracteristicDimension : Coordinates not set !");
 465   const double *coords=_coords->getConstPointer();
 466   int nbOfValues=_coords->getNbOfElems();
 467   return std::abs(*std::max_element(coords,coords+nbOfValues,MEDCouplingCompAbs()));
 468 }
 469
 470 /*!
 471  * 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
 472  * around origin of 'radius' 1.
 473  *
 474  * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
 475  *
 476  * \warning this method is non const and alterates coordinates in \b this without modifying.
 477  */
 478 void MEDCouplingPointSet::recenterForMaxPrecision(double eps) throw(INTERP_KERNEL::Exception)
 479 {
 480   if(!_coords)
 481     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::recenterForMaxPrecision : Coordinates not set !");
 482   _coords->recenterForMaxPrecision(eps);
 483   updateTime();
 484 }
 485
 486 /*!
 487  * Rotates \a this set of nodes by \a angle around either an axis (in 3D) or a point
 488  * (in 2D).
 489  *  \param [in] center - coordinates either of an origin of rotation axis (in 3D) or
 490  *         of center of rotation (in 2D). This array is to be of size \a
 491  *         this->getSpaceDimension() at least.
 492  *  \param [in] vector - 3 components of a vector defining direction of the rotation
 493  *         axis in 3D. In 2D this parameter is not used.
 494  *  \param [in] angle - the rotation angle in radians.
 495  *  \throw If the coordinates array is not set.
 496  *  \throw If \a this->getSpaceDimension() != 2 && \a this->getSpaceDimension() != 3.
 497  *  \throw If \a center == NULL
 498  *  \throw If \a vector == NULL && \a this->getSpaceDimension() == 3.
 499  *  \throw If Magnitude of \a vector is zero.
 500  *
 501  *  \ref cpp_mcpointset_rotate "Here is a C++ example".<br>
 502  *  \ref  py_mcpointset_rotate "Here is a Python example".
 503  */
 504 void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
 505 {
 506   int spaceDim=getSpaceDimension();
 507   if(spaceDim==3)
 508     rotate3D(center,vector,angle);
 509   else if(spaceDim==2)
 510     rotate2D(center,angle);
 511   else
 512     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
 513   _coords->declareAsNew();
 514   updateTime();
 515 }
 516
 517 /*!
 518  * Translates \a this set of nodes.
 519  *  \param [in] vector - components of a translation vector. This array is to be of
 520  *         size \a this->getSpaceDimension() at least.
 521  *  \throw If the coordinates array is not set.
 522  *  \throw If \a vector == NULL.
 523  *
 524  *  \ref cpp_mcpointset_translate "Here is a C++ example".<br>
 525  *  \ref  py_mcpointset_translate "Here is a Python example".
 526  */
 527 void MEDCouplingPointSet::translate(const double *vector)
 528 {
 529   if(!vector)
 530     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::translate : NULL input vector !");
 531   if(!_coords)
 532     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::translate : no coordinates set !");
 533   double *coords=_coords->getPointer();
 534   int nbNodes=getNumberOfNodes();
 535   int dim=getSpaceDimension();
 536   for(int i=0; i<nbNodes; i++)
 537     for(int idim=0; idim<dim;idim++)
 538       coords[i*dim+idim]+=vector[idim];
 539   _coords->declareAsNew();
 540   updateTime();
 541 }
 542
 543
 544 /*!
 545  * Applies scaling transformation to \a this set of nodes.
 546  *  \param [in] point - coordinates of a scaling center. This array is to be of
 547  *         size \a this->getSpaceDimension() at least.
 548  *  \param [in] factor - a scale factor.
 549  *  \throw If the coordinates array is not set.
 550  *  \throw If \a point == NULL.
 551  *
 552  *  \ref cpp_mcpointset_scale "Here is a C++ example".<br>
 553  *  \ref  py_mcpointset_scale "Here is a Python example".
 554  */
 555 void MEDCouplingPointSet::scale(const double *point, double factor)
 556 {
 557   if(!point)
 558     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::scale : NULL input point !");
 559   if(!_coords)
 560     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::scale : no coordinates set !");
 561   double *coords=_coords->getPointer();
 562   int nbNodes=getNumberOfNodes();
 563   int dim=getSpaceDimension();
 564   for(int i=0;i<nbNodes;i++)
 565     {
 566       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
 567       std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
 568       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
 569     }
 570   _coords->declareAsNew();
 571   updateTime();
 572 }
 573
 574 /*!
 575  * Converts \a this set of points to an other dimension by changing number of
 576  * components of point coordinates. If the dimension increases, added components
 577  * are filled with \a dftValue. If the dimension decreases, last components are lost.
 578  * If the new dimension is same as \a this->getSpaceDimension(), nothing is done.
 579  *  \param [in] newSpaceDim - the new space dimension.
 580  *  \param [in] dftValue - the value to assign to added components of point coordinates
 581  *         (if the dimension increases).
 582  *  \throw If the coordinates array is not set.
 583  *  \throw If \a newSpaceDim < 1.
 584  */
 585 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue) throw(INTERP_KERNEL::Exception)
 586 {
 587   if(getCoords()==0)
 588     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
 589   if(newSpaceDim<1)
 590     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
 591   int oldSpaceDim=getSpaceDimension();
 592   if(newSpaceDim==oldSpaceDim)
 593     return ;
 594   DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
 595   setCoords(newCoords);
 596   newCoords->decrRef();
 597   updateTime();
 598 }
 599
 600 /*!
 601  * Substitutes \a this->_coords with \a other._coords provided that coordinates of
 602  * the two point sets match with a specified precision, else an exception is thrown.
 603  *  \param [in] other - the other point set whose coordinates array will be used by
 604  *         \a this point set in case of their equality.
 605  *  \param [in] epsilon - the precision used to compare coordinates.
 606  *  \throw If the coordinates array of \a this is not set.
 607  *  \throw If the coordinates array of \a other is not set.
 608  *  \throw If the coordinates of \a this and \a other do not match.
 609  */
 610 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
 611 {
 612   if(_coords==other._coords)
 613     return ;
 614   if(!_coords)
 615     throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
 616   if(!other._coords)
 617     throw INTERP_KERNEL::Exception("Other instance has no coords whereas current has !");
 618   if(!_coords->isEqualWithoutConsideringStr(*other._coords,epsilon))
 619     throw INTERP_KERNEL::Exception("Coords are not the same !");
 620   setCoords(other._coords);
 621 }
 622
 623 /*!
 624  * This method duplicates the nodes whose ids are in [\b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd) and put the result of their duplication at the end
 625  * of existing node ids.
 626  *
 627  * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
 628  * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
 629  */
 630 void MEDCouplingPointSet::duplicateNodesInCoords(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd) throw(INTERP_KERNEL::Exception)
 631 {
 632   if(!_coords)
 633     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::duplicateNodesInCoords : no coords set !");
 634   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->selectByTupleIdSafe(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd);
 635   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords2=DataArrayDouble::Aggregate(_coords,newCoords);
 636   setCoords(newCoords2);
 637 }
 638
 639 /*!
 640  * Finds nodes located at distance lower that \a eps from a specified plane.
 641  *  \param [in] pt - 3 components of a point defining location of the plane.
 642  *  \param [in] vec - 3 components of a normal vector to the plane. Vector magnitude
 643  *         must be greater than 10*\a eps.
 644  *  \param [in] eps - maximal distance of a node from the plane at which the node is
 645  *         considered to lie on the plane.
 646  *  \param [in,out] nodes - a vector returning ids of found nodes. This vector is not
 647  *         cleared before filling in.
 648  *  \throw If the coordinates array is not set.
 649  *  \throw If \a pt == NULL.
 650  *  \throw If \a vec == NULL.
 651  *  \throw If the magnitude of \a vec is zero.
 652  *  \throw If \a this->getSpaceDimension() != 3.
 653  */
 654 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
 655 {
 656   if(getSpaceDimension()!=3)
 657     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
 658   if(!pt)
 659     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : NULL point pointer specified !");
 660   if(!vec)
 661     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : NULL vector pointer specified !");
 662   int nbOfNodes=getNumberOfNodes();
 663   double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
 664   double deno=sqrt(a*a+b*b+c*c);
 665   if(deno<std::numeric_limits<double>::min())
 666     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : vector pointer specified has norm equal to 0. !");
 667   const double *work=_coords->getConstPointer();
 668   for(int i=0;i<nbOfNodes;i++)
 669     {
 670       if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
 671         nodes.push_back(i);
 672       work+=3;
 673     }
 674 }
 675
 676 /*!
 677  * Finds nodes located at distance lower that \a eps from a specified line in 2D and 3D.
 678  *  \param [in] pt - components of coordinates of an initial point of the line. This
 679  *         array is to be of size \a this->getSpaceDimension() at least.
 680  *  \param [in] vec - components of a vector defining the line direction. This array
 681  *         is to be of size \a this->getSpaceDimension() at least. Vector magnitude
 682  *         must be greater than 10*\a eps.
 683  *  \param [in] eps - maximal distance of a node from the line at which the node is
 684  *         considered to lie on the line.
 685  *  \param [in,out] nodes - a vector returning ids of found nodes. This vector is not
 686  *         cleared before filling in.
 687  *  \throw If the coordinates array is not set.
 688  *  \throw If \a pt == NULL.
 689  *  \throw If \a vec == NULL.
 690  *  \throw If the magnitude of \a vec is zero.
 691  *  \throw If ( \a this->getSpaceDimension() != 3 && \a this->getSpaceDimension() != 2 ).
 692  */
 693 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
 694 {
 695   int spaceDim=getSpaceDimension();
 696   if(spaceDim!=2 && spaceDim!=3)
 697     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
 698   if(!pt)
 699     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : NULL point pointer specified !");
 700   if(!vec)
 701     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : NULL vector pointer specified !");
 702   int nbOfNodes=getNumberOfNodes();
 703   double den=0.;
 704   for(int i=0;i<spaceDim;i++)
 705     den+=vec[i]*vec[i];
 706   double deno=sqrt(den);
 707   if(deno<10.*eps)
 708     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
 709   INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
 710   for(int i=0;i<spaceDim;i++)
 711     vecn[i]=vec[i]/deno;
 712   const double *work=_coords->getConstPointer();
 713   if(spaceDim==2)
 714     {
 715       for(int i=0;i<nbOfNodes;i++)
 716         {
 717           if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
 718             nodes.push_back(i);
 719           work+=2;
 720         }
 721     }
 722   else
 723     {
 724       for(int i=0;i<nbOfNodes;i++)
 725         {
 726           double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
 727           double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
 728           double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
 729           if(std::sqrt(a*a+b*b+c*c)<eps)
 730             nodes.push_back(i);
 731           work+=3;
 732         }
 733     }
 734 }
 735
 736 /*!
 737  * Returns a new array of node coordinates by concatenating node coordinates of two
 738  * given point sets, so that (1) the number of nodes in the result array is a sum of the
 739  * number of nodes of given point sets and (2) the number of component in the result array
 740  * is same as that of each of given point sets. Info on components is copied from the first
 741  * of the given point set. Space dimension of the given point sets must be the same.
 742  *  \param [in] m1 - a point set whose coordinates will be included in the result array.
 743  *  \param [in] m2 - another point set whose coordinates will be included in the
 744  *         result array.
 745  *  \return DataArrayDouble * - the new instance of DataArrayDouble.
 746  *          The caller is to delete this result array using decrRef() as it is no more
 747  *          needed.
 748  *  \throw If both \a m1 and \a m2 are NULL.
 749  *  \throw If \a m1->getSpaceDimension() != \a m2->getSpaceDimension().
 750  */
 751 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception)
 752 {
 753   int spaceDim=m1->getSpaceDimension();
 754   if(spaceDim!=m2->getSpaceDimension())
 755     throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
 756   return DataArrayDouble::Aggregate(m1->getCoords(),m2->getCoords());
 757 }
 758
 759 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception)
 760 {
 761   if(ms.empty())
 762     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : input array must be NON EMPTY !");
 763   std::vector<const MEDCouplingPointSet *>::const_iterator it=ms.begin();
 764   std::vector<const DataArrayDouble *> coo(ms.size());
 765   int spaceDim=(*it)->getSpaceDimension();
 766   coo[0]=(*it++)->getCoords();
 767   for(int i=1;it!=ms.end();it++,i++)
 768     {
 769       const DataArrayDouble *tmp=(*it)->getCoords();
 770       if(tmp)
 771         {
 772           if((*it)->getSpaceDimension()==spaceDim)
 773             coo[i]=tmp;
 774           else
 775             throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
 776         }
 777       else
 778         throw INTERP_KERNEL::Exception("Empty coords detected during call of MergeNodesArray !");
 779     }
 780   return DataArrayDouble::Aggregate(coo);
 781 }
 782
 783 /*!
 784  * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
 785  * This method is used during unserialization process.
 786  */
 787 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
 788 {
 789   switch(type)
 790     {
 791     case UNSTRUCTURED:
 792       return MEDCouplingUMesh::New();
 793     case UNSTRUCTURED_DESC:
 794       return MEDCouplingUMeshDesc::New();
 795     default:
 796       throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
 797     }
 798 }
 799
 800 /*!
 801  * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
 802  */
 803 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 804 {
 805   int it,order;
 806   double time=getTime(it,order);
 807   if(_coords)
 808     {
 809       int spaceDim=getSpaceDimension();
 810       littleStrings.resize(spaceDim+4);
 811       littleStrings[0]=getName();
 812       littleStrings[1]=getDescription();
 813       littleStrings[2]=_coords->getName();
 814       littleStrings[3]=getTimeUnit();
 815       for(int i=0;i<spaceDim;i++)
 816         littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
 817       tinyInfo.clear();
 818       tinyInfo.push_back(getType());
 819       tinyInfo.push_back(spaceDim);
 820       tinyInfo.push_back(getNumberOfNodes());
 821       tinyInfo.push_back(it);
 822       tinyInfo.push_back(order);
 823       tinyInfoD.push_back(time);
 824     }
 825   else
 826     {
 827       littleStrings.resize(3);
 828       littleStrings[0]=getName();
 829       littleStrings[1]=getDescription();
 830       littleStrings[2]=getTimeUnit();
 831       tinyInfo.clear();
 832       tinyInfo.push_back(getType());
 833       tinyInfo.push_back(-1);
 834       tinyInfo.push_back(-1);
 835       tinyInfo.push_back(it);
 836       tinyInfo.push_back(order);
 837       tinyInfoD.push_back(time);
 838     }
 839 }
 840
 841 /*!
 842  * Third and final step of serialization process.
 843  */
 844 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 845 {
 846   if(_coords)
 847     {
 848       a2=const_cast<DataArrayDouble *>(getCoords());
 849       a2->incrRef();
 850     }
 851   else
 852     a2=0;
 853 }
 854
 855 /*!
 856  * Second step of serialization process.
 857  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 858  */
 859 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 860 {
 861   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 862     {
 863       a2->alloc(tinyInfo[2],tinyInfo[1]);
 864       littleStrings.resize(tinyInfo[1]+4);
 865     }
 866   else
 867     {
 868       littleStrings.resize(3);
 869     }
 870 }
 871
 872 /*!
 873  * Second and final unserialization process.
 874  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 875  */
 876 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
 877 {
 878   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 879     {
 880       setCoords(a2);
 881       setName(littleStrings[0].c_str());
 882       setDescription(littleStrings[1].c_str());
 883       a2->setName(littleStrings[2].c_str());
 884       setTimeUnit(littleStrings[3].c_str());
 885       for(int i=0;i<tinyInfo[1];i++)
 886         getCoords()->setInfoOnComponent(i,littleStrings[i+4].c_str());
 887       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 888     }
 889   else
 890     {
 891       setName(littleStrings[0].c_str());
 892       setDescription(littleStrings[1].c_str());
 893       setTimeUnit(littleStrings[2].c_str());
 894       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 895     }
 896 }
 897
 898 void MEDCouplingPointSet::checkCoherency() const throw(INTERP_KERNEL::Exception)
 899 {
 900   if(!_coords)
 901     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkCoherency : no coordinates set !");
 902 }
 903
 904 /*!
 905  * Intersect Bounding Box given 2 Bounding Boxes.
 906  */
 907 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
 908 {
 909   double* bbtemp = new double[2*dim];
 910   double deltamax=0.0;
 911
 912   for (int i=0; i< dim; i++)
 913     {
 914       double delta = bb1[2*i+1]-bb1[2*i];
 915       if ( delta > deltamax )
 916         {
 917           deltamax = delta ;
 918         }
 919     }
 920   for (int i=0; i<dim; i++)
 921     {
 922       bbtemp[i*2]=bb1[i*2]-deltamax*eps;
 923       bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
 924     }
 925
 926   for (int idim=0; idim < dim; idim++)
 927     {
 928       bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
 929         && (bb2[idim*2]<bbtemp[idim*2+1]) ;
 930       if (!intersects)
 931         {
 932           delete [] bbtemp;
 933           return false;
 934         }
 935     }
 936   delete [] bbtemp;
 937   return true;
 938 }
 939
 940 /*!
 941  * Intersect 2 given Bounding Boxes.
 942  */
 943 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
 944 {
 945   double* bbtemp = new double[2*dim];
 946   double deltamax=0.0;
 947
 948   for (int i=0; i< dim; i++)
 949     {
 950       double delta = bb2[2*i+1]-bb2[2*i];
 951       if ( delta > deltamax )
 952         {
 953           deltamax = delta ;
 954         }
 955     }
 956   for (int i=0; i<dim; i++)
 957     {
 958       bbtemp[i*2]=bb2[i*2]-deltamax*eps;
 959       bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
 960     }
 961
 962   bool intersects = !bb1.isDisjointWith( bbtemp );
 963   delete [] bbtemp;
 964   return intersects;
 965 }
 966
 967 /*!
 968  * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
 969  */
 970 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
 971 {
 972   double *coords=_coords->getPointer();
 973   int nbNodes=getNumberOfNodes();
 974   Rotate3DAlg(center,vect,angle,nbNodes,coords);
 975 }
 976
 977 /*!
 978  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
 979  * around an axe ('center','vect') and with angle 'angle'.
 980  */
 981 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
 982 {
 983   if(!center || !vect)
 984     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : null vector in input !");
 985   double sina=sin(angle);
 986   double cosa=cos(angle);
 987   double vectorNorm[3];
 988   double matrix[9];
 989   double matrixTmp[9];
 990   double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
 991   if(norm<std::numeric_limits<double>::min())
 992     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : magnitude of input vector is too close of 0. !");
 993   std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
 994   //rotation matrix computation
 995   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;
 996   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
 997   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
 998   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
 999   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
1000   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
1001   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
1002   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
1003   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
1004   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
1005   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
1006   //rotation matrix computed.
1007   double tmp[3];
1008   for(int i=0; i<nbNodes; i++)
1009     {
1010       std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
1011       coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
1012       coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
1013       coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
1014     }
1015 }
1016
1017 /*!
1018  * Creates a new MEDCouplingMesh containing a part of cells of \a this mesh. The new
1019  * mesh shares a coordinates array with \a this one. The cells to include to the
1020  * result mesh are specified by an array of cell ids.
1021  *  \param [in] start - an array of cell ids to include to the result mesh.
1022  *  \param [in] end - specifies the end of the array \a start, so that
1023  *              the last value of \a start is \a end[ -1 ].
1024  *  \return MEDCouplingMesh * - a new instance of MEDCouplingMesh. The caller is to
1025  *         delete this mesh using decrRef() as it is no more needed.
1026  */
1027 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
1028 {
1029   return buildPartOfMySelf(start,end,true);
1030 }
1031
1032 /*!
1033  * Creates a new MEDCouplingMesh containing a part of cells of \a this mesh. The
1034  * cells to include to the result mesh are specified by an array of cell ids.
1035  * <br> This method additionally returns a renumbering map in "Old to New" mode
1036  * which allows the caller to know the mapping between nodes in \a this and the result mesh.
1037  *  \param [in] start - an array of cell ids to include to the result mesh.
1038  *  \param [in] end - specifies the end of the array \a start, so that
1039  *              the last value of \a start is \a end[ -1 ].
1040  *  \param [out] arr - a new DataArrayInt that is the "Old to New" renumbering
1041  *         map. The caller is to delete this array using decrRef() as it is no more needed.
1042  *  \return MEDCouplingMesh * - a new instance of MEDCouplingMesh. The caller is to
1043  *         delete this mesh using decrRef() as it is no more needed.
1044  */
1045 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
1046 {
1047   MEDCouplingPointSet *ret=buildPartOfMySelf(start,end,true);
1048   arr=ret->zipCoordsTraducer();
1049   return ret;
1050 }
1051
1052
1053 /*!
1054  * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
1055  */
1056 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
1057 {
1058   double *coords=_coords->getPointer();
1059   int nbNodes=getNumberOfNodes();
1060   Rotate2DAlg(center,angle,nbNodes,coords);
1061 }
1062
1063 /*!
1064  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
1065  * around the center point 'center' and with angle 'angle'.
1066  */
1067 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
1068 {
1069   double cosa=cos(angle);
1070   double sina=sin(angle);
1071   double matrix[4];
1072   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
1073   double tmp[2];
1074   for(int i=0; i<nbNodes; i++)
1075     {
1076       std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
1077       coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
1078       coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
1079     }
1080 }
1081
1082 /// @cond INTERNAL
1083
1084 class DummyClsMCPS
1085 {
1086 public:
1087   static const int MY_SPACEDIM=3;
1088   static const int MY_MESHDIM=2;
1089   typedef int MyConnType;
1090   static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
1091 };
1092
1093 /// @endcond
1094
1095 /*!
1096  * res should be an empty vector before calling this method.
1097  * This method returns all the node coordinates included in _coords which ids are in [startConn;endConn) and put it into 'res' vector.
1098  * If spaceDim==3 a projection will be done for each nodes on the middle plane containing these all nodes in [startConn;endConn).
1099  * And after each projected nodes are moved to Oxy plane in order to consider these nodes as 2D nodes.
1100  */
1101 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
1102 {
1103   const double *coords=_coords->getConstPointer();
1104   int spaceDim=getSpaceDimension();
1105   for(const int *it=startConn;it!=endConn;it++)
1106     res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
1107   if(spaceDim==2)
1108     return ;
1109   if(spaceDim==3)
1110     {
1111       std::vector<double> cpy(res);
1112       int nbNodes=(int)std::distance(startConn,endConn);
1113       INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0.,0.,true);
1114       res.resize(2*nbNodes);
1115       for(int i=0;i<nbNodes;i++)
1116         {
1117           res[2*i]=cpy[3*i];
1118           res[2*i+1]=cpy[3*i+1];
1119         }
1120       return ;
1121     }
1122   throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
1123 }
1124
1125 /*!
1126  * low level method that checks that the 2D cell is not a butterfly cell.
1127  */
1128 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
1129 {
1130   std::size_t nbOfNodes=res.size()/2;
1131   std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
1132   for(std::size_t i=0;i<nbOfNodes;i++)
1133     {
1134       INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
1135       nodes[i]=tmp;
1136     }
1137   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1138   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1139   INTERP_KERNEL::QuadraticPolygon *pol=0;
1140   if(isQuad)
1141     pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
1142   else
1143     pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
1144   bool ret=pol->isButterflyAbs();
1145   delete pol;
1146   return ret;
1147 }