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