src/ParaMEDMEM/ParaFIELD.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 "Topology.hxx"
  21 #include "BlockTopology.hxx"
  22 #include "ComponentTopology.hxx"
  23 #include "ExplicitCoincidentDEC.hxx"
  24 #include "StructuredCoincidentDEC.hxx"
  25 #include "CommInterface.hxx"
  26 #include "ProcessorGroup.hxx"
  27 #include "MPIProcessorGroup.hxx"
  28 #include "ParaFIELD.hxx"
  29 #include "ParaMESH.hxx"
  30 #include "InterpKernelUtilities.hxx"
  31 #include "InterpolationMatrix.hxx"
  32
  33 #include <numeric>
  34
  35 namespace ParaMEDMEM
  36 {
  37   /*!
  38     \defgroup parafield ParaFIELD
  39     This class encapsulates parallel fields. It basically encapsulates
  40     a MEDCouplingField with extra information related to parallel
  41     topology.
  42     It is most conveniently created by giving a pointer to a MEDCouplingField
  43     object and a \c ProcessorGroup.
  44     By default, a ParaFIELD object will be constructed with all field components
  45     located on the same processors. In some specific cases, it might be necessary to scatter components over several processors. In this case, the constructor
  46     using a ComponentTopology is required.
  47
  48     @{ */
  49
  50   /*!
  51
  52     \brief  Constructing a \c ParaFIELD from a \c ParaSUPPORT and a \c ComponentTopology.
  53
  54     This constructor creates an empty field based on the ParaSUPPORT description
  55     and the partitioning of components described in \a component_topology.
  56     It takes ownership over the \c _field object that it creates.
  57
  58     Here come the three ComponentTopology constructors :
  59     \verbatim
  60     ComponentTopology c; // one component in the field
  61     ComponentTopology c(6); //six components, all of them on the same processor
  62     ComponentTopology c(6, proc_group); // six components, evenly distributed over the processors of procgroup
  63     \endverbatim
  64
  65   */
  66   ParaFIELD::ParaFIELD(TypeOfField type, TypeOfTimeDiscretization td, ParaMESH* para_support, const ComponentTopology& component_topology)
  67     :_field(0),
  68      _component_topology(component_topology),_topology(0),_own_support(false),
  69      _support(para_support)
  70   {
  71     if (para_support->isStructured() || (para_support->getTopology()->getProcGroup()->size()==1 && component_topology.nbBlocks()!=1))
  72       {
  73         const BlockTopology* source_topo = dynamic_cast<const BlockTopology*>(para_support->getTopology());
  74         _topology=new BlockTopology(*source_topo,component_topology);
  75       }
  76     else
  77       {
  78         if (component_topology.nbBlocks()!=1 &&  para_support->getTopology()->getProcGroup()->size()!=1)
  79           throw INTERP_KERNEL::Exception(LOCALIZED("ParaFIELD constructor : Unstructured Support not taken into account with component topology yet"));
  80         else
  81           {
  82             const BlockTopology* source_topo=dynamic_cast<const BlockTopology*> (para_support->getTopology());
  83             int nb_local_comp=component_topology.nbLocalComponents();
  84             _topology=new BlockTopology(*source_topo,nb_local_comp);
  85           }
  86       }
  87     int nb_components = component_topology.nbLocalComponents();
  88     if (nb_components!=0)
  89       {
  90         _field=MEDCouplingFieldDouble::New(type,td);
  91         _field->setMesh(_support->getCellMesh());
  92         DataArrayDouble *array=DataArrayDouble::New();
  93         array->alloc(_field->getNumberOfTuples(),nb_components);
  94         _field->setArray(array);
  95         array->decrRef();
  96       }
  97     else return;
  98
  99     _field->setName("Default ParaFIELD name");
 100     _field->setDescription("Default ParaFIELD description");
 101   }
 102
 103   /*! \brief Constructor creating the ParaFIELD
 104     from a given FIELD and a processor group.
 105
 106     This constructor supposes that support underlying \a subdomain_field has no ParaSUPPORT
 107     attached and it therefore recreates one. It therefore takes ownership over _support. The component topology associated with the field is a basic one (all components on the same processor).
 108   */
 109   ParaFIELD::ParaFIELD(MEDCouplingFieldDouble* subdomain_field, ParaMESH *sup, const ProcessorGroup& proc_group):
 110     _field(subdomain_field),
 111     _component_topology(ComponentTopology(_field->getNumberOfComponents())),_topology(0),_own_support(false),
 112     _support(sup)
 113   {
 114     if(_field)
 115       _field->incrRef();
 116     const BlockTopology* source_topo=dynamic_cast<const BlockTopology*> (_support->getTopology());
 117     _topology=new BlockTopology(*source_topo,_component_topology.nbLocalComponents());
 118   }
 119
 120   ParaFIELD::~ParaFIELD()
 121   {
 122     if(_field)
 123       _field->decrRef();
 124     if(_own_support)
 125       delete _support;
 126     delete _topology;
 127   }
 128
 129   void ParaFIELD::synchronizeTarget(ParaFIELD* source_field)
 130   {
 131     DisjointDEC* data_channel;
 132     if (dynamic_cast<BlockTopology*>(_topology)!=0)
 133       {
 134         data_channel=new StructuredCoincidentDEC;
 135       }
 136     else
 137       {
 138         data_channel=new ExplicitCoincidentDEC;
 139       }
 140     data_channel->attachLocalField(this);
 141     data_channel->synchronize();
 142     data_channel->prepareTargetDE();
 143     data_channel->recvData();
 144
 145     delete data_channel;
 146   }
 147
 148   void ParaFIELD::synchronizeSource(ParaFIELD* target_field)
 149   {
 150     DisjointDEC* data_channel;
 151     if (dynamic_cast<BlockTopology*>(_topology)!=0)
 152       {
 153         data_channel=new StructuredCoincidentDEC;
 154       }
 155     else
 156       {
 157         data_channel=new ExplicitCoincidentDEC;
 158       }
 159     data_channel->attachLocalField(this);
 160     data_channel->synchronize();
 161     data_channel->prepareSourceDE();
 162     data_channel->sendData();
 163
 164     delete data_channel;
 165   }
 166
 167   /*!
 168    * This method returns, if it exists, an array with only one component and as many as tuples as _field has.
 169    * This array gives for every element on which this->_field lies, its global number, if this->_field is nodal.
 170    * For example if _field is a nodal field : returned array will be the nodal global numbers.
 171    * The content of this method is used to inform Working side to accumulate data recieved by lazy side.
 172    */
 173   DataArrayInt* ParaFIELD::returnCumulativeGlobalNumbering() const
 174   {
 175     if(!_field)
 176       return 0;
 177     TypeOfField type=_field->getTypeOfField();
 178     switch(type)
 179       {
 180       case ON_CELLS:
 181         return 0;
 182       case ON_NODES:
 183         return _support->getGlobalNumberingNodeDA();
 184       default:
 185         return 0;
 186       }
 187   }
 188
 189   DataArrayInt* ParaFIELD::returnGlobalNumbering() const
 190   {
 191     if(!_field)
 192       return 0;
 193     TypeOfField type=_field->getTypeOfField();
 194     switch(type)
 195       {
 196       case ON_CELLS:
 197         return _support->getGlobalNumberingCellDA();
 198       case ON_NODES:
 199         return _support->getGlobalNumberingNodeDA();
 200       default:
 201         return 0;
 202       }
 203   }
 204
 205   int ParaFIELD::nbComponents() const
 206   {
 207     return _component_topology.nbComponents();
 208   }
 209
 210
 211   /*! This method retrieves the integral of component \a icomp
 212     over the all domain. */
 213   double ParaFIELD::getVolumeIntegral(int icomp, bool isWAbs) const
 214   {
 215     CommInterface comm_interface = _topology->getProcGroup()->getCommInterface();
 216     double integral=_field->integral(icomp,isWAbs);
 217     double total=0.;
 218     const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*>(_topology->getProcGroup()))->getComm();
 219     comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, *comm);
 220
 221     return total;
 222   }
 223 }