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 void MEDCouplingPointSet::setCoords(const DataArrayDouble *coords)
80 if( coords != _coords )
84 _coords=const_cast<DataArrayDouble *>(coords);
91 DataArrayDouble *MEDCouplingPointSet::getCoordinatesAndOwner() const
99 * This method copyies all tiny strings from other (name and components name).
100 * @throw if other and this have not same mesh type.
102 void MEDCouplingPointSet::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
104 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
106 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::copyTinyStringsFrom : meshes have not same type !");
107 MEDCouplingMesh::copyTinyStringsFrom(other);
108 if(_coords && otherC->_coords)
109 _coords->copyStringInfoFrom(*otherC->_coords);
112 bool MEDCouplingPointSet::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
115 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::isEqualIfNotWhy : null mesh instance in input !");
116 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
119 reason="mesh given in input is not castable in MEDCouplingPointSet !";
122 if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
124 if(!areCoordsEqualIfNotWhy(*otherC,prec,reason))
129 bool MEDCouplingPointSet::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
131 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
134 if(!areCoordsEqualWithoutConsideringStr(*otherC,prec))
139 bool MEDCouplingPointSet::areCoordsEqualIfNotWhy(const MEDCouplingPointSet& other, double prec, std::string& reason) const
141 if(_coords==0 && other._coords==0)
143 if(_coords==0 || other._coords==0)
145 reason="Only one PointSet between the two this and other has coordinate defined !";
148 if(_coords==other._coords)
150 bool ret=_coords->isEqualIfNotWhy(*other._coords,prec,reason);
152 reason.insert(0,"Coordinates DataArray do not match : ");
156 bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const
159 return areCoordsEqualIfNotWhy(other,prec,tmp);
162 bool MEDCouplingPointSet::areCoordsEqualWithoutConsideringStr(const MEDCouplingPointSet& other, double prec) const
164 if(_coords==0 && other._coords==0)
166 if(_coords==0 || other._coords==0)
168 if(_coords==other._coords)
170 return _coords->isEqualWithoutConsideringStr(*other._coords,prec);
174 * Returns coordinates of node with id 'nodeId' and append it in 'coo'.
176 void MEDCouplingPointSet::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
179 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCoordinatesOfNode : no coordinates array set !");
180 int nbNodes=getNumberOfNodes();
181 if(nodeId>=0 && nodeId<nbNodes)
183 const double *cooPtr=_coords->getConstPointer();
184 int spaceDim=getSpaceDimension();
185 coo.insert(coo.end(),cooPtr+spaceDim*nodeId,cooPtr+spaceDim*(nodeId+1));
189 std::ostringstream oss; oss << "MEDCouplingPointSet::getCoordinatesOfNode : request of nodeId \"" << nodeId << "\" but it should be in [0,"<< nbNodes << ") !";
190 throw INTERP_KERNEL::Exception(oss.str().c_str());
195 * This method is typically the base method used for implementation of mergeNodes. This method computes this permutation array using as input,
196 * This method is const ! So this method simply computes the array, no permutation of nodes is done.
197 * 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.
198 * To desactivate this advanced feature put -1 to this argument.
199 * @param areNodesMerged output parameter that states if some nodes have been "merged" in returned array
200 * @param newNbOfNodes output parameter too this is the maximal id in returned array to avoid to recompute it.
202 DataArrayInt *MEDCouplingPointSet::buildPermArrayForMergeNode(double precision, int limitNodeId, bool& areNodesMerged, int& newNbOfNodes) const
204 DataArrayInt *comm,*commI;
205 findCommonNodes(precision,limitNodeId,comm,commI);
206 int oldNbOfNodes=getNumberOfNodes();
207 DataArrayInt *ret=buildNewNumberingFromCommonNodesFormat(comm,commI,newNbOfNodes);
208 areNodesMerged=(oldNbOfNodes!=newNbOfNodes);
215 * 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
216 * comm and commIndex.
217 * @param limitNodeId is the limit node id. All nodes which id is strictly lower than 'limitNodeId' will not be merged each other.
218 * @param comm out parameter (not inout)
219 * @param commIndex out parameter (not inout)
221 void MEDCouplingPointSet::findCommonNodes(double prec, int limitNodeId, DataArrayInt *&comm, DataArrayInt *&commIndex) const
224 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findCommonNodes : no coords specified !");
225 _coords->findCommonTuples(prec,limitNodeId,comm,commIndex);
228 std::vector<int> MEDCouplingPointSet::getNodeIdsNearPoint(const double *pos, double eps) const throw(INTERP_KERNEL::Exception)
230 std::vector<int> c,cI;
231 getNodeIdsNearPoints(pos,1,eps,c,cI);
236 * Given a point given by its position 'pos' this method finds the set of node ids that are a a distance lower than eps.
237 * Position 'pos' is expected to be of size getSpaceDimension()*nbOfNodes. If not the behabiour is not warranted.
238 * This method throws an exception if no coordiantes are set.
240 void MEDCouplingPointSet::getNodeIdsNearPoints(const double *pos, int nbOfNodes, double eps, std::vector<int>& c, std::vector<int>& cI) const throw(INTERP_KERNEL::Exception)
243 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : no coordiantes set !");
244 int spaceDim=getSpaceDimension();
245 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> points=DataArrayDouble::New();
246 points->useArray(pos,false,CPP_DEALLOC,nbOfNodes,spaceDim);
247 _coords->computeTupleIdsNearTuples(points,eps,c,cI);
251 * @param comm in param in the same format than one returned by findCommonNodes method.
252 * @param commI in param in the same format than one returned by findCommonNodes method.
253 * @return the old to new correspondance array.
255 DataArrayInt *MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat(const DataArrayInt *comm, const DataArrayInt *commIndex,
256 int& newNbOfNodes) const
259 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat : no coords specified !");
260 return DataArrayInt::BuildOld2NewArrayFromSurjectiveFormat2(getNumberOfNodes(),comm,commIndex,newNbOfNodes);
264 * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
265 * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
266 * This value is asked because often known by the caller of this method.
267 * @param newNodeNumbers array specifying the new numbering in old2New convention..
268 * @param newNbOfNodes the new number of nodes.
270 void MEDCouplingPointSet::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
273 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::renumberNodes : no coords specified !");
274 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->renumberAndReduce(newNodeNumbers,newNbOfNodes);
275 setCoords(newCoords);
279 * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
280 * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
281 * This value is asked because often known by the caller of this method.
282 * Contrary to ParaMEDMEM::MEDCouplingPointSet::renumberNodes method for merged nodes the barycenter of them is computed here.
284 * @param newNodeNumbers array specifying the new numbering.
285 * @param newNbOfNodes the new number of nodes.
287 void MEDCouplingPointSet::renumberNodes2(const int *newNodeNumbers, int newNbOfNodes)
289 DataArrayDouble *newCoords=DataArrayDouble::New();
290 std::vector<int> div(newNbOfNodes);
291 int spaceDim=getSpaceDimension();
292 newCoords->alloc(newNbOfNodes,spaceDim);
293 newCoords->copyStringInfoFrom(*_coords);
294 newCoords->fillWithZero();
295 int oldNbOfNodes=getNumberOfNodes();
296 double *ptToFill=newCoords->getPointer();
297 const double *oldCoordsPtr=_coords->getConstPointer();
298 for(int i=0;i<oldNbOfNodes;i++)
300 std::transform(oldCoordsPtr+i*spaceDim,oldCoordsPtr+(i+1)*spaceDim,ptToFill+newNodeNumbers[i]*spaceDim,
301 ptToFill+newNodeNumbers[i]*spaceDim,std::plus<double>());
302 div[newNodeNumbers[i]]++;
304 for(int i=0;i<newNbOfNodes;i++)
305 ptToFill=std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)div[i]));
306 setCoords(newCoords);
307 newCoords->decrRef();
311 * This method fills bbox params like that : bbox[0]=XMin, bbox[1]=XMax, bbox[2]=YMin...
312 * The returned bounding box is arranged along trihedron.
313 * @param bbox out array of size 2*this->getSpaceDimension().
315 void MEDCouplingPointSet::getBoundingBox(double *bbox) const throw(INTERP_KERNEL::Exception)
318 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getBoundingBox : Coordinates not set !");
319 _coords->getMinMaxPerComponent(bbox);
323 * This method removes useless nodes in coords.
325 void MEDCouplingPointSet::zipCoords()
328 DataArrayInt *traducer=zipCoordsTraducer();
332 struct MEDCouplingCompAbs
334 bool operator()(double x, double y) { return std::abs(x)<std::abs(y);}
338 * This method expects that _coords attribute is set.
339 * @return the carateristic dimension of point set. This caracteristic dimension is the max of difference
340 * @exception If _coords attribute not set.
342 double MEDCouplingPointSet::getCaracteristicDimension() const
345 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCaracteristicDimension : Coordinates not set !");
346 const double *coords=_coords->getConstPointer();
347 int nbOfValues=_coords->getNbOfElems();
348 return std::abs(*std::max_element(coords,coords+nbOfValues,MEDCouplingCompAbs()));
352 * 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
353 * around origin of 'radius' 1.
355 * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
357 * \warning this method is non const and alterates coordinates in \b this without modifying.
359 void MEDCouplingPointSet::recenterForMaxPrecision(double eps) throw(INTERP_KERNEL::Exception)
362 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::recenterForMaxPrecision : Coordinates not set !");
363 _coords->recenterForMaxPrecision(eps);
368 * Non const method that operates a rotation of 'this'.
369 * If spaceDim==2 'vector' parameter is ignored (and could be 0) and the rotation is done around 'center' with angle specified by 'angle'.
370 * If spaceDim==3 the rotation axe is defined by ('center','vector') and the angle is 'angle'.
371 * @param center an array of size getSpaceDimension().
372 * @param vector in array of size getSpaceDimension().
373 * @param angle angle of rotation in radian.
375 void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
377 int spaceDim=getSpaceDimension();
379 rotate3D(center,vector,angle);
381 rotate2D(center,angle);
383 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
384 _coords->declareAsNew();
389 * Non const method that operates a translation of 'this'.
390 * @param vector in array of size getSpaceDimension().
392 void MEDCouplingPointSet::translate(const double *vector)
394 double *coords=_coords->getPointer();
395 int nbNodes=getNumberOfNodes();
396 int dim=getSpaceDimension();
397 for(int i=0; i<nbNodes; i++)
398 for(int idim=0; idim<dim;idim++)
399 coords[i*dim+idim]+=vector[idim];
400 _coords->declareAsNew();
405 * Non const method that operates a scale on 'this' with 'point' as reference point of scale and with factor 'factor'.
406 * @param point in array of size getSpaceDimension().
407 * @param factor factor of the scaling
409 void MEDCouplingPointSet::scale(const double *point, double factor)
411 double *coords=_coords->getPointer();
412 int nbNodes=getNumberOfNodes();
413 int dim=getSpaceDimension();
414 double *tmp=new double[dim];
415 for(int i=0;i<nbNodes;i++)
417 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
418 std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
419 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
422 _coords->declareAsNew();
427 * This method is only available for already defined coordinates.
428 * If not an INTERP_KERNEL::Exception is thrown. The 'newSpaceDim' input must be greater or equal to 1.
429 * This method simply convert this to newSpaceDim space :
430 * - by putting a 0. for each \f$ i^{th} \f$ components of each coord of nodes so that i>=getSpaceDim(), if 'newSpaceDim'>getSpaceDimsion()
431 * - by ignoring each \f$ i^{th} \f$ components of each coord of nodes so that i >= 'newSpaceDim', 'newSpaceDim'<getSpaceDimension()
432 * If newSpaceDim==getSpaceDim() nothing is done by this method.
434 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue) throw(INTERP_KERNEL::Exception)
437 throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
439 throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
440 int oldSpaceDim=getSpaceDimension();
441 if(newSpaceDim==oldSpaceDim)
443 DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
444 setCoords(newCoords);
445 newCoords->decrRef();
450 * This method try to substitute this->_coords with other._coords if arrays match.
451 * This method potentially modifies 'this' if it succeeds, otherway an exception is thrown.
453 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
455 if(_coords==other._coords)
458 throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
460 throw INTERP_KERNEL::Exception("Other instance has no coords whereas current has !");
461 if(!_coords->isEqualWithoutConsideringStr(*other._coords,epsilon))
462 throw INTERP_KERNEL::Exception("Coords are not the same !");
463 setCoords(other._coords);
467 * This method duplicates the nodes whose ids are in [\b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd) and put the result of their duplication at the end
468 * of existing node ids.
470 * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
471 * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
473 void MEDCouplingPointSet::duplicateNodesInCoords(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd) throw(INTERP_KERNEL::Exception)
476 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::duplicateNodesInCoords : no coords set !");
477 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->selectByTupleIdSafe(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd);
478 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords2=DataArrayDouble::Aggregate(_coords,newCoords);
479 setCoords(newCoords2);
483 * This method is expecting to be called for meshes so that getSpaceDimension() returns 3.
484 * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from plane
485 * defined by the point 'pt' and the vector 'vec'.
486 * @param pt points to an array of size 3 and represents a point that owns to plane.
487 * @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)
488 * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
489 * @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.
491 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
493 if(getSpaceDimension()!=3)
494 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
495 int nbOfNodes=getNumberOfNodes();
496 double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
497 double deno=sqrt(a*a+b*b+c*c);
498 const double *work=_coords->getConstPointer();
499 for(int i=0;i<nbOfNodes;i++)
501 if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
508 * This method is expecting to be called for meshes so that getSpaceDimension() returns 2 or 3.
509 * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from a line
510 * defined by the point 'pt' and the vector 'vec'. 'pt' and 'vec' are expected to have a dimension equal to space dimension of 'this'
511 * @param pt points to an array of size this->getSpaceDimension and represents a point that owns to plane.
512 * @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.
513 * But norm must be greater than 10*abs(eps)
514 * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
515 * @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.
517 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
519 int spaceDim=getSpaceDimension();
520 if(spaceDim!=2 && spaceDim!=3)
521 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
522 int nbOfNodes=getNumberOfNodes();
524 for(int i=0;i<spaceDim;i++)
526 double deno=sqrt(den);
528 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
529 INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
530 for(int i=0;i<spaceDim;i++)
532 const double *work=_coords->getConstPointer();
535 for(int i=0;i<nbOfNodes;i++)
537 if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
544 for(int i=0;i<nbOfNodes;i++)
546 double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
547 double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
548 double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
549 if(std::sqrt(a*a+b*b+c*c)<eps)
557 * merge _coords arrays of m1 and m2 and returns the union. The returned instance is newly created with ref count == 1.
559 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception)
561 int spaceDim=m1->getSpaceDimension();
562 if(spaceDim!=m2->getSpaceDimension())
563 throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
564 return DataArrayDouble::Aggregate(m1->getCoords(),m2->getCoords());
567 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception)
570 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : input array must be NON EMPTY !");
571 std::vector<const MEDCouplingPointSet *>::const_iterator it=ms.begin();
572 std::vector<const DataArrayDouble *> coo(ms.size());
573 int spaceDim=(*it)->getSpaceDimension();
574 coo[0]=(*it++)->getCoords();
575 for(int i=1;it!=ms.end();it++,i++)
577 const DataArrayDouble *tmp=(*it)->getCoords();
580 if((*it)->getSpaceDimension()==spaceDim)
583 throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
586 throw INTERP_KERNEL::Exception("Empty coords detected during call of MergeNodesArray !");
588 return DataArrayDouble::Aggregate(coo);
592 * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
593 * This method is used during unserialization process.
595 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
600 return MEDCouplingUMesh::New();
601 case UNSTRUCTURED_DESC:
602 return MEDCouplingUMeshDesc::New();
604 throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
609 * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
611 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
614 double time=getTime(it,order);
617 int spaceDim=getSpaceDimension();
618 littleStrings.resize(spaceDim+4);
619 littleStrings[0]=getName();
620 littleStrings[1]=getDescription();
621 littleStrings[2]=_coords->getName();
622 littleStrings[3]=getTimeUnit();
623 for(int i=0;i<spaceDim;i++)
624 littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
626 tinyInfo.push_back(getType());
627 tinyInfo.push_back(spaceDim);
628 tinyInfo.push_back(getNumberOfNodes());
629 tinyInfo.push_back(it);
630 tinyInfo.push_back(order);
631 tinyInfoD.push_back(time);
635 littleStrings.resize(3);
636 littleStrings[0]=getName();
637 littleStrings[1]=getDescription();
638 littleStrings[2]=getTimeUnit();
640 tinyInfo.push_back(getType());
641 tinyInfo.push_back(-1);
642 tinyInfo.push_back(-1);
643 tinyInfo.push_back(it);
644 tinyInfo.push_back(order);
645 tinyInfoD.push_back(time);
650 * Third and final step of serialization process.
652 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
656 a2=const_cast<DataArrayDouble *>(getCoords());
664 * Second step of serialization process.
665 * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
667 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
669 if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
671 a2->alloc(tinyInfo[2],tinyInfo[1]);
672 littleStrings.resize(tinyInfo[1]+4);
676 littleStrings.resize(3);
681 * Second and final unserialization process.
682 * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
684 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
686 if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
689 setName(littleStrings[0].c_str());
690 setDescription(littleStrings[1].c_str());
691 a2->setName(littleStrings[2].c_str());
692 setTimeUnit(littleStrings[3].c_str());
693 for(int i=0;i<tinyInfo[1];i++)
694 getCoords()->setInfoOnComponent(i,littleStrings[i+4].c_str());
695 setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
699 setName(littleStrings[0].c_str());
700 setDescription(littleStrings[1].c_str());
701 setTimeUnit(littleStrings[2].c_str());
702 setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
707 * Intersect Bounding Box given 2 Bounding Boxes.
709 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
711 double* bbtemp = new double[2*dim];
714 for (int i=0; i< dim; i++)
716 double delta = bb1[2*i+1]-bb1[2*i];
717 if ( delta > deltamax )
722 for (int i=0; i<dim; i++)
724 bbtemp[i*2]=bb1[i*2]-deltamax*eps;
725 bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
728 for (int idim=0; idim < dim; idim++)
730 bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
731 && (bb2[idim*2]<bbtemp[idim*2+1]) ;
743 * Intersect 2 given Bounding Boxes.
745 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
747 double* bbtemp = new double[2*dim];
750 for (int i=0; i< dim; i++)
752 double delta = bb2[2*i+1]-bb2[2*i];
753 if ( delta > deltamax )
758 for (int i=0; i<dim; i++)
760 bbtemp[i*2]=bb2[i*2]-deltamax*eps;
761 bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
764 bool intersects = !bb1.isDisjointWith( bbtemp );
770 * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
772 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
774 double *coords=_coords->getPointer();
775 int nbNodes=getNumberOfNodes();
776 Rotate3DAlg(center,vect,angle,nbNodes,coords);
780 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
781 * around an axe ('center','vect') and with angle 'angle'.
783 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
786 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : null vector in input !");
787 double sina=sin(angle);
788 double cosa=cos(angle);
789 double vectorNorm[3];
792 double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
793 std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
794 //rotation matrix computation
795 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;
796 matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
797 matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
798 matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
799 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
800 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
801 matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
802 matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
803 matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
804 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
805 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
806 //rotation matrix computed.
808 for(int i=0; i<nbNodes; i++)
810 std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
811 coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
812 coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
813 coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
818 * This method implements pure virtual method MEDCouplingMesh::buildPart.
819 * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end').
820 * The coords are kept unchanged contrary to pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
821 * The returned mesh has to be managed by the caller.
823 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
825 return buildPartOfMySelf(start,end,true);
829 * This method implements pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
830 * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end') \b and potentially reduces the nodes set
831 * behind returned mesh. This cause an overhead but it is lesser in memory.
832 * This method returns an array too. This array allows to the caller to know the mapping between nodeids in 'this' and nodeids in
833 * returned mesh. This is quite usefull for MEDCouplingFieldDouble on nodes for example...
834 * 'arr' is in old2New format of size ret->getNumberOfCells like MEDCouplingUMesh::zipCoordsTraducer is.
835 * The returned mesh has to be managed by the caller.
837 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
839 MEDCouplingPointSet *ret=buildPartOfMySelf(start,end,true);
840 arr=ret->zipCoordsTraducer();
846 * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
848 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
850 double *coords=_coords->getPointer();
851 int nbNodes=getNumberOfNodes();
852 Rotate2DAlg(center,angle,nbNodes,coords);
856 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
857 * around the center point 'center' and with angle 'angle'.
859 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
861 double cosa=cos(angle);
862 double sina=sin(angle);
864 matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
866 for(int i=0; i<nbNodes; i++)
868 std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
869 coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
870 coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
879 static const int MY_SPACEDIM=3;
880 static const int MY_MESHDIM=2;
881 typedef int MyConnType;
882 static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
888 * res should be an empty vector before calling this method.
889 * This method returns all the node coordinates included in _coords which ids are in [startConn;endConn) and put it into 'res' vector.
890 * If spaceDim==3 a projection will be done for each nodes on the middle plane containing these all nodes in [startConn;endConn).
891 * And after each projected nodes are moved to Oxy plane in order to consider these nodes as 2D nodes.
893 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
895 const double *coords=_coords->getConstPointer();
896 int spaceDim=getSpaceDimension();
897 for(const int *it=startConn;it!=endConn;it++)
898 res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
903 std::vector<double> cpy(res);
904 int nbNodes=(int)std::distance(startConn,endConn);
905 INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0.,0.,true);
906 res.resize(2*nbNodes);
907 for(int i=0;i<nbNodes;i++)
910 res[2*i+1]=cpy[3*i+1];
914 throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
918 * low level method that checks that the 2D cell is not a butterfly cell.
920 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
922 std::size_t nbOfNodes=res.size()/2;
923 std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
924 for(std::size_t i=0;i<nbOfNodes;i++)
926 INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
929 INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
930 INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
931 INTERP_KERNEL::QuadraticPolygon *pol=0;
933 pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
935 pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
936 bool ret=pol->isButterflyAbs();