src/ParaMEDMEM/ParaFIELD.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 "Topology.hxx"
  20 #include "BlockTopology.hxx"
  21 #include "ComponentTopology.hxx"
  22 #include "ExplicitCoincidentDEC.hxx"
  23 #include "StructuredCoincidentDEC.hxx"
  24 #include "CommInterface.hxx"
  25 #include "ProcessorGroup.hxx"
  26 #include "MPIProcessorGroup.hxx"
  27 #include "ParaFIELD.hxx"
  28 #include "ParaMESH.hxx"
  29 #include "InterpKernelUtilities.hxx"
  30 #include "InterpolationMatrix.hxx"
  31
  32 #include <numeric>
  33
  34 namespace ParaMEDMEM
  35 {
  36   /*!
  37     \defgroup parafield ParaFIELD
  38     This class encapsulates parallel fields. It basically encapsulates
  39     a MEDCouplingField with extra information related to parallel
  40     topology.
  41     It is most conveniently created by giving a pointer to a MEDCouplingField
  42     object and a \c ProcessorGroup.
  43     By default, a ParaFIELD object will be constructed with all field components
  44     located on the same processors. In some specific cases, it might be necessary to scatter components over several processors. In this case, the constructor
  45     using a ComponentTopology is required.
  46
  47     @{ */
  48
  49   /*!
  50
  51     \brief  Constructing a \c ParaFIELD from a \c ParaSUPPORT and a \c ComponentTopology.
  52
  53     This constructor creates an empty field based on the ParaSUPPORT description
  54     and the partitioning of components described in \a component_topology.
  55     It takes ownership over the \c _field object that it creates.
  56
  57     Here come the three ComponentTopology constructors :
  58     \verbatim
  59     ComponentTopology c; // one component in the field
  60     ComponentTopology c(6); //six components, all of them on the same processor
  61     ComponentTopology c(6, proc_group); // six components, evenly distributed over the processors of procgroup
  62     \endverbatim
  63
  64   */
  65   ParaFIELD::ParaFIELD(TypeOfField type, ParaMESH* para_support, const ComponentTopology& component_topology)
  66     :_support(para_support),
  67      _component_topology(component_topology),_topology(0),
  68      _field(0),
  69      _has_field_ownership(true),
  70      _has_support_ownership(false)
  71   {
  72     if (para_support->isStructured() || (para_support->getTopology()->getProcGroup()->size()==1 && component_topology.nbBlocks()!=1))
  73       {
  74         const BlockTopology* source_topo = dynamic_cast<const BlockTopology*>(para_support->getTopology());
  75         _topology=new BlockTopology(*source_topo,component_topology);
  76       }
  77     else
  78       {
  79         if (component_topology.nbBlocks()!=1 &&  para_support->getTopology()->getProcGroup()->size()!=1)
  80           throw INTERP_KERNEL::Exception(LOCALIZED("ParaFIELD constructor : Unstructured Support not taken into account with component topology yet"));
  81         else
  82           {
  83             const BlockTopology* source_topo=dynamic_cast<const BlockTopology*> (para_support->getTopology());
  84             int nb_local_comp=component_topology.nbLocalComponents();
  85             _topology=new BlockTopology(*source_topo,nb_local_comp);
  86           }
  87       }
  88     int nb_components = component_topology.nbLocalComponents();
  89     if (nb_components!=0)
  90       {
  91         _field=MEDCouplingFieldDouble::New(type);
  92         _field->setMesh(_support->getCellMesh());
  93         DataArrayDouble *array=DataArrayDouble::New();
  94         array->alloc(_field->getNumberOfTuples(),nb_components);
  95         _field->setArray(array);
  96         array->decrRef();
  97       }
  98     else return;
  99
 100     _field->setName("Default ParaFIELD name");
 101     _field->setDescription("Default ParaFIELD description");
 102     _field->setDtIt(0,0);
 103     _field->setTime(0.0);
 104   }
 105
 106   /*! \brief Constructor creating the ParaFIELD
 107     from a given FIELD and a processor group.
 108
 109     This constructor supposes that support underlying \a subdomain_field has no ParaSUPPORT
 110     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).
 111   */
 112   ParaFIELD::ParaFIELD(MEDCouplingFieldDouble* subdomain_field, const ProcessorGroup& proc_group):
 113     _field(subdomain_field),
 114     _support(),
 115     _component_topology(ComponentTopology(_field->getNumberOfComponents())),_topology(0),
 116     _has_field_ownership(false),
 117     _has_support_ownership(true)
 118   {
 119     if(_field)
 120       _field->incrRef();
 121     const ExplicitTopology* source_topo=dynamic_cast<const ExplicitTopology*> (_support->getTopology());
 122     _topology=new ExplicitTopology(*source_topo,_component_topology.nbLocalComponents());
 123   }
 124
 125   ParaFIELD::~ParaFIELD()
 126   {
 127     if(_field)
 128       _field->decrRef();
 129     delete _topology;
 130   }
 131
 132   void ParaFIELD::synchronizeTarget(ParaFIELD* source_field)
 133   {
 134     DEC* data_channel;
 135     if (dynamic_cast<BlockTopology*>(_topology)!=0)
 136       {
 137         data_channel=new StructuredCoincidentDEC;
 138       }
 139     else
 140       {
 141         data_channel=new ExplicitCoincidentDEC;
 142       }
 143     data_channel->attachLocalField(this);
 144     data_channel->synchronize();
 145     data_channel->prepareTargetDE();
 146     data_channel->recvData();
 147
 148     delete data_channel;
 149   }
 150
 151   void ParaFIELD::synchronizeSource(ParaFIELD* target_field)
 152   {
 153     DEC* data_channel;
 154     if (dynamic_cast<BlockTopology*>(_topology)!=0)
 155       {
 156         data_channel=new StructuredCoincidentDEC;
 157       }
 158     else
 159       {
 160         data_channel=new ExplicitCoincidentDEC;
 161       }
 162     data_channel->attachLocalField(this);
 163     data_channel->synchronize();
 164     data_channel->prepareSourceDE();
 165     data_channel->sendData();
 166
 167     delete data_channel;
 168   }
 169
 170   /*! This method retrieves the integral of component \a icomp
 171     over the all domain. */
 172   double ParaFIELD::getVolumeIntegral(int icomp) const
 173   {
 174     CommInterface comm_interface = _topology->getProcGroup()->getCommInterface();
 175     MEDCouplingFieldDouble *volume=InterpolationMatrix::getSupportVolumes(_field->getMesh());
 176     double *ptr=volume->getArray()->getPointer();
 177     int nbOfValues=volume->getArray()->getNbOfElems();
 178     double integral=0.;
 179     for (int i=0; i<nbOfValues; i++)
 180       integral+=_field->getIJ(i,icomp)*ptr[i];
 181
 182     volume->decrRef();
 183
 184     double total=0.;
 185     const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*>(_topology->getProcGroup()))->getComm();
 186     comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, *comm);
 187
 188     return total;
 189   }
 190 }
 191