From 0f3787c331d2154939053507c786977c933e6ba8 Mon Sep 17 00:00:00 2001 From: ageay Date: Mon, 28 Jun 2010 12:34:25 +0000 Subject: [PATCH] Management of 1D interpolation in ParaMEDMEM. --- src/ParaMEDMEM/ElementLocator.cxx | 2 +- src/ParaMEDMEM/InterpolationMatrix.cxx | 24 ++ src/ParaMEDMEMTest/ParaMEDMEMTest.hxx | 4 + .../ParaMEDMEMTest_InterpKernelDEC.cxx | 381 +++++++++++++++++- 4 files changed, 409 insertions(+), 2 deletions(-) diff --git a/src/ParaMEDMEM/ElementLocator.cxx b/src/ParaMEDMEM/ElementLocator.cxx index e694c4e1d..29e2e26f5 100644 --- a/src/ParaMEDMEM/ElementLocator.cxx +++ b/src/ParaMEDMEM/ElementLocator.cxx @@ -35,7 +35,7 @@ using namespace std; -#define USE_DIRECTED_BB +//#define USE_DIRECTED_BB namespace ParaMEDMEM { diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index 5f3a424ea..c7db0371c 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -24,6 +24,8 @@ #include "InterpolationMatrix.hxx" #include "TranslationRotationMatrix.hxx" #include "Interpolation.hxx" +#include "Interpolation1D.txx" +#include "Interpolation2DCurve.txx" #include "Interpolation2D.txx" #include "Interpolation3DSurf.txx" #include "Interpolation3D.txx" @@ -160,6 +162,28 @@ namespace ParaMEDMEM { throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions"); } + else if( distant_support.getMeshDimension() == 1 + && distant_support.getSpaceDimension() == 1 ) + { + MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(source_supportC); + + INTERP_KERNEL::Interpolation1D interpolation(*this); + colSize=interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); + } + else if( distant_support.getMeshDimension() == 1 + && distant_support.getSpaceDimension() == 2 ) + { + MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(distant_supportC); + MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(source_supportC); + + INTERP_KERNEL::Interpolation2DCurve interpolation(*this); + colSize=interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str()); + target_wrapper.releaseTempArrays(); + source_wrapper.releaseTempArrays(); + } else if ( distant_support.getMeshDimension() == 2 && distant_support.getSpaceDimension() == 3 ) { diff --git a/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx b/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx index da4e4873c..2fe285702 100644 --- a/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx +++ b/src/ParaMEDMEMTest/ParaMEDMEMTest.hxx @@ -36,6 +36,8 @@ class ParaMEDMEMTest : public CppUnit::TestFixture CPPUNIT_TEST(testMPIProcessorGroup_rank); CPPUNIT_TEST(testBlockTopology_constructor); CPPUNIT_TEST(testBlockTopology_serialize); + CPPUNIT_TEST(testInterpKernelDEC_1D); + CPPUNIT_TEST(testInterpKernelDEC_2DCurve); CPPUNIT_TEST(testInterpKernelDEC_2D); CPPUNIT_TEST(testInterpKernelDEC2_2D); CPPUNIT_TEST(testInterpKernelDEC_2DP0P1); @@ -87,6 +89,8 @@ public: void testMPIProcessorGroup_rank(); void testBlockTopology_constructor(); void testBlockTopology_serialize(); + void testInterpKernelDEC_1D(); + void testInterpKernelDEC_2DCurve(); void testInterpKernelDEC_2D(); void testInterpKernelDEC2_2D(); void testInterpKernelDEC_2DP0P1(); diff --git a/src/ParaMEDMEMTest/ParaMEDMEMTest_InterpKernelDEC.cxx b/src/ParaMEDMEMTest/ParaMEDMEMTest_InterpKernelDEC.cxx index c41981608..513792c60 100644 --- a/src/ParaMEDMEMTest/ParaMEDMEMTest_InterpKernelDEC.cxx +++ b/src/ParaMEDMEMTest/ParaMEDMEMTest_InterpKernelDEC.cxx @@ -47,7 +47,7 @@ using namespace std; using namespace ParaMEDMEM; - + void ParaMEDMEMTest::testInterpKernelDEC_2D() { testInterpKernelDEC_2D_("P0","P0"); @@ -68,6 +68,385 @@ void ParaMEDMEMTest::testInterpKernelDEC_2DP0P1() //testInterpKernelDEC_2D_("P0","P1"); } +void ParaMEDMEMTest::testInterpKernelDEC_1D() +{ + int size; + int rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + // + if(size!=5) + return ; + int nproc_source = 3; + set self_procs; + set procs_source; + set procs_target; + + for (int i=0; icontainsMyRank()) + { + if(rank==0) + { + double coords[4]={0.3,0.7, 0.9,1.0}; + int conn[4]={0,1,2,3}; + mesh=MEDCouplingUMesh::New("Source mesh Proc0",1); + mesh->allocateCells(2); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+2); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(4,1); + std::copy(coords,coords+4,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + if(rank==1) + { + double coords[2]={0.7,0.9}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Source mesh Proc1",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,1); + std::copy(coords,coords+2,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + if(rank==2) + { + double coords[2]={1.,1.12}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Source mesh Proc2",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,1); + std::copy(coords,coords+2,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + paramesh=new ParaMESH(mesh,*source_group,"source mesh"); + ParaMEDMEM::ComponentTopology comptopo; + parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + double *valueP0=parafieldP0->getField()->getArray()->getPointer(); + parafieldP0->getField()->setNature(ConservativeVolumic); + if(rank==0) + { + valueP0[0]=7.; valueP0[1]=8.; + } + if(rank==1) + { + valueP0[0]=9.; + } + if(rank==2) + { + valueP0[0]=10.; + } + } + else + { + const char targetMeshName[]="target mesh"; + if(rank==3) + { + double coords[2]={0.5,0.75}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Target mesh Proc3",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,1); + std::copy(coords,coords+2,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + paramesh=new ParaMESH(mesh,*target_group,targetMeshName); + } + if(rank==4) + { + double coords[2]={0.75,1.2}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Target mesh Proc4",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,1); + std::copy(coords,coords+2,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + paramesh=new ParaMESH(mesh,*target_group,targetMeshName); + } + ParaMEDMEM::ComponentTopology comptopo; + parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafieldP0->getField()->setNature(ConservativeVolumic); + } + // test 1 + ParaMEDMEM::InterpKernelDEC dec(*source_group,*target_group); + if (source_group->containsMyRank()) + { + dec.setMethod("P0"); + dec.attachLocalField(parafieldP0); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.sendData(); + dec.recvData(); + const double *valueP0=parafieldP0->getField()->getArray()->getPointer(); + if(rank==0) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(7.4,valueP0[0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0540540540540544,valueP0[1],1e-7); + } + if(rank==1) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(8.64054054054054,valueP0[0],1e-7); + } + if(rank==2) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0540540540540544,valueP0[0],1e-7); + } + } + else + { + dec.setMethod("P0"); + dec.attachLocalField(parafieldP0); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.recvData(); + const double *res=parafieldP0->getField()->getArray()->getConstPointer(); + if(rank==3) + { + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfComponents()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(7.4,res[0],1e-12); + } + if(rank==4) + { + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfComponents()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0540540540540526,res[0],1e-12); + } + dec.sendData(); + } + // + delete parafieldP0; + mesh->decrRef(); + delete paramesh; + delete self_group; + delete target_group; + delete source_group; + // + MPI_Barrier(MPI_COMM_WORLD); +} + +void ParaMEDMEMTest::testInterpKernelDEC_2DCurve() +{ + int size; + int rank; + MPI_Comm_size(MPI_COMM_WORLD,&size); + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + // + if(size!=5) + return ; + int nproc_source = 3; + set self_procs; + set procs_source; + set procs_target; + + for (int i=0; icontainsMyRank()) + { + if(rank==0) + { + double coords[8]={0.3,0.3,0.7,0.7, 0.9,0.9,1.0,1.0}; + int conn[4]={0,1,2,3}; + mesh=MEDCouplingUMesh::New("Source mesh Proc0",1); + mesh->allocateCells(2); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+2); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(4,2); + std::copy(coords,coords+8,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + if(rank==1) + { + double coords[4]={0.7,0.7,0.9,0.9}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Source mesh Proc1",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,2); + std::copy(coords,coords+4,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + if(rank==2) + { + double coords[4]={1.,1.,1.12,1.12}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Source mesh Proc2",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,2); + std::copy(coords,coords+4,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + } + paramesh=new ParaMESH(mesh,*source_group,"source mesh"); + ParaMEDMEM::ComponentTopology comptopo; + parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + double *valueP0=parafieldP0->getField()->getArray()->getPointer(); + parafieldP0->getField()->setNature(ConservativeVolumic); + if(rank==0) + { + valueP0[0]=7.; valueP0[1]=8.; + } + if(rank==1) + { + valueP0[0]=9.; + } + if(rank==2) + { + valueP0[0]=10.; + } + } + else + { + const char targetMeshName[]="target mesh"; + if(rank==3) + { + double coords[4]={0.5,0.5,0.75,0.75}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Target mesh Proc3",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,2); + std::copy(coords,coords+4,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + paramesh=new ParaMESH(mesh,*target_group,targetMeshName); + } + if(rank==4) + { + double coords[4]={0.75,0.75,1.2,1.2}; + int conn[2]={0,1}; + mesh=MEDCouplingUMesh::New("Target mesh Proc4",1); + mesh->allocateCells(1); + mesh->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(2,2); + std::copy(coords,coords+4,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + paramesh=new ParaMESH(mesh,*target_group,targetMeshName); + } + ParaMEDMEM::ComponentTopology comptopo; + parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo); + parafieldP0->getField()->setNature(ConservativeVolumic); + } + // test 1 + ParaMEDMEM::InterpKernelDEC dec(*source_group,*target_group); + if (source_group->containsMyRank()) + { + dec.setMethod("P0"); + dec.attachLocalField(parafieldP0); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.sendData(); + dec.recvData(); + const double *valueP0=parafieldP0->getField()->getArray()->getPointer(); + if(rank==0) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(7.4,valueP0[0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0540540540540544,valueP0[1],1e-7); + } + if(rank==1) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(8.64054054054054,valueP0[0],1e-7); + } + if(rank==2) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0540540540540544,valueP0[0],1e-7); + } + } + else + { + dec.setMethod("P0"); + dec.attachLocalField(parafieldP0); + dec.synchronize(); + dec.setForcedRenormalization(false); + dec.recvData(); + const double *res=parafieldP0->getField()->getArray()->getConstPointer(); + if(rank==3) + { + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfComponents()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(7.4,res[0],1e-12); + } + if(rank==4) + { + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,parafieldP0->getField()->getNumberOfComponents()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(9.0540540540540526,res[0],1e-12); + } + dec.sendData(); + } + // + delete parafieldP0; + mesh->decrRef(); + delete paramesh; + delete self_group; + delete target_group; + delete source_group; + // + MPI_Barrier(MPI_COMM_WORLD); +} + + /* * Check methods defined in InterpKernelDEC.hxx * -- 2.39.2