src/ParaMEDMEM/InterpKernelDEC.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 <mpi.h>
  21 #include "CommInterface.hxx"
  22 #include "Topology.hxx"
  23 #include "BlockTopology.hxx"
  24 #include "ComponentTopology.hxx"
  25 #include "ParaFIELD.hxx"
  26 #include "MPIProcessorGroup.hxx"
  27 #include "ParaMESH.hxx"
  28 #include "DEC.hxx"
  29 #include "InterpolationMatrix.hxx"
  30 #include "InterpKernelDEC.hxx"
  31 #include "ElementLocator.hxx"
  32
  33 namespace ParaMEDMEM
  34 {
  35
  36   /*!
  37     \anchor InterpKernelDEC-det
  38     \class InterpKernelDEC
  39
  40     \section InterpKernelDEC-over Overview
  41
  42     The InterpKernelDEC enables the \ref InterpKerRemapGlobal "remapping" (or interpolation) of fields between
  43     two parallel codes.
  44
  45     The projection
  46     methodology is based on the algorithms of %INTERP_KERNEL, that is to say, they work in a similar fashion than
  47     what the \ref remapper "sequential remapper" does. The following \ref discretization "projection methods"
  48     are supported: P0->P0 (the most common case), P1->P0, P0->P1.
  49
  50     The computation is possible for 3D meshes, 2D meshes, and 3D-surface
  51     meshes. Dimensions must be identical for code A and code B (for instance, though it could be
  52     desirable, it is not yet possible to couple 3D surfaces with 2D surfaces).
  53
  54     The name "InterpKernelDEC" comes from the fact that this class uses exactly the same algorithms
  55     as the sequential remapper. Both this class and the sequential
  56     \ref ParaMEDMEM::MEDCouplingRemapper "MEDCouplingRemapper" are built on top of the %INTERP_KERNEL
  57     algorithms (notably the computation of the intersection volumes).
  58
  59     Among the important properties inherited from the parent abstract class \ref DisjointDEC-det "DisjointDEC",
  60     the two \ref MPIProcessorGroup-det "processor groups" (source and target) must have a void intersection.
  61
  62     \image html NonCoincident_small.png "Transfer of a field supported by a quadrangular mesh to a triangular mesh".
  63
  64     \image latex NonCoincident_small.eps "Transfer of a field supported by a quadrangular mesh to a triangular mesh"
  65
  66     In the figure above we see the transfer of a field based on a quadrangular mesh to a new field supported by
  67     a triangular mesh. In a P0-P0 interpolation, to obtain the value on a triangle, the values on the
  68     quadrangles are weighted by their intersection area and summed.
  69
  70     A typical use of InterpKernelDEC encompasses two distinct phases :
  71     - A setup phase during which the intersection volumes are computed and the communication structures are
  72     setup. This corresponds to calling the InterpKernelDEC::synchronize() method.
  73     - A running phase during which the projections are actually performed. This corresponds to the calls to
  74     sendData() and recvData() which actually trigger the data exchange. The data exchange are synchronous
  75     in the current version of the library so that recvData() and sendData() calls must be synchronized
  76     on code A and code B processor groups.
  77
  78     The following code excerpt illustrates a typical use of the InterpKernelDEC class.
  79
  80     \code
  81     ...
  82     InterpKernelDEC dec(groupA, groupB);
  83     dec.attachLocalField(field);
  84     dec.synchronize();
  85     if (groupA.containsMyRank())
  86     dec.recvData();
  87     else if (groupB.containsMyRank())
  88     dec.sendData();
  89     ...
  90     \endcode
  91     A \ref InterpKerRemapGlobal "remapping" of the field from the source mesh to the target mesh is performed by
  92     the function synchronise(), which computes the interpolation matrix.
  93
  94     Computing the field on the receiving side can be expressed in terms of a matrix-vector product :
  95     \f$ \phi_t=W.\phi_s\f$, with \f$ \phi_t \f$ the field on the target side and \f$ \phi_s \f$ the field
  96     on the source side.
  97     When remapping a 3D surface to another 3D surface, a projection phase is necessary to match elements
  98     from both sides. Care must be taken when defining this projection to obtain a
  99     \ref InterpKerRemapGlobal "conservative remapping".
 100
 101     In the P0-P0 case, this matrix is a plain rectangular matrix with coefficients equal to the
 102     intersection areas between triangle and quadrangles. For instance, in the above figure, the matrix
 103     is :
 104
 105     \f[
 106     \begin{tabular}{|cccc|}
 107     0.72 & 0 & 0.2 & 0 \\
 108     0.46 & 0 & 0.51 & 0.03\\
 109     0.42 & 0.53 & 0 & 0.05\\
 110     0 & 0 & 0.92 & 0.05 \\
 111     \end{tabular}
 112     \f]
 113
 114     \section InterpKernelDEC-options Options
 115     On top of the usual \ref ParaMEDMEM::DECOptions "DEC options", the options supported by %InterpKernelDEC objects are
 116     related to the underlying \ref InterpKerIntersectors "intersector class".
 117     All the options available in the intersector objects are
 118     available for the %InterpKernelDEC object. The various options available for  intersectors can
 119     be reviewed in \ref InterpKerIntersectors.
 120
 121     For instance :
 122     \verbatim
 123     InterpKernelDEC dec(source_group, target_group);
 124     dec.attachLocalField(field);
 125     dec.setDoRotate(false);
 126     dec.setPrecision(1e-12);
 127     dec.synchronize();
 128     \endverbatim
 129
 130     \warning{  Options must be set before calling the synchronize method. }
 131   */
 132
 133   InterpKernelDEC::InterpKernelDEC():
 134     DisjointDEC(),
 135     _nb_distant_points(0), _distant_coords(0),
 136     _distant_locations(0), _interpolation_matrix(0)
 137   {
 138   }
 139
 140   /*!
 141     This constructor creates an InterpKernelDEC which has \a source_group as a working side
 142     and  \a target_group as an idle side. All the processors will actually participate, but intersection computations will be performed on the working side during the \a synchronize() phase.
 143     The constructor must be called synchronously on all processors of both processor groups.
 144
 145     \param source_group working side ProcessorGroup
 146     \param target_group lazy side ProcessorGroup
 147
 148   */
 149   InterpKernelDEC::InterpKernelDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
 150     DisjointDEC(source_group, target_group),
 151     _nb_distant_points(0), _distant_coords(0),
 152     _distant_locations(0), _interpolation_matrix(0)
 153   {
 154
 155   }
 156
 157   InterpKernelDEC::InterpKernelDEC(const std::set<int>& src_ids, const std::set<int>& trg_ids,
 158                                    const MPI_Comm& world_comm):
 159     DisjointDEC(src_ids,trg_ids,world_comm),
 160     _nb_distant_points(0), _distant_coords(0),
 161     _distant_locations(0), _interpolation_matrix(0)
 162   {
 163   }
 164
 165   InterpKernelDEC::~InterpKernelDEC()
 166   {
 167     if (_interpolation_matrix !=0)
 168       delete _interpolation_matrix;
 169   }
 170
 171   /*!
 172     \brief Synchronization process for exchanging topologies.
 173
 174     This method prepares all the structures necessary for sending data from a processor group to the other. It uses the mesh
 175     underlying the fields that have been set with attachLocalField method.
 176     It works in four steps :
 177     -# Bounding boxes are computed for each sub-domain,
 178     -# The lazy side mesh parts that are likely to intersect the working side local processor are sent to the working side,
 179     -# The working side calls the interpolation kernel to compute the intersection between local and imported mesh.
 180     -# The lazy side is updated so that it knows the structure of the data that will be sent by
 181     the working side during a \a sendData() call.
 182
 183   */
 184   void InterpKernelDEC::synchronize()
 185   {
 186     if(!isInUnion())
 187       return ;
 188     delete _interpolation_matrix;
 189     _interpolation_matrix = new InterpolationMatrix (_local_field, *_source_group,*_target_group,*this,*this);
 190
 191     //setting up the communication DEC on both sides
 192     if (_source_group->containsMyRank())
 193       {
 194         //locate the distant meshes
 195         ElementLocator locator(*_local_field, *_target_group, *_source_group);
 196         //transfering option from InterpKernelDEC to ElementLocator
 197         locator.copyOptions(*this);
 198         MEDCouplingPointSet* distant_mesh=0;
 199         int* distant_ids=0;
 200         std::string distantMeth;
 201         for (int i=0; i<_target_group->size(); i++)
 202           {
 203             //        int idistant_proc = (i+_source_group->myRank())%_target_group->size();
 204             int idistant_proc=i;
 205
 206             //gathers pieces of the target meshes that can intersect the local mesh
 207             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
 208             if (distant_mesh !=0)
 209               {
 210                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
 211                 //adds the contribution of the distant mesh on the local one
 212                 int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc);
 213                 //std::cout <<"add contribution from proc "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<std::endl;
 214                 _interpolation_matrix->addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth);
 215                 distant_mesh->decrRef();
 216                 delete [] distant_ids;
 217                 distant_mesh=0;
 218                 distant_ids=0;
 219               }
 220           }
 221        _interpolation_matrix->finishContributionW(locator);
 222       }
 223
 224     if (_target_group->containsMyRank())
 225       {
 226         ElementLocator locator(*_local_field, *_source_group, *_target_group);
 227         //transfering option from InterpKernelDEC to ElementLocator
 228         locator.copyOptions(*this);
 229         MEDCouplingPointSet* distant_mesh=0;
 230         int* distant_ids=0;
 231         for (int i=0; i<_source_group->size(); i++)
 232           {
 233             //        int idistant_proc = (i+_target_group->myRank())%_source_group->size();
 234             int  idistant_proc=i;
 235             //gathers pieces of the target meshes that can intersect the local mesh
 236             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
 237             //std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<std::endl;
 238             if (distant_mesh!=0)
 239               {
 240                 std::string distantMeth;
 241                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
 242                 distant_mesh->decrRef();
 243                 delete [] distant_ids;
 244                 distant_mesh=0;
 245                 distant_ids=0;
 246               }
 247           }
 248         _interpolation_matrix->finishContributionL(locator);
 249       }
 250     _interpolation_matrix->prepare();
 251   }
 252
 253
 254   /*!
 255     Receives the data whether the processor is on the working side or on the lazy side. It must match a \a sendData() call on the other side.
 256   */
 257   void InterpKernelDEC::recvData()
 258   {
 259     if (_source_group->containsMyRank())
 260       _interpolation_matrix->transposeMultiply(*_local_field->getField());
 261     else if (_target_group->containsMyRank())
 262       {
 263         _interpolation_matrix->multiply(*_local_field->getField());
 264         if (getForcedRenormalization())
 265           renormalizeTargetField(getMeasureAbsStatus());
 266       }
 267   }
 268
 269
 270   /*!
 271     Receives the data at time \a time in asynchronous mode. The value of the field
 272     will be time-interpolated from the field values received.
 273     \param time time at which the value is desired
 274   */
 275   void InterpKernelDEC::recvData( double time )
 276   {
 277     _interpolation_matrix->getAccessDEC()->setTime(time);
 278     recvData() ;
 279   }
 280
 281   /*!
 282     Sends the data whether the processor is on the working side or on the lazy side.
 283     It must match a recvData() call on the other side.
 284   */
 285   void InterpKernelDEC::sendData()
 286   {
 287     if (_source_group->containsMyRank())
 288       {
 289
 290         _interpolation_matrix->multiply(*_local_field->getField());
 291         if (getForcedRenormalization())
 292           renormalizeTargetField(getMeasureAbsStatus());
 293
 294       }
 295     else if (_target_group->containsMyRank())
 296       _interpolation_matrix->transposeMultiply(*_local_field->getField());
 297   }
 298
 299   /*!
 300     Sends the data available at time \a time in asynchronous mode.
 301     \param time time at which the value is available
 302     \param deltatime time interval between the value presently sent and the next one.
 303   */
 304   void InterpKernelDEC::sendData( double time , double deltatime )
 305   {
 306     _interpolation_matrix->getAccessDEC()->setTime(time,deltatime);
 307     sendData() ;
 308   }
 309
 310 }