src/ParaMEDMEM/NonCoincidentDEC.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 "DEC.hxx"
  27 #include "NonCoincidentDEC.hxx"
  28
  29 extern "C" {
  30 #include <fvm_parall.h>
  31 #include <fvm_nodal.h>
  32 #include <fvm_nodal_append.h>
  33 #include <fvm_locator.h>
  34 }
  35
  36 namespace ParaMEDMEM
  37 {
  38
  39   /*!
  40     \defgroup noncoincidentdec NonCoincidentDEC
  41
  42     \section overview Overview
  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, based on the
  48     FVM library. The computation 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 intersectiondec which gives a \ref conservativeremapping.
  51     It is particularly true for the initialisation phase (synchronize)
  52     which is very computationnaly intensive in \ref intersectiondec.
  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     fvm_nodal_t * fvm_nodal = fvm_nodal_create(mesh->getName().c_str(),mesh->getMeshDimension());
 102
 103     //loop on cell types
 104     int nbtypes = mesh->getNumberOfTypes(MED_EN::MED_CELL);
 105     const MED_EN::medGeometryElement* types = mesh->getTypes(MED_EN::MED_CELL);
 106     for (int itype=0; itype<nbtypes; itype++)
 107       {
 108         fvm_element_t fvm_type;
 109         switch (types[itype])
 110           {
 111           case MED_EN::MED_TRIA3 :
 112             fvm_type=FVM_FACE_TRIA;
 113             break;
 114           case MED_EN::MED_QUAD4 :
 115             fvm_type=FVM_FACE_QUAD;
 116             break;
 117           case MED_EN::MED_TETRA4 :
 118             fvm_type=FVM_CELL_TETRA;
 119             break;
 120           case MED_EN::MED_HEXA8 :
 121             fvm_type=FVM_CELL_HEXA;
 122             break;
 123           default:
 124             throw MEDEXCEPTION(" MED type  conversion to fvm is not handled yet.");
 125             break;
 126
 127           }
 128
 129         fvm_lnum_t nbelems = mesh->getNumberOfElements(MED_EN::MED_CELL, types[itype]);
 130         fvm_lnum_t* conn = new fvm_lnum_t[nbelems*(types[itype]%100)];
 131         const int* mesh_conn =mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL, MED_EN::MED_CELL, types[itype]);
 132         for (int i=0; i<nbelems*(types[itype]%100); i++)
 133           conn[i]=mesh_conn[i];
 134         //swapping trias
 135         if (types[itype]==MED_EN::MED_TRIA3)
 136           {
 137             for (int i=0; i<nbelems;i++)
 138               {
 139                 int tmp=conn[3*i];
 140                 conn[3*i]=mesh_conn[3*i+1];
 141                 conn[3*i+1]=tmp;
 142               }
 143           }
 144         //swapping tetras
 145         if (types[itype]==MED_EN::MED_TETRA4)
 146           {
 147             for (int i=0; i<nbelems;i++)
 148               {
 149                 int tmp=conn[4*i];
 150                 conn[4*i]=mesh_conn[4*i+1];
 151                 conn[4*i+1]=tmp;
 152               }
 153           }
 154         fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,0);
 155
 156         int nbnodes= mesh->getNumberOfNodes();
 157         int spacedim=mesh->getSpaceDimension();
 158         fvm_coord_t* coords = new fvm_coord_t[nbnodes*spacedim];
 159         const double* mesh_coords=mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
 160         for (int i=0; i<nbnodes*spacedim; i++)
 161           coords[i]=mesh_coords[i];
 162         fvm_nodal_transfer_vertices(fvm_nodal,coords);
 163       }
 164     return fvm_nodal;
 165   }
 166
 167   fvm_nodal_t*  medmemSupportToFVMMesh(const MEDMEM::SUPPORT* support)
 168   {
 169
 170     // create an FVM structure from the paramesh structure
 171     fvm_nodal_t * fvm_nodal = fvm_nodal_create(support->getName().c_str(),1);
 172
 173     const MEDMEM::MESH* mesh= support->getMesh();
 174
 175     //loop on cell types
 176     MED_EN::medEntityMesh entity = support->getEntity();
 177
 178     int nbtypes = support->getNumberOfTypes();
 179     const MED_EN::medGeometryElement* types = support->getTypes();
 180     int ioffset=0;
 181     const int* type_offset = support->getNumberIndex();
 182
 183     //browsing through all types
 184     for (int itype=0; itype<nbtypes; itype++)
 185       {
 186         fvm_element_t fvm_type;
 187         switch (types[itype])
 188           {
 189           case MED_EN::MED_TRIA3 :
 190             fvm_type=FVM_FACE_TRIA;
 191             break;
 192           case MED_EN::MED_QUAD4 :
 193             fvm_type=FVM_FACE_QUAD;
 194             break;
 195           case MED_EN::MED_TETRA4 :
 196             fvm_type=FVM_CELL_TETRA;
 197             break;
 198           case MED_EN::MED_HEXA8 :
 199             fvm_type=FVM_CELL_HEXA;
 200             break;
 201           default:
 202             throw MEDEXCEPTION(" MED type  conversion to fvm is not handled yet.");
 203             break;
 204
 205           }
 206         fvm_lnum_t nbelems = support->getNumberOfElements(types[itype]);
 207
 208         //for a partial support, defining the element numbers that are taken into
 209         //account in the support
 210         fvm_lnum_t* elem_numbers=0;
 211         if (!support->isOnAllElements())
 212           {
 213             elem_numbers = const_cast<fvm_lnum_t*> (support->getNumber(types[itype]));
 214
 215             //creating work arrays to store list of elems for partial suports
 216             if (itype>0)
 217               {
 218                 fvm_lnum_t* temp = new int[nbelems];
 219                 for (int i=0; i< nbelems; i++)
 220                   temp[i] = elem_numbers [i]-ioffset;
 221                 ioffset+=type_offset[itype];
 222                 elem_numbers = temp;
 223               }
 224           }
 225         //retrieving original mesh connectivity
 226         fvm_lnum_t* conn = const_cast<fvm_lnum_t*> (mesh->getConnectivity(MED_EN::MED_FULL_INTERLACE,MED_EN::MED_NODAL,entity, types[itype]));
 227
 228         // adding the elements to the FVM structure
 229         fvm_nodal_append_by_transfer(fvm_nodal, nbelems, fvm_type,0,0,0,conn,elem_numbers);
 230
 231         //cleaning work arrays (for partial supports)
 232         if (!support->isOnAllElements() && itype>0)
 233           delete[] elem_numbers;
 234
 235       }
 236     return fvm_nodal;
 237   }
 238
 239   NonCoincidentDEC::NonCoincidentDEC()
 240   {
 241   }
 242
 243   /*!
 244     \addtogroup noncoincidentdec
 245     @{
 246   */
 247
 248   /*! Constructor of a non coincident 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 = ParaMEDMEM::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 = ParaMEDMEM::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 transfered 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   /*!
 394     @}
 395   */
 396 }