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
20 #include "MEDCouplingPointSet.hxx"
21 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
22 #include "MEDCouplingUMesh.hxx"
23 #include "MEDCouplingUMeshDesc.hxx"
24 #include "MEDCouplingMemArray.hxx"
25 #include "PlanarIntersector.txx"
26 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
27 #include "InterpKernelGeo2DNode.hxx"
28 #include "DirectedBoundingBox.hxx"
29 #include "InterpKernelAutoPtr.hxx"
35 using namespace ParaMEDMEM;
37 MEDCouplingPointSet::MEDCouplingPointSet():_coords(0)
41 MEDCouplingPointSet::MEDCouplingPointSet(const MEDCouplingPointSet& other, bool deepCopy):MEDCouplingMesh(other),_coords(0)
44 _coords=other._coords->performCpy(deepCopy);
47 MEDCouplingPointSet::~MEDCouplingPointSet()
53 int MEDCouplingPointSet::getNumberOfNodes() const
56 return _coords->getNumberOfTuples();
58 throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !");
61 int MEDCouplingPointSet::getSpaceDimension() const
64 return _coords->getNumberOfComponents();
66 throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !");
69 void MEDCouplingPointSet::updateTime() const
73 updateTimeWith(*_coords);
77 void MEDCouplingPointSet::setCoords(const DataArrayDouble *coords)
79 if( coords != _coords )
83 _coords=const_cast<DataArrayDouble *>(coords);
90 DataArrayDouble *MEDCouplingPointSet::getCoordinatesAndOwner() const
98 * This method copyies all tiny strings from other (name and components name).
99 * @throw if other and this have not same mesh type.
101 void MEDCouplingPointSet::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
103 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
105 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::copyTinyStringsFrom : meshes have not same type !");
106 MEDCouplingMesh::copyTinyStringsFrom(other);
107 if(_coords && otherC->_coords)
108 _coords->copyStringInfoFrom(*otherC->_coords);
111 bool MEDCouplingPointSet::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
114 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::isEqualIfNotWhy : null mesh instance in input !");
115 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
118 reason="mesh given in input is not castable in MEDCouplingPointSet !";
121 if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
123 if(!areCoordsEqualIfNotWhy(*otherC,prec,reason))
128 bool MEDCouplingPointSet::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
130 const MEDCouplingPointSet *otherC=dynamic_cast<const MEDCouplingPointSet *>(other);
133 if(!areCoordsEqualWithoutConsideringStr(*otherC,prec))
138 bool MEDCouplingPointSet::areCoordsEqualIfNotWhy(const MEDCouplingPointSet& other, double prec, std::string& reason) const
140 if(_coords==0 && other._coords==0)
142 if(_coords==0 || other._coords==0)
144 reason="Only one PointSet between the two this and other has coordinate defined !";
147 if(_coords==other._coords)
149 bool ret=_coords->isEqualIfNotWhy(*other._coords,prec,reason);
151 reason.insert(0,"Coordinates DataArray do not match : ");
155 bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const
158 return areCoordsEqualIfNotWhy(other,prec,tmp);
161 bool MEDCouplingPointSet::areCoordsEqualWithoutConsideringStr(const MEDCouplingPointSet& other, double prec) const
163 if(_coords==0 && other._coords==0)
165 if(_coords==0 || other._coords==0)
167 if(_coords==other._coords)
169 return _coords->isEqualWithoutConsideringStr(*other._coords,prec);
173 * Returns coordinates of node with id 'nodeId' and append it in 'coo'.
175 void MEDCouplingPointSet::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
178 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCoordinatesOfNode : no coordinates array set !");
179 int nbNodes=getNumberOfNodes();
180 if(nodeId>=0 && nodeId<nbNodes)
182 const double *cooPtr=_coords->getConstPointer();
183 int spaceDim=getSpaceDimension();
184 coo.insert(coo.end(),cooPtr+spaceDim*nodeId,cooPtr+spaceDim*(nodeId+1));
188 std::ostringstream oss; oss << "MEDCouplingPointSet::getCoordinatesOfNode : request of nodeId \"" << nodeId << "\" but it should be in [0,"<< nbNodes << ") !";
189 throw INTERP_KERNEL::Exception(oss.str().c_str());
194 * This method is typically the base method used for implementation of mergeNodes. This method computes this permutation array using as input,
195 * This method is const ! So this method simply computes the array, no permutation of nodes is done.
196 * 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.
197 * To desactivate this advanced feature put -1 to this argument.
198 * @param areNodesMerged output parameter that states if some nodes have been "merged" in returned array
199 * @param newNbOfNodes output parameter too this is the maximal id in returned array to avoid to recompute it.
201 DataArrayInt *MEDCouplingPointSet::buildPermArrayForMergeNode(double precision, int limitNodeId, bool& areNodesMerged, int& newNbOfNodes) const
203 DataArrayInt *comm,*commI;
204 findCommonNodes(precision,limitNodeId,comm,commI);
205 int oldNbOfNodes=getNumberOfNodes();
206 DataArrayInt *ret=buildNewNumberingFromCommonNodesFormat(comm,commI,newNbOfNodes);
207 areNodesMerged=(oldNbOfNodes!=newNbOfNodes);
214 * 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
215 * comm and commIndex.
216 * @param limitNodeId is the limit node id. All nodes which id is strictly lower than 'limitNodeId' will not be merged each other.
217 * @param comm out parameter (not inout)
218 * @param commIndex out parameter (not inout)
220 void MEDCouplingPointSet::findCommonNodes(double prec, int limitNodeId, DataArrayInt *&comm, DataArrayInt *&commIndex) const
223 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findCommonNodes : no coords specified !");
224 _coords->findCommonTuples(prec,limitNodeId,comm,commIndex);
227 std::vector<int> MEDCouplingPointSet::getNodeIdsNearPoint(const double *pos, double eps) const throw(INTERP_KERNEL::Exception)
229 std::vector<int> c,cI;
230 getNodeIdsNearPoints(pos,1,eps,c,cI);
235 * Given a point given by its position 'pos' this method finds the set of node ids that are a a distance lower than eps.
236 * Position 'pos' is expected to be of size getSpaceDimension()*nbOfNodes. If not the behabiour is not warranted.
237 * This method throws an exception if no coordiantes are set.
239 void MEDCouplingPointSet::getNodeIdsNearPoints(const double *pos, int nbOfNodes, double eps, std::vector<int>& c, std::vector<int>& cI) const throw(INTERP_KERNEL::Exception)
242 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : no coordiantes set !");
243 int spaceDim=getSpaceDimension();
244 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> points=DataArrayDouble::New();
245 points->useArray(pos,false,CPP_DEALLOC,nbOfNodes,spaceDim);
246 _coords->computeTupleIdsNearTuples(points,eps,c,cI);
250 * @param comm in param in the same format than one returned by findCommonNodes method.
251 * @param commI in param in the same format than one returned by findCommonNodes method.
252 * @return the old to new correspondance array.
254 DataArrayInt *MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat(const DataArrayInt *comm, const DataArrayInt *commIndex,
255 int& newNbOfNodes) const
258 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat : no coords specified !");
259 return DataArrayInt::BuildOld2NewArrayFromSurjectiveFormat2(getNumberOfNodes(),comm,commIndex,newNbOfNodes);
263 * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
264 * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
265 * This value is asked because often known by the caller of this method.
266 * @param newNodeNumbers array specifying the new numbering in old2New convention..
267 * @param newNbOfNodes the new number of nodes.
269 void MEDCouplingPointSet::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
272 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::renumberNodes : no coords specified !");
273 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->renumberAndReduce(newNodeNumbers,newNbOfNodes);
274 setCoords(newCoords);
278 * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
279 * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
280 * This value is asked because often known by the caller of this method.
281 * Contrary to ParaMEDMEM::MEDCouplingPointSet::renumberNodes method for merged nodes the barycenter of them is computed here.
283 * @param newNodeNumbers array specifying the new numbering.
284 * @param newNbOfNodes the new number of nodes.
286 void MEDCouplingPointSet::renumberNodes2(const int *newNodeNumbers, int newNbOfNodes)
288 DataArrayDouble *newCoords=DataArrayDouble::New();
289 std::vector<int> div(newNbOfNodes);
290 int spaceDim=getSpaceDimension();
291 newCoords->alloc(newNbOfNodes,spaceDim);
292 newCoords->copyStringInfoFrom(*_coords);
293 newCoords->fillWithZero();
294 int oldNbOfNodes=getNumberOfNodes();
295 double *ptToFill=newCoords->getPointer();
296 const double *oldCoordsPtr=_coords->getConstPointer();
297 for(int i=0;i<oldNbOfNodes;i++)
299 std::transform(oldCoordsPtr+i*spaceDim,oldCoordsPtr+(i+1)*spaceDim,ptToFill+newNodeNumbers[i]*spaceDim,
300 ptToFill+newNodeNumbers[i]*spaceDim,std::plus<double>());
301 div[newNodeNumbers[i]]++;
303 for(int i=0;i<newNbOfNodes;i++)
304 ptToFill=std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)div[i]));
305 setCoords(newCoords);
306 newCoords->decrRef();
310 * This method fills bbox params like that : bbox[0]=XMin, bbox[1]=XMax, bbox[2]=YMin...
311 * The returned bounding box is arranged along trihedron.
312 * @param bbox out array of size 2*this->getSpaceDimension().
314 void MEDCouplingPointSet::getBoundingBox(double *bbox) const throw(INTERP_KERNEL::Exception)
317 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getBoundingBox : Coordinates not set !");
318 _coords->getMinMaxPerComponent(bbox);
322 * This method removes useless nodes in coords.
324 void MEDCouplingPointSet::zipCoords()
327 DataArrayInt *traducer=zipCoordsTraducer();
331 struct MEDCouplingCompAbs
333 bool operator()(double x, double y) { return std::abs(x)<std::abs(y);}
337 * This method expects that _coords attribute is set.
338 * @return the carateristic dimension of point set. This caracteristic dimension is the max of difference
339 * @exception If _coords attribute not set.
341 double MEDCouplingPointSet::getCaracteristicDimension() const
344 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getCaracteristicDimension : Coordinates not set !");
345 const double *coords=_coords->getConstPointer();
346 int nbOfValues=_coords->getNbOfElems();
347 return std::abs(*std::max_element(coords,coords+nbOfValues,MEDCouplingCompAbs()));
351 * 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
352 * around origin of 'radius' 1.
354 * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
356 * \warning this method is non const and alterates coordinates in \b this without modifying.
358 void MEDCouplingPointSet::recenterForMaxPrecision(double eps) throw(INTERP_KERNEL::Exception)
361 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::recenterForMaxPrecision : Coordinates not set !");
362 _coords->recenterForMaxPrecision(eps);
367 * Non const method that operates a rotation of 'this'.
368 * If spaceDim==2 'vector' parameter is ignored (and could be 0) and the rotation is done around 'center' with angle specified by 'angle'.
369 * If spaceDim==3 the rotation axe is defined by ('center','vector') and the angle is 'angle'.
370 * @param center an array of size getSpaceDimension().
371 * @param vector in array of size getSpaceDimension().
372 * @param angle angle of rotation in radian.
374 void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
376 int spaceDim=getSpaceDimension();
378 rotate3D(center,vector,angle);
380 rotate2D(center,angle);
382 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
383 _coords->declareAsNew();
388 * Non const method that operates a translation of 'this'.
389 * @param vector in array of size getSpaceDimension().
391 void MEDCouplingPointSet::translate(const double *vector)
393 double *coords=_coords->getPointer();
394 int nbNodes=getNumberOfNodes();
395 int dim=getSpaceDimension();
396 for(int i=0; i<nbNodes; i++)
397 for(int idim=0; idim<dim;idim++)
398 coords[i*dim+idim]+=vector[idim];
399 _coords->declareAsNew();
404 * Non const method that operates a scale on 'this' with 'point' as reference point of scale and with factor 'factor'.
405 * @param point in array of size getSpaceDimension().
406 * @param factor factor of the scaling
408 void MEDCouplingPointSet::scale(const double *point, double factor)
410 double *coords=_coords->getPointer();
411 int nbNodes=getNumberOfNodes();
412 int dim=getSpaceDimension();
413 double *tmp=new double[dim];
414 for(int i=0;i<nbNodes;i++)
416 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
417 std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
418 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
421 _coords->declareAsNew();
426 * This method is only available for already defined coordinates.
427 * If not an INTERP_KERNEL::Exception is thrown. The 'newSpaceDim' input must be greater or equal to 1.
428 * This method simply convert this to newSpaceDim space :
429 * - by putting a 0. for each \f$ i^{th} \f$ components of each coord of nodes so that i>=getSpaceDim(), if 'newSpaceDim'>getSpaceDimsion()
430 * - by ignoring each \f$ i^{th} \f$ components of each coord of nodes so that i >= 'newSpaceDim', 'newSpaceDim'<getSpaceDimension()
431 * If newSpaceDim==getSpaceDim() nothing is done by this method.
433 void MEDCouplingPointSet::changeSpaceDimension(int newSpaceDim, double dftValue) throw(INTERP_KERNEL::Exception)
436 throw INTERP_KERNEL::Exception("changeSpaceDimension must be called on an MEDCouplingPointSet instance with coordinates set !");
438 throw INTERP_KERNEL::Exception("changeSpaceDimension must be called a newSpaceDim >=1 !");
439 int oldSpaceDim=getSpaceDimension();
440 if(newSpaceDim==oldSpaceDim)
442 DataArrayDouble *newCoords=getCoords()->changeNbOfComponents(newSpaceDim,dftValue);
443 setCoords(newCoords);
444 newCoords->decrRef();
449 * This method try to substitute this->_coords with other._coords if arrays match.
450 * This method potentially modifies 'this' if it succeeds, otherway an exception is thrown.
452 void MEDCouplingPointSet::tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
454 if(_coords==other._coords)
457 throw INTERP_KERNEL::Exception("Current instance has no coords whereas other has !");
459 throw INTERP_KERNEL::Exception("Other instance has no coords whereas current has !");
460 if(!_coords->isEqualWithoutConsideringStr(*other._coords,epsilon))
461 throw INTERP_KERNEL::Exception("Coords are not the same !");
462 setCoords(other._coords);
466 * This method duplicates the nodes whose ids are in [\b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd) and put the result of their duplication at the end
467 * of existing node ids.
469 * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
470 * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
472 void MEDCouplingPointSet::duplicateNodesInCoords(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd) throw(INTERP_KERNEL::Exception)
475 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::duplicateNodesInCoords : no coords set !");
476 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=_coords->selectByTupleIdSafe(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd);
477 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords2=DataArrayDouble::Aggregate(_coords,newCoords);
478 setCoords(newCoords2);
482 * This method is expecting to be called for meshes so that getSpaceDimension() returns 3.
483 * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from plane
484 * defined by the point 'pt' and the vector 'vec'.
485 * @param pt points to an array of size 3 and represents a point that owns to plane.
486 * @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)
487 * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
488 * @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.
490 void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
492 if(getSpaceDimension()!=3)
493 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
494 int nbOfNodes=getNumberOfNodes();
495 double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
496 double deno=sqrt(a*a+b*b+c*c);
497 const double *work=_coords->getConstPointer();
498 for(int i=0;i<nbOfNodes;i++)
500 if(std::abs(a*work[0]+b*work[1]+c*work[2]+d)/deno<eps)
507 * This method is expecting to be called for meshes so that getSpaceDimension() returns 2 or 3.
508 * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from a line
509 * defined by the point 'pt' and the vector 'vec'. 'pt' and 'vec' are expected to have a dimension equal to space dimension of 'this'
510 * @param pt points to an array of size this->getSpaceDimension and represents a point that owns to plane.
511 * @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.
512 * But norm must be greater than 10*abs(eps)
513 * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
514 * @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.
516 void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
518 int spaceDim=getSpaceDimension();
519 if(spaceDim!=2 && spaceDim!=3)
520 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
521 int nbOfNodes=getNumberOfNodes();
523 for(int i=0;i<spaceDim;i++)
525 double deno=sqrt(den);
527 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
528 INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
529 for(int i=0;i<spaceDim;i++)
531 const double *work=_coords->getConstPointer();
534 for(int i=0;i<nbOfNodes;i++)
536 if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
543 for(int i=0;i<nbOfNodes;i++)
545 double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
546 double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
547 double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
548 if(std::sqrt(a*a+b*b+c*c)<eps)
556 * merge _coords arrays of m1 and m2 and returns the union. The returned instance is newly created with ref count == 1.
558 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception)
560 int spaceDim=m1->getSpaceDimension();
561 if(spaceDim!=m2->getSpaceDimension())
562 throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
563 return DataArrayDouble::Aggregate(m1->getCoords(),m2->getCoords());
566 DataArrayDouble *MEDCouplingPointSet::MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception)
569 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::MergeNodesArray : input array must be NON EMPTY !");
570 std::vector<const MEDCouplingPointSet *>::const_iterator it=ms.begin();
571 std::vector<const DataArrayDouble *> coo(ms.size());
572 int spaceDim=(*it)->getSpaceDimension();
573 coo[0]=(*it++)->getCoords();
574 for(int i=1;it!=ms.end();it++,i++)
576 const DataArrayDouble *tmp=(*it)->getCoords();
579 if((*it)->getSpaceDimension()==spaceDim)
582 throw INTERP_KERNEL::Exception("Mismatch in SpaceDim during call of MergeNodesArray !");
585 throw INTERP_KERNEL::Exception("Empty coords detected during call of MergeNodesArray !");
587 return DataArrayDouble::Aggregate(coo);
591 * Factory to build new instance of instanciable subclasses of MEDCouplingPointSet.
592 * This method is used during unserialization process.
594 MEDCouplingPointSet *MEDCouplingPointSet::BuildInstanceFromMeshType(MEDCouplingMeshType type)
599 return MEDCouplingUMesh::New();
600 case UNSTRUCTURED_DESC:
601 return MEDCouplingUMeshDesc::New();
603 throw INTERP_KERNEL::Exception("Invalid type of mesh specified");
608 * First step of serialization process. Used by ParaMEDMEM and MEDCouplingCorba to transfert data between process.
610 void MEDCouplingPointSet::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
613 double time=getTime(it,order);
616 int spaceDim=getSpaceDimension();
617 littleStrings.resize(spaceDim+4);
618 littleStrings[0]=getName();
619 littleStrings[1]=getDescription();
620 littleStrings[2]=_coords->getName();
621 littleStrings[3]=getTimeUnit();
622 for(int i=0;i<spaceDim;i++)
623 littleStrings[i+4]=getCoords()->getInfoOnComponent(i);
625 tinyInfo.push_back(getType());
626 tinyInfo.push_back(spaceDim);
627 tinyInfo.push_back(getNumberOfNodes());
628 tinyInfo.push_back(it);
629 tinyInfo.push_back(order);
630 tinyInfoD.push_back(time);
634 littleStrings.resize(3);
635 littleStrings[0]=getName();
636 littleStrings[1]=getDescription();
637 littleStrings[2]=getTimeUnit();
639 tinyInfo.push_back(getType());
640 tinyInfo.push_back(-1);
641 tinyInfo.push_back(-1);
642 tinyInfo.push_back(it);
643 tinyInfo.push_back(order);
644 tinyInfoD.push_back(time);
649 * Third and final step of serialization process.
651 void MEDCouplingPointSet::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
655 a2=const_cast<DataArrayDouble *>(getCoords());
663 * Second step of serialization process.
664 * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
666 void MEDCouplingPointSet::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
668 if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
670 a2->alloc(tinyInfo[2],tinyInfo[1]);
671 littleStrings.resize(tinyInfo[1]+4);
675 littleStrings.resize(3);
680 * Second and final unserialization process.
681 * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
683 void MEDCouplingPointSet::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
685 if(tinyInfo[2]>=0 && tinyInfo[1]>=1)
688 setName(littleStrings[0].c_str());
689 setDescription(littleStrings[1].c_str());
690 a2->setName(littleStrings[2].c_str());
691 setTimeUnit(littleStrings[3].c_str());
692 for(int i=0;i<tinyInfo[1];i++)
693 getCoords()->setInfoOnComponent(i,littleStrings[i+4].c_str());
694 setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
698 setName(littleStrings[0].c_str());
699 setDescription(littleStrings[1].c_str());
700 setTimeUnit(littleStrings[2].c_str());
701 setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
706 * Intersect Bounding Box given 2 Bounding Boxes.
708 bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
710 double* bbtemp = new double[2*dim];
713 for (int i=0; i< dim; i++)
715 double delta = bb1[2*i+1]-bb1[2*i];
716 if ( delta > deltamax )
721 for (int i=0; i<dim; i++)
723 bbtemp[i*2]=bb1[i*2]-deltamax*eps;
724 bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
727 for (int idim=0; idim < dim; idim++)
729 bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
730 && (bb2[idim*2]<bbtemp[idim*2+1]) ;
742 * Intersect 2 given Bounding Boxes.
744 bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps)
746 double* bbtemp = new double[2*dim];
749 for (int i=0; i< dim; i++)
751 double delta = bb2[2*i+1]-bb2[2*i];
752 if ( delta > deltamax )
757 for (int i=0; i<dim; i++)
759 bbtemp[i*2]=bb2[i*2]-deltamax*eps;
760 bbtemp[i*2+1]=bb2[i*2+1]+deltamax*eps;
763 bool intersects = !bb1.isDisjointWith( bbtemp );
769 * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
771 void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
773 double *coords=_coords->getPointer();
774 int nbNodes=getNumberOfNodes();
775 Rotate3DAlg(center,vect,angle,nbNodes,coords);
779 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
780 * around an axe ('center','vect') and with angle 'angle'.
782 void MEDCouplingPointSet::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
785 throw INTERP_KERNEL::Exception("MEDCouplingPointSet::Rotate3DAlg : null vector in input !");
786 double sina=sin(angle);
787 double cosa=cos(angle);
788 double vectorNorm[3];
791 double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
792 std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
793 //rotation matrix computation
794 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;
795 matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
796 matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
797 matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
798 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
799 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
800 matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
801 matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
802 matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
803 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
804 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
805 //rotation matrix computed.
807 for(int i=0; i<nbNodes; i++)
809 std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
810 coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
811 coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
812 coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
817 * This method implements pure virtual method MEDCouplingMesh::buildPart.
818 * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end').
819 * The coords are kept unchanged contrary to pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
820 * The returned mesh has to be managed by the caller.
822 MEDCouplingMesh *MEDCouplingPointSet::buildPart(const int *start, const int *end) const
824 return buildPartOfMySelf(start,end,true);
828 * This method implements pure virtual method MEDCouplingMesh::buildPartAndReduceNodes.
829 * This method build a part of 'this' by simply keeping cells whose ids are in ['start','end') \b and potentially reduces the nodes set
830 * behind returned mesh. This cause an overhead but it is lesser in memory.
831 * This method returns an array too. This array allows to the caller to know the mapping between nodeids in 'this' and nodeids in
832 * returned mesh. This is quite usefull for MEDCouplingFieldDouble on nodes for example...
833 * 'arr' is in old2New format of size ret->getNumberOfCells like MEDCouplingUMesh::zipCoordsTraducer is.
834 * The returned mesh has to be managed by the caller.
836 MEDCouplingMesh *MEDCouplingPointSet::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
838 MEDCouplingPointSet *ret=buildPartOfMySelf(start,end,true);
839 arr=ret->zipCoordsTraducer();
845 * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
847 void MEDCouplingPointSet::rotate2D(const double *center, double angle)
849 double *coords=_coords->getPointer();
850 int nbNodes=getNumberOfNodes();
851 Rotate2DAlg(center,angle,nbNodes,coords);
855 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in 'coords'
856 * around the center point 'center' and with angle 'angle'.
858 void MEDCouplingPointSet::Rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
860 double cosa=cos(angle);
861 double sina=sin(angle);
863 matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
865 for(int i=0; i<nbNodes; i++)
867 std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
868 coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
869 coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
878 static const int MY_SPACEDIM=3;
879 static const int MY_MESHDIM=2;
880 typedef int MyConnType;
881 static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
887 * res should be an empty vector before calling this method.
888 * This method returns all the node coordinates included in _coords which ids are in [startConn;endConn) and put it into 'res' vector.
889 * If spaceDim==3 a projection will be done for each nodes on the middle plane containing these all nodes in [startConn;endConn).
890 * And after each projected nodes are moved to Oxy plane in order to consider these nodes as 2D nodes.
892 void MEDCouplingPointSet::project2DCellOnXY(const int *startConn, const int *endConn, std::vector<double>& res) const
894 const double *coords=_coords->getConstPointer();
895 int spaceDim=getSpaceDimension();
896 for(const int *it=startConn;it!=endConn;it++)
897 res.insert(res.end(),coords+spaceDim*(*it),coords+spaceDim*(*it+1));
902 std::vector<double> cpy(res);
903 int nbNodes=(int)std::distance(startConn,endConn);
904 INTERP_KERNEL::PlanarIntersector<DummyClsMCPS,int>::projection(&res[0],&cpy[0],nbNodes,nbNodes,1.e-12,0.,0.,true);
905 res.resize(2*nbNodes);
906 for(int i=0;i<nbNodes;i++)
909 res[2*i+1]=cpy[3*i+1];
913 throw INTERP_KERNEL::Exception("Invalid spacedim for project2DCellOnXY !");
917 * low level method that checks that the 2D cell is not a butterfly cell.
919 bool MEDCouplingPointSet::isButterfly2DCell(const std::vector<double>& res, bool isQuad, double eps)
921 std::size_t nbOfNodes=res.size()/2;
922 std::vector<INTERP_KERNEL::Node *> nodes(nbOfNodes);
923 for(std::size_t i=0;i<nbOfNodes;i++)
925 INTERP_KERNEL::Node *tmp=new INTERP_KERNEL::Node(res[2*i],res[2*i+1]);
928 INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
929 INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
930 INTERP_KERNEL::QuadraticPolygon *pol=0;
932 pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
934 pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
935 bool ret=pol->isButterflyAbs();