src/ParaMEDMEM/DisjointDEC.cxx

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