-// Copyright (C) 2007-2015 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2021 CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
#include "CurveIntersectorP0P1.txx"
#include "CurveIntersectorP1P1.txx"
#include "CurveIntersectorP1P1PL.txx"
-#include "BBTree.txx"
#include <time.h>
+#include <memory>
namespace INTERP_KERNEL
{
- /**
- * \defgroup interpolationCurve InterpolationCurve
- *
- * \class InterpolationCurve
- * \brief Class used to compute the coefficients of the interpolation matrix between
- * two local meshes in two dimensions.
- */
+
template<class RealCurve>
InterpolationCurve<RealCurve>::InterpolationCurve()
{
*/
template<class RealCurve>
template<class MyMeshType, class MatrixType>
- int InterpolationCurve<RealCurve>::interpolateMeshes (const MyMeshType& myMeshS,
- const MyMeshType& myMeshT,
- MatrixType& result,
- const std::string& method)
+ typename MyMeshType::MyConnType InterpolationCurve<RealCurve>::interpolateMeshesInternal (const MyMeshType& myMeshS,
+ const MyMeshType& myMeshT,
+ MatrixType& result,
+ const std::string& method,
+ std::function< void(const BBTree< MyMeshType::MY_SPACEDIM ,
+ typename MyMeshType::MyConnType>&, const double*,
+ std::vector<typename MyMeshType::MyConnType>&) > bbtreeMethod
+ )
{
static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
typedef typename MyMeshType::MyConnType ConnType;
static const NumberingPolicy numPol = MyMeshType::My_numPol;
long global_start = clock();
- int counter=0;
+ std::size_t counter=0;
- long nbMailleS = myMeshS.getNumberOfElements();
- long nbMailleT = myMeshT.getNumberOfElements();
+ ConnType nbMailleS = myMeshS.getNumberOfElements();
+ ConnType nbMailleT = myMeshT.getNumberOfElements();
- CurveIntersector<MyMeshType,MatrixType>* intersector=0;
+ std::unique_ptr< CurveIntersector<MyMeshType,MatrixType> > intersector;
if(method=="P0P0")
{
switch (InterpolationOptions::getIntersectionType())
{
case Triangulation:
{
- intersector = new CurveIntersectorP0P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP0P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
}
default:
{
case Triangulation:
{
- intersector = new CurveIntersectorP0P1<MyMeshType,MatrixType>(myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP0P1<MyMeshType,MatrixType>(myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
}
default:
{
case Triangulation:
{
- intersector = new CurveIntersectorP1P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP1P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
}
- default:
- throw INTERP_KERNEL::Exception("For P1P0 in 1D or 2D curve only Triangulation supported for the moment !");
+ default:
+ throw INTERP_KERNEL::Exception("For P1P0 in 1D or 2D curve only Triangulation supported for the moment !");
}
}
else if(method=="P1P1")
switch (InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector = new CurveIntersectorP1P1<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP1P1<MyMeshType,MatrixType>
+ (myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
case PointLocator:
- intersector = new CurveIntersectorP1P1PL<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP1P1PL<MyMeshType,MatrixType>
+ (myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
default:
throw INTERP_KERNEL::Exception("For P1P1 in 1D or 2D curve only Triangulation and PointLocator supported !");
throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
/****************************************************************/
/* Create a search tree based on the bounding boxes */
- /* Instanciate the intersector and initialise the result vector */
+ /* Instantiate the intersector and initialise the result vector */
/****************************************************************/
long start_filtering=clock();
std::vector<double> bbox;
intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
- intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustmentAbs());
+ intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustment(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs());
BBTree<SPACEDIM,ConnType> my_tree(&bbox[0], 0, 0, nbMailleS);//creating the search structure
long end_filtering = clock();
/****************************************************/
long start_intersection = clock();
const ConnType *connIndxT = myMeshT.getConnectivityIndexPtr();
- for(int iT=0; iT<nbMailleT; iT++)
+ for(ConnType iT=0; iT<nbMailleT; iT++)
{
- int nb_nodesT = connIndxT[iT+1] - connIndxT[iT];
- std::vector<int> intersecting_elems;
+ ConnType nb_nodesT = connIndxT[iT+1] - connIndxT[iT];
+ std::vector<ConnType> intersecting_elems;
double bb[2*SPACEDIM];
intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
- my_tree.getIntersectingElems(bb, intersecting_elems);
+ bbtreeMethod(my_tree,bb,intersecting_elems);
intersector->intersectCells(iT,intersecting_elems,result);
counter += intersecting_elems.size();
}
- int ret = intersector->getNumberOfColsOfResMatrix();
- delete intersector;
if (InterpolationOptions::getPrintLevel() >= 1)
{
std::cout << "Number of computed intersections = " << counter << std::endl;
std::cout << "Global time= " << global_end - global_start << std::endl;
}
- return ret;
+ return intersector->getNumberOfColsOfResMatrix();
}
}
