src/ParaMEDMEM/InterpolationMatrix.cxx

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