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   if(!vector)
 404     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::translate : NULL input vector !");
 405   if(!_coords)
 406     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::translate : no coordinates set !");
 407   double *coords=_coords->getPointer();
 408   int nbNodes=getNumberOfNodes();
 409   int dim=getSpaceDimension();
 410   for(int i=0; i<nbNodes; i++)
 411     for(int idim=0; idim<dim;idim++)
 412       coords[i*dim+idim]+=vector[idim];
 413   _coords->declareAsNew();
 414   updateTime();
 415 }
 416
 417 /*!
 418  * Non const method that operates a scale on 'this' with 'point' as reference point of scale and with factor 'factor'.
 419  * @param point in array of size getSpaceDimension().
 420  * @param factor factor of the scaling
 421  */
 422 void MEDCouplingPointSet::scale(const double *point, double factor)
 423 {
 424   if(!point)
 425     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::scale : NULL input point !");
 426   if(!_coords)
 427     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::scale : no coordinates set !");
 428   double *coords=_coords->getPointer();
 429   int nbNodes=getNumberOfNodes();
 430   int dim=getSpaceDimension();
 431   for(int i=0;i<nbNodes;i++)
 432     {
 433       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
 434       std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
 435       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
 436     }
 437   _coords->declareAsNew();
 438   updateTime();
 439 }
 440
 441 /*!
 442  * This method is only available for already defined coordinates.
 443  * If not an INTERP_KERNEL::Exception is thrown. The 'newSpaceDim' input must be greater or equal to 1.
 444  * This method simply convert this to newSpaceDim space :
 445  * - by putting a 0. for each \f$ i^{th} \f$ components of each coord of nodes so that i>=getSpaceDim(), if 'newSpaceDim'>getSpaceDimsion()
 446  * - by ignoring each \f$ i^{th} \f$ components of each coord of nodes so that i >= 'newSpaceDim', 'newSpaceDim'<getSpaceDimension()
 447  * If newSpaceDim==getSpaceDim() nothing is done by this method.
 448  */
 449 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue) throw(INTERP_KERNEL::Exception)
 450 {
 451   if(getCoords()==0)
 452     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
 453   if(newSpaceDim<1)
 454     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
 455   int oldSpaceDim=getSpaceDimension();
 456   if(newSpaceDim==oldSpaceDim)
 457     return ;
 458   DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
 459   setCoords(newCoords);
 460   newCoords->decrRef();
 461   updateTime();
 462 }
 463
 464 /*!
 465  * This method try to substitute this->_coords with other._coords if arrays match.
 466  * This method potentially modifies 'this' if it succeeds, otherway an exception is thrown.
 467  */
 468 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
 469 {
 470   if(_coords==other._coords)
 471     return ;
 472   if(!_coords)
 473     throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
 474   if(!other._coords)
 475     throw INTERP_KERNEL::Exception("Other instance has no coords whereas current has !");
 476   if(!_coords->isEqualWithoutConsideringStr(*other._coords,epsilon))
 477     throw INTERP_KERNEL::Exception("Coords are not the same !");
 478   setCoords(other._coords);
 479 }
 480
 481 /*!
 482  * This method duplicates the nodes whose ids are in [\b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd) and put the result of their duplication at the end
 483  * of existing node ids.
 484  *
 485  * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
 486  * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
 487  */
 488 void MEDCouplingPointSet::duplicateNodesInCoords(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd) throw(INTERP_KERNEL::Exception)
 489 {
 490   if(!_coords)
 491     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::duplicateNodesInCoords : no coords set !");
 492   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->selectByTupleIdSafe(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd);
 493   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords2=DataArrayDouble::Aggregate(_coords,newCoords);
 494   setCoords(newCoords2);
 495 }
 496
 497 /*!
 498  * This method is expecting to be called for meshes so that getSpaceDimension() returns 3.
 499  * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from plane
 500  * defined by the point 'pt' and the vector 'vec'.
 501  * @param pt points to an array of size 3 and represents a point that owns to plane.
 502  * @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)
 503  * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
 504  * @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.
 505  */
 506 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
 507 {
 508   if(getSpaceDimension()!=3)
 509     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
 510   int nbOfNodes=getNumberOfNodes();
 511   double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
 512   double deno=sqrt(a*a+b*b+c*c);
 513   const double *work=_coords->getConstPointer();
 514   for(int i=0;i<nbOfNodes;i++)
 515     {
 516       if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
 517         nodes.push_back(i);
 518       work+=3;
 519     }
 520 }
 521
 522 /*!
 523  * This method is expecting to be called for meshes so that getSpaceDimension() returns 2 or 3.
 524  * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from a line
 525  * defined by the point 'pt' and the vector 'vec'. 'pt' and 'vec' are expected to have a dimension equal to space dimension of 'this'
 526  * @param pt points to an array of size this->getSpaceDimension and represents a point that owns to plane.
 527  * @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.
 528  *            But norm must be greater than 10*abs(eps)
 529  * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
 530  * @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.
 531  */
 532 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
 533 {
 534   int spaceDim=getSpaceDimension();
 535   if(spaceDim!=2 && spaceDim!=3)
 536     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
 537   int nbOfNodes=getNumberOfNodes();
 538   double den=0.;
 539   for(int i=0;i<spaceDim;i++)
 540     den+=vec[i]*vec[i];
 541   double deno=sqrt(den);
 542   if(deno<10.*eps)
 543     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
 544   INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
 545   for(int i=0;i<spaceDim;i++)
 546     vecn[i]=vec[i]/deno;
 547   const double *work=_coords->getConstPointer();
 548   if(spaceDim==2)
 549     {
 550       for(int i=0;i<nbOfNodes;i++)
 551         {
 552           if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
 553             nodes.push_back(i);
 554           work+=2;
 555         }
 556     }
 557   else
 558     {
 559       for(int i=0;i<nbOfNodes;i++)
 560         {
 561           double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
 562           double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
 563           double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
 564           if(std::sqrt(a*a+b*b+c*c)<eps)
 565             nodes.push_back(i);
 566           work+=3;
 567         }
 568     }
 569 }
 570
 571 /*!
 572  * merge _coords arrays of m1 and m2 and returns the union. The returned instance is newly created with ref count == 1.
 573  */
 574 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception)
 575 {
 576   int spaceDim=m1->getSpaceDimension();
 577   if(spaceDim!=m2->getSpaceDimension())
 578     throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
 579   return DataArrayDouble::Aggregate(m1->getCoords(),m2->getCoords());
 580 }
 581
 582 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception)
 583 {
 584   if(ms.empty())
 585     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : input array must be NON EMPTY !");
 586   std::vector<const MEDCouplingPointSet *>::const_iterator it=ms.begin();
 587   std::vector<const DataArrayDouble *> coo(ms.size());
 588   int spaceDim=(*it)->getSpaceDimension();
 589   coo[0]=(*it++)->getCoords();
 590   for(int i=1;it!=ms.end();it++,i++)
 591     {
 592       const DataArrayDouble *tmp=(*it)->getCoords();
 593       if(tmp)
 594         {
 595           if((*it)->getSpaceDimension()==spaceDim)
 596             coo[i]=tmp;
 597           else
 598             throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
 599         }
 600       else
 601         throw INTERP_KERNEL::Exception("Empty coords detected during call of MergeNodesArray !");
 602     }
 603   return DataArrayDouble::Aggregate(coo);
 604 }
 605
 606 /*!
 607  * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
 608  * This method is used during unserialization process.
 609  */
 610 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
 611 {
 612   switch(type)
 613     {
 614     case UNSTRUCTURED:
 615       return MEDCouplingUMesh::New();
 616     case UNSTRUCTURED_DESC:
 617       return MEDCouplingUMeshDesc::New();
 618     default:
 619       throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
 620     }
 621 }
 622
 623 /*!
 624  * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
 625  */
 626 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 627 {
 628   int it,order;
 629   double time=getTime(it,order);
 630   if(_coords)
 631     {
 632       int spaceDim=getSpaceDimension();
 633       littleStrings.resize(spaceDim+4);
 634       littleStrings[0]=getName();
 635       littleStrings[1]=getDescription();
 636       littleStrings[2]=_coords->getName();
 637       littleStrings[3]=getTimeUnit();
 638       for(int i=0;i<spaceDim;i++)
 639         littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
 640       tinyInfo.clear();
 641       tinyInfo.push_back(getType());
 642       tinyInfo.push_back(spaceDim);
 643       tinyInfo.push_back(getNumberOfNodes());
 644       tinyInfo.push_back(it);
 645       tinyInfo.push_back(order);
 646       tinyInfoD.push_back(time);
 647     }
 648   else
 649     {
 650       littleStrings.resize(3);
 651       littleStrings[0]=getName();
 652       littleStrings[1]=getDescription();
 653       littleStrings[2]=getTimeUnit();
 654       tinyInfo.clear();
 655       tinyInfo.push_back(getType());
 656       tinyInfo.push_back(-1);
 657       tinyInfo.push_back(-1);
 658       tinyInfo.push_back(it);
 659       tinyInfo.push_back(order);
 660       tinyInfoD.push_back(time);
 661     }
 662 }
 663
 664 /*!
 665  * Third and final step of serialization process.
 666  */
 667 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 668 {
 669   if(_coords)
 670     {
 671       a2=const_cast<DataArrayDouble *>(getCoords());
 672       a2->incrRef();
 673     }
 674   else
 675     a2=0;
 676 }
 677
 678 /*!
 679  * Second step of serialization process.
 680  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 681  */
 682 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 683 {
 684   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 685     {
 686       a2->alloc(tinyInfo[2],tinyInfo[1]);
 687       littleStrings.resize(tinyInfo[1]+4);
 688     }
 689   else
 690     {
 691       littleStrings.resize(3);
 692     }
 693 }
 694
 695 /*!
 696  * Second and final unserialization process.
 697  * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
 698  */
 699 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
 700 {
 701   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
 702     {
 703       setCoords(a2);
 704       setName(littleStrings[0].c_str());
 705       setDescription(littleStrings[1].c_str());
 706       a2->setName(littleStrings[2].c_str());
 707       setTimeUnit(littleStrings[3].c_str());
 708       for(int i=0;i<tinyInfo[1];i++)
 709         getCoords()->setInfoOnComponent(i,littleStrings[i+4].c_str());
 710       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 711     }
 712   else
 713     {
 714       setName(littleStrings[0].c_str());
 715       setDescription(littleStrings[1].c_str());
 716       setTimeUnit(littleStrings[2].c_str());
 717       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 718     }
 719 }
 720
 721 void MEDCouplingPointSet::checkCoherency() const throw(INTERP_KERNEL::Exception)
 722 {
 723   if(!_coords)
 724     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkCoherency : no coordinates set !");
 725 }
 726
 727 /*!
 728  * Intersect Bounding Box given 2 Bounding Boxes.
 729  */
 730 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
 731 {
 732   double* bbtemp = new double[2*dim];
 733   double deltamax=0.0;
 734
 735   for (int i=0; i< dim; i++)
 736     {
 737       double delta = bb1[2*i+1]-bb1[2*i];
 738       if ( delta > deltamax )
 739         {
 740           deltamax = delta ;
 741         }
 742     }
 743   for (int i=0; i<dim; i++)
 744     {
 745       bbtemp[i*2]=bb1[i*2]-deltamax*eps;
 746       bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
 747     }
 748
 749   for (int idim=0; idim < dim; idim++)
 750     {
 751       bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
 752         && (bb2[idim*2]<bbtemp[idim*2+1]) ;
 753       if (!intersects)
 754         {
 755           delete [] bbtemp;
 756           return false;
 757         }
 758     }
 759   delete [] bbtemp;
 760   return true;
 761 }
 762
 763 /*!
 764  * Intersect 2 given Bounding Boxes.
 765  */
 766 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
 767 {
 768   double* bbtemp = new double[2*dim];
 769   double deltamax=0.0;
 770
 771   for (int i=0; i< dim; i++)
 772     {
 773       double delta = bb2[2*i+1]-bb2[2*i];
 774       if ( delta > deltamax )
 775         {
 776           deltamax = delta ;
 777         }
 778     }
 779   for (int i=0; i<dim; i++)
 780     {
 781       bbtemp[i*2]=bb2[i*2]-deltamax*eps;
 782       bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
 783     }
 784
 785   bool intersects = !bb1.isDisjointWith( bbtemp );
 786   delete [] bbtemp;
 787   return intersects;
 788 }
 789
 790 /*!
 791  * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
 792  */
 793 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
 794 {
 795   double *coords=_coords->getPointer();
 796   int nbNodes=getNumberOfNodes();
 797   Rotate3DAlg(center,vect,angle,nbNodes,coords);
 798 }
 799
 800 /*!
 801  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
 802  * around an axe ('center','vect') and with angle 'angle'.
 803  */
 804 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
 805 {
 806   if(!center || !vect)
 807     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : null vector in input !");
 808   double sina=sin(angle);
 809   double cosa=cos(angle);
 810   double vectorNorm[3];
 811   double matrix[9];
 812   double matrixTmp[9];
 813   double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
 814   if(norm<std::numeric_limits<double>::min())
 815     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : magnitude of input vector is too close of 0. !");
 816   std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
 817   //rotation matrix computation
 818   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;
 819   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
 820   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
 821   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
 822   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
 823   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
 824   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
 825   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
 826   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
 827   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
 828   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
 829   //rotation matrix computed.
 830   double tmp[3];
 831   for(int i=0; i<nbNodes; i++)
 832     {
 833       std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
 834       coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
 835       coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
 836       coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
 837     }
 838 }
 839
 840 /*!
 841  * This method implements pure virtual method MEDCouplingMesh::buildPart.
 842  * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end').
 843  * The coords are kept unchanged contrary to pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
 844  * The returned mesh has to be managed by the caller.
 845  */
 846 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
 847 {
 848   return buildPartOfMySelf(start,end,true);
 849 }
 850
 851 /*!
 852  * This method implements pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
 853  * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end') \b and potentially reduces the nodes set
 854  * behind returned mesh. This cause an overhead but it is lesser in memory.
 855  * This method returns an array too. This array allows to the caller to know the mapping between nodeids in 'this' and nodeids in
 856  * returned mesh. This is quite useful for MEDCouplingFieldDouble on nodes for example...
 857  * 'arr' is in old2New format of size ret->getNumberOfCells like MEDCouplingUMesh::zipCoordsTraducer is.
 858  * The returned mesh has to be managed by the caller.
 859  */
 860 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
 861 {
 862   MEDCouplingPointSet *ret=buildPartOfMySelf(start,end,true);
 863   arr=ret->zipCoordsTraducer();
 864   return ret;
 865 }
 866
 867
 868 /*!
 869  * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
 870  */
 871 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
 872 {
 873   double *coords=_coords->getPointer();
 874   int nbNodes=getNumberOfNodes();
 875   Rotate2DAlg(center,angle,nbNodes,coords);
 876 }
 877
 878 /*!
 879  * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
 880  * around the center point 'center' and with angle 'angle'.
 881  */
 882 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
 883 {
 884   double cosa=cos(angle);
 885   double sina=sin(angle);
 886   double matrix[4];
 887   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
 888   double tmp[2];
 889   for(int i=0; i<nbNodes; i++)
 890     {
 891       std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
 892       coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
 893       coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
 894     }
 895 }
 896
 897 /// @cond INTERNAL
 898
 899 class DummyClsMCPS
 900 {
 901 public:
 902   static const int MY_SPACEDIM=3;
 903   static const int MY_MESHDIM=2;
 904   typedef int MyConnType;
 905   static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
 906 };
 907
 908 /// @endcond
 909
 910 /*!
 911  * res should be an empty vector before calling this method.
 912  * This method returns all the node coordinates included in _coords which ids are in [startConn;endConn) and put it into 'res' vector.
 913  * If spaceDim==3 a projection will be done for each nodes on the middle plane containing these all nodes in [startConn;endConn).
 914  * And after each projected nodes are moved to Oxy plane in order to consider these nodes as 2D nodes.
 915  */
 916 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
 917 {
 918   const double *coords=_coords->getConstPointer();
 919   int spaceDim=getSpaceDimension();
 920   for(const int *it=startConn;it!=endConn;it++)
 921     res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
 922   if(spaceDim==2)
 923     return ;
 924   if(spaceDim==3)
 925     {
 926       std::vector<double> cpy(res);
 927       int nbNodes=(int)std::distance(startConn,endConn);
 928       INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0.,0.,true);
 929       res.resize(2*nbNodes);
 930       for(int i=0;i<nbNodes;i++)
 931         {
 932           res[2*i]=cpy[3*i];
 933           res[2*i+1]=cpy[3*i+1];
 934         }
 935       return ;
 936     }
 937   throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
 938 }
 939
 940 /*!
 941  * low level method that checks that the 2D cell is not a butterfly cell.
 942  */
 943 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
 944 {
 945   std::size_t nbOfNodes=res.size()/2;
 946   std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
 947   for(std::size_t i=0;i<nbOfNodes;i++)
 948     {
 949       INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
 950       nodes[i]=tmp;
 951     }
 952   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
 953   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
 954   INTERP_KERNEL::QuadraticPolygon *pol=0;
 955   if(isQuad)
 956     pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
 957   else
 958     pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
 959   bool ret=pol->isButterflyAbs();
 960   delete pol;
 961   return ret;
 962 }