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