From 4c76b05991a6ab0ae3f08298869b9d2195e89f25 Mon Sep 17 00:00:00 2001
From: geay <anthony.geay@cea.fr>
Date: Thu, 24 Apr 2014 15:36:06 +0200
Subject: [PATCH] open kriging API

---
 .../MEDCouplingFieldDiscretization.cxx        | 64 +++++++++++++++----
 .../MEDCouplingFieldDiscretization.hxx        |  8 ++-
 .../MEDCouplingFieldDiscretization.i          | 41 ++++++++++++
 3 files changed, 99 insertions(+), 14 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 08aa48cdf..83720b7ab 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -2877,21 +2877,12 @@ void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimens
   {
     case 1:
       {
-        for(int i=0;i<nbOfElems;i++)
-          {
-            double val=matrixPtr[i];
-            matrixPtr[i]=val*val*val;
-          }
+        OperateOnDenseMatrixH3(nbOfElems,matrixPtr);
         break;
       }
     case 2:
       {
-        for(int i=0;i<nbOfElems;i++)
-          {
-            double val=matrixPtr[i];
-            if(val!=0.)
-              matrixPtr[i]=val*val*log(val);
-          }
+        OperateOnDenseMatrixH2Ln(nbOfElems,matrixPtr);
         break;
       }
     case 3:
@@ -2904,6 +2895,56 @@ void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimens
   }
 }
 
+void MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH3(int nbOfElems, double *matrixPtr)
+{
+  for(int i=0;i<nbOfElems;i++)
+    {
+      double val=matrixPtr[i];
+      matrixPtr[i]=val*val*val;
+    }
+}
+
+void MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH2Ln(int nbOfElems, double *matrixPtr)
+{
+  for(int i=0;i<nbOfElems;i++)
+    {
+      double val=matrixPtr[i];
+      if(val!=0.)
+        matrixPtr[i]=val*val*log(val);
+    }
+}
+
+/*!
+ * Performs a drift to the rectangular input matrix \a matr.
+ * This method generate a dense matrix starting from an input dense matrix \a matr and input array \a arr.
+ * \param [in] matr The rectangular dense matrix (with only one component). The number of rows of \a matr must be equal to the number of tuples of \a arr
+ * \param [in] arr The array of coords to be appended in the input dense matrix \a matr. Typically arr is an array of coordinates.
+ * \param [out] delta the delta of number of columns between returned dense matrix and input dense matrix \a matr. \a delta is equal to number of components of \a arr + 1.
+ * \sa performDrift
+ */
+DataArrayDouble *MEDCouplingFieldDiscretizationKriging::PerformDriftRect(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta)
+{
+  if(!matr || !matr->isAllocated() || matr->getNumberOfComponents()!=1)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftRect : invalid input dense matrix ! Must be allocated not NULL and with exactly one component !");
+  if(!arr || !arr->isAllocated())
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftRect : invalid input array of coordiantes ! Must be allocated and not NULL !");
+  int spaceDimension(arr->getNumberOfComponents()),nbOfPts(arr->getNumberOfTuples()),nbOfEltInMatrx(matr->getNumberOfTuples());
+  delta=spaceDimension+1;
+  int nbOfCols(nbOfEltInMatrx/nbOfPts);
+  if(nbOfEltInMatrx%nbOfPts!=0)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::PerformDriftRect : size of input dense matrix and input arrays mismatch ! NbOfElems in matrix % nb of tuples in array must be equal to 0 !");
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfPts*(nbOfCols+delta));
+  double *retPtr(ret->getPointer());
+  const double *mPtr(matr->begin()),*aPtr(arr->begin());
+  for(int i=0;i<nbOfPts;i++,aPtr+=spaceDimension,mPtr+=nbOfCols)
+    {
+      retPtr=std::copy(mPtr,mPtr+nbOfCols,retPtr);
+      *retPtr++=1.;
+      retPtr=std::copy(aPtr,aPtr+spaceDimension,retPtr);
+    }
+  return ret.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.
@@ -2913,6 +2954,7 @@ void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimens
  * \param [in] matr input matrix whose drift part will be added
  * \param [out] delta the difference between the size of the output matrix and the input matrix \a matr.
  * \return a newly allocated matrix bigger than input matrix \a matr.
+ * \sa MEDCouplingFieldDiscretizationKriging::PerformDriftRect
  */
 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const
 {
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
index bcb1a809c..40afef593 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
@@ -405,13 +405,15 @@ 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 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;
-  protected:
-    void operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const;
-    DataArrayDouble *performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const;
   public:
     static const char REPR[];
     static const TypeOfField TYPE;
diff --git a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i
index 5a7dee310..dc9da4538 100644
--- a/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i
+++ b/src/MEDCoupling_Swig/MEDCouplingFieldDiscretization.i
@@ -413,6 +413,47 @@ namespace ParaMEDMEM
         PyTuple_SetItem(ret,1,PyInt_FromLong(ret1));
         return ret;
       }
+
+      void operateOnDenseMatrix(int spaceDimension, DataArrayDouble *myMatrix) const throw(INTERP_KERNEL::Exception)
+      {
+        if(!myMatrix || !myMatrix->isAllocated() || myMatrix->getNumberOfComponents()!=1)
+          throw INTERP_KERNEL::Exception("Wrap of MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix : invalid input matrix as DataArrayDouble ! Must be allocated with one component !");
+        self->operateOnDenseMatrix(spaceDimension,myMatrix->getNumberOfTuples(),myMatrix->getPointer());
+      }
+
+      PyObject *performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr) const throw(INTERP_KERNEL::Exception)
+      {
+        int ret1(-1);
+        DataArrayDouble *ret0(self->performDrift(matr,arr,ret1));
+        PyObject *res(PyTuple_New(2));
+        PyTuple_SetItem(res,0,SWIG_NewPointerObj((void*)ret0,SWIGTYPE_p_ParaMEDMEM__DataArrayDouble,SWIG_POINTER_OWN | 0));
+        PyTuple_SetItem(res,1,PyInt_FromLong(ret1));
+        return res;
+      }
+
+      static PyObject *PerformDriftRect(const DataArrayDouble *matr, const DataArrayDouble *arr) throw(INTERP_KERNEL::Exception)
+      {
+        int ret1(-1);
+        DataArrayDouble *ret0(MEDCouplingFieldDiscretizationKriging::PerformDriftRect(matr,arr,ret1));
+        PyObject *res(PyTuple_New(2));
+        PyTuple_SetItem(res,0,SWIG_NewPointerObj((void*)ret0,SWIGTYPE_p_ParaMEDMEM__DataArrayDouble,SWIG_POINTER_OWN | 0));
+        PyTuple_SetItem(res,1,PyInt_FromLong(ret1));
+        return res;
+      }
+
+      static void OperateOnDenseMatrixH3(DataArrayDouble *myMatrix) throw(INTERP_KERNEL::Exception)
+      {
+        if(!myMatrix || !myMatrix->isAllocated() || myMatrix->getNumberOfComponents()!=1)
+          throw INTERP_KERNEL::Exception("Wrap of MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH3 : invalid input matrix as DataArrayDouble ! Must be allocated with one component !");
+        MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH3(myMatrix->getNumberOfTuples(),myMatrix->getPointer());
+      }
+      
+      static void OperateOnDenseMatrixH2Ln(DataArrayDouble *myMatrix) throw(INTERP_KERNEL::Exception)
+      {
+        if(!myMatrix || !myMatrix->isAllocated() || myMatrix->getNumberOfComponents()!=1)
+          throw INTERP_KERNEL::Exception("Wrap of MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH2Ln : invalid input matrix as DataArrayDouble ! Must be allocated with one component !");
+        MEDCouplingFieldDiscretizationKriging::OperateOnDenseMatrixH2Ln(myMatrix->getNumberOfTuples(),myMatrix->getPointer());
+      }
     }
   };
 }
-- 
2.39.2