#include "IntersectorCU1D.txx"
#include "IntersectorCU2D.txx"
#include "IntersectorCU3D.txx"
+#include "BBTree.txx"
#include <map>
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 ];
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<MeshElement<int>*> 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<int>(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<dim, int> 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<int> intersectElems;
+
+ tree.getIntersectingElems(targetBox, intersectElems);
+
+ if ( !intersectElems.empty() )
+ for (auto elemIdx : intersectElems)
+ {
+ unsigned long ii, jj, kk;
+ std::vector<int> 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;
}
typedef SplitterTetra2<MyUMeshType, _Cartesian3D2UnstructHexMesh > TSplitter;
typedef SplitterTetra <_Cartesian3D2UnstructHexMesh > TTetra;
+
+ void prepareTetrasForSourceCell(UConnType icellT);
+
_Cartesian3D2UnstructHexMesh* _uHexMesh;
TSplitter* _split;
+
+ std::map< UConnType, std::vector< TTetra* > > _tetraCache;
};
}
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;
+ }
}
//================================================================================
{
delete _uHexMesh; _uHexMesh=0;
delete _split; _split=0;
+
+ for ( auto tetraIntVec : _tetraCache)
+ for (auto t : tetraIntVec.second)
+ delete t;
}
//================================================================================
const std::vector<CConnType>& 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<TTetra *>& 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;
}
}
~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,
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,
// 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;
*
*/
template<class MyMeshType>
- inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
+ inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum) const
{
const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
* @param key key associated with the face
*/
template<class MyMeshType>
- inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
+ inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key) const
{
const double vol = tri.calculateIntersectionVolume();
_volumes.insert(std::make_pair(key, vol));
_volumes.clear();
}
+ template<class MyMeshType>
+ void SplitterTetra<MyMeshType>::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.
*/
template<class MyMeshType>
double SplitterTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element,
- double* baryCentre)
+ double* baryCentre) const
{
typedef typename MyMeshType::MyConnType ConnType;
const NumberingPolicy numPol=MyMeshType::My_numPol;
//std::cout << std::endl << "*** " << globalNodeNum << std::endl;
calculateNode(globalNodeNum);
}
- CheckIsOutside(_nodes[globalNodeNum], isOutside);
+ const double * nn = _nodes[globalNodeNum];
+ CheckIsOutside(nn, isOutside);
}
// halfspace filtering check
}
}
_nodes.clear();
+ _node_ids.clear();
}
/*!
* 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())
~TransformedTriangle();
double calculateIntersectionVolume();
- double calculateIntersectionSurface(TetraAffineTransform* tat);
+ double calculateIntersectionSurface(const TetraAffineTransform* tat);
void dumpCoords() const;
