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 void MEDCouplingPointSet::setCoords(const DataArrayDouble *coords)
  79 {
  80   if( coords != _coords )
  81     {
  82       if (_coords)
  83         _coords->decrRef();
  84       _coords=const_cast<DataArrayDouble *>(coords);
  85       if(_coords)
  86         _coords->incrRef();
  87       declareAsNew();
  88     }
  89 }
  90
  91 DataArrayDouble *MEDCouplingPointSet::getCoordinatesAndOwner() const
  92 {
  93   if(_coords)
  94     _coords->incrRef();
  95   return _coords;
  96 }
  97
  98 /*!
  99  * This method copyies all tiny strings from other (name and components name).
 100  * @throw if other and this have not same mesh type.
 101  */
 102 void MEDCouplingPointSet::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
 103 {
 104   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
 105   if(!otherC)
 106     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::copyTinyStringsFrom : meshes have not same type !");
 107   MEDCouplingMesh::copyTinyStringsFrom(other);
 108   if(_coords && otherC->_coords)
 109     _coords->copyStringInfoFrom(*otherC->_coords);
 110 }
 111
 112 bool MEDCouplingPointSet::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
 113 {
 114   if(!other)
 115     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::isEqualIfNotWhy : null mesh instance in input !");
 116   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
 117   if(!otherC)
 118     {
 119       reason="mesh given in input is not castable in MEDCouplingPointSet !";
 120       return false;
 121     }
 122   if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
 123     return false;
 124   if(!areCoordsEqualIfNotWhy(*otherC,prec,reason))
 125     return false;
 126   return true;
 127 }
 128
 129 bool MEDCouplingPointSet::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 130 {
 131   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
 132   if(!otherC)
 133     return false;
 134   if(!areCoordsEqualWithoutConsideringStr(*otherC,prec))
 135     return false;
 136   return true;
 137 }
 138
 139 bool MEDCouplingPointSet::areCoordsEqualIfNotWhy(const MEDCouplingPointSet& other, double prec, std::string& reason) const
 140 {
 141   if(_coords==0 && other._coords==0)
 142     return true;
 143   if(_coords==0 || other._coords==0)
 144     {
 145       reason="Only one PointSet between the two this and other has coordinate defined !";
 146       return false;
 147     }
 148   if(_coords==other._coords)
 149     return true;
 150   bool ret=_coords->isEqualIfNotWhy(*other._coords,prec,reason);
 151   if(!ret)
 152     reason.insert(0,"Coordinates DataArray do not match : ");
 153   return ret;
 154 }
 155
 156 bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const
 157 {
 158   std::string tmp;
 159   return areCoordsEqualIfNotWhy(other,prec,tmp);
 160 }
 161
 162 bool MEDCouplingPointSet::areCoordsEqualWithoutConsideringStr(const MEDCouplingPointSet& other, double prec) const
 163 {
 164   if(_coords==0 && other._coords==0)
 165     return true;
 166   if(_coords==0 || other._coords==0)
 167     return false;
 168   if(_coords==other._coords)
 169     return true;
 170   return _coords->isEqualWithoutConsideringStr(*other._coords,prec);
 171 }
 172
 173 /*!
 174  * Returns coordinates of node with id 'nodeId' and append it in 'coo'.
 175  */
 176 void MEDCouplingPointSet::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
 177 {
 178   if(!_coords)
 179     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCoordinatesOfNode : no coordinates array set !");
 180   int nbNodes=getNumberOfNodes();
 181   if(nodeId>=0 && nodeId<nbNodes)
 182     {
 183       const double *cooPtr=_coords->getConstPointer();
 184       int spaceDim=getSpaceDimension();
 185       coo.insert(coo.end(),cooPtr+spaceDim*nodeId,cooPtr+spaceDim*(nodeId+1));
 186     }
 187   else
 188     {
 189       std::ostringstream oss; oss << "MEDCouplingPointSet::getCoordinatesOfNode : request of nodeId \"" << nodeId << "\" but it should be in [0,"<< nbNodes << ") !";
 190       throw INTERP_KERNEL::Exception(oss.str().c_str());
 191     }
 192 }
 193
 194 /*!
 195  * This method is typically the base method used for implementation of mergeNodes. This method computes this permutation array using as input,
 196  * This method is const ! So this method simply computes the array, no permutation of nodes is done.
 197  * a precision 'precision' and a 'limitNodeId' that is the node id so that every nodes which id is strictly lower than 'limitNodeId' will not be merged.
 198  * To desactivate this advanced feature put -1 to this argument.
 199  * @param areNodesMerged output parameter that states if some nodes have been "merged" in returned array
 200  * @param newNbOfNodes output parameter too this is the maximal id in returned array to avoid to recompute it.
 201  */
 202 DataArrayInt *MEDCouplingPointSet::buildPermArrayForMergeNode(double precision, int limitNodeId, bool& areNodesMerged, int& newNbOfNodes) const
 203 {
 204   DataArrayInt *comm,*commI;
 205   findCommonNodes(precision,limitNodeId,comm,commI);
 206   int oldNbOfNodes=getNumberOfNodes();
 207   DataArrayInt *ret=buildNewNumberingFromCommonNodesFormat(comm,commI,newNbOfNodes);
 208   areNodesMerged=(oldNbOfNodes!=newNbOfNodes);
 209   comm->decrRef();
 210   commI->decrRef();
 211   return ret;
 212 }
 213
 214 /*!
 215  * This methods searches for each node if there are any nodes in _coords that are less far than 'prec' from n1. if any, these nodes are stored in out params
 216  * comm and commIndex.
 217  * @param limitNodeId is the limit node id. All nodes which id is strictly lower than 'limitNodeId' will not be merged each other.
 218  * @param comm out parameter (not inout)
 219  * @param commIndex out parameter (not inout)
 220  */
 221 void MEDCouplingPointSet::findCommonNodes(double prec, int limitNodeId, DataArrayInt *&comm, DataArrayInt *&commIndex) const
 222 {
 223   if(!_coords)
 224     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findCommonNodes : no coords specified !");
 225   _coords->findCommonTuples(prec,limitNodeId,comm,commIndex);
 226 }
 227
 228 std::vector<int> MEDCouplingPointSet::getNodeIdsNearPoint(const double *pos, double eps) const throw(INTERP_KERNEL::Exception)
 229 {
 230   std::vector<int> c,cI;
 231   getNodeIdsNearPoints(pos,1,eps,c,cI);
 232   return c;
 233 }
 234
 235 /*!
 236  * Given a point given by its position 'pos' this method finds the set of node ids that are a a distance lower than eps.
 237  * Position 'pos' is expected to be of size getSpaceDimension()*nbOfNodes. If not the behabiour is not warranted.
 238  * This method throws an exception if no coordiantes are set.
 239  */
 240 void MEDCouplingPointSet::getNodeIdsNearPoints(const double *pos, int nbOfNodes, double eps, std::vector<int>& c, std::vector<int>& cI) const throw(INTERP_KERNEL::Exception)
 241 {
 242   if(!_coords)
 243     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : no coordiantes set !");
 244   int spaceDim=getSpaceDimension();
 245   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> points=DataArrayDouble::New();
 246   points->useArray(pos,false,CPP_DEALLOC,nbOfNodes,spaceDim);
 247   _coords->computeTupleIdsNearTuples(points,eps,c,cI);
 248 }
 249
 250 /*!
 251  * @param comm in param in the same format than one returned by findCommonNodes method.
 252  * @param commI in param in the same format than one returned by findCommonNodes method.
 253  * @return the old to new correspondance array.
 254  */
 255 DataArrayInt *MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat(const DataArrayInt *comm, const DataArrayInt *commIndex,
 256                                                                           int& newNbOfNodes) const
 257 {
 258   if(!_coords)
 259     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat : no coords specified !");
 260   return DataArrayInt::BuildOld2NewArrayFromSurjectiveFormat2(getNumberOfNodes(),comm,commIndex,newNbOfNodes);
 261 }
 262
 263 /*
 264  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
 265  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
 266  * This value is asked because often known by the caller of this method.
 267  * @param newNodeNumbers array specifying the new numbering in old2New convention..
 268  * @param newNbOfNodes the new number of nodes.
 269  */
 270 void MEDCouplingPointSet::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
 271 {
 272   if(!_coords)
 273     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::renumberNodes : no coords specified !");
 274   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->renumberAndReduce(newNodeNumbers,newNbOfNodes);
 275   setCoords(newCoords);
 276 }
 277
 278 /*
 279  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
 280  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
 281  * This value is asked because often known by the caller of this method.
 282  * Contrary to ParaMEDMEM::MEDCouplingPointSet::renumberNodes method for merged nodes the barycenter of them is computed here.
 283  *
 284  * @param newNodeNumbers array specifying the new numbering.
 285  * @param newNbOfNodes the new number of nodes.
 286  */
 287 void MEDCouplingPointSet::renumberNodes2(const int *newNodeNumbers, int newNbOfNodes)
 288 {
 289   DataArrayDouble *newCoords=DataArrayDouble::New();
 290   std::vector<int> div(newNbOfNodes);
 291   int spaceDim=getSpaceDimension();
 292   newCoords->alloc(newNbOfNodes,spaceDim);
 293   newCoords->copyStringInfoFrom(*_coords);
 294   newCoords->fillWithZero();
 295   int oldNbOfNodes=getNumberOfNodes();
 296   double *ptToFill=newCoords->getPointer();
 297   const double *oldCoordsPtr=_coords->getConstPointer();
 298   for(int i=0;i<oldNbOfNodes;i++)
 299     {
 300       std::transform(oldCoordsPtr+i*spaceDim,oldCoordsPtr+(i+1)*spaceDim,ptToFill+newNodeNumbers[i]*spaceDim,
 301                      ptToFill+newNodeNumbers[i]*spaceDim,std::plus<double>());
 302       div[newNodeNumbers[i]]++;
 303     }
 304   for(int i=0;i<newNbOfNodes;i++)
 305     ptToFill=std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)div[i]));
 306   setCoords(newCoords);
 307   newCoords->decrRef();
 308 }
 309
 310 /*!
 311  * This method fills bbox params like that : bbox[0]=XMin, bbox[1]=XMax, bbox[2]=YMin...
 312  * The returned bounding box is arranged along trihedron.
 313  * @param bbox out array of size 2*this->getSpaceDimension().
 314  */
 315 void MEDCouplingPointSet::getBoundingBox(double *bbox) const throw(INTERP_KERNEL::Exception)
 316 {
 317   if(!_coords)
 318     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getBoundingBox : Coordinates not set !");
 319   _coords->getMinMaxPerComponent(bbox);
 320 }
 321
 322 /*!
 323  * This method removes useless nodes in coords.
 324  */
 325 void MEDCouplingPointSet::zipCoords()
 326 {
 327   checkFullyDefined();
 328   DataArrayInt *traducer=zipCoordsTraducer();
 329   traducer->decrRef();
 330 }
 331
 332 struct MEDCouplingCompAbs
 333 {
 334   bool operator()(double x, double y) { return std::abs(x)<std::abs(y);}
 335 };
 336
 337 /*!
 338  * This method expects that _coords attribute is set.
 339  * @return the carateristic dimension of point set. This caracteristic dimension is the max of difference
 340  * @exception If _coords attribute not set.
 341  */
 342 double MEDCouplingPointSet::getCaracteristicDimension() const
 343 {
 344   if(!_coords)
 345     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCaracteristicDimension : Coordinates not set !");
 346   const double *coords=_coords->getConstPointer();
 347   int nbOfValues=_coords->getNbOfElems();
 348   return std::abs(*std::max_element(coords,coords+nbOfValues,MEDCouplingCompAbs()));
 349 }
 350
 351 /*!
 352  * 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
 353  * around origin of 'radius' 1.
 354  *
 355  * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
 356  *
 357  * \warning this method is non const and alterates coordinates in \b this without modifying.
 358  */
 359 void MEDCouplingPointSet::recenterForMaxPrecision(double eps) throw(INTERP_KERNEL::Exception)
 360 {
 361   if(!_coords)
 362     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::recenterForMaxPrecision : Coordinates not set !");
 363   _coords->recenterForMaxPrecision(eps);
 364   updateTime();
 365 }
 366
 367 /*!
 368  * Non const method that operates a rotation of 'this'.
 369  * If spaceDim==2 'vector' parameter is ignored (and could be 0) and the rotation is done around 'center' with angle specified by 'angle'.
 370  * If spaceDim==3 the rotation axe is defined by ('center','vector') and the angle is 'angle'.
 371  * @param center an array of size getSpaceDimension().
 372  * @param vector in array of size getSpaceDimension().
 373  * @param angle angle of rotation in radian.
 374  */
 375 void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
 376 {
 377   int spaceDim=getSpaceDimension();
 378   if(spaceDim==3)
 379     rotate3D(center,vector,angle);
 380   else if(spaceDim==2)
 381     rotate2D(center,angle);
 382   else
 383     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
 384   _coords->declareAsNew();
 385   updateTime();
 386 }
 387
 388 /*!
 389  * Non const method that operates a translation of 'this'.
 390  * @param vector in array of size getSpaceDimension().
 391  */
 392 void MEDCouplingPointSet::translate(const double *vector)
 393 {
 394   double *coords=_coords->getPointer();
 395   int nbNodes=getNumberOfNodes();
 396   int dim=getSpaceDimension();
 397   for(int i=0; i<nbNodes; i++)
 398     for(int idim=0; idim<dim;idim++)
 399       coords[i*dim+idim]+=vector[idim];
 400   _coords->declareAsNew();
 401   updateTime();
 402 }
 403
 404 /*!
 405  * Non const method that operates a scale on 'this' with 'point' as reference point of scale and with factor 'factor'.
 406  * @param point in array of size getSpaceDimension().
 407  * @param factor factor of the scaling
 408  */
 409 void MEDCouplingPointSet::scale(const double *point, double factor)
 410 {
 411   double *coords=_coords->getPointer();
 412   int nbNodes=getNumberOfNodes();
 413   int dim=getSpaceDimension();
 414   double *tmp=new double[dim];
 415   for(int i=0;i<nbNodes;i++)
 416     {
 417       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
 418       std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
 419       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
 420     }
 421   delete [] tmp;
 422   _coords->declareAsNew();
 423   updateTime();
 424 }
 425
 426 /*!
 427  * This method is only available for already defined coordinates.
 428  * If not an INTERP_KERNEL::Exception is thrown. The 'newSpaceDim' input must be greater or equal to 1.
 429  * This method simply convert this to newSpaceDim space :
 430  * - by putting a 0. for each \f$ i^{th} \f$ components of each coord of nodes so that i>=getSpaceDim(), if 'newSpaceDim'>getSpaceDimsion()
 431  * - by ignoring each \f$ i^{th} \f$ components of each coord of nodes so that i >= 'newSpaceDim', 'newSpaceDim'<getSpaceDimension()
 432  * If newSpaceDim==getSpaceDim() nothing is done by this method.
 433  */
 434 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue) throw(INTERP_KERNEL::Exception)
 435 {
 436   if(getCoords()==0)
 437     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
 438   if(newSpaceDim<1)
 439     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
 440   int oldSpaceDim=getSpaceDimension();
 441   if(newSpaceDim==oldSpaceDim)
 442     return ;
 443   DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
 444   setCoords(newCoords);
 445   newCoords->decrRef();
 446   updateTime();
 447 }
 448
 449 /*!
 450  * This method try to substitute this->_coords with other._coords if arrays match.
 451  * This method potentially modifies 'this' if it succeeds, otherway an exception is thrown.
 452  */
 453 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
 454 {
 455   if(_coords==other._coords)
 456     return ;
 457   if(!_coords)
 458     throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
 459   if(!other._coords)
 460     throw INTERP_KERNEL::Exception("Other instance has no coords whereas current has !");
 461   if(!_coords->isEqualWithoutConsideringStr(*other._coords,epsilon))
 462     throw INTERP_KERNEL::Exception("Coords are not the same !");
 463   setCoords(other._coords);
 464 }
 465
 466 /*!
 467  * This method duplicates the nodes whose ids are in [\b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd) and put the result of their duplication at the end
 468  * of existing node ids.
 469  *
 470  * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
 471  * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
 472  */
 473 void MEDCouplingPointSet::duplicateNodesInCoords(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd) throw(INTERP_KERNEL::Exception)
 474 {
 475   if(!_coords)
 476     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::duplicateNodesInCoords : no coords set !");
 477   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->selectByTupleIdSafe(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd);
 478   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords2=DataArrayDouble::Aggregate(_coords,newCoords);
 479   setCoords(newCoords2);
 480 }
 481
 482 /*!
 483  * This method is expecting to be called for meshes so that getSpaceDimension() returns 3.
 484  * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from plane
 485  * defined by the point 'pt' and the vector 'vec'.
 486  * @param pt points to an array of size 3 and represents a point that owns to plane.
 487  * @param vec points to an array of size 3 and represents the normal vector of the plane. The norm of the vector is not compulsory equal to 1. But norm must be greater than 10*abs(eps)
 488  * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
 489  * @param nodes is the output of the method. The vector is not compulsory empty before call. The nodes that fulfills the condition will be added at the end of the nodes.
 490  */
 491 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
 492 {
 493   if(getSpaceDimension()!=3)
 494     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
 495   int nbOfNodes=getNumberOfNodes();
 496   double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
 497   double deno=sqrt(a*a+b*b+c*c);
 498   const double *work=_coords->getConstPointer();
 499   for(int i=0;i<nbOfNodes;i++)
 500     {
 501       if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
 502         nodes.push_back(i);
 503       work+=3;
 504     }
 505 }
 506
 507 /*!
 508  * This method is expecting to be called for meshes so that getSpaceDimension() returns 2 or 3.
 509  * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from a line
 510  * defined by the point 'pt' and the vector 'vec'. 'pt' and 'vec' are expected to have a dimension equal to space dimension of 'this'
 511  * @param pt points to an array of size this->getSpaceDimension and represents a point that owns to plane.
 512  * @param vec points to an array of size this->getSpaceDimension and represents the direction vector of the line. The norm of the vector is not compulsory equal to 1.
 513  *            But norm must be greater than 10*abs(eps)
 514  * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
 515  * @param nodes is the output of the method. The vector is not compulsory empty before call. The nodes that fulfills the condition will be added at the end of the nodes.
 516  */
 517 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
 518 {
 519   int spaceDim=getSpaceDimension();
 520   if(spaceDim!=2 && spaceDim!=3)
 521     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
 522   int nbOfNodes=getNumberOfNodes();
 523   double den=0.;
 524   for(int i=0;i<spaceDim;i++)
 525     den+=vec[i]*vec[i];
 526   double deno=sqrt(den);
 527   if(deno<10.*eps)
 528     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
 529   INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
 530   for(int i=0;i<spaceDim;i++)
 531     vecn[i]=vec[i]/deno;
 532   const double *work=_coords->getConstPointer();
 533   if(spaceDim==2)
 534     {
 535       for(int i=0;i<nbOfNodes;i++)
 536         {
 537           if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
 538             nodes.push_back(i);
 539           work+=2;
 540         }
 541     }
 542   else
 543     {
 544       for(int i=0;i<nbOfNodes;i++)
 545         {
 546           double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
 547           double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
 548           double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
 549           if(std::sqrt(a*a+b*b+c*c)<eps)
 550             nodes.push_back(i);
 551           work+=3;
 552         }
 553     }
 554 }
 555
 556 /*!
 557  * merge _coords arrays of m1 and m2 and returns the union. The returned instance is newly created with ref count == 1.
 558  */
 559 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception)
 560 {
 561   int spaceDim=m1->getSpaceDimension();
 562   if(spaceDim!=m2->getSpaceDimension())
 563     throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
 564   return DataArrayDouble::Aggregate(m1->getCoords(),m2->getCoords());
 565 }
 566
 567 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception)
 568 {
 569   if(ms.empty())
 570     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : input array must be NON EMPTY !");
 571   std::vector<const MEDCouplingPointSet *>::const_iterator it=ms.begin();
 572   std::vector<const DataArrayDouble *> coo(ms.size());
 573   int spaceDim=(*it)->getSpaceDimension();
 574   coo[0]=(*it++)->getCoords();
 575   for(int i=1;it!=ms.end();it++,i++)
 576     {
 577       const DataArrayDouble *tmp=(*it)->getCoords();
 578       if(tmp)
 579         {
 580           if((*it)->getSpaceDimension()==spaceDim)
 581             coo[i]=tmp;
 582           else
 583             throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
 584         }
 585       else
 586         throw INTERP_KERNEL::Exception("Empty coords detected during call of MergeNodesArray !");
 587     }
 588   return DataArrayDouble::Aggregate(coo);
 589 }
 590
 591 /*!
 592  * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
 593  * This method is used during unserialization process.
 594  */
 595 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
 596 {
 597   switch(type)
 598     {
 599     case UNSTRUCTURED:
 600       return MEDCouplingUMesh::New();
 601     case UNSTRUCTURED_DESC:
 602       return MEDCouplingUMeshDesc::New();
 603     default:
 604       throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
 605     }
 606 }
 607
 608 /*!
 609  * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
 610  */
 611 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 612 {
 613   int it,order;
 614   double time=getTime(it,order);
 615   if(_coords)
 616     {
 617       int spaceDim=getSpaceDimension();
 618       littleStrings.resize(spaceDim+4);
 619       littleStrings[0]=getName();
 620       littleStrings[1]=getDescription();
 621       littleStrings[2]=_coords->getName();
 622       littleStrings[3]=getTimeUnit();
 623       for(int i=0;i<spaceDim;i++)
 624         littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
 625       tinyInfo.clear();
 626       tinyInfo.push_back(getType());
 627       tinyInfo.push_back(spaceDim);
 628       tinyInfo.push_back(getNumberOfNodes());
 629       tinyInfo.push_back(it);
 630       tinyInfo.push_back(order);
 631       tinyInfoD.push_back(time);
 632     }
 633   else
 634     {
 635       littleStrings.resize(3);
 636       littleStrings[0]=getName();
 637       littleStrings[1]=getDescription();
 638       littleStrings[2]=getTimeUnit();
 639       tinyInfo.clear();
 640       tinyInfo.push_back(getType());
 641       tinyInfo.push_back(-1);
 642       tinyInfo.push_back(-1);
 643       tinyInfo.push_back(it);
 644       tinyInfo.push_back(order);
 645       tinyInfoD.push_back(time);
 646     }
 647 }
 648
 649 /*!
 650  * Third and final step of serialization process.
 651  */
 652 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 653 {
 654   if(_coords)
 655     {
 656       a2=const_cast<DataArrayDouble *>(getCoords());
 657       a2->incrRef();
 658     }
 659   else
 660     a2=0;
 661 }
 662
 663 /*!
 664  * Second step of serialization process.
 665  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 666  */
 667 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 668 {
 669   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 670     {
 671       a2->alloc(tinyInfo[2],tinyInfo[1]);
 672       littleStrings.resize(tinyInfo[1]+4);
 673     }
 674   else
 675     {
 676       littleStrings.resize(3);
 677     }
 678 }
 679
 680 /*!
 681  * Second and final unserialization process.
 682  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 683  */
 684 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
 685 {
 686   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 687     {
 688       setCoords(a2);
 689       setName(littleStrings[0].c_str());
 690       setDescription(littleStrings[1].c_str());
 691       a2->setName(littleStrings[2].c_str());
 692       setTimeUnit(littleStrings[3].c_str());
 693       for(int i=0;i<tinyInfo[1];i++)
 694         getCoords()->setInfoOnComponent(i,littleStrings[i+4].c_str());
 695       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 696     }
 697   else
 698     {
 699       setName(littleStrings[0].c_str());
 700       setDescription(littleStrings[1].c_str());
 701       setTimeUnit(littleStrings[2].c_str());
 702       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 703     }
 704 }
 705
 706 /*!
 707  * Intersect Bounding Box given 2 Bounding Boxes.
 708  */
 709 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
 710 {
 711   double* bbtemp = new double[2*dim];
 712   double deltamax=0.0;
 713
 714   for (int i=0; i< dim; i++)
 715     {
 716       double delta = bb1[2*i+1]-bb1[2*i];
 717       if ( delta > deltamax )
 718         {
 719           deltamax = delta ;
 720         }
 721     }
 722   for (int i=0; i<dim; i++)
 723     {
 724       bbtemp[i*2]=bb1[i*2]-deltamax*eps;
 725       bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
 726     }
 727
 728   for (int idim=0; idim < dim; idim++)
 729     {
 730       bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
 731         && (bb2[idim*2]<bbtemp[idim*2+1]) ;
 732       if (!intersects)
 733         {
 734           delete [] bbtemp;
 735           return false;
 736         }
 737     }
 738   delete [] bbtemp;
 739   return true;
 740 }
 741
 742 /*!
 743  * Intersect 2 given Bounding Boxes.
 744  */
 745 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
 746 {
 747   double* bbtemp = new double[2*dim];
 748   double deltamax=0.0;
 749
 750   for (int i=0; i< dim; i++)
 751     {
 752       double delta = bb2[2*i+1]-bb2[2*i];
 753       if ( delta > deltamax )
 754         {
 755           deltamax = delta ;
 756         }
 757     }
 758   for (int i=0; i<dim; i++)
 759     {
 760       bbtemp[i*2]=bb2[i*2]-deltamax*eps;
 761       bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
 762     }
 763
 764   bool intersects = !bb1.isDisjointWith( bbtemp );
 765   delete [] bbtemp;
 766   return intersects;
 767 }
 768
 769 /*!
 770  * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
 771  */
 772 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
 773 {
 774   double *coords=_coords->getPointer();
 775   int nbNodes=getNumberOfNodes();
 776   Rotate3DAlg(center,vect,angle,nbNodes,coords);
 777 }
 778
 779 /*!
 780  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
 781  * around an axe ('center','vect') and with angle 'angle'.
 782  */
 783 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
 784 {
 785   if(!center || !vect)
 786     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : null vector in input !");
 787   double sina=sin(angle);
 788   double cosa=cos(angle);
 789   double vectorNorm[3];
 790   double matrix[9];
 791   double matrixTmp[9];
 792   double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
 793   std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
 794   //rotation matrix computation
 795   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;
 796   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
 797   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
 798   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
 799   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
 800   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
 801   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
 802   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
 803   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
 804   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
 805   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
 806   //rotation matrix computed.
 807   double tmp[3];
 808   for(int i=0; i<nbNodes; i++)
 809     {
 810       std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
 811       coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
 812       coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
 813       coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
 814     }
 815 }
 816
 817 /*!
 818  * This method implements pure virtual method MEDCouplingMesh::buildPart.
 819  * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end').
 820  * The coords are kept unchanged contrary to pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
 821  * The returned mesh has to be managed by the caller.
 822  */
 823 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
 824 {
 825   return buildPartOfMySelf(start,end,true);
 826 }
 827
 828 /*!
 829  * This method implements pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
 830  * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end') \b and potentially reduces the nodes set
 831  * behind returned mesh. This cause an overhead but it is lesser in memory.
 832  * This method returns an array too. This array allows to the caller to know the mapping between nodeids in 'this' and nodeids in
 833  * returned mesh. This is quite usefull for MEDCouplingFieldDouble on nodes for example...
 834  * 'arr' is in old2New format of size ret->getNumberOfCells like MEDCouplingUMesh::zipCoordsTraducer is.
 835  * The returned mesh has to be managed by the caller.
 836  */
 837 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
 838 {
 839   MEDCouplingPointSet *ret=buildPartOfMySelf(start,end,true);
 840   arr=ret->zipCoordsTraducer();
 841   return ret;
 842 }
 843
 844
 845 /*!
 846  * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
 847  */
 848 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
 849 {
 850   double *coords=_coords->getPointer();
 851   int nbNodes=getNumberOfNodes();
 852   Rotate2DAlg(center,angle,nbNodes,coords);
 853 }
 854
 855 /*!
 856  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
 857  * around the center point 'center' and with angle 'angle'.
 858  */
 859 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
 860 {
 861   double cosa=cos(angle);
 862   double sina=sin(angle);
 863   double matrix[4];
 864   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
 865   double tmp[2];
 866   for(int i=0; i<nbNodes; i++)
 867     {
 868       std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
 869       coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
 870       coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
 871     }
 872 }
 873
 874 /// @cond INTERNAL
 875
 876 class DummyClsMCPS
 877 {
 878 public:
 879   static const int MY_SPACEDIM=3;
 880   static const int MY_MESHDIM=2;
 881   typedef int MyConnType;
 882   static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
 883 };
 884
 885 /// @endcond
 886
 887 /*!
 888  * res should be an empty vector before calling this method.
 889  * This method returns all the node coordinates included in _coords which ids are in [startConn;endConn) and put it into 'res' vector.
 890  * If spaceDim==3 a projection will be done for each nodes on the middle plane containing these all nodes in [startConn;endConn).
 891  * And after each projected nodes are moved to Oxy plane in order to consider these nodes as 2D nodes.
 892  */
 893 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
 894 {
 895   const double *coords=_coords->getConstPointer();
 896   int spaceDim=getSpaceDimension();
 897   for(const int *it=startConn;it!=endConn;it++)
 898     res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
 899   if(spaceDim==2)
 900     return ;
 901   if(spaceDim==3)
 902     {
 903       std::vector<double> cpy(res);
 904       int nbNodes=(int)std::distance(startConn,endConn);
 905       INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0.,0.,true);
 906       res.resize(2*nbNodes);
 907       for(int i=0;i<nbNodes;i++)
 908         {
 909           res[2*i]=cpy[3*i];
 910           res[2*i+1]=cpy[3*i+1];
 911         }
 912       return ;
 913     }
 914   throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
 915 }
 916
 917 /*!
 918  * low level method that checks that the 2D cell is not a butterfly cell.
 919  */
 920 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
 921 {
 922   std::size_t nbOfNodes=res.size()/2;
 923   std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
 924   for(std::size_t i=0;i<nbOfNodes;i++)
 925     {
 926       INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
 927       nodes[i]=tmp;
 928     }
 929   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
 930   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
 931   INTERP_KERNEL::QuadraticPolygon *pol=0;
 932   if(isQuad)
 933     pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
 934   else
 935     pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
 936   bool ret=pol->isButterflyAbs();
 937   delete pol;
 938   return ret;
 939 }