src/ParaMEDMEM/InterpolationMatrix.cxx

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