From a6c4fa941f2f989606ab1ccf0d89c5992a0193a9 Mon Sep 17 00:00:00 2001 From: abn Date: Tue, 7 Aug 2018 16:17:35 +0200 Subject: [PATCH] getCellsContainingPoints is sensitive to 2D quadratic static cells. --- src/INTERP_KERNEL/PointLocatorAlgos.txx | 3 + src/MEDCoupling/MEDCouplingMesh.cxx | 12 +++ src/MEDCoupling/MEDCouplingMesh.hxx | 1 + src/MEDCoupling/MEDCouplingRemapper.cxx | 2 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 90 +++++++++++-------- src/MEDCoupling/MEDCouplingUMesh.hxx | 8 +- src/MEDCoupling/MEDCouplingUMesh_internal.hxx | 5 +- .../MEDCouplingBasicsTest3.py | 1 + .../MEDCouplingBasicsTest6.py | 21 +++++ src/MEDCoupling_Swig/MEDCouplingCommon.i | 51 +++++------ src/MEDCoupling_Swig/MEDCouplingTypemaps.i | 37 ++++++++ 11 files changed, 160 insertions(+), 71 deletions(-) diff --git a/src/INTERP_KERNEL/PointLocatorAlgos.txx b/src/INTERP_KERNEL/PointLocatorAlgos.txx index b02403c3e..b3b4ba7c2 100644 --- a/src/INTERP_KERNEL/PointLocatorAlgos.txx +++ b/src/INTERP_KERNEL/PointLocatorAlgos.txx @@ -163,6 +163,9 @@ namespace INTERP_KERNEL return ret; } + /*! + * Precondition : spacedim==meshdim. To be checked upstream to this call. + */ static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords, const typename MyMeshType::MyConnType *conn_elem, int conn_elem_sz, double eps) { const int SPACEDIM=MyMeshType::MY_SPACEDIM; diff --git a/src/MEDCoupling/MEDCouplingMesh.cxx b/src/MEDCoupling/MEDCouplingMesh.cxx index 8dc00b5ef..f0799497e 100644 --- a/src/MEDCoupling/MEDCouplingMesh.cxx +++ b/src/MEDCoupling/MEDCouplingMesh.cxx @@ -662,6 +662,18 @@ void MEDCouplingMesh::getCellsContainingPoints(const double *pos, int nbOfPoints } } +/*! + * Behaves like MEDCouplingMesh::getCellsContainingPoints for cells in \a this that are linear. + * For quadratic cells in \a this, this method behaves by just considering linear part of cells. + * This method is here only for backward compatibility (interpolation GaussPoints to GaussPoints). + * + * \sa MEDCouplingMesh::getCellsContainingPoints, MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss + */ +void MEDCouplingMesh::getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const +{ + this->getCellsContainingPoints(pos,nbOfPoints,eps,elts,eltsIndex); +} + /*! * Writes \a this mesh into a VTK format file named as specified. * \param [in] fileName - the name of the file to write in. If the extension is OK the fileName will be used directly. diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx index a45a41749..537cbd536 100644 --- a/src/MEDCoupling/MEDCouplingMesh.hxx +++ b/src/MEDCoupling/MEDCouplingMesh.hxx @@ -116,6 +116,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT virtual int getCellContainingPoint(const double *pos, double eps) const = 0; MEDCOUPLING_EXPORT virtual void getCellsContainingPoint(const double *pos, double eps, std::vector& elts) const = 0; MEDCOUPLING_EXPORT virtual void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const; + MEDCOUPLING_EXPORT virtual void getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const; MEDCOUPLING_EXPORT virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const; MEDCOUPLING_EXPORT virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, const std::string& func) const; MEDCOUPLING_EXPORT virtual MEDCouplingFieldDouble *fillFromAnalyticCompo(TypeOfField t, int nbOfComp, const std::string& func) const; diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index e9f3eeb66..2aa3f3bf0 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -878,7 +878,7 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss() MCAuto eltsArr,eltsIndexArr; int trgNbOfGaussPts=trgLoc->getNumberOfTuples(); _matrix.resize(trgNbOfGaussPts); - _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr); + _src_ft->getMesh()->getCellsContainingPointsLinearPartOnlyOnNonDynType(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr); const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin()); MCAuto nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex()); MCAuto ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0); diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index ca95dc944..4a8344d8c 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -4081,6 +4081,45 @@ void MEDCouplingUMesh::getCellsContainingPoint(const double *pos, double eps, st elts.clear(); elts.insert(elts.end(),eltsUg->begin(),eltsUg->end()); } +void MEDCouplingUMesh::getCellsContainingPointsZeAlg(const double *pos, int nbOfPoints, double eps, + MCAuto& elts, MCAuto& eltsIndex, + std::function sensibilityTo2DQuadraticLinearCellsFunc) const +{ + int spaceDim(getSpaceDimension()),mDim(getMeshDimension()); + if(spaceDim==3) + { + if(mDim==3) + { + const double *coords=_coords->getConstPointer(); + getCellsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex,sensibilityTo2DQuadraticLinearCellsFunc); + } + else + throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !"); + } + else if(spaceDim==2) + { + if(mDim==2) + { + const double *coords=_coords->getConstPointer(); + getCellsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex,sensibilityTo2DQuadraticLinearCellsFunc); + } + else + throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !"); + } + else if(spaceDim==1) + { + if(mDim==1) + { + const double *coords=_coords->getConstPointer(); + getCellsContainingPointsAlg<1>(coords,pos,nbOfPoints,eps,elts,eltsIndex,sensibilityTo2DQuadraticLinearCellsFunc); + } + else + throw INTERP_KERNEL::Exception("For spaceDim==1 only meshDim==1 implemented for getelementscontainingpoints !"); + } + else + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPoints : not managed for mdim not in [1,2,3] !"); +} + /*! * Finds cells in contact with several balls (i.e. points with precision). * This method is an extension of getCellContainingPoint() and @@ -4113,44 +4152,21 @@ void MEDCouplingUMesh::getCellsContainingPoint(const double *pos, double eps, st void MEDCouplingUMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const { - int spaceDim=getSpaceDimension(); - int mDim=getMeshDimension(); - if(spaceDim==3) - { - if(mDim==3) - { - const double *coords=_coords->getConstPointer(); - getCellsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex); - } - /*else if(mDim==2) - { + auto yesImSensibleTo2DQuadraticLinearCellsFunc([](INTERP_KERNEL::NormalizedCellType ct, int mdim) { return INTERP_KERNEL::CellModel::GetCellModel(ct).isQuadratic() && mdim == 2; } ); + this->getCellsContainingPointsZeAlg(pos,nbOfPoints,eps,elts,eltsIndex,yesImSensibleTo2DQuadraticLinearCellsFunc); +} - }*/ - else - throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !"); - } - else if(spaceDim==2) - { - if(mDim==2) - { - const double *coords=_coords->getConstPointer(); - getCellsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex); - } - else - throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !"); - } - else if(spaceDim==1) - { - if(mDim==1) - { - const double *coords=_coords->getConstPointer(); - getCellsContainingPointsAlg<1>(coords,pos,nbOfPoints,eps,elts,eltsIndex); - } - else - throw INTERP_KERNEL::Exception("For spaceDim==1 only meshDim==1 implemented for getelementscontainingpoints !"); - } - else - throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPoints : not managed for mdim not in [1,2,3] !"); +/*! + * Behaves like MEDCouplingMesh::getCellsContainingPoints for cells in \a this that are linear. + * For quadratic cells in \a this, this method behaves by just considering linear part of cells. + * This method is here only for backward compatibility (interpolation GaussPoints to GaussPoints). + * + * \sa MEDCouplingUMesh::getCellsContainingPoints, MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss + */ +void MEDCouplingUMesh::getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const +{ + auto noImNotSensibleTo2DQuadraticLinearCellsFunc([](INTERP_KERNEL::NormalizedCellType,int) { return false; } ); + this->getCellsContainingPointsZeAlg(pos,nbOfPoints,eps,elts,eltsIndex,noImNotSensibleTo2DQuadraticLinearCellsFunc); } /*! diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index a7be645d8..f5fae2a31 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -170,7 +170,8 @@ namespace MEDCoupling MEDCOUPLING_EXPORT DataArrayDouble *distanceToPoints(const DataArrayDouble *pts, DataArrayInt *& cellIds) const; MEDCOUPLING_EXPORT int getCellContainingPoint(const double *pos, double eps) const; MEDCOUPLING_EXPORT void getCellsContainingPoint(const double *pos, double eps, std::vector& elts) const; - MEDCOUPLING_EXPORT void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const; + MEDCOUPLING_EXPORT void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const override; + MEDCOUPLING_EXPORT void getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, int nbOfPoints, double eps, MCAuto& elts, MCAuto& eltsIndex) const override; MEDCOUPLING_EXPORT void checkButterflyCells(std::vector& cells, double eps=1e-12) const; MEDCOUPLING_EXPORT DataArrayInt *convexEnvelop2D(); MEDCOUPLING_EXPORT DataArrayInt *findAndCorrectBadOriented3DExtrudedCells(); @@ -322,7 +323,10 @@ namespace MEDCoupling DataArrayInt *buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayInt *n2o) const; template void getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, - double eps, MCAuto& elts, MCAuto& eltsIndex) const; + double eps, MCAuto& elts, MCAuto& eltsIndex, std::function sensibilityTo2DQuadraticLinearCellsFunc) const; + void getCellsContainingPointsZeAlg(const double *pos, int nbOfPoints, double eps, + MCAuto& elts, MCAuto& eltsIndex, + std::function sensibilityTo2DQuadraticLinearCellsFunc) const; /// @cond INTERNAL static MEDCouplingUMesh *MergeUMeshesLL(const std::vector& a); typedef int (*DimM1DescNbrer)(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2); diff --git a/src/MEDCoupling/MEDCouplingUMesh_internal.hxx b/src/MEDCoupling/MEDCouplingUMesh_internal.hxx index 59a306753..47e6ac7ae 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_internal.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh_internal.hxx @@ -100,7 +100,7 @@ namespace MEDCoupling template void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, - double eps, MCAuto& elts, MCAuto& eltsIndex) const + double eps, MCAuto& elts, MCAuto& eltsIndex, std::function sensibilityTo2DQuadraticLinearCellsFunc) const { // Override precision for this method only: INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); @@ -129,7 +129,8 @@ void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const d int sz(connI[(*iter)+1]-connI[*iter]-1); INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]]); bool status(false); - if(ct!=INTERP_KERNEL::NORM_POLYGON && ct!=INTERP_KERNEL::NORM_QPOLYG) + // [ABN] : point locator algorithms are only properly working for linear cells. + if(ct!=INTERP_KERNEL::NORM_POLYGON && !sensibilityTo2DQuadraticLinearCellsFunc(ct,_mesh_dim)) status=INTERP_KERNEL::PointLocatorAlgos >::isElementContainsPoint(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps); else { diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py index 267589138..a5d5e16fc 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest3.py @@ -1259,6 +1259,7 @@ class MEDCouplingBasicsTest3(unittest.TestCase): vals1=[1.1,2.1,3.1,4.1,5.2,6.2,7.2,8.2,9.2,10.2] da.setValues(vals1,10,1); f.setArray(da); + f.checkConsistencyLight() # loc=[0.64637931739890486, -0.16185896817550552, 0.22678966365273748] locs=f.getValueOnMulti(loc); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest6.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest6.py index 30b4af304..e7cc79e23 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest6.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest6.py @@ -210,6 +210,27 @@ class MEDCouplingBasicsTest6(unittest.TestCase): ptsPosExp=DataArrayDouble([6.+a,3.+b,3.+a,6.+a,3.,3.+b,6.+b,3.+b,3.+b,7.,3.+b,3.+b,6.+a,6.+a,3.+a,6.+b,6.+a,3.+b,7.,6.+a,3.+b,6.+a,7.,3.+b,6.+a,3.+b,3.,6.+a,6.+a,3.],10,3) self.assertTrue(m.getCoords()[ptsExpToBeModified].isEqual(ptsPosExp,1e-12)) pass + + def testUMeshGetCellsContainingPtOn2DNonDynQuadraticCells(self): + """getCellsContainingPoint is now dealing curves of quadratic 2D elements. +This test is a mesh containing 2 QUAD8 cells. The input point is located at a special loc. +If true geometry (with curve as edges) is considered the result of getCellsContainingPoint is not the same as if only linear part of cells is considered.""" + coords=DataArrayDouble([-0.9428090415820631,0.9428090415820631,-1.06066017177982,1.06066017177982,-1.1785113019775801,1.1785113019775801,-1.2963624321753402,1.2963624321753402,-1.4142135623731,1.41421356237309,-0.7653668647301801,1.8477590650225701,-0.6378057206084831,1.53979922085214,-0.510244576486786,1.23183937668172,-0.701586292669331,1.6937791429373599,-0.574025148547635,1.38581929876693,-0.9259503883660041,1.38578268717091,-0.740760310692803,1.10862614973673,-1.1111404660392,1.66293922460509],13,2) + m=MEDCouplingUMesh("mesh",2) + m.setCoords(coords) + m.allocateCells() + m.insertNextCell(NORM_QUAD8,[4,2,6,5,3,10,8,12]) + m.insertNextCell(NORM_QUAD8,[2,0,7,6,1,11,9,10]) + # + zePt=DataArrayDouble([-0.85863751450784975,1.4203162316045934],1,2) + a,b=m.getCellsContainingPoints(zePt,1e-12) + self.assertTrue(b.isEqual(DataArrayInt([0,1]))) + self.assertTrue(a.isEqual(DataArrayInt([1]))) # <- test is here. 0 if only linear parts are considered. + # + a,b=m.getCellsContainingPointsLinearPartOnlyOnNonDynType(zePt,1e-12) + self.assertTrue(b.isEqual(DataArrayInt([0,1]))) + self.assertTrue(a.isEqual(DataArrayInt([0]))) # <- like before + pass pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 45d2fc793..eb3d38b30 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -714,43 +714,36 @@ namespace MEDCoupling return ret; } - PyObject *getCellsContainingPoints(PyObject *p, double eps) const throw(INTERP_KERNEL::Exception) + PyObject *getCellsContainingPointsLinearPartOnlyOnNonDynType(PyObject *p, int nbOfPoints, double eps) const throw(INTERP_KERNEL::Exception) { - MCAuto elts,eltsIndex; + double val; + DataArrayDouble *a; + DataArrayDoubleTuple *aa; + std::vector bb; + int sw; int spaceDim=self->getSpaceDimension(); - void *da=0; - int res1=SWIG_ConvertPtr(p,&da,SWIGTYPE_p_MEDCoupling__DataArrayDouble, 0 | 0 ); - if (!SWIG_IsOK(res1)) - { - int size; - INTERP_KERNEL::AutoCPtr tmp=convertPyToNewDblArr2(p,&size); - int nbOfPoints=size/spaceDim; - if(size%spaceDim!=0) - { - throw INTERP_KERNEL::Exception("MEDCouplingMesh::getCellsContainingPoints : Invalid list length ! Must be a multiple of self.getSpaceDimension() !"); - } - self->getCellsContainingPoints(tmp,nbOfPoints,eps,elts,eltsIndex); - } - else - { - DataArrayDouble *da2=reinterpret_cast< DataArrayDouble * >(da); - if(!da2) - throw INTERP_KERNEL::Exception("MEDCouplingMesh::getCellsContainingPoints : Not null DataArrayDouble instance expected !"); - da2->checkAllocated(); - int size=da2->getNumberOfTuples(); - int nbOfCompo=da2->getNumberOfComponents(); - if(nbOfCompo!=spaceDim) - { - throw INTERP_KERNEL::Exception("MEDCouplingMesh::getCellsContainingPoints : Invalid DataArrayDouble nb of components ! Expected same as self.getSpaceDimension() !"); - } - self->getCellsContainingPoints(da2->getConstPointer(),size,eps,elts,eltsIndex); - } + const char msg[]="Python wrap of MEDCouplingMesh::getCellsContainingPointsLinearPartOnlyOnNonDynType : "; + const double *pos=convertObjToPossibleCpp5_Safe(p,sw,val,a,aa,bb,msg,nbOfPoints,spaceDim,true); + MCAuto elts,eltsIndex; + self->getCellsContainingPointsLinearPartOnlyOnNonDynType(pos,nbOfPoints,eps,elts,eltsIndex); PyObject *ret=PyTuple_New(2); PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(elts.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 )); PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(eltsIndex.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 )); return ret; } + PyObject *getCellsContainingPoints(PyObject *p, double eps) const throw(INTERP_KERNEL::Exception) + { + auto getCellsContainingPointsFunc=[self](const double *a, int b,double c, MCAuto& d, MCAuto& e) { self->getCellsContainingPoints(a,b,c,d,e); }; + return Mesh_getCellsContainingPointsLike(p,eps,self,getCellsContainingPointsFunc); + } + + PyObject *getCellsContainingPointsLinearPartOnlyOnNonDynType(PyObject *p, double eps) const throw(INTERP_KERNEL::Exception) + { + auto getCellsContainingPointsFunc=[self](const double *a, int b,double c, MCAuto& d, MCAuto& e) { self->getCellsContainingPointsLinearPartOnlyOnNonDynType(a,b,c,d,e); }; + return Mesh_getCellsContainingPointsLike(p,eps,self,getCellsContainingPointsFunc); + } + PyObject *getCellsContainingPoint(PyObject *p, double eps) const throw(INTERP_KERNEL::Exception) { double val; diff --git a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i index b7a4ade9d..09545d6dd 100644 --- a/src/MEDCoupling_Swig/MEDCouplingTypemaps.i +++ b/src/MEDCoupling_Swig/MEDCouplingTypemaps.i @@ -689,4 +689,41 @@ void field__setstate__(typename MEDCoupling::Traits::FieldType *self, PyObjec self->finishUnserialization(a1,a0,a2); } +PyObject *Mesh_getCellsContainingPointsLike(PyObject *p, double eps, const MEDCoupling::MEDCouplingMesh *self, std::function&,MEDCoupling::MCAuto&)> func) +{ + MEDCoupling::MCAuto elts,eltsIndex; + int spaceDim=self->getSpaceDimension(); + void *da=0; + int res1=SWIG_ConvertPtr(p,&da,SWIGTYPE_p_MEDCoupling__DataArrayDouble, 0 | 0 ); + if (!SWIG_IsOK(res1)) + { + int size; + INTERP_KERNEL::AutoCPtr tmp=convertPyToNewDblArr2(p,&size); + int nbOfPoints=size/spaceDim; + if(size%spaceDim!=0) + { + throw INTERP_KERNEL::Exception("MEDCouplingMesh::getCellsContainingPoints : Invalid list length ! Must be a multiple of self.getSpaceDimension() !"); + } + func(tmp,nbOfPoints,eps,elts,eltsIndex); + } + else + { + MEDCoupling::DataArrayDouble *da2=reinterpret_cast< MEDCoupling::DataArrayDouble * >(da); + if(!da2) + throw INTERP_KERNEL::Exception("MEDCouplingMesh::getCellsContainingPoints : Not null DataArrayDouble instance expected !"); + da2->checkAllocated(); + int size=da2->getNumberOfTuples(); + int nbOfCompo=da2->getNumberOfComponents(); + if(nbOfCompo!=spaceDim) + { + throw INTERP_KERNEL::Exception("MEDCouplingMesh::getCellsContainingPoints : Invalid DataArrayDouble nb of components ! Expected same as self.getSpaceDimension() !"); + } + func(da2->getConstPointer(),size,eps,elts,eltsIndex); + } + PyObject *ret=PyTuple_New(2); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(elts.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(eltsIndex.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 )); + return ret; +} + #endif -- 2.39.2