BoundingBox::BoundingBox(const double** pts, const unsigned numPts)
:_coords(new double[6])
{
- assert(numPts > 1);
+ assert(numPts > 0);
// initialize with first two points
- const double* pt1 = pts[0];
- const double* pt2 = pts[1];
+ const double *pt0(pts[0]);
for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1))
{
- _coords[c] = std::min(pt1[c], pt2[c]);
- _coords[c + 3] = std::max(pt1[c], pt2[c]);
+ _coords[c] = pt0[c];
+ _coords[c + 3] = pt0[c];
}
- for(unsigned i = 2 ; i < numPts ; ++i)
+ for(unsigned i = 1 ; i < numPts ; ++i)
{
updateWithPoint(pts[i]);
}
*/
BoundingBox::~BoundingBox()
{
- delete[] _coords;
+ delete [] _coords;
}
/**
*
*/
- Interpolation3D1D::Interpolation3D1D()
- {}
+ Interpolation3D1D::Interpolation3D1D() { }
- Interpolation3D1D::Interpolation3D1D(const InterpolationOptions& io):Interpolation<Interpolation3D1D>(io)
- {}
+ Interpolation3D1D::Interpolation3D1D(const InterpolationOptions& io):Interpolation<Interpolation3D1D>(io) { }
/**
* Inspired from PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes
//\r
// Author : A Bruneton (CEA/DEN)\r
\r
-#ifndef __INTERPOLATION3D1D_HXX__\r
-#define __INTERPOLATION3D1D_HXX__\r
+#pragma once\r
\r
#include "INTERPKERNELDefines.hxx"\r
#include "Interpolation.hxx"\r
};\r
}\r
\r
-#endif\r
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-// Author : Anthony Geay (CEA/DEN)
+// Author : Anthony Geay (EDF R&D)
-#ifndef __INTERPOLATION3D1D_TXX__
-#define __INTERPOLATION3D1D_TXX__
+#pragma once
#include "Interpolation3D1D.hxx"
#include "Interpolation.txx"
#include "BBTree.txx"
+#include <memory>
+
namespace INTERP_KERNEL
{
/**
typename MyMeshType::MyConnType Interpolation3D1D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
{
if(InterpolationOptions::getIntersectionType() != PointLocator)
- INTERP_KERNEL::Exception("Invalid 3D/1D intersection type specified : must be PointLocator.");
+ INTERP_KERNEL::Exception("Invalid 3D/1D-0D intersection type specified : must be PointLocator.");
typedef typename MyMeshType::MyConnType ConnType;
// create MeshElement objects corresponding to each element of the two meshes
LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
- std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
- std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
+ 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] = new MeshElement<ConnType>(i, srcMesh);
+ srcElems[i].reset( new MeshElement<ConnType>(i, srcMesh) );
for(ConnType i = 0 ; i < numTargetElems ; ++i)
- targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
+ targetElems[i].reset( new MeshElement<ConnType>(i, targetMesh) );
- Intersector3D<MyMeshType,MatrixType>* intersector=0;
+ std::unique_ptr< Intersector3D<MyMeshType,MatrixType> > intersector;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
if(methC=="P0P0")
- { intersector=new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ { intersector.reset( new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
}
else if(methC=="P0P1")
- { intersector=new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ { intersector.reset( new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
}
else if(methC=="P1P0")
- { intersector=new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ { intersector.reset( new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
}
else if(methC=="P1P1")
- { intersector=new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+ { intersector.reset( new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
}
else
throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\".");
// create BBTree structure
// - get bounding boxes
std::vector<double> bboxes(6*numSrcElems);
- ConnType* srcElemIdx = new ConnType[numSrcElems];
+ std::unique_ptr<ConnType[]> srcElemIdx{ new ConnType[numSrcElems] };
for(ConnType i = 0; i < numSrcElems ; ++i)
{
// get source bboxes in right order
}
adjustBoundingBoxes(bboxes);
- const double *bboxPtr=0;
+ const double *bboxPtr = nullptr;
if(numSrcElems>0)
- bboxPtr=&bboxes[0];
- BBTree<3,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems);
+ bboxPtr=bboxes.data();
+ BBTree<3,ConnType> tree(bboxPtr, srcElemIdx.get(), 0, numSrcElems);
// for each target element, get source elements with which to calculate intersection
// - calculate intersection by calling intersectCells
if ( !intersectElems.empty() )
intersector->intersectCells(targetIdx,intersectElems,result);
}
-
- // free allocated memory
- delete [] srcElemIdx;
-
- 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;
-
+ return intersector->getNumberOfColsOfResMatrix();
}
}
-
-#endif
INTERP_KERNEL::Interpolation3D1D interpolation(*this);
nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
}
+ else if(srcMeshDim==3 && trgMeshDim==0 && srcSpaceDim==3)
+ {
+ if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
+ throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 0D ! Select PointLocator as intersection type !");
+ MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
+ INTERP_KERNEL::Interpolation3D1D interpolation(*this);//Not a bug : 3D1D deal with 3D0D
+ nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
+ }
else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
{
if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
%{
#include "MEDCouplingRemapper.hxx"
+#include <memory>
%}
%include "InterpolationOptions.hxx"
const std::map<mcIdType,double>& row=m[i];
PyObject *ret0=PyDict_New();
for(std::map<mcIdType,double>::const_iterator it=row.begin();it!=row.end();it++)
- PyDict_SetItem(ret0,PyInt_FromLong((*it).first),PyFloat_FromDouble((*it).second));
+ {
+ std::unique_ptr<PyObject,std::function<void(PyObject*)>> k(PyInt_FromLong((*it).first),[](PyObject *obj) { Py_XDECREF(obj); } ),v(PyFloat_FromDouble((*it).second),[](PyObject *obj) { Py_XDECREF(obj); } );
+ PyDict_SetItem(ret0,k.get(),v.get());
+ }
PyList_SetItem(ret,i,ret0);
}
return ret;
self.assertTrue(abs(res-ref)/ref<1e-12)
pass
+ def test3D0DP1P1(self):
+ """
+ For pointlocator fans, Remapper support following intersection
+ IntersectionType == PointLocator
+ - source == 3D
+ - target == 0D
+ """
+ src = MEDCouplingUMesh("src",3)
+ src.allocateCells()
+ src.setCoords( DataArrayDouble([(0,0,0),(1,0,0),(0,1,0),(0,0,1)]) )
+ src.insertNextCell(NORM_TETRA4,[0,1,2,3])
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([(0.4,0.3,0.07)]) )
+ # P1P1
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.23,1:0.4,2:0.3,3:0.07}],src.getNumberOfNodes(),1e-12)
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.23,1:0.4,2:0.3,3:0.07}],src.getNumberOfNodes(),1e-12)
+ # P1P0
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P0")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.23,1:0.4,2:0.3,3:0.07}],src.getNumberOfNodes(),1e-12)
+ # P0P1
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P0P1")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:1.0}],src.getNumberOfCells(),1e-12)
+ # P0P0
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P0P0")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:1.0}],src.getNumberOfCells(),1e-12)
+ pass
+
def checkMatrix(self,mat1,mat2,nbCols,eps):
self.assertEqual(len(mat1),len(mat2))
for i in range(len(mat1)):
