From 39160b669e265ae326fe6be6a264048405d8c228 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Tue, 31 Dec 2019 10:16:34 +0100 Subject: [PATCH] Deal with 3D/0D point locator intersection in Remapper + fix memory leak on getCrudeMatrix --- src/INTERP_KERNEL/BoundingBox.cxx | 13 ++--- src/INTERP_KERNEL/Interpolation3D1D.cxx | 6 +- src/INTERP_KERNEL/Interpolation3D1D.hxx | 4 +- src/INTERP_KERNEL/Interpolation3D1D.txx | 56 +++++++------------ src/MEDCoupling/MEDCouplingRemapper.cxx | 9 +++ .../MEDCouplingRemapperCommon.i | 6 +- .../MEDCouplingRemapperTest.py | 35 ++++++++++++ 7 files changed, 77 insertions(+), 52 deletions(-) diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index 032369073..3057e6349 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -38,19 +38,18 @@ namespace INTERP_KERNEL 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]); } @@ -84,7 +83,7 @@ namespace INTERP_KERNEL */ BoundingBox::~BoundingBox() { - delete[] _coords; + delete [] _coords; } /** diff --git a/src/INTERP_KERNEL/Interpolation3D1D.cxx b/src/INTERP_KERNEL/Interpolation3D1D.cxx index faa65ec42..a08f6f29c 100644 --- a/src/INTERP_KERNEL/Interpolation3D1D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D1D.cxx @@ -30,11 +30,9 @@ namespace INTERP_KERNEL * */ - Interpolation3D1D::Interpolation3D1D() - {} + Interpolation3D1D::Interpolation3D1D() { } - Interpolation3D1D::Interpolation3D1D(const InterpolationOptions& io):Interpolation(io) - {} + Interpolation3D1D::Interpolation3D1D(const InterpolationOptions& io):Interpolation(io) { } /** * Inspired from PlanarIntersector::adjustBoundingBoxes diff --git a/src/INTERP_KERNEL/Interpolation3D1D.hxx b/src/INTERP_KERNEL/Interpolation3D1D.hxx index 59a5cec7c..f04172173 100755 --- a/src/INTERP_KERNEL/Interpolation3D1D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D1D.hxx @@ -18,8 +18,7 @@ // // Author : A Bruneton (CEA/DEN) -#ifndef __INTERPOLATION3D1D_HXX__ -#define __INTERPOLATION3D1D_HXX__ +#pragma once #include "INTERPKERNELDefines.hxx" #include "Interpolation.hxx" @@ -42,4 +41,3 @@ namespace INTERP_KERNEL }; } -#endif diff --git a/src/INTERP_KERNEL/Interpolation3D1D.txx b/src/INTERP_KERNEL/Interpolation3D1D.txx index 9bc67b338..cddf62c4f 100755 --- a/src/INTERP_KERNEL/Interpolation3D1D.txx +++ b/src/INTERP_KERNEL/Interpolation3D1D.txx @@ -16,10 +16,9 @@ // // 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" @@ -32,6 +31,8 @@ #include "BBTree.txx" +#include + namespace INTERP_KERNEL { /** @@ -42,7 +43,7 @@ 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 @@ -51,30 +52,30 @@ namespace INTERP_KERNEL LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); - std::vector*> srcElems(numSrcElems); - std::vector*> targetElems(numTargetElems); + std::vector< std::unique_ptr< MeshElement > > srcElems(numSrcElems); + std::vector< std::unique_ptr< MeshElement > > targetElems(numTargetElems); std::map*, ConnType> indices; for(ConnType i = 0 ; i < numSrcElems ; ++i) - srcElems[i] = new MeshElement(i, srcMesh); + srcElems[i].reset( new MeshElement(i, srcMesh) ); for(ConnType i = 0 ; i < numTargetElems ; ++i) - targetElems[i] = new MeshElement(i, targetMesh); + targetElems[i].reset( new MeshElement(i, targetMesh) ); - Intersector3D* intersector=0; + std::unique_ptr< Intersector3D > intersector; std::string methC = InterpolationOptions::filterInterpolationMethod(method); if(methC=="P0P0") - { intersector=new PointLocator3DIntersectorP0P0(targetMesh, srcMesh, getPrecision()); + { intersector.reset( new PointLocator3DIntersectorP0P0(targetMesh, srcMesh, getPrecision()) ); } else if(methC=="P0P1") - { intersector=new PointLocator3DIntersectorP0P1(targetMesh, srcMesh, getPrecision()); + { intersector.reset( new PointLocator3DIntersectorP0P1(targetMesh, srcMesh, getPrecision()) ); } else if(methC=="P1P0") - { intersector=new PointLocator3DIntersectorP1P0(targetMesh, srcMesh, getPrecision()); + { intersector.reset( new PointLocator3DIntersectorP1P0(targetMesh, srcMesh, getPrecision()) ); } else if(methC=="P1P1") - { intersector=new PointLocator3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()); + { intersector.reset( new PointLocator3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()) ); } else throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\"."); @@ -84,7 +85,7 @@ namespace INTERP_KERNEL // create BBTree structure // - get bounding boxes std::vector bboxes(6*numSrcElems); - ConnType* srcElemIdx = new ConnType[numSrcElems]; + std::unique_ptr srcElemIdx{ new ConnType[numSrcElems] }; for(ConnType i = 0; i < numSrcElems ; ++i) { // get source bboxes in right order @@ -100,10 +101,10 @@ namespace INTERP_KERNEL } 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 @@ -128,25 +129,6 @@ namespace INTERP_KERNEL 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 diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index bc18d38fe..b52977152 100755 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -464,6 +464,15 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() 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) diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i b/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i index 4c0c0188b..e77d21d84 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i @@ -26,6 +26,7 @@ %{ #include "MEDCouplingRemapper.hxx" +#include %} %include "InterpolationOptions.hxx" @@ -77,7 +78,10 @@ namespace MEDCoupling const std::map& row=m[i]; PyObject *ret0=PyDict_New(); for(std::map::const_iterator it=row.begin();it!=row.end();it++) - PyDict_SetItem(ret0,PyInt_FromLong((*it).first),PyFloat_FromDouble((*it).second)); + { + std::unique_ptr> 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; diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index 58054d519..63bf0ee7f 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -1302,6 +1302,41 @@ class MEDCouplingBasicsTest(unittest.TestCase): 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)): -- 2.39.2