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     \anchor NonCoincidentDEC-det
  42     \class NonCoincidentDEC
  43
  44     \c NonCoincidentDEC enables nonconservative 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. The computation enables fast parallel localization, and is based on a point in element search, followed
  48     by a field evaluation at the point location. Thus, it is typically
  49     faster than the \ref InterpKernelDEC-det "InterpKernelDEC" which gives a
  50     \ref InterpKerRemapGlobal "conservative remapping".
  51     It is particularly true for the initialisation phase (synchronize)
  52     which is very computationnaly intensive in \ref InterpKernelDEC-det.
  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   /*! Constructor of a non coincident \ref para-dec "DEC" with
 246    * a source group on which lies a field lying on a mesh and a
 247    * target group on which lies a mesh.
 248    *
 249    * \param source_group ProcessorGroup on the source side
 250    * \param target_group ProcessorGroup on the target side
 251    */
 252
 253   NonCoincidentDEC::NonCoincidentDEC(ProcessorGroup& source_group,
 254                                      ProcessorGroup& target_group)
 255     :DEC(source_group, target_group)
 256   {}
 257
 258   NonCoincidentDEC::~NonCoincidentDEC()
 259   {
 260   }
 261
 262   /*! Synchronization process. Calling this method
 263    * synchronizes the topologies so that the target side
 264    * gets the information which enable it to fetch the field value
 265    * from the source side.
 266    * A typical call is :
 267    \verbatim
 268    NonCoincidentDEC dec(source_group,target_group);
 269    dec.attachLocalField(field);
 270    dec.synchronize();
 271    \endverbatim
 272   */
 273   void NonCoincidentDEC::synchronize()
 274   {
 275
 276     //initializing FVM parallel environment
 277     const MPI_Comm* comm=dynamic_cast<const MPIProcessorGroup*> (_union_group)->getComm();
 278     fvm_parall_set_mpi_comm(*const_cast<MPI_Comm*> (comm));
 279
 280
 281     //setting up the communication DEC on both sides
 282
 283     if (_source_group->containsMyRank())
 284       {
 285         MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
 286         fvm_nodal_t* source_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh);
 287
 288         int target_size = _target_group->size()  ;
 289         int start_rank=  _source_group->size();
 290         const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
 291
 292         _locator =  fvm_locator_create(1e-6,
 293                                        *comm,
 294                                        target_size,
 295                                        start_rank);
 296
 297         fvm_locator_set_nodal(_locator,
 298                               source_nodal,
 299                               mesh->getSpaceDimension(),
 300                               0,
 301                               NULL,
 302                               0);
 303
 304
 305         _nb_distant_points = fvm_locator_get_n_dist_points(_locator);
 306         _distant_coords = fvm_locator_get_dist_coords(_locator);
 307         _distant_locations = fvm_locator_get_dist_locations(_locator);
 308
 309       }
 310     if (_target_group->containsMyRank())
 311       {
 312         MEDMEM::MESH* mesh = _local_field->getField()->getSupport()->getMesh();
 313
 314         fvm_nodal_t* target_nodal = ParaMEDMEM::medmemMeshToFVMMesh(mesh);
 315         int source_size = _source_group->size();
 316         int start_rank=  0 ;
 317         const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*> (_union_group))->getComm();
 318
 319         _locator = fvm_locator_create(1e-6,
 320                                       *comm,
 321                                       source_size,
 322                                       start_rank);
 323         int nbcells = mesh->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
 324         const MEDMEM::SUPPORT* support=_local_field->getField()->getSupport();
 325         MEDMEM::FIELD<double>* barycenter_coords = mesh->getBarycenter(support);
 326         const double* coords = barycenter_coords->getValue();
 327         fvm_locator_set_nodal(_locator,
 328                               target_nodal,
 329                               mesh->getSpaceDimension(),
 330                               nbcells,
 331                               NULL,
 332                               coords);
 333         delete barycenter_coords;
 334       }
 335   }
 336
 337
 338   /*! This method is called on the target group in order to
 339    * trigger the retrieveal of field data. It must
 340    * be called synchronously with a sendData() call on
 341    * the source group.
 342    */
 343   void NonCoincidentDEC::recvData()
 344   {
 345     int nbelems = _local_field->getField()->getSupport()->getMesh()->getNumberOfElements(MED_EN::MED_CELL, MED_EN::MED_ALL_ELEMENTS);
 346     int nbcomp =  _local_field->getField()->getNumberOfComponents();
 347     double* values = new double [nbelems*nbcomp];
 348     fvm_locator_exchange_point_var(_locator,
 349                                    0,
 350                                    values,
 351                                    0,
 352                                    sizeof(double),
 353                                    nbcomp,
 354                                    0);
 355     _local_field->getField()->setValue(values);
 356     if (_forced_renormalization_flag)
 357       renormalizeTargetField();
 358     delete[]values;
 359   }
 360
 361   /*! This method is called on the source group in order to
 362    * send field data. It must be called synchronously with
 363    * a recvData() call on
 364    * the target group.
 365    */
 366   void NonCoincidentDEC::sendData()
 367   {
 368     const double* values=_local_field->getField()->getValue();
 369     int nbcomp = _local_field->getField()->getNumberOfComponents();
 370     double* distant_values = new double [_nb_distant_points*nbcomp];
 371
 372     //cheap interpolation :  the value of the cell is transfered to the point
 373     for (int i=0; i<_nb_distant_points; i++)
 374       for (int j=0; j <nbcomp; j++)
 375         distant_values[i*nbcomp+j]=values[(_distant_locations[i]-1)*nbcomp+j];
 376
 377     fvm_locator_exchange_point_var(_locator,
 378                                    distant_values,
 379                                    0,
 380                                    0,
 381                                    sizeof(double),
 382                                    nbcomp,
 383                                    0);
 384
 385     delete [] distant_values;
 386     if (_forced_renormalization_flag)
 387       renormalizeTargetField();
 388
 389   }
 390 }