src/ParaMEDMEM/IntersectionDEC.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 <mpi.h>
  20 #include "CommInterface.hxx"
  21 #include "Topology.hxx"
  22 #include "BlockTopology.hxx"
  23 #include "ComponentTopology.hxx"
  24 #include "ParaFIELD.hxx"
  25 #include "MPIProcessorGroup.hxx"
  26 #include "ParaMESH.hxx"
  27 #include "DEC.hxx"
  28 #include "InterpolationMatrix.hxx"
  29 #include "IntersectionDEC.hxx"
  30 #include "ElementLocator.hxx"
  31
  32
  33
  34 namespace ParaMEDMEM
  35 {
  36
  37   /*!
  38     \defgroup intersectiondec IntersectionDEC
  39
  40     \section overview Overview
  41
  42     The IntersectionDEC 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).
  43
  44     In the present version, only fields lying on elements are considered.
  45
  46     \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."
  47
  48     \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."
  49
  50     A typical use of IntersectionDEC encompasses two distinct phases :
  51     - A setup phase during which the intersection volumes are computed and the communication structures are setup. This corresponds to calling the IntersectionDEC::synchronize() method.
  52     - 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.
  53
  54     The following code excerpt illutrates a typical use of the IntersectionDEC class.
  55
  56     \code
  57     ...
  58     IntersectionDEC dec(groupA, groupB);
  59     dec.attachLocalField(field);
  60     dec.synchronize();
  61     if (groupA.containsMyRank())
  62     dec.recvData();
  63     else if (groupB.containsMyRank())
  64     dec.sendData();
  65     ...
  66     \endcode
  67     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.
  68
  69     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.
  70     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.
  71
  72     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 :
  73
  74     \f[
  75     \begin{tabular}{|cccc|}
  76     0.72 & 0 & 0.2 & 0 \\
  77     0.46 & 0 & 0.51 & 0.03\\
  78     0.42 & 0.53 & 0 & 0.05\\
  79     0 & 0 & 0.92 & 0.05 \\
  80     \end{tabular}
  81     \f]
  82
  83
  84
  85     \section intersectiondec_options Options
  86     On top of \ref dec_options, options supported by %IntersectionDEC objects are
  87     related to the underlying Intersector class.
  88     All the options available in the intersector objects are
  89     available for the %IntersectionDEC object. The various options available for  * intersectors can be reviewed in \ref InterpKerIntersectors.
  90
  91     For instance :
  92     \verbatim
  93     IntersectionDEC dec(source_group, target_group);
  94     dec.attachLocalField(field);
  95     dec.setOptions("DoRotate",false);
  96     dec.setOptions("Precision",1e-12);
  97     dec.synchronize();
  98     \endverbatim
  99
 100     \warning{  Options must be set before calling the synchronize method. }
 101   */
 102
 103   /*!
 104     \addtogroup intersectiondec
 105     @{
 106   */
 107
 108   IntersectionDEC::IntersectionDEC()
 109   {
 110   }
 111
 112   /*!
 113     This constructor creates an IntersectionDEC which has \a source_group as a working side
 114     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.
 115     The constructor must be called synchronously on all processors of both processor groups.
 116
 117     \param source_group working side ProcessorGroup
 118     \param target_group lazy side ProcessorGroup
 119
 120   */
 121   IntersectionDEC::IntersectionDEC(ProcessorGroup& source_group, ProcessorGroup& target_group):
 122     DEC(source_group, target_group),_interpolation_matrix(0)
 123   {
 124
 125   }
 126
 127   IntersectionDEC::~IntersectionDEC()
 128   {
 129     if (_interpolation_matrix !=0)
 130       delete _interpolation_matrix;
 131   }
 132
 133   /*!
 134     \brief Synchronization process for exchanging topologies.
 135
 136     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.
 137     It works in four steps :
 138     -# Bounding boxes are computed for each subdomain,
 139     -# The lazy side mesh parts that are likely to intersect the working side local processor are sent to the working side,
 140     -# The working side calls the interpolation kernel to compute the intersection between local and imported mesh.
 141     -# The lazy side is updated so that it knows the structure of the data that will be sent by
 142     the working side during a \a sendData() call.
 143
 144   */
 145   void IntersectionDEC::synchronize()
 146   {
 147     ParaMEDMEM::ParaMESH* para_mesh = _local_field->getSupport();
 148     //cout <<"size of Interpolation Matrix"<<sizeof(InterpolationMatrix)<<endl;
 149     _interpolation_matrix = new InterpolationMatrix (para_mesh, *_source_group,*_target_group,*this,*this);
 150
 151     //setting up the communication DEC on both sides
 152     if (_source_group->containsMyRank())
 153       {
 154         //locate the distant meshes
 155         ElementLocator locator(*para_mesh, *_target_group);
 156
 157         //transfering option from IntersectionDEC to ElementLocator
 158         locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment());
 159
 160         MEDCouplingUMesh* distant_mesh=0;
 161         int* distant_ids=0;
 162         for (int i=0; i<_target_group->size(); i++)
 163           {
 164             //        int idistant_proc = (i+_source_group->myRank())%_target_group->size();
 165             int       idistant_proc=i;
 166
 167             //gathers pieces of the target meshes that can intersect the local mesh
 168             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
 169             std::string distantMeth;
 170             locator.exchangeMethod(_method,idistant_proc,distantMeth);
 171             if (distant_mesh !=0)
 172               {
 173                 //adds the contribution of the distant mesh on the local one
 174                 int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc);
 175                 std::cout <<"add contribution from proc "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<std::endl;
 176                 _interpolation_matrix->addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth);
 177
 178                 distant_mesh->decrRef();
 179                 delete[] distant_ids;
 180                 distant_mesh=0;
 181                 distant_ids=0;
 182               }
 183           }
 184       }
 185
 186     if (_target_group->containsMyRank())
 187       {
 188         ElementLocator locator(*para_mesh, *_source_group);
 189         //transfering option from IntersectionDEC to ElementLocator
 190         locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment());
 191
 192         MEDCouplingUMesh* distant_mesh=0;
 193         int* distant_ids=0;
 194         for (int i=0; i<_source_group->size(); i++)
 195           {
 196             //        int idistant_proc = (i+_target_group->myRank())%_source_group->size();
 197             int  idistant_proc=i;
 198             //gathers pieces of the target meshes that can intersect the local mesh
 199             locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
 200             std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<std::endl;
 201             std::string distantMeth;
 202             locator.exchangeMethod(_method,idistant_proc,distantMeth);
 203             if (distant_mesh!=0)
 204               {
 205                 distant_mesh->decrRef();
 206                 delete[] distant_ids;
 207                 distant_mesh=0;
 208                 distant_ids=0;
 209               }
 210           }
 211       }
 212     _interpolation_matrix->prepare();
 213   }
 214
 215
 216   /*!
 217     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.
 218   */
 219   void IntersectionDEC::recvData()
 220   {
 221     if (_source_group->containsMyRank())
 222       _interpolation_matrix->transposeMultiply(*_local_field->getField());
 223     else if (_target_group->containsMyRank())
 224       {
 225         _interpolation_matrix->multiply(*_local_field->getField());
 226         if (getForcedRenormalization())
 227           renormalizeTargetField();
 228       }
 229   }
 230
 231
 232   /*!
 233     Receives the data at time \a time in asynchronous mode. The value of the field
 234     will be time-interpolated from the field values received.
 235     \param time time at which the value is desired
 236   */
 237   void IntersectionDEC::recvData( double time )
 238   {
 239     _interpolation_matrix->getAccessDEC()->setTime(time);
 240     recvData() ;
 241   }
 242
 243   /*!
 244     Sends the data whether the processor is on the working side or on the lazy side.
 245     It must match a recvData() call on the other side.
 246   */
 247   void IntersectionDEC::sendData()
 248   {
 249     if (_source_group->containsMyRank())
 250       {
 251
 252         _interpolation_matrix->multiply(*_local_field->getField());
 253         if (getForcedRenormalization())
 254           renormalizeTargetField();
 255
 256       }
 257     else if (_target_group->containsMyRank())
 258       _interpolation_matrix->transposeMultiply(*_local_field->getField());
 259   }
 260
 261   /*!
 262     Sends the data available at time \a time in asynchronous mode.
 263     \param time time at which the value is available
 264     \param deltatime time interval between the value presently sent and the next one.
 265   */
 266   void IntersectionDEC::sendData( double time , double deltatime )
 267   {
 268     _interpolation_matrix->getAccessDEC()->setTime(time,deltatime);
 269     sendData() ;
 270   }
 271
 272   /*!
 273     @}
 274   */
 275
 276 }