{
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:
}
}
+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.
* \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
{
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;
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());
+ }
}
};
}
