src/ParaMEDMEM/ParaFIELD.cxx

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