src/ParaMEDMEM/DisjointDEC.cxx

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