+ return ret.retn();
+}
+
+/*!
+ * This methods returns the minimal distance between the two set of points \a this and \a other.
+ * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
+ * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
+ *
+ * \param [out] thisTupleId the tuple id in \a this corresponding to the returned minimal distance
+ * \param [out] otherTupleId the tuple id in \a other corresponding to the returned minimal distance
+ * \return the minimal distance between the two set of points \a this and \a other.
+ * \sa DataArrayDouble::findClosestTupleId
+ */
+double DataArrayDouble::minimalDistanceTo(const DataArrayDouble *other, int& thisTupleId, int& otherTupleId) const
+{
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> part1=findClosestTupleId(other);
+ int nbOfCompo(getNumberOfComponents());
+ int otherNbTuples(other->getNumberOfTuples());
+ const double *thisPt(begin()),*otherPt(other->begin());
+ const int *part1Pt(part1->begin());
+ double ret=std::numeric_limits<double>::max();
+ for(int i=0;i<otherNbTuples;i++,part1Pt++,otherPt+=nbOfCompo)
+ {
+ double tmp(0.);
+ for(int j=0;j<nbOfCompo;j++)
+ tmp+=(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j])*(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j]);
+ if(tmp<ret)
+ { ret=tmp; thisTupleId=*part1Pt; otherTupleId=i; }
+ }
+ return sqrt(ret);
+}
+
+/*!
+ * This methods returns for each tuple in \a other which tuple in \a this is the closest.
+ * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
+ * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
+ *
+ * \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components.
+ * \sa DataArrayDouble::minimalDistanceTo
+ */
+DataArrayInt *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const
+{
+ if(!other)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::findClosestTupleId : other instance is NULL !");
+ checkAllocated(); other->checkAllocated();
+ int nbOfCompo=getNumberOfComponents();
+ if(nbOfCompo!=other->getNumberOfComponents())
+ {
+ std::ostringstream oss; oss << "DataArrayDouble::findClosestTupleId : number of components in this is " << nbOfCompo;
+ oss << ", whereas number of components in other is " << other->getNumberOfComponents() << "! Should be equal !";
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
+ int nbOfTuples=other->getNumberOfTuples();
+ int thisNbOfTuples=getNumberOfTuples();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(nbOfTuples,1);
+ double bounds[6];
+ getMinMaxPerComponent(bounds);
+ switch(nbOfCompo)
+ {
+ case 3:
+ {
+ double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4]));
+ double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta);
+ double characSize=pow((delta*delta*delta)/((double)thisNbOfTuples),1./3.);
+ BBTreePts<3,int> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
+ FindClosestTupleIdAlg<3>(myTree,3.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
+ break;
+ }
+ case 2:
+ {
+ double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2]));
+ double delta=std::max(xDelta,yDelta);
+ double characSize=sqrt(delta/(double)thisNbOfTuples);
+ BBTreePts<2,int> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
+ FindClosestTupleIdAlg<2>(myTree,2.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
+ break;
+ }
+ case 1:
+ {
+ double characSize=fabs(bounds[1]-bounds[0])/thisNbOfTuples;
+ BBTreePts<1,int> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
+ FindClosestTupleIdAlg<1>(myTree,1.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for findClosestTupleId. Must be 1, 2 or 3.");
+ }
+ return ret.retn();
+}
+
+/*!
+ * This method expects that \a this and \a otherBBoxFrmt arrays are bounding box arrays ( as the output of MEDCouplingPointSet::getBoundingBoxForBBTree method ).
+ * This method will return a DataArrayInt array having the same number of tuples than \a this. This returned array tells for each cell in \a this
+ * how many bounding boxes in \a otherBBoxFrmt.
+ * So, this method expects that \a this and \a otherBBoxFrmt have the same number of components.
+ *
+ * \param [in] otherBBoxFrmt - It is an array .
+ * \param [in] eps - the absolute precision of the detection. when eps < 0 the bboxes are enlarged so more interactions are detected. Inversely when > 0 the bboxes are stretched.
+ * \sa MEDCouplingPointSet::getBoundingBoxForBBTree
+ * \throw If \a this and \a otherBBoxFrmt have not the same number of components.
+ * \throw If \a this and \a otherBBoxFrmt number of components is not even (BBox format).
+ */
+DataArrayInt *DataArrayDouble::computeNbOfInteractionsWith(const DataArrayDouble *otherBBoxFrmt, double eps) const
+{
+ if(!otherBBoxFrmt)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : input array is NULL !");
+ if(!isAllocated() || !otherBBoxFrmt->isAllocated())
+ throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : this and input array must be allocated !");
+ int nbOfComp(getNumberOfComponents()),nbOfTuples(getNumberOfTuples());
+ if(nbOfComp!=otherBBoxFrmt->getNumberOfComponents())
+ {
+ std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : this number of components (" << nbOfComp << ") must be equal to the number of components of input array (" << otherBBoxFrmt->getNumberOfComponents() << ") !";
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
+ if(nbOfComp%2!=0)
+ {
+ std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : Number of components (" << nbOfComp << ") is not even ! It should be to be compatible with bbox format !";
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(nbOfTuples,1);
+ const double *thisBBPtr(begin());
+ int *retPtr(ret->getPointer());
+ switch(nbOfComp/2)
+ {
+ case 3:
+ {
+ BBTree<3,int> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
+ for(int i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
+ *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
+ break;
+ }
+ case 2:
+ {
+ BBTree<2,int> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
+ for(int i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
+ *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
+ break;
+ }
+ case 1:
+ {
+ BBTree<1,int> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
+ for(int i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
+ *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : space dimension supported are [1,2,3] !");
+ }
+
+ return ret.retn();