open kriging API (2)

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 83720b7abc7bd5bcb896dc87c0974308302f195b..76344523e13a62855fb58a21f731efcf256ddb80 100644 (file)
--- 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.
  * \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
 {
  */
 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const
 {

-  if(!mesh)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients : NULL input mesh !");
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
-  int nbOfPts=coords->getNumberOfTuples();
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
-  operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
-  // Drift
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
-  matSz=nbOfPts+isDrift;
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift(computeMatrix(mesh,isDrift,matSz));
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(DataArrayDouble::New());

   matrixInv->alloc(matSz*matSz,1);
   INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),matSz,matrixInv->getPointer());
   return matrixInv.retn();
 }
 
   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<DataArrayDouble> coords(getLocalizationOfDiscValues(mesh));
+    int nbOfPts(coords->getNumberOfTuples());
+    MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix(coords->buildEuclidianDistanceDenseMatrix());
+    operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
+    // Drift
+    MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> 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.
 /*!
  * 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<DataArrayDouble> matrixInv(computeInverseMatrix(mesh,isDrift,nbRows));
 {
   int nbRows(-1);
   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(computeInverseMatrix(mesh,isDrift,nbRows));

-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=DataArrayDouble::New();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK(DataArrayDouble::New());

   KnewiK->alloc(nbRows*1,1);
   KnewiK->alloc(nbRows*1,1);

-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> 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<DataArrayDouble> arr2(PerformDriftOfVec(arr,isDrift));
+  INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbRows,nbRows,arr2->getConstPointer(),arr2->getNumberOfTuples(),1,KnewiK->getPointer());

   return KnewiK.retn();
 }
 
   return KnewiK.retn();
 }
 


@@ -2945,6 +2954,23 @@ DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftRect(const D
   return ret.retn();
 }
 
   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<DataArrayDouble> 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.
 /*!
  * 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 40afef5937719783297f05408fd49326d2dbac9c..5bb752e0a9e2ea26a8c9bac545646ad4a8fca50c 100644 (file)
--- 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;
     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);
     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;
   public:
     static const char REPR[];
     static const TypeOfField TYPE;

diff --git a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i
index dc9da45388ca3f62332dd04f59998e66a97998a5..5394a2cc2f3324e3a6ab0d4ff95994d9d58efd5a 100644 (file)
--- 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::MEDCouplingFieldDiscretization::clonePart;
 %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::getValueOnMulti;
 %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::computeTupleIdsToSelectFromCellIds;

+%newobject ParaMEDMEM::MEDCouplingFieldDiscretizationKriging::PerformDriftOfVec;

 
 namespace ParaMEDMEM
 {
 
 namespace ParaMEDMEM
 {


@@ -374,6 +375,7 @@ namespace ParaMEDMEM
   class MEDCouplingFieldDiscretizationKriging : public MEDCouplingFieldDiscretizationOnNodes
   {
   public:
   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)
     %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;
       }
         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)
       {
       
       PyObject *computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, PyObject *locs) const throw(INTERP_KERNEL::Exception)
       {