Back to index

salome-med  6.5.0
MEDCouplingPointSet.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDCouplingPointSet.hxx"
00021 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
00022 #include "MEDCouplingUMesh.hxx"
00023 #include "MEDCouplingUMeshDesc.hxx"
00024 #include "MEDCouplingMemArray.hxx"
00025 #include "MEDCouplingPointSet.txx"
00026 #include "PlanarIntersector.txx"
00027 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
00028 #include "InterpKernelGeo2DNode.hxx"
00029 #include "DirectedBoundingBox.hxx"
00030 #include "InterpKernelAutoPtr.hxx"
00031 
00032 #include <cmath>
00033 #include <limits>
00034 #include <numeric>
00035 
00036 using namespace ParaMEDMEM;
00037 
00038 MEDCouplingPointSet::MEDCouplingPointSet():_coords(0)
00039 {
00040 }
00041 
00042 MEDCouplingPointSet::MEDCouplingPointSet(const MEDCouplingPointSet& other, bool deepCopy):MEDCouplingMesh(other),_coords(0)
00043 {
00044   if(other._coords)
00045     _coords=other._coords->performCpy(deepCopy);
00046 }
00047 
00048 MEDCouplingPointSet::~MEDCouplingPointSet()
00049 {
00050   if(_coords)
00051     _coords->decrRef();
00052 }
00053 
00054 int MEDCouplingPointSet::getNumberOfNodes() const
00055 {
00056   if(_coords)
00057     return _coords->getNumberOfTuples();
00058   else
00059     throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !");
00060 }
00061 
00062 int MEDCouplingPointSet::getSpaceDimension() const
00063 {
00064   if(_coords)
00065     return _coords->getNumberOfComponents();
00066   else
00067     throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !");
00068 }
00069 
00070 void MEDCouplingPointSet::updateTime() const
00071 {
00072   if(_coords)
00073     {
00074       updateTimeWith(*_coords);
00075     }
00076 }
00077 
00078 void MEDCouplingPointSet::setCoords(const DataArrayDouble *coords)
00079 {
00080   if( coords != _coords )
00081     {
00082       if (_coords)
00083         _coords->decrRef();
00084       _coords=const_cast<DataArrayDouble *>(coords);
00085       if(_coords)
00086         _coords->incrRef();
00087       declareAsNew();
00088     }
00089 }
00090 
00091 DataArrayDouble *MEDCouplingPointSet::getCoordinatesAndOwner() const
00092 {
00093   if(_coords)
00094     _coords->incrRef();
00095   return _coords;
00096 }
00097 
00102 void MEDCouplingPointSet::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
00103 {
00104   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
00105   if(!otherC)
00106     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::copyTinyStringsFrom : meshes have not same type !");
00107   MEDCouplingMesh::copyTinyStringsFrom(other);
00108   if(_coords && otherC->_coords)
00109     _coords->copyStringInfoFrom(*otherC->_coords);
00110 }
00111 
00112 bool MEDCouplingPointSet::isEqual(const MEDCouplingMesh *other, double prec) const
00113 {
00114   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
00115   if(!otherC)
00116     return false;
00117   if(!MEDCouplingMesh::isEqual(other,prec))
00118     return false;
00119   if(!areCoordsEqual(*otherC,prec))
00120     return false;
00121   return true;
00122 }
00123 
00124 bool MEDCouplingPointSet::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
00125 {
00126   const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
00127   if(!otherC)
00128     return false;
00129   if(!areCoordsEqualWithoutConsideringStr(*otherC,prec))
00130     return false;
00131   return true;
00132 }
00133 
00134 bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const
00135 {
00136   if(_coords==0 && other._coords==0)
00137     return true;
00138   if(_coords==0 || other._coords==0)
00139     return false;
00140   if(_coords==other._coords)
00141     return true;
00142   return _coords->isEqual(*other._coords,prec);
00143 }
00144 
00145 bool MEDCouplingPointSet::areCoordsEqualWithoutConsideringStr(const MEDCouplingPointSet& other, double prec) const
00146 {
00147   if(_coords==0 && other._coords==0)
00148     return true;
00149   if(_coords==0 || other._coords==0)
00150     return false;
00151   if(_coords==other._coords)
00152     return true;
00153   return _coords->isEqualWithoutConsideringStr(*other._coords,prec);
00154 }
00155 
00159 void MEDCouplingPointSet::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
00160 {
00161   if(!_coords)
00162     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCoordinatesOfNode : no coordinates array set !");
00163   int nbNodes=getNumberOfNodes();
00164   if(nodeId>=0 && nodeId<nbNodes)
00165     {
00166       const double *cooPtr=_coords->getConstPointer();
00167       int spaceDim=getSpaceDimension();
00168       coo.insert(coo.end(),cooPtr+spaceDim*nodeId,cooPtr+spaceDim*(nodeId+1));
00169     }
00170   else
00171     {
00172       std::ostringstream oss; oss << "MEDCouplingPointSet::getCoordinatesOfNode : request of nodeId \"" << nodeId << "\" but it should be in [0,"<< nbNodes << ") !";
00173       throw INTERP_KERNEL::Exception(oss.str().c_str());
00174     }
00175 }
00176 
00185 DataArrayInt *MEDCouplingPointSet::buildPermArrayForMergeNode(double precision, int limitNodeId, bool& areNodesMerged, int& newNbOfNodes) const
00186 {
00187   DataArrayInt *comm,*commI;
00188   findCommonNodes(precision,limitNodeId,comm,commI);
00189   int oldNbOfNodes=getNumberOfNodes();
00190   DataArrayInt *ret=buildNewNumberingFromCommonNodesFormat(comm,commI,newNbOfNodes);
00191   areNodesMerged=(oldNbOfNodes!=newNbOfNodes);
00192   comm->decrRef();
00193   commI->decrRef();
00194   return ret;
00195 }
00196 
00204 void MEDCouplingPointSet::findCommonNodes(double prec, int limitNodeId, DataArrayInt *&comm, DataArrayInt *&commIndex) const
00205 {
00206   if(!_coords)
00207     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findCommonNodes : no coords specified !");
00208   _coords->findCommonTuples(prec,limitNodeId,comm,commIndex);
00209 }
00210 
00211 std::vector<int> MEDCouplingPointSet::getNodeIdsNearPoint(const double *pos, double eps) const throw(INTERP_KERNEL::Exception)
00212 {
00213   std::vector<int> c,cI;
00214   getNodeIdsNearPoints(pos,1,eps,c,cI);
00215   return c;
00216 }
00217 
00223 void MEDCouplingPointSet::getNodeIdsNearPoints(const double *pos, int nbOfNodes, double eps, std::vector<int>& c, std::vector<int>& cI) const throw(INTERP_KERNEL::Exception)
00224 {
00225   if(!_coords)
00226     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : no coordiantes set !");
00227   const double *coordsPtr=_coords->getConstPointer();
00228   if(!coordsPtr)
00229     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : coordiante array set but no data inside it !");
00230   int spaceDim=getSpaceDimension();
00231   int nbNodes=getNumberOfNodes();
00232   std::vector<double> bbox(2*nbNodes*spaceDim);
00233   for(int i=0;i<nbNodes;i++)
00234     {
00235       for(int j=0;j<spaceDim;j++)
00236         {
00237           bbox[2*spaceDim*i+2*j]=coordsPtr[spaceDim*i+j];
00238           bbox[2*spaceDim*i+2*j+1]=coordsPtr[spaceDim*i+j];
00239         }
00240     }
00241   std::vector<int> ret;
00242   c.clear();
00243   cI.resize(1); cI[0]=0;
00244   switch(spaceDim)
00245     {
00246     case 3:
00247       findNodeIdsNearPointAlg<3>(bbox,pos,nbOfNodes,eps,c,cI);
00248       break;
00249     case 2:
00250       findNodeIdsNearPointAlg<2>(bbox,pos,nbOfNodes,eps,c,cI);
00251       break;
00252     case 1:
00253       findNodeIdsNearPointAlg<1>(bbox,pos,nbOfNodes,eps,c,cI);
00254       break;
00255     default:
00256       throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for getNodeIdsNearPoint. Must be 1, 2 or 3.");
00257     }
00258 }
00259 
00265 DataArrayInt *MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat(const DataArrayInt *comm, const DataArrayInt *commIndex,
00266                                                                           int& newNbOfNodes) const
00267 {
00268   if(!_coords)
00269     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat : no coords specified !");
00270   return DataArrayInt::BuildOld2NewArrayFromSurjectiveFormat2(getNumberOfNodes(),comm,commIndex,newNbOfNodes);
00271 }
00272 
00273 /*
00274  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
00275  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
00276  * This value is asked because often known by the caller of this method.
00277  * @param newNodeNumbers array specifying the new numbering in old2New convention..
00278  * @param newNbOfNodes the new number of nodes.
00279  */
00280 void MEDCouplingPointSet::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
00281 {
00282   if(!_coords)
00283     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::renumberNodes : no coords specified !");
00284   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->renumberAndReduce(newNodeNumbers,newNbOfNodes);
00285   setCoords(newCoords);
00286 }
00287 
00288 /*
00289  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
00290  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
00291  * This value is asked because often known by the caller of this method.
00292  * Contrary to ParaMEDMEM::MEDCouplingPointSet::renumberNodes method for merged nodes the barycenter of them is computed here.
00293  *
00294  * @param newNodeNumbers array specifying the new numbering.
00295  * @param newNbOfNodes the new number of nodes.
00296  */
00297 void MEDCouplingPointSet::renumberNodes2(const int *newNodeNumbers, int newNbOfNodes)
00298 {
00299   DataArrayDouble *newCoords=DataArrayDouble::New();
00300   std::vector<int> div(newNbOfNodes);
00301   int spaceDim=getSpaceDimension();
00302   newCoords->alloc(newNbOfNodes,spaceDim);
00303   newCoords->copyStringInfoFrom(*_coords);
00304   newCoords->fillWithZero();
00305   int oldNbOfNodes=getNumberOfNodes();
00306   double *ptToFill=newCoords->getPointer();
00307   const double *oldCoordsPtr=_coords->getConstPointer();
00308   for(int i=0;i<oldNbOfNodes;i++)
00309     {
00310       std::transform(oldCoordsPtr+i*spaceDim,oldCoordsPtr+(i+1)*spaceDim,ptToFill+newNodeNumbers[i]*spaceDim,
00311                      ptToFill+newNodeNumbers[i]*spaceDim,std::plus<double>());
00312       div[newNodeNumbers[i]]++;
00313     }
00314   for(int i=0;i<newNbOfNodes;i++)
00315     ptToFill=std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)div[i]));
00316   setCoords(newCoords);
00317   newCoords->decrRef();
00318 }
00319 
00325 void MEDCouplingPointSet::getBoundingBox(double *bbox) const throw(INTERP_KERNEL::Exception)
00326 {
00327   if(!_coords)
00328     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getBoundingBox : Coordinates not set !");
00329   _coords->getMinMaxPerComponent(bbox);
00330 }
00331 
00335 void MEDCouplingPointSet::zipCoords()
00336 {
00337   checkFullyDefined();
00338   DataArrayInt *traducer=zipCoordsTraducer();
00339   traducer->decrRef();
00340 }
00341 
00342 struct MEDCouplingCompAbs
00343 {
00344   bool operator()(double x, double y) { return std::abs(x)<std::abs(y);}
00345 };
00346 
00352 double MEDCouplingPointSet::getCaracteristicDimension() const
00353 {
00354   if(!_coords)
00355     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCaracteristicDimension : Coordinates not set !");
00356   const double *coords=_coords->getConstPointer();
00357   int nbOfValues=_coords->getNbOfElems();
00358   return std::abs(*std::max_element(coords,coords+nbOfValues,MEDCouplingCompAbs()));
00359 }
00360 
00369 void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
00370 {
00371   int spaceDim=getSpaceDimension();
00372   if(spaceDim==3)
00373     rotate3D(center,vector,angle);
00374   else if(spaceDim==2)
00375     rotate2D(center,angle);
00376   else
00377     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
00378   _coords->declareAsNew();
00379   updateTime();
00380 }
00381 
00386 void MEDCouplingPointSet::translate(const double *vector)
00387 {
00388   double *coords=_coords->getPointer();
00389   int nbNodes=getNumberOfNodes();
00390   int dim=getSpaceDimension();
00391   for(int i=0; i<nbNodes; i++)
00392     for(int idim=0; idim<dim;idim++)
00393       coords[i*dim+idim]+=vector[idim];
00394   _coords->declareAsNew();
00395   updateTime();
00396 }
00397 
00403 void MEDCouplingPointSet::scale(const double *point, double factor)
00404 {
00405   double *coords=_coords->getPointer();
00406   int nbNodes=getNumberOfNodes();
00407   int dim=getSpaceDimension();
00408   double *tmp=new double[dim];
00409   for(int i=0;i<nbNodes;i++)
00410     {
00411       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
00412       std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
00413       std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
00414     }
00415   delete [] tmp;
00416   _coords->declareAsNew();
00417   updateTime();
00418 }
00419 
00428 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue) throw(INTERP_KERNEL::Exception)
00429 {
00430   if(getCoords()==0)
00431     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
00432   if(newSpaceDim<1)
00433     throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
00434   int oldSpaceDim=getSpaceDimension();
00435   if(newSpaceDim==oldSpaceDim)
00436     return ;
00437   DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
00438   setCoords(newCoords);
00439   newCoords->decrRef();
00440   updateTime();
00441 }
00442 
00447 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
00448 {
00449   if(_coords==other._coords)
00450     return ;
00451   if(!_coords)
00452     throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
00453   if(!other._coords)
00454     throw INTERP_KERNEL::Exception("Other instance has no coords whereas current has !");
00455   if(!_coords->isEqualWithoutConsideringStr(*other._coords,epsilon))
00456     throw INTERP_KERNEL::Exception("Coords are not the same !");
00457   setCoords(other._coords);
00458 }
00459 
00469 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
00470 {
00471   if(getSpaceDimension()!=3)
00472     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
00473   int nbOfNodes=getNumberOfNodes();
00474   double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
00475   double deno=sqrt(a*a+b*b+c*c);
00476   const double *work=_coords->getConstPointer();
00477   for(int i=0;i<nbOfNodes;i++)
00478     {
00479       if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
00480         nodes.push_back(i);
00481       work+=3;
00482     }
00483 }
00484 
00495 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
00496 {
00497   int spaceDim=getSpaceDimension();
00498   if(spaceDim!=2 && spaceDim!=3)
00499     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
00500   int nbOfNodes=getNumberOfNodes();
00501   double den=0.;
00502   for(int i=0;i<spaceDim;i++)
00503     den+=vec[i]*vec[i];
00504   double deno=sqrt(den);
00505   if(deno<10.*eps)
00506     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
00507   INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
00508   for(int i=0;i<spaceDim;i++)
00509     vecn[i]=vec[i]/deno;
00510   const double *work=_coords->getConstPointer();
00511   if(spaceDim==2)
00512     {
00513       for(int i=0;i<nbOfNodes;i++)
00514         {
00515           if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
00516             nodes.push_back(i);
00517           work+=2;
00518         }
00519     }
00520   else
00521     {
00522       for(int i=0;i<nbOfNodes;i++)
00523         {
00524           double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
00525           double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
00526           double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
00527           if(std::sqrt(a*a+b*b+c*c)<eps)
00528             nodes.push_back(i);
00529           work+=3;
00530         }
00531     }
00532 }
00533 
00537 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception)
00538 {
00539   int spaceDim=m1->getSpaceDimension();
00540   if(spaceDim!=m2->getSpaceDimension())
00541     throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
00542   return DataArrayDouble::Aggregate(m1->getCoords(),m2->getCoords());
00543 }
00544 
00545 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception)
00546 {
00547   if(ms.empty())
00548     throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : input array must be NON EMPTY !");
00549   std::vector<const MEDCouplingPointSet *>::const_iterator it=ms.begin();
00550   std::vector<const DataArrayDouble *> coo(ms.size());
00551   int spaceDim=(*it)->getSpaceDimension();
00552   coo[0]=(*it++)->getCoords();
00553   for(int i=1;it!=ms.end();it++,i++)
00554     {
00555       const DataArrayDouble *tmp=(*it)->getCoords();
00556       if(tmp)
00557         {
00558           if((*it)->getSpaceDimension()==spaceDim)
00559             coo[i]=tmp;
00560           else
00561             throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
00562         }
00563       else
00564         throw INTERP_KERNEL::Exception("Empty coords detected during call of MergeNodesArray !");
00565     }
00566   return DataArrayDouble::Aggregate(coo);
00567 }
00568 
00573 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
00574 {
00575   switch(type)
00576     {
00577     case UNSTRUCTURED:
00578       return MEDCouplingUMesh::New();
00579     case UNSTRUCTURED_DESC:
00580       return MEDCouplingUMeshDesc::New();
00581     default:
00582       throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
00583     }
00584 }
00585 
00589 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
00590 {
00591   int it,order;
00592   double time=getTime(it,order);
00593   if(_coords)
00594     {
00595       int spaceDim=getSpaceDimension();
00596       littleStrings.resize(spaceDim+4);
00597       littleStrings[0]=getName();
00598       littleStrings[1]=getDescription();
00599       littleStrings[2]=_coords->getName();
00600       littleStrings[3]=getTimeUnit();
00601       for(int i=0;i<spaceDim;i++)
00602         littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
00603       tinyInfo.clear();
00604       tinyInfo.push_back(getType());
00605       tinyInfo.push_back(spaceDim);
00606       tinyInfo.push_back(getNumberOfNodes());
00607       tinyInfo.push_back(it);
00608       tinyInfo.push_back(order);
00609       tinyInfoD.push_back(time);
00610     }
00611   else
00612     {
00613       littleStrings.resize(3);
00614       littleStrings[0]=getName();
00615       littleStrings[1]=getDescription();
00616       littleStrings[2]=getTimeUnit();
00617       tinyInfo.clear();
00618       tinyInfo.push_back(getType());
00619       tinyInfo.push_back(-1);
00620       tinyInfo.push_back(-1);
00621       tinyInfo.push_back(it);
00622       tinyInfo.push_back(order);
00623       tinyInfoD.push_back(time);
00624     }
00625 }
00626 
00630 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
00631 {
00632   if(_coords)
00633     {
00634       a2=const_cast<DataArrayDouble *>(getCoords());
00635       a2->incrRef();
00636     }
00637   else
00638     a2=0;
00639 }
00640 
00645 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
00646 {
00647   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
00648     {
00649       a2->alloc(tinyInfo[2],tinyInfo[1]);
00650       littleStrings.resize(tinyInfo[1]+4);
00651     }
00652   else
00653     {
00654       littleStrings.resize(3);
00655     }
00656 }
00657 
00662 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
00663 {
00664   if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
00665     {
00666       setCoords(a2);
00667       setName(littleStrings[0].c_str());
00668       setDescription(littleStrings[1].c_str());
00669       a2->setName(littleStrings[2].c_str());
00670       setTimeUnit(littleStrings[3].c_str());
00671       for(int i=0;i<tinyInfo[1];i++)
00672         getCoords()->setInfoOnComponent(i,littleStrings[i+4].c_str());
00673       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
00674     }
00675   else
00676     {
00677       setName(littleStrings[0].c_str());
00678       setDescription(littleStrings[1].c_str());
00679       setTimeUnit(littleStrings[2].c_str());
00680       setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
00681     }
00682 }
00683 
00687 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
00688 {
00689   double* bbtemp = new double[2*dim];
00690   double deltamax=0.0;
00691 
00692   for (int i=0; i< dim; i++)
00693     {
00694       double delta = bb1[2*i+1]-bb1[2*i];
00695       if ( delta > deltamax )
00696         {
00697           deltamax = delta ;
00698         }
00699     }
00700   for (int i=0; i<dim; i++)
00701     {
00702       bbtemp[i*2]=bb1[i*2]-deltamax*eps;
00703       bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
00704     }
00705   
00706   for (int idim=0; idim < dim; idim++)
00707     {
00708       bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
00709         && (bb2[idim*2]<bbtemp[idim*2+1]) ;
00710       if (!intersects)
00711         {
00712           delete [] bbtemp;
00713           return false; 
00714         }
00715     }
00716   delete [] bbtemp;
00717   return true;
00718 }
00719 
00723 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
00724 {
00725   double* bbtemp = new double[2*dim];
00726   double deltamax=0.0;
00727 
00728   for (int i=0; i< dim; i++)
00729     {
00730       double delta = bb2[2*i+1]-bb2[2*i];
00731       if ( delta > deltamax )
00732         {
00733           deltamax = delta ;
00734         }
00735     }
00736   for (int i=0; i<dim; i++)
00737     {
00738       bbtemp[i*2]=bb2[i*2]-deltamax*eps;
00739       bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
00740     }
00741   
00742   bool intersects = !bb1.isDisjointWith( bbtemp );
00743   delete [] bbtemp;
00744   return intersects;
00745 }
00746 
00750 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
00751 {
00752   double *coords=_coords->getPointer();
00753   int nbNodes=getNumberOfNodes();
00754   Rotate3DAlg(center,vect,angle,nbNodes,coords);
00755 }
00756 
00761 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
00762 {
00763   double sina=sin(angle);
00764   double cosa=cos(angle);
00765   double vectorNorm[3];
00766   double matrix[9];
00767   double matrixTmp[9];
00768   double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
00769   std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
00770   //rotation matrix computation
00771   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;
00772   matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
00773   matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
00774   matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
00775   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
00776   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
00777   matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
00778   matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
00779   matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
00780   std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
00781   std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
00782   //rotation matrix computed.
00783   double tmp[3];
00784   for(int i=0; i<nbNodes; i++)
00785     {
00786       std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
00787       coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
00788       coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
00789       coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
00790     }
00791 }
00792 
00799 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
00800 {
00801   return buildPartOfMySelf(start,end,true);
00802 }
00803 
00813 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
00814 {
00815   MEDCouplingPointSet *ret=buildPartOfMySelf(start,end,true);
00816   arr=ret->zipCoordsTraducer();
00817   return ret;
00818 }
00819 
00820 
00824 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
00825 {
00826   double *coords=_coords->getPointer();
00827   int nbNodes=getNumberOfNodes();
00828   Rotate2DAlg(center,angle,nbNodes,coords);
00829 }
00830 
00835 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
00836 {
00837   double cosa=cos(angle);
00838   double sina=sin(angle);
00839   double matrix[4];
00840   matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
00841   double tmp[2];
00842   for(int i=0; i<nbNodes; i++)
00843     {
00844       std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
00845       coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
00846       coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
00847     }
00848 }
00849 
00851 
00852 class DummyClsMCPS
00853 {
00854 public:
00855   static const int MY_SPACEDIM=3;
00856   static const int MY_MESHDIM=2;
00857   typedef int MyConnType;
00858   static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
00859 };
00860 
00862 
00869 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
00870 {
00871   const double *coords=_coords->getConstPointer();
00872   int spaceDim=getSpaceDimension();
00873   for(const int *it=startConn;it!=endConn;it++)
00874     res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
00875   if(spaceDim==2)
00876     return ;
00877   if(spaceDim==3)
00878     {
00879       std::vector<double> cpy(res);
00880       int nbNodes=(int)std::distance(startConn,endConn);
00881       INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0.,0.,true);
00882       res.resize(2*nbNodes);
00883       for(int i=0;i<nbNodes;i++)
00884         {
00885           res[2*i]=cpy[3*i];
00886           res[2*i+1]=cpy[3*i+1];
00887         }
00888       return ;
00889     }
00890   throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
00891 }
00892 
00896 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
00897 {
00898   std::size_t nbOfNodes=res.size()/2;
00899   std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
00900   for(std::size_t i=0;i<nbOfNodes;i++)
00901     {
00902       INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
00903       nodes[i]=tmp;
00904     }
00905   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
00906   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
00907   INTERP_KERNEL::QuadraticPolygon *pol=0;
00908   if(isQuad)
00909     pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
00910   else
00911     pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
00912   bool ret=pol->isButterflyAbs();
00913   delete pol;
00914   return ret;
00915 }