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