From dde8f52747f20c9fb67a781eb79c31f6c3678e12 Mon Sep 17 00:00:00 2001 From: ageay Date: Fri, 5 Mar 2010 15:40:13 +0000 Subject: [PATCH] getElementContainingPoint-getValueOn P0 --- src/MEDCoupling/MEDCouplingCMesh.cxx | 6 + src/MEDCoupling/MEDCouplingCMesh.hxx | 1 + src/MEDCoupling/MEDCouplingExtrudedMesh.cxx | 6 + src/MEDCoupling/MEDCouplingExtrudedMesh.hxx | 1 + .../MEDCouplingFieldDiscretization.cxx | 16 ++ .../MEDCouplingFieldDiscretization.hxx | 10 ++ src/MEDCoupling/MEDCouplingFieldDouble.cxx | 13 +- src/MEDCoupling/MEDCouplingMesh.hxx | 1 + src/MEDCoupling/MEDCouplingPointSet.cxx | 18 +- src/MEDCoupling/MEDCouplingPointSet.hxx | 2 + .../MEDCouplingTimeDiscretization.cxx | 33 +++- .../MEDCouplingTimeDiscretization.hxx | 12 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 170 ++++++++++++++++++ src/MEDCoupling/MEDCouplingUMesh.hxx | 8 + src/MEDCoupling/MEDCouplingUMeshDesc.cxx | 6 + src/MEDCoupling/MEDCouplingUMeshDesc.hxx | 1 + .../Test/MEDCouplingBasicsTest.cxx | 96 ++++++++++ .../Test/MEDCouplingBasicsTest.hxx | 4 + 18 files changed, 385 insertions(+), 19 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingCMesh.cxx b/src/MEDCoupling/MEDCouplingCMesh.cxx index d1b6aec3c..4f4b5150f 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCMesh.cxx @@ -200,6 +200,12 @@ MEDCouplingFieldDouble *MEDCouplingCMesh::buildOrthogonalField() const return ret; } +int MEDCouplingCMesh::getElementContainingPoint(const double *pos, double eps) const +{ + //not implemented yet ! + return -1; +} + void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle) { throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !"); diff --git a/src/MEDCoupling/MEDCouplingCMesh.hxx b/src/MEDCoupling/MEDCouplingCMesh.hxx index 4f4096a44..201a828d5 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.hxx +++ b/src/MEDCoupling/MEDCouplingCMesh.hxx @@ -48,6 +48,7 @@ namespace ParaMEDMEM MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const; MEDCouplingFieldDouble *buildOrthogonalField() const; + int getElementContainingPoint(const double *pos, double eps) const; void rotate(const double *center, const double *vector, double angle); void translate(const double *vector); MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx index a286dc69b..39e440fc2 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx @@ -132,6 +132,12 @@ MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::buildOrthogonalField() const throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::buildOrthogonalField not implemented yet !"); } +int MEDCouplingExtrudedMesh::getElementContainingPoint(const double *pos, double eps) const +{ + //not implemented yet + return -1; +} + MEDCouplingExtrudedMesh::~MEDCouplingExtrudedMesh() { if(_mesh2D) diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx index 8480bcf2d..ee3abeeef 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx @@ -49,6 +49,7 @@ namespace ParaMEDMEM MEDCouplingFieldDouble *getMeasureField(bool) const; MEDCouplingFieldDouble *getMeasureFieldOnNode(bool) const; MEDCouplingFieldDouble *buildOrthogonalField() const; + int getElementContainingPoint(const double *pos, double eps) const; static int findCorrespCellByNodalConn(const std::vector& nodalConnec, const int *revNodalPtr, const int *revNodalIndxPtr) throw(INTERP_KERNEL::Exception); void rotate(const double *center, const double *vector, double angle); diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index f3756291c..266ee4bdd 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -26,6 +26,8 @@ using namespace ParaMEDMEM; +const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12; + const char MEDCouplingFieldDiscretizationP0::REPR[]="P0"; const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS; @@ -34,6 +36,10 @@ const char MEDCouplingFieldDiscretizationP1::REPR[]="P1"; const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES; +MEDCouplingFieldDiscretization::MEDCouplingFieldDiscretization():_precision(DFLT_PRECISION) +{ +} + MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type) { switch(type) @@ -108,6 +114,12 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(cons return mesh->getMeasureField(isAbs); } +void MEDCouplingFieldDiscretizationP0::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const +{ + int id=mesh->getElementContainingPoint(loc,_precision); + arr->getTuple(id,res); +} + /*! * Nothing to do. It's not a bug. */ @@ -183,6 +195,10 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(cons return mesh->getMeasureFieldOnNode(isAbs); } +void MEDCouplingFieldDiscretizationP1::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const +{ +} + void MEDCouplingFieldDiscretizationP1::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const { int oldNbOfElems=old2New->getNbOfElems(); diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index fdb2fd115..0d9b90c26 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -35,6 +35,8 @@ namespace ParaMEDMEM { public: static MEDCouplingFieldDiscretization *New(TypeOfField type); + double getPrecision() const { return _precision; } + void setPrecision(double val) { _precision=val; } static TypeOfField getTypeOfFieldFromStringRepr(const char *repr) throw(INTERP_KERNEL::Exception); virtual TypeOfField getEnum() const = 0; virtual bool isEqual(const MEDCouplingFieldDiscretization *other) const = 0; @@ -45,8 +47,14 @@ namespace ParaMEDMEM virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0; virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0; virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const = 0; + virtual void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const = 0; virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0; virtual void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const = 0; + protected: + MEDCouplingFieldDiscretization(); + protected: + double _precision; + static const double DFLT_PRECISION; }; class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP0 : public MEDCouplingFieldDiscretization @@ -61,6 +69,7 @@ namespace ParaMEDMEM void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const; + void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; public: @@ -80,6 +89,7 @@ namespace ParaMEDMEM void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const; + void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const; static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di); diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index dc2d27652..45d04a2e0 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -165,12 +165,21 @@ void MEDCouplingFieldDouble::measureAccumulate(bool isWAbs, double *res) const void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception) { - _time_discr->checkNoTimePresence(); + const DataArrayDouble *arr=_time_discr->getArray(); + _type->getValueOn(arr,_mesh,spaceLoc,res); } void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception) { - _time_discr->checkTimePresence(time); + std::vector< const DataArrayDouble *> arrs=_time_discr->getArraysForTime(time); + std::vector res2; + for(std::vector< const DataArrayDouble *>::const_iterator iter=arrs.begin();iter!=arrs.end();iter++) + { + int sz=res2.size(); + res2.resize(sz+(*iter)->getNumberOfComponents()); + _type->getValueOn(*iter,_mesh,spaceLoc,&res2[sz]); + } + _time_discr->getValueForTime(time,res2,res); } void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId) diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx index 3bbdec6e2..bd5491fb0 100644 --- a/src/MEDCoupling/MEDCouplingMesh.hxx +++ b/src/MEDCoupling/MEDCouplingMesh.hxx @@ -56,6 +56,7 @@ namespace ParaMEDMEM virtual void getBoundingBox(double *bbox) const = 0; virtual MEDCouplingFieldDouble *getMeasureField(bool isAbs) const = 0; virtual MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const = 0; + virtual int getElementContainingPoint(const double *pos, double eps) const = 0; virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const; virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const; virtual MEDCouplingFieldDouble *buildOrthogonalField() const = 0; diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index 0cac14c02..3a4a4738b 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -401,6 +401,13 @@ bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect' */ void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle) +{ + double *coords=_coords->getPointer(); + int nbNodes=getNumberOfNodes(); + rotate3DAlg(center,vect,angle,nbNodes,coords); +} + +void MEDCouplingPointSet::rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords) { double sina=sin(angle); double cosa=cos(angle); @@ -422,8 +429,6 @@ void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, dou std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies(),sina)); std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus()); //rotation matrix computed. - double *coords=_coords->getPointer(); - int nbNodes=getNumberOfNodes(); double tmp[3]; for(int i=0; igetPointer(); + int nbNodes=getNumberOfNodes(); + rotate2DAlg(center,angle,nbNodes,coords); +} + +void MEDCouplingPointSet::rotate2DAlg(const double *center, double angle, int nbNodes, double *coords) { double cosa=cos(angle); double sina=sin(angle); double matrix[4]; matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa; - double *coords=_coords->getPointer(); - int nbNodes=getNumberOfNodes(); double tmp[2]; for(int i=0; i& nodes) const = 0; diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx index fa900e7a4..0f3337323 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx @@ -395,7 +395,12 @@ void MEDCouplingNoTimeLabel::checkTimePresence(double time) const throw(INTERP_K throw INTERP_KERNEL::Exception(EXCEPTION_MSG); } -DataArrayDouble *MEDCouplingNoTimeLabel::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) +std::vector< const DataArrayDouble *> MEDCouplingNoTimeLabel::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception(EXCEPTION_MSG); +} + +void MEDCouplingNoTimeLabel::getValueForTime(double time, const std::vector& vals, double *res) const { throw INTERP_KERNEL::Exception(EXCEPTION_MSG); } @@ -590,18 +595,23 @@ void MEDCouplingWithTimeStep::checkTimePresence(double time) const throw(INTERP_ } } -DataArrayDouble *MEDCouplingWithTimeStep::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) +std::vector< const DataArrayDouble *> MEDCouplingWithTimeStep::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception) { if(std::fabs(time-_time)<=_time_tolerance) { - if(_array) - _array->incrRef(); - return _array; + std::vector< const DataArrayDouble *> ret(1); + ret[0]=_array; + return ret; } else throw INTERP_KERNEL::Exception(EXCEPTION_MSG); } +void MEDCouplingWithTimeStep::getValueForTime(double time, const std::vector& vals, double *res) const +{ + std::copy(vals.begin(),vals.end(),res); +} + void MEDCouplingWithTimeStep::getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception) { if(std::fabs(time-_time)<=_time_tolerance) @@ -666,18 +676,23 @@ MEDCouplingTimeDiscretization *MEDCouplingConstOnTimeInterval::performCpy(bool d return new MEDCouplingConstOnTimeInterval(*this,deepCpy); } -DataArrayDouble *MEDCouplingConstOnTimeInterval::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) +std::vector< const DataArrayDouble *> MEDCouplingConstOnTimeInterval::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception) { if(time>_start_time-_time_tolerance && time<_end_time+_time_tolerance) { - if(_array) - _array->incrRef(); - return _array; + std::vector< const DataArrayDouble *> ret(1); + ret[0]=_array; + return ret; } else throw INTERP_KERNEL::Exception(EXCEPTION_MSG); } +void MEDCouplingConstOnTimeInterval::getValueForTime(double time, const std::vector& vals, double *res) const +{ + std::copy(vals.begin(),vals.end(),res); +} + bool MEDCouplingConstOnTimeInterval::areCompatible(const MEDCouplingTimeDiscretization *other) const { if(!MEDCouplingTimeDiscretization::areCompatible(other)) diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx index 158c6bb51..07aea3bfb 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx @@ -61,8 +61,8 @@ namespace ParaMEDMEM virtual void setArrays(const std::vector& arrays, TimeLabel *owner) throw(INTERP_KERNEL::Exception); DataArrayDouble *getArray() const { return _array; } virtual DataArrayDouble *getEndArray() const { return _array; } - //! Warning contrary to getArray method this method returns an object to deal with. - virtual DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) = 0; + virtual std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception) = 0; + virtual void getValueForTime(double time, const std::vector& vals, double *res) const = 0; virtual void getArrays(std::vector& arrays) const; virtual bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); virtual bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); @@ -104,7 +104,8 @@ namespace ParaMEDMEM MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const; void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) { } void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception); - DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception); + std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception); + void getValueForTime(double time, const std::vector& vals, double *res) const; bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception); double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception); @@ -144,6 +145,8 @@ namespace ParaMEDMEM double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; } double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; } DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception); + std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception); + void getValueForTime(double time, const std::vector& vals, double *res) const; void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception); void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception); public: @@ -168,7 +171,8 @@ namespace ParaMEDMEM MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const; bool areCompatible(const MEDCouplingTimeDiscretization *other) const; bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const; - DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception); + std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception); + void getValueForTime(double time, const std::vector& vals, double *res) const; void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception); void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception); TypeOfTimeDiscretization getEnum() const { return DISCRETIZATION; } diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 86ad3e08e..e76f933e5 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -21,6 +21,8 @@ #include "CellModel.hxx" #include "VolSurfUser.txx" #include "InterpolationUtils.hxx" +#include "PointLocatorAlgos.txx" +#include "BBTree.txx" #include #include @@ -907,6 +909,140 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const return ret; } +/*{ + const double *coordsPtr=_coords->getConstPointer(); + BBTree myTree(&bbox[0],0,0,nbNodes,-prec); + double bb[2*SPACEDIM]; + double prec2=prec*prec; + for(int i=0;i intersectingElems; + myTree.getIntersectingElems(bb,intersectingElems); + if(intersectingElems.size()>1) + { + std::vector commonNodes; + for(std::vector::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++) + if(*it!=i) + if(INTERP_KERNEL::distance2(coordsPtr+SPACEDIM*i,coordsPtr+SPACEDIM*(*it)) elts; + getElementsContainingPoint(pos,eps,elts); + if(elts.empty()) + return -1; + return elts.front(); +} + +void MEDCouplingUMesh::getElementsContainingPoint(const double *pos, double eps, std::vector& elts) const +{ + std::vector eltsIndex; + getElementsContainingPoints(pos,1,eps,elts,eltsIndex); +} + +namespace ParaMEDMEM +{ + template + class DummyClsMCUG + { + public: + static const int MY_SPACEDIM=SPACEDIMM; + static const int MY_MESHDIM=8; + typedef int MyConnType; + static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE; + }; +} + +template +void MEDCouplingUMesh::getElementsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, + double eps, std::vector& elts, std::vector& eltsIndex) const +{ + std::vector bbox; + eltsIndex.resize(nbOfPoints+1); + eltsIndex[0]=0; + elts.clear(); + getBoundingBoxForBBTree(bbox); + int nbOfCells=getNumberOfCells(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connI=_nodal_connec_index->getConstPointer(); + double bb[2*SPACEDIM]; + for(int j=0;j::max(); bb[2*j+1]=-std::numeric_limits::max(); } + BBTree myTree(&bbox[0],0,0,nbOfCells,-eps); + for(int i=0;i candidates; + myTree.getIntersectingElems(bb,candidates); + for(std::vector::const_iterator iter=candidates.begin();iter!=candidates.end();iter++) + { + int sz=connI[(*iter)+1]-connI[*iter]-1; + if(INTERP_KERNEL::PointLocatorAlgos >::isElementContainsPoint(pos+i*SPACEDIM, + (INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]], + coords,conn+connI[*iter]+1,sz,eps)) + { + eltsIndex[i+1]++; + elts.push_back(*iter); + } + } + } +} + +void MEDCouplingUMesh::getElementsContainingPoints(const double *pos, int nbOfPoints, double eps, + std::vector& elts, std::vector& eltsIndex) const +{ + int spaceDim=getSpaceDimension(); + int mDim=getMeshDimension(); + if(spaceDim==3) + { + if(mDim==3) + { + const double *coords=_coords->getConstPointer(); + getElementsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex); + } + /*else if(mDim==2) + { + + }*/ + 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(); + getElementsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex); + } + else + throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !"); + } + + + +} + /*! * This method is only available for a mesh with meshDim==2 and spaceDim==2||spaceDim==3. * This method returns a vector 'cells' where all detected butterfly cells have been added to cells. @@ -937,6 +1073,40 @@ void MEDCouplingUMesh::checkButterflyCells(std::vector& cells) const } } +/*! + * This method aggregate the bbox of each cell and put it into bbox parameter. + * @param bbox out parameter of size 2*spacedim*nbOfcells. + */ +void MEDCouplingUMesh::getBoundingBoxForBBTree(std::vector& bbox) const +{ + int spaceDim=getSpaceDimension(); + int nbOfCells=getNumberOfCells(); + bbox.resize(2*nbOfCells*spaceDim); + for(int i=0;i::max(); + bbox[2*i+1]=-std::numeric_limits::max(); + } + const double *coordsPtr=_coords->getConstPointer(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connI=_nodal_connec_index->getConstPointer(); + for(int i=0;i=0) + for(int k=0;k& elts) const; + void getElementsContainingPoints(const double *pos, int nbOfPoints, double eps, + std::vector& elts, std::vector& eltsIndex) const; void checkButterflyCells(std::vector& cells) const; + void getBoundingBoxForBBTree(std::vector& bbox) const; bool checkConsecutiveCellTypes() const; MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; DataArrayDouble *getBarycenterAndOwner() const; @@ -85,6 +90,9 @@ namespace ParaMEDMEM void checkFullyDefined() const throw(INTERP_KERNEL::Exception); //tools MEDCouplingUMesh *buildPartOfMySelfKeepCoords(const int *start, const int *end) const; + template + void getElementsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, + double eps, std::vector& elts, std::vector& eltsIndex) const; private: //! this iterator stores current position in _nodal_connec array. mutable int _iterator; diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index 3a34d1f84..feca3dc3c 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -312,3 +312,9 @@ DataArrayDouble *MEDCouplingUMeshDesc::getBarycenterAndOwner() const //not implemented yet. return 0; } + +int MEDCouplingUMeshDesc::getElementContainingPoint(const double *pos, double eps) const +{ + //not implemented yet. + return -1; +} diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index ebd34115a..94876129b 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -57,6 +57,7 @@ namespace ParaMEDMEM void renumberNodes(const int *newNodeNumbers, int newNbOfNodes); MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const; + int getElementContainingPoint(const double *pos, double eps) const; MEDCouplingFieldDouble *buildOrthogonalField() const; DataArrayInt *zipCoordsTraducer(); MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx index c53720bed..0b8f9b70d 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx @@ -1473,6 +1473,102 @@ void MEDCouplingBasicsTest::testBuildOrthogonalField() targetMesh->decrRef(); } +void MEDCouplingBasicsTest::testGetElementsContainingPoint() +{ + MEDCouplingUMesh *targetMesh=build2DTargetMesh_1(); + double pos[12]={0.,0.,0.4,0.4,0.,0.4,0.1,0.1,0.25,0.,0.65,0.}; + std::vector t1,t2; + //2D basic + targetMesh->getElementsContainingPoints(pos,6,1e-12,t1,t2); + CPPUNIT_ASSERT_EQUAL(6,(int)t1.size()); + CPPUNIT_ASSERT_EQUAL(7,(int)t2.size()); + const int expectedValues1[6]={0,4,3,0,1,2}; + const int expectedValues2[7]={0,1,2,3,4,5,6}; + CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues1)); + CPPUNIT_ASSERT(std::equal(t2.begin(),t2.end(),expectedValues2)); + //2D with no help of bounding box. + double center[2]={0.2,0.2}; + MEDCouplingPointSet::rotate2DAlg(center,0.78539816339744830962,6,pos); + targetMesh->rotate(center,0,0.78539816339744830962); + t1.clear(); t2.clear(); + targetMesh->getElementsContainingPoints(pos,6,1e-12,t1,t2); + CPPUNIT_ASSERT_EQUAL(6,(int)t1.size()); + CPPUNIT_ASSERT_EQUAL(7,(int)t2.size()); + CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues1)); + CPPUNIT_ASSERT(std::equal(t2.begin(),t2.end(),expectedValues2)); + //2D outside + const double pos1bis[2]={-0.3303300858899107,-0.11819805153394641}; + CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos1bis,1e-12)); + targetMesh->decrRef(); + //test limits 2D + targetMesh=build2DTargetMesh_1(); + const double pos2[2]={0.2,-0.05}; + t1.clear(); + targetMesh->getElementsContainingPoint(pos2,1e-12,t1); + CPPUNIT_ASSERT_EQUAL(2,(int)t1.size()); + const int expectedValues3[2]={0,1}; + CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues3)); + const double pos3[2]={0.2,0.2}; + t1.clear(); + targetMesh->getElementsContainingPoint(pos3,1e-12,t1); + CPPUNIT_ASSERT_EQUAL(5,(int)t1.size()); + const int expectedValues4[5]={0,1,2,3,4}; + CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues4)); + CPPUNIT_ASSERT_EQUAL(0,targetMesh->getElementContainingPoint(pos3,1e-12)); + targetMesh->decrRef(); + //3D + targetMesh=build3DTargetMesh_1(); + const double pos4[3]={25.,25.,25.}; + CPPUNIT_ASSERT_EQUAL(0,targetMesh->getElementContainingPoint(pos4,1e-12)); + const double pos5[3]={50.,50.,50.}; + t1.clear(); + targetMesh->getElementsContainingPoint(pos5,1e-12,t1); + CPPUNIT_ASSERT_EQUAL(8,(int)t1.size()); + const int expectedValues5[8]={0,1,2,3,4,5,6,7}; + CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues5)); + const double pos6[3]={0., 50., 0.}; + t1.clear(); + targetMesh->getElementsContainingPoint(pos6,1e-12,t1); + CPPUNIT_ASSERT_EQUAL(2,(int)t1.size()); + const int expectedValues6[2]={0,2}; + CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues6)); + //3D outside + const double pos7[3]={-1.0,-1.0,0.}; + CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos7,1e-12)); + //3D outside 2 + const double center2[3]={0.,0.,0.}; + const double vec2[3]={0.,-1.,0.}; + targetMesh->rotate(center2,vec2,0.78539816339744830962); + const double pos8[3]={-25,25.,12.}; + CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos8,1e-12)); + // + targetMesh->decrRef(); +} + +void MEDCouplingBasicsTest::testGetValueOn1() +{ + MEDCouplingUMesh *targetMesh=build2DTargetMesh_1(); + MEDCouplingFieldDouble *fieldOnCells=MEDCouplingFieldDouble::New(ON_CELLS); + int nbOfCells=targetMesh->getNumberOfCells(); + fieldOnCells->setMesh(targetMesh); + DataArrayDouble *array=DataArrayDouble::New(); + array->alloc(nbOfCells,2); + fieldOnCells->setArray(array); + double *tmp=array->getPointer(); + for(int i=0;idecrRef(); + // + const double pos1[2]={0.25,0.}; + double res[2]; + fieldOnCells->getValueOn(pos1,res); + CPPUNIT_ASSERT_DOUBLES_EQUAL(8.,res[0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(18.,res[1],1e-12); + // + fieldOnCells->decrRef(); + targetMesh->decrRef(); +} + void MEDCouplingBasicsTest::test2DInterpP0P0_1() { MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index f2c054659..d1a50bfea 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -62,6 +62,8 @@ namespace ParaMEDMEM CPPUNIT_TEST( testMergeNodesOnField ); CPPUNIT_TEST( testCheckConsecutiveCellTypes ); CPPUNIT_TEST( testBuildOrthogonalField ); + CPPUNIT_TEST( testGetElementsContainingPoint ); + CPPUNIT_TEST( testGetValueOn1 ); CPPUNIT_TEST( test2DInterpP0P0_1 ); CPPUNIT_TEST( test2DInterpP0P0PL_1 ); CPPUNIT_TEST( test2DInterpP0P0PL_2 ); @@ -143,6 +145,8 @@ namespace ParaMEDMEM void testMergeNodesOnField(); void testCheckConsecutiveCellTypes(); void testBuildOrthogonalField(); + void testGetElementsContainingPoint(); + void testGetValueOn1(); void test2DInterpP0P0_1(); void test2DInterpP0P0PL_1(); void test2DInterpP0P0PL_2(); -- 2.39.2