medtool/src/ParaMEDMEM/ParaFIELD.cxx

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