src/ParaMEDMEM/InterpKernelDEC.cxx

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