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 "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(static_cast<MEDCouplingPointSet *>(const_cast<MEDCouplingMesh *>(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         /* Strange part of code localgroup not used !
 279         ProcessorGroup* localgroup;
 280         if (_source_group->containsMyRank())
 281           localgroup=_source_group;
 282         else
 283         localgroup=_target_group;*/
 284         delete _icoco_field;
 285
 286         _icoco_field=new ICoCo::MEDField(*const_cast<ICoCo::TrioField* >(triofield));
 287         attachLocalField(_icoco_field);
 288         return;
 289       }
 290     throw INTERP_KERNEL::Exception("incompatible field type");
 291   }
 292
 293   /*!
 294     Computes the field norm over its support
 295     on the source side and renormalizes the field on the target side
 296     so that the norms match.
 297
 298     \f[
 299     I_{source}=\sum_{i=1}^{n_{source}}V_{i}.|\Phi^{source}_{i}|^2,
 300     \f]
 301
 302     \f[
 303     I_{target}=\sum_{i=1}^{n_{target}}V_{i}.|\Phi^{target}_{i}|^2,
 304     \f]
 305
 306     \f[
 307     \Phi^{target}:=\Phi^{target}.\sqrt{I_{source}/I_{target}}.
 308     \f]
 309
 310   */
 311   void DisjointDEC::renormalizeTargetField(bool isWAbs)
 312   {
 313     if (_source_group->containsMyRank())
 314       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 315         {
 316           double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
 317           double source_norm = total_norm;
 318           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 319
 320         }
 321     if (_target_group->containsMyRank())
 322       {
 323         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 324           {
 325             double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
 326             double source_norm=total_norm;
 327             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 328
 329             if (fabs(total_norm)>1e-100)
 330               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
 331           }
 332       }
 333   }
 334   /*! @} */
 335
 336   bool DisjointDEC::isInSourceSide() const
 337   {
 338     if(!_source_group)
 339       return false;
 340     return _source_group->containsMyRank();
 341   }
 342
 343   bool DisjointDEC::isInTargetSide() const
 344   {
 345     if(!_target_group)
 346       return false;
 347     return _target_group->containsMyRank();
 348   }
 349
 350   bool DisjointDEC::isInUnion() const
 351   {
 352     if(!_union_group)
 353       return false;
 354     return _union_group->containsMyRank();
 355   }
 356
 357   void DisjointDEC::compareFieldAndMethod() const throw(INTERP_KERNEL::Exception)
 358   {
 359     if (_local_field)
 360       {
 361         TypeOfField entity = _local_field->getField()->getTypeOfField();
 362         if ( getMethod() == "P0" )
 363           {
 364             if ( entity != ON_CELLS )
 365               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 366                                              " For P0 interpolation, field must be on MED_CELL's");
 367           }
 368         else if ( getMethod() == "P1" )
 369           {
 370             if ( entity != ON_NODES )
 371               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 372                                              " For P1 interpolation, field must be on MED_NODE's");
 373           }
 374         else if ( getMethod() == "P1d" )
 375           {
 376             if ( entity != ON_CELLS )
 377               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 378                                              " For P1d interpolation, field must be on MED_CELL's");
 379             if ( _target_group->containsMyRank() )
 380               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
 381           }
 382         else
 383           {
 384             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
 385           }
 386       }
 387   }
 388
 389   /*!
 390     If way==true, source procs call sendData() and target procs call recvData().
 391     if way==false, it's the other way round.
 392   */
 393   void DisjointDEC::sendRecvData(bool way)
 394   {
 395     if(!isInUnion())
 396       return;
 397     if(isInSourceSide())
 398       {
 399         if(way)
 400           sendData();
 401         else
 402           recvData();
 403       }
 404     else if(isInTargetSide())
 405       {
 406         if(way)
 407           recvData();
 408         else
 409           sendData();
 410       }
 411   }
 412 }