From: ageay <ageay>
Date: Mon, 28 Jun 2010 12:34:25 +0000 (+0000)
Subject: Management of 1D interpolation in ParaMEDMEM.
X-Git-Tag: V5_1_main_FINAL~111
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=0f3787c331d2154939053507c786977c933e6ba8;p=tools%2Fmedcoupling.git

Management of 1D interpolation in ParaMEDMEM.
---

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<int> self_procs;
+  set<int> procs_source;
+  set<int> procs_target;
+  
+  for (int i=0; i<nproc_source; i++)
+    procs_source.insert(i);
+  for (int i=nproc_source; i<size; i++)
+    procs_target.insert(i);
+  self_procs.insert(rank);
+  //
+  ParaMEDMEM::MEDCouplingUMesh *mesh=0;
+  ParaMEDMEM::ParaMESH *paramesh=0;
+  ParaMEDMEM::ParaFIELD *parafieldP0=0;
+  //
+  ParaMEDMEM::CommInterface interface;
+  //
+  ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
+  ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
+  ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
+  //
+  MPI_Barrier(MPI_COMM_WORLD);
+  if(source_group->containsMyRank())
+    {
+      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<int> self_procs;
+  set<int> procs_source;
+  set<int> procs_target;
+  
+  for (int i=0; i<nproc_source; i++)
+    procs_source.insert(i);
+  for (int i=nproc_source; i<size; i++)
+    procs_target.insert(i);
+  self_procs.insert(rank);
+  //
+  ParaMEDMEM::MEDCouplingUMesh *mesh=0;
+  ParaMEDMEM::ParaMESH *paramesh=0;
+  ParaMEDMEM::ParaFIELD *parafieldP0=0;
+  //
+  ParaMEDMEM::CommInterface interface;
+  //
+  ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
+  ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
+  ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
+  //
+  MPI_Barrier(MPI_COMM_WORLD);
+  if(source_group->containsMyRank())
+    {
+      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
  *