src/ParaMEDMEMTest/ParaMEDMEMTest_OverlapDEC.cxx

   1 // Copyright (C) 2007-2019  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 "MCAuto.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 MEDCoupling;
  46
  47 typedef  MCAuto<MEDCouplingUMesh> MUMesh;
  48 typedef  MCAuto<MEDCouplingFieldDouble> MFDouble;
  49 typedef  MCAuto<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=ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
  63 //  MUMesh tgt_mesh=ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
  64 ////  MUMesh tgt_mesh=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(IntensiveMaximum);
  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 //  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 MEDCoupling;
  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=ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
 132 //      MUMesh tgt_mesh=ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
 133 //    //  MUMesh tgt_mesh=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(IntensiveMaximum);
 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(IntensiveMaximum);//ExtensiveConservation
 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(IntensiveMaximum);//ExtensiveConservation
 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 //      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 void prepareData2_buildOneSquare(MEDCouplingUMesh* & meshS_0, MEDCouplingUMesh* & meshT_0)
 315 {
 316   const double coords[10] = {0.0,0.0,  0.0,1.0,  1.0,1.0,  1.0,0.0, 0.5,0.5};
 317   meshS_0 = MEDCouplingUMesh::New("source", 2);
 318   DataArrayDouble *myCoords=DataArrayDouble::New();
 319   myCoords->alloc(5,2);
 320   std::copy(coords,coords+10,myCoords->getPointer());
 321   meshS_0->setCoords(myCoords);  myCoords->decrRef();
 322   int connS[4]={0,1,2,3};
 323   meshS_0->allocateCells(2);
 324   meshS_0->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS);
 325   //
 326   meshT_0 = MEDCouplingUMesh::New("target", 2);
 327   myCoords=DataArrayDouble::New();
 328   myCoords->alloc(5,2);
 329   std::copy(coords,coords+10,myCoords->getPointer());
 330   meshT_0->setCoords(myCoords);
 331   myCoords->decrRef();
 332   int connT[12]={0,1,4,  1,2,4,  2,3,4,  3,0,4};
 333   meshT_0->allocateCells(4);
 334   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
 335   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+3);
 336   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+6);
 337   meshT_0->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT+9);
 338 }
 339
 340 /**
 341  * Prepare five (detached) QUAD4 disposed like this:
 342  *   (0)  (1)  (2)
 343  *   (3)  (4)
 344  *
 345  * On the target side the global mesh is identical except that each QUAD4 is split in 4 TRI3 (along the diagonals).
 346  * This is a case for two procs:
 347  *    - proc #0 has source squares 0,1,2 and target squares 0,3 (well, sets of TRI3s actually)
 348  *    - proc #1 has source squares 3,4 and target squares 1,2,4
 349  */
 350 void prepareData2(int rank, ProcessorGroup * grp, NatureOfField nature,
 351                   MEDCouplingUMesh *& meshS, MEDCouplingUMesh *& meshT,
 352                   ParaMESH*& parameshS, ParaMESH*& parameshT,
 353                   ParaFIELD*& parafieldS, ParaFIELD*& parafieldT,
 354                   bool stripPartOfSource=false,
 355                   int fieldCompoNum=1)
 356 {
 357   MEDCouplingUMesh *meshS_0 = 0, *meshT_0 = 0;
 358   prepareData2_buildOneSquare(meshS_0, meshT_0);
 359
 360   if(rank==0)
 361     {
 362       const double tr1[] = {1.5, 0.0};
 363       MEDCouplingUMesh *meshS_1 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCopy());
 364       meshS_1->translate(tr1);
 365       const double tr2[] = {3.0, 0.0};
 366       MEDCouplingUMesh *meshS_2 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCopy());
 367       meshS_2->translate(tr2);
 368
 369       std::vector<const MEDCouplingUMesh*> vec;
 370       vec.push_back(meshS_0);vec.push_back(meshS_1);
 371       if (!stripPartOfSource)
 372         vec.push_back(meshS_2);
 373       meshS = MEDCouplingUMesh::MergeUMeshes(vec);
 374       meshS_1->decrRef(); meshS_2->decrRef();
 375
 376       ComponentTopology comptopo(fieldCompoNum);
 377       parameshS=new ParaMESH(meshS, *grp,"source mesh");
 378       parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 379       parafieldS->getField()->setNature(nature);
 380       double *valsS=parafieldS->getField()->getArray()->getPointer();
 381       for(int i=0; i < fieldCompoNum; i++)
 382         {
 383           valsS[i] = 1. * (10^i);
 384           valsS[fieldCompoNum+i] = 2. * (10^i);
 385           if (!stripPartOfSource)
 386             {
 387               valsS[2*fieldCompoNum+i] = 3. * (10^i);
 388             }
 389         }
 390
 391       //
 392       const double tr3[] = {0.0, -1.5};
 393       MEDCouplingUMesh *meshT_3 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCopy());
 394       meshT_3->translate(tr3);
 395       vec.clear();
 396       vec.push_back(meshT_0);vec.push_back(meshT_3);
 397       meshT = MEDCouplingUMesh::MergeUMeshes(vec);
 398       meshT_3->decrRef();
 399
 400       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 401       parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 402       parafieldT->getField()->setNature(nature);
 403     }
 404   //
 405   if(rank==1)
 406     {
 407       const double tr3[] = {0.0, -1.5};
 408       MEDCouplingUMesh *meshS_3 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCopy());
 409       meshS_3->translate(tr3);
 410       const double tr4[] = {1.5, -1.5};
 411       MEDCouplingUMesh *meshS_4 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCopy());
 412       meshS_4->translate(tr4);
 413
 414       std::vector<const MEDCouplingUMesh*> vec;
 415       vec.push_back(meshS_3);vec.push_back(meshS_4);
 416       meshS = MEDCouplingUMesh::MergeUMeshes(vec);
 417       meshS_3->decrRef(); meshS_4->decrRef();
 418
 419       ComponentTopology comptopo(fieldCompoNum);
 420       parameshS=new ParaMESH(meshS, *grp,"source mesh");
 421       parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 422       parafieldS->getField()->setNature(nature);
 423       double *valsS=parafieldS->getField()->getArray()->getPointer();
 424       for(int i=0; i < fieldCompoNum; i++)
 425         {
 426           valsS[i] = 4. * (10^i);
 427           valsS[fieldCompoNum+i] = 5. * (10^i);
 428         }
 429
 430       //
 431       const double tr5[] = {1.5, 0.0};
 432       MEDCouplingUMesh *meshT_1 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCopy());
 433       meshT_1->translate(tr5);
 434       const double tr6[] = {3.0, 0.0};
 435       MEDCouplingUMesh *meshT_2 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCopy());
 436       meshT_2->translate(tr6);
 437       const double tr7[] = {1.5, -1.5};
 438       MEDCouplingUMesh *meshT_4 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCopy());
 439       meshT_4->translate(tr7);
 440
 441       vec.clear();
 442       vec.push_back(meshT_1);vec.push_back(meshT_2);vec.push_back(meshT_4);
 443       meshT = MEDCouplingUMesh::MergeUMeshes(vec);
 444       meshT_1->decrRef(); meshT_2->decrRef(); meshT_4->decrRef();
 445
 446       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 447       parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 448       parafieldT->getField()->setNature(nature);
 449     }
 450   meshS_0->decrRef();
 451   meshT_0->decrRef();
 452 }
 453
 454 /*! Test case from the official doc of the OverlapDEC.
 455  *  WARNING: bounding boxes might be tweaked here to make the case more interesting (i.e. to avoid an all to all exchange
 456  *  between all procs).
 457  */
 458 void testOverlapDEC_generic(int workSharingAlgo, double bbAdj)
 459 {
 460   int size, rank;
 461   MPI_Comm_size(MPI_COMM_WORLD,&size);
 462   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 463   //  char hostname[256];
 464   //  printf("(%d) PID %d on localhost ready for attach\n", rank, getpid());
 465   //  fflush(stdout);
 466
 467 //    if (rank == 0)
 468 //      {
 469 //        int i=1, j=0;
 470 //        while (i!=0)
 471 //          j=2;
 472 //      }
 473
 474   if (size != 3) return ;
 475   int nproc = 3;
 476   std::set<int> procs;
 477   for (int i=0; i<nproc; i++)
 478     procs.insert(i);
 479
 480   CommInterface interface;
 481   OverlapDEC dec(procs);
 482   MEDCouplingFieldDouble * mcfieldS=0, *mcfieldT=0;
 483
 484   prepareData1(rank, IntensiveMaximum, mcfieldS, mcfieldT);
 485
 486   // See comment in the caller:
 487   dec.setBoundingBoxAdjustmentAbs(bbAdj);
 488   dec.setWorkSharingAlgo(workSharingAlgo);  // just to ease debugging
 489
 490   dec.attachSourceLocalField(mcfieldS);
 491   dec.attachTargetLocalField(mcfieldT);
 492   dec.synchronize();
 493 //  dec.debugPrintWorkSharing(std::cout);
 494   dec.sendRecvData(true);
 495   //
 496   MPI_Barrier(MPI_COMM_WORLD);
 497   if(rank==0)
 498     {
 499       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.75,mcfieldT->getArray()->getIJ(0,0),1e-12);
 500     }
 501   if(rank==1)
 502     {
 503       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,mcfieldT->getArray()->getIJ(0,0),1e-12);
 504     }
 505   if(rank==2)
 506     {
 507       CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,mcfieldT->getArray()->getIJ(0,0),1e-12);
 508     }
 509
 510   mcfieldS->decrRef();
 511   mcfieldT->decrRef();
 512
 513   MPI_Barrier(MPI_COMM_WORLD);
 514 }
 515
 516 void ParaMEDMEMTest::testOverlapDEC1()
 517 {
 518   /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 519    *  HACK ON BOUNDING BOX TO MAKE THIS CASE SIMPLE AND USABLE IN DEBUG
 520    * Bounding boxes are slightly smaller than should be, thus localizing the work to be done
 521    * and avoiding every proc talking to everyone else.
 522    * Obviously this is NOT a good idea to do this in production code :-)
 523    * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 524    */
 525   testOverlapDEC_generic(0,-1.0e-12);
 526 }
 527
 528 void ParaMEDMEMTest::testOverlapDEC1_bis()
 529 {
 530    // Same BB hack as above
 531   testOverlapDEC_generic(1,-1.0e-12);
 532 }
 533
 534 void ParaMEDMEMTest::testOverlapDEC1_ter()
 535 {
 536    // Same BB hack as above
 537   testOverlapDEC_generic(2, -1.0e-12);
 538 }
 539
 540
 541 /*!
 542  * Same as testOverlapDEC1() but with regular bounding boxes. If you're looking for a nice debug case,
 543  * testOverlapDEC1() is identical in terms of geometry and field values, and more appropriate.
 544  */
 545 void ParaMEDMEMTest::testOverlapDEC2()
 546 {
 547   testOverlapDEC_generic(0,1.0e-12);
 548 }
 549
 550 void ParaMEDMEMTest::testOverlapDEC2_bis()
 551 {
 552   testOverlapDEC_generic(1,1.0e-12);
 553 }
 554
 555 void ParaMEDMEMTest::testOverlapDEC2_ter()
 556 {
 557   testOverlapDEC_generic(2,1.0e-12);
 558 }
 559
 560
 561 /*! Test focused on the mapping of cell IDs.
 562  * (i.e. when only part of the source/target mesh is transmitted)
 563  */
 564 void ParaMEDMEMTest::testOverlapDEC3()
 565 {
 566   int size, rank;
 567   MPI_Comm_size(MPI_COMM_WORLD,&size);
 568   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 569
 570   int nproc = 2;
 571   if (size != nproc) return ;
 572   std::set<int> procs;
 573   for (int i=0; i<nproc; i++)
 574     procs.insert(i);
 575
 576   CommInterface interface;
 577   OverlapDEC dec(procs);
 578   ProcessorGroup * grp = dec.getGroup();
 579   MEDCouplingUMesh* meshS=0, *meshT=0;
 580   ParaMESH* parameshS=0, *parameshT=0;
 581   ParaFIELD* parafieldS=0, *parafieldT=0;
 582
 583   prepareData2(rank, grp, IntensiveMaximum, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT);
 584
 585   dec.attachSourceLocalField(parafieldS);
 586   dec.attachTargetLocalField(parafieldT);
 587   dec.synchronize();
 588   dec.sendRecvData(true);
 589   //
 590   MEDCouplingFieldDouble * resField = parafieldT->getField();
 591   if(rank==0)
 592     {
 593       CPPUNIT_ASSERT_EQUAL(8, (int)resField->getNumberOfTuples());
 594       for(int i=0;i<4;i++)
 595         CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,resField->getArray()->getIJ(i,0),1e-12);
 596       for(int i=4;i<8;i++)
 597         CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0,resField->getArray()->getIJ(i,0),1e-12);
 598     }
 599   if(rank==1)
 600     {
 601       CPPUNIT_ASSERT_EQUAL(12, (int)resField->getNumberOfTuples());
 602       for(int i=0;i<4;i++)
 603         CPPUNIT_ASSERT_DOUBLES_EQUAL(2.0,resField->getArray()->getIJ(i,0),1e-12);
 604       for(int i=4;i<8;i++)
 605         CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0,resField->getArray()->getIJ(i,0),1e-12);
 606       for(int i=8;i<12;i++)
 607         CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0,resField->getArray()->getIJ(i,0),1e-12);
 608     }
 609   delete parafieldS;
 610   delete parafieldT;
 611   delete parameshS;
 612   delete parameshT;
 613   meshS->decrRef();
 614   meshT->decrRef();
 615
 616   MPI_Barrier(MPI_COMM_WORLD);
 617 }
 618
 619 /*!
 620  * Tests:
 621  *  - default value
 622  *  - multi-component fields
 623  */
 624 void ParaMEDMEMTest::testOverlapDEC4()
 625 {
 626   int size, rank;
 627   MPI_Comm_size(MPI_COMM_WORLD,&size);
 628   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 629
 630   int nproc = 2;
 631   if (size != nproc) return ;
 632   std::set<int> procs;
 633   for (int i=0; i<nproc; i++)
 634     procs.insert(i);
 635
 636   CommInterface interface;
 637   OverlapDEC dec(procs);
 638   ProcessorGroup * grp = dec.getGroup();
 639   MEDCouplingUMesh* meshS=0, *meshT=0;
 640   ParaMESH* parameshS=0, *parameshT=0;
 641   ParaFIELD* parafieldS=0, *parafieldT=0;
 642
 643   // As before, except than one of the source cell is removed, and that the field now has 2 components
 644   prepareData2(rank, grp, IntensiveMaximum, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT,
 645                true, 2);
 646 //  if (rank == 1)
 647 //    {
 648 //      int i=1, j=0;
 649 //      while (i!=0)
 650 //        j=2;
 651 //    }
 652
 653   dec.attachSourceLocalField(parafieldS);
 654   dec.attachTargetLocalField(parafieldT);
 655   double defVal = -300.0;
 656   dec.setDefaultValue(defVal);
 657   dec.synchronize();
 658   dec.sendRecvData(true);
 659   //
 660   MEDCouplingFieldDouble * resField = parafieldT->getField();
 661   if(rank==0)
 662     {
 663       CPPUNIT_ASSERT_EQUAL(8, (int)resField->getNumberOfTuples());
 664       for(int i=0;i<4;i++)
 665         {
 666           CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,resField->getArray()->getIJ(i*2,0),1e-12);
 667           CPPUNIT_ASSERT_DOUBLES_EQUAL(10.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 668         }
 669       for(int i=4;i<8;i++)
 670         {
 671           CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0,resField->getArray()->getIJ(i*2,0),1e-12);
 672           CPPUNIT_ASSERT_DOUBLES_EQUAL(40.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 673         }
 674     }
 675   if(rank==1)
 676     {
 677       CPPUNIT_ASSERT_EQUAL(12, (int)resField->getNumberOfTuples());
 678       for(int i=0;i<4;i++)
 679         {
 680           CPPUNIT_ASSERT_DOUBLES_EQUAL(2.0,resField->getArray()->getIJ(i*2,0),1e-12);
 681           CPPUNIT_ASSERT_DOUBLES_EQUAL(20.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 682         }
 683       // Default value should be here:
 684       for(int i=4;i<8;i++)
 685         {
 686           CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2,0),1e-12);
 687           CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2+1,0),1e-12);
 688         }
 689       for(int i=8;i<12;i++)
 690         {
 691           CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0,resField->getArray()->getIJ(i*2,0),1e-12);
 692           CPPUNIT_ASSERT_DOUBLES_EQUAL(50.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 693         }
 694     }
 695   delete parafieldS;
 696   delete parafieldT;
 697   delete parameshS;
 698   delete parameshT;
 699   meshS->decrRef();
 700   meshT->decrRef();
 701
 702   MPI_Barrier(MPI_COMM_WORLD);
 703 }
 704