src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "ParaMEDMEMTest.hxx"
  21 #include <cppunit/TestAssert.h>
  22
  23 #include "CommInterface.hxx"
  24 #include "ProcessorGroup.hxx"
  25 #include "MPIProcessorGroup.hxx"
  26 #include "Topology.hxx"
  27 #include "OverlapDEC.hxx"
  28 #include "ParaMESH.hxx"
  29 #include "ParaFIELD.hxx"
  30 #include "ComponentTopology.hxx"
  31
  32 #include "MEDCouplingUMesh.hxx"
  33
  34 #include <set>
  35
  36 using namespace std;
  37
  38 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
  39 #include "MEDLoader.hxx"
  40 #include "MEDLoaderBase.hxx"
  41 #include "MEDCouplingFieldDouble.hxx"
  42 #include "MEDCouplingMemArray.hxx"
  43 #include "MEDCouplingRemapper.hxx"
  44
  45 using namespace ParaMEDMEM;
  46
  47 typedef  MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> MUMesh;
  48 typedef  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> MFDouble;
  49 typedef  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> DADouble;
  50
  51 //void ParaMEDMEMTest::testOverlapDEC_LMEC_seq()
  52 //{
  53 //  //  T_SC_Trio_src.med  -- "SupportOf_"
  54 //  //  T_SC_Trio_dst.med  -- "SupportOf_T_SC_Trio"
  55 //  //  h_TH_Trio_src.med  -- "SupportOf_"
  56 //  //  h_TH_Trio_dst.med  -- "SupportOf_h_TH_Trio"
  57 //  string rep("/export/home/adrien/support/antoine_LMEC/");
  58 //  string src_mesh_nam(rep + string("T_SC_Trio_src.med"));
  59 //  string tgt_mesh_nam(rep + string("T_SC_Trio_dst.med"));
  60 ////  string src_mesh_nam(rep + string("h_TH_Trio_src.med"));
  61 ////  string tgt_mesh_nam(rep + string("h_TH_Trio_dst.med"));
  62 //  MUMesh src_mesh=MEDLoader::ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
  63 //  MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
  64 ////  MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_h_TH_Trio",0);
  65 //
  66 //  MFDouble srcField = MEDCouplingFieldDouble::New(ON_CELLS, ONE_TIME);
  67 //  srcField->setMesh(src_mesh);
  68 //  DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1);
  69 //  dad->fillWithValue(1.0);
  70 //  srcField->setArray(dad);
  71 //  srcField->setNature(ConservativeVolumic);
  72 //
  73 //  MEDCouplingRemapper remap;
  74 //  remap.setOrientation(2); // always consider surface intersections as absolute areas.
  75 //  remap.prepare(src_mesh, tgt_mesh, "P0P0");
  76 //  MFDouble tgtField = remap.transferField(srcField, 1.0e+300);
  77 //  tgtField->setName("result");
  78 //  string out_nam(rep + string("adrien.med"));
  79 //  MEDLoader::WriteField(out_nam,tgtField, true);
  80 //  cout << "wrote: " << out_nam << "\n";
  81 //  double integ1 = 0.0, integ2 = 0.0;
  82 //  srcField->integral(true, &integ1);
  83 //  tgtField->integral(true, &integ2);
  84 ////  tgtField->reprQuickOverview(cout);
  85 //  CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8);
  86 //
  87 //  dad->decrRef();
  88 //}
  89 //
  90 //void ParaMEDMEMTest::testOverlapDEC_LMEC_para()
  91 //{
  92 //  using namespace ParaMEDMEM;
  93 //
  94 //  int size;
  95 //  int rank;
  96 //  MPI_Comm_size(MPI_COMM_WORLD,&size);
  97 //  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
  98 //
  99 //  if (size != 1) return ;
 100 //
 101 //  int nproc = 1;
 102 //  std::set<int> procs;
 103 //
 104 //  for (int i=0; i<nproc; i++)
 105 //    procs.insert(i);
 106 //
 107 //  CommInterface interface;
 108 //  OverlapDEC dec(procs);
 109 //
 110 //  ParaMESH* parameshS=0;
 111 //  ParaMESH* parameshT=0;
 112 //  ParaFIELD* parafieldS=0;
 113 //  ParaFIELD* parafieldT=0;
 114 //  MFDouble srcField;
 115 //
 116 //  // **** FILE LOADING
 117 //  //  T_SC_Trio_src.med  -- "SupportOf_"
 118 //  //  T_SC_Trio_dst.med  -- "SupportOf_T_SC_Trio"
 119 //  //  h_TH_Trio_src.med  -- "SupportOf_"
 120 //  //  h_TH_Trio_dst.med  -- "SupportOf_h_TH_Trio"
 121 //  string rep("/export/home/adrien/support/antoine_LMEC/");
 122 //  string src_mesh_nam(rep + string("T_SC_Trio_src.med"));
 123 //  string tgt_mesh_nam(rep + string("T_SC_Trio_dst.med"));
 124 //
 125 //
 126 //  MPI_Barrier(MPI_COMM_WORLD);
 127 //  if(rank==0)
 128 //    {
 129 //    //  string src_mesh_nam(rep + string("h_TH_Trio_src.med"));
 130 //    //  string tgt_mesh_nam(rep + string("h_TH_Trio_dst.med"));
 131 //      MUMesh src_mesh=MEDLoader::ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
 132 //      MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
 133 //    //  MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_h_TH_Trio",0);
 134 //
 135 //      // **** SOURCE
 136 //      srcField = MEDCouplingFieldDouble::New(ON_CELLS, ONE_TIME);
 137 //      srcField->setMesh(src_mesh);
 138 //      DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1);
 139 //      dad->fillWithValue(1.0);
 140 //      srcField->setArray(dad);
 141 //      srcField->setNature(ConservativeVolumic);
 142 //
 143 //      ComponentTopology comptopo;
 144 //      parameshS = new ParaMESH(src_mesh,*dec.getGroup(),"source mesh");
 145 //      parafieldS = new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 146 //      parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 147 //      parafieldS->getField()->setArray(dad);
 148 //
 149 //      // **** TARGET
 150 //      parameshT=new ParaMESH(tgt_mesh,*dec.getGroup(),"target mesh");
 151 //      parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 152 //      parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 153 //      parafieldT->getField()->getArray()->fillWithValue(1.0e300);
 154 ////      valsT[0]=7.;
 155 //    }
 156 //  dec.setOrientation(2);
 157 //  dec.attachSourceLocalField(parafieldS);
 158 //  dec.attachTargetLocalField(parafieldT);
 159 //  dec.synchronize();
 160 //  dec.sendRecvData(true);
 161 //  //
 162 //  if(rank==0)
 163 //    {
 164 //      double integ1 = 0.0, integ2 = 0.0;
 165 //      MEDCouplingFieldDouble * tgtField;
 166 //
 167 //      srcField->integral(true, &integ1);
 168 //      tgtField = parafieldT->getField();
 169 ////      tgtField->reprQuickOverview(cout);
 170 //      tgtField->integral(true, &integ2);
 171 //      tgtField->setName("result");
 172 //      string out_nam(rep + string("adrien_para.med"));
 173 //      MEDLoader::WriteField(out_nam,tgtField, true);
 174 //      cout << "wrote: " << out_nam << "\n";
 175 //      CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8);
 176 //    }
 177 //  delete parafieldS;
 178 //  delete parafieldT;
 179 //  delete parameshS;
 180 //  delete parameshT;
 181 //
 182 //  MPI_Barrier(MPI_COMM_WORLD);
 183 //}
 184 //
 185 void prepareData1(int rank, NatureOfField nature,
 186                   MEDCouplingFieldDouble *& fieldS, MEDCouplingFieldDouble *& fieldT)
 187 {
 188   if(rank==0)
 189     {
 190       const double coordsS[10]={0.,0.,0.5,0.,1.,0.,0.,0.5,0.5,0.5};
 191       const double coordsT[6]={0.,0.,1.,0.,1.,1.};
 192       MUMesh meshS=MEDCouplingUMesh::New();
 193       meshS->setMeshDimension(2);
 194       DataArrayDouble *myCoords=DataArrayDouble::New();
 195       myCoords->alloc(5,2);
 196       std::copy(coordsS,coordsS+10,myCoords->getPointer());
 197       meshS->setCoords(myCoords);
 198       myCoords->decrRef();
 199       int connS[7]={0,3,4,1, 1,4,2};
 200       meshS->allocateCells(2);
 201       meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS);
 202       meshS->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connS+4);
 203       meshS->finishInsertingCells();
 204       fieldS = MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 205       DADouble arr = DataArrayDouble::New(); arr->alloc(meshS->getNumberOfCells(), 1);
 206       fieldS->setMesh(meshS); fieldS->setArray(arr);
 207       fieldS->setNature(nature);
 208       double *valsS=fieldS->getArray()->getPointer();
 209       valsS[0]=7.; valsS[1]=8.;
 210       //
 211       MUMesh meshT=MEDCouplingUMesh::New();
 212       meshT->setMeshDimension(2);
 213       myCoords=DataArrayDouble::New();
 214       myCoords->alloc(3,2);
 215       std::copy(coordsT,coordsT+6,myCoords->getPointer());
 216       meshT->setCoords(myCoords);
 217       myCoords->decrRef();
 218       int connT[3]={0,2,1};
 219       meshT->allocateCells(1);
 220       meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
 221       meshT->finishInsertingCells();
 222       fieldT = MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 223       DADouble arr2 = DataArrayDouble::New(); arr2->alloc(meshT->getNumberOfCells(), 1);
 224       fieldT->setMesh(meshT);  fieldT->setArray(arr2);
 225       fieldT->setNature(nature);
 226       double *valsT=fieldT->getArray()->getPointer();
 227       valsT[0]=7.;
 228     }
 229   //
 230   if(rank==1)
 231     {
 232       const double coordsS[10]={1.,0.,0.5,0.5,1.,0.5,0.5,1.,1.,1.};
 233       const double coordsT[6]={0.,0.,0.5,0.5,0.,1.};
 234       MUMesh meshS=MEDCouplingUMesh::New();
 235       meshS->setMeshDimension(2);
 236       DataArrayDouble *myCoords=DataArrayDouble::New();
 237       myCoords->alloc(5,2);
 238       std::copy(coordsS,coordsS+10,myCoords->getPointer());
 239       meshS->setCoords(myCoords);
 240       myCoords->decrRef();
 241       int connS[7]={0,1,2, 1,3,4,2};
 242       meshS->allocateCells(2);
 243       meshS->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connS);
 244       meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS+3);
 245       meshS->finishInsertingCells();
 246       fieldS = MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 247       DADouble arr = DataArrayDouble::New(); arr->alloc(meshS->getNumberOfCells(), 1);
 248       fieldS->setMesh(meshS); fieldS->setArray(arr);
 249       fieldS->setNature(nature);
 250       double *valsS=fieldS->getArray()->getPointer();
 251       valsS[0]=9.; valsS[1]=11.;
 252       //
 253       MUMesh meshT=MEDCouplingUMesh::New();
 254       meshT->setMeshDimension(2);
 255       myCoords=DataArrayDouble::New();
 256       myCoords->alloc(3,2);
 257       std::copy(coordsT,coordsT+6,myCoords->getPointer());
 258       meshT->setCoords(myCoords);
 259       myCoords->decrRef();
 260       int connT[3]={0,2,1};
 261       meshT->allocateCells(1);
 262       meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
 263       meshT->finishInsertingCells();
 264       fieldT = MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 265       DADouble arr2 = DataArrayDouble::New(); arr2->alloc(meshT->getNumberOfCells(), 1);
 266       fieldT->setMesh(meshT);  fieldT->setArray(arr2);
 267       fieldT->setNature(nature);
 268       double *valsT=fieldT->getArray()->getPointer();
 269       valsT[0]=8.;
 270     }
 271   //
 272   if(rank==2)
 273     {
 274       const double coordsS[8]={0.,0.5, 0.5,0.5, 0.,1., 0.5,1.};
 275       const double coordsT[6]={0.5,0.5,0.,1.,1.,1.};
 276       MUMesh meshS=MEDCouplingUMesh::New();
 277       meshS->setMeshDimension(2);
 278       DataArrayDouble *myCoords=DataArrayDouble::New();
 279       myCoords->alloc(4,2);
 280       std::copy(coordsS,coordsS+8,myCoords->getPointer());
 281       meshS->setCoords(myCoords);
 282       myCoords->decrRef();
 283       int connS[4]={0,2,3,1};
 284       meshS->allocateCells(1);
 285       meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS);
 286       meshS->finishInsertingCells();
 287       fieldS = MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 288       DADouble arr = DataArrayDouble::New(); arr->alloc(meshS->getNumberOfCells(), 1);
 289       fieldS->setMesh(meshS); fieldS->setArray(arr);
 290       fieldS->setNature(nature);
 291       double *valsS=fieldS->getArray()->getPointer();
 292       valsS[0]=10.;
 293       //
 294       MUMesh meshT=MEDCouplingUMesh::New();
 295       meshT->setMeshDimension(2);
 296       myCoords=DataArrayDouble::New();
 297       myCoords->alloc(3,2);
 298       std::copy(coordsT,coordsT+6,myCoords->getPointer());
 299       meshT->setCoords(myCoords);
 300       myCoords->decrRef();
 301       int connT[3]={0,1,2};
 302       meshT->allocateCells(1);
 303       meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
 304       meshT->finishInsertingCells();
 305       fieldT = MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 306       DADouble arr2 = DataArrayDouble::New(); arr2->alloc(meshT->getNumberOfCells(), 1);
 307       fieldT->setMesh(meshT); fieldT->setArray(arr2);
 308       fieldT->setNature(nature);
 309       double *valsT=fieldT->getArray()->getPointer();
 310       valsT[0]=9.;
 311     }
 312 }
 313
 314 /*! Test case from the official doc of the OverlapDEC.
 315  *  WARNING: bounding boxes are tweaked here to make the case more interesting (i.e. to avoid an all to all exchange
 316  *  between all procs).
 317  */
 318 void ParaMEDMEMTest::testOverlapDEC1()
 319 {
 320   int size, rank;
 321   MPI_Comm_size(MPI_COMM_WORLD,&size);
 322   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 323   //  char hostname[256];
 324   //  printf("(%d) PID %d on localhost ready for attach\n", rank, getpid());
 325   //  fflush(stdout);
 326
 327 //    if (rank == 1)
 328 //      {
 329 //        int i=1, j=0;
 330 //        while (i!=0)
 331 //          j=2;
 332 //      }
 333
 334   if (size != 3) return ;
 335   int nproc = 3;
 336   std::set<int> procs;
 337   for (int i=0; i<nproc; i++)
 338     procs.insert(i);
 339
 340   CommInterface interface;
 341   OverlapDEC dec(procs);
 342   MEDCouplingFieldDouble * mcfieldS=0, *mcfieldT=0;
 343
 344   prepareData1(rank, ConservativeVolumic, mcfieldS, mcfieldT);
 345
 346   /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 347    *  HACK ON BOUNDING BOX TO MAKE THIS CASE SIMPLE AND USABLE IN DEBUG
 348    * Bounding boxes are slightly smaller than should be, thus localizing the work to be done
 349    * and avoiding every proc talking to everyone else.
 350    * Obviously this is NOT a good idea to do this in production code :-)
 351    * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 352    */
 353   dec.setBoundingBoxAdjustmentAbs(-1.0e-12);
 354
 355   dec.attachSourceLocalField(mcfieldS);
 356   dec.attachTargetLocalField(mcfieldT);
 357   dec.synchronize();
 358   dec.sendRecvData(true);
 359   //
 360   MPI_Barrier(MPI_COMM_WORLD);
 361   if(rank==0)
 362     {
 363       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.75,mcfieldT->getArray()->getIJ(0,0),1e-12);
 364     }
 365   if(rank==1)
 366     {
 367       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,mcfieldT->getArray()->getIJ(0,0),1e-12);
 368     }
 369   if(rank==2)
 370     {
 371       CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,mcfieldT->getArray()->getIJ(0,0),1e-12);
 372     }
 373
 374   mcfieldS->decrRef();
 375   mcfieldT->decrRef();
 376
 377   MPI_Barrier(MPI_COMM_WORLD);
 378 }
 379
 380 /*!
 381  * Same as testOverlapDEC1() but with regular bounding boxes. If you're looking for a nice debug case,
 382  * testOverlapDEC1() is identical in terms of geometry and field values, and more appropriate.
 383  */
 384 void ParaMEDMEMTest::testOverlapDEC2()
 385 {
 386   int size, rank;
 387   MPI_Comm_size(MPI_COMM_WORLD,&size);
 388   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 389
 390   if (size != 3) return ;
 391   int nproc = 3;
 392   std::set<int> procs;
 393   for (int i=0; i<nproc; i++)
 394     procs.insert(i);
 395
 396 //      if (rank == 0)
 397 //        {
 398 //          int i=1, j=0;
 399 //          while (i!=0)
 400 //            j=2;
 401 //        }
 402
 403
 404   CommInterface interface;
 405   OverlapDEC dec(procs);
 406   MEDCouplingFieldDouble * mcfieldS=0, *mcfieldT=0;
 407   prepareData1(rank, ConservativeVolumic, mcfieldS, mcfieldT);
 408
 409   /* Normal bounding boxes here: */
 410   dec.setBoundingBoxAdjustmentAbs(+1.0e-12);
 411
 412   dec.attachSourceLocalField(mcfieldS);
 413   dec.attachTargetLocalField(mcfieldT);
 414   dec.synchronize();
 415   dec.sendRecvData(true);
 416   //
 417   if(rank==0)
 418     {
 419       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.75,mcfieldT->getArray()->getIJ(0,0),1e-12);
 420     }
 421   if(rank==1)
 422     {
 423       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,mcfieldT->getArray()->getIJ(0,0),1e-12);
 424     }
 425   if(rank==2)
 426     {
 427       CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,mcfieldT->getArray()->getIJ(0,0),1e-12);
 428     }
 429
 430   mcfieldS->decrRef();
 431   mcfieldT->decrRef();
 432
 433   MPI_Barrier(MPI_COMM_WORLD);
 434 }
 435
 436 void prepareData2_buildOneSquare(MEDCouplingUMesh* & meshS_0, MEDCouplingUMesh* & meshT_0)
 437 {
 438   const double coords[10] = {0.0,0.0,  0.0,1.0,  1.0,1.0,  1.0,0.0, 0.5,0.5};
 439   meshS_0 = MEDCouplingUMesh::New("source", 2);
 440   DataArrayDouble *myCoords=DataArrayDouble::New();
 441   myCoords->alloc(5,2);
 442   std::copy(coords,coords+10,myCoords->getPointer());
 443   meshS_0->setCoords(myCoords);  myCoords->decrRef();
 444   int connS[4]={0,1,2,3};
 445   meshS_0->allocateCells(2);
 446   meshS_0->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS);
 447   //
 448   meshT_0 = MEDCouplingUMesh::New("target", 2);
 449   myCoords=DataArrayDouble::New();
 450   myCoords->alloc(5,2);
 451   std::copy(coords,coords+10,myCoords->getPointer());
 452   meshT_0->setCoords(myCoords);
 453   myCoords->decrRef();
 454   int connT[12]={0,1,4,  1,2,4,  2,3,4,  3,0,4};
 455   meshT_0->allocateCells(4);
 456   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
 457   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+3);
 458   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+6);
 459   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+9);
 460
 461 }
 462
 463 /**
 464  * Prepare five (detached) QUAD4 disposed like this:
 465  *   (0)  (1)  (2)
 466  *   (3)  (4)
 467  *
 468  * On the target side the global mesh is identical except that each QUAD4 is split in 4 TRI3 (along the diagonals).
 469  * This is a case for two procs:
 470  *    - proc #0 has source squares 0,1,2 and target squares 0,3 (well, sets of TRI3s actually)
 471  *    - proc #1 has source squares 3,4 and target squares 1,2,4
 472  */
 473 void prepareData2(int rank, ProcessorGroup * grp, NatureOfField nature,
 474                   MEDCouplingUMesh *& meshS, MEDCouplingUMesh *& meshT,
 475                   ParaMESH*& parameshS, ParaMESH*& parameshT,
 476                   ParaFIELD*& parafieldS, ParaFIELD*& parafieldT,
 477                   bool stripPartOfSource=false,
 478                   int fieldCompoNum=1)
 479 {
 480   MEDCouplingUMesh *meshS_0 = 0, *meshT_0 = 0;
 481   prepareData2_buildOneSquare(meshS_0, meshT_0);
 482
 483   if(rank==0)
 484     {
 485       const double tr1[] = {1.5, 0.0};
 486       MEDCouplingUMesh *meshS_1 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 487       meshS_1->translate(tr1);
 488       const double tr2[] = {3.0, 0.0};
 489       MEDCouplingUMesh *meshS_2 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 490       meshS_2->translate(tr2);
 491
 492       std::vector<const MEDCouplingUMesh*> vec;
 493       vec.push_back(meshS_0);vec.push_back(meshS_1);
 494       if (!stripPartOfSource)
 495         vec.push_back(meshS_2);
 496       meshS = MEDCouplingUMesh::MergeUMeshes(vec);
 497       meshS_1->decrRef(); meshS_2->decrRef();
 498
 499       ComponentTopology comptopo(fieldCompoNum);
 500       parameshS=new ParaMESH(meshS, *grp,"source mesh");
 501       parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 502       parafieldS->getField()->setNature(nature);
 503       double *valsS=parafieldS->getField()->getArray()->getPointer();
 504       for(int i=0; i < fieldCompoNum; i++)
 505         {
 506           valsS[i] = 1. * (10^i);
 507           valsS[fieldCompoNum+i] = 2. * (10^i);
 508           if (!stripPartOfSource)
 509             {
 510               valsS[2*fieldCompoNum+i] = 3. * (10^i);
 511             }
 512         }
 513
 514       //
 515       const double tr3[] = {0.0, -1.5};
 516       MEDCouplingUMesh *meshT_3 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 517       meshT_3->translate(tr3);
 518       vec.clear();
 519       vec.push_back(meshT_0);vec.push_back(meshT_3);
 520       meshT = MEDCouplingUMesh::MergeUMeshes(vec);
 521       meshT_3->decrRef();
 522
 523       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 524       parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 525       parafieldT->getField()->setNature(nature);
 526     }
 527   //
 528   if(rank==1)
 529     {
 530       const double tr3[] = {0.0, -1.5};
 531       MEDCouplingUMesh *meshS_3 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 532       meshS_3->translate(tr3);
 533       const double tr4[] = {1.5, -1.5};
 534       MEDCouplingUMesh *meshS_4 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 535       meshS_4->translate(tr4);
 536
 537       std::vector<const MEDCouplingUMesh*> vec;
 538       vec.push_back(meshS_3);vec.push_back(meshS_4);
 539       meshS = MEDCouplingUMesh::MergeUMeshes(vec);
 540       meshS_3->decrRef(); meshS_4->decrRef();
 541
 542       ComponentTopology comptopo(fieldCompoNum);
 543       parameshS=new ParaMESH(meshS, *grp,"source mesh");
 544       parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 545       parafieldS->getField()->setNature(nature);
 546       double *valsS=parafieldS->getField()->getArray()->getPointer();
 547       for(int i=0; i < fieldCompoNum; i++)
 548         {
 549           valsS[i] = 4. * (10^i);
 550           valsS[fieldCompoNum+i] = 5. * (10^i);
 551         }
 552
 553       //
 554       const double tr5[] = {1.5, 0.0};
 555       MEDCouplingUMesh *meshT_1 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 556       meshT_1->translate(tr5);
 557       const double tr6[] = {3.0, 0.0};
 558       MEDCouplingUMesh *meshT_2 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 559       meshT_2->translate(tr6);
 560       const double tr7[] = {1.5, -1.5};
 561       MEDCouplingUMesh *meshT_4 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 562       meshT_4->translate(tr7);
 563
 564       vec.clear();
 565       vec.push_back(meshT_1);vec.push_back(meshT_2);vec.push_back(meshT_4);
 566       meshT = MEDCouplingUMesh::MergeUMeshes(vec);
 567       meshT_1->decrRef(); meshT_2->decrRef(); meshT_4->decrRef();
 568
 569       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 570       parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 571       parafieldT->getField()->setNature(nature);
 572     }
 573   meshS_0->decrRef();
 574   meshT_0->decrRef();
 575 }
 576
 577 /*! Test focused on the mapping of cell IDs.
 578  * (i.e. when only part of the source/target mesh is transmitted)
 579  */
 580 void ParaMEDMEMTest::testOverlapDEC3()
 581 {
 582   int size, rank;
 583   MPI_Comm_size(MPI_COMM_WORLD,&size);
 584   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 585
 586   int nproc = 2;
 587   if (size != nproc) return ;
 588   std::set<int> procs;
 589   for (int i=0; i<nproc; i++)
 590     procs.insert(i);
 591
 592   CommInterface interface;
 593   OverlapDEC dec(procs);
 594   ProcessorGroup * grp = dec.getGroup();
 595   MEDCouplingUMesh* meshS=0, *meshT=0;
 596   ParaMESH* parameshS=0, *parameshT=0;
 597   ParaFIELD* parafieldS=0, *parafieldT=0;
 598
 599   prepareData2(rank, grp, ConservativeVolumic, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT);
 600
 601   dec.attachSourceLocalField(parafieldS);
 602   dec.attachTargetLocalField(parafieldT);
 603   dec.synchronize();
 604   dec.sendRecvData(true);
 605   //
 606   MEDCouplingFieldDouble * resField = parafieldT->getField();
 607   if(rank==0)
 608     {
 609       CPPUNIT_ASSERT_EQUAL(8, resField->getNumberOfTuples());
 610       for(int i=0;i<4;i++)
 611         CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,resField->getArray()->getIJ(i,0),1e-12);
 612       for(int i=4;i<8;i++)
 613         CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0,resField->getArray()->getIJ(i,0),1e-12);
 614     }
 615   if(rank==1)
 616     {
 617       CPPUNIT_ASSERT_EQUAL(12, resField->getNumberOfTuples());
 618       for(int i=0;i<4;i++)
 619         CPPUNIT_ASSERT_DOUBLES_EQUAL(2.0,resField->getArray()->getIJ(i,0),1e-12);
 620       for(int i=4;i<8;i++)
 621         CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0,resField->getArray()->getIJ(i,0),1e-12);
 622       for(int i=8;i<12;i++)
 623         CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0,resField->getArray()->getIJ(i,0),1e-12);
 624     }
 625   delete parafieldS;
 626   delete parafieldT;
 627   delete parameshS;
 628   delete parameshT;
 629   meshS->decrRef();
 630   meshT->decrRef();
 631
 632   MPI_Barrier(MPI_COMM_WORLD);
 633 }
 634
 635 /*!
 636  * Tests:
 637  *  - default value
 638  *  - multi-component fields
 639  */
 640 void ParaMEDMEMTest::testOverlapDEC4()
 641 {
 642   int size, rank;
 643   MPI_Comm_size(MPI_COMM_WORLD,&size);
 644   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 645
 646   int nproc = 2;
 647   if (size != nproc) return ;
 648   std::set<int> procs;
 649   for (int i=0; i<nproc; i++)
 650     procs.insert(i);
 651
 652   CommInterface interface;
 653   OverlapDEC dec(procs);
 654   ProcessorGroup * grp = dec.getGroup();
 655   MEDCouplingUMesh* meshS=0, *meshT=0;
 656   ParaMESH* parameshS=0, *parameshT=0;
 657   ParaFIELD* parafieldS=0, *parafieldT=0;
 658
 659   // As before, except than one of the source cell is removed, and that the field now has 2 components
 660   prepareData2(rank, grp, ConservativeVolumic, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT,
 661                true, 2);
 662 //  if (rank == 1)
 663 //    {
 664 //      int i=1, j=0;
 665 //      while (i!=0)
 666 //        j=2;
 667 //    }
 668
 669   dec.attachSourceLocalField(parafieldS);
 670   dec.attachTargetLocalField(parafieldT);
 671   double defVal = -300.0;
 672   dec.setDefaultValue(defVal);
 673   dec.synchronize();
 674   dec.sendRecvData(true);
 675   //
 676   MEDCouplingFieldDouble * resField = parafieldT->getField();
 677   if(rank==0)
 678     {
 679       CPPUNIT_ASSERT_EQUAL(8, resField->getNumberOfTuples());
 680       for(int i=0;i<4;i++)
 681         {
 682           CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,resField->getArray()->getIJ(i*2,0),1e-12);
 683           CPPUNIT_ASSERT_DOUBLES_EQUAL(10.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 684         }
 685       for(int i=4;i<8;i++)
 686         {
 687           CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0,resField->getArray()->getIJ(i*2,0),1e-12);
 688           CPPUNIT_ASSERT_DOUBLES_EQUAL(40.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 689         }
 690     }
 691   if(rank==1)
 692     {
 693       CPPUNIT_ASSERT_EQUAL(12, resField->getNumberOfTuples());
 694       for(int i=0;i<4;i++)
 695         {
 696           CPPUNIT_ASSERT_DOUBLES_EQUAL(2.0,resField->getArray()->getIJ(i*2,0),1e-12);
 697           CPPUNIT_ASSERT_DOUBLES_EQUAL(20.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 698         }
 699       // Default value should be here:
 700       for(int i=4;i<8;i++)
 701         {
 702           CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2,0),1e-12);
 703           CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2+1,0),1e-12);
 704         }
 705       for(int i=8;i<12;i++)
 706         {
 707           CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0,resField->getArray()->getIJ(i*2,0),1e-12);
 708           CPPUNIT_ASSERT_DOUBLES_EQUAL(50.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 709         }
 710     }
 711   delete parafieldS;
 712   delete parafieldT;
 713   delete parameshS;
 714   delete parameshT;
 715   meshS->decrRef();
 716   meshT->decrRef();
 717
 718   MPI_Barrier(MPI_COMM_WORLD);
 719
 720 }
 721