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;
}
}
+/*!
+ * 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<DataArrayInt>& elts, MCAuto<DataArrayInt>& 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.
MEDCOUPLING_EXPORT virtual int getCellContainingPoint(const double *pos, double eps) const = 0;
MEDCOUPLING_EXPORT virtual void getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const = 0;
MEDCOUPLING_EXPORT virtual void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const;
+ MEDCOUPLING_EXPORT virtual void getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, int nbOfPoints, double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& 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;
MCAuto<DataArrayInt> 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<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
elts.clear(); elts.insert(elts.end(),eltsUg->begin(),eltsUg->end());
}
+void MEDCouplingUMesh::getCellsContainingPointsZeAlg(const double *pos, int nbOfPoints, double eps,
+ MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex,
+ std::function<bool(INTERP_KERNEL::NormalizedCellType,int)> 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
void MEDCouplingUMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps,
MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& 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<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const
+{
+ auto noImNotSensibleTo2DQuadraticLinearCellsFunc([](INTERP_KERNEL::NormalizedCellType,int) { return false; } );
+ this->getCellsContainingPointsZeAlg(pos,nbOfPoints,eps,elts,eltsIndex,noImNotSensibleTo2DQuadraticLinearCellsFunc);
}
/*!
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<int>& elts) const;
- MEDCOUPLING_EXPORT void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const;
+ MEDCOUPLING_EXPORT void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const override;
+ MEDCOUPLING_EXPORT void getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, int nbOfPoints, double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const override;
MEDCOUPLING_EXPORT void checkButterflyCells(std::vector<int>& cells, double eps=1e-12) const;
MEDCOUPLING_EXPORT DataArrayInt *convexEnvelop2D();
MEDCOUPLING_EXPORT DataArrayInt *findAndCorrectBadOriented3DExtrudedCells();
DataArrayInt *buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayInt *n2o) const;
template<int SPACEDIM>
void getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
- double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const;
+ double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex, std::function<bool(INTERP_KERNEL::NormalizedCellType,int)> sensibilityTo2DQuadraticLinearCellsFunc) const;
+ void getCellsContainingPointsZeAlg(const double *pos, int nbOfPoints, double eps,
+ MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex,
+ std::function<bool(INTERP_KERNEL::NormalizedCellType,int)> sensibilityTo2DQuadraticLinearCellsFunc) const;
/// @cond INTERNAL
static MEDCouplingUMesh *MergeUMeshesLL(const std::vector<const MEDCouplingUMesh *>& a);
typedef int (*DimM1DescNbrer)(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2);
template<int SPACEDIM>
void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
- double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const
+ double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex, std::function<bool(INTERP_KERNEL::NormalizedCellType,int)> sensibilityTo2DQuadraticLinearCellsFunc) const
{
// Override precision for this method only:
INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
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<DummyClsMCUG<SPACEDIM> >::isElementContainsPoint(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps);
else
{
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);
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
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<DataArrayInt> elts,eltsIndex;
+ double val;
+ DataArrayDouble *a;
+ DataArrayDoubleTuple *aa;
+ std::vector<double> 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<double> 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<DataArrayInt> 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<DataArrayInt>& d, MCAuto<DataArrayInt>& 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<DataArrayInt>& d, MCAuto<DataArrayInt>& 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;
self->finishUnserialization(a1,a0,a2);
}
+PyObject *Mesh_getCellsContainingPointsLike(PyObject *p, double eps, const MEDCoupling::MEDCouplingMesh *self, std::function<void(const double *,int,double,MEDCoupling::MCAuto<MEDCoupling::DataArrayInt>&,MEDCoupling::MCAuto<MEDCoupling::DataArrayInt>&)> func)
+{
+ MEDCoupling::MCAuto<MEDCoupling::DataArrayInt> 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<double> 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
