src/ParaMEDMEM/DEC.cxx

   1 //  Copyright (C) 2007-2008  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 #include "CommInterface.hxx"
  20 #include "Topology.hxx"
  21 #include "BlockTopology.hxx"
  22 #include "ComponentTopology.hxx"
  23 #include "ParaFIELD.hxx"
  24 #include "DEC.hxx"
  25 #include "ICoCoField.hxx"
  26 #include "ICoCoMEDField.hxx"
  27 #include "ICoCoTrioField.hxx"
  28 #include "MPIProcessorGroup.hxx"
  29
  30 #include <cmath>
  31
  32 /*! \defgroup dec DEC
  33  *
  34  * \section decintroduction Introduction
  35  *
  36  * Interface class for creation of a link between two
  37  * processor groups for exhanging mesh or field data.
  38  * The DEC is defined by attaching a field on the receiving or on the
  39  * sending side.
  40  * On top of attaching a ParaMEDMEM::FIELD, it is possible to
  41  * attach a ICoCo::Field. This class is an abstract class that enables
  42  * coupling of codes that respect the ICoCo interface \ref icoco. It has two implementations:
  43  * one for codes that express their fields as MEDCoupling fields (ICoCo::MEDField) and one
  44  * for codes that express their fields as Trio/U fields.
  45  *
  46  * \section dec_options DEC Options
  47  * Options supported by DEC objects are
  48  *
  49  * <TABLE BORDER=1 >
  50  * <TR><TD>Option</TD><TD>Description</TD><TD>Default value</TD></TR>
  51  * <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>
  52  *</TABLE>
  53
  54
  55  The following code excerpt shows how to set options for an object that inherits from DEC :
  56
  57  \code
  58  IntersectionDEC dec(source_group,target_group);
  59  dec.setOptions("ForcedRenormalization",true);
  60  dec.attachLocalField(field);
  61  dec.synchronize();
  62  if (source_group.containsMyRank())
  63  dec.sendData();
  64  else
  65  dec.recvData();
  66  \endcode
  67 */
  68
  69 namespace ParaMEDMEM
  70 {
  71
  72
  73   /*! \addtogroup dec
  74     @{
  75   */
  76   DEC::DEC(ProcessorGroup& source_group, ProcessorGroup& target_group):_local_field(0),
  77                                                                        _source_group(&source_group),
  78                                                                        _target_group(&target_group),
  79                                                                        _owns_field(false),
  80                                                                        _icoco_field(0)
  81   {
  82     _union_group = source_group.fuse(target_group);
  83   }
  84
  85   DEC::~DEC()
  86   {
  87     //    delete _union_group;
  88     if(_owns_field)
  89       delete _local_field;
  90     delete _icoco_field;
  91     delete _union_group;
  92   }
  93   /*! Attaches a local field to a DEC.
  94     If the processor is on the receiving end of the DEC, the field
  95     will be updated by a recvData() call.
  96     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
  97   */
  98   void DEC::attachLocalField(const ParaFIELD* field)
  99   {
 100     _local_field=field;
 101     _comm_interface=&(field->getTopology()->getProcGroup()->getCommInterface());
 102     compareFieldAndMethod();
 103   }
 104
 105   /*! Attaches a local field to a DEC. The method will test whether the processor
 106     is on the source or the target side and will associate the mesh underlying the
 107     field to the local side.
 108
 109     If the processor is on the receiving end of the DEC, the field
 110     will be updated by a recvData() call.
 111     Reversely, if the processor is on the sending end, the field will be read, possibly transformed,
 112     and sent appropriately to the other side.
 113   */
 114
 115   void DEC::attachLocalField(MEDCouplingFieldDouble* field)
 116   {
 117     ProcessorGroup* local_group;
 118     if (_source_group->containsMyRank())
 119       local_group=_source_group;
 120     else if (_target_group->containsMyRank())
 121       local_group=_target_group;
 122     else
 123       throw INTERP_KERNEL::Exception("Invalid procgroup for field attachment to DEC");
 124
 125     _local_field= new ParaFIELD(field, *local_group);
 126     _owns_field=true;
 127     _comm_interface=&(local_group->getCommInterface());
 128   }
 129
 130   /*!
 131     Attaches a local field to a DEC.
 132     If the processor is on the receiving end of the DEC, the field
 133     will be updated by a recvData() call.
 134     Reversely, if the processor is on the sending end, the field will be read, possibly transformed, and sent appropriately to the other side.
 135     The field type is a generic ICoCo Field, so that the DEC can couple a number of different fields :
 136     - a ICoCo::MEDField, that is created from a MEDCoupling structure
 137     - a ICOCo::TrioField, that is created from tables extracted from a TRIO-U structure.
 138
 139   */
 140
 141   void DEC::attachLocalField(const ICoCo::Field* field)
 142   {
 143     const ICoCo::MEDField* medfield=dynamic_cast<const ICoCo::MEDField*> (field);
 144     if(medfield !=0)
 145       {
 146         attachLocalField(medfield->getField());
 147         return;
 148       }
 149     const ICoCo::TrioField* triofield=dynamic_cast<const ICoCo::TrioField*> (field);
 150     if (triofield !=0)
 151       {
 152         ProcessorGroup* localgroup;
 153         if (_source_group->containsMyRank())
 154           localgroup=_source_group;
 155         else
 156           localgroup=_target_group;
 157         _icoco_field=new ICoCo::MEDField(*const_cast<ICoCo::TrioField* >(triofield), *localgroup);
 158         attachLocalField(_icoco_field);
 159         return;
 160       }
 161     throw INTERP_KERNEL::Exception("incompatible field type");
 162   }
 163
 164   /*!
 165     Computes the field norm over its support
 166     on the source side and renormalizes the field on the target side
 167     so that the norms match.
 168
 169     \f[
 170     I_{source}=\sum_{i=1}^{n_{source}}V_{i}.|\Phi^{source}_{i}|^2,
 171     \f]
 172
 173     \f[
 174     I_{target}=\sum_{i=1}^{n_{target}}V_{i}.|\Phi^{target}_{i}|^2,
 175     \f]
 176
 177     \f[
 178     \Phi^{target}:=\Phi^{target}.\sqrt{I_{source}/I_{target}}.
 179     \f]
 180
 181   */
 182   void DEC::renormalizeTargetField()
 183   {
 184     if (_source_group->containsMyRank())
 185       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 186         {
 187           double total_norm = _local_field->getVolumeIntegral(icomp+1);
 188           double source_norm = total_norm;
 189           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 190
 191         }
 192     if (_target_group->containsMyRank())
 193       {
 194         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
 195           {
 196             double total_norm = _local_field->getVolumeIntegral(icomp+1);
 197             double source_norm=total_norm;
 198             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 199
 200             if (fabs(total_norm)>1e-100)
 201               _local_field->getField()->applyLin(source_norm/total_norm,0.0,icomp+1);
 202           }
 203       }
 204   }
 205   /*! @} */
 206
 207   void DEC::compareFieldAndMethod() const throw(INTERP_KERNEL::Exception)
 208   {
 209     if (_local_field)
 210       {
 211         TypeOfField entity = _local_field->getField()->getEntity();
 212         if ( getMethod() == "P0" )
 213           {
 214             if ( entity != ON_CELLS )
 215               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 216                                              " For P0 interpolation, field must be on MED_CELL's");
 217           }
 218         else if ( getMethod() == "P1" )
 219           {
 220             if ( entity != ON_NODES )
 221               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 222                                              " For P1 interpolation, field must be on MED_NODE's");
 223           }
 224         else if ( getMethod() == "P1d" )
 225           {
 226             if ( entity != ON_CELLS )
 227               throw INTERP_KERNEL::Exception("Field support and interpolation method mismatch."
 228                                              " For P1d interpolation, field must be on MED_CELL's");
 229             if ( _target_group->containsMyRank() )
 230               throw INTERP_KERNEL::Exception("Projection to P1d field not supported");
 231           }
 232         else
 233           {
 234             throw INTERP_KERNEL::Exception("Unknown interpolation method. Possible methods: P0, P1, P1d");
 235           }
 236       }
 237   }
 238 }