From 5e9836f2e98cc70c437b2da1acec483fc8af21a7 Mon Sep 17 00:00:00 2001 From: geay Date: Fri, 25 Apr 2014 17:26:38 +0200 Subject: [PATCH] open kriging API (2) --- .../MEDCouplingFieldDiscretization.cxx | 58 ++++++++++++++----- .../MEDCouplingFieldDiscretization.hxx | 10 ++-- .../MEDCouplingFieldDiscretization.i | 13 +++++ 3 files changed, 61 insertions(+), 20 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index 83720b7ab..76344523e 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -2822,24 +2822,36 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeEvaluationMatrixO * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients. If different from 0 there is presence of drift coefficients. * \param [out] matSz the size of returned square matrix * \return the new result matrix to be deallocated by the caller. + * \sa computeMatrix */ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const { - if(!mesh) - throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients : NULL input mesh !"); - MEDCouplingAutoRefCountObjectPtr coords=getLocalizationOfDiscValues(mesh); - int nbOfPts=coords->getNumberOfTuples(); - MEDCouplingAutoRefCountObjectPtr matrix=coords->buildEuclidianDistanceDenseMatrix(); - operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer()); - // Drift - MEDCouplingAutoRefCountObjectPtr matrixWithDrift=performDrift(matrix,coords,isDrift); - MEDCouplingAutoRefCountObjectPtr matrixInv=DataArrayDouble::New(); - matSz=nbOfPts+isDrift; + MEDCouplingAutoRefCountObjectPtr matrixWithDrift(computeMatrix(mesh,isDrift,matSz)); + MEDCouplingAutoRefCountObjectPtr matrixInv(DataArrayDouble::New()); matrixInv->alloc(matSz*matSz,1); INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),matSz,matrixInv->getPointer()); return matrixInv.retn(); } +/*! + * This method computes the kriging matrix. + * \return the new result matrix to be deallocated by the caller. + * \sa computeInverseMatrix + */ +DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const +{ + if(!mesh) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeMatrix : NULL input mesh !"); + MEDCouplingAutoRefCountObjectPtr coords(getLocalizationOfDiscValues(mesh)); + int nbOfPts(coords->getNumberOfTuples()); + MEDCouplingAutoRefCountObjectPtr matrix(coords->buildEuclidianDistanceDenseMatrix()); + operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer()); + // Drift + MEDCouplingAutoRefCountObjectPtr matrixWithDrift(performDrift(matrix,coords,isDrift)); + matSz=nbOfPts+isDrift; + return matrixWithDrift.retn(); +} + /*! * This method computes coefficients to apply to each representing points of \a mesh, that is to say the nodes of \a mesh given a field array \a arr whose * number of tuples should be equal to the number of representing points in \a mesh. @@ -2854,13 +2866,10 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficie { int nbRows(-1); MEDCouplingAutoRefCountObjectPtr matrixInv(computeInverseMatrix(mesh,isDrift,nbRows)); - MEDCouplingAutoRefCountObjectPtr KnewiK=DataArrayDouble::New(); + MEDCouplingAutoRefCountObjectPtr KnewiK(DataArrayDouble::New()); KnewiK->alloc(nbRows*1,1); - MEDCouplingAutoRefCountObjectPtr arr2=DataArrayDouble::New(); - arr2->alloc(nbRows*1,1); - double *work=std::copy(arr->begin(),arr->end(),arr2->getPointer()); - std::fill(work,work+isDrift,0.); - INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbRows,nbRows,arr2->getConstPointer(),nbRows,1,KnewiK->getPointer()); + MEDCouplingAutoRefCountObjectPtr arr2(PerformDriftOfVec(arr,isDrift)); + INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbRows,nbRows,arr2->getConstPointer(),arr2->getNumberOfTuples(),1,KnewiK->getPointer()); return KnewiK.retn(); } @@ -2945,6 +2954,23 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftRect(const D return ret.retn(); } +/*! + * \return a newly allocated array having \a isDrift more tuples than \a arr. + * \sa computeVectorOfCoefficients + */ +DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec(const DataArrayDouble *arr, int isDrift) +{ + if(!arr || !arr->isAllocated() || arr->getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec : input array must be not NULL allocated and with one component !"); + if(isDrift<0) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec : isDrift parameter must be >=0 !"); + MEDCouplingAutoRefCountObjectPtr arr2(DataArrayDouble::New()); + arr2->alloc((arr->getNumberOfTuples()+isDrift)*1,1); + double *work(std::copy(arr->begin(),arr->end(),arr2->getPointer())); + std::fill(work,work+isDrift,0.); + return arr2.retn(); +} + /*! * Starting from a square matrix \a matr, this method returns a newly allocated dense square matrix whose \a matr is included in returned matrix * in the top left corner, and in the remaining returned matrix the parameters to take into account about the kriging drift. diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index 40afef593..5bb752e0a 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -405,15 +405,17 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; MEDCOUPLING_EXPORT DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const; MEDCOUPLING_EXPORT void reprQuickOverview(std::ostream& stream) const; + public://specific part + MEDCOUPLING_EXPORT DataArrayDouble *computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints, int& nbCols) const; + MEDCOUPLING_EXPORT DataArrayDouble *computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const; + MEDCOUPLING_EXPORT DataArrayDouble *computeMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const; + MEDCOUPLING_EXPORT DataArrayDouble *computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const; MEDCOUPLING_EXPORT void operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const; MEDCOUPLING_EXPORT DataArrayDouble *performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const; MEDCOUPLING_EXPORT static void OperateOnDenseMatrixH3(int nbOfElems, double *matrixPtr); MEDCOUPLING_EXPORT static void OperateOnDenseMatrixH2Ln(int nbOfElems, double *matrixPtr); MEDCOUPLING_EXPORT static DataArrayDouble *PerformDriftRect(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta); - public://specific part - MEDCOUPLING_EXPORT DataArrayDouble *computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints, int& nbCols) const; - MEDCOUPLING_EXPORT DataArrayDouble *computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const; - MEDCOUPLING_EXPORT DataArrayDouble *computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const; + MEDCOUPLING_EXPORT static DataArrayDouble *PerformDriftOfVec(const DataArrayDouble *arr, int isDrift); public: static const char REPR[]; static const TypeOfField TYPE; diff --git a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i index dc9da4538..5394a2cc2 100644 --- a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i +++ b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i @@ -28,6 +28,7 @@ %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::clonePart; %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::getValueOnMulti; %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::computeTupleIdsToSelectFromCellIds; +%newobject ParaMEDMEM::MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec; namespace ParaMEDMEM { @@ -374,6 +375,7 @@ namespace ParaMEDMEM class MEDCouplingFieldDiscretizationKriging : public MEDCouplingFieldDiscretizationOnNodes { public: + static DataArrayDouble *PerformDriftOfVec(const DataArrayDouble *arr, int isDrift) throw(INTERP_KERNEL::Exception); %extend { PyObject *computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr) const throw(INTERP_KERNEL::Exception) @@ -396,6 +398,17 @@ namespace ParaMEDMEM PyTuple_SetItem(ret,2,PyInt_FromLong(ret2)); return ret; } + + PyObject *computeMatrix(const MEDCouplingMesh *mesh) const throw(INTERP_KERNEL::Exception) + { + int ret1(-1),ret2(-1); + DataArrayDouble *ret0=self->computeMatrix(mesh,ret1,ret2); + PyObject *ret=PyTuple_New(3); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(ret0),SWIGTYPE_p_ParaMEDMEM__DataArrayDouble, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,PyInt_FromLong(ret1)); + PyTuple_SetItem(ret,2,PyInt_FromLong(ret2)); + return ret; + } PyObject *computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, PyObject *locs) const throw(INTERP_KERNEL::Exception) { -- 2.39.2