From 9c65a3b1c6eee41fed156c5b27ff245175c72f6c Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Wed, 29 Aug 2018 10:22:11 +0200 Subject: [PATCH] End --- src/INTERP_KERNEL/PointLocatorAlgos.txx | 5 +- 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 | 4 +- .../MEDCouplingBasicsTest6.py | 4 + src/MEDCoupling_Swig/MEDCouplingCommon.i | 51 +++++------ src/MEDCoupling_Swig/MEDCouplingTypemaps.i | 37 ++++++++ 10 files changed, 141 insertions(+), 73 deletions(-) diff --git a/src/INTERP_KERNEL/PointLocatorAlgos.txx b/src/INTERP_KERNEL/PointLocatorAlgos.txx index 07abf8a6b..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; @@ -173,8 +176,6 @@ namespace INTERP_KERNEL // if (SPACEDIM==2) { - if (cmType.isQuadratic()) - throw INTERP_KERNEL::Exception("PointLocatorAlgos.txx::isElementContainsPoint(): Invalid cell type detected! Quadratic types not supported!"); int nbEdges=cmType.getNumberOfSons(); double *pts = new double[nbEdges*SPACEDIM]; for (int iedge=0; iedge& 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 9b37ee294..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); @@ -130,7 +130,7 @@ void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const d INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]]); bool status(false); // [ABN] : point locator algorithms are only properly working for linear cells. - if(ct!=INTERP_KERNEL::NORM_POLYGON && ! (INTERP_KERNEL::CellModel::GetCellModel(ct).isQuadratic() && _mesh_dim == 2)) + 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/MEDCouplingBasicsTest6.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest6.py index 9475b27e9..e7cc79e23 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest6.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest6.py @@ -226,6 +226,10 @@ If true geometry (with curve as edges) is considered the result of getCellsConta 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