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
  50 //void ParaMEDMEMTest::testOverlapDEC_LMEC_seq()
  51 //{
  52 //  //  T_SC_Trio_src.med  -- "SupportOf_"
  53 //  //  T_SC_Trio_dst.med  -- "SupportOf_T_SC_Trio"
  54 //  //  h_TH_Trio_src.med  -- "SupportOf_"
  55 //  //  h_TH_Trio_dst.med  -- "SupportOf_h_TH_Trio"
  56 //  string rep("/export/home/adrien/support/antoine_LMEC/");
  57 //  string src_mesh_nam(rep + string("T_SC_Trio_src.med"));
  58 //  string tgt_mesh_nam(rep + string("T_SC_Trio_dst.med"));
  59 ////  string src_mesh_nam(rep + string("h_TH_Trio_src.med"));
  60 ////  string tgt_mesh_nam(rep + string("h_TH_Trio_dst.med"));
  61 //  MUMesh src_mesh=MEDLoader::ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
  62 //  MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
  63 ////  MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_h_TH_Trio",0);
  64 //
  65 //  MFDouble srcField = MEDCouplingFieldDouble::New(ON_CELLS, ONE_TIME);
  66 //  srcField->setMesh(src_mesh);
  67 //  DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1);
  68 //  dad->fillWithValue(1.0);
  69 //  srcField->setArray(dad);
  70 //  srcField->setNature(ConservativeVolumic);
  71 //
  72 //  MEDCouplingRemapper remap;
  73 //  remap.setOrientation(2); // always consider surface intersections as absolute areas.
  74 //  remap.prepare(src_mesh, tgt_mesh, "P0P0");
  75 //  MFDouble tgtField = remap.transferField(srcField, 1.0e+300);
  76 //  tgtField->setName("result");
  77 //  string out_nam(rep + string("adrien.med"));
  78 //  MEDLoader::WriteField(out_nam,tgtField, true);
  79 //  cout << "wrote: " << out_nam << "\n";
  80 //  double integ1 = 0.0, integ2 = 0.0;
  81 //  srcField->integral(true, &integ1);
  82 //  tgtField->integral(true, &integ2);
  83 ////  tgtField->reprQuickOverview(cout);
  84 //  CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8);
  85 //
  86 //  dad->decrRef();
  87 //}
  88 //
  89 //void ParaMEDMEMTest::testOverlapDEC_LMEC_para()
  90 //{
  91 //  using namespace ParaMEDMEM;
  92 //
  93 //  int size;
  94 //  int rank;
  95 //  MPI_Comm_size(MPI_COMM_WORLD,&size);
  96 //  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
  97 //
  98 //  if (size != 1) return ;
  99 //
 100 //  int nproc = 1;
 101 //  std::set<int> procs;
 102 //
 103 //  for (int i=0; i<nproc; i++)
 104 //    procs.insert(i);
 105 //
 106 //  CommInterface interface;
 107 //  OverlapDEC dec(procs);
 108 //
 109 //  ParaMESH* parameshS=0;
 110 //  ParaMESH* parameshT=0;
 111 //  ParaFIELD* parafieldS=0;
 112 //  ParaFIELD* parafieldT=0;
 113 //  MFDouble srcField;
 114 //
 115 //  // **** FILE LOADING
 116 //  //  T_SC_Trio_src.med  -- "SupportOf_"
 117 //  //  T_SC_Trio_dst.med  -- "SupportOf_T_SC_Trio"
 118 //  //  h_TH_Trio_src.med  -- "SupportOf_"
 119 //  //  h_TH_Trio_dst.med  -- "SupportOf_h_TH_Trio"
 120 //  string rep("/export/home/adrien/support/antoine_LMEC/");
 121 //  string src_mesh_nam(rep + string("T_SC_Trio_src.med"));
 122 //  string tgt_mesh_nam(rep + string("T_SC_Trio_dst.med"));
 123 //
 124 //
 125 //  MPI_Barrier(MPI_COMM_WORLD);
 126 //  if(rank==0)
 127 //    {
 128 //    //  string src_mesh_nam(rep + string("h_TH_Trio_src.med"));
 129 //    //  string tgt_mesh_nam(rep + string("h_TH_Trio_dst.med"));
 130 //      MUMesh src_mesh=MEDLoader::ReadUMeshFromFile(src_mesh_nam,"SupportOf_",0);
 131 //      MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_T_SC_Trio",0);
 132 //    //  MUMesh tgt_mesh=MEDLoader::ReadUMeshFromFile(tgt_mesh_nam,"SupportOf_h_TH_Trio",0);
 133 //
 134 //      // **** SOURCE
 135 //      srcField = MEDCouplingFieldDouble::New(ON_CELLS, ONE_TIME);
 136 //      srcField->setMesh(src_mesh);
 137 //      DataArrayDouble * dad = DataArrayDouble::New(); dad->alloc(src_mesh->getNumberOfCells(),1);
 138 //      dad->fillWithValue(1.0);
 139 //      srcField->setArray(dad);
 140 //      srcField->setNature(ConservativeVolumic);
 141 //
 142 //      ComponentTopology comptopo;
 143 //      parameshS = new ParaMESH(src_mesh,*dec.getGroup(),"source mesh");
 144 //      parafieldS = new ParaFIELD(ON_CELLS,ONE_TIME,parameshS,comptopo);
 145 //      parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 146 //      parafieldS->getField()->setArray(dad);
 147 //
 148 //      // **** TARGET
 149 //      parameshT=new ParaMESH(tgt_mesh,*dec.getGroup(),"target mesh");
 150 //      parafieldT=new ParaFIELD(ON_CELLS,ONE_TIME,parameshT,comptopo);
 151 //      parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 152 //      parafieldT->getField()->getArray()->fillWithValue(1.0e300);
 153 ////      valsT[0]=7.;
 154 //    }
 155 //  dec.setOrientation(2);
 156 //  dec.attachSourceLocalField(parafieldS);
 157 //  dec.attachTargetLocalField(parafieldT);
 158 //  dec.synchronize();
 159 //  dec.sendRecvData(true);
 160 //  //
 161 //  if(rank==0)
 162 //    {
 163 //      double integ1 = 0.0, integ2 = 0.0;
 164 //      MEDCouplingFieldDouble * tgtField;
 165 //
 166 //      srcField->integral(true, &integ1);
 167 //      tgtField = parafieldT->getField();
 168 ////      tgtField->reprQuickOverview(cout);
 169 //      tgtField->integral(true, &integ2);
 170 //      tgtField->setName("result");
 171 //      string out_nam(rep + string("adrien_para.med"));
 172 //      MEDLoader::WriteField(out_nam,tgtField, true);
 173 //      cout << "wrote: " << out_nam << "\n";
 174 //      CPPUNIT_ASSERT_DOUBLES_EQUAL(integ1,integ2,1e-8);
 175 //    }
 176 //  delete parafieldS;
 177 //  delete parafieldT;
 178 //  delete parameshS;
 179 //  delete parameshT;
 180 //
 181 //  MPI_Barrier(MPI_COMM_WORLD);
 182 //}
 183
 184 void prepareData1(int rank, ProcessorGroup * grp,
 185                                   MEDCouplingUMesh *& meshS, MEDCouplingUMesh *& meshT,
 186                                   ParaMESH*& parameshS, ParaMESH*& parameshT,
 187                                   ParaFIELD*& parafieldS, ParaFIELD*& parafieldT)
 188 {
 189   if(rank==0)
 190     {
 191       const double coordsS[10]={0.,0.,0.5,0.,1.,0.,0.,0.5,0.5,0.5};
 192       const double coordsT[6]={0.,0.,1.,0.,1.,1.};
 193       meshS=MEDCouplingUMesh::New();
 194       meshS->setMeshDimension(2);
 195       DataArrayDouble *myCoords=DataArrayDouble::New();
 196       myCoords->alloc(5,2);
 197       std::copy(coordsS,coordsS+10,myCoords->getPointer());
 198       meshS->setCoords(myCoords);
 199       myCoords->decrRef();
 200       int connS[7]={0,3,4,1, 1,4,2};
 201       meshS->allocateCells(2);
 202       meshS->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,connS);
 203       meshS->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connS+4);
 204       meshS->finishInsertingCells();
 205       ComponentTopology comptopo;
 206       parameshS=new ParaMESH(meshS, *grp,"source mesh");
 207       parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo);
 208       parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 209       double *valsS=parafieldS->getField()->getArray()->getPointer();
 210       valsS[0]=7.; valsS[1]=8.;
 211       //
 212       meshT=MEDCouplingUMesh::New();
 213       meshT->setMeshDimension(2);
 214       myCoords=DataArrayDouble::New();
 215       myCoords->alloc(3,2);
 216       std::copy(coordsT,coordsT+6,myCoords->getPointer());
 217       meshT->setCoords(myCoords);
 218       myCoords->decrRef();
 219       int connT[3]={0,2,1};
 220       meshT->allocateCells(1);
 221       meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
 222       meshT->finishInsertingCells();
 223       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 224       parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo);
 225       parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 226       double *valsT=parafieldT->getField()->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       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       ComponentTopology comptopo;
 247       parameshS=new ParaMESH(meshS,*grp,"source mesh");
 248       parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo);
 249       parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 250       double *valsS=parafieldS->getField()->getArray()->getPointer();
 251       valsS[0]=9.;
 252       valsS[1]=11.;
 253       //
 254       meshT=MEDCouplingUMesh::New();
 255       meshT->setMeshDimension(2);
 256       myCoords=DataArrayDouble::New();
 257       myCoords->alloc(3,2);
 258       std::copy(coordsT,coordsT+6,myCoords->getPointer());
 259       meshT->setCoords(myCoords);
 260       myCoords->decrRef();
 261       int connT[3]={0,2,1};
 262       meshT->allocateCells(1);
 263       meshT->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,connT);
 264       meshT->finishInsertingCells();
 265       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 266       parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo);
 267       parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 268       double *valsT=parafieldT->getField()->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       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       ComponentTopology comptopo;
 288       parameshS=new ParaMESH(meshS,*grp,"source mesh");
 289       parafieldS=new ParaFIELD(ON_CELLS,NO_TIME,parameshS,comptopo);
 290       parafieldS->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 291       double *valsS=parafieldS->getField()->getArray()->getPointer();
 292       valsS[0]=10.;
 293       //
 294       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       parameshT=new ParaMESH(meshT,*grp,"target mesh");
 306       parafieldT=new ParaFIELD(ON_CELLS,NO_TIME,parameshT,comptopo);
 307       parafieldT->getField()->setNature(ConservativeVolumic);//IntegralGlobConstraint
 308       double *valsT=parafieldT->getField()->getArray()->getPointer();
 309       valsT[0]=9.;
 310     }
 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 (size != 3) return ;
 328
 329   int nproc = 3;
 330   std::set<int> procs;
 331
 332   for (int i=0; i<nproc; i++)
 333     procs.insert(i);
 334
 335   CommInterface interface;
 336   OverlapDEC dec(procs);
 337   ProcessorGroup * grp = dec.getGroup();
 338   MEDCouplingUMesh* meshS=0, *meshT=0;
 339   ParaMESH* parameshS=0, *parameshT=0;
 340   ParaFIELD* parafieldS=0, *parafieldT=0;
 341
 342   MPI_Barrier(MPI_COMM_WORLD);
 343   prepareData1(rank, grp, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT);
 344
 345   /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 346    *  HACK ON BOUNDING BOX TO MAKE THIS CASE SIMPLE AND USABLE IN DEBUG
 347    * Bounding boxes are slightly smaller than should be thus localising the work to be done
 348    * and avoiding every proc talking to everyone else.
 349    * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 350    */
 351   dec.setBoundingBoxAdjustmentAbs(-1.0e-12);
 352
 353   dec.attachSourceLocalField(parafieldS);
 354   dec.attachTargetLocalField(parafieldT);
 355   dec.synchronize();
 356   dec.sendRecvData(true);
 357   //
 358   if(rank==0)
 359     {
 360       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.75,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
 361     }
 362   if(rank==1)
 363     {
 364       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
 365     }
 366   if(rank==2)
 367     {
 368       CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
 369     }
 370   delete parafieldS;
 371   delete parafieldT;
 372   delete parameshS;
 373   delete parameshT;
 374   meshS->decrRef();
 375   meshT->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 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
 392   int nproc = 3;
 393   std::set<int> procs;
 394
 395   for (int i=0; i<nproc; i++)
 396     procs.insert(i);
 397
 398   CommInterface interface;
 399   OverlapDEC dec(procs);
 400   ProcessorGroup * grp = dec.getGroup();
 401   MEDCouplingUMesh* meshS=0, *meshT=0;
 402   ParaMESH* parameshS=0, *parameshT=0;
 403   ParaFIELD* parafieldS=0, *parafieldT=0;
 404
 405   MPI_Barrier(MPI_COMM_WORLD);
 406   prepareData1(rank, grp, meshS, meshT, parameshS, parameshT, parafieldS, parafieldT);
 407
 408   /* Normal bounding boxes here: */
 409   dec.setBoundingBoxAdjustmentAbs(+1.0e-12);
 410
 411   dec.attachSourceLocalField(parafieldS);
 412   dec.attachTargetLocalField(parafieldT);
 413   dec.synchronize();
 414   dec.sendRecvData(true);
 415   //
 416   if(rank==0)
 417     {
 418       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.75,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
 419     }
 420   if(rank==1)
 421     {
 422       CPPUNIT_ASSERT_DOUBLES_EQUAL(8.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
 423     }
 424   if(rank==2)
 425     {
 426       CPPUNIT_ASSERT_DOUBLES_EQUAL(10.5,parafieldT->getField()->getArray()->getIJ(0,0),1e-12);
 427     }
 428   delete parafieldS;
 429   delete parafieldT;
 430   delete parameshS;
 431   delete parameshT;
 432   meshS->decrRef();
 433   meshT->decrRef();
 434
 435   MPI_Barrier(MPI_COMM_WORLD);
 436 }
 437