From 27137ea210f48ee6e5d53d743843a056a4324886 Mon Sep 17 00:00:00 2001 From: abn Date: Thu, 26 Jul 2018 16:02:24 +0200 Subject: [PATCH] Improved CU and UC remapper by caching tetra splitting. Speed gain factor x6. --- src/INTERP_KERNEL/InterpolationCU.txx | 97 +++++++++++++++++++++++ src/INTERP_KERNEL/IntersectorCU3D.hxx | 5 ++ src/INTERP_KERNEL/IntersectorCU3D.txx | 35 ++++++-- src/INTERP_KERNEL/SplitterTetra.hxx | 17 ++-- src/INTERP_KERNEL/SplitterTetra.txx | 14 +++- src/INTERP_KERNEL/TransformedTriangle.cxx | 2 +- src/INTERP_KERNEL/TransformedTriangle.hxx | 2 +- 7 files changed, 153 insertions(+), 19 deletions(-) diff --git a/src/INTERP_KERNEL/InterpolationCU.txx b/src/INTERP_KERNEL/InterpolationCU.txx index 69d0c6ade..57a9093b7 100644 --- a/src/INTERP_KERNEL/InterpolationCU.txx +++ b/src/INTERP_KERNEL/InterpolationCU.txx @@ -29,6 +29,7 @@ #include "IntersectorCU1D.txx" #include "IntersectorCU2D.txx" #include "IntersectorCU3D.txx" +#include "BBTree.txx" #include @@ -119,6 +120,7 @@ namespace INTERP_KERNEL result.resize( intersector->getNumberOfRowsOfResMatrix() ); const int ret = intersector->getNumberOfColsOfResMatrix(); +#ifndef USE_NEW_BB_FOR_CU const double* src_coords[ dim ]; int src_nb_coords[ dim ]; std::map< double, int> src_coord_to_index[ dim ]; @@ -185,6 +187,101 @@ namespace INTERP_KERNEL for ( unsigned int iInd = 0; iInd < structIndices.size(); ++iInd ) intersector->intersectCells( iT, structIndices[iInd], result ); } +#else + // Create BBTree structure + // - get bounding boxes + int numSrcElems = src_mesh.getNumberOfElements(); + int numTargetElems = tgt_mesh.getNumberOfElements(); + + std::vector*> targetElems(numTargetElems); + + // The bounding boxes of the CMesh are built manually: + const double* src_coords[ dim ]; + int src_nb_cells[ dim ]; + for ( int j = 0; j < dim; ++j ) + { + src_coords [j] = src_mesh.getCoordsAlongAxis(j); + src_nb_cells[j] = src_mesh.nbCellsAlongAxis(j); + } + + for(unsigned long i = 0 ; i < numTargetElems ; ++i) + targetElems[i] = new MeshElement(i, tgt_mesh); + + const int bboxStride = dim*2; + double* bboxes = new double[bboxStride * numSrcElems]; + int* srcElemIdx = new int[numSrcElems]; + const int Nx = src_nb_cells[0]; + const int Ny = dim > 1 ? src_nb_cells[1] : 0; + const int Nz = dim > 2 ? src_nb_cells[2] : 0; + for(unsigned long kk = 0; kk < Nx; kk++) + for(unsigned long jj = 0; jj < Ny; jj++) + for(unsigned long ii = 0; ii < Nz; ii++) + { + // The bounding boxes of the CMesh are built manually: + unsigned long i = ii+jj*Nx+kk*Nx*Ny; + bboxes[bboxStride*i+0] = src_coords[0][ii]; // BoundingBox::XMIN + bboxes[bboxStride*i+1] = src_coords[0][ii+1]; // BoundingBox::XMAX + if (dim > 1) + { + bboxes[bboxStride*i+2] = src_coords[1][jj]; // BoundingBox::YMIN + bboxes[bboxStride*i+3] = src_coords[1][jj+1]; // BoundingBox::YMAX + if (dim > 2) + { + bboxes[bboxStride*i+4] = src_coords[2][kk]; // BoundingBox::ZMIN + bboxes[bboxStride*i+5] = src_coords[2][kk+1]; // BoundingBox::ZMAX + } + } + srcElemIdx[i] = i; + } + + BBTree tree(bboxes, srcElemIdx, 0, numSrcElems); + + // for each target element, get source elements with which to calculate intersection + // - calculate intersection by calling intersectCells + for(unsigned long tgtIdx = 0; tgtIdx < numTargetElems; ++tgtIdx) + { + const BoundingBox* box = targetElems[tgtIdx]->getBoundingBox(); + + // get target bbox in right order + double targetBox[bboxStride]; + targetBox[0] = box->getCoordinate(BoundingBox::XMIN); + targetBox[1] = box->getCoordinate(BoundingBox::XMAX); + if (dim > 1) + { + targetBox[2] = box->getCoordinate(BoundingBox::YMIN); + targetBox[3] = box->getCoordinate(BoundingBox::YMAX); + if (dim > 2) + { + targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); + targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + } + } + + std::vector intersectElems; + + tree.getIntersectingElems(targetBox, intersectElems); + + if ( !intersectElems.empty() ) + for (auto elemIdx : intersectElems) + { + unsigned long ii, jj, kk; + std::vector conn(dim); + conn[0] = elemIdx % src_nb_cells[0]; + if (dim > 1) + conn[1] = (elemIdx / src_nb_cells[0]) % src_nb_cells[1]; + if (dim > 2) + conn[2] = elemIdx / (src_nb_cells[0]*src_nb_cells[1]); + intersector->intersectCells(tgtIdx,conn,result); + } + } + + delete [] bboxes; + delete [] srcElemIdx; + + for(unsigned long i = 0 ; i < numTargetElems ; ++i) + delete targetElems[i]; +#endif + delete intersector; return ret; } diff --git a/src/INTERP_KERNEL/IntersectorCU3D.hxx b/src/INTERP_KERNEL/IntersectorCU3D.hxx index 09d746355..63edf77be 100644 --- a/src/INTERP_KERNEL/IntersectorCU3D.hxx +++ b/src/INTERP_KERNEL/IntersectorCU3D.hxx @@ -46,8 +46,13 @@ namespace INTERP_KERNEL typedef SplitterTetra2 TSplitter; typedef SplitterTetra <_Cartesian3D2UnstructHexMesh > TTetra; + + void prepareTetrasForSourceCell(UConnType icellT); + _Cartesian3D2UnstructHexMesh* _uHexMesh; TSplitter* _split; + + std::map< UConnType, std::vector< TTetra* > > _tetraCache; }; } diff --git a/src/INTERP_KERNEL/IntersectorCU3D.txx b/src/INTERP_KERNEL/IntersectorCU3D.txx index 8fc6749e4..3ac9a6ebe 100644 --- a/src/INTERP_KERNEL/IntersectorCU3D.txx +++ b/src/INTERP_KERNEL/IntersectorCU3D.txx @@ -108,7 +108,22 @@ namespace INTERP_KERNEL throw Exception("IntersectorCU3D(): Invalid mesh dimension, it must be 3"); _uHexMesh = new _Cartesian3D2UnstructHexMesh( _INTERSECTOR_CU::_coordsC ); - _split = new TSplitter( meshT, *_uHexMesh, splitting_policy ); + _split = new TSplitter( meshT, *_uHexMesh, splitting_policy ); + } + + IntersectorCU3D_TEMPLATE + void INTERSECTOR_CU3D::prepareTetrasForSourceCell(UConnType icellT) + { + auto it = _tetraCache.find(icellT); + if (it == _tetraCache.end()) + { + UConnType nb_nodes = + _INTERSECTOR_CU::_connIndexU[icellT+1] - _INTERSECTOR_CU::_connIndexU[icellT]; + std::vector< TTetra* > tetra; + _split->releaseArrays(); + _split->splitTargetCell( icellT, nb_nodes, tetra); + _tetraCache[icellT] = tetra; + } } //================================================================================ @@ -122,6 +137,10 @@ namespace INTERP_KERNEL { delete _uHexMesh; _uHexMesh=0; delete _split; _split=0; + + for ( auto tetraIntVec : _tetraCache) + for (auto t : tetraIntVec.second) + delete t; } //================================================================================ @@ -136,20 +155,22 @@ namespace INTERP_KERNEL const std::vector& icellS) { // split an unstructured cell into tetra - std::vector< TTetra* > tetra; - UConnType nb_nodes = - _INTERSECTOR_CU::_connIndexU[icellT+1] - _INTERSECTOR_CU::_connIndexU[icellT]; - _split->releaseArrays(); - _split->splitTargetCell( icellT, nb_nodes, tetra); + prepareTetrasForSourceCell(icellT); + auto pr = _tetraCache.find(icellT); + if (pr == _tetraCache.end()) + throw Exception("INTERSECTOR_CU3D:intersectGeometry() - internal error"); + const std::vector& tetra = pr->second; // intersect a cartesian 3d cell with tetra _uHexMesh->setHexa( _FMIC(icellS[0]),_FMIC(icellS[1]),_FMIC(icellS[2])); // set cell at i,j,k double res = 0; for ( unsigned int t = 0; t < tetra.size(); ++t ) { + tetra[t]->clearVolumesCache(); + tetra[t]->clearNodesCache(); res += tetra[t]->intersectSourceCell( 0 ); - delete tetra[t]; } + return res; } } diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index 5a011c752..7a1f80aba 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -368,7 +368,7 @@ namespace INTERP_KERNEL ~SplitterTetra(); - double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0); + double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0) const; double intersectSourceFace(const NormalizedCellType polyType, const int polyNodesNbr, const int *const polyNodes, @@ -387,13 +387,14 @@ namespace INTERP_KERNEL void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId); void clearVolumesCache(); + void clearNodesCache(); private: inline static void CheckIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL); inline static void CheckIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL); - inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum); + inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum) const; inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node); - inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); + inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key) const; inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key); static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant, @@ -412,13 +413,13 @@ namespace INTERP_KERNEL // member variables /// affine transform associated with this target element - TetraAffineTransform* _t; + const TetraAffineTransform* _t; /// HashMap relating node numbers to transformed nodes, used for caching - HashMap< int , double* > _nodes; + mutable HashMap< int , double* > _nodes; /// HashMap relating triangular faces to calculated volume contributions, used for caching - HashMap< TriangleFaceKey, double > _volumes; + mutable HashMap< TriangleFaceKey, double > _volumes; /// reference to the source mesh const MyMeshType& _src_mesh; @@ -478,7 +479,7 @@ namespace INTERP_KERNEL * */ template - inline void SplitterTetra::calculateNode(typename MyMeshType::MyConnType globalNodeNum) + inline void SplitterTetra::calculateNode(typename MyMeshType::MyConnType globalNodeNum) const { const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; double* transformedNode = new double[MyMeshType::MY_SPACEDIM]; @@ -515,7 +516,7 @@ namespace INTERP_KERNEL * @param key key associated with the face */ template - inline void SplitterTetra::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key) + inline void SplitterTetra::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key) const { const double vol = tri.calculateIntersectionVolume(); _volumes.insert(std::make_pair(key, vol)); diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index d4d0bae74..b7dd40bc5 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -134,6 +134,14 @@ namespace INTERP_KERNEL _volumes.clear(); } + template + void SplitterTetra::clearNodesCache() + { + for(HashMap< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter) + delete[] iter->second; + _nodes.clear(); + } + /*! * This method destroys the 4 pointers pointed by tetraCorners[0],tetraCorners[1],tetraCorners[2] and tetraCorners[3] * @param i is in 0..23 included. @@ -175,7 +183,7 @@ namespace INTERP_KERNEL */ template double SplitterTetra::intersectSourceCell(typename MyMeshType::MyConnType element, - double* baryCentre) + double* baryCentre) const { typedef typename MyMeshType::MyConnType ConnType; const NumberingPolicy numPol=MyMeshType::My_numPol; @@ -210,7 +218,8 @@ namespace INTERP_KERNEL //std::cout << std::endl << "*** " << globalNodeNum << std::endl; calculateNode(globalNodeNum); } - CheckIsOutside(_nodes[globalNodeNum], isOutside); + const double * nn = _nodes[globalNodeNum]; + CheckIsOutside(nn, isOutside); } // halfspace filtering check @@ -920,6 +929,7 @@ namespace INTERP_KERNEL } } _nodes.clear(); + _node_ids.clear(); } /*! diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 693c2d5d7..0fdb7b504 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -238,7 +238,7 @@ namespace INTERP_KERNEL * as described in Grandy * */ - double TransformedTriangle::calculateIntersectionSurface(TetraAffineTransform* tat) + double TransformedTriangle::calculateIntersectionSurface(const TetraAffineTransform* tat) { // check first that we are not below z - plane if(isTriangleBelowTetraeder()) diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index b0e771f7e..170535b06 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -138,7 +138,7 @@ namespace INTERP_KERNEL ~TransformedTriangle(); double calculateIntersectionVolume(); - double calculateIntersectionSurface(TetraAffineTransform* tat); + double calculateIntersectionSurface(const TetraAffineTransform* tat); void dumpCoords() const; -- 2.39.2