template <int dim, class ConnType = int>
class BBTree
{
-
private:
BBTree* _left;
BBTree* _right;
static const int MIN_NB_ELEMS=15;
static const int MAX_LEVEL=20;
public:
-
+ BBTree() = default;
/*!
Constructor of the bounding box tree
\param bbs pointer to the [xmin1 xmax1 ymin1 ymax1 xmin2 xmax2 ...] array containing the bounding boxes that are to be indexed.
BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12);
\endcode
*/
-
BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON):
_left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
{
_right=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_right.size(),_epsilon);
}
-
-
~BBTree()
{
}
-
/*! returns in \a elems the list of elements potentially intersecting the bounding box pointed to by \a bb
\param bb pointer to query bounding box
\param elems list of elements (given in 0-indexing that is to say in \b C \b mode) intersecting the bounding box
*/
-
void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const
{
// terminal node : return list of elements intersecting bb
\param xx pointer to query point coords
\param elems list of elements (given in 0-indexing) intersecting the bounding box
*/
-
void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
{
// terminal node : return list of elements intersecting bb
_right->getElementsAroundPoint(xx,elems);
}
-
-
ConnType size()
{
if (_terminal) return _nbelems;
--- /dev/null
+// Copyright (C) 2022 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
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// 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
+//
+
+#pragma once
+
+#include "BBTree.txx"
+#include <memory>
+
+/*!
+ * Wrapper over BBTree to deal with ownership of bbox double array.
+ */
+template <int dim, class ConnType>
+class BBTreeStandAlone
+{
+private:
+ std::unique_ptr<double[]> _bbox;
+ BBTree<dim,ConnType> _effective;
+public:
+ BBTreeStandAlone(std::unique_ptr<double[]>&& bbs, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON):_bbox(std::move(bbs)),_effective(_bbox.get(),nullptr,0,nbelems,epsilon) { }
+ void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const { _effective.getIntersectingElems(bb,elems); }
+};
inline void dumpCoords() const;
void toCompactData(double data[6]) const;
+
+ BoundingBox& operator=(const BoundingBox& box) = delete;
private:
/// disallow copying
BoundingBox(const BoundingBox& box);
- /// disallow assignment
- BoundingBox& operator=(const BoundingBox& box);
-
/// Vector containing the coordinates of the box
/// interlaced in the order XMIN, YMIN, ZMIN, XMAX, YMAX, ZMAX
double _coords[6];
protected:
template<class MyMeshType, class MyMatrixType>
- void performAdjustmentOfBB(Intersector3D<MyMeshType,MyMatrixType>* intersector, std::vector<double>& bbox) const
+ void performAdjustmentOfBB(Intersector3D<MyMeshType,MyMatrixType>* intersector, double *bbox, std::size_t sz) const
{
- intersector->adjustBoundingBoxes(bbox,InterpolationOptions::getBoundingBoxAdjustment(),InterpolationOptions::getBoundingBoxAdjustmentAbs());
+ intersector->adjustBoundingBoxes(bbox,sz,InterpolationOptions::getBoundingBoxAdjustment(),InterpolationOptions::getBoundingBoxAdjustmentAbs());
}
private:
#include "PointLocator3DIntersectorP1P0.txx"
#include "PolyhedronIntersectorP1P1.txx"
#include "PointLocator3DIntersectorP1P1.txx"
-#include "Log.hxx"
+#include "InterpolationHelper.txx"
#include "BBTree.txx"
{
typedef typename MyMeshType::MyConnType ConnType;
// create MeshElement objects corresponding to each element of the two meshes
- const ConnType numSrcElems = srcMesh.getNumberOfElements();
const ConnType numTargetElems = targetMesh.getNumberOfElements();
- LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
+ LOG(2, "Target mesh has " << numTargetElems << " elements ");
- std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
- std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
-
- std::map<MeshElement<ConnType>*, ConnType> indices;
- DuplicateFacesType intersectFaces;
-
- for(ConnType i = 0 ; i < numSrcElems ; ++i)
- srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
+ std::vector< MeshElement<ConnType> > targetElems(numTargetElems);
+ DuplicateFacesType intersectFaces;
for(ConnType i = 0 ; i < numTargetElems ; ++i)
- targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
+ targetElems[i].assign(i, targetMesh);
- Intersector3D<MyMeshType,MyMatrixType>* intersector=0;
+ std::unique_ptr< Intersector3D<MyMeshType,MyMatrixType> > intersector;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel());
if(methC=="P0P0")
switch(InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,
+ intersector.reset( new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,
srcMesh,
dimCaracteristic,
getPrecision(),
intersectFaces,
- getSplittingPolicy());
+ getSplittingPolicy()) );
break;
default:
throw INTERP_KERNEL::Exception("Invalid 2D to 3D intersection type for P0P0 interp specified : must be Triangulation.");
matrix.resize(intersector->getNumberOfRowsOfResMatrix());
// create BBTree structure
- // - get bounding boxes
- std::vector<double> bboxes(6 * numSrcElems);
- for(ConnType i = 0; i < numSrcElems ; ++i)
- {
- // get source bboxes in right order
- const BoundingBox* box = srcElems[i]->getBoundingBox();
- box->fillInXMinXmaxYminYmaxZminZmaxFormat(bboxes.data()+6*i);
- }
-
// [ABN] Adjust 2D bounding box (those might be flat in the cases where the 2D surf are perfectly aligned with the axis)
- performAdjustmentOfBB(intersector, bboxes);
-
- BBTree<3,ConnType> tree(bboxes.data(), nullptr, 0, numSrcElems, 0.);
+ BBTreeStandAlone<3,ConnType> tree( BuildBBTreeWithAdjustment(srcMesh,[this,&intersector](double *bbox, typename MyMeshType::MyConnType sz){ this->performAdjustmentOfBB(intersector.get(),bbox,sz); }) );
// for each target element, get source elements with which to calculate intersection
// - calculate intersection by calling intersectCells
for(ConnType i = 0; i < numTargetElems; ++i)
{
- const BoundingBox* box = targetElems[i]->getBoundingBox();
+ const BoundingBox* box = targetElems[i].getBoundingBox();
// get target bbox in right order
double targetBox[6];
// free allocated memory
ConnType ret=intersector->getNumberOfColsOfResMatrix();
- delete intersector;
-
- for(ConnType i = 0 ; i < numSrcElems ; ++i)
- {
- delete srcElems[i];
- }
- for(ConnType i = 0 ; i < numTargetElems ; ++i)
- {
- delete targetElems[i];
- }
return ret;
}
#else // use BBTree class
-#include "BBTree.txx"
+#include "InterpolationHelper.txx"
#endif
template<class MyMeshType, class MatrixType>
typename MyMeshType::MyConnType Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
{
- typedef typename MyMeshType::MyConnType ConnType;
+ using ConnType = typename MyMeshType::MyConnType;
// create MeshElement objects corresponding to each element of the two meshes
- const ConnType numSrcElems = srcMesh.getNumberOfElements();
const ConnType numTargetElems = targetMesh.getNumberOfElements();
- LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
+ LOG(2, "Target mesh has " << numTargetElems << " elements ");
- std::vector<std::unique_ptr<MeshElement<ConnType>>> srcElems(numSrcElems);
- std::vector<std::unique_ptr<MeshElement<ConnType>>> targetElems(numTargetElems);
-
- std::map<MeshElement<ConnType>*, ConnType> indices;
-
- for(ConnType i = 0 ; i < numSrcElems ; ++i)
- srcElems[i].reset( new MeshElement<ConnType>(i, srcMesh) );
+ std::vector< MeshElement<ConnType> > targetElems(numTargetElems);
for(ConnType i = 0 ; i < numTargetElems ; ++i)
- targetElems[i].reset( new MeshElement<ConnType>(i, targetMesh) );
+ {
+ targetElems[i].assign(i, targetMesh);
+ }
std::unique_ptr<Intersector3D<MyMeshType,MatrixType>> intersector;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
}
#else // Use BBTree
-
- // create BBTree structure
- // - get bounding boxes
- std::unique_ptr<double[]> bboxes( new double[6 * numSrcElems] );
- for(ConnType i = 0; i < numSrcElems ; ++i)
- {
- // get source bboxes in right order
- const BoundingBox *box = srcElems[i]->getBoundingBox();
- box->fillInXMinXmaxYminYmaxZminZmaxFormat(bboxes.get()+6*i);
- }
-
- BBTree<3,ConnType> tree(bboxes.get(), nullptr, 0, numSrcElems);
+ // create BBTree structure
+ BBTreeStandAlone<3,ConnType> tree( BuildBBTree(srcMesh) );
// for each target element, get source elements with which to calculate intersection
// - calculate intersection by calling intersectCells
for(ConnType i = 0; i < numTargetElems; ++i)
{
- const BoundingBox* box = targetElems[i]->getBoundingBox();
+ const BoundingBox* box = targetElems[i].getBoundingBox();
// get target bbox in right order
double targetBox[6];
Interpolation3D1D::Interpolation3D1D(const InterpolationOptions& io):Interpolation<Interpolation3D1D>(io) { }
- /**
- * Inspired from PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes
- */
- void Interpolation3D1D::adjustBoundingBoxes(std::vector<double>& bbox)
+ void Interpolation3D1D::adjustBoundingBoxes(double *bbox, std::size_t sz)
{
const int SPACE_DIM = 3;
const double adj = getBoundingBoxAdjustmentAbs();
const double adjRel = getBoundingBoxAdjustment();
- std::size_t size = bbox.size()/(2*SPACE_DIM);
+ std::size_t size = sz/(2*SPACE_DIM);
for (std::size_t i=0; i<size; i++)
{
double max=- std::numeric_limits<double>::max();
}
}
}
+
+ /**
+ * Inspired from PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes
+ */
+ void Interpolation3D1D::adjustBoundingBoxes(std::vector<double>& bbox)
+ {
+ adjustBoundingBoxes(bbox.data(),bbox.size());
+ }
}
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-// Author : A Bruneton (CEA/DEN)\r
-\r
-#pragma once\r
-\r
-#include "INTERPKERNELDefines.hxx"\r
-#include "Interpolation.hxx"\r
-#include "NormalizedUnstructuredMesh.hxx"\r
-#include "InterpolationOptions.hxx"\r
-\r
-#include <vector>\r
-\r
-namespace INTERP_KERNEL\r
+// Author : A Bruneton (CEA/DEN)
+
+#pragma once
+
+#include "INTERPKERNELDefines.hxx"
+#include "Interpolation.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "InterpolationOptions.hxx"
+
+#include <vector>
+
+namespace INTERP_KERNEL
{
/**
* \class Interpolation3D1D
* Can be seen as a specialization of Interpolation3D, and allows notably the adjustment of bounding boxes.
*/
- class INTERPKERNEL_EXPORT Interpolation3D1D : public Interpolation<Interpolation3D1D>\r
- {\r
- public:\r
- Interpolation3D1D();\r
- Interpolation3D1D(const InterpolationOptions& io);\r
- template<class MyMeshType, class MatrixType>\r
- typename MyMeshType::MyConnType interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method);\r
- private:\r
- void adjustBoundingBoxes(std::vector<double>& bbox);\r
- };\r
-}\r
-\r
+ class INTERPKERNEL_EXPORT Interpolation3D1D : public Interpolation<Interpolation3D1D>
+ {
+ public:
+ Interpolation3D1D();
+ Interpolation3D1D(const InterpolationOptions& io);
+ template<class MyMeshType, class MatrixType>
+ typename MyMeshType::MyConnType interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method);
+ private:
+ void adjustBoundingBoxes(double *bbox, std::size_t sz);
+ void adjustBoundingBoxes(std::vector<double>& bbox);
+ };
+}
+
#include "PointLocator3DIntersectorP1P1.txx"
#include "Log.hxx"
-#include "BBTree.txx"
-
-#include <memory>
+#include "InterpolationHelper.txx"
namespace INTERP_KERNEL
{
typedef typename MyMeshType::MyConnType ConnType;
// create MeshElement objects corresponding to each element of the two meshes
- const ConnType numSrcElems = srcMesh.getNumberOfElements();
const ConnType numTargetElems = targetMesh.getNumberOfElements();
- LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
-
- std::vector< std::unique_ptr< MeshElement<ConnType> > > srcElems(numSrcElems);
- std::vector< std::unique_ptr< MeshElement<ConnType> > > targetElems(numTargetElems);
-
- std::map<MeshElement<ConnType>*, ConnType> indices;
+ LOG(2, "Target mesh has " << numTargetElems << " elements ");
- for(ConnType i = 0 ; i < numSrcElems ; ++i)
- srcElems[i].reset( new MeshElement<ConnType>(i, srcMesh) );
+ std::vector< MeshElement<ConnType> > targetElems(numTargetElems);
for(ConnType i = 0 ; i < numTargetElems ; ++i)
- targetElems[i].reset( new MeshElement<ConnType>(i, targetMesh) );
+ targetElems[i].assign(i, targetMesh);
std::unique_ptr< Intersector3D<MyMeshType,MatrixType> > intersector;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
// create empty maps for all source elements
result.resize(intersector->getNumberOfRowsOfResMatrix());
- // create BBTree structure
- // - get bounding boxes
- std::vector<double> bboxes(6*numSrcElems);
- for(ConnType i = 0; i < numSrcElems ; ++i)
- {
- // get source bboxes in right order
- const BoundingBox* box = srcElems[i]->getBoundingBox();
- box->fillInXMinXmaxYminYmaxZminZmaxFormat(bboxes.data()+6*i);
- }
-
- adjustBoundingBoxes(bboxes);
- const double *bboxPtr = nullptr;
- if(numSrcElems>0)
- bboxPtr=bboxes.data();
- BBTree<3,ConnType> tree(bboxPtr, nullptr, 0, numSrcElems);
+ BBTreeStandAlone<3,ConnType> tree( BuildBBTreeWithAdjustment(srcMesh,[this](double *bbox, typename MyMeshType::MyConnType sz){ this->adjustBoundingBoxes(bbox,sz); }) );
// for each target element, get source elements with which to calculate intersection
// - calculate intersection by calling intersectCells
for(ConnType i = 0; i < numTargetElems; ++i)
{
- const BoundingBox* box = targetElems[i]->getBoundingBox();
+ const BoundingBox* box = targetElems[i].getBoundingBox();
// get target bbox in right order
double targetBox[6];
box->fillInXMinXmaxYminYmaxZminZmaxFormat(targetBox);
--- /dev/null
+// Copyright (C) 2022 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
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// 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
+//
+// Author : Anthony Geay (EDF R&D)
+
+#pragma once
+
+#include "BBTreeStandAlone.txx"
+#include "Log.hxx"
+
+#include <memory>
+#include <functional>
+
+namespace INTERP_KERNEL
+{
+ template<class MyMeshType>
+ BBTreeStandAlone<3,typename MyMeshType::MyConnType> BuildBBTree(const MyMeshType& srcMesh)
+ {
+ return BuildBBTreeWithAdjustment(srcMesh,[](double *,typename MyMeshType::MyConnType){});
+ }
+
+ template<class MyMeshType>
+ BBTreeStandAlone<3,typename MyMeshType::MyConnType> BuildBBTreeWithAdjustment(const MyMeshType& srcMesh, std::function<void(double *,typename MyMeshType::MyConnType)> bboxAdjuster)
+ {
+ using ConnType = typename MyMeshType::MyConnType;
+ const ConnType numSrcElems = srcMesh.getNumberOfElements();
+ LOG(2, "Source mesh has " << numSrcElems << " elements");
+ std::vector< MeshElement<ConnType> > srcElems(numSrcElems);
+
+ for(ConnType i = 0 ; i < numSrcElems ; ++i)
+ {
+ srcElems[i].assign(i,srcMesh);
+ }
+ // create BBTree structure
+ // - get bounding boxes
+ const ConnType nbElts = 6 * numSrcElems;
+ std::unique_ptr<double[]> bboxes( new double[nbElts] );
+ for(ConnType i = 0; i < numSrcElems ; ++i)
+ {
+ // get source bboxes in right order
+ const BoundingBox *box = srcElems[i].getBoundingBox();
+ box->fillInXMinXmaxYminYmaxZminZmaxFormat(bboxes.get()+6*i);
+ }
+ bboxAdjuster(bboxes.get(),nbElts);
+ return BBTreeStandAlone<3,ConnType>(std::move(bboxes),numSrcElems);
+ }
+}
template<class MyMeshType>
MeshElement(const ConnType index, const MyMeshType& mesh);
+ MeshElement() = default;
+
+ template<class MyMeshType>
+ void assign(const ConnType index, const MyMeshType& mesh);
+
~MeshElement() { }
nbnodesincelltype getNumberOfNodes() const { return _number; }
const BoundingBox *getBoundingBox() const { return &_box; }
+ MeshElement& operator=(const MeshElement& elem) = delete;
+
private:
/// disallow copying
MeshElement(const MeshElement& elem);
- /// disallow assignment
- MeshElement& operator=(const MeshElement& elem);
-
nbnodesincelltype _number;
/// bounding box of the element - does not change after having been initialised
*/
template<class ConnType>
template<class MyMeshType>
- MeshElement<ConnType>::MeshElement(const ConnType index, const MyMeshType& mesh)
- : _number( 0 )
+ MeshElement<ConnType>::MeshElement(const ConnType index, const MyMeshType& mesh): _number( 0 )
+ {
+ this->assign(index,mesh);
+ }
+
+ template<class ConnType>
+ template<class MyMeshType>
+ void MeshElement<ConnType>::assign(const ConnType index, const MyMeshType& mesh)
{
auto numberCore = mesh.getNumberOfNodesOfElement(OTT<typename MyMeshType::MyConnType,MyMeshType::My_numPol>::indFC(index));
if(numberCore < std::numeric_limits<nbnodesincelltype>::max())
std::unique_ptr<const double*[]> vertices( new const double*[_number] );
for( nbnodesincelltype i = 0 ; i < _number ; ++i)
vertices[i] = getCoordsOfNode(i , OTT<typename MyMeshType::MyConnType,MyMeshType::My_numPol>::indFC(index), mesh);
-
// create bounding box
_box.initializeWith(vertices.get(),_number);
}
virtual ~TargetIntersector() { }
void adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEps, double adjustmentEpsAbs);
+ void adjustBoundingBoxes(double *bbox, std::size_t sz, double adjustmentEps, double adjustmentEpsAbs);
};
}
namespace INTERP_KERNEL
{
+ template<class MyMeshType, class MyMatrix>
+ void TargetIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEps, double adjustmentEpsAbs)
+ {
+ this->adjustBoundingBoxes(bbox.data(),bbox.size(),adjustmentEps,adjustmentEpsAbs);
+ }
+
/*! Readjusts a set of bounding boxes so that they are extended
in all dimensions for avoiding missing interesting intersections
@param adjustmentEpsAbs absolute adjustment value (added on each side of the BBox in each dimension)
*/
template<class MyMeshType, class MyMatrix>
- void TargetIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEps, double adjustmentEpsAbs)
+ void TargetIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(double *bbox, std::size_t sz, double adjustmentEps, double adjustmentEpsAbs)
{
- std::size_t size = bbox.size()/(2*SPACEDIM);
+ std::size_t size = sz/(2*SPACEDIM);
for (std::size_t i=0; i<size; i++)
{
double max=- std::numeric_limits<double>::max();
