1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
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"
36 using namespace ParaMEDMEM;
38 MEDCouplingPointSet::MEDCouplingPointSet():_coords(0)
42 MEDCouplingPointSet::MEDCouplingPointSet(const MEDCouplingPointSet& other, bool deepCopy):MEDCouplingMesh(other),_coords(0)
45 _coords=other._coords->performCpy(deepCopy);
48 MEDCouplingPointSet::~MEDCouplingPointSet()
54 int MEDCouplingPointSet::getNumberOfNodes() const
57 return _coords->getNumberOfTuples();
59 throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !");
62 int MEDCouplingPointSet::getSpaceDimension() const
65 return _coords->getNumberOfComponents();
67 throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !");
70 void MEDCouplingPointSet::updateTime() const
74 updateTimeWith(*_coords);
78 std::size_t MEDCouplingPointSet::getHeapMemorySize() const
82 ret+=_coords->getHeapMemorySize();
83 return MEDCouplingMesh::getHeapMemorySize()+ret;
86 void MEDCouplingPointSet::setCoords(const DataArrayDouble *coords)
88 if( coords != _coords )
92 _coords=const_cast<DataArrayDouble *>(coords);
99 DataArrayDouble *MEDCouplingPointSet::getCoordinatesAndOwner() const
107 * This method copyies all tiny strings from other (name and components name).
108 * @throw if other and this have not same mesh type.
110 void MEDCouplingPointSet::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
112 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
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);
120 bool MEDCouplingPointSet::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
123 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::isEqualIfNotWhy : null mesh instance in input !");
124 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
127 reason="mesh given in input is not castable in MEDCouplingPointSet !";
130 if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
132 if(!areCoordsEqualIfNotWhy(*otherC,prec,reason))
137 bool MEDCouplingPointSet::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
139 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
142 if(!areCoordsEqualWithoutConsideringStr(*otherC,prec))
147 bool MEDCouplingPointSet::areCoordsEqualIfNotWhy(const MEDCouplingPointSet& other, double prec, std::string& reason) const
149 if(_coords==0 && other._coords==0)
151 if(_coords==0 || other._coords==0)
153 reason="Only one PointSet between the two this and other has coordinate defined !";
156 if(_coords==other._coords)
158 bool ret=_coords->isEqualIfNotWhy(*other._coords,prec,reason);
160 reason.insert(0,"Coordinates DataArray do not match : ");
164 bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const
167 return areCoordsEqualIfNotWhy(other,prec,tmp);
170 bool MEDCouplingPointSet::areCoordsEqualWithoutConsideringStr(const MEDCouplingPointSet& other, double prec) const
172 if(_coords==0 && other._coords==0)
174 if(_coords==0 || other._coords==0)
176 if(_coords==other._coords)
178 return _coords->isEqualWithoutConsideringStr(*other._coords,prec);
182 * Returns coordinates of node with id 'nodeId' and append it in 'coo'.
184 void MEDCouplingPointSet::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
187 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCoordinatesOfNode : no coordinates array set !");
188 int nbNodes=getNumberOfNodes();
189 if(nodeId>=0 && nodeId<nbNodes)
191 const double *cooPtr=_coords->getConstPointer();
192 int spaceDim=getSpaceDimension();
193 coo.insert(coo.end(),cooPtr+spaceDim*nodeId,cooPtr+spaceDim*(nodeId+1));
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());
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.
210 DataArrayInt *MEDCouplingPointSet::buildPermArrayForMergeNode(double precision, int limitNodeId, bool& areNodesMerged, int& newNbOfNodes) const
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);
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)
229 void MEDCouplingPointSet::findCommonNodes(double prec, int limitNodeId, DataArrayInt *&comm, DataArrayInt *&commIndex) const
232 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findCommonNodes : no coords specified !");
233 _coords->findCommonTuples(prec,limitNodeId,comm,commIndex);
236 DataArrayInt *MEDCouplingPointSet::getNodeIdsNearPoint(const double *pos, double eps) const throw(INTERP_KERNEL::Exception)
238 DataArrayInt *c=0,*cI=0;
239 getNodeIdsNearPoints(pos,1,eps,c,cI);
240 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cITmp(cI);
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.
249 void MEDCouplingPointSet::getNodeIdsNearPoints(const double *pos, int nbOfNodes, double eps, DataArrayInt *& c, DataArrayInt *& cI) const throw(INTERP_KERNEL::Exception)
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);
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.
264 DataArrayInt *MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat(const DataArrayInt *comm, const DataArrayInt *commIndex,
265 int& newNbOfNodes) const
268 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat : no coords specified !");
269 return DataArrayInt::BuildOld2NewArrayFromSurjectiveFormat2(getNumberOfNodes(),comm->begin(),commIndex->begin(),commIndex->end(),newNbOfNodes);
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.
279 void MEDCouplingPointSet::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
282 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::renumberNodes : no coords specified !");
283 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->renumberAndReduce(newNodeNumbers,newNbOfNodes);
284 setCoords(newCoords);
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.
293 * @param newNodeNumbers array specifying the new numbering.
294 * @param newNbOfNodes the new number of nodes.
296 void MEDCouplingPointSet::renumberNodes2(const int *newNodeNumbers, int newNbOfNodes)
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++)
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]]++;
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();
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().
324 void MEDCouplingPointSet::getBoundingBox(double *bbox) const throw(INTERP_KERNEL::Exception)
327 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getBoundingBox : Coordinates not set !");
328 _coords->getMinMaxPerComponent(bbox);
332 * This method removes useless nodes in coords.
334 void MEDCouplingPointSet::zipCoords()
337 DataArrayInt *traducer=zipCoordsTraducer();
341 struct MEDCouplingCompAbs
343 bool operator()(double x, double y) { return std::abs(x)<std::abs(y);}
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.
351 double MEDCouplingPointSet::getCaracteristicDimension() const
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()));
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.
364 * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
366 * \warning this method is non const and alterates coordinates in \b this without modifying.
368 void MEDCouplingPointSet::recenterForMaxPrecision(double eps) throw(INTERP_KERNEL::Exception)
371 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::recenterForMaxPrecision : Coordinates not set !");
372 _coords->recenterForMaxPrecision(eps);
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.
384 void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
386 int spaceDim=getSpaceDimension();
388 rotate3D(center,vector,angle);
390 rotate2D(center,angle);
392 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
393 _coords->declareAsNew();
398 * Non const method that operates a translation of 'this'.
399 * @param vector in array of size getSpaceDimension().
401 void MEDCouplingPointSet::translate(const double *vector)
404 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::translate : NULL input vector !");
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();
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
422 void MEDCouplingPointSet::scale(const double *point, double factor)
425 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::scale : NULL input point !");
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++)
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>());
437 _coords->declareAsNew();
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.
449 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue) throw(INTERP_KERNEL::Exception)
452 throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
454 throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
455 int oldSpaceDim=getSpaceDimension();
456 if(newSpaceDim==oldSpaceDim)
458 DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
459 setCoords(newCoords);
460 newCoords->decrRef();
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.
468 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
470 if(_coords==other._coords)
473 throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
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);
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.
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
488 void MEDCouplingPointSet::duplicateNodesInCoords(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd) throw(INTERP_KERNEL::Exception)
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);
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.
506 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
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++)
516 if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
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.
532 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
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();
539 for(int i=0;i<spaceDim;i++)
541 double deno=sqrt(den);
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++)
547 const double *work=_coords->getConstPointer();
550 for(int i=0;i<nbOfNodes;i++)
552 if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
559 for(int i=0;i<nbOfNodes;i++)
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)
572 * merge _coords arrays of m1 and m2 and returns the union. The returned instance is newly created with ref count == 1.
574 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception)
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());
582 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception)
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++)
592 const DataArrayDouble *tmp=(*it)->getCoords();
595 if((*it)->getSpaceDimension()==spaceDim)
598 throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
601 throw INTERP_KERNEL::Exception("Empty coords detected during call of MergeNodesArray !");
603 return DataArrayDouble::Aggregate(coo);
607 * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
608 * This method is used during unserialization process.
610 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
615 return MEDCouplingUMesh::New();
616 case UNSTRUCTURED_DESC:
617 return MEDCouplingUMeshDesc::New();
619 throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
624 * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
626 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
629 double time=getTime(it,order);
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);
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);
650 littleStrings.resize(3);
651 littleStrings[0]=getName();
652 littleStrings[1]=getDescription();
653 littleStrings[2]=getTimeUnit();
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);
665 * Third and final step of serialization process.
667 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
671 a2=const_cast<DataArrayDouble *>(getCoords());
679 * Second step of serialization process.
680 * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
682 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
684 if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
686 a2->alloc(tinyInfo[2],tinyInfo[1]);
687 littleStrings.resize(tinyInfo[1]+4);
691 littleStrings.resize(3);
696 * Second and final unserialization process.
697 * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
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)
701 if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
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]);
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]);
721 void MEDCouplingPointSet::checkCoherency() const throw(INTERP_KERNEL::Exception)
724 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::checkCoherency : no coordinates set !");
728 * Intersect Bounding Box given 2 Bounding Boxes.
730 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
732 double* bbtemp = new double[2*dim];
735 for (int i=0; i< dim; i++)
737 double delta = bb1[2*i+1]-bb1[2*i];
738 if ( delta > deltamax )
743 for (int i=0; i<dim; i++)
745 bbtemp[i*2]=bb1[i*2]-deltamax*eps;
746 bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
749 for (int idim=0; idim < dim; idim++)
751 bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
752 && (bb2[idim*2]<bbtemp[idim*2+1]) ;
764 * Intersect 2 given Bounding Boxes.
766 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
768 double* bbtemp = new double[2*dim];
771 for (int i=0; i< dim; i++)
773 double delta = bb2[2*i+1]-bb2[2*i];
774 if ( delta > deltamax )
779 for (int i=0; i<dim; i++)
781 bbtemp[i*2]=bb2[i*2]-deltamax*eps;
782 bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
785 bool intersects = !bb1.isDisjointWith( bbtemp );
791 * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
793 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
795 double *coords=_coords->getPointer();
796 int nbNodes=getNumberOfNodes();
797 Rotate3DAlg(center,vect,angle,nbNodes,coords);
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'.
804 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
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];
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.
831 for(int i=0; i<nbNodes; i++)
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];
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.
846 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
848 return buildPartOfMySelf(start,end,true);
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.
860 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
862 MEDCouplingPointSet *ret=buildPartOfMySelf(start,end,true);
863 arr=ret->zipCoordsTraducer();
869 * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
871 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
873 double *coords=_coords->getPointer();
874 int nbNodes=getNumberOfNodes();
875 Rotate2DAlg(center,angle,nbNodes,coords);
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'.
882 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
884 double cosa=cos(angle);
885 double sina=sin(angle);
887 matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
889 for(int i=0; i<nbNodes; i++)
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];
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;
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.
916 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
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));
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++)
933 res[2*i+1]=cpy[3*i+1];
937 throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
941 * low level method that checks that the 2D cell is not a butterfly cell.
943 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
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++)
949 INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
952 INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
953 INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
954 INTERP_KERNEL::QuadraticPolygon *pol=0;
956 pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
958 pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
959 bool ret=pol->isButterflyAbs();