From 67310dbac78707d934fcdc4d9f51b70cf6d3e6cd Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Mon, 26 Mar 2018 15:23:26 +0200 Subject: [PATCH] MEDCouplingRemapper::SetCrudeMatrix implementation --- src/MEDCoupling/MEDCouplingRemapper.cxx | 10 ++- src/MEDCoupling/MEDCouplingRemapper.hxx | 4 +- .../MEDCouplingDataArrayTypemaps.i | 42 ++++++++++ .../MEDCouplingRemapperCommon.i | 27 ++++++- .../MEDCouplingRemapperTest.py | 79 +++++++++++++++++++ 5 files changed, 155 insertions(+), 7 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index f7c3924d4..1e8143833 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -94,14 +94,14 @@ int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const ME return prepareNotInterpKernelOnly(); } -void MEDCouplingRemapper::setMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector >& m) +void MEDCouplingRemapper::setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector >& m) { MCAuto src,target; BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target); - setMatrixEx(src,target,m); + setCrudeMatrixEx(src,target,m); } -void MEDCouplingRemapper::setMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector >& m) +void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector >& m) { #if __cplusplus >= 201103L restartUsing(src,target); @@ -124,6 +124,10 @@ void MEDCouplingRemapper::setMatrixEx(const MEDCouplingFieldTemplate *src, const } } _matrix=m; + _deno_multiply.clear(); + _deno_multiply.resize(_matrix.size()); + _deno_reverse_multiply.clear(); + _deno_reverse_multiply.resize(srcNbElem); #else throw INTERP_KERNEL::Exception("Breaking news : 10% off for C++11 compiler :)"); #endif diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx index afb75bb90..5d7ac9da8 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.hxx +++ b/src/MEDCoupling/MEDCouplingRemapper.hxx @@ -56,8 +56,8 @@ namespace MEDCoupling MEDCOUPLINGREMAPPER_EXPORT ~MEDCouplingRemapper(); MEDCOUPLINGREMAPPER_EXPORT int prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method); MEDCOUPLINGREMAPPER_EXPORT int prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target); - MEDCOUPLINGREMAPPER_EXPORT void setMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector >& m); - MEDCOUPLINGREMAPPER_EXPORT void setMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector >& m); + MEDCOUPLINGREMAPPER_EXPORT void setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector >& m); + MEDCOUPLINGREMAPPER_EXPORT void setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector >& m); MEDCOUPLINGREMAPPER_EXPORT void transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue); MEDCOUPLINGREMAPPER_EXPORT void partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField); MEDCOUPLINGREMAPPER_EXPORT void reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue); diff --git a/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i b/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i index 5f44c75ca..6fd671e1b 100644 --- a/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i +++ b/src/MEDCoupling_Swig/MEDCouplingDataArrayTypemaps.i @@ -3032,6 +3032,48 @@ PyObject *DataArrayT__getitem__internal(const typename MEDCoupling::Traits::A } } +bool isCSRMatrix(PyObject *m) +{ +#if defined(WITH_NUMPY) && defined(WITH_SCIPY) + PyObject* pdict(PyDict_New()); + PyDict_SetItemString(pdict, "__builtins__", PyEval_GetBuiltins()); + PyObject *tmp(PyRun_String("from scipy.sparse import csr_matrix", Py_single_input, pdict, pdict)); + if(!tmp) + throw INTERP_KERNEL::Exception("Problem during loading csr_matrix in scipy.sparse ! Is Scipy module available in present ?"); + PyObject *csrMatrixCls=PyDict_GetItemString(pdict,"csr_matrix"); + if(!csrMatrixCls) + throw INTERP_KERNEL::Exception("csr_matrix not found in scipy.sparse ! Is Scipy module available in present ?"); + bool ret(PyObject_IsInstance(m,csrMatrixCls)); + Py_DECREF(pdict); Py_XDECREF(tmp); + return ret; +#else + return false; +#endif +} + +void convertCSR_MCDataToVectMapIntDouble(const MEDCoupling::DataArrayInt *indptrPtr, const MEDCoupling::DataArrayInt *indicesPtr, const MEDCoupling::DataArrayDouble *dataPtr, std::vector >& mCpp) +{ +#if __cplusplus >= 201103L + auto nbOfRows(indptrPtr->getNumberOfTuples()-1); + if(nbOfRows<0) + throw INTERP_KERNEL::Exception("pywrap of MEDCouplingRemapper::setMatrix : input CSR matrix looks bad regarding indptr array !"); + mCpp.resize(nbOfRows); + auto indPtrCPtr(indptrPtr->begin()); + auto indicesCPtr(indicesPtr->begin()); + auto dataCPtr(dataPtr->begin()); + for(auto i=0;i >& mCpp) { if(!PyList_Check(pyobj)) diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i b/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i index a88c0b468..7f3ce4239 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperCommon.i @@ -88,11 +88,34 @@ namespace MEDCoupling return ToCSRMatrix(self->getCrudeMatrix(),self->getNumberOfColsOfMatrix()); } #endif - void setMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, PyObject *m) throw(INTERP_KERNEL::Exception) + void setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, PyObject *m) throw(INTERP_KERNEL::Exception) + { + std::vector > mCpp; + if(isCSRMatrix(m)) + { +#if defined(WITH_NUMPY) && defined(WITH_SCIPY) + PyObject *indptr(PyObject_GetAttrString(m,"indptr")); + PyObject *indices(PyObject_GetAttrString(m,"indices")); + PyObject *data(PyObject_GetAttrString(m,"data")); + MCAuto indptrPtr(MEDCoupling_DataArrayInt_New__SWIG_1(indptr,NULL,NULL)); + MCAuto indicesPtr(MEDCoupling_DataArrayInt_New__SWIG_1(indices,NULL,NULL)); + MCAuto dataPtr(MEDCoupling_DataArrayDouble_New__SWIG_1(data,NULL,NULL)); + convertCSR_MCDataToVectMapIntDouble(indptrPtr,indicesPtr,dataPtr,mCpp); + Py_XDECREF(data); Py_XDECREF(indptr); Py_XDECREF(indices); +#else + throw INTERP_KERNEL::Exception("pywrap of MEDCouplingRemapper::setCrudeMatrix : unexpected situation regarding numpy/scipy !"); +#endif + } + else + convertToVectMapIntDouble(m,mCpp); + self->setCrudeMatrix(srcMesh,targetMesh,method,mCpp); + } + + void setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, PyObject *m) throw(INTERP_KERNEL::Exception) { std::vector > mCpp; convertToVectMapIntDouble(m,mCpp); - self->setMatrix(srcMesh,targetMesh,method,mCpp); + self->setCrudeMatrixEx(src,target,mCpp); } } }; diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index dc0d1ea79..e5c50edb1 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -1200,6 +1200,85 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertAlmostEqual(delta2.sum(),0.,14) pass + def testSetMatrix1(self): + """ Remapper has now setCrudeMatrix method to reload matrix to skip prepare phase """ + cooS=DataArrayDouble([1,1, 7,1, 7,2, 1,2],4,2) + cooT=DataArrayDouble([0,0, 3,0, 3,3, 0,3, 6,0, 12,0, 12,3, 6,3],8,2) + ms=MEDCouplingUMesh("source",2) ; ms.allocateCells(1) ; ms.insertNextCell(NORM_QUAD4,[0,1,2,3]) ; ms.setCoords(cooS) + mt=MEDCouplingUMesh("target",2) ; mt.allocateCells(2) ; mt.insertNextCell(NORM_QUAD4,[0,1,2,3]) ; mt.insertNextCell(NORM_QUAD4,[4,5,6,7]) ; mt.setCoords(cooT) + rem=MEDCouplingRemapper() + self.assertEqual(rem.prepare(ms,mt,"P0P0"),1) # [{0: 2.0}, {0: 1.0}] + fs=MEDCouplingFieldDouble(ON_CELLS) + fs.setMesh(ms) + fs.setArray(DataArrayDouble([10])) + fs.checkConsistencyLight() + # + fs.setNature(ExtensiveConservation) + self.assertTrue(rem.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([20./3,10./3.]),1e-12))# sum is equal to 10. First value is twice than second value + # + fs.setNature(ExtensiveMaximum) + self.assertTrue(rem.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([20./6.,10./6.]),1e-12))#sum is equal to 5 (10/2. because only half part on input cell is intercepted by the target cells). First value is twice than second value + # + fs.setNature(IntensiveConservation) + self.assertTrue(rem.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([2./9.*10.,1./18.*10.]),1e-12))# + # + fs.setNature(IntensiveMaximum) + self.assertTrue(rem.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([10.,10.]),1e-12))# + #### + rem2=MEDCouplingRemapper() + rem2.setCrudeMatrix(ms,mt,"P0P0",rem.getCrudeMatrix()) + fs.setNature(ExtensiveConservation) + self.assertTrue(rem2.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([20./3,10./3.]),1e-12)) + # + fs.setNature(ExtensiveMaximum) + self.assertTrue(rem2.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([20./6.,10./6.]),1e-12)) + # + fs.setNature(IntensiveConservation) + self.assertTrue(rem2.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([2./9.*10.,1./18.*10.]),1e-12)) + # + fs.setNature(IntensiveMaximum) + self.assertTrue(rem2.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([10.,10.]),1e-12)) + # + srcFt=MEDCouplingFieldTemplate.New(ON_CELLS); + trgFt=MEDCouplingFieldTemplate.New(ON_CELLS); + srcFt.setMesh(ms); + trgFt.setMesh(mt); + rem3=MEDCouplingRemapper() + rem3.setCrudeMatrixEx(srcFt,trgFt,rem.getCrudeMatrix()) + fs.setNature(ExtensiveConservation) + self.assertTrue(rem3.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([20./3,10./3.]),1e-12)) + pass + + @unittest.skipUnless(MEDCouplingHasNumPyBindings() and MEDCouplingHasSciPyBindings(),"requires numpy AND scipy") + def testSetMatrix2(self): + """ Remapper has now setCrudeMatrix method to reload matrix to skip prepare phase. Same as testSetMatrix1 but with CSR scipy matrix """ + arrx_s=DataArrayDouble(6) ; arrx_s.iota() + arry_s=DataArrayDouble(6) ; arry_s.iota() + ms=MEDCouplingCMesh() ; ms.setCoords(arrx_s,arry_s) + ms=ms.buildUnstructured() + # + arrx_t=DataArrayDouble([2.5,4.5,5.5]) + arry_t=DataArrayDouble([2.5,3.5,5.5]) + mt=MEDCouplingCMesh() ; mt.setCoords(arrx_t,arry_t) + mt=mt.buildUnstructured() + # + rem=MEDCouplingRemapper() + self.assertEqual(rem.prepare(ms,mt,"P0P0"),1) + # + fs=MEDCouplingFieldDouble(ON_CELLS) + fs.setMesh(ms) + arr=DataArrayDouble(25) ; arr.iota() + fs.setArray(arr) + fs.checkConsistencyLight() + # + fs.setNature(ExtensiveConservation) + self.assertTrue(rem.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([54.25,11.75,79.25,16.75]),1e-12)) + mat=rem.getCrudeCSRMatrix() + rem2=MEDCouplingRemapper() + rem2.setCrudeMatrix(ms,mt,"P0P0",mat) + self.assertTrue(rem2.transferField(fs,1e300).getArray().isEqual(DataArrayDouble([54.25,11.75,79.25,16.75]),1e-12)) + pass + def checkMatrix(self,mat1,mat2,nbCols,eps): self.assertEqual(len(mat1),len(mat2)) for i in range(len(mat1)): -- 2.39.2