From 23c72a3fb303a8f6c9277ac1e49be1210bce125d Mon Sep 17 00:00:00 2001 From: ageay Date: Mon, 11 Mar 2013 16:23:32 +0000 Subject: [PATCH] DataArrayDouble::findClosestTupleId (for Gauss<->Gauss interpolation) --- src/INTERP_KERNEL/Interpolation.hxx | 2 +- src/INTERP_KERNEL/Interpolation.txx | 4 +- src/MEDCoupling/MEDCouplingMemArray.cxx | 116 +++++++++++++++++++-- src/MEDCoupling/MEDCouplingMemArray.hxx | 5 + src/MEDCoupling_Swig/MEDCouplingCommon.i | 1 + src/MEDCoupling_Swig/MEDCouplingRemapper.i | 18 +++- 6 files changed, 129 insertions(+), 17 deletions(-) diff --git a/src/INTERP_KERNEL/Interpolation.hxx b/src/INTERP_KERNEL/Interpolation.hxx index 3c5eab131..92c38349d 100644 --- a/src/INTERP_KERNEL/Interpolation.hxx +++ b/src/INTERP_KERNEL/Interpolation.hxx @@ -43,7 +43,7 @@ namespace INTERP_KERNEL int fromIntegralUniform(const MyMeshType& meshT, MatrixType& result, const char *method) { return fromToIntegralUniform(false,meshT,result,method); } template int toIntegralUniform(const MyMeshType& meshS, MatrixType& result, const char *method) { return fromToIntegralUniform(true,meshS,result,method); } - static void checkAndSplitInterpolationMethod(const char *method, std::string& srcMeth, std::string& trgMeth) throw(INTERP_KERNEL::Exception); + static void CheckAndSplitInterpolationMethod(const char *method, std::string& srcMeth, std::string& trgMeth) throw(INTERP_KERNEL::Exception); template static double CalculateCharacteristicSizeOfMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, const int printLevel); protected: diff --git a/src/INTERP_KERNEL/Interpolation.txx b/src/INTERP_KERNEL/Interpolation.txx index 201b44671..befbd4714 100644 --- a/src/INTERP_KERNEL/Interpolation.txx +++ b/src/INTERP_KERNEL/Interpolation.txx @@ -56,7 +56,7 @@ namespace INTERP_KERNEL } template - void Interpolation::checkAndSplitInterpolationMethod(const char *method, std::string& srcMeth, std::string& trgMeth) throw(INTERP_KERNEL::Exception) + void Interpolation::CheckAndSplitInterpolationMethod(const char *method, std::string& srcMeth, std::string& trgMeth) throw(INTERP_KERNEL::Exception) { const int NB_OF_METH_MANAGED=4; const char *METH_MANAGED[NB_OF_METH_MANAGED]={"P0P0","P0P1","P1P0","P1P1"}; @@ -66,7 +66,7 @@ namespace INTERP_KERNEL found=(methodC==METH_MANAGED[i]); if(!found) { - std::string msg("The interpolation method : \'"); msg+=method; msg+="\' not managed !"; + std::string msg("The interpolation method : \'"); msg+=method; msg+="\' not managed by INTERP_KERNEL interpolators ! Supported are \"P0P0\", \"P0P1\", \"P1P0\" and \"P1P1\"."; throw INTERP_KERNEL::Exception(msg.c_str()); } srcMeth=methodC.substr(0,2); diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 0bb80568b..04e3cff2a 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -84,6 +84,47 @@ void DataArrayDouble::FindTupleIdsNearTuplesAlg(const BBTree& myTr } } +template +int DataArrayDouble::ComputeBasicClosestTupleIdAlg(const std::vector& elems, const double *thisPt, const double *zePt) +{ + double val=std::numeric_limits::max(); + int ret=-1; + for(std::vector::const_iterator it=elems.begin();it!=elems.end();it++) + { + double tmp=0.; + for(int j=0;j +void DataArrayDouble::FindClosestTupleIdAlg(const BBTree& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res) +{ + double distOpt(dist); + const double *p(pos); + int *r(res); + double bbox[2*SPACEDIM]; + for(int i=0;i::max()); + int nbOfTurn=15; + while(nbOfTurn>0) + { + for(int j=0;j elems; + myTree.getIntersectingElems(bbox,elems); nbOfTurn--; + if(elems.empty()) + { failVal=distOpt; distOpt=okVal==std::numeric_limits::max()?2*distOpt:(okVal+failVal)/2.; continue; } + if( elems.size()<=15 || nbOfTurn>0) + { *r=ComputeBasicClosestTupleIdAlg(elems,thisPt,p); break; } + else + { okVal=distOpt; distOpt=(distOpt+failVal)/2.; continue; } + } + } +} + std::size_t DataArray::getHeapMemorySize() const { std::size_t sz1=_name.capacity(); @@ -1311,6 +1352,65 @@ DataArrayDouble *DataArrayDouble::duplicateEachTupleNTimes(int nbTimes) const th return ret.retn(); } +/*! + * This methods returns for each tuple in \a other which tuple in \a this is the closest. + * So \a this and \a other have to have the same number of components. + * + * \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components. + */ +DataArrayInt *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const throw(INTERP_KERNEL::Exception) +{ + if(!other) + throw INTERP_KERNEL::Exception("DataArrayDouble::findClosestTupleId : other instance is NULL !"); + checkAllocated(); other->checkAllocated(); + int nbOfCompo=getNumberOfComponents(); + if(nbOfCompo!=other->getNumberOfComponents()) + { + std::ostringstream oss; oss << "DataArrayDouble::findClosestTupleId : number of components in this is " << nbOfCompo; + oss << ", whereas number of components in other is " << other->getNumberOfComponents() << "! Should be equal !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + int nbOfTuples=other->getNumberOfTuples(); + int thisNbOfTuples=getNumberOfTuples(); + MEDCouplingAutoRefCountObjectPtr ret=DataArrayInt::New(); ret->alloc(nbOfTuples,1); + double bounds[6]; + getMinMaxPerComponent(bounds); + switch(nbOfCompo) + { + case 3: + { + double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4])); + double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta); + double characSize=pow((delta*delta*delta)/thisNbOfTuples,1./3.); + MEDCouplingAutoRefCountObjectPtr bbox=computeBBoxPerTuple(characSize*1e-12); + BBTree<3,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12); + FindClosestTupleIdAlg<3>(myTree,2.*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer()); + break; + } + case 2: + { + double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])); + double delta=std::max(xDelta,yDelta); + double characSize=sqrt((delta*delta)/thisNbOfTuples); + MEDCouplingAutoRefCountObjectPtr bbox=computeBBoxPerTuple(characSize*1e-12); + BBTree<2,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12); + FindClosestTupleIdAlg<2>(myTree,2.*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer()); + break; + } + case 1: + { + double characSize=fabs(bounds[1]-bounds[0])/thisNbOfTuples; + MEDCouplingAutoRefCountObjectPtr bbox=computeBBoxPerTuple(characSize*1e-12); + BBTree<1,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12); + FindClosestTupleIdAlg<1>(myTree,2.*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer()); + break; + } + default: + throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for findClosestTupleId. Must be 1, 2 or 3."); + } + return ret.retn(); +} + /*! * This method returns a newly allocated object the user should deal with. * This method works for arrays which have number of components into [1,2,3]. If not an exception will be thrown. @@ -5345,13 +5445,11 @@ DataArrayInt *DataArrayInt::getIdsEqualList(const int *valsBg, const int *valsEn const int *cptr=getConstPointer(); std::vector res; int nbOfTuples=getNumberOfTuples(); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); ret->alloc(0,1); for(int i=0;ialloc((int)res.size(),1); - std::copy(res.begin(),res.end(),ret->getPointer()); - return ret; + ret->pushBackSilent(i); + return ret.retn(); } DataArrayInt *DataArrayInt::getIdsNotEqualList(const int *valsBg, const int *valsEnd) const throw(INTERP_KERNEL::Exception) @@ -5362,13 +5460,11 @@ DataArrayInt *DataArrayInt::getIdsNotEqualList(const int *valsBg, const int *val const int *cptr=getConstPointer(); std::vector res; int nbOfTuples=getNumberOfTuples(); + MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); ret->alloc(0,1); for(int i=0;ialloc((int)res.size(),1); - std::copy(res.begin(),res.end(),ret->getPointer()); - return ret; + ret->pushBackSilent(i); + return ret.retn(); } /*! diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index e09e9e8c6..b7205f15f 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -215,6 +215,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void findCommonTuples(double prec, int limitTupleId, DataArrayInt *&comm, DataArrayInt *&commIndex) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *duplicateEachTupleNTimes(int nbTimes) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *getDifferentValues(double prec, int limitTupleId=-1) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayInt *findClosestTupleId(const DataArrayDouble *other) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void setSelectedComponents(const DataArrayDouble *a, const std::vector& compoIds) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void setPartOfValues1(const DataArrayDouble *a, int bgTuples, int endTuples, int stepTuples, int bgComp, int endComp, int stepComp, bool strictCompoCompare=true) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void setPartOfValuesSimple1(double a, int bgTuples, int endTuples, int stepTuples, int bgComp, int endComp, int stepComp) throw(INTERP_KERNEL::Exception); @@ -316,6 +317,10 @@ namespace ParaMEDMEM template void findCommonTuplesAlg(const double *bbox, int nbNodes, int limitNodeId, double prec, DataArrayInt *c, DataArrayInt *cI) const; template + static void FindClosestTupleIdAlg(const BBTree& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res); + template + static int ComputeBasicClosestTupleIdAlg(const std::vector& elems, const double *thisPt, const double *zePt); + template static void FindTupleIdsNearTuplesAlg(const BBTree& myTree, const double *pos, int nbOfTuples, double eps, DataArrayInt *c, DataArrayInt *cI); private: diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 804e96ab0..72ddf3688 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -245,6 +245,7 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::DataArrayDouble::fromCylToCart; %newobject ParaMEDMEM::DataArrayDouble::fromSpherToCart; %newobject ParaMEDMEM::DataArrayDouble::getDifferentValues; +%newobject ParaMEDMEM::DataArrayDouble::findClosestTupleId; %newobject ParaMEDMEM::DataArrayDouble::duplicateEachTupleNTimes; %newobject ParaMEDMEM::DataArrayDouble::__neg__; %newobject ParaMEDMEM::DataArrayDouble::__add__; diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapper.i b/src/MEDCoupling_Swig/MEDCouplingRemapper.i index 3b2359a61..489c9e565 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapper.i +++ b/src/MEDCoupling_Swig/MEDCouplingRemapper.i @@ -42,6 +42,14 @@ using namespace INTERP_KERNEL; namespace ParaMEDMEM { + typedef enum + { + IK_ONLY_PREFERED = 0, + NOT_IK_ONLY_PREFERED = 1, + IK_ONLY_FORCED = 2, + NOT_IK_ONLY_FORCED =3 + } InterpolationMatrixPolicy; + class MEDCouplingRemapper : public TimeLabel, public INTERP_KERNEL::InterpolationOptions { private: @@ -56,16 +64,18 @@ namespace ParaMEDMEM void reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception); - bool setOptionInt(const std::string& key, int value); - bool setOptionDouble(const std::string& key, double value); - bool setOptionString(const std::string& key, const std::string& value); + bool setOptionInt(const std::string& key, int value) throw(INTERP_KERNEL::Exception); + bool setOptionDouble(const std::string& key, double value) throw(INTERP_KERNEL::Exception); + bool setOptionString(const std::string& key, const std::string& value) throw(INTERP_KERNEL::Exception); + int getInterpolationMatrixPolicy() const throw(INTERP_KERNEL::Exception); + void setInterpolationMatrixPolicy(int newInterpMatPol) throw(INTERP_KERNEL::Exception); // int nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs) throw(INTERP_KERNEL::Exception); int nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception); double getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception); %extend { - PyObject *getCrudeMatrix() const + PyObject *getCrudeMatrix() const throw(INTERP_KERNEL::Exception) { const std::vector >& m=self->getCrudeMatrix(); std::size_t sz=m.size(); -- 2.39.2