From eac4a49c5dc4e9ffe046b60e798e7cbe46d1bf89 Mon Sep 17 00:00:00 2001 From: eap Date: Thu, 22 Apr 2010 09:28:29 +0000 Subject: [PATCH] EHPOC * DirectedBoundingBox has been created --- src/INTERP_KERNEL/DirectedBoundingBox.cxx | 740 ++++++++++++++++++++++ src/INTERP_KERNEL/DirectedBoundingBox.hxx | 118 ++++ src/INTERP_KERNEL/Makefile.am | 6 +- src/INTERP_KERNELTest/BBTreeTest.cxx | 222 +++++++ src/INTERP_KERNELTest/BBTreeTest.hxx | 6 + src/MEDCoupling/MEDCouplingPointSet.cxx | 28 + src/MEDCoupling/MEDCouplingPointSet.hxx | 7 + src/MEDCoupling/MEDCouplingUMesh.cxx | 53 ++ src/MEDCoupling/MEDCouplingUMesh.hxx | 1 + src/MEDCoupling/MEDCouplingUMeshDesc.cxx | 45 ++ src/MEDCoupling/MEDCouplingUMeshDesc.hxx | 1 + src/ParaMEDMEM/ElementLocator.cxx | 38 +- 12 files changed, 1259 insertions(+), 6 deletions(-) create mode 100644 src/INTERP_KERNEL/DirectedBoundingBox.cxx create mode 100644 src/INTERP_KERNEL/DirectedBoundingBox.hxx diff --git a/src/INTERP_KERNEL/DirectedBoundingBox.cxx b/src/INTERP_KERNEL/DirectedBoundingBox.cxx new file mode 100644 index 000000000..4d0ac2bff --- /dev/null +++ b/src/INTERP_KERNEL/DirectedBoundingBox.cxx @@ -0,0 +1,740 @@ +// Copyright (C) 2009-2010 OPEN CASCADE +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : DirectedBoundingBox.cxx +// Created : Mon Apr 12 14:41:22 2010 +// Author : Edward AGAPOV (eap) + + +#include "DirectedBoundingBox.hxx" + +#include "InterpolationUtils.hxx" + +#define __TENSOR(i,j) tensor[(i)*_dim+(j)] +#define __AXIS(i) (&_axes[(i)*_dim]) +#define __MIN(i) _minmax[i*2] +#define __MAX(i) _minmax[i*2+1] +#define __MYID (long(this)%10000) +#define __DMP(msg) \ + cout << msg << endl + +using namespace std; + +namespace +{ + //================================================================================ + /*! + * \brief Add point coordinates to inertia tensor in 3D space + */ + //================================================================================ + + inline void addPointToInertiaTensor3D(const double* coord, + const double* gc, + vector& tensor) + { + // we fill the upper triangle of tensor only + const int _dim = 3; + double x = coord[0] - gc[0], y = coord[1] - gc[1], z = coord[2] - gc[2]; + __TENSOR(0,0) += y*y + z*z; + __TENSOR(1,1) += x*x + z*z; + __TENSOR(2,2) += x*x + y*y; + __TENSOR(0,1) -= x*y; + __TENSOR(0,2) -= x*z; + __TENSOR(1,2) -= y*z; + } + //================================================================================ + /*! + * \brief Add point coordinates to inertia tensor in 2D space + */ + //================================================================================ + + inline void addPointToInertiaTensor2D(const double* coord, + const double* gc, + vector& tensor) + { + // we fill the upper triangle of tensor only + const int _dim = 2; + double x = coord[0] - gc[0], y = coord[1] - gc[1]; + __TENSOR(0,0) += y*y; + __TENSOR(1,1) += x*x; + __TENSOR(0,1) -= x*y; + } + + //================================================================================ + /*! + * \brief Find eigenvectors of tensor using Jacobi's method + */ + //================================================================================ + + bool JacobiEigenvectorsSearch( const int _dim, vector& tensor, vector& _axes) + { + if ( _dim == 3 ) + __DMP( "Tensor : {" + << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << ", "<<__TENSOR(0,2) << "} " + << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << ", "<<__TENSOR(1,2) << "} " + << "{ "<<__TENSOR(2,0) << ", "<<__TENSOR(2,1) << ", "<<__TENSOR(2,2) << "}} "); + else + __DMP( "Tensor : {" + << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << "} " + << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << "}} "); + + const int maxRot = 5*_dim*_dim; // limit on number of rotations + const double tol = 1e-9; + + // set _axes to identity + int i,j; + for ( i = 0; i < _dim; ++i ) + for ( j = 0; j < _dim; ++j ) + __AXIS(i)[j] = ( i==j ? 1. : 0 ); + + bool solved = false; + for ( int iRot = 0; iRot < maxRot; ++ iRot ) + { + // find max off-diagonal element of the tensor + int k = 0, l = 0; + double max = 0.; + for ( i = 0; i < _dim-1; ++i ) + for ( j = i+1; j < _dim; ++j ) + if ( fabs( __TENSOR(i,j)) > max ) + max = fabs( __TENSOR(i,j) ), k = i, l = j; + solved = ( max < tol ); + if ( solved ) + break; + + // Rotate to make __TENSOR(k,l) == 0 + + double diff = __TENSOR(l,l) - __TENSOR(k,k); + double t; // tangent of rotation angle + if ( fabs(__TENSOR(k,l)) < abs(diff)*1.0e-36) + { + t = __TENSOR(k,l)/diff; + } + else + { + double phi = diff/(2.0*__TENSOR(k,l)); + t = 1.0/(abs(phi) + sqrt(phi*phi + 1.0)); + if ( phi < 0.0) t = -t; + } + double c = 1.0/sqrt(t*t + 1.0); // cosine of rotation angle + double s = t*c; // sine of rotation angle + double tau = s/(1.0 + c); + __TENSOR(k,k) -= t*__TENSOR(k,l); + __TENSOR(l,l) += t*__TENSOR(k,l); + __TENSOR(k,l) = 0.0; + +#define __ROTATE(T,r1,c1,r2,c2) \ +{ \ +int i1 = r1*_dim+c1, i2 = r2*_dim+c2; \ +double t1 = T[i1], t2 = T[i2]; \ +T[i1] -= s * ( t2 + tau * t1);\ +T[i2] += s * ( t1 - tau * t2);\ +} + for ( i = 0; i < k; ++i ) // Case of i < k + __ROTATE(tensor, i,k,i,l); + + for ( i = k+1; i < l; ++i ) // Case of k < i < l + __ROTATE(tensor, k,i,i,l); + + for ( i = l + 1; i < _dim; ++i ) // Case of i > l + __ROTATE(tensor, k,i,l,i); + + for ( i = 0; i < _dim; ++i ) // Update transformation matrix + __ROTATE(_axes, i,k,i,l); + } + + __DMP( "Solved = " << solved ); + if ( _dim == 3 ) { + __DMP( " Eigen " << __TENSOR(0,0)<<", "<<__TENSOR(1,1)<<", "<<__TENSOR(2,2) ); + for ( int i=0; i <3; ++i ) + __DMP( i << ": " << __AXIS(i)[0] << ", " << __AXIS(i)[1] << ", " << __AXIS(i)[2] ); + } + else { + __DMP( " Eigen " << __TENSOR(0,0) << ", " << __TENSOR(1,1) ); + for ( int i=0; i <2; ++i ) + __DMP( i << ": " << __AXIS(i)[0] << ", " << __AXIS(i)[1] ); + } + + return solved; + } + + //================================================================================ + /*! + * \brief Return true if two minmaxes do not intersect + */ + //================================================================================ + + inline bool isMinMaxOut(const double* minmax1, + const double* minmax2, + int dim) + { + for ( int i = 0; i < dim; ++i ) + { + if ( minmax1[i*2] > minmax2[i*2+1] || + minmax1[i*2+1] < minmax2[i*2] ) + return true; + } + return false; + } + +} // noname namespace + +namespace INTERP_KERNEL +{ + + //================================================================================ + /*! + * \brief Creates empty box intended to further initalization via setData() + */ + //================================================================================ + + DirectedBoundingBox::DirectedBoundingBox():_dim(-1) + { + } + + //================================================================================ + /*! + * \brief Creates bounding box of a mesh + * \param pts - coordinates of points in full interlace + * \param numPts - number of points in the mesh + * \param dim - space dimension + */ + //================================================================================ + + DirectedBoundingBox::DirectedBoundingBox(const double* pts, + const unsigned numPts, + const unsigned dim) + : _dim(dim), _axes(dim*dim), _minmax(2*dim) + { + // init box extremities + for ( unsigned i = 0; i < _dim; ++i ) + _minmax[1+i*2] = -numeric_limits::max(), + _minmax[i*2] = numeric_limits::max(); + + if ( numPts < 1 ) return; + + __DMP( "DirectedBoundingBox " << __MYID ); + + const double* coord = pts; + const double* coordEnd = coord + numPts * dim; + + // compute gravity center of points + double gc[3] = {0,0,0}; + if ( dim > 1 ) + { + for ( coord = pts; coord < coordEnd; ) + for ( int i = 0; i < dim; ++i ) + gc[i] += *coord++; + for ( int j = 0; j < dim; ++j ) + gc[j] /= numPts; + + } + + // compute axes and box extremities + vector tensor( dim * dim, 0.); + switch ( dim ) + { + case 3: + for ( coord = pts; coord < coordEnd; coord += dim ) + addPointToInertiaTensor3D( coord, gc, tensor ); + + //computeAxes3D(tensor); + JacobiEigenvectorsSearch(_dim, tensor, _axes); + + for ( coord = pts; coord < coordEnd; coord += dim ) + addPointToBox( coord ); + + break; + + case 2: + for ( coord = pts; coord < coordEnd; coord += dim ) + addPointToInertiaTensor2D( coord, gc, tensor ); + + //computeAxes2D(tensor); + JacobiEigenvectorsSearch(_dim, tensor, _axes); + + for ( coord = pts; coord < coordEnd; coord += dim ) + addPointToBox( coord ); + + break; + + default: + for ( coord = pts; coord < coordEnd; coord += dim ) + { + if ( *coord < _minmax[0] ) _minmax[0] = *coord; + if ( *coord > _minmax[1] ) _minmax[1] = *coord; + } + } + } + + //================================================================================ + /*! + * \brief Creates bounding box of an element + * \param pts - coordinates of points of element + * \param numPts - number of points in the element + * \param dim - space dimension + */ + //================================================================================ + + DirectedBoundingBox::DirectedBoundingBox(const double** pts, + const unsigned numPts, + const unsigned dim) + : _dim(dim), _axes(dim*dim), _minmax(2*dim) + { + // init box extremities + for ( unsigned i = 0; i < _dim; ++i ) + _minmax[1+i*2] = -numeric_limits::max(), + _minmax[i*2] = numeric_limits::max(); + + if ( numPts < 1 ) return; + + __DMP( "DirectedBoundingBox " << __MYID ); + + // compute gravity center of points + double gc[3] = {0,0,0}; + if ( dim > 1 ) + { + for ( unsigned i = 0; i < numPts; ++i ) + for ( int j = 0; j < dim; ++j ) + gc[j] += pts[i][j]; + for ( int j = 0; j < dim; ++j ) + gc[j] /= numPts; + } + + // compute axes and box extremities + vector tensor( dim * dim, 0.); + switch ( dim ) + { + case 3: + for ( unsigned i = 0; i < numPts; ++i ) + addPointToInertiaTensor3D( pts[i], gc, tensor ); + + //computeAxes3D(tensor); + JacobiEigenvectorsSearch(_dim, tensor, _axes); + + for ( unsigned i = 0; i < numPts; ++i ) + addPointToBox( pts[i] ); + + break; + case 2: + for ( unsigned i = 0; i < numPts; ++i ) + addPointToInertiaTensor2D( pts[i], gc, tensor ); + + //computeAxes2D(tensor); + JacobiEigenvectorsSearch(_dim, tensor, _axes); + + for ( unsigned i = 0; i < numPts; ++i ) + addPointToBox( pts[i] ); + + break; + default: + for ( unsigned i = 0; i < numPts; ++i ) + { + if ( pts[i][0] < _minmax[0] ) _minmax[0] = pts[i][0]; + if ( pts[i][0] > _minmax[1] ) _minmax[1] = pts[i][0]; + } + } + } + + //================================================================================ + /*! + * \brief Compute eigenvectors of inertia tensor + */ + //================================================================================ + + // void DirectedBoundingBox::computeAxes3D(const std::vector& tensor) +// { +// // compute principal moments of inertia which are eigenvalues of the tensor +// double eig[3]; +// { +// // coefficients of polynomial equation det(tensor-eig*I) = 0 +// double a = -1; +// double b = __TENSOR(0,0)+__TENSOR(1,1)+__TENSOR(2,2); +// double c = +// __TENSOR(0,1)*__TENSOR(0,1) + +// __TENSOR(0,2)*__TENSOR(0,2) + +// __TENSOR(1,2)*__TENSOR(1,2) - +// __TENSOR(0,0)*__TENSOR(1,1) - +// __TENSOR(0,0)*__TENSOR(2,2) - +// __TENSOR(1,1)*__TENSOR(2,2); +// double d = +// __TENSOR(0,0)*__TENSOR(1,1)*__TENSOR(2,2) - +// __TENSOR(0,0)*__TENSOR(1,2)*__TENSOR(1,2) - +// __TENSOR(1,1)*__TENSOR(0,2)*__TENSOR(0,2) - +// __TENSOR(2,2)*__TENSOR(0,1)*__TENSOR(0,1) + +// __TENSOR(0,1)*__TENSOR(0,2)*__TENSOR(1,2)*2; + +// // find eigenvalues which are roots of characteristic polynomial +// double x = (3*c/a - b*b/(a*a))/3; +// double y = (2*b*b*b/(a*a*a) - 9*b*c/(a*a) + 27*d/a)/27; +// double z = y*y/4 + x*x*x/27; + +// double i = sqrt(y*y/4 - z) + 1e-300; +// double j = -pow(i,1/3.); +// double y2 = -y/(2*i); +// if ( y2 > 1.0) y2 = 1.; else if ( y2 < -1.0) y2 = -1.; +// double k = acos(y2); +// double m = cos(k/3); +// double n = sqrt(3)*sin(k/3); +// double p = -b/(3*a); + +// eig[0] = -2*j*m + p; +// eig[1] = j *(m + n) + p; +// eig[2] = j *(m - n) + p; +// } +// // compute eigenvector of the tensor at each eigenvalue +// // by solving system [tensor-eig*I]*[axis] = 0 +// bool ok = true; +// __DMP( "Tensor : {" +// << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << ", "<<__TENSOR(0,2) << "} " +// << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << ", "<<__TENSOR(1,2) << "} " +// << "{ "<<__TENSOR(2,0) << ", "<<__TENSOR(2,1) << ", "<<__TENSOR(2,2) << "}} "); +// for ( int i = 0; i < 3 && ok; ++i ) // loop on 3 eigenvalues +// { +// // [tensor-eig*I] +// double T[3][3]= +// {{ __TENSOR(0,0)-eig[i],__TENSOR(0,1), __TENSOR(0,2), }, +// { __TENSOR(0,1), __TENSOR(1,1)-eig[i],__TENSOR(1,2), }, +// { __TENSOR(0,2), __TENSOR(1,2), __TENSOR(2,2)-eig[i]}}; +// // The determinant of T is zero, so that the equations are not linearly independent. +// // Therefore, we assign an arbitrary value (1.) to i-th component of eigenvector +// // and use two of the equations to compute the other two components +// double M[2][3], sol[2]; +// for ( int j = 0, c = 0; j < 3; ++j ) +// if ( i == j ) +// M[0][2] = -T[0][j], M[1][2] = -T[1][j]; +// else +// M[0][c] = T[0][j], M[1][c] = T[1][j], c++; + +// ok = solveSystemOfEquations<2>( M, sol ); + +// double* eigenVec = __AXIS(i); +// for ( int j = 0, c = 0; j < 3; ++j ) +// eigenVec[j] = ( i == j ) ? 1. : sol[c++]; + +// // normilize +// double size = sqrt(eigenVec[0]*eigenVec[0] + +// eigenVec[1]*eigenVec[1] + +// eigenVec[2]*eigenVec[2] ); +// if ((ok = (size > numeric_limits::min() ))) +// { +// eigenVec[0] /= size; +// eigenVec[1] /= size; +// eigenVec[2] /= size; +// } +// } +// if ( !ok ) +// { +// __DMP( " solve3EquationSystem() - KO " ); +// _axes = vector( _dim*_dim, 0); +// __AXIS(0)[0] = __AXIS(1)[1] = __AXIS(2)[2] = 1.; +// } +// __DMP( " Eigen " << eig[0] << ", " << eig[1] << ", " << eig[2] ); +// for ( int i=0; i <3; ++i ) +// __DMP( i << ": " << __AXIS(i)[0] << ", " << __AXIS(i)[1] << ", " << __AXIS(i)[2] ); + +// double* a0 = __AXIS(0), *a1 = __AXIS(1); +// double cross[3] = { a0[1]*a1[2]-a1[1]*a0[2], +// a0[2]*a1[0]-a1[2]*a0[0], +// a0[0]*a1[1]-a1[0]*a0[1] }; +// __DMP( " Cross a1^a2 " << cross[0] << ", " << cross[1] << ", " << cross[2] ); +// } + + //================================================================================ + /*! + * \brief Compute eigenvectors of inertia tensor + */ + //================================================================================ + + // void DirectedBoundingBox::computeAxes2D(const std::vector& tensor) +// { +// // compute principal moments of inertia which are eigenvalues of the tensor +// // by solving square equation det(tensor-eig*I) +// double X = (__TENSOR(0,0)+__TENSOR(1,1))/2; +// double Y = sqrt(4*__TENSOR(0,1)*__TENSOR(0,1) + +// (__TENSOR(0,0)-__TENSOR(1,1)) * (__TENSOR(0,0)-__TENSOR(1,1)))/2; +// double eig[2] = +// { +// X + Y, +// X - Y +// }; +// // compute eigenvector of the tensor at each eigenvalue +// // by solving system [tensor-eig*I]*[axis] = 0 +// bool ok = true; +// for ( int i = 0; i < 2 && ok; ++i ) +// { +// // [tensor-eig*I] +// double T[2][2]= +// {{ __TENSOR(0,0)-eig[i],__TENSOR(0,1) }, +// { __TENSOR(0,1), __TENSOR(1,1)-eig[i] }}; + +// // The determinant of T is zero, so that the equations are not linearly independent. +// // Therefore, we assign an arbitrary value (1.) to i-th component of eigenvector +// // and use one equation to compute the other component +// double* eigenVec = __AXIS(i); +// eigenVec[i] = 1.; +// int j = 1-i; +// if ((ok = ( fabs( T[j][j] ) > numeric_limits::min() ))) +// eigenVec[j] = -T[j][i] / T[j][j]; +// } +// if ( !ok ) +// { +// _axes = vector( _dim*_dim, 0); +// __AXIS(0)[0] = __AXIS(1)[1] = 1.; +// } +// } + + //================================================================================ + /*! + * \brief Convert point coordinates into local coordinate system of the box + */ + //================================================================================ + + void DirectedBoundingBox::toLocalCS(const double* p, double* pLoc) const + { + switch ( _dim ) + { + case 3: + pLoc[0] = dotprod<3>( p, __AXIS(0)); + pLoc[1] = dotprod<3>( p, __AXIS(1)); + pLoc[2] = dotprod<3>( p, __AXIS(2)); + break; + case 2: + pLoc[0] = dotprod<2>( p, __AXIS(0)); + pLoc[1] = dotprod<2>( p, __AXIS(1)); + break; + default: + pLoc[0] = p[0]; + } + } + + //================================================================================ + /*! + * \brief Convert point coordinates from local coordinate system of the box to global CS + */ + //================================================================================ + + void DirectedBoundingBox::fromLocalCS(const double* p, double* pGlob) const + { + switch ( _dim ) + { + case 3: + pGlob[0] = p[0] * __AXIS(0)[0] + p[1] * __AXIS(1)[0] + p[2] * __AXIS(2)[0]; + pGlob[1] = p[0] * __AXIS(0)[1] + p[1] * __AXIS(1)[1] + p[2] * __AXIS(2)[1]; + pGlob[2] = p[0] * __AXIS(0)[2] + p[1] * __AXIS(1)[2] + p[2] * __AXIS(2)[2]; + break; + case 2: + pGlob[0] = p[0] * __AXIS(0)[0] + p[1] * __AXIS(1)[0]; + pGlob[1] = p[0] * __AXIS(0)[1] + p[1] * __AXIS(1)[1]; + break; + default: + pGlob[0] = p[0]; + } + } + + //================================================================================ + /*! + * \brief Enlarge box size by given value + */ + //================================================================================ + + void DirectedBoundingBox::enlarge(const double tol) + { + for ( unsigned i = 0; i < _dim; ++i ) + __MIN(i) -= tol, __MAX(i) += tol; + } + + //================================================================================ + /*! + * \brief Return coordinates of corners of bounding box + */ + //================================================================================ + + void DirectedBoundingBox::getCorners(std::vector& corners, + const double* minmax) const + { + int iC, nbCorners = 1; + for ( int i=0;i<_dim;++i ) nbCorners *= 2; + corners.resize( nbCorners * _dim ); + // each coordinate is filled with either min or max, nbSwap is number of corners + // after which min and max swap + int nbSwap = nbCorners/2; + for ( unsigned i = 0; i < _dim; ++i ) + { + iC = 0; + while ( iC < nbCorners ) + { + for (int j = 0; j < nbSwap; ++j, ++iC ) corners[iC*_dim+i] = minmax[i*2]; + for (int j = 0; j < nbSwap; ++j, ++iC ) corners[iC*_dim+i] = minmax[i*2+1]; + } + nbSwap /= 2; + } + } + + //================================================================================ + /*! + * \brief Test if this box intersects with the other + * \retval bool - true if there is no intersection + */ + //================================================================================ + + bool DirectedBoundingBox::isDisjointWith(const DirectedBoundingBox& box) const + { + // boxes are disjoined if their minmaxes in local CS of either of boxes do not intersect + for ( int isThisCS = 0; isThisCS < 2; ++isThisCS ) + { + const DirectedBoundingBox* axisBox = isThisCS ? this : &box; + const DirectedBoundingBox* cornerBox = isThisCS ? &box : this; + + // find minmax of cornerBox in the CS of axisBox + + DirectedBoundingBox mmBox((double*)0,0,_dim); //!< empty box with CS == axisBox->_axes + mmBox._axes = axisBox->_axes; + + vector corners; + getCorners( corners, &cornerBox->_minmax[0] ); + + double globCorner[3]; + for ( int iC = 0, nC = corners.size()/_dim; iC < nC; ++iC) + { + cornerBox->fromLocalCS( &corners[iC*_dim], globCorner ); + mmBox.addPointToBox( globCorner ); + } + if ( isMinMaxOut( &mmBox._minmax[0], &axisBox->_minmax[0], _dim )) + return true; + } + return false; + } + + //================================================================================ + /*! + * \brief Test if this box intersects with an non-directed box + * \retval bool - true if there is no intersection + */ + //================================================================================ + + bool DirectedBoundingBox::isDisjointWith(const double* box) const + { + // boxes are disjoined if their minmaxes in local CS of either of boxes do not intersect + + // compare minmaxes in locals CS of this directed box + { + vector cornersOther; + getCorners( cornersOther, box ); + DirectedBoundingBox mmBox((double*)0,0,_dim); //!< empty box with CS == this->_axes + mmBox._axes = this->_axes; + for ( int iC = 0, nC = cornersOther.size()/_dim; iC < nC; ++iC) + mmBox.addPointToBox( &cornersOther[iC*_dim] ); + + if ( isMinMaxOut( &mmBox._minmax[0], &this->_minmax[0], _dim )) + return true; + } + + // compare minmaxes in global CS + { + vector cornersThis; + getCorners( cornersThis, &_minmax[0] ); + DirectedBoundingBox mmBox((double*)0,0,_dim); //!< initailized _minmax + double globCorner[3]; + for ( int iC = 0, nC = cornersThis.size()/_dim; iC < nC; ++iC) + { + fromLocalCS( &cornersThis[iC*_dim], globCorner ); + for ( int i = 0; i < _dim; ++i ) + { + if ( globCorner[i] < mmBox._minmax[i*2] ) mmBox._minmax[i*2] = globCorner[i]; + if ( globCorner[i] > mmBox._minmax[i*2+1] ) mmBox._minmax[i*2+1] = globCorner[i]; + } + } + if ( isMinMaxOut( &mmBox._minmax[0], box, _dim )) + return true; + } + return false; + } + + //================================================================================ + /*! + * \brief Return true if given point is out of this box + */ + //================================================================================ + + bool DirectedBoundingBox::isOut(const double* point) const + { + double pLoc[3]; + toLocalCS( point, pLoc ); + bool out = isLocalOut( pLoc ); +#ifdef _DEBUG_ + switch (_dim) + { + case 3: + __DMP(__MYID<<": "< DirectedBoundingBox::getData() const + { + vector data(1, _dim); + data.insert( data.end(), &_axes[0], &_axes[0] + _axes.size()); + data.insert( data.end(), &_minmax[0], &_minmax[0] + _minmax.size()); + if ( data.size() < dataSize( _dim )) + data.resize( dataSize( _dim ), 0 ); + return data; + } + + //================================================================================ + /*! + * \brief Initializes self with data retrieved via getData() + */ + //================================================================================ + + void DirectedBoundingBox::setData(const double* data) + { + _dim = unsigned( *data++ ); + if ( _dim > 0 ) + { + _axes.assign( data, data+_dim*_dim ); data += _dim*_dim; + _minmax.assign( data, data+2*_dim ); + } + else + { + _axes.clear(); + _minmax.clear(); + } + } + + //================================================================================ + /*! + * \brief Return size of internal data returned by getData() depending on space dim + */ + //================================================================================ + + int DirectedBoundingBox::dataSize(int dim) + { + return 1 + dim*dim + 2*dim; // : _dim + _axes + _minmax + } +} diff --git a/src/INTERP_KERNEL/DirectedBoundingBox.hxx b/src/INTERP_KERNEL/DirectedBoundingBox.hxx new file mode 100644 index 000000000..93874acdb --- /dev/null +++ b/src/INTERP_KERNEL/DirectedBoundingBox.hxx @@ -0,0 +1,118 @@ +// Copyright (C) 2009-2010 OPEN CASCADE +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef __DirectedBoundingBox_HXX__ +#define __DirectedBoundingBox_HXX__ + +#include "INTERPKERNELDefines.hxx" + +#include + +namespace INTERP_KERNEL +{ + + /** + * \brief Class representing the bounding box of a number of points + * with box axes parallel to principal axes of inertia of points + */ + class INTERPKERNEL_EXPORT DirectedBoundingBox + { + public: + + DirectedBoundingBox(); + + DirectedBoundingBox(const double* pts, const unsigned numPts, const unsigned dim); + + DirectedBoundingBox(const double** pts, const unsigned numPts, const unsigned dim); + + //~DirectedBoundingBox(); + + void enlarge(const double tol); + + bool isDisjointWith(const DirectedBoundingBox& box) const; + + bool isDisjointWith(const double* box) const; + + bool isOut(const double* point) const; + + + // return internal data + std::vector getData() const; + + // initialize with data returned by getData() + void setData(const double* data); + + // return size of internal data + static int dataSize(int dim); + + private: + + //void computeAxes3D(const std::vector& tensor); + + //void computeAxes2D(const std::vector& tensor); + + inline void addPointToBox(const double* coord); + + void toLocalCS(const double* p, double* pLoc) const; + + void fromLocalCS(const double* p, double* pGlob) const; + + inline bool isLocalOut(const double* pLoc) const; + + void getCorners(std::vector& corners, const double* minmax) const; + + unsigned _dim; + + std::vector _axes; //!< principal axes of inertia in full interlace + std::vector _minmax; //!< pairs of min an max coordinates along the axes + + }; + + //================================================================================ + /*! + * \brief Test point in local CS against box extremities + * + */ + //================================================================================ + + inline bool DirectedBoundingBox::isLocalOut(const double* pLoc) const + { + for ( int i = 0; i < _dim; ++i ) + if ( pLoc[i] < _minmax[i*2] || pLoc[i] > _minmax[i*2+1] ) + return true; + return false; + } + + //================================================================================ + /*! + * \brief Update box extremities + */ + //================================================================================ + + inline void DirectedBoundingBox::addPointToBox(const double* coord) + { + for ( int i = 0; i < _dim; ++i ) + { + double c = 0; + for ( int j = 0; j < _dim; ++j ) c += coord[j]*_axes[i*_dim+j]; + if ( c < _minmax[i*2] ) _minmax[i*2] = c; + if ( c > _minmax[i*2+1] ) _minmax[i*2+1] = c; + } + } +} +#endif diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index e1fcbdd7d..0cf9b6692 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -155,7 +155,8 @@ Interpolation1D.txx \ Interpolation2DCurve.hxx \ Interpolation2DCurve.txx \ InterpolationCurve.hxx \ -InterpolationCurve.txx +InterpolationCurve.txx \ +DirectedBoundingBox.hxx EXTRA_DIST += \ @@ -180,7 +181,8 @@ dist_libinterpkernel_la_SOURCES = \ TetraAffineTransform.cxx\ CellModel.cxx\ UnitTetraIntersectionBary.cxx \ - InterpolationOptions.cxx + InterpolationOptions.cxx \ + DirectedBoundingBox.cxx libinterpkernel_la_CPPFLAGS=-I$(srcdir)/Geometric2D -I$(srcdir)/Bases diff --git a/src/INTERP_KERNELTest/BBTreeTest.cxx b/src/INTERP_KERNELTest/BBTreeTest.cxx index ccff041b0..11a5a78d2 100644 --- a/src/INTERP_KERNELTest/BBTreeTest.cxx +++ b/src/INTERP_KERNELTest/BBTreeTest.cxx @@ -20,6 +20,8 @@ #include "BBTreeTest.hxx" #include #include +#include "DirectedBoundingBox.hxx" + namespace INTERP_TEST { @@ -81,5 +83,225 @@ namespace INTERP_TEST delete[] bbox; } + void BBTreeTest::test_DirectedBB_3D() + { + // a rectangle 1x2 extruded along vector (10,0,10) + const int nbP = 8, dim = 3; + double coords[nbP*dim] = + { + 0,0,0, 2,0,0, 2,1,0, 0,1,0, + 10,0,10, 12,0,10, 12,1,10, 10,1,10 + }; + INTERP_KERNEL::DirectedBoundingBox bb( coords, nbP, dim); + bb.enlarge( 1e-12 ); + + // corners of extrusion are IN + for ( int i = 0; i < nbP*dim; i+=dim ) + CPPUNIT_ASSERT( !bb.isOut( coords + i )); + + // points near corners of extrusion are OUT + double p[nbP*dim] = + { + 0,0,3, 6,0,3, 5,1,2, 0,1,2, + 8,0,9, 11,0,8, 11,0.5,8, 8,0.5,9 + }; + for ( int i = 0; i < nbP*dim; i+=dim ) + CPPUNIT_ASSERT( bb.isOut( p + i )); + + // the extrusions shifted by 3 in XOY plane are OUT + double shifted_X[nbP*dim] = + { + 3,0,0, 5,0,0, 5,1,0, 3,1,0, + 13,0,10, 15,0,10, 15,1,10, 13,1,10 + }; + double shifted_x[nbP*dim] = + { + -3,0,0, -1,0,0, -1,1,0, -3,1,0, + 7,0,10, 9,0,10, 9,1,10, 7,1,10 + }; + double shifted_Y[nbP*dim] = + { + 0,3,0, 2,3,0, 2,4,0, 0,4,0, + 10,3,10, 12,3,10, 12,4,10, 10,4,10 + }; + double shifted_y[nbP*dim] = + { + 0,-3,0, 2,-3,0, 2,-2,0, 0,-2,0, + 10,-3,10, 12,-3,10, 12,-2,10, 10,-2,10 + }; + INTERP_KERNEL::DirectedBoundingBox shiftedBB_x( shifted_x, nbP, dim); + INTERP_KERNEL::DirectedBoundingBox shiftedBB_X( shifted_X, nbP, dim); + INTERP_KERNEL::DirectedBoundingBox shiftedBB_y( shifted_y, nbP, dim); + INTERP_KERNEL::DirectedBoundingBox shiftedBB_Y( shifted_Y, nbP, dim); + + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_x )); + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_X )); + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_y )); + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_Y )); + + // intersecting box is IN + double inters_coords[nbP*dim] = + { + 0,0,0, 2,0,0, 2,1,0, 0,1,0, + 0,0,2, 2,0,2, 2,1,2, 0,1,2 + }; + INTERP_KERNEL::DirectedBoundingBox ibb( inters_coords, nbP, dim); + CPPUNIT_ASSERT( !bb.isDisjointWith( ibb )); + + // overlapping non-directed BB + double overlappingBB[2*dim] = + { + 5,6, 0, 1, -5,4 + }; + CPPUNIT_ASSERT( !bb.isDisjointWith( overlappingBB )); + + // non-overlapping non-directed BB + double nonoverlappingBB_1[2*dim] = + { + 5,6, 0,1, -5,2 + }; + CPPUNIT_ASSERT( bb.isDisjointWith( nonoverlappingBB_1 )); + double nonoverlappingBB_2[2*dim] = + { + 5,6, 0,1, 7,20 + }; + CPPUNIT_ASSERT( bb.isDisjointWith( nonoverlappingBB_2 )); + } + + void BBTreeTest::test_DirectedBB_2D() + { + // a segment of length 2 extruded along vector (10,10) + const int nbP = 4, dim = 2; + double coords[nbP*dim] = + { + 0,0, 2,0, + 10,10, 12,10, + }; + INTERP_KERNEL::DirectedBoundingBox bb( coords, nbP, dim); + bb.enlarge( 1e-12 ); + + // corners of extrusion are IN + for ( int i = 0; i < nbP*dim; i+=dim ) + CPPUNIT_ASSERT( !bb.isOut( coords + i )); + + // points near corners of extrusion are OUT + double p[nbP*dim] = + { + 1,2, 4,1, + 11,8, 8,9, + }; + for ( int i = 0; i < nbP*dim; i+=dim ) + CPPUNIT_ASSERT( bb.isOut( p + i )); + + // the extrusions shifted by 3 along OX are OUT + double shifted_X[nbP*dim] = + { + 3,0, 5,0, + 13,10, 15,10, + }; + double shifted_x[nbP*dim] = + { + -3,0, -1,0, + 7,10, 9,10, + }; + INTERP_KERNEL::DirectedBoundingBox shiftedBB_x( shifted_x, nbP, dim); + INTERP_KERNEL::DirectedBoundingBox shiftedBB_X( shifted_X, nbP, dim); + + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_x )); + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_X )); + + // intersecting box is IN + double inters_coords[nbP*dim] = + { + 0,0, 2,0, + 0,2, 2,2 + }; + INTERP_KERNEL::DirectedBoundingBox ibb( inters_coords, nbP, dim); + CPPUNIT_ASSERT( !bb.isDisjointWith( ibb )); + + // overlapping non-directed BB + double overlappingBB[2*dim] = + { + 5,6, -5,4 + }; + CPPUNIT_ASSERT( !bb.isDisjointWith( overlappingBB )); + + // non-overlapping non-directed BB + double nonoverlappingBB_1[2*dim] = + { + 5,6, -5,2 + }; + CPPUNIT_ASSERT( bb.isDisjointWith( nonoverlappingBB_1 )); + double nonoverlappingBB_2[2*dim] = + { + 5,6, 7,20 + }; + CPPUNIT_ASSERT( bb.isDisjointWith( nonoverlappingBB_2 )); + } + + void BBTreeTest::test_DirectedBB_1D() + { + const int nbP = 2, dim = 1; + double coords[nbP*dim] = + { + 0, 10 + }; + INTERP_KERNEL::DirectedBoundingBox bb( coords, nbP, dim); + bb.enlarge( 1e-12 ); + + // coords are IN + for ( int i = 0; i < nbP*dim; i+=dim ) + CPPUNIT_ASSERT( !bb.isOut( coords + i )); + + // points near ends are OUT + double p[nbP*dim] = + { + -0.0001, 10.1 + }; + for ( int i = 0; i < nbP*dim; i+=dim ) + CPPUNIT_ASSERT( bb.isOut( p + i )); + + // shifted boxes are OUT + double shifted_X[nbP*dim] = + { + 10.1, 11 + }; + double shifted_x[nbP*dim] = + { + -3.0, -0.001 + }; + INTERP_KERNEL::DirectedBoundingBox shiftedBB_x( shifted_x, nbP, dim); + INTERP_KERNEL::DirectedBoundingBox shiftedBB_X( shifted_X, nbP, dim); + + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_x )); + CPPUNIT_ASSERT( bb.isDisjointWith( shiftedBB_X )); + + // intersecting box is IN + double inters_coords[nbP*dim] = + { + -2,2 + }; + INTERP_KERNEL::DirectedBoundingBox ibb( inters_coords, nbP, dim); + CPPUNIT_ASSERT( !bb.isDisjointWith( ibb )); + + // overlapping non-directed BB + double overlappingBB[2*dim] = + { + -5,4 + }; + CPPUNIT_ASSERT( !bb.isDisjointWith( overlappingBB )); + + // non-overlapping non-directed BB + double nonoverlappingBB_1[2*dim] = + { + -5,-2 + }; + CPPUNIT_ASSERT( bb.isDisjointWith( nonoverlappingBB_1 )); + double nonoverlappingBB_2[2*dim] = + { + 11,16 + }; + CPPUNIT_ASSERT( bb.isDisjointWith( nonoverlappingBB_2 )); + } } diff --git a/src/INTERP_KERNELTest/BBTreeTest.hxx b/src/INTERP_KERNELTest/BBTreeTest.hxx index 7b7393114..7db6208d9 100644 --- a/src/INTERP_KERNELTest/BBTreeTest.hxx +++ b/src/INTERP_KERNELTest/BBTreeTest.hxx @@ -37,6 +37,9 @@ namespace INTERP_TEST CPPUNIT_TEST_SUITE( BBTreeTest ); CPPUNIT_TEST( test_BBTree ); + CPPUNIT_TEST( test_DirectedBB_1D ); + CPPUNIT_TEST( test_DirectedBB_2D ); + CPPUNIT_TEST( test_DirectedBB_3D ); CPPUNIT_TEST_SUITE_END(); @@ -47,6 +50,9 @@ namespace INTERP_TEST // tests void test_BBTree(); + void test_DirectedBB_1D(); + void test_DirectedBB_2D(); + void test_DirectedBB_3D(); }; diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index 607fbb817..0c0ce1887 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -25,6 +25,7 @@ #include "PlanarIntersector.txx" #include "InterpKernelGeo2DQuadraticPolygon.hxx" #include "InterpKernelGeo2DNode.hxx" +#include "DirectedBoundingBox.hxx" #include #include @@ -550,6 +551,33 @@ bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* return true; } +/*! + * Intersect 2 given Bounding Boxes. + */ +bool MEDCouplingPointSet::intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps) +{ + double* bbtemp = new double[2*dim]; + double deltamax=0.0; + + for (int i=0; i< dim; i++) + { + double delta = bb2[2*i+1]-bb2[2*i]; + if ( delta > deltamax ) + { + deltamax = delta ; + } + } + for (int i=0; i +namespace INTERP_KERNEL +{ + class DirectedBoundingBox; +} + namespace ParaMEDMEM { class DataArrayInt; @@ -76,10 +81,12 @@ namespace ParaMEDMEM virtual void unserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); virtual void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems) = 0; + virtual void giveElemsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector& elems) = 0; virtual DataArrayInt *zipCoordsTraducer() = 0; protected: virtual void checkFullyDefined() const throw(INTERP_KERNEL::Exception) = 0; static bool intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps); + static bool intersectsBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bb1, const double* bb2, int dim, double eps); void rotate2D(const double *center, double angle); void rotate3D(const double *center, const double *vect, double angle); void project2DCellOnXY(const int *startConn, const int *endConn, std::vector& res) const; diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 2261f3486..62e64ce7e 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -24,6 +24,7 @@ #include "InterpolationUtils.hxx" #include "PointLocatorAlgos.txx" #include "BBTree.txx" +#include "DirectedBoundingBox.hxx" #include #include @@ -709,6 +710,58 @@ void MEDCouplingUMesh::giveElemsInBoundingBox(const double *bbox, double eps, st delete [] elem_bb; } +/*! + * Given a boundary box 'bbox' returns elements 'elems' contained in this 'bbox'. + * Warning 'elems' is incremented during the call so if elems is not empty before call returned elements will be + * added in 'elems' parameter. + */ +void MEDCouplingUMesh::giveElemsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector& elems) +{ + if(getMeshDimension()==-1) + { + elems.push_back(0); + return; + } + int dim=getSpaceDimension(); + double* elem_bb=new double[2*dim]; + const int* conn = getNodalConnectivity()->getConstPointer(); + const int* conn_index= getNodalConnectivityIndex()->getConstPointer(); + const double* coords = getCoords()->getConstPointer(); + int nbOfCells=getNumberOfCells(); + for ( int ielem=0; ielem::max(); + elem_bb[i*2+1]=-std::numeric_limits::max(); + } + + for (int inode=conn_index[ielem]+1; inode=0)//avoid polyhedron separator + { + for (int idim=0; idim elem_bb[idim*2+1] ) + { + elem_bb[idim*2+1] = coords[node*dim+idim] ; + } + } + } + } + if (intersectsBoundingBox(bbox, elem_bb, dim, eps)) + { + elems.push_back(ielem); + } + } + delete [] elem_bb; +} + /*! * Returns the cell type of cell with id 'cellId'. */ diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 664fbd28c..074de0b7b 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -73,6 +73,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void renumberNodes(const int *newNodeNumbers, int newNbOfNodes); MEDCOUPLING_EXPORT void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check); MEDCOUPLING_EXPORT void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); + MEDCOUPLING_EXPORT void giveElemsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector& elems); MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildOrthogonalField() const; diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index c077447a0..29f1e58e3 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -238,6 +238,51 @@ void MEDCouplingUMeshDesc::giveElemsInBoundingBox(const double *bbox, double eps delete [] elem_bb; } +void MEDCouplingUMeshDesc::giveElemsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox &bbox, double eps, std::vector& elems) +{ + int dim=getSpaceDimension(); + double* elem_bb=new double[2*dim]; + const int* conn = _desc_connec->getConstPointer(); + const int* conn_index= _desc_connec_index->getConstPointer(); + const int* face = _nodal_connec_face->getConstPointer(); + const int* face_index= _nodal_connec_face_index->getConstPointer(); + const double* coords = getCoords()->getConstPointer(); + int nbOfCells=getNumberOfCells(); + for ( int ielem=0; ielem::max(); + elem_bb[i*2+1]=-std::numeric_limits::max(); + } + + for (int jface=conn_index[ielem]+1; jface elem_bb[idim*2+1] ) + { + elem_bb[idim*2+1] = coords[node*dim+idim] ; + } + } + } + } + if (intersectsBoundingBox(bbox, elem_bb, dim, eps)) + { + elems.push_back(ielem); + } + } + delete [] elem_bb; +} + DataArrayInt *MEDCouplingUMeshDesc::mergeNodes(double precision, bool& areNodesMerged) { //not implemented yet. diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index b4c2d3fb1..91ce649ed 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -52,6 +52,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const; MEDCOUPLING_EXPORT void unserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); MEDCOUPLING_EXPORT void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); + MEDCOUPLING_EXPORT void giveElemsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox &bbox, double eps, std::vector& elems); MEDCOUPLING_EXPORT DataArrayInt *mergeNodes(double precision, bool& areNodesMerged); MEDCOUPLING_EXPORT MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const; MEDCOUPLING_EXPORT MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const; diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index ed638da80..1ae99850c 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -27,6 +27,7 @@ #include "ProcessorGroup.hxx" #include "MPIProcessorGroup.hxx" #include "MEDCouplingFieldDouble.hxx" +#include "DirectedBoundingBox.hxx" #include #include @@ -34,6 +35,8 @@ using namespace std; +#define USE_DIRECTED_BB + namespace ParaMEDMEM { ElementLocator::ElementLocator(const ParaFIELD& sourceField, @@ -85,8 +88,15 @@ namespace ParaMEDMEM return; vector elems; +#ifdef USE_DIRECTED_BB + INTERP_KERNEL::DirectedBoundingBox dbb; + double* distant_bb = _domain_bounding_boxes+rank*dbb.dataSize(_local_cell_mesh_space_dim); + dbb.setData(distant_bb); + _local_cell_mesh->giveElemsInBoundingBox(dbb,getBoundingBoxAdjustment(),elems); +#else double* distant_bb = _domain_bounding_boxes+rank*2*_local_cell_mesh_space_dim; _local_cell_mesh->giveElemsInBoundingBox(distant_bb,getBoundingBoxAdjustment(),elems); +#endif DataArrayInt *distant_ids_send; MEDCouplingPointSet *send_mesh = (MEDCouplingPointSet *)_local_para_field.getField()->buildSubMeshData(&elems[0],&elems[elems.size()],distant_ids_send); @@ -138,8 +148,20 @@ namespace ParaMEDMEM if(spaceDimForAll[i]!=_local_cell_mesh_space_dim && spaceDimForAll[i]!=-1) throw INTERP_KERNEL::Exception("Spacedim not matches !"); delete [] spaceDimForAll; - _domain_bounding_boxes = new double[2*_local_cell_mesh_space_dim*_union_group->size()]; - double * minmax=new double [2*_local_cell_mesh_space_dim]; +#ifdef USE_DIRECTED_BB + INTERP_KERNEL::DirectedBoundingBox dbb; + int bbSize = dbb.dataSize(_local_cell_mesh_space_dim); + _domain_bounding_boxes = new double[bbSize*_union_group->size()]; + if(_local_cell_mesh->getMeshDimension() != -1) + dbb = INTERP_KERNEL::DirectedBoundingBox(_local_cell_mesh->getCoords()->getPointer(), + _local_cell_mesh->getNumberOfNodes(), + _local_cell_mesh_space_dim); + std::vector dbbData = dbb.getData(); + double * minmax= &dbbData[0]; +#else + int bbSize = 2*_local_cell_mesh_space_dim; + _domain_bounding_boxes = new double[bbSize*_union_group->size()]; + double * minmax=new double [bbSize]; if(_local_cell_mesh->getMeshDimension() != -1) _local_cell_mesh->getBoundingBox(minmax); else @@ -148,9 +170,10 @@ namespace ParaMEDMEM minmax[i*2]=-std::numeric_limits::max(); minmax[i*2+1]=std::numeric_limits::max(); } +#endif - comm_interface.allGather(minmax, 2*_local_cell_mesh_space_dim, MPI_DOUBLE, - _domain_bounding_boxes,2*_local_cell_mesh_space_dim, MPI_DOUBLE, + comm_interface.allGather(minmax, bbSize, MPI_DOUBLE, + _domain_bounding_boxes,bbSize, MPI_DOUBLE, *comm); for (int i=0; i< _distant_group.size(); i++) @@ -171,6 +194,12 @@ namespace ParaMEDMEM // ============================================= bool ElementLocator::_intersectsBoundingBox(int irank) { +#ifdef USE_DIRECTED_BB + INTERP_KERNEL::DirectedBoundingBox local_dbb, distant_dbb; + local_dbb.setData( _domain_bounding_boxes+_union_group->myRank()*local_dbb.dataSize( _local_cell_mesh_space_dim )); + distant_dbb.setData( _domain_bounding_boxes+irank*distant_dbb.dataSize( _local_cell_mesh_space_dim )); + return !local_dbb.isDisjointWith( distant_dbb ); +#else double* local_bb = _domain_bounding_boxes+_union_group->myRank()*2*_local_cell_mesh_space_dim; double* distant_bb = _domain_bounding_boxes+irank*2*_local_cell_mesh_space_dim; @@ -182,6 +211,7 @@ namespace ParaMEDMEM if (!intersects) return false; } return true; +#endif } // ====================== -- 2.39.2