From 0402b119ef271b422db415ddad33365fe82acce2 Mon Sep 17 00:00:00 2001 From: geay Date: Wed, 26 Feb 2014 12:37:30 +0100 Subject: [PATCH] debug for 3D1D interpolation in P0P1 --- src/MEDCoupling/MEDCouplingRemapper.cxx | 21 +++++++--- src/MEDCoupling/MEDCouplingRemapper.hxx | 1 + src/MEDCoupling_Swig/MEDCouplingRemapper.i | 1 + .../MEDCouplingRemapperTest.py | 38 +++++++++++++++++++ 4 files changed, 55 insertions(+), 6 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index 8700012b0..3327043a4 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -304,7 +304,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() const MEDCouplingPointSet *src_mesh=static_cast(_src_ft->getMesh()); const MEDCouplingPointSet *target_mesh=static_cast(_target_ft->getMesh()); std::string srcMeth,trgMeth; - std::string method=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth); + std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth)); const int srcMeshDim=src_mesh->getMeshDimension(); int srcSpaceDim=-1; if(srcMeshDim!=-1) @@ -369,7 +369,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation3D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); } @@ -388,7 +389,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation2D1D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces(); @@ -412,7 +414,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation2D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); } @@ -459,7 +462,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation3D2D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces(); @@ -869,7 +873,12 @@ std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !"); srcMeth=_src_ft->getDiscretization()->getRepr(); trgMeth=_target_ft->getDiscretization()->getRepr(); - std::string method(srcMeth); method+=trgMeth; + return BuildMethodFrom(srcMeth,trgMeth); +} + +std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2) +{ + std::string method(meth1); method+=meth2; return method; } diff --git a/src/MEDCoupling/MEDCouplingRemapper.hxx b/src/MEDCoupling/MEDCouplingRemapper.hxx index 55df5ab9b..5325745b9 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.hxx +++ b/src/MEDCoupling/MEDCouplingRemapper.hxx @@ -74,6 +74,7 @@ namespace ParaMEDMEM MEDCOUPLINGREMAPPER_EXPORT const std::vector >& getCrudeMatrix() const; MEDCOUPLINGREMAPPER_EXPORT int getNumberOfColsOfMatrix() const; MEDCOUPLINGREMAPPER_EXPORT static void PrintMatrix(const std::vector >& m); + MEDCOUPLINGREMAPPER_EXPORT static std::string BuildMethodFrom(const std::string& meth1, const std::string& meth2); private: int prepareInterpKernelOnly(); int prepareInterpKernelOnlyUU(); diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapper.i b/src/MEDCoupling_Swig/MEDCouplingRemapper.i index c77a2d09f..5154fc43e 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapper.i +++ b/src/MEDCoupling_Swig/MEDCouplingRemapper.i @@ -74,6 +74,7 @@ namespace ParaMEDMEM int nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception); double getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception); int getNumberOfColsOfMatrix() const throw(INTERP_KERNEL::Exception); + static std::string BuildMethodFrom(const std::string& meth1, const std::string& meth2) throw(INTERP_KERNEL::Exception); %extend { PyObject *getCrudeMatrix() const throw(INTERP_KERNEL::Exception) diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index c65f85d99..ab4e726b5 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -722,6 +722,44 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(cellIdsOfSkin.isEqual(DataArrayInt([1,2,3,5,6,7,8,9,10,11,12,13,14,15,16,17,19,20,21,23]))) self.assertTrue(cellIdsOfNonConformCells.isEqual(DataArrayInt([0,4,18,22]))) pass + + def test3D1DOnP1P0_1(self): + """ This test focused on P1P0 interpolation with a source with meshDim=1 spaceDim=3 and a target with meshDim=3. + This test has revealed a bug in remapper. A reverse matrix is computed so a reverse method should be given in input. + """ + target=MEDCouplingCMesh() + arrX=DataArrayDouble([0,1]) ; arrY=DataArrayDouble([0,1]) ; arrZ=DataArrayDouble(11) ; arrZ.iota() + target.setCoords(arrX,arrY,arrZ) + target=target.buildUnstructured() ; target.setName("TargetSecondaire") + # + sourceCoo=DataArrayDouble([(0.5,0.5,0.1),(0.5,0.5,1.2),(0.5,0.5,1.6),(0.5,0.5,1.8),(0.5,0.5,2.43),(0.5,0.5,2.55),(0.5,0.5,4.1),(0.5,0.5,4.4),(0.5,0.5,4.9),(0.5,0.5,5.1),(0.5,0.5,7.6),(0.5,0.5,7.7),(0.5,0.5,8.2),(0.5,0.5,8.4),(0.5,0.5,8.6),(0.5,0.5,8.8),(0.5,0.5,9.2),(0.5,0.5,9.6),(0.5,0.5,11.5)]) + source=MEDCoupling1SGTUMesh("SourcePrimaire",NORM_SEG2) + source.setCoords(sourceCoo) + source.allocateCells() + for i in xrange(len(sourceCoo)-1): + source.insertNextCell([i,i+1]) + pass + source=source.buildUnstructured() + fsource=MEDCouplingFieldDouble(ON_NODES) ; fsource.setName("field") + fsource.setMesh(source) + arr=DataArrayDouble(len(sourceCoo)) ; arr.iota(0.7) ; arr*=arr + fsource.setArray(arr) + fsource.setNature(ConservativeVolumic) + # + rem=MEDCouplingRemapper() + rem.setIntersectionType(PointLocator) + rem.prepare(source,target,"P1P0") + f2Test=rem.transferField(fsource,-27) + self.assertEqual(f2Test.getName(),fsource.getName()) + self.assertEqual(f2Test.getMesh().getHiddenCppPointer(),target.getHiddenCppPointer()) + expArr=DataArrayDouble([0.49,7.956666666666667,27.29,-27,59.95666666666667,94.09,-27,125.69,202.89,296.09]) + self.assertTrue(f2Test.getArray().isEqual(expArr,1e-12)) + f2Test=rem.reverseTransferField(f2Test,-36) + self.assertEqual(f2Test.getName(),fsource.getName()) + self.assertEqual(f2Test.getMesh().getHiddenCppPointer(),source.getHiddenCppPointer()) + expArr2=DataArrayDouble([0.49,7.956666666666667,7.956666666666667,7.956666666666667,27.29,27.29,59.95666666666667,59.95666666666667,59.95666666666667,94.09,125.69,125.69,202.89,202.89,202.89,202.89,296.09,296.09,-36.]) + self.assertTrue(f2Test.getArray().isEqual(expArr2,1e-12)) + pass def build2DSourceMesh_1(self): sourceCoords=[-0.3,-0.3, 0.7,-0.3, -0.3,0.7, 0.7,0.7] -- 2.39.2