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