From fc6c67329ec7967143b8e16efb66dcc6f7a6bf61 Mon Sep 17 00:00:00 2001 From: ageay Date: Thu, 29 Sep 2011 16:08:08 +0000 Subject: [PATCH] Introduction of 2D/1D interpolation into MEDCouplingRemapper. --- src/MEDCoupling/MEDCouplingRemapper.cxx | 50 ++++++++++----- .../Test/MEDCouplingRemapperTest.cxx | 62 ++++++++++++++++++- 2 files changed, 96 insertions(+), 16 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index 56d8c819a..8c9f3ee40 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -30,6 +30,7 @@ #include "Interpolation2D.txx" #include "Interpolation3D.txx" #include "Interpolation3DSurf.hxx" +#include "Interpolation2D1D.txx" using namespace ParaMEDMEM; @@ -234,24 +235,43 @@ int MEDCouplingRemapper::prepareUU(const char *method) throw(INTERP_KERNEL::Exce } else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2) { - if(getIntersectionType()!=INTERP_KERNEL::PointLocator) - throw INTERP_KERNEL::Exception("Invalid interpolation requested between 2D and 1D ! Select PointLocator as intersection type !"); - MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh); - MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh); - INTERP_KERNEL::Interpolation2D interpolation(*this); - nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method); + if(getIntersectionType()==INTERP_KERNEL::PointLocator) + { + MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh); + MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh); + INTERP_KERNEL::Interpolation2D interpolation(*this); + nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method); + } + else + { + MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh); + 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); + reverseMatrix(matrixTmp,nbCols,_matrix); + nbCols=matrixTmp.size(); + } } else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2) { - if(getIntersectionType()!=INTERP_KERNEL::PointLocator) - throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 2D ! Select PointLocator as intersection type !"); - MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh); - 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); - reverseMatrix(matrixTmp,nbCols,_matrix); - nbCols=matrixTmp.size(); + if(getIntersectionType()==INTERP_KERNEL::PointLocator) + { + MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh); + 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); + reverseMatrix(matrixTmp,nbCols,_matrix); + nbCols=matrixTmp.size(); + } + else + { + MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh); + MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh); + INTERP_KERNEL::Interpolation2D1D interpolation(*this); + nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method); + } } else if(trgMeshDim==-1) { diff --git a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx index 930b561ea..58d2de4f6 100644 --- a/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingRemapperTest.cxx @@ -622,7 +622,67 @@ void MEDCouplingRemapperTest::testMultiDimCombi() values=srcField->getArray()->getConstPointer(); for(int i0=0;i0<5;i0++) CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesExpected7[i0],values[i0],1e-14); - // + srcField->decrRef(); + trgField->decrRef(); + sourceMesh->decrRef(); + targetMesh->decrRef(); + //------------- 1D -> 2D + const int conn[8]={0,1,1,2,2,3,3,0}; + const int conn2[12]={6,7,5,4,2,7,6,3,0,4,5,1}; + const double coords1[]={0.17,0.93,0.56,0.93,0.56,0.25,0.17,0.52}; + const double coords2[]={0.,0.,1.,0.,1.,1.,0.,1.,0.,0.5,1.,0.5,0.,0.8,1.,0.8}; + sourceMesh=MEDCouplingUMesh::New("src1D",1); + sourceMesh->allocateCells(4); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+2); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+4); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+6); + sourceMesh->finishInsertingCells(); + array=DataArrayDouble::New(); array->alloc(4,2); + std::copy(coords1,coords1+8,array->getPointer()); + sourceMesh->setCoords(array); array->decrRef(); + targetMesh=MEDCouplingUMesh::New("trg2D",2); + targetMesh->allocateCells(3); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn2); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn2+4); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn2+8); + targetMesh->finishInsertingCells(); + array=DataArrayDouble::New(); array->alloc(8,2); + std::copy(coords2,coords2+16,array->getPointer()); + targetMesh->setCoords(array); array->decrRef(); + remapper.setPrecision(1e-12); + remapper.setIntersectionType(INTERP_KERNEL::Geometric2D); + CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(sourceMesh,targetMesh,"P0P0")); + srcField=MEDCouplingFieldDouble::New(ON_CELLS); + srcField->setNature(ConservativeVolumic); + srcField->setMesh(sourceMesh); + array=DataArrayDouble::New(); + array->alloc(4,1); array->iota(2.); + srcField->setArray(array); array->decrRef(); + trgField=remapper.transferField(srcField,4.57); + const double valuesExpected10[3]={3.9674868868103834, 2.8, 3.6372633449255796}; + CPPUNIT_ASSERT_EQUAL(3,trgField->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,trgField->getNumberOfComponents()); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesExpected10[i],trgField->getIJ(i,0),1e-13); + srcField->decrRef(); + trgField->decrRef(); + //------------- 2D -> 1D + remapper.setPrecision(1e-12); + remapper.setIntersectionType(INTERP_KERNEL::Geometric2D); + CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(targetMesh,sourceMesh,"P0P0")); + srcField=MEDCouplingFieldDouble::New(ON_CELLS); + srcField->setNature(ConservativeVolumic); + srcField->setMesh(targetMesh); + array=DataArrayDouble::New(); + array->alloc(3,1); array->iota(2.); + srcField->setArray(array); array->decrRef(); + trgField=remapper.transferField(srcField,4.57); + const double valuesExpected11[4]={3., 2.9264705882352944, 3.8518518518518516, 2.3170731707317076}; + CPPUNIT_ASSERT_EQUAL(4,trgField->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,trgField->getNumberOfComponents()); + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesExpected11[i],trgField->getIJ(i,0),1e-13); srcField->decrRef(); trgField->decrRef(); sourceMesh->decrRef(); -- 2.39.2