src/ParaMEDMEM/InterpolationMatrix.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 "ParaMESH.hxx"
  20 #include "ProcessorGroup.hxx"
  21 #include "MxN_Mapping.hxx"
  22 #include "InterpolationMatrix.hxx"
  23 #include "TranslationRotationMatrix.hxx"
  24 #include "Interpolation.hxx"
  25 #include "Interpolation2D.txx"
  26 #include "Interpolation3DSurf.txx"
  27 #include "Interpolation3D.txx"
  28 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
  29 #include "InterpolationOptions.hxx"
  30 #include "VolSurfFormulae.hxx"
  31 #include "NormalizedUnstructuredMesh.hxx"
  32
  33 // class InterpolationMatrix
  34 // This class enables the storage of an interpolation matrix Wij mapping
  35 // source field Sj to target field Ti via Ti=Vi^(-1).Wij.Sj.
  36 // The matrix is built and stored on the processors belonging to the source
  37 // group.
  38
  39 using namespace std;
  40
  41 namespace ParaMEDMEM
  42 {
  43
  44   //   ====================================================================
  45   //   Creates an empty matrix structure linking two distributed supports.
  46   //   The method must be called by all processors belonging to source
  47   //   and target groups.
  48   //   param source_support local support
  49   //   param source_group processor group containing the local processors
  50   //   param target_group processor group containing the distant processors
  51   //   param method interpolation method
  52   //   ====================================================================
  53
  54   InterpolationMatrix::InterpolationMatrix(ParaMEDMEM::ParaMESH *source_support,
  55                                            const ProcessorGroup& source_group,
  56                                            const ProcessorGroup& target_group,
  57                                            const DECOptions& dec_options,
  58                                            const INTERP_KERNEL::InterpolationOptions& interp_options):
  59     _source_support(source_support->getCellMesh()),
  60     _mapping(source_group, target_group, dec_options),
  61     _source_group(source_group),
  62     _target_group(target_group),
  63     _source_volume(0),
  64     DECOptions(dec_options),
  65     INTERP_KERNEL::InterpolationOptions(interp_options)
  66   {
  67     int nbelems = _source_support->getNumberOfCells();
  68
  69     _row_offsets.resize(nbelems+1);
  70     for (int i=0; i<nbelems+1; i++)
  71       {
  72         _row_offsets[i]=0;
  73       }
  74
  75     _coeffs.resize(nbelems);
  76   }
  77
  78   InterpolationMatrix::~InterpolationMatrix()
  79   {
  80   }
  81
  82
  83   //   ======================================================================
  84   //   \brief Adds the contribution of a distant subdomain to the*
  85   //   interpolation matrix.
  86   //   The method adds contribution to the interpolation matrix.
  87   //   For each row of the matrix, elements are addded as
  88   //   a (column, coeff) pair in the _coeffs array. This column number refers
  89   //   to an element on the target side via the _col_offsets array.
  90   //   It is made of a series of (iproc, ielem) pairs.
  91   //   The number of elements per row is stored in the row_offsets array.
  92
  93   //   param distant_support local representation of the distant subdomain
  94   //   param iproc_distant id of the distant subdomain (in the distant group)
  95   //   param distant_elems mapping between the local representation of
  96   //   the subdomain and the actual elem ids on the distant subdomain
  97   //   ======================================================================
  98
  99   void InterpolationMatrix::addContribution ( MEDCouplingUMesh& distant_support,
 100                                               int iproc_distant,
 101                                               int* distant_elems,
 102                                               const std::string& srcMeth,
 103                                               const std::string& targetMeth)
 104   {
 105     if (distant_support.getMeshDimension() != _source_support->getMeshDimension())
 106       {
 107         throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
 108       }
 109     std::string interpMethod(srcMeth);
 110     interpMethod+=targetMeth;
 111     //creating the interpolator structure
 112     vector<map<int,double> > surfaces;
 113     int colSize=0;
 114     //computation of the intersection volumes between source and target elements
 115
 116     if ( distant_support.getMeshDimension() == 2
 117          && distant_support.getSpaceDimension() == 3 )
 118       {
 119         MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(&distant_support);
 120         MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(_source_support);
 121
 122         INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
 123         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
 124         target_wrapper.ReleaseTempArrays();
 125         source_wrapper.ReleaseTempArrays();
 126       }
 127     else if ( distant_support.getMeshDimension() == 2
 128               && distant_support.getSpaceDimension() == 2)
 129       {
 130         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(&distant_support);
 131         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(_source_support);
 132
 133         INTERP_KERNEL::Interpolation2D interpolator (*this);
 134         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
 135         target_wrapper.ReleaseTempArrays();
 136         source_wrapper.ReleaseTempArrays();
 137       }
 138     else if ( distant_support.getMeshDimension() == 3
 139               && distant_support.getSpaceDimension() ==3 )
 140       {
 141         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(&distant_support);
 142         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(_source_support);
 143
 144         INTERP_KERNEL::Interpolation3D interpolator (*this);
 145         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
 146         target_wrapper.ReleaseTempArrays();
 147         source_wrapper.ReleaseTempArrays();
 148       }
 149     else
 150       {
 151         throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
 152       }
 153
 154     int source_size=surfaces.size();
 155
 156     MEDCouplingFieldDouble *target_triangle_surf =
 157       getSupportVolumes(&distant_support);
 158     MEDCouplingFieldDouble *source_triangle_surf =
 159       getSupportVolumes(_source_support) ;
 160
 161     // Storing the source volumes
 162     _source_volume.resize(source_size);
 163     for (int i=0; i<source_size; i++)
 164       {
 165         _source_volume[i] = source_triangle_surf->getIJ(i,0);
 166       }
 167
 168     source_triangle_surf->decrRef();
 169
 170     //loop over the elements to build the interpolation
 171     //matrix structures
 172     for (int ielem=0; ielem < surfaces.size(); ielem++)
 173       {
 174         _row_offsets[ielem+1] += surfaces[ielem].size();
 175         //_source_indices.push_back(make_pair(iproc_distant, ielem));
 176
 177         for (map<int,double>::const_iterator iter = surfaces[ielem].begin();
 178              iter != surfaces[ielem].end();
 179              iter++)
 180           {
 181             //surface of the target triangle
 182             double surf = target_triangle_surf->getIJ(iter->first,0);
 183
 184             //locating the (iproc, itriangle) pair in the list of columns
 185             vector<pair<int,int> >::iterator iter2 =
 186               find(_col_offsets.begin(), _col_offsets.end(),
 187                    make_pair(iproc_distant,iter->first));
 188             int col_id;
 189
 190             if (iter2 == _col_offsets.end())
 191               {
 192                 //(iproc, itriangle) is not registered in the list
 193                 //of distant elements
 194
 195                 _col_offsets.push_back(make_pair(iproc_distant,iter->first));
 196                 col_id =_col_offsets.size();
 197                 _mapping.addElementFromSource(iproc_distant,
 198                                               distant_elems[iter->first]);
 199                 _target_volume.push_back(surf);
 200               }
 201             else
 202               {
 203                 col_id = iter2 - _col_offsets.begin() + 1;
 204               }
 205
 206             //the non zero coefficient is stored
 207             //ielem is the row,
 208             //col_id is the number of the column
 209             //iter->second is the value of the coefficient
 210
 211             _coeffs[ielem].push_back(make_pair(col_id,iter->second));
 212           }
 213       }
 214     target_triangle_surf->decrRef();
 215   }
 216
 217
 218   // ==================================================================
 219   // The call to this method updates the arrays on the target side
 220   //   so that they know which amount of data from which processor they
 221   //   should expect.
 222   //   That call makes actual interpolations via multiply method
 223   //   available.
 224   // ==================================================================
 225
 226   void InterpolationMatrix::prepare()
 227   {
 228     int nbelems = _source_support->getNumberOfCells();
 229     for (int ielem=0; ielem < nbelems; ielem++)
 230       {
 231         _row_offsets[ielem+1]+=_row_offsets[ielem];
 232       }
 233     _mapping.prepareSendRecv();
 234   }
 235
 236
 237   //   =======================================================================
 238   //   brief performs t=Ws, where t is the target field, s is the source field
 239
 240   //   The call to this method must be called both on the working side
 241   //   and on the idle side. On the working side, the vector  T=VT^(-1).(W.S)
 242   //   is computed and sent. On the idle side, no computation is done, but the
 243   //   result from the working side is received and the field is updated.
 244
 245   //   param field source field on processors involved on the source side,
 246   //   target field on processors on the target side
 247   //   =======================================================================
 248
 249   void InterpolationMatrix::multiply(MEDCouplingFieldDouble& field) const
 250   {
 251     int nbcomp = field.getArray()->getNumberOfComponents();
 252     vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
 253
 254     //computing the matrix multiply on source side
 255     if (_source_group.containsMyRank())
 256       {
 257         int nbrows = _source_support->getNumberOfCells();
 258
 259         // performing W.S
 260         // W is the intersection matrix
 261         // S is the source vector
 262
 263         for (int irow=0; irow<nbrows; irow++)
 264           {
 265             for (int icomp=0; icomp< nbcomp; icomp++)
 266               {
 267                 double coeff_row = field.getIJ(irow,icomp);
 268                 for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
 269                   {
 270                     int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
 271                     double value = _coeffs[irow][icol-_row_offsets[irow]].second;
 272                     target_value[(colid-1)*nbcomp+icomp]+=value*coeff_row;
 273                   }
 274               }
 275           }
 276
 277         // performing VT^(-1).(W.S)
 278         // where VT^(-1) is the inverse of the diagonal matrix containing
 279         // the volumes of target cells
 280
 281         for (int i=0; i<_col_offsets.size();i++)
 282           {
 283             for (int icomp=0; icomp<nbcomp; icomp++)
 284               {
 285                 target_value[i*nbcomp+icomp] /= _target_volume[i];
 286               }
 287           }
 288
 289       }
 290
 291     if (_target_group.containsMyRank())
 292       {
 293         int nbelems = field.getArray()->getNumberOfTuples() ;
 294         double* value = const_cast<double*> (field.getArray()->getPointer());
 295         for (int i=0; i<nbelems*nbcomp; i++)
 296           {
 297             value[i]=0.0;
 298           }
 299       }
 300
 301     //on source side : sending  T=VT^(-1).(W.S)
 302     //on target side :: receiving T and storing it in field
 303     _mapping.sendRecv(&target_value[0],field);
 304   }
 305
 306
 307   // =========================================================================
 308   // brief performs s=WTt, where t is the target field, s is the source field,
 309   // WT is the transpose matrix from W
 310
 311   //   The call to this method must be called both on the working side
 312   //   and on the idle side. On the working side, the target vector T is
 313   //   received and the vector  S=VS^(-1).(WT.T) is computed to update
 314   //   the field.
 315   //   On the idle side, no computation is done, but the field is sent.
 316
 317   //   param field source field on processors involved on the source side,
 318   //   target field on processors on the target side
 319   // =========================================================================
 320
 321   void InterpolationMatrix::transposeMultiply(MEDCouplingFieldDouble& field) const
 322   {
 323     //  int nbcomp = field.getNumberOfComponents();
 324     int nbcomp = field.getArray()->getNumberOfComponents();
 325     vector<double> source_value(_col_offsets.size()* nbcomp,0.0);
 326     _mapping.reverseSendRecv(&source_value[0],field);
 327
 328     //treatment of the transpose matrix multiply on the source side
 329     if (_source_group.containsMyRank())
 330       {
 331         int nbrows    = _source_support->getNumberOfCells();
 332         double *array = field.getArray()->getPointer() ;
 333
 334         // Initialization
 335         std::fill(array, array+nbrows*nbcomp, 0.0) ;
 336
 337         //performing WT.T
 338         //WT is W transpose
 339         //T is the target vector
 340         for (int irow = 0; irow < nbrows; irow++)
 341           {
 342             for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
 343               {
 344                 int colid    = _coeffs[irow][icol-_row_offsets[irow]].first;
 345                 double value = _coeffs[irow][icol-_row_offsets[irow]].second;
 346
 347                 for (int icomp=0; icomp<nbcomp; icomp++)
 348                   {
 349                     double coeff_row = source_value[(colid-1)*nbcomp+icomp];
 350                     array[irow*nbcomp+icomp] += value*coeff_row;
 351                   }
 352               }
 353           }
 354
 355         //performing VS^(-1).(WT.T)
 356         //VS^(-1) is the inverse of the diagonal matrix storing
 357         //volumes of the source cells
 358
 359         for (int irow=0; irow<nbrows; irow++)
 360           {
 361             for (int icomp=0; icomp<nbcomp; icomp++)
 362               {
 363                 array[irow*nbcomp+icomp] /= _source_volume[irow];
 364               }
 365           }
 366
 367       }
 368   }
 369
 370   MEDCouplingFieldDouble* InterpolationMatrix::getSupportVolumes(MEDCouplingMesh * mesh)
 371   {
 372     if(!mesh->isStructured())
 373       return getSupportUnstructuredVolumes((MEDCouplingUMesh *)mesh);
 374     else
 375       throw INTERP_KERNEL::Exception("Not implemented yet !!!");
 376   }
 377
 378   //   ====================================================================
 379   //   brief returns the volumes of the cells underlying the field \a field
 380
 381   //   For 2D geometries, the returned field contains the areas.
 382   //   For 3D geometries, the returned field contains the volumes.
 383
 384   //   param field field on which cells the volumes are required
 385   //   return field containing the volumes
 386   //   ====================================================================
 387
 388   MEDCouplingFieldDouble* InterpolationMatrix::getSupportUnstructuredVolumes(MEDCouplingUMesh * mesh)
 389   {
 390     int ipt, type ;
 391     int nbelem       = mesh->getNumberOfCells() ;
 392     int dim_mesh     = mesh->getMeshDimension();
 393     int dim_space    = mesh->getSpaceDimension() ;
 394     double *coords    = mesh->getCoords()->getPointer() ;
 395     int *connec       = mesh->getNodalConnectivity()->getPointer() ;
 396     int *connec_index = mesh->getNodalConnectivityIndex()->getPointer() ;
 397
 398
 399     MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS);
 400     DataArrayDouble* array = DataArrayDouble::New() ;
 401     array->alloc(nbelem, 1) ;
 402     double *area_vol = array->getPointer() ;
 403
 404     switch (dim_mesh)
 405       {
 406       case 2: // getting the areas
 407         for ( int iel=0 ; iel<nbelem ; iel++ )
 408           {
 409             ipt = connec_index[iel] ;
 410             type = connec[ipt] ;
 411
 412             switch ( type )
 413               {
 414               case INTERP_KERNEL::NORM_TRI3 :
 415               case INTERP_KERNEL::NORM_TRI6 :
 416                 {
 417                   int N1 = connec[ipt+1];
 418                   int N2 = connec[ipt+2];
 419                   int N3 = connec[ipt+3];
 420
 421                   area_vol[iel]=INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1),
 422                                                                     coords+(dim_space*N2),
 423                                                                     coords+(dim_space*N3),
 424                                                                     dim_space);
 425                 }
 426                 break ;
 427
 428               case INTERP_KERNEL::NORM_QUAD4 :
 429               case INTERP_KERNEL::NORM_QUAD8 :
 430                 {
 431                   int N1 = connec[ipt+1];
 432                   int N2 = connec[ipt+2];
 433                   int N3 = connec[ipt+3];
 434                   int N4 = connec[ipt+4];
 435
 436                   area_vol[iel]=INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1,
 437                                                                     coords+dim_space*N2,
 438                                                                     coords+dim_space*N3,
 439                                                                     coords+dim_space*N4,
 440                                                                     dim_space) ;
 441                 }
 442                 break ;
 443
 444               case INTERP_KERNEL::NORM_POLYGON :
 445                 {
 446                   // We must remember that the first item is the type. That's
 447                   // why we substract 1 to get the number of nodes of this polygon
 448                   int size = connec_index[iel+1] - connec_index[iel] - 1 ;
 449
 450                   double **pts = new double *[size] ;
 451
 452                   for ( int inod=0 ; inod<size ; inod++ )
 453                     {
 454                       // Remember the first item is the type
 455                       pts[inod] = coords+dim_space*connec[ipt+inod+1] ;
 456                     }
 457
 458                   area_vol[iel]=INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,
 459                                                                      size,dim_space);
 460                   delete [] pts;
 461                 }
 462                 break ;
 463
 464               default :
 465                 throw INTERP_KERNEL::Exception("Bad Support to get Areas on it !");
 466
 467               } // End of switch
 468
 469           } // End of the loop over the cells
 470         break;
 471       case 3: // getting the volumes
 472         for ( int iel=0 ; iel<nbelem ; iel++ )
 473           {
 474             ipt = connec_index[iel] ;
 475             type = connec[ipt] ;
 476
 477             switch ( type )
 478               {
 479               case INTERP_KERNEL::NORM_TETRA4 :
 480               case INTERP_KERNEL::NORM_TETRA10 :
 481                 {
 482                   int N1 = connec[ipt+1];
 483                   int N2 = connec[ipt+2];
 484                   int N3 = connec[ipt+3];
 485                   int N4 = connec[ipt+4];
 486
 487                   area_vol[iel]=INTERP_KERNEL::calculateVolumeForTetra(coords+dim_space*N1,
 488                                                                        coords+dim_space*N2,
 489                                                                        coords+dim_space*N3,
 490                                                                        coords+dim_space*N4) ;
 491                 }
 492                 break ;
 493
 494               case INTERP_KERNEL::NORM_PYRA5 :
 495               case INTERP_KERNEL::NORM_PYRA13 :
 496                 {
 497                   int N1 = connec[ipt+1];
 498                   int N2 = connec[ipt+2];
 499                   int N3 = connec[ipt+3];
 500                   int N4 = connec[ipt+4];
 501                   int N5 = connec[ipt+5];
 502
 503                   area_vol[iel]=INTERP_KERNEL::calculateVolumeForPyra(coords+dim_space*N1,
 504                                                                       coords+dim_space*N2,
 505                                                                       coords+dim_space*N3,
 506                                                                       coords+dim_space*N4,
 507                                                                       coords+dim_space*N5) ;
 508                 }
 509                 break ;
 510
 511               case INTERP_KERNEL::NORM_PENTA6 :
 512               case INTERP_KERNEL::NORM_PENTA15 :
 513                 {
 514                   int N1 = connec[ipt+1];
 515                   int N2 = connec[ipt+2];
 516                   int N3 = connec[ipt+3];
 517                   int N4 = connec[ipt+4];
 518                   int N5 = connec[ipt+5];
 519                   int N6 = connec[ipt+6];
 520
 521                   area_vol[iel]=INTERP_KERNEL::calculateVolumeForPenta(coords+dim_space*N1,
 522                                                                        coords+dim_space*N2,
 523                                                                        coords+dim_space*N3,
 524                                                                        coords+dim_space*N4,
 525                                                                        coords+dim_space*N5,
 526                                                                        coords+dim_space*N6) ;
 527                 }
 528                 break ;
 529
 530               case INTERP_KERNEL::NORM_HEXA8 :
 531               case INTERP_KERNEL::NORM_HEXA20 :
 532                 {
 533                   int N1 = connec[ipt+1];
 534                   int N2 = connec[ipt+2];
 535                   int N3 = connec[ipt+3];
 536                   int N4 = connec[ipt+4];
 537                   int N5 = connec[ipt+5];
 538                   int N6 = connec[ipt+6];
 539                   int N7 = connec[ipt+7];
 540                   int N8 = connec[ipt+8];
 541
 542                   area_vol[iel]=INTERP_KERNEL::calculateVolumeForHexa(coords+dim_space*N1,
 543                                                                       coords+dim_space*N2,
 544                                                                       coords+dim_space*N3,
 545                                                                       coords+dim_space*N4,
 546                                                                       coords+dim_space*N5,
 547                                                                       coords+dim_space*N6,
 548                                                                       coords+dim_space*N7,
 549                                                                       coords+dim_space*N8) ;
 550                 }
 551                 break ;
 552
 553               case INTERP_KERNEL::NORM_POLYHED :
 554                 {
 555                   throw INTERP_KERNEL::Exception("Not yet implemented !");
 556                 }
 557                 break ;
 558
 559               default:
 560                 throw INTERP_KERNEL::Exception("Bad Support to get Volume on it !");
 561               }
 562           }
 563         break;
 564       default:
 565         throw INTERP_KERNEL::Exception("interpolation is not available for this dimension");
 566       }
 567
 568     field->setArray(array) ;
 569     array->decrRef();
 570     field->setMesh(mesh) ;
 571
 572     return field ;
 573   }
 574 }