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