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