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 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 /**
 342  * Prepare five (detached) QUAD4 disposed like this:
 343  *   (0)  (1)  (2)
 344  *   (3)  (4)
 345  *
 346  * On the target side the global mesh is identical except that each QUAD4 is split in 4 TRI3 (along the diagonals).
 347  * This is a case for two procs:
 348  *    - proc #0 has source squares 0,1,2 and target squares 0,3 (well, sets of TRI3s actually)
 349  *    - proc #1 has source squares 3,4 and target squares 1,2,4
 350  */
 351 void prepareData2(int rank, ProcessorGroup * grp, NatureOfField nature,
 352                   MEDCouplingUMesh *& meshS, MEDCouplingUMesh *& meshT,
 353                   ParaMESH*& parameshS, ParaMESH*& parameshT,
 354                   ParaFIELD*& parafieldS, ParaFIELD*& parafieldT,
 355                   bool stripPartOfSource=false,
 356                   int fieldCompoNum=1)
 357 {
 358   MEDCouplingUMesh *meshS_0 = 0, *meshT_0 = 0;
 359   prepareData2_buildOneSquare(meshS_0, meshT_0);
 360
 361   if(rank==0)
 362     {
 363       const double tr1[] = {1.5, 0.0};
 364       MEDCouplingUMesh *meshS_1 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 365       meshS_1->translate(tr1);
 366       const double tr2[] = {3.0, 0.0};
 367       MEDCouplingUMesh *meshS_2 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 368       meshS_2->translate(tr2);
 369
 370       std::vector<const MEDCouplingUMesh*> vec;
 371       vec.push_back(meshS_0);vec.push_back(meshS_1);
 372       if (!stripPartOfSource)
 373         vec.push_back(meshS_2);
 374       meshS = MEDCouplingUMesh::MergeUMeshes(vec);
 375       meshS_1->decrRef(); meshS_2->decrRef();
 376
 377       ComponentTopology comptopo(fieldCompoNum);
 378       parameshS=new ParaMESH(meshS, *grp,"source mesh");
 379       parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 380       parafieldS->getField()->setNature(nature);
 381       double *valsS=parafieldS->getField()->getArray()->getPointer();
 382       for(int i=0; i < fieldCompoNum; i++)
 383         {
 384           valsS[i] = 1. * (10^i);
 385           valsS[fieldCompoNum+i] = 2. * (10^i);
 386           if (!stripPartOfSource)
 387             {
 388               valsS[2*fieldCompoNum+i] = 3. * (10^i);
 389             }
 390         }
 391
 392       //
 393       const double tr3[] = {0.0, -1.5};
 394       MEDCouplingUMesh *meshT_3 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 395       meshT_3->translate(tr3);
 396       vec.clear();
 397       vec.push_back(meshT_0);vec.push_back(meshT_3);
 398       meshT = MEDCouplingUMesh::MergeUMeshes(vec);
 399       meshT_3->decrRef();
 400
 401       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 402       parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 403       parafieldT->getField()->setNature(nature);
 404     }
 405   //
 406   if(rank==1)
 407     {
 408       const double tr3[] = {0.0, -1.5};
 409       MEDCouplingUMesh *meshS_3 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 410       meshS_3->translate(tr3);
 411       const double tr4[] = {1.5, -1.5};
 412       MEDCouplingUMesh *meshS_4 = static_cast<MEDCouplingUMesh*>(meshS_0->deepCpy());
 413       meshS_4->translate(tr4);
 414
 415       std::vector<const MEDCouplingUMesh*> vec;
 416       vec.push_back(meshS_3);vec.push_back(meshS_4);
 417       meshS = MEDCouplingUMesh::MergeUMeshes(vec);
 418       meshS_3->decrRef(); meshS_4->decrRef();
 419
 420       ComponentTopology comptopo(fieldCompoNum);
 421       parameshS=new ParaMESH(meshS, *grp,"source mesh");
 422       parafieldS=new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 423       parafieldS->getField()->setNature(nature);
 424       double *valsS=parafieldS->getField()->getArray()->getPointer();
 425       for(int i=0; i < fieldCompoNum; i++)
 426         {
 427           valsS[i] = 4. * (10^i);
 428           valsS[fieldCompoNum+i] = 5. * (10^i);
 429         }
 430
 431       //
 432       const double tr5[] = {1.5, 0.0};
 433       MEDCouplingUMesh *meshT_1 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 434       meshT_1->translate(tr5);
 435       const double tr6[] = {3.0, 0.0};
 436       MEDCouplingUMesh *meshT_2 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 437       meshT_2->translate(tr6);
 438       const double tr7[] = {1.5, -1.5};
 439       MEDCouplingUMesh *meshT_4 = static_cast<MEDCouplingUMesh*>(meshT_0->deepCpy());
 440       meshT_4->translate(tr7);
 441
 442       vec.clear();
 443       vec.push_back(meshT_1);vec.push_back(meshT_2);vec.push_back(meshT_4);
 444       meshT = MEDCouplingUMesh::MergeUMeshes(vec);
 445       meshT_1->decrRef(); meshT_2->decrRef(); meshT_4->decrRef();
 446
 447       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 448       parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 449       parafieldT->getField()->setNature(nature);
 450     }
 451   meshS_0->decrRef();
 452   meshT_0->decrRef();
 453 }
 454
 455 /*! Test case from the official doc of the OverlapDEC.
 456  *  WARNING: bounding boxes might be tweaked here to make the case more interesting (i.e. to avoid an all to all exchange
 457  *  between all procs).
 458  */
 459 void testOverlapDEC_generic(int workSharingAlgo, double bbAdj)
 460 {
 461   int size, rank;
 462   MPI_Comm_size(MPI_COMM_WORLD,&size);
 463   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 464   //  char hostname[256];
 465   //  printf("(%d) PID %d on localhost ready for attach\n", rank, getpid());
 466   //  fflush(stdout);
 467
 468 //    if (rank == 0)
 469 //      {
 470 //        int i=1, j=0;
 471 //        while (i!=0)
 472 //          j=2;
 473 //      }
 474
 475   if (size != 3) return ;
 476   int nproc = 3;
 477   std::set<int> procs;
 478   for (int i=0; i<nproc; i++)
 479     procs.insert(i);
 480
 481   CommInterface interface;
 482   OverlapDEC dec(procs);
 483   MEDCouplingFieldDouble * mcfieldS=0, *mcfieldT=0;
 484
 485   prepareData1(rank, ConservativeVolumic, mcfieldS, mcfieldT);
 486
 487   // See comment in the caller:
 488   dec.setBoundingBoxAdjustmentAbs(bbAdj);
 489   dec.setWorkSharingAlgo(workSharingAlgo);  // just to ease debugging
 490
 491   dec.attachSourceLocalField(mcfieldS);
 492   dec.attachTargetLocalField(mcfieldT);
 493   dec.synchronize();
 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   /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 531    * Same BB hack as above
 532    * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 533    */
 534   testOverlapDEC_generic(1,-1.0e-12);
 535 }
 536
 537 /*!
 538  * Same as testOverlapDEC1() but with regular bounding boxes. If you're looking for a nice debug case,
 539  * testOverlapDEC1() is identical in terms of geometry and field values, and more appropriate.
 540  */
 541 void ParaMEDMEMTest::testOverlapDEC2()
 542 {
 543   testOverlapDEC_generic(0,1.0e-12);
 544 }
 545
 546 void ParaMEDMEMTest::testOverlapDEC2_bis()
 547 {
 548   testOverlapDEC_generic(1,1.0e-12);
 549 }
 550
 551
 552 /*! Test focused on the mapping of cell IDs.
 553  * (i.e. when only part of the source/target mesh is transmitted)
 554  */
 555 void ParaMEDMEMTest::testOverlapDEC3()
 556 {
 557   int size, rank;
 558   MPI_Comm_size(MPI_COMM_WORLD,&size);
 559   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 560
 561   int nproc = 2;
 562   if (size != nproc) return ;
 563   std::set<int> procs;
 564   for (int i=0; i<nproc; i++)
 565     procs.insert(i);
 566
 567   CommInterface interface;
 568   OverlapDEC dec(procs);
 569   ProcessorGroup * grp = dec.getGroup();
 570   MEDCouplingUMesh* meshS=0, *meshT=0;
 571   ParaMESH* parameshS=0, *parameshT=0;
 572   ParaFIELD* parafieldS=0, *parafieldT=0;
 573
 574   prepareData2(rank, grp, ConservativeVolumic, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT);
 575
 576   dec.attachSourceLocalField(parafieldS);
 577   dec.attachTargetLocalField(parafieldT);
 578   dec.synchronize();
 579   dec.sendRecvData(true);
 580   //
 581   MEDCouplingFieldDouble * resField = parafieldT->getField();
 582   if(rank==0)
 583     {
 584       CPPUNIT_ASSERT_EQUAL(8, resField->getNumberOfTuples());
 585       for(int i=0;i<4;i++)
 586         CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,resField->getArray()->getIJ(i,0),1e-12);
 587       for(int i=4;i<8;i++)
 588         CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0,resField->getArray()->getIJ(i,0),1e-12);
 589     }
 590   if(rank==1)
 591     {
 592       CPPUNIT_ASSERT_EQUAL(12, resField->getNumberOfTuples());
 593       for(int i=0;i<4;i++)
 594         CPPUNIT_ASSERT_DOUBLES_EQUAL(2.0,resField->getArray()->getIJ(i,0),1e-12);
 595       for(int i=4;i<8;i++)
 596         CPPUNIT_ASSERT_DOUBLES_EQUAL(3.0,resField->getArray()->getIJ(i,0),1e-12);
 597       for(int i=8;i<12;i++)
 598         CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0,resField->getArray()->getIJ(i,0),1e-12);
 599     }
 600   delete parafieldS;
 601   delete parafieldT;
 602   delete parameshS;
 603   delete parameshT;
 604   meshS->decrRef();
 605   meshT->decrRef();
 606
 607   MPI_Barrier(MPI_COMM_WORLD);
 608 }
 609
 610 /*!
 611  * Tests:
 612  *  - default value
 613  *  - multi-component fields
 614  */
 615 void ParaMEDMEMTest::testOverlapDEC4()
 616 {
 617   int size, rank;
 618   MPI_Comm_size(MPI_COMM_WORLD,&size);
 619   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 620
 621   int nproc = 2;
 622   if (size != nproc) return ;
 623   std::set<int> procs;
 624   for (int i=0; i<nproc; i++)
 625     procs.insert(i);
 626
 627   CommInterface interface;
 628   OverlapDEC dec(procs);
 629   ProcessorGroup * grp = dec.getGroup();
 630   MEDCouplingUMesh* meshS=0, *meshT=0;
 631   ParaMESH* parameshS=0, *parameshT=0;
 632   ParaFIELD* parafieldS=0, *parafieldT=0;
 633
 634   // As before, except than one of the source cell is removed, and that the field now has 2 components
 635   prepareData2(rank, grp, ConservativeVolumic, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT,
 636                true, 2);
 637 //  if (rank == 1)
 638 //    {
 639 //      int i=1, j=0;
 640 //      while (i!=0)
 641 //        j=2;
 642 //    }
 643
 644   dec.attachSourceLocalField(parafieldS);
 645   dec.attachTargetLocalField(parafieldT);
 646   double defVal = -300.0;
 647   dec.setDefaultValue(defVal);
 648   dec.synchronize();
 649   dec.sendRecvData(true);
 650   //
 651   MEDCouplingFieldDouble * resField = parafieldT->getField();
 652   if(rank==0)
 653     {
 654       CPPUNIT_ASSERT_EQUAL(8, resField->getNumberOfTuples());
 655       for(int i=0;i<4;i++)
 656         {
 657           CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,resField->getArray()->getIJ(i*2,0),1e-12);
 658           CPPUNIT_ASSERT_DOUBLES_EQUAL(10.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 659         }
 660       for(int i=4;i<8;i++)
 661         {
 662           CPPUNIT_ASSERT_DOUBLES_EQUAL(4.0,resField->getArray()->getIJ(i*2,0),1e-12);
 663           CPPUNIT_ASSERT_DOUBLES_EQUAL(40.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 664         }
 665     }
 666   if(rank==1)
 667     {
 668       CPPUNIT_ASSERT_EQUAL(12, resField->getNumberOfTuples());
 669       for(int i=0;i<4;i++)
 670         {
 671           CPPUNIT_ASSERT_DOUBLES_EQUAL(2.0,resField->getArray()->getIJ(i*2,0),1e-12);
 672           CPPUNIT_ASSERT_DOUBLES_EQUAL(20.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 673         }
 674       // Default value should be here:
 675       for(int i=4;i<8;i++)
 676         {
 677           CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2,0),1e-12);
 678           CPPUNIT_ASSERT_DOUBLES_EQUAL(defVal,resField->getArray()->getIJ(i*2+1,0),1e-12);
 679         }
 680       for(int i=8;i<12;i++)
 681         {
 682           CPPUNIT_ASSERT_DOUBLES_EQUAL(5.0,resField->getArray()->getIJ(i*2,0),1e-12);
 683           CPPUNIT_ASSERT_DOUBLES_EQUAL(50.0,resField->getArray()->getIJ(i*2+1,0),1e-12);
 684         }
 685     }
 686   delete parafieldS;
 687   delete parafieldT;
 688   delete parameshS;
 689   delete parameshT;
 690   meshS->decrRef();
 691   meshT->decrRef();
 692
 693   MPI_Barrier(MPI_COMM_WORLD);
 694 }
 695