return prepareNotInterpKernelOnly();
}
-void MEDCouplingRemapper::setMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector<std::map<int,double> >& m)
+void MEDCouplingRemapper::setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector<std::map<int,double> >& m)
{
MCAuto<MEDCouplingFieldTemplate> 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<std::map<int,double> >& m)
+void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector<std::map<int,double> >& m)
{
#if __cplusplus >= 201103L
restartUsing(src,target);
}
}
_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
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<std::map<int,double> >& m);
- MEDCOUPLINGREMAPPER_EXPORT void setMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector<std::map<int,double> >& m);
+ MEDCOUPLINGREMAPPER_EXPORT void setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector<std::map<int,double> >& m);
+ MEDCOUPLINGREMAPPER_EXPORT void setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector<std::map<int,double> >& 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);
}
}
+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<std::map<int,double> >& 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<nbOfRows;i++)
+ {
+ auto& line(mCpp[i]);
+ for(auto j=indPtrCPtr[i];j<indPtrCPtr[i+1];j++)
+ {
+ line[indicesCPtr[j]]=dataCPtr[j];
+ }
+ }
+#else
+ throw INTERP_KERNEL::Exception("Breaking news : 10% off for C++11 compiler :)");
+#endif
+}
+
void convertToVectMapIntDouble(PyObject *pyobj, std::vector<std::map<int,double> >& mCpp)
{
if(!PyList_Check(pyobj))
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<std::map<int,double> > 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<DataArrayInt> indptrPtr(MEDCoupling_DataArrayInt_New__SWIG_1(indptr,NULL,NULL));
+ MCAuto<DataArrayInt> indicesPtr(MEDCoupling_DataArrayInt_New__SWIG_1(indices,NULL,NULL));
+ MCAuto<DataArrayDouble> 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<std::map<int,double> > mCpp;
convertToVectMapIntDouble(m,mCpp);
- self->setMatrix(srcMesh,targetMesh,method,mCpp);
+ self->setCrudeMatrixEx(src,target,mCpp);
}
}
};
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)):
