src/ParaMEDMEM/InterpolationMatrix.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   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 "ParaMESH.hxx"
  21 #include "ParaFIELD.hxx"
  22 #include "ProcessorGroup.hxx"
  23 #include "MxN_Mapping.hxx"
  24 #include "InterpolationMatrix.hxx"
  25 #include "TranslationRotationMatrix.hxx"
  26 #include "Interpolation.hxx"
  27 #include "Interpolation1D.txx"
  28 #include "Interpolation2DCurve.hxx"
  29 #include "Interpolation2D.txx"
  30 #include "Interpolation3DSurf.hxx"
  31 #include "Interpolation3D.txx"
  32 #include "Interpolation2D3D.txx"
  33 #include "Interpolation2D1D.txx"
  34 #include "MEDCouplingUMesh.hxx"
  35 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
  36 #include "InterpolationOptions.hxx"
  37 #include "NormalizedUnstructuredMesh.hxx"
  38 #include "ElementLocator.hxx"
  39
  40 #include <algorithm>
  41
  42 using namespace std;
  43
  44 namespace MEDCoupling
  45 {
  46
  47   /**!
  48      Creates an empty matrix structure linking two distributed supports.
  49      The method must be called by all processors belonging to source
  50      and target groups.
  51      \param source_field local source field
  52      \param source_group processor group containing the local processors
  53      \param target_group processor group containing the distant processors
  54      \param dec_options DEC options (including projection method P0, P1)
  55      \param interp_options interpolation options
  56   */
  57   InterpolationMatrix::InterpolationMatrix(const MEDCoupling::ParaFIELD *source_field,
  58                                            const ProcessorGroup& source_group,
  59                                            const ProcessorGroup& target_group,
  60                                            const DECOptions& dec_options,
  61                                            const INTERP_KERNEL::InterpolationOptions& interp_options):
  62     INTERP_KERNEL::InterpolationOptions(interp_options),
  63     DECOptions(dec_options),
  64     _source_field(source_field),
  65     _source_support(source_field->getSupport()->getCellMesh()),
  66     _mapping(source_group, target_group, dec_options),
  67     _source_group(source_group),
  68     _target_group(target_group)
  69   {
  70     mcIdType nbelems = source_field->getField()->getNumberOfTuples();
  71     _row_offsets.resize(nbelems+1);
  72     _coeffs.resize(nbelems);
  73     _target_volume.resize(nbelems);
  74   }
  75
  76   InterpolationMatrix::~InterpolationMatrix()
  77   {
  78   }
  79
  80
  81   /*!
  82      \brief Adds the contribution of a distant subdomain to the*
  83      interpolation matrix.
  84      The method adds contribution to the interpolation matrix.
  85      For each row of the matrix, elements are added as
  86      a (column, coeff) pair in the _coeffs array. This column number refers
  87      to an element on the target side via the _col_offsets array.
  88      It is made of a series of (iproc, ielem) pairs.
  89      The number of elements per row is stored in the row_offsets array.
  90
  91      param distant_support local representation of the distant subdomain
  92      param iproc_distant id of the distant subdomain (in the distant group)
  93      param distant_elems mapping between the local representation of
  94      the subdomain and the actual elem ids on the distant subdomain
  95    */
  96   void InterpolationMatrix::addContribution ( MEDCouplingPointSet& distant_support,
  97                                               int iproc_distant,
  98                                               const mcIdType* distant_elems,
  99                                               const std::string& srcMeth,
 100                                               const std::string& targetMeth)
 101   {
 102     std::string interpMethod(targetMeth);
 103     interpMethod+=srcMeth;
 104     //creating the interpolator structure
 105     vector<map<mcIdType,double> > surfaces;
 106     //computation of the intersection volumes between source and target elements
 107     MEDCouplingUMesh *distant_supportC=dynamic_cast<MEDCouplingUMesh *>(&distant_support);
 108     MEDCouplingUMesh *source_supportC=dynamic_cast<MEDCouplingUMesh *>(_source_support);
 109     if ( distant_support.getMeshDimension() == -1 )
 110       {
 111         if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==2)
 112           {
 113             MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(source_supportC);
 114             INTERP_KERNEL::Interpolation2D interpolation(*this);
 115             interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
 116           }
 117         else if(source_supportC->getMeshDimension()==3 && source_supportC->getSpaceDimension()==3)
 118           {
 119             MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(source_supportC);
 120             INTERP_KERNEL::Interpolation3D interpolation(*this);
 121             interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
 122           }
 123         else if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==3)
 124           {
 125             MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(source_supportC);
 126             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 127             interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
 128           }
 129         else
 130           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 131       }
 132     else if ( source_supportC->getMeshDimension() == -1 )
 133       {
 134         if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==2)
 135           {
 136             MEDCouplingNormalizedUnstructuredMesh<2,2> distant_mesh_wrapper(distant_supportC);
 137             INTERP_KERNEL::Interpolation2D interpolation(*this);
 138             interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
 139           }
 140         else if(distant_supportC->getMeshDimension()==3 && distant_supportC->getSpaceDimension()==3)
 141           {
 142             MEDCouplingNormalizedUnstructuredMesh<3,3> distant_mesh_wrapper(distant_supportC);
 143             INTERP_KERNEL::Interpolation3D interpolation(*this);
 144             interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
 145           }
 146         else if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==3)
 147           {
 148             MEDCouplingNormalizedUnstructuredMesh<3,2> distant_mesh_wrapper(distant_supportC);
 149             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 150             interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
 151           }
 152         else
 153           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
 154       }
 155     else if ( distant_support.getMeshDimension() == 2
 156               && _source_support->getMeshDimension() == 3
 157               && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
 158       {
 159         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
 160         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
 161         INTERP_KERNEL::Interpolation2D3D interpolator (*this);
 162         interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 163         target_wrapper.releaseTempArrays();
 164         source_wrapper.releaseTempArrays();
 165       }
 166     else if ( distant_support.getMeshDimension() == 1
 167               && _source_support->getMeshDimension() == 2
 168               && distant_support.getSpaceDimension() == 2 && _source_support->getSpaceDimension() == 2)
 169       {
 170         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
 171         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
 172         INTERP_KERNEL::Interpolation2D1D interpolator (*this);
 173         interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 174         target_wrapper.releaseTempArrays();
 175         source_wrapper.releaseTempArrays();
 176       }
 177     else if ( distant_support.getMeshDimension() == 3
 178               && _source_support->getMeshDimension() == 1
 179               && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
 180       {
 181         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
 182         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
 183         INTERP_KERNEL::Interpolation3D interpolator (*this);
 184         interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 185         target_wrapper.releaseTempArrays();
 186         source_wrapper.releaseTempArrays();
 187       }
 188     else if (distant_support.getMeshDimension() != _source_support->getMeshDimension())
 189       {
 190         throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
 191       }
 192     else if( distant_support.getMeshDimension() == 1
 193              && distant_support.getSpaceDimension() == 1 )
 194       {
 195         MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(distant_supportC);
 196         MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(source_supportC);
 197
 198         INTERP_KERNEL::Interpolation1D interpolation(*this);
 199         interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 200         target_wrapper.releaseTempArrays();
 201         source_wrapper.releaseTempArrays();
 202       }
 203     else if( distant_support.getMeshDimension() == 1
 204              && distant_support.getSpaceDimension() == 2 )
 205       {
 206         MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(distant_supportC);
 207         MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(source_supportC);
 208
 209         INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
 210         interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 211         target_wrapper.releaseTempArrays();
 212         source_wrapper.releaseTempArrays();
 213       }
 214     else if ( distant_support.getMeshDimension() == 2
 215               && distant_support.getSpaceDimension() == 3 )
 216       {
 217         MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC);
 218         MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC);
 219
 220         INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
 221         interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 222         target_wrapper.releaseTempArrays();
 223         source_wrapper.releaseTempArrays();
 224       }
 225     else if ( distant_support.getMeshDimension() == 2
 226               && distant_support.getSpaceDimension() == 2)
 227       {
 228         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
 229         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
 230
 231         INTERP_KERNEL::Interpolation2D interpolator (*this);
 232         interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 233         target_wrapper.releaseTempArrays();
 234         source_wrapper.releaseTempArrays();
 235       }
 236     else if ( distant_support.getMeshDimension() == 3
 237               && distant_support.getSpaceDimension() == 3 )
 238       {
 239         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
 240         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
 241
 242         INTERP_KERNEL::Interpolation3D interpolator (*this);
 243         interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
 244         target_wrapper.releaseTempArrays();
 245         source_wrapper.releaseTempArrays();
 246       }
 247     else
 248       {
 249         throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
 250       }
 251     bool needTargetSurf=isSurfaceComputationNeeded(targetMeth);
 252
 253     MEDCouplingFieldDouble *target_triangle_surf=0;
 254     if(needTargetSurf)
 255       target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus());
 256     fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf);
 257
 258     if(needTargetSurf)
 259       target_triangle_surf->decrRef();
 260   }
 261
 262   void InterpolationMatrix::fillDSFromVM(int iproc_distant, const mcIdType* distant_elems, const std::vector< std::map<mcIdType,double> >& values, MEDCouplingFieldDouble *surf)
 263   {
 264     //loop over the elements to build the interpolation
 265     //matrix structures
 266     std::size_t source_size=values.size();
 267     for (std::size_t ielem=0; ielem < source_size; ielem++)
 268       {
 269         _row_offsets[ielem+1] += ToIdType(values[ielem].size());
 270         for(map<mcIdType,double>::const_iterator iter=values[ielem].begin();iter!=values[ielem].end();iter++)
 271           {
 272             mcIdType localId;
 273             if(distant_elems)
 274               localId=distant_elems[iter->first];
 275             else
 276               localId=iter->first;
 277             //locating the (iproc, itriangle) pair in the list of columns
 278             map<pair<int,mcIdType>,mcIdType >::iterator iter2 = _col_offsets.find(make_pair(iproc_distant,localId));
 279             mcIdType col_id;
 280
 281             if (iter2 == _col_offsets.end())
 282               {
 283                 //(iproc, itriangle) is not registered in the list
 284                 //of distant elements
 285                 col_id =ToIdType(_col_offsets.size());
 286                 _col_offsets.insert(make_pair(make_pair(iproc_distant,localId),col_id));
 287                 _mapping.addElementFromSource(iproc_distant,localId);
 288               }
 289             else
 290               {
 291                 col_id = iter2->second;
 292               }
 293             //the non zero coefficient is stored
 294             //ielem is the row,
 295             //col_id is the number of the column
 296             //iter->second is the value of the coefficient
 297             if(surf)
 298               {
 299                 double surface = surf->getIJ(iter->first,0);
 300                 _target_volume[ielem].push_back(surface);
 301               }
 302             _coeffs[ielem].push_back(make_pair(col_id,iter->second));
 303           }
 304       }
 305   }
 306
 307   void InterpolationMatrix::serializeMe(std::vector< std::vector< std::map<mcIdType,double> > >& data1, std::vector<int>& data2) const
 308   {
 309     data1.clear();
 310     data2.clear();
 311     const std::vector<std::pair<int,mcIdType> >& sendingIds=_mapping.getSendingIds();
 312     std::set<int> procsS;
 313     for(std::vector<std::pair<int,mcIdType> >::const_iterator iter1=sendingIds.begin();iter1!=sendingIds.end();iter1++)
 314       procsS.insert((*iter1).first);
 315     data1.resize(procsS.size());
 316     data2.resize(procsS.size());
 317     std::copy(procsS.begin(),procsS.end(),data2.begin());
 318     std::map<int,int> fastProcAcc;
 319     int id=0;
 320     for(std::set<int>::const_iterator iter2=procsS.begin();iter2!=procsS.end();iter2++,id++)
 321       fastProcAcc[*iter2]=id;
 322     mcIdType nbOfSrcElt=ToIdType(_coeffs.size());
 323     for(std::vector< std::vector< std::map<mcIdType,double> > >::iterator iter3=data1.begin();iter3!=data1.end();iter3++)
 324       (*iter3).resize(nbOfSrcElt);
 325     id=0;
 326     for(std::vector< std::vector< std::pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,id++)
 327       {
 328         for(std::vector< std::pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++)
 329           {
 330             const std::pair<int,mcIdType>& elt=sendingIds[(*iter5).first];
 331             data1[fastProcAcc[elt.first]][id][elt.second]=(*iter5).second;
 332           }
 333       }
 334   }
 335
 336   void InterpolationMatrix::initialize()
 337   {
 338     mcIdType lgth=ToIdType(_coeffs.size());
 339     _row_offsets.clear(); _row_offsets.resize(lgth+1);
 340     _coeffs.clear(); _coeffs.resize(lgth);
 341     _target_volume.clear(); _target_volume.resize(lgth);
 342     _col_offsets.clear();
 343     _mapping.initialize();
 344   }
 345
 346   void InterpolationMatrix::finishContributionW(ElementLocator& elementLocator)
 347   {
 348     NatureOfField nature=elementLocator.getLocalNature();
 349     switch(nature)
 350       {
 351       case IntensiveMaximum:
 352         computeConservVolDenoW(elementLocator);
 353         break;
 354       case ExtensiveMaximum:
 355         {
 356           if(!elementLocator.isM1DCorr())
 357             computeIntegralDenoW(elementLocator);
 358           else
 359             computeGlobConstraintDenoW(elementLocator);
 360           break;
 361         }
 362       case ExtensiveConservation:
 363         computeGlobConstraintDenoW(elementLocator);
 364         break;
 365       case IntensiveConservation:
 366         {
 367           if(!elementLocator.isM1DCorr())
 368             computeRevIntegralDenoW(elementLocator);
 369           else
 370             computeConservVolDenoW(elementLocator);
 371           break;
 372         }
 373       default:
 374         throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
 375         break;
 376       }
 377   }
 378
 379   void InterpolationMatrix::finishContributionL(ElementLocator& elementLocator)
 380   {
 381     NatureOfField nature=elementLocator.getLocalNature();
 382     switch(nature)
 383       {
 384       case IntensiveMaximum:
 385         computeConservVolDenoL(elementLocator);
 386         break;
 387       case ExtensiveMaximum:
 388         {
 389           if(!elementLocator.isM1DCorr())
 390             computeIntegralDenoL(elementLocator);
 391           else
 392             computeConservVolDenoL(elementLocator);
 393           break;
 394         }
 395       case ExtensiveConservation:
 396         //this is not a bug doing like IntensiveMaximum
 397         computeConservVolDenoL(elementLocator);
 398         break;
 399       case IntensiveConservation:
 400         {
 401           if(!elementLocator.isM1DCorr())
 402             computeRevIntegralDenoL(elementLocator);
 403           else
 404             computeConservVolDenoL(elementLocator);
 405           break;
 406         }
 407       default:
 408         throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
 409         break;
 410       }
 411   }
 412
 413   void InterpolationMatrix::computeConservVolDenoW(ElementLocator& elementLocator)
 414   {
 415     computeGlobalColSum(_deno_reverse_multiply);
 416     computeGlobalRowSum(elementLocator,_deno_multiply,_deno_reverse_multiply);
 417   }
 418
 419   void InterpolationMatrix::computeConservVolDenoL(ElementLocator& elementLocator)
 420   {
 421     int pol1=elementLocator.sendPolicyToWorkingSideL();
 422     if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
 423       {
 424         elementLocator.recvFromWorkingSideL();
 425         elementLocator.sendToWorkingSideL();
 426       }
 427     else if(ElementLocator::CUMULATIVE_POLICY)
 428       {
 429         //ask for lazy side to deduce ids eventually missing on working side and to send it back.
 430         elementLocator.recvLocalIdsFromWorkingSideL();
 431         elementLocator.sendCandidatesGlobalIdsToWorkingSideL();
 432         elementLocator.recvCandidatesForAddElementsL();
 433         elementLocator.sendAddElementsToWorkingSideL();
 434         //Working side has updated its eventually missing ids updates its global ids with lazy side procs contribution
 435         elementLocator.recvLocalIdsFromWorkingSideL();
 436         elementLocator.sendGlobalIdsToWorkingSideL();
 437         //like no post treatment
 438         elementLocator.recvFromWorkingSideL();
 439         elementLocator.sendToWorkingSideL();
 440       }
 441     else
 442       throw INTERP_KERNEL::Exception("Not managed policy detected on lazy side : not implemented !");
 443   }
 444
 445   void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator)
 446   {
 447     MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
 448     _deno_multiply.resize(_coeffs.size());
 449     vector<vector<double> >::iterator iter6=_deno_multiply.begin();
 450     const double *values=source_triangle_surf->getArray()->getConstPointer();
 451     for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
 452       {
 453         (*iter6).resize((*iter4).size());
 454         std::fill((*iter6).begin(),(*iter6).end(),*values);
 455       }
 456     source_triangle_surf->decrRef();
 457     _deno_reverse_multiply=_target_volume;
 458   }
 459
 460   void InterpolationMatrix::computeRevIntegralDenoW(ElementLocator& elementLocator)
 461   {
 462     _deno_multiply=_target_volume;
 463     MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
 464     _deno_reverse_multiply.resize(_coeffs.size());
 465     vector<vector<double> >::iterator iter6=_deno_reverse_multiply.begin();
 466     const double *values=source_triangle_surf->getArray()->getConstPointer();
 467     for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
 468       {
 469         (*iter6).resize((*iter4).size());
 470         std::fill((*iter6).begin(),(*iter6).end(),*values);
 471       }
 472     source_triangle_surf->decrRef();
 473   }
 474
 475   /*!
 476    * Nothing to do because surface computation is on working side.
 477    */
 478   void InterpolationMatrix::computeIntegralDenoL(ElementLocator& elementLocator)
 479   {
 480   }
 481
 482   /*!
 483    * Nothing to do because surface computation is on working side.
 484    */
 485   void InterpolationMatrix::computeRevIntegralDenoL(ElementLocator& elementLocator)
 486   {
 487   }
 488
 489
 490   void InterpolationMatrix::computeGlobConstraintDenoW(ElementLocator& elementLocator)
 491   {
 492     computeGlobalColSum(_deno_multiply);
 493     computeGlobalRowSum(elementLocator,_deno_reverse_multiply,_deno_multiply);
 494   }
 495
 496   void InterpolationMatrix::computeGlobalRowSum(ElementLocator& elementLocator, std::vector<std::vector<double> >& denoStrorage, std::vector<std::vector<double> >& denoStrorageInv)
 497   {
 498     //stores id in distant procs sorted by lazy procs connected with
 499     vector< vector<mcIdType> > rowsPartialSumI;
 500     //stores for each lazy procs connected with, if global info is available and if it's the case the policy
 501     vector<int> policyPartial;
 502     //stores the corresponding values.
 503     vector< vector<double> > rowsPartialSumD;
 504     elementLocator.recvPolicyFromLazySideW(policyPartial);
 505     int pol1=mergePolicies(policyPartial);
 506     if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
 507       {
 508         computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
 509         elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
 510         elementLocator.recvSumFromLazySideW(rowsPartialSumD);
 511       }
 512     else if(pol1==ElementLocator::CUMULATIVE_POLICY)
 513       {
 514         //updateWithNewAdditionnalElements(addingElements);
 515         //stores for each lazy procs connected with, the ids in global mode if it exists (regarding policyPartial). This array has exactly the size of  rowsPartialSumI,
 516         //if policyPartial has CUMALATIVE_POLICY in each.
 517         vector< vector<mcIdType> > globalIdsPartial;
 518         computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
 519         elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
 520         elementLocator.recvCandidatesGlobalIdsFromLazyProcsW(globalIdsPartial);
 521         std::vector< std::vector<mcIdType> > addingElements;
 522         findAdditionnalElements(elementLocator,addingElements,rowsPartialSumI,globalIdsPartial);
 523         addGhostElements(elementLocator.getDistantProcIds(),addingElements);
 524         rowsPartialSumI.clear();
 525         globalIdsPartial.clear();
 526         computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
 527         elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
 528         elementLocator.recvGlobalIdsFromLazyProcsW(rowsPartialSumI,globalIdsPartial);
 529         //
 530         elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
 531         elementLocator.recvSumFromLazySideW(rowsPartialSumD);
 532         mergeRowSum3(globalIdsPartial,rowsPartialSumD);
 533         mergeCoeffs(elementLocator.getDistantProcIds(),rowsPartialSumI,globalIdsPartial,denoStrorageInv);
 534       }
 535     else
 536       throw INTERP_KERNEL::Exception("Not managed policy detected : not implemented !");
 537     divideByGlobalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD,denoStrorage);
 538   }
 539
 540   /*!
 541    * @param distantProcs in parameter that indicates which lazy procs are concerned.
 542    * @param resPerProcI out parameter that must be cleared before calling this method. The size of 1st dimension is equal to the size of 'distantProcs'.
 543    *                    It contains the element ids (2nd dimension) of the corresponding lazy proc.
 544    * @param  resPerProcD out parameter with the same format than 'resPerProcI'. It contains corresponding sum values.
 545    */
 546   void InterpolationMatrix::computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<mcIdType> >& resPerProcI,
 547                                                std::vector<std::vector<double> >& resPerProcD) const
 548   {
 549     resPerProcI.resize(distantProcs.size());
 550     resPerProcD.resize(distantProcs.size());
 551     std::vector<double> res(_col_offsets.size());
 552     for(vector<vector<pair<int,double> > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++)
 553       for(vector<pair<int,double> >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++)
 554         res[(*iter3).first]+=(*iter3).second;
 555     set<int> procsSet;
 556     std::size_t id=-1;
 557     const vector<std::pair<int,mcIdType> >& mapping=_mapping.getSendingIds();
 558     for(vector<std::pair<int,mcIdType> >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++)
 559       {
 560         std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first);
 561         if(isIns.second)
 562           id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin();
 563         resPerProcI[id].push_back((*iter2).second);
 564         resPerProcD[id].push_back(res[iter2-mapping.begin()]);
 565       }
 566   }
 567
 568   /*!
 569    * This method is only usable when CUMULATIVE_POLICY detected. This method finds elements ids (typically nodes) lazy side that
 570    * are not present in columns of 'this' and that should regarding cumulative merge of elements regarding their global ids.
 571    */
 572   void InterpolationMatrix::findAdditionnalElements(ElementLocator& elementLocator, std::vector<std::vector<mcIdType> >& elementsToAdd,
 573                                                     const std::vector<std::vector<mcIdType> >& resPerProcI, const std::vector<std::vector<mcIdType> >& globalIdsPartial)
 574   {
 575     std::set<mcIdType> globalIds;
 576     std::size_t nbLazyProcs=globalIdsPartial.size();
 577     for(std::size_t i=0;i<nbLazyProcs;i++)
 578       globalIds.insert(globalIdsPartial[i].begin(),globalIdsPartial[i].end());
 579     std::vector<mcIdType> tmp(globalIds.size());
 580     std::copy(globalIds.begin(),globalIds.end(),tmp.begin());
 581     globalIds.clear();
 582     elementLocator.sendCandidatesForAddElementsW(tmp);
 583     elementLocator.recvAddElementsFromLazyProcsW(elementsToAdd);
 584   }
 585
 586   void InterpolationMatrix::addGhostElements(const std::vector<int>& distantProcs, const std::vector<std::vector<mcIdType> >& elementsToAdd)
 587   {
 588     std::vector< std::vector< std::map<mcIdType,double> > > data1;
 589     std::vector<int> data2;
 590     serializeMe(data1,data2);
 591     initialize();
 592     std::size_t nbOfDistProcs=distantProcs.size();
 593     for(std::size_t i=0;i<nbOfDistProcs;i++)
 594       {
 595         int procId=distantProcs[i];
 596         const std::vector<mcIdType>& eltsForThisProc=elementsToAdd[i];
 597         if(!eltsForThisProc.empty())
 598           {
 599             std::vector<int>::iterator iter1=std::find(data2.begin(),data2.end(),procId);
 600             std::map<mcIdType,double> *toFeed=0;
 601             if(iter1!=data2.end())
 602               {//to test
 603                 std::size_t rank=iter1-data2.begin();
 604                 toFeed=&(data1[rank].back());
 605               }
 606             else
 607               {
 608                 iter1=std::lower_bound(data2.begin(),data2.end(),procId);
 609                 std::size_t rank=iter1-data2.begin();
 610                 data2.insert(iter1,procId);
 611                 std::vector< std::map<mcIdType,double> > tmp(data1.front().size());
 612                 data1.insert(data1.begin()+rank,tmp);
 613                 toFeed=&(data1[rank].back());
 614               }
 615             for(std::vector<mcIdType>::const_iterator iter2=eltsForThisProc.begin();iter2!=eltsForThisProc.end();iter2++)
 616               (*toFeed)[*iter2]=0.;
 617           }
 618       }
 619     //
 620     nbOfDistProcs=data2.size();
 621     for(std::size_t j=0;j<nbOfDistProcs;j++)
 622       fillDSFromVM(data2[j],0,data1[j],0);
 623   }
 624
 625   int InterpolationMatrix::mergePolicies(const std::vector<int>& policyPartial)
 626   {
 627     if(policyPartial.empty())
 628       return ElementLocator::NO_POST_TREATMENT_POLICY;
 629     int ref=policyPartial[0];
 630      std::vector<int>::const_iterator iter1=std::find_if(policyPartial.begin(),policyPartial.end(),std::bind2nd(std::not_equal_to<int>(),ref));
 631     if(iter1!=policyPartial.end())
 632       {
 633         std::ostringstream msg; msg << "Incompatible policies between lazy procs each other : proc # " << iter1-policyPartial.begin();
 634         throw INTERP_KERNEL::Exception(msg.str().c_str());
 635       }
 636     return ref;
 637   }
 638
 639   /*!
 640    * This method introduce global ids aspects in computed 'rowsPartialSumD'.
 641    * As precondition rowsPartialSumD.size()==policyPartial.size()==globalIdsPartial.size(). Foreach i in [0;rowsPartialSumD.size() ) rowsPartialSumD[i].size()==globalIdsPartial[i].size()
 642    * @param rowsPartialSumD : in parameter, Partial row sum computed for each lazy procs connected with.
 643    * @param globalIdsPartial : in parameter, the global numbering, of elements connected with.
 644    * @param globalIdsLazySideInteraction : out parameter, constituted from all global ids of lazy procs connected with.
 645    * @para sumCorresponding : out parameter, relative to 'globalIdsLazySideInteraction'
 646    */
 647   void InterpolationMatrix::mergeRowSum(const std::vector< std::vector<double> >& rowsPartialSumD, const std::vector< std::vector<mcIdType> >& globalIdsPartial,
 648                                         std::vector<mcIdType>& globalIdsLazySideInteraction, std::vector<double>& sumCorresponding)
 649   {
 650     std::map<mcIdType,double> sumToReturn;
 651     std::size_t nbLazyProcs=rowsPartialSumD.size();
 652     for(std::size_t i=0;i<nbLazyProcs;i++)
 653       {
 654         const std::vector<double>& rowSumOfP=rowsPartialSumD[i];
 655         const std::vector<mcIdType>& globalIdsOfP=globalIdsPartial[i];
 656         std::vector<double>::const_iterator iter1=rowSumOfP.begin();
 657         std::vector<mcIdType>::const_iterator iter2=globalIdsOfP.begin();
 658         for(;iter1!=rowSumOfP.end();iter1++,iter2++)
 659           sumToReturn[*iter2]+=*iter1;
 660       }
 661     //
 662     std::size_t lgth=sumToReturn.size();
 663     globalIdsLazySideInteraction.resize(lgth);
 664     sumCorresponding.resize(lgth);
 665     std::vector<mcIdType>::iterator iter3=globalIdsLazySideInteraction.begin();
 666     std::vector<double>::iterator iter4=sumCorresponding.begin();
 667     for(std::map<mcIdType,double>::const_iterator iter5=sumToReturn.begin();iter5!=sumToReturn.end();iter5++,iter3++,iter4++)
 668       {
 669         *iter3=(*iter5).first;
 670         *iter4=(*iter5).second;
 671       }
 672   }
 673
 674   /*!
 675    * This method simply reorganize the result contained in 'sumCorresponding' computed by lazy side into 'rowsPartialSumD' with help of 'globalIdsPartial' and 'globalIdsLazySideInteraction'
 676    *
 677    * @param globalIdsPartial : in parameter, global ids sorted by lazy procs
 678    * @param rowsPartialSumD : in/out parameter, with exactly the same size as 'globalIdsPartial'
 679    * @param globalIdsLazySideInteraction : in parameter that represents ALL the global ids of every lazy procs in interaction
 680    * @param sumCorresponding : in parameter with same size as 'globalIdsLazySideInteraction' that stores the corresponding sum of 'globalIdsLazySideInteraction'
 681    */
 682   void InterpolationMatrix::mergeRowSum2(const std::vector< std::vector<mcIdType> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD,
 683                                          const std::vector<mcIdType>& globalIdsLazySideInteraction, const std::vector<double>& sumCorresponding)
 684   {
 685     std::map<mcIdType,double> acc;
 686     std::vector<mcIdType>::const_iterator iter1=globalIdsLazySideInteraction.begin();
 687     std::vector<double>::const_iterator iter2=sumCorresponding.begin();
 688     for(;iter1!=globalIdsLazySideInteraction.end();iter1++,iter2++)
 689       acc[*iter1]=*iter2;
 690     //
 691     std::size_t nbLazyProcs=globalIdsPartial.size();
 692     for(std::size_t i=0;i<nbLazyProcs;i++)
 693       {
 694         const std::vector<mcIdType>& tmp1=globalIdsPartial[i];
 695         std::vector<double>& tmp2=rowsPartialSumD[i];
 696         std::vector<mcIdType>::const_iterator iter3=tmp1.begin();
 697         std::vector<double>::iterator iter4=tmp2.begin();
 698         for(;iter3!=tmp1.end();iter3++,iter4++)
 699           *iter4=acc[*iter3];
 700       }
 701   }
 702
 703   void InterpolationMatrix::mergeRowSum3(const std::vector< std::vector<mcIdType> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD)
 704   {
 705     std::map<mcIdType,double> sum;
 706     std::vector< std::vector<mcIdType> >::const_iterator iter1=globalIdsPartial.begin();
 707     std::vector< std::vector<double> >::iterator iter2=rowsPartialSumD.begin();
 708     for(;iter1!=globalIdsPartial.end();iter1++,iter2++)
 709       {
 710         std::vector<mcIdType>::const_iterator iter3=(*iter1).begin();
 711         std::vector<double>::const_iterator iter4=(*iter2).begin();
 712         for(;iter3!=(*iter1).end();iter3++,iter4++)
 713           sum[*iter3]+=*iter4;
 714       }
 715     iter2=rowsPartialSumD.begin();
 716     for(iter1=globalIdsPartial.begin();iter1!=globalIdsPartial.end();iter1++,iter2++)
 717       {
 718         std::vector<mcIdType>::const_iterator iter3=(*iter1).begin();
 719         std::vector<double>::iterator iter4=(*iter2).begin();
 720         for(;iter3!=(*iter1).end();iter3++,iter4++)
 721           *iter4=sum[*iter3];
 722       }
 723   }
 724
 725   /*!
 726    * This method updates this->_coeffs attribute in order to take into account hidden (because having the same global number) similar nodes in _coeffs array.
 727    * If in this->_coeffs two distant element id have the same global id their values will be replaced for each by the sum of the two.
 728    * @param procsInInteraction input parameter : specifies the procId in absolute of distant lazy procs in interaction with
 729    * @param rowsPartialSumI input parameter : local ids of distant lazy procs elements in interaction with
 730    * @param globalIdsPartial input parameter : global ids of distant lazy procs elements in interaction with
 731    */
 732   void InterpolationMatrix::mergeCoeffs(const std::vector<int>& procsInInteraction, const std::vector< std::vector<mcIdType> >& rowsPartialSumI,
 733                                         const std::vector<std::vector<mcIdType> >& globalIdsPartial, std::vector<std::vector<double> >& denoStrorageInv)
 734   {
 735     //preparing fast access structures
 736     std::map<int,int> procT;
 737     int localProcId=0;
 738     for(std::vector<int>::const_iterator iter1=procsInInteraction.begin();iter1!=procsInInteraction.end();iter1++,localProcId++)
 739       procT[*iter1]=localProcId;
 740     std::size_t size=procsInInteraction.size();
 741     std::vector<std::map<mcIdType,mcIdType> > localToGlobal(size);
 742     for(std::size_t i=0;i<size;i++)
 743       {
 744         std::map<mcIdType,mcIdType>& myLocalToGlobal=localToGlobal[i];
 745         const std::vector<mcIdType>& locals=rowsPartialSumI[i];
 746         const std::vector<mcIdType>& globals=globalIdsPartial[i];
 747         std::vector<mcIdType>::const_iterator iter3=locals.begin();
 748         std::vector<mcIdType>::const_iterator iter4=globals.begin();
 749         for(;iter3!=locals.end();iter3++,iter4++)
 750           myLocalToGlobal[*iter3]=*iter4;
 751       }
 752     //
 753     const vector<std::pair<int,mcIdType> >& mapping=_mapping.getSendingIds();
 754     std::map<mcIdType,double> globalIdVal;
 755     //accumulate for same global id on lazy part.
 756     for(vector<vector<pair<int,double> > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++)
 757       for(vector<pair<int,double> >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 758         {
 759           const std::pair<int,mcIdType>& distantLocalLazyId=mapping[(*iter2).first];
 760           int localLazyProcId=procT[distantLocalLazyId.first];
 761           mcIdType globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second];
 762           globalIdVal[globalDistantLazyId]+=(*iter2).second;
 763         }
 764     //perform merge
 765     std::vector<std::vector<double> >::iterator iter3=denoStrorageInv.begin();
 766     for(vector<vector<pair<int,double> > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter3++)
 767       {
 768         double val=(*iter3).back();
 769         (*iter3).resize((*iter1).size());
 770         std::vector<double>::iterator iter4=(*iter3).begin();
 771         for(vector<pair<int,double> >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter4++)
 772           {
 773             const std::pair<int,mcIdType>& distantLocalLazyId=mapping[(*iter2).first];
 774             int localLazyProcId=procT[distantLocalLazyId.first];
 775             mcIdType globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second];
 776             double newVal=globalIdVal[globalDistantLazyId];
 777             if((*iter2).second!=0.)
 778               (*iter4)=val*newVal/(*iter2).second;
 779             else
 780               (*iter4)=std::numeric_limits<double>::max();
 781             (*iter2).second=newVal;
 782           }
 783       }
 784   }
 785
 786   void InterpolationMatrix::divideByGlobalRowSum(const std::vector<int>& distantProcs, const std::vector<std::vector<mcIdType> >& resPerProcI,
 787                                                  const std::vector<std::vector<double> >& resPerProcD, std::vector<std::vector<double> >& deno)
 788   {
 789     map<mcIdType,double> fastSums;
 790     int procId=0;
 791     for(vector<int>::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++)
 792       {
 793         const std::vector<mcIdType>& currentProcI=resPerProcI[procId];
 794         const std::vector<double>& currentProcD=resPerProcD[procId];
 795         vector<double>::const_iterator iter3=currentProcD.begin();
 796         for(vector<mcIdType>::const_iterator iter2=currentProcI.begin();iter2!=currentProcI.end();iter2++,iter3++)
 797           fastSums[_col_offsets[std::make_pair(*iter1,*iter2)]]=*iter3;
 798       }
 799     deno.resize(_coeffs.size());
 800     vector<vector<double> >::iterator iter6=deno.begin();
 801     for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++)
 802       {
 803         (*iter6).resize((*iter4).size());
 804         vector<double>::iterator iter7=(*iter6).begin();
 805         for(vector<pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++,iter7++)
 806           *iter7=fastSums[(*iter5).first];
 807       }
 808   }
 809
 810   void InterpolationMatrix::computeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
 811   {
 812     denoStrorage.resize(_coeffs.size());
 813     vector<vector<double> >::iterator iter2=denoStrorage.begin();
 814     for(vector<vector<pair<int,double> > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++)
 815       {
 816         (*iter2).resize((*iter1).size());
 817         double sumOfCurrentRow=0.;
 818         for(vector<pair<int,double> >::const_iterator iter3=(*iter1).begin();iter3!=(*iter1).end();iter3++)
 819           sumOfCurrentRow+=(*iter3).second;
 820         std::fill((*iter2).begin(),(*iter2).end(),sumOfCurrentRow);
 821       }
 822   }
 823
 824   void InterpolationMatrix::resizeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
 825   {
 826     vector<vector<double> >::iterator iter2=denoStrorage.begin();
 827     for(vector<vector<pair<int,double> > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++)
 828       {
 829         double val=(*iter2).back();
 830         (*iter2).resize((*iter1).size());
 831         std::fill((*iter2).begin(),(*iter2).end(),val);
 832       }
 833   }
 834
 835
 836  /**!  The call to this method updates the arrays on the target side
 837      so that they know which amount of data from which processor they
 838      should expect.
 839      That call makes actual interpolations via multiply method
 840      available.
 841      */
 842   void InterpolationMatrix::prepare()
 843   {
 844     mcIdType nbelems = _source_field->getField()->getNumberOfTuples();
 845     for (mcIdType ielem=0; ielem < nbelems; ielem++)
 846       {
 847         _row_offsets[ielem+1]+=_row_offsets[ielem];
 848       }
 849     _mapping.prepareSendRecv();
 850   }
 851
 852   MCAuto<DataArrayIdType> InterpolationMatrix::retrieveNonFetchedIdsTarget(mcIdType nbTuples) const
 853   {
 854     return _mapping.retrieveNonFetchedIdsTarget(nbTuples);
 855   }
 856
 857   /*!
 858      \brief performs t=Ws, where t is the target field, s is the source field
 859
 860      The call to this method must be called both on the working side
 861      and on the idle side. On the working side, the vector  T=VT^(-1).(W.S)
 862      is computed and sent. On the idle side, no computation is done, but the
 863      result from the working side is received and the field is updated.
 864
 865      \param field source field on processors involved on the source side,
 866      target field on processors on the target side
 867    */
 868   void InterpolationMatrix::multiply(MEDCouplingFieldDouble& field) const
 869   {
 870     mcIdType nbcomp = ToIdType(field.getArray()->getNumberOfComponents());
 871     vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
 872     //computing the matrix multiply on source side
 873     if (_source_group.containsMyRank())
 874       {
 875         mcIdType nbrows = ToIdType(_coeffs.size());
 876         // performing W.S
 877         // W is the intersection matrix
 878         // S is the source vector
 879
 880         for (mcIdType irow=0; irow<nbrows; irow++)
 881           {
 882             for (mcIdType icomp=0; icomp< nbcomp; icomp++)
 883               {
 884                 double coeff_row = field.getIJ(irow,icomp);
 885                 for (mcIdType icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
 886                   {
 887                     int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
 888                     double value = _coeffs[irow][icol-_row_offsets[irow]].second;
 889                     double deno = _deno_multiply[irow][icol-_row_offsets[irow]];
 890                     target_value[colid*nbcomp+icomp]+=value*coeff_row/deno;
 891                   }
 892               }
 893           }
 894       }
 895
 896     if (_target_group.containsMyRank())
 897       {
 898         field.getArray()->fillWithZero();
 899       }
 900
 901     //on source side : sending  T=VT^(-1).(W.S)
 902     //on target side :: receiving T and storing it in field
 903     _mapping.sendRecv(target_value.data(),field);
 904
 905     if( _target_group.containsMyRank() )
 906     {
 907       if( this->_presence_dft_value )
 908       {
 909         const MCAuto<DataArrayIdType> nonFetchedEntities = _mapping.retrieveNonFetchedIdsTarget(field.getArray()->getNumberOfTuples());
 910         double *fieldPtr( field.getArray()->getPointerSilent() );
 911         for( const mcIdType *eltId = nonFetchedEntities->begin() ; eltId != nonFetchedEntities->end() ; ++eltId)
 912           std::fill( fieldPtr + (*eltId)*nbcomp, fieldPtr + ((*eltId)+1)*nbcomp, this->_dft_value );
 913       }
 914     }
 915   }
 916
 917   MCAuto<DataArrayIdType> InterpolationMatrix::retrieveNonFetchedIdsSource() const
 918   {
 919     MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
 920     if (_source_group.containsMyRank())
 921       {
 922         mcIdType nbrows = ToIdType( _coeffs.size() );
 923         for (mcIdType irow = 0; irow < nbrows; irow++)
 924           {
 925             if( _row_offsets[irow+1] == _row_offsets[irow] )
 926             {
 927               ret->pushBackSilent( irow );
 928             }
 929           }
 930       }
 931     else
 932       THROW_IK_EXCEPTION("InterpolationMatrix::retrieveNonFetchedIdsSource : not supposed to be called target side !");
 933     return ret;
 934   }
 935
 936
 937   /**!
 938    \brief performs s=WTt, where t is the target field, s is the source field,
 939    WT is the transpose matrix from W
 940
 941      The call to this method must be called both on the working side
 942      and on the idle side. On the working side, the target vector T is
 943      received and the vector  S=VS^(-1).(WT.T) is computed to update
 944      the field.
 945      On the idle side, no computation is done, but the field is sent.
 946
 947      param field source field on processors involved on the source side,
 948      target field on processors on the target side
 949      */
 950   void InterpolationMatrix::transposeMultiply(MEDCouplingFieldDouble& field) const
 951   {
 952     std::size_t nbcomp = field.getArray()->getNumberOfComponents();
 953     vector<double> source_value(_col_offsets.size()* nbcomp,0.0);
 954     _mapping.reverseSendRecv(&source_value[0],field);
 955
 956     //treatment of the transpose matrix multiply on the source side
 957     if (_source_group.containsMyRank())
 958       {
 959         mcIdType nbrows = ToIdType( _coeffs.size() );
 960         double   *array = field.getArray()->getPointer() ;
 961
 962         // Initialization
 963         std::fill(array, array+nbrows*nbcomp, 0.0) ;
 964
 965         //performing WT.T
 966         //WT is W transpose
 967         //T is the target vector
 968         for (mcIdType irow = 0; irow < nbrows; irow++)
 969           {
 970             if( _row_offsets[irow+1] > _row_offsets[irow] )
 971             {
 972               for (mcIdType icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
 973                 {
 974                   int colid    = _coeffs[irow][icol-_row_offsets[irow]].first;
 975                   double value = _coeffs[irow][icol-_row_offsets[irow]].second;
 976                   double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
 977                   for (std::size_t icomp=0; icomp<nbcomp; icomp++)
 978                     {
 979                       double coeff_row = source_value[colid*nbcomp+icomp];
 980                       array[irow*nbcomp+icomp] += value*coeff_row/deno;
 981                     }
 982                 }
 983             }
 984             else
 985             {
 986               if( _presence_dft_value )
 987                 std::fill(array+irow*nbcomp,array+(irow+1)*nbcomp,this->_dft_value);
 988             }
 989           }
 990       }
 991   }
 992
 993   bool InterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const
 994   {
 995     return method=="P0";
 996   }
 997 }