src/ParaMEDMEM/DisjointDEC.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 "DisjointDEC.hxx"
  21 #include "CommInterface.hxx"
  22 #include "Topology.hxx"
  23 #include "BlockTopology.hxx"
  24 #include "ComponentTopology.hxx"
  25 #include "ParaFIELD.hxx"
  26 #include "ParaMESH.hxx"
  27 #include "ICoCoField.hxx"
  28 #include "ICoCoMEDField.hxx"
  29 #include "MPIProcessorGroup.hxx"
  30
  31 #include <cmath>
  32 #include <iostream>
  33
  34
  35 namespace ParaMEDMEM
  36 {
  37
  38   /*!
  39    * \anchor DisjointDEC-det
  40    * \class DisjointDEC
  41    *
  42    * \section DisjointDEC-over Overview
  43    *
  44    * Abstract interface class representing a link between two
  45    * processor groups for exchanging mesh or field data. The two processor groups must
  46    * have a void intersection (\ref ParaMEDMEM::OverlapDEC "OverlapDEC" is to be considered otherwise).
  47    * The %DEC is initialized by attaching a field on the receiving or on the
  48    * sending side.
  49    *
  50    * The data is sent or received through calls to the (abstract) methods recvData() and sendData().
  51    *
  52    * One can attach either a \c ParaMEDMEM::ParaFIELD, or a
  53    * \c ICoCo::Field, or directly a \c ParaMEDMEM::MEDCouplingFieldDouble instance.
  54    * See the various signatures of the method DisjointDEC::attachLocalField()
  55    *
  56    * The derivations of this class should be considered for practical instanciation:
  57    * - \ref InterpKernelDEC-det "InterpKernelDEC"
  58    * - \ref ExplicitCoincidentDEC-det "ExplicitCoincidentDEC"
  59    * - \ref StructuredCoincidentDEC-det "StructuredCoincidentDEC"
  60    *
  61    * \section DisjointDEC-options DisjointDEC options
  62    * The options supported by %DisjointDEC objects are the same that the ones supported for all
  63    * DECs in general and are all inherited from the class \ref ParaMEDMEM::DECOptions "DECOptions"
  64    *
  65   */
  66
  67   DisjointDEC::DisjointDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
  68       _local_field(0),
  69       _source_group(&source_group),
  70       _target_group(&target_group),
  71       _comm_interface(0),
  72       _owns_field(false),
  73       _owns_groups(false),
  74       _union_comm(MPI_COMM_NULL)
  75   {
  76     _union_group = source_group.fuse(target_group);
  77   }
  78
  79   DisjointDEC::DisjointDEC(const DisjointDEC& s):
  80       DEC(s),
  81       _local_field(0),
  82       _union_group(0),
  83       _source_group(0),
  84       _target_group(0),
  85       _comm_interface(0),
  86       _owns_field(false),
  87       _owns_groups(false),
  88       _union_comm(MPI_COMM_NULL)
  89   {
  90     copyInstance(s);
  91   }
  92
  93   DisjointDEC & DisjointDEC::operator=(const DisjointDEC& s)
  94   {
  95     cleanInstance();
  96     copyInstance(s);
  97     return *this;
  98
  99   }
 100
 101   void DisjointDEC::copyInstance(const DisjointDEC& other)
 102   {
 103     DEC::copyFrom(other);
 104     if(other._target_group)
 105       {
 106         _target_group=other._target_group->deepCpy();
 107         _owns_groups=true;
 108       }
 109     if(other._source_group)
 110       {
 111         _source_group=other._source_group->deepCpy();
 112         _owns_groups=true;
 113       }
 114     if (_source_group && _target_group)
 115       _union_group = _source_group->fuse(*_target_group);
 116   }
 117
 118   DisjointDEC::DisjointDEC(const std::set<int>& source_ids,
 119                            const std::set<int>& target_ids,
 120                            const MPI_Comm& world_comm):
 121      _local_field(0),
 122      _owns_field(false),
 123      _owns_groups(true),
 124      _comm_interface(0),
 125      _union_comm(MPI_COMM_NULL)
 126   {
 127     ParaMEDMEM::CommInterface comm;
 128     // Create the list of procs including source and target
 129     std::set<int> union_ids; // source and target ids in world_comm
 130     union_ids.insert(source_ids.begin(),source_ids.end());
 131     union_ids.insert(target_ids.begin(),target_ids.end());
 132     if(union_ids.size()!=(source_ids.size()+target_ids.size()))
 133       throw INTERP_KERNEL::Exception("DisjointDEC constructor : source_ids and target_ids overlap partially or fully. This type of DEC does not support it! OverlapDEC class could be the solution!");
 134     int* union_ranks_world=new int[union_ids.size()]; // ranks of sources and targets in world_comm
 135     std::copy(union_ids.begin(), union_ids.end(), union_ranks_world);
 136
 137     // Create a communicator on these procs
 138     MPI_Group union_group,world_group;
 139     comm.commGroup(world_comm,&world_group);
 140     comm.groupIncl(world_group,union_ids.size(),union_ranks_world,&union_group);
 141     comm.commCreate(world_comm,union_group,&_union_comm);
 142     delete[] union_ranks_world;
 143     if (_union_comm==MPI_COMM_NULL)
 144       { // This process is not in union
 145         _source_group=0;
 146         _target_group=0;
 147         _union_group=0;
 148         comm.groupFree(&union_group);
 149         comm.groupFree(&world_group);
 150         return;
 151       }
 152
 153     // Translate source_ids and target_ids from world_comm to union_comm
 154     int* source_ranks_world=new int[source_ids.size()]; // ranks of sources in world_comm
 155     std::copy(source_ids.begin(), source_ids.end(),source_ranks_world);
 156     int* source_ranks_union=new int[source_ids.size()]; // ranks of sources in union_comm
 157     int* target_ranks_world=new int[target_ids.size()]; // ranks of targets in world_comm
 158     std::copy(target_ids.begin(), target_ids.end(),target_ranks_world);
 159     int* target_ranks_union=new int[target_ids.size()]; // ranks of targets in union_comm
 160     MPI_Group_translate_ranks(world_group,source_ids.size(),source_ranks_world,union_group,source_ranks_union);
 161     MPI_Group_translate_ranks(world_group,target_ids.size(),target_ranks_world,union_group,target_ranks_union);
 162     std::set<int> source_ids_union;
 163     for (int i=0;i<(int)source_ids.size();i++)
 164       source_ids_union.insert(source_ranks_union[i]);
 165     std::set<int> target_ids_union;
 166     for (int i=0;i<(int)target_ids.size();i++)
 167       target_ids_union.insert(target_ranks_union[i]);
 168     delete [] source_ranks_world;
 169     delete [] source_ranks_union;
 170     delete [] target_ranks_world;
 171     delete [] target_ranks_union;
 172
 173     // Create the MPIProcessorGroups
 174     _source_group = new MPIProcessorGroup(comm,source_ids_union,_union_comm);
 175     _target_group = new MPIProcessorGroup(comm,target_ids_union,_union_comm);
 176     _union_group = _source_group->fuse(*_target_group);
 177     comm.groupFree(&union_group);
 178     comm.groupFree(&world_group);
 179   }
 180
 181   DisjointDEC::~DisjointDEC()
 182   {
 183     cleanInstance();
 184   }
 185
 186   void DisjointDEC::cleanInstance()
 187   {
 188     if(_owns_field)
 189       {
 190         delete _local_field;
 191       }
 192     _local_field=0;
 193     _owns_field=false;
 194     if(_owns_groups)
 195       {
 196         delete _source_group;
 197         delete _target_group;
 198       }
 199     _owns_groups=false;
 200     _source_group=0;
 201     _target_group=0;
 202     delete _union_group;
 203     _union_group=0;
 204     if (_union_comm != MPI_COMM_NULL)
 205       _comm_interface->commFree(&_union_comm);
 206   }
 207
 208   void DisjointDEC::setNature(NatureOfField nature)
 209   {
 210     if(_local_field)
 211       _local_field->getField()->setNature(nature);
 212   }
 213
 214   /*! Attaches a local field to a DEC.
 215     If the processor is on the receiving end of the DEC, the field
 216     will be updated by a recvData() call.
 217     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
 218   */
 219   void DisjointDEC::attachLocalField(const ParaFIELD *field, bool ownPt)
 220   {
 221     if(!isInUnion())
 222       return ;
 223     if(_owns_field)
 224       delete _local_field;
 225     _local_field=field;
 226     _owns_field=ownPt;
 227     _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
 228     compareFieldAndMethod();
 229   }
 230
 231   /*! Attaches a local field to a DEC. The method will test whether the processor
 232     is on the source or the target side and will associate the mesh underlying the
 233     field to the local side.
 234
 235     If the processor is on the receiving end of the DEC, the field
 236     will be updated by a recvData() call.
 237     Reversely, if the processor is on the sending end, the field will be read, possibly transformed,
 238     and sent appropriately to the other side.
 239   */
 240
 241   void DisjointDEC::attachLocalField(MEDCouplingFieldDouble *field)
 242   {
 243     if(!isInUnion())
 244       return ;
 245     ProcessorGroup* local_group;
 246     if (_source_group->containsMyRank())
 247       local_group=_source_group;
 248     else if (_target_group->containsMyRank())
 249       local_group=_target_group;
 250     else
 251       throw INTERP_KERNEL::Exception("Invalid procgroup for field attachment to DEC");
 252     ParaMESH *paramesh=new ParaMESH(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(field->getMesh())),*local_group,field->getMesh()->getName());
 253     ParaFIELD *tmp=new ParaFIELD(field, paramesh, *local_group);
 254     tmp->setOwnSupport(true);
 255     attachLocalField(tmp,true);
 256     //_comm_interface=&(local_group->getCommInterface());
 257   }
 258
 259   /*!
 260     Attaches a local field to a DEC.
 261     If the processor is on the receiving end of the DEC, the field
 262     will be updated by a recvData() call.
 263     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
 264     The field type is a generic ICoCo Field, so that the DEC can couple a number of different fields :
 265     - a ICoCo::MEDField, that is created from a MEDCoupling structure
 266
 267   */
 268   void DisjointDEC::attachLocalField(const ICoCo::MEDField *field)
 269   {
 270     if(!isInUnion())
 271       return ;
 272     if(!field)
 273       throw INTERP_KERNEL::Exception("DisjointDEC::attachLocalField : ICoCo::MEDField pointer is NULL !");
 274     attachLocalField(field->getField());
 275   }
 276
 277   /*!
 278     Computes the field norm over its support
 279     on the source side and renormalizes the field on the target side
 280     so that the norms match.
 281
 282     \f[
 283     I_{source}=\sum_{i=1}^{n_{source}}V_{i}.|\Phi^{source}_{i}|^2,
 284     \f]
 285
 286     \f[
 287     I_{target}=\sum_{i=1}^{n_{target}}V_{i}.|\Phi^{target}_{i}|^2,
 288     \f]
 289
 290     \f[
 291     \Phi^{target}:=\Phi^{target}.\sqrt{I_{source}/I_{target}}.
 292     \f]
 293
 294   */
 295   void DisjointDEC::renormalizeTargetField(bool isWAbs)
 296   {
 297     if (_source_group->containsMyRank())
 298       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 299         {
 300           double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
 301           double source_norm = total_norm;
 302           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 303
 304         }
 305     if (_target_group->containsMyRank())
 306       {
 307         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 308           {
 309             double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
 310             double source_norm=total_norm;
 311             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 312
 313             if (fabs(total_norm)>1e-100)
 314               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
 315           }
 316       }
 317   }
 318
 319   bool DisjointDEC::isInSourceSide() const
 320   {
 321     if(!_source_group)
 322       return false;
 323     return _source_group->containsMyRank();
 324   }
 325
 326   bool DisjointDEC::isInTargetSide() const
 327   {
 328     if(!_target_group)
 329       return false;
 330     return _target_group->containsMyRank();
 331   }
 332
 333   bool DisjointDEC::isInUnion() const
 334   {
 335     if(!_union_group)
 336       return false;
 337     return _union_group->containsMyRank();
 338   }
 339
 340   void DisjointDEC::compareFieldAndMethod() const throw(INTERP_KERNEL::Exception)
 341   {
 342     if (_local_field)
 343       {
 344         TypeOfField entity = _local_field->getField()->getTypeOfField();
 345         if ( getMethod() == "P0" )
 346           {
 347             if ( entity != ON_CELLS )
 348               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 349                                              " For P0 interpolation, field must be on MED_CELL's");
 350           }
 351         else if ( getMethod() == "P1" )
 352           {
 353             if ( entity != ON_NODES )
 354               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 355                                              " For P1 interpolation, field must be on MED_NODE's");
 356           }
 357         else if ( getMethod() == "P1d" )
 358           {
 359             if ( entity != ON_CELLS )
 360               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 361                                              " For P1d interpolation, field must be on MED_CELL's");
 362             if ( _target_group->containsMyRank() )
 363               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
 364           }
 365         else
 366           {
 367             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
 368           }
 369       }
 370   }
 371
 372   /*!
 373     If way==true, source procs call sendData() and target procs call recvData().
 374     if way==false, it's the other way round.
 375   */
 376   void DisjointDEC::sendRecvData(bool way)
 377   {
 378     if(!isInUnion())
 379       return;
 380     if(isInSourceSide())
 381       {
 382         if(way)
 383           sendData();
 384         else
 385           recvData();
 386       }
 387     else if(isInTargetSide())
 388       {
 389         if(way)
 390           recvData();
 391         else
 392           sendData();
 393       }
 394   }
 395 }