src/ParaMEDMEM/InterpKernelDEC.cxx

   1 // Copyright (C) 2007-2020  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 MEDCoupling
  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 MEDCoupling::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 MEDCoupling::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     The source group and target group MUST form a partition of all the procs within the communicator passed as 'world_comm'
 145     when building the group.
 146
 147     \param source_group working side ProcessorGroup
 148     \param target_group lazy side ProcessorGroup
 149
 150   */
 151   InterpKernelDEC::InterpKernelDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
 152     DisjointDEC(source_group, target_group),
 153     _nb_distant_points(0), _distant_coords(0),
 154     _distant_locations(0), _interpolation_matrix(0)
 155   {
 156
 157   }
 158
 159
 160   /*!
 161    * Creates an InterpKernelDEC from a set of source procs IDs and target group IDs.
 162    * The difference with the ctor using groups is that the set of procs might not cover entirely MPI_COMM_WORLD
 163    * (a sub-communicator holding the union of source and target procs is recreated internally).
 164    */
 165   InterpKernelDEC::InterpKernelDEC(const std::set<int>& src_ids, const std::set<int>& trg_ids,
 166                                    const MPI_Comm& world_comm):
 167     DisjointDEC(src_ids,trg_ids,world_comm),
 168     _nb_distant_points(0), _distant_coords(0),
 169     _distant_locations(0), _interpolation_matrix(0)
 170   {
 171   }
 172
 173   InterpKernelDEC::~InterpKernelDEC()
 174   {
 175     if (_interpolation_matrix !=0)
 176       delete _interpolation_matrix;
 177   }
 178
 179   /*!
 180     \brief Synchronization process for exchanging topologies.
 181
 182     This method prepares all the structures necessary for sending data from a processor group to the other. It uses the mesh
 183     underlying the fields that have been set with attachLocalField method.
 184     It works in four steps :
 185     -# Bounding boxes are computed for each sub-domain,
 186     -# The lazy side mesh parts that are likely to intersect the working side local processor are sent to the working side,
 187     -# The working side calls the interpolation kernel to compute the intersection between local and imported mesh.
 188     -# The lazy side is updated so that it knows the structure of the data that will be sent by
 189     the working side during a \a sendData() call.
 190
 191   */
 192   void InterpKernelDEC::synchronize()
 193   {
 194     if(!isInUnion())
 195       return ;
 196     delete _interpolation_matrix;
 197     _interpolation_matrix = new InterpolationMatrix (_local_field, *_source_group,*_target_group,*this,*this);
 198
 199     //setting up the communication DEC on both sides
 200     if (_source_group->containsMyRank())
 201       {
 202         //locate the distant meshes
 203         ElementLocator locator(*_local_field, *_target_group, *_source_group);
 204         //transferring option from InterpKernelDEC to ElementLocator
 205         locator.copyOptions(*this);
 206         MEDCouplingPointSet* distant_mesh=0;
 207         mcIdType* distant_ids=0;
 208         std::string distantMeth;
 209         for (int i=0; i<_target_group->size(); i++)
 210           {
 211             //        int idistant_proc = (i+_source_group->myRank())%_target_group->size();
 212             int idistant_proc=i;
 213
 214             //gathers pieces of the target meshes that can intersect the local mesh
 215             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
 216             if (distant_mesh !=0)
 217               {
 218                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
 219                 //adds the contribution of the distant mesh on the local one
 220                 int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc);
 221                 //std::cout <<"add contribution from proc "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<std::endl;
 222                 _interpolation_matrix->addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth);
 223                 distant_mesh->decrRef();
 224                 delete [] distant_ids;
 225                 distant_mesh=0;
 226                 distant_ids=0;
 227               }
 228           }
 229        _interpolation_matrix->finishContributionW(locator);
 230       }
 231
 232     if (_target_group->containsMyRank())
 233       {
 234         ElementLocator locator(*_local_field, *_source_group, *_target_group);
 235         //transferring option from InterpKernelDEC to ElementLocator
 236         locator.copyOptions(*this);
 237         MEDCouplingPointSet* distant_mesh=0;
 238         mcIdType* distant_ids=0;
 239         for (int i=0; i<_source_group->size(); i++)
 240           {
 241             //        int idistant_proc = (i+_target_group->myRank())%_source_group->size();
 242             int  idistant_proc=i;
 243             //gathers pieces of the target meshes that can intersect the local mesh
 244             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
 245             //std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<std::endl;
 246             if (distant_mesh!=0)
 247               {
 248                 std::string distantMeth;
 249                 locator.exchangeMethod(_method,idistant_proc,distantMeth);
 250                 distant_mesh->decrRef();
 251                 delete [] distant_ids;
 252                 distant_mesh=0;
 253                 distant_ids=0;
 254               }
 255           }
 256         _interpolation_matrix->finishContributionL(locator);
 257       }
 258     _interpolation_matrix->prepare();
 259   }
 260
 261
 262   /*!
 263     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.
 264   */
 265   void InterpKernelDEC::recvData()
 266   {
 267     if (_source_group->containsMyRank())
 268       _interpolation_matrix->transposeMultiply(*_local_field->getField());
 269     else if (_target_group->containsMyRank())
 270       {
 271         _interpolation_matrix->multiply(*_local_field->getField());
 272         if (getForcedRenormalization())
 273           renormalizeTargetField(getMeasureAbsStatus());
 274       }
 275   }
 276
 277
 278   /*!
 279     Receives the data at time \a time in asynchronous mode. The value of the field
 280     will be time-interpolated from the field values received.
 281     \param time time at which the value is desired
 282   */
 283   void InterpKernelDEC::recvData( double time )
 284   {
 285     _interpolation_matrix->getAccessDEC()->setTime(time);
 286     recvData() ;
 287   }
 288
 289   /*!
 290     Sends the data whether the processor is on the working side or on the lazy side.
 291     It must match a recvData() call on the other side.
 292   */
 293   void InterpKernelDEC::sendData()
 294   {
 295     if (_source_group->containsMyRank())
 296       {
 297
 298         _interpolation_matrix->multiply(*_local_field->getField());
 299         if (getForcedRenormalization())
 300           renormalizeTargetField(getMeasureAbsStatus());
 301
 302       }
 303     else if (_target_group->containsMyRank())
 304       _interpolation_matrix->transposeMultiply(*_local_field->getField());
 305   }
 306
 307   /*!
 308     Sends the data available at time \a time in asynchronous mode.
 309     \param time time at which the value is available
 310     \param deltatime time interval between the value presently sent and the next one.
 311   */
 312   void InterpKernelDEC::sendData( double time , double deltatime )
 313   {
 314     _interpolation_matrix->getAccessDEC()->setTime(time,deltatime);
 315     sendData() ;
 316   }
 317
 318 }