src/ParaMEDMEM/NonCoincidentDEC.cxx

   1 // Copyright (C) 2007-2016  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 "DEC.hxx"
  28 #include "NonCoincidentDEC.hxx"
  29
  30 extern "C" {
  31 #include <fvm_parall.h>
  32 #include <fvm_nodal.h>
  33 #include <fvm_nodal_append.h>
  34 #include <fvm_locator.h>
  35 }
  36
  37 namespace MEDCoupling
  38 {
  39
  40   /*!
  41     \anchor NonCoincidentDEC-det
  42     \class NonCoincidentDEC
  43
  44     \c NonCoincidentDEC enables non-conservative remapping of fields
  45     between two parallel codes.
  46     The computation is possible for 3D meshes and 2D meshes.
  47     It is not available for 3D surfaces.
  48
  49     The computation enables fast parallel localization, and is based on a point in element search, followed
  50     by a field evaluation at the point location. Thus, it is typically
  51     faster than the \ref InterpKernelDEC-det "InterpKernelDEC" which uses a
  52     \ref InterpKerRemapGlobal "conservative remapping" (i.e. the same algorithms of volume
  53     intersection as in the \ref remapper "sequential remapper")
  54     It is particularly true for the initialisation phase (synchronize() method)
  55     which has a significant computation cost in \ref InterpKernelDEC-det.
  56
  57     In the present version, only fields lying on elements ("P0") are considered.
  58     The value is estimated by locating the barycenter of the target
  59     side cell in a source cell and sending the value of this source cell
  60     as the value of the target cell.
  61
  62     \image html NonCoincident_small.png "Example showing the transfer from a field based on a quadrangular mesh to a triangular mesh. The triangle barycenters are computed and located in the quadrangles. In a P0-P0 interpolation, the value on the quadrangle is then applied to the triangles whose barycenter lies within."
  63
  64     \image latex NonCoincident_small.eps "Example showing the transfer from a field based on a quadrangular mesh to a triangular mesh. The triangle barycenters are computed and located in the quadrangles. In a P0-P0 interpolation, the value on the quadrangle is then applied to the triangles whose barycenter lies within."
  65
  66     A typical use of NonCoincidentDEC encompasses two distinct phases :
  67     - A setup phase during which the intersection volumes are computed and the communication structures are setup. This corresponds to calling the NonCoincidentDEC::synchronize() method.
  68     - 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.
  69
  70     The following code excerpt illutrates a typical use of the NonCoincidentDEC class.
  71
  72     \code
  73     ...
  74     NonCoincidentDEC dec(groupA, groupB);
  75     dec.attachLocalField(field);
  76     dec.synchronize();
  77     if (groupA.containsMyRank())
  78     dec.recvData();
  79     else if (groupB.containsMyRank())
  80     dec.sendData();
  81     ...
  82     \endcode
  83
  84     Computing the field on the receiving side can be expressed in terms
  85     of a matrix-vector product : \f$ \phi_t=W.\phi_s\f$, with \f$ \phi_t
  86     \f$ the field on the target side and \f$ \phi_s \f$ the field on
  87     the source side.
  88     In the P0-P0 case, this matrix is a plain rectangular matrix with one
  89     non-zero element per row (with value 1). For instance, in the above figure, the matrix is :
  90     \f[
  91
  92     \begin{tabular}{|cccc|}
  93     1 & 0 & 0 & 0\\
  94     0 & 0 & 1 & 0\\
  95     1 & 0 & 0 & 0\\
  96     0 & 0 & 1 & 0\\
  97     \end{tabular}
  98     \f]
  99   */
 100
 101   fvm_nodal_t*  medmemMeshToFVMMesh(const MEDMEM::MESH* mesh)
 102   {
 103     // create an FVM structure from the paramesh structure
 104     std::string meshName(mesh->getName());//this line avoid that mesh->getName() object killed before fvm_nodal_create read the const char *.
 105     fvm_nodal_t * fvm_nodal = fvm_nodal_create(meshName.c_str(),mesh->getMeshDimension());
 106
 107     //loop on cell types
 108     int nbtypes = mesh->getNumberOfTypes(MED_EN::MED_CELL);
 109     const MED_EN::medGeometryElement* types = mesh->getTypes(MED_EN::MED_CELL);
 110     for (int itype=0; itype<nbtypes; itype++)
 111       {
 112         fvm_element_t fvm_type;
 113         switch (types[itype])
 114           {
 115           case MED_EN::MED_TRIA3 :
 116             fvm_type=FVM_FACE_TRIA;
 117             break;
 118           case MED_EN::MED_QUAD4 :
 119             fvm_type=FVM_FACE_QUAD;
 120             break;
 121           case MED_EN::MED_TETRA4 :
 122             fvm_type=FVM_CELL_TETRA;
 123             break;
 124           case MED_EN::MED_HEXA8 :
 125             fvm_type=FVM_CELL_HEXA;
 126             break;
 127           default:
 128             throw MEDEXCEPTION(" MED type  conversion to fvm is not handled yet.");
 129             break;
 130
 131           }
 132
 133         fvm_lnum_t nbelems = mesh->getNumberOfElements(MED_EN::MED_CELL, types[itype]);
 134         fvm_lnum_t* conn = new fvm_lnum_t[nbelems*(types[itype]%100)];
 135         const int* mesh_conn =mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, types[itype]);
 136         for (int i=0; i<nbelems*(types[itype]%100); i++)
 137           conn[i]=mesh_conn[i];
 138         //swapping trias
 139         if (types[itype]==MED_EN::MED_TRIA3)
 140           {
 141             for (int i=0; i<nbelems;i++)
 142               {
 143                 int tmp=conn[3*i];
 144                 conn[3*i]=mesh_conn[3*i+1];
 145                 conn[3*i+1]=tmp;
 146               }
 147           }
 148         //swapping tetras
 149         if (types[itype]==MED_EN::MED_TETRA4)
 150           {
 151             for (int i=0; i<nbelems;i++)
 152               {
 153                 int tmp=conn[4*i];
 154                 conn[4*i]=mesh_conn[4*i+1];
 155                 conn[4*i+1]=tmp;
 156               }
 157           }
 158         fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,0);
 159
 160         int nbnodes= mesh->getNumberOfNodes();
 161         int spacedim=mesh->getSpaceDimension();
 162         fvm_coord_t* coords = new fvm_coord_t[nbnodes*spacedim];
 163         const double* mesh_coords=mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
 164         for (int i=0; i<nbnodes*spacedim; i++)
 165           coords[i]=mesh_coords[i];
 166         fvm_nodal_transfer_vertices(fvm_nodal,coords);
 167       }
 168     return fvm_nodal;
 169   }
 170
 171   fvm_nodal_t*  medmemSupportToFVMMesh(const MEDMEM::SUPPORT* support)
 172   {
 173
 174     // create an FVM structure from the paramesh structure
 175     std::string supportName(support->getName());//this line avoid that support->getName() object killed before fvm_nodal_create read the const char *.
 176     fvm_nodal_t * fvm_nodal = fvm_nodal_create(supportName.c_str(),1);
 177
 178     const MEDMEM::MESH* mesh= support->getMesh();
 179
 180     //loop on cell types
 181     MED_EN::medEntityMesh entity = support->getEntity();
 182
 183     int nbtypes = support->getNumberOfTypes();
 184     const MED_EN::medGeometryElement* types = support->getTypes();
 185     int ioffset=0;
 186     const int* type_offset = support->getNumberIndex();
 187
 188     //browsing through all types
 189     for (int itype=0; itype<nbtypes; itype++)
 190       {
 191         fvm_element_t fvm_type;
 192         switch (types[itype])
 193           {
 194           case MED_EN::MED_TRIA3 :
 195             fvm_type=FVM_FACE_TRIA;
 196             break;
 197           case MED_EN::MED_QUAD4 :
 198             fvm_type=FVM_FACE_QUAD;
 199             break;
 200           case MED_EN::MED_TETRA4 :
 201             fvm_type=FVM_CELL_TETRA;
 202             break;
 203           case MED_EN::MED_HEXA8 :
 204             fvm_type=FVM_CELL_HEXA;
 205             break;
 206           default:
 207             throw MEDEXCEPTION(" MED type  conversion to fvm is not handled yet.");
 208             break;
 209
 210           }
 211         fvm_lnum_t nbelems = support->getNumberOfElements(types[itype]);
 212
 213         //for a partial support, defining the element numbers that are taken into
 214         //account in the support
 215         fvm_lnum_t* elem_numbers=0;
 216         if (!support->isOnAllElements())
 217           {
 218             elem_numbers = const_cast<fvm_lnum_t*> (support->getNumber(types[itype]));
 219
 220             //creating work arrays to store list of elems for partial supports
 221             if (itype>0)
 222               {
 223                 fvm_lnum_t* temp = new int[nbelems];
 224                 for (int i=0; i< nbelems; i++)
 225                   temp[i] = elem_numbers [i]-ioffset;
 226                 ioffset+=type_offset[itype];
 227                 elem_numbers = temp;
 228               }
 229           }
 230         //retrieving original mesh connectivity
 231         fvm_lnum_t* conn = const_cast<fvm_lnum_t*> (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,entity, types[itype]));
 232
 233         // adding the elements to the FVM structure
 234         fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,elem_numbers);
 235
 236         //cleaning work arrays (for partial supports)
 237         if (!support->isOnAllElements() && itype>0)
 238           delete[] elem_numbers;
 239
 240       }
 241     return fvm_nodal;
 242   }
 243
 244   NonCoincidentDEC::NonCoincidentDEC()
 245   {
 246   }
 247
 248   /*! Constructor of a non coincident \ref para-dec "DEC" with
 249    * a source group on which lies a field lying on a mesh and a
 250    * target group on which lies a mesh.
 251    *
 252    * \param source_group ProcessorGroup on the source side
 253    * \param target_group ProcessorGroup on the target side
 254    */
 255
 256   NonCoincidentDEC::NonCoincidentDEC(ProcessorGroup& source_group,
 257                                      ProcessorGroup& target_group)
 258     :DEC(source_group, target_group)
 259   {}
 260
 261   NonCoincidentDEC::~NonCoincidentDEC()
 262   {
 263   }
 264
 265   /*! Synchronization process. Calling this method
 266    * synchronizes the topologies so that the target side
 267    * gets the information which enable it to fetch the field value
 268    * from the source side.
 269    * A typical call is :
 270    \verbatim
 271    NonCoincidentDEC dec(source_group,target_group);
 272    dec.attachLocalField(field);
 273    dec.synchronize();
 274    \endverbatim
 275   */
 276   void NonCoincidentDEC::synchronize()
 277   {
 278
 279     //initializing FVM parallel environment
 280     const MPI_Comm* comm=dynamic_cast<const MPIProcessorGroup*> (_union_group)->getComm();
 281     fvm_parall_set_mpi_comm(*const_cast<MPI_Comm*> (comm));
 282
 283
 284     //setting up the communication DEC on both sides
 285
 286     if (_source_group->containsMyRank())
 287       {
 288         MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
 289         fvm_nodal_t* source_nodal = MEDCoupling::medmemMeshToFVMMesh(mesh);
 290
 291         int target_size = _target_group->size()  ;
 292         int start_rank=  _source_group->size();
 293         const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
 294
 295         _locator =  fvm_locator_create(1e-6,
 296                                        *comm,
 297                                        target_size,
 298                                        start_rank);
 299
 300         fvm_locator_set_nodal(_locator,
 301                               source_nodal,
 302                               mesh->getSpaceDimension(),
 303                               0,
 304                               NULL,
 305                               0);
 306
 307
 308         _nb_distant_points = fvm_locator_get_n_dist_points(_locator);
 309         _distant_coords = fvm_locator_get_dist_coords(_locator);
 310         _distant_locations = fvm_locator_get_dist_locations(_locator);
 311
 312       }
 313     if (_target_group->containsMyRank())
 314       {
 315         MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
 316
 317         fvm_nodal_t* target_nodal = MEDCoupling::medmemMeshToFVMMesh(mesh);
 318         int source_size = _source_group->size();
 319         int start_rank=  0 ;
 320         const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
 321
 322         _locator = fvm_locator_create(1e-6,
 323                                       *comm,
 324                                       source_size,
 325                                       start_rank);
 326         int nbcells = mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
 327         const MEDMEM::SUPPORT* support=_local_field->getField()->getSupport();
 328         MEDMEM::FIELD<double>* barycenter_coords = mesh->getBarycenter(support);
 329         const double* coords = barycenter_coords->getValue();
 330         fvm_locator_set_nodal(_locator,
 331                               target_nodal,
 332                               mesh->getSpaceDimension(),
 333                               nbcells,
 334                               NULL,
 335                               coords);
 336         delete barycenter_coords;
 337       }
 338   }
 339
 340
 341   /*! This method is called on the target group in order to
 342    * trigger the retrieveal of field data. It must
 343    * be called synchronously with a sendData() call on
 344    * the source group.
 345    */
 346   void NonCoincidentDEC::recvData()
 347   {
 348     int nbelems = _local_field->getField()->getSupport()->getMesh()->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
 349     int nbcomp =  _local_field->getField()->getNumberOfComponents();
 350     double* values = new double [nbelems*nbcomp];
 351     fvm_locator_exchange_point_var(_locator,
 352                                    0,
 353                                    values,
 354                                    0,
 355                                    sizeof(double),
 356                                    nbcomp,
 357                                    0);
 358     _local_field->getField()->setValue(values);
 359     if (_forced_renormalization_flag)
 360       renormalizeTargetField();
 361     delete[]values;
 362   }
 363
 364   /*! This method is called on the source group in order to
 365    * send field data. It must be called synchronously with
 366    * a recvData() call on
 367    * the target group.
 368    */
 369   void NonCoincidentDEC::sendData()
 370   {
 371     const double* values=_local_field->getField()->getValue();
 372     int nbcomp = _local_field->getField()->getNumberOfComponents();
 373     double* distant_values = new double [_nb_distant_points*nbcomp];
 374
 375     //cheap interpolation :  the value of the cell is transferred to the point
 376     for (int i=0; i<_nb_distant_points; i++)
 377       for (int j=0; j <nbcomp; j++)
 378         distant_values[i*nbcomp+j]=values[(_distant_locations[i]-1)*nbcomp+j];
 379
 380     fvm_locator_exchange_point_var(_locator,
 381                                    distant_values,
 382                                    0,
 383                                    0,
 384                                    sizeof(double),
 385                                    nbcomp,
 386                                    0);
 387
 388     delete [] distant_values;
 389     if (_forced_renormalization_flag)
 390       renormalizeTargetField();
 391
 392   }
 393 }