Back to index

salome-med  6.5.0
InterpolationMatrix.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "ParaMESH.hxx"
00021 #include "ParaFIELD.hxx"
00022 #include "ProcessorGroup.hxx"
00023 #include "MxN_Mapping.hxx"
00024 #include "InterpolationMatrix.hxx"
00025 #include "TranslationRotationMatrix.hxx"
00026 #include "Interpolation.hxx"
00027 #include "Interpolation1D.txx"
00028 #include "Interpolation2DCurve.hxx"
00029 #include "Interpolation2D.txx"
00030 #include "Interpolation3DSurf.hxx"
00031 #include "Interpolation3D.txx"
00032 #include "Interpolation3D2D.txx"
00033 #include "Interpolation2D1D.txx"
00034 #include "MEDCouplingUMesh.hxx"
00035 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
00036 #include "InterpolationOptions.hxx"
00037 #include "NormalizedUnstructuredMesh.hxx"
00038 #include "ElementLocator.hxx"
00039 
00040 #include <algorithm>
00041 
00042 // class InterpolationMatrix
00043 // This class enables the storage of an interpolation matrix Wij mapping 
00044 // source field Sj to target field Ti via Ti=Vi^(-1).Wij.Sj.
00045 // The matrix is built and stored on the processors belonging to the source
00046 // group. 
00047 
00048 using namespace std;
00049 
00050 namespace ParaMEDMEM
00051 {
00052 
00053   //   ====================================================================
00054   //   Creates an empty matrix structure linking two distributed supports.
00055   //   The method must be called by all processors belonging to source
00056   //   and target groups.
00057   //   param source_support local support
00058   //   param source_group processor group containing the local processors
00059   //   param target_group processor group containing the distant processors
00060   //   param method interpolation method
00061   //   ====================================================================
00062 
00063   InterpolationMatrix::InterpolationMatrix(const ParaMEDMEM::ParaFIELD *source_field, 
00064                                            const ProcessorGroup& source_group,
00065                                            const ProcessorGroup& target_group,
00066                                            const DECOptions& dec_options,
00067                                            const INTERP_KERNEL::InterpolationOptions& interp_options):
00068     INTERP_KERNEL::InterpolationOptions(interp_options),
00069     DECOptions(dec_options),
00070     _source_field(source_field),
00071     _source_support(source_field->getSupport()->getCellMesh()),
00072     _mapping(source_group, target_group, dec_options),
00073     _source_group(source_group),
00074     _target_group(target_group)
00075   {
00076     int nbelems = source_field->getField()->getNumberOfTuples();
00077     _row_offsets.resize(nbelems+1);
00078     _coeffs.resize(nbelems);
00079     _target_volume.resize(nbelems);
00080   }
00081 
00082   InterpolationMatrix::~InterpolationMatrix()
00083   {
00084   }
00085 
00086 
00087   //   ======================================================================
00088   //   \brief Adds the contribution of a distant subdomain to the*
00089   //   interpolation matrix.
00090   //   The method adds contribution to the interpolation matrix.
00091   //   For each row of the matrix, elements are addded as
00092   //   a (column, coeff) pair in the _coeffs array. This column number refers
00093   //   to an element on the target side via the _col_offsets array.
00094   //   It is made of a series of (iproc, ielem) pairs. 
00095   //   The number of elements per row is stored in the row_offsets array.
00096 
00097   //   param distant_support local representation of the distant subdomain
00098   //   param iproc_distant id of the distant subdomain (in the distant group)
00099   //   param distant_elems mapping between the local representation of
00100   //   the subdomain and the actual elem ids on the distant subdomain
00101   //   ======================================================================
00102 
00103   void InterpolationMatrix::addContribution ( MEDCouplingPointSet& distant_support,
00104                                               int iproc_distant,
00105                                               const int* distant_elems,
00106                                               const std::string& srcMeth,
00107                                               const std::string& targetMeth)
00108   {
00109     std::string interpMethod(targetMeth);
00110     interpMethod+=srcMeth;
00111     //creating the interpolator structure
00112     vector<map<int,double> > surfaces;
00113     int colSize=0;
00114     //computation of the intersection volumes between source and target elements
00115     MEDCouplingUMesh *distant_supportC=dynamic_cast<MEDCouplingUMesh *>(&distant_support);
00116     MEDCouplingUMesh *source_supportC=dynamic_cast<MEDCouplingUMesh *>(_source_support);
00117     if ( distant_support.getMeshDimension() == -1 )
00118       {
00119         if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==2)
00120           {
00121             MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(source_supportC);
00122             INTERP_KERNEL::Interpolation2D interpolation(*this);
00123             colSize=interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth.c_str());
00124           }
00125         else if(source_supportC->getMeshDimension()==3 && source_supportC->getSpaceDimension()==3)
00126           {
00127             MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(source_supportC);
00128             INTERP_KERNEL::Interpolation3D interpolation(*this);
00129             colSize=interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth.c_str());
00130           }
00131         else if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==3)
00132           {
00133             MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(source_supportC);
00134             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00135             colSize=interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth.c_str());
00136           }
00137         else
00138           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
00139       }
00140     else if ( source_supportC->getMeshDimension() == -1 )
00141       {
00142         if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==2)
00143           {
00144             MEDCouplingNormalizedUnstructuredMesh<2,2> distant_mesh_wrapper(distant_supportC);
00145             INTERP_KERNEL::Interpolation2D interpolation(*this);
00146             colSize=interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth.c_str());
00147           }
00148         else if(distant_supportC->getMeshDimension()==3 && distant_supportC->getSpaceDimension()==3)
00149           {
00150             MEDCouplingNormalizedUnstructuredMesh<3,3> distant_mesh_wrapper(distant_supportC);
00151             INTERP_KERNEL::Interpolation3D interpolation(*this);
00152             colSize=interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth.c_str());
00153           }
00154         else if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==3)
00155           {
00156             MEDCouplingNormalizedUnstructuredMesh<3,2> distant_mesh_wrapper(distant_supportC);
00157             INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00158             colSize=interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth.c_str());
00159           }
00160         else
00161           throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
00162       }
00163     else if ( distant_support.getMeshDimension() == 2
00164               && _source_support->getMeshDimension() == 3
00165               && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
00166       {
00167         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
00168         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
00169         INTERP_KERNEL::Interpolation3D2D interpolator (*this);
00170         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
00171         target_wrapper.releaseTempArrays();
00172         source_wrapper.releaseTempArrays();
00173       }
00174     else if ( distant_support.getMeshDimension() == 1
00175               && _source_support->getMeshDimension() == 2
00176               && distant_support.getSpaceDimension() == 2 && _source_support->getSpaceDimension() == 2)
00177       {
00178         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
00179         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
00180         INTERP_KERNEL::Interpolation2D1D interpolator (*this);
00181         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
00182         target_wrapper.releaseTempArrays();
00183         source_wrapper.releaseTempArrays();
00184       }
00185     else if (distant_support.getMeshDimension() != _source_support->getMeshDimension())
00186       {
00187         throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
00188       }
00189     else if( distant_support.getMeshDimension() == 1
00190              && distant_support.getSpaceDimension() == 1 )
00191       {
00192         MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(distant_supportC);
00193         MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(source_supportC);
00194 
00195         INTERP_KERNEL::Interpolation1D interpolation(*this);
00196         colSize=interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
00197         target_wrapper.releaseTempArrays();
00198         source_wrapper.releaseTempArrays();
00199       }
00200     else if( distant_support.getMeshDimension() == 1
00201              && distant_support.getSpaceDimension() == 2 )
00202       {
00203         MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(distant_supportC);
00204         MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(source_supportC);
00205 
00206         INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
00207         colSize=interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
00208         target_wrapper.releaseTempArrays();
00209         source_wrapper.releaseTempArrays();
00210       }
00211     else if ( distant_support.getMeshDimension() == 2
00212               && distant_support.getSpaceDimension() == 3 )
00213       {
00214         MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC);
00215         MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC);
00216 
00217         INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
00218         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
00219         target_wrapper.releaseTempArrays();
00220         source_wrapper.releaseTempArrays();
00221       }
00222     else if ( distant_support.getMeshDimension() == 2
00223               && distant_support.getSpaceDimension() == 2)
00224       {
00225         MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
00226         MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
00227 
00228         INTERP_KERNEL::Interpolation2D interpolator (*this);
00229         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
00230         target_wrapper.releaseTempArrays();
00231         source_wrapper.releaseTempArrays();
00232       }
00233     else if ( distant_support.getMeshDimension() == 3
00234               && distant_support.getSpaceDimension() == 3 )
00235       {
00236         MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
00237         MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
00238 
00239         INTERP_KERNEL::Interpolation3D interpolator (*this);
00240         colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
00241         target_wrapper.releaseTempArrays();
00242         source_wrapper.releaseTempArrays();
00243       }
00244     else
00245       {
00246         throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
00247       }
00248     bool needTargetSurf=isSurfaceComputationNeeded(targetMeth);
00249 
00250     MEDCouplingFieldDouble *target_triangle_surf=0;
00251     if(needTargetSurf)
00252       target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus());
00253     fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf);
00254 
00255     if(needTargetSurf)
00256       target_triangle_surf->decrRef();
00257   }
00258 
00259   void InterpolationMatrix::fillDSFromVM(int iproc_distant, const int* distant_elems, const std::vector< std::map<int,double> >& values, MEDCouplingFieldDouble *surf)
00260   {
00261     //loop over the elements to build the interpolation
00262     //matrix structures
00263     int source_size=values.size();
00264     for (int ielem=0; ielem < source_size; ielem++) 
00265       {
00266         _row_offsets[ielem+1] += values[ielem].size();
00267         for(map<int,double>::const_iterator iter=values[ielem].begin();iter!=values[ielem].end();iter++)
00268           {
00269             int localId;
00270             if(distant_elems)
00271               localId=distant_elems[iter->first];
00272             else
00273               localId=iter->first;
00274             //locating the (iproc, itriangle) pair in the list of columns
00275             map<pair<int,int>,int >::iterator iter2 = _col_offsets.find(make_pair(iproc_distant,localId));
00276             int col_id;
00277 
00278             if (iter2 == _col_offsets.end())
00279               {
00280                 //(iproc, itriangle) is not registered in the list
00281                 //of distant elements
00282                 col_id =_col_offsets.size();
00283                 _col_offsets.insert(make_pair(make_pair(iproc_distant,localId),col_id));
00284                 _mapping.addElementFromSource(iproc_distant,localId);
00285               }
00286             else 
00287               {
00288                 col_id = iter2->second;
00289               }
00290             //the non zero coefficient is stored 
00291             //ielem is the row,
00292             //col_id is the number of the column
00293             //iter->second is the value of the coefficient
00294             if(surf)
00295               {
00296                 double surface = surf->getIJ(iter->first,0);
00297                 _target_volume[ielem].push_back(surface);
00298               }
00299             _coeffs[ielem].push_back(make_pair(col_id,iter->second));
00300           }
00301       }
00302   }
00303 
00304   void InterpolationMatrix::serializeMe(std::vector< std::vector< std::map<int,double> > >& data1, std::vector<int>& data2) const
00305   {
00306     data1.clear();
00307     data2.clear();
00308     const std::vector<std::pair<int,int> >& sendingIds=_mapping.getSendingIds();
00309     std::set<int> procsS;
00310     for(std::vector<std::pair<int,int> >::const_iterator iter1=sendingIds.begin();iter1!=sendingIds.end();iter1++)
00311       procsS.insert((*iter1).first);
00312     data1.resize(procsS.size());
00313     data2.resize(procsS.size());
00314     std::copy(procsS.begin(),procsS.end(),data2.begin());
00315     std::map<int,int> fastProcAcc;
00316     int id=0;
00317     for(std::set<int>::const_iterator iter2=procsS.begin();iter2!=procsS.end();iter2++,id++)
00318       fastProcAcc[*iter2]=id;
00319     int nbOfSrcElt=_coeffs.size();
00320     for(std::vector< std::vector< std::map<int,double> > >::iterator iter3=data1.begin();iter3!=data1.end();iter3++)
00321       (*iter3).resize(nbOfSrcElt);
00322     id=0;
00323     for(std::vector< std::vector< std::pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,id++)
00324       {
00325         for(std::vector< std::pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++)
00326           {
00327             const std::pair<int,int>& elt=sendingIds[(*iter5).first];
00328             data1[fastProcAcc[elt.first]][id][elt.second]=(*iter5).second;
00329           }
00330       }
00331   }
00332 
00333   void InterpolationMatrix::initialize()
00334   {
00335     int lgth=_coeffs.size();
00336     _row_offsets.clear(); _row_offsets.resize(lgth+1);
00337     _coeffs.clear(); _coeffs.resize(lgth);
00338     _target_volume.clear(); _target_volume.resize(lgth);
00339     _col_offsets.clear();
00340     _mapping.initialize();
00341   }
00342 
00343   void InterpolationMatrix::finishContributionW(ElementLocator& elementLocator)
00344   {
00345     NatureOfField nature=elementLocator.getLocalNature();
00346     switch(nature)
00347       {
00348       case ConservativeVolumic:
00349         computeConservVolDenoW(elementLocator);
00350         break;
00351       case Integral:
00352         {
00353           if(!elementLocator.isM1DCorr())
00354             computeIntegralDenoW(elementLocator);
00355           else
00356             computeGlobConstraintDenoW(elementLocator);
00357           break;
00358         }
00359       case IntegralGlobConstraint:
00360         computeGlobConstraintDenoW(elementLocator);
00361         break;
00362       case RevIntegral:
00363         {
00364           if(!elementLocator.isM1DCorr())
00365             computeRevIntegralDenoW(elementLocator);
00366           else
00367             computeConservVolDenoW(elementLocator);
00368           break;
00369         }
00370       default:
00371         throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
00372         break;
00373       }
00374   }
00375 
00376   void InterpolationMatrix::finishContributionL(ElementLocator& elementLocator)
00377   {
00378     NatureOfField nature=elementLocator.getLocalNature();
00379     switch(nature)
00380       {
00381       case ConservativeVolumic:
00382         computeConservVolDenoL(elementLocator);
00383         break;
00384       case Integral:
00385         {
00386           if(!elementLocator.isM1DCorr())
00387             computeIntegralDenoL(elementLocator);
00388           else
00389             computeConservVolDenoL(elementLocator);
00390           break;
00391         }
00392       case IntegralGlobConstraint:
00393         //this is not a bug doing like ConservativeVolumic
00394         computeConservVolDenoL(elementLocator);
00395         break;
00396       case RevIntegral:
00397         {
00398           if(!elementLocator.isM1DCorr())
00399             computeRevIntegralDenoL(elementLocator);
00400           else
00401             computeConservVolDenoL(elementLocator);
00402           break;
00403         }
00404       default:
00405         throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
00406         break;
00407       }
00408   }
00409   
00410   void InterpolationMatrix::computeConservVolDenoW(ElementLocator& elementLocator)
00411   {
00412     computeGlobalColSum(_deno_reverse_multiply);
00413     computeGlobalRowSum(elementLocator,_deno_multiply,_deno_reverse_multiply);
00414   }
00415   
00416   void InterpolationMatrix::computeConservVolDenoL(ElementLocator& elementLocator)
00417   {
00418     int pol1=elementLocator.sendPolicyToWorkingSideL();
00419     if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
00420       {
00421         elementLocator.recvFromWorkingSideL();
00422         elementLocator.sendToWorkingSideL();
00423       }
00424     else if(ElementLocator::CUMULATIVE_POLICY)
00425       {
00426         //ask for lazy side to deduce ids eventually missing on working side and to send it back.
00427         elementLocator.recvLocalIdsFromWorkingSideL();
00428         elementLocator.sendCandidatesGlobalIdsToWorkingSideL();
00429         elementLocator.recvCandidatesForAddElementsL();
00430         elementLocator.sendAddElementsToWorkingSideL();
00431         //Working side has updated its eventually missing ids updates its global ids with lazy side procs contribution
00432         elementLocator.recvLocalIdsFromWorkingSideL();
00433         elementLocator.sendGlobalIdsToWorkingSideL();
00434         //like no post treatment
00435         elementLocator.recvFromWorkingSideL();
00436         elementLocator.sendToWorkingSideL();
00437       }
00438     else
00439       throw INTERP_KERNEL::Exception("Not managed policy detected on lazy side : not implemented !");
00440   }
00441 
00442   void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator)
00443   {
00444     MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
00445     _deno_multiply.resize(_coeffs.size());
00446     vector<vector<double> >::iterator iter6=_deno_multiply.begin();
00447     const double *values=source_triangle_surf->getArray()->getConstPointer();
00448     for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
00449       {
00450         (*iter6).resize((*iter4).size());
00451         std::fill((*iter6).begin(),(*iter6).end(),*values);
00452       }
00453     source_triangle_surf->decrRef();
00454     _deno_reverse_multiply=_target_volume;
00455   }
00456 
00457   void InterpolationMatrix::computeRevIntegralDenoW(ElementLocator& elementLocator)
00458   {
00459     _deno_multiply=_target_volume;
00460     MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
00461     _deno_reverse_multiply.resize(_coeffs.size());
00462     vector<vector<double> >::iterator iter6=_deno_reverse_multiply.begin();
00463     const double *values=source_triangle_surf->getArray()->getConstPointer();
00464     for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
00465       {
00466         (*iter6).resize((*iter4).size());
00467         std::fill((*iter6).begin(),(*iter6).end(),*values);
00468       }
00469     source_triangle_surf->decrRef();
00470   }
00471   
00475   void InterpolationMatrix::computeIntegralDenoL(ElementLocator& elementLocator)
00476   {
00477   }
00478 
00482   void InterpolationMatrix::computeRevIntegralDenoL(ElementLocator& elementLocator)
00483   {
00484   }
00485 
00486 
00487   void InterpolationMatrix::computeGlobConstraintDenoW(ElementLocator& elementLocator)
00488   {
00489     computeGlobalColSum(_deno_multiply);
00490     computeGlobalRowSum(elementLocator,_deno_reverse_multiply,_deno_multiply);
00491   }
00492 
00493   void InterpolationMatrix::computeGlobalRowSum(ElementLocator& elementLocator, std::vector<std::vector<double> >& denoStrorage, std::vector<std::vector<double> >& denoStrorageInv)
00494   {
00495     //stores id in distant procs sorted by lazy procs connected with
00496     vector< vector<int> > rowsPartialSumI;
00497     //stores for each lazy procs connected with, if global info is available and if it's the case the policy
00498     vector<int> policyPartial;
00499     //stores the corresponding values.
00500     vector< vector<double> > rowsPartialSumD;
00501     elementLocator.recvPolicyFromLazySideW(policyPartial);
00502     int pol1=mergePolicies(policyPartial);
00503     if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
00504       {
00505         computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
00506         elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
00507         elementLocator.recvSumFromLazySideW(rowsPartialSumD);
00508       }
00509     else if(pol1==ElementLocator::CUMULATIVE_POLICY)
00510       {
00511         //updateWithNewAdditionnalElements(addingElements);
00512         //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,
00513         //if policyPartial has CUMALATIVE_POLICY in each.
00514         vector< vector<int> > globalIdsPartial;
00515         computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
00516         elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
00517         elementLocator.recvCandidatesGlobalIdsFromLazyProcsW(globalIdsPartial);
00518         std::vector< std::vector<int> > addingElements;
00519         findAdditionnalElements(elementLocator,addingElements,rowsPartialSumI,globalIdsPartial);
00520         addGhostElements(elementLocator.getDistantProcIds(),addingElements);
00521         rowsPartialSumI.clear();
00522         globalIdsPartial.clear();
00523         computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
00524         elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
00525         elementLocator.recvGlobalIdsFromLazyProcsW(rowsPartialSumI,globalIdsPartial);
00526         //
00527         elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
00528         elementLocator.recvSumFromLazySideW(rowsPartialSumD);
00529         mergeRowSum3(globalIdsPartial,rowsPartialSumD);
00530         mergeCoeffs(elementLocator.getDistantProcIds(),rowsPartialSumI,globalIdsPartial,denoStrorageInv);
00531       }
00532     else
00533       throw INTERP_KERNEL::Exception("Not managed policy detected : not implemented !");
00534     divideByGlobalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD,denoStrorage);
00535   }
00536 
00543   void InterpolationMatrix::computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<int> >& resPerProcI,
00544                                                std::vector<std::vector<double> >& resPerProcD) const
00545   {
00546     resPerProcI.resize(distantProcs.size());
00547     resPerProcD.resize(distantProcs.size());
00548     std::vector<double> res(_col_offsets.size());
00549     for(vector<vector<pair<int,double> > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++)
00550       for(vector<pair<int,double> >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++)
00551         res[(*iter3).first]+=(*iter3).second;
00552     set<int> procsSet;
00553     int id=-1;
00554     const vector<std::pair<int,int> >& mapping=_mapping.getSendingIds();
00555     for(vector<std::pair<int,int> >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++)
00556       {
00557         std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first);
00558         if(isIns.second)
00559           id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin();
00560         resPerProcI[id].push_back((*iter2).second);
00561         resPerProcD[id].push_back(res[iter2-mapping.begin()]);
00562       }
00563   }
00564 
00569   void InterpolationMatrix::findAdditionnalElements(ElementLocator& elementLocator, std::vector<std::vector<int> >& elementsToAdd,
00570                                                     const std::vector<std::vector<int> >& resPerProcI, const std::vector<std::vector<int> >& globalIdsPartial)
00571   {
00572     std::set<int> globalIds;
00573     int nbLazyProcs=globalIdsPartial.size();
00574     for(int i=0;i<nbLazyProcs;i++)
00575       globalIds.insert(globalIdsPartial[i].begin(),globalIdsPartial[i].end());
00576     std::vector<int> tmp(globalIds.size());
00577     std::copy(globalIds.begin(),globalIds.end(),tmp.begin());
00578     globalIds.clear();
00579     elementLocator.sendCandidatesForAddElementsW(tmp);
00580     elementLocator.recvAddElementsFromLazyProcsW(elementsToAdd);
00581   }
00582 
00583   void InterpolationMatrix::addGhostElements(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& elementsToAdd)
00584   {
00585     std::vector< std::vector< std::map<int,double> > > data1;
00586     std::vector<int> data2;
00587     serializeMe(data1,data2);
00588     initialize();
00589     int nbOfDistProcs=distantProcs.size();
00590     for(int i=0;i<nbOfDistProcs;i++)
00591       {
00592         int procId=distantProcs[i];
00593         const std::vector<int>& eltsForThisProc=elementsToAdd[i];
00594         if(!eltsForThisProc.empty())
00595           {
00596             std::vector<int>::iterator iter1=std::find(data2.begin(),data2.end(),procId);
00597             std::map<int,double> *toFeed=0;
00598             if(iter1!=data2.end())
00599               {//to test
00600                 int rank=iter1-data2.begin();
00601                 toFeed=&(data1[rank].back());
00602               }
00603             else
00604               {
00605                 iter1=std::lower_bound(data2.begin(),data2.end(),procId);
00606                 int rank=iter1-data2.begin();
00607                 data2.insert(iter1,procId);
00608                 std::vector< std::map<int,double> > tmp(data1.front().size());
00609                 data1.insert(data1.begin()+rank,tmp);
00610                 toFeed=&(data1[rank].back());
00611               }
00612             for(std::vector<int>::const_iterator iter2=eltsForThisProc.begin();iter2!=eltsForThisProc.end();iter2++)
00613               (*toFeed)[*iter2]=0.;
00614           }
00615       }
00616     //
00617     nbOfDistProcs=data2.size();
00618     for(int j=0;j<nbOfDistProcs;j++)
00619       fillDSFromVM(data2[j],0,data1[j],0);
00620   }
00621 
00622   int InterpolationMatrix::mergePolicies(const std::vector<int>& policyPartial)
00623   {
00624     if(policyPartial.empty())
00625       return ElementLocator::NO_POST_TREATMENT_POLICY;
00626     int ref=policyPartial[0];
00627      std::vector<int>::const_iterator iter1=std::find_if(policyPartial.begin(),policyPartial.end(),std::bind2nd(std::not_equal_to<int>(),ref));
00628     if(iter1!=policyPartial.end())
00629       {
00630         std::ostringstream msg; msg << "Incompatible policies between lazy procs each other : proc # " << iter1-policyPartial.begin();
00631         throw INTERP_KERNEL::Exception(msg.str().c_str());
00632       }
00633     return ref;
00634   }
00635 
00645   void InterpolationMatrix::mergeRowSum(const std::vector< std::vector<double> >& rowsPartialSumD, const std::vector< std::vector<int> >& globalIdsPartial,
00646                                         std::vector<int>& globalIdsLazySideInteraction, std::vector<double>& sumCorresponding)
00647   {
00648     std::map<int,double> sumToReturn;
00649     int nbLazyProcs=rowsPartialSumD.size();
00650     for(int i=0;i<nbLazyProcs;i++)
00651       {
00652         const std::vector<double>& rowSumOfP=rowsPartialSumD[i];
00653         const std::vector<int>& globalIdsOfP=globalIdsPartial[i];
00654         std::vector<double>::const_iterator iter1=rowSumOfP.begin();
00655         std::vector<int>::const_iterator iter2=globalIdsOfP.begin();
00656         for(;iter1!=rowSumOfP.end();iter1++,iter2++)
00657           sumToReturn[*iter2]+=*iter1;
00658       }
00659     //
00660     int lgth=sumToReturn.size();
00661     globalIdsLazySideInteraction.resize(lgth);
00662     sumCorresponding.resize(lgth);
00663     std::vector<int>::iterator iter3=globalIdsLazySideInteraction.begin();
00664     std::vector<double>::iterator iter4=sumCorresponding.begin();
00665     for(std::map<int,double>::const_iterator iter5=sumToReturn.begin();iter5!=sumToReturn.end();iter5++,iter3++,iter4++)
00666       {
00667         *iter3=(*iter5).first;
00668         *iter4=(*iter5).second;
00669       }
00670   }
00671 
00680   void InterpolationMatrix::mergeRowSum2(const std::vector< std::vector<int> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD,
00681                                          const std::vector<int>& globalIdsLazySideInteraction, const std::vector<double>& sumCorresponding)
00682   {
00683     std::map<int,double> acc;
00684     std::vector<int>::const_iterator iter1=globalIdsLazySideInteraction.begin();
00685     std::vector<double>::const_iterator iter2=sumCorresponding.begin();
00686     for(;iter1!=globalIdsLazySideInteraction.end();iter1++,iter2++)
00687       acc[*iter1]=*iter2;
00688     //
00689     int nbLazyProcs=globalIdsPartial.size();
00690     for(int i=0;i<nbLazyProcs;i++)
00691       {
00692         const std::vector<int>& tmp1=globalIdsPartial[i];
00693         std::vector<double>& tmp2=rowsPartialSumD[i];
00694         std::vector<int>::const_iterator iter3=tmp1.begin();
00695         std::vector<double>::iterator iter4=tmp2.begin();
00696         for(;iter3!=tmp1.end();iter3++,iter4++)
00697           *iter4=acc[*iter3];
00698       }
00699   }
00700   
00701   void InterpolationMatrix::mergeRowSum3(const std::vector< std::vector<int> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD)
00702   {
00703     std::map<int,double> sum;
00704     std::vector< std::vector<int> >::const_iterator iter1=globalIdsPartial.begin();
00705     std::vector< std::vector<double> >::iterator iter2=rowsPartialSumD.begin();
00706     for(;iter1!=globalIdsPartial.end();iter1++,iter2++)
00707       {
00708         std::vector<int>::const_iterator iter3=(*iter1).begin();
00709         std::vector<double>::const_iterator iter4=(*iter2).begin();
00710         for(;iter3!=(*iter1).end();iter3++,iter4++)
00711           sum[*iter3]+=*iter4;
00712       }
00713     iter2=rowsPartialSumD.begin();
00714     for(iter1=globalIdsPartial.begin();iter1!=globalIdsPartial.end();iter1++,iter2++)
00715       {
00716         std::vector<int>::const_iterator iter3=(*iter1).begin();
00717         std::vector<double>::iterator iter4=(*iter2).begin();
00718         for(;iter3!=(*iter1).end();iter3++,iter4++)
00719           *iter4=sum[*iter3];
00720       }
00721   }
00722 
00730   void InterpolationMatrix::mergeCoeffs(const std::vector<int>& procsInInteraction, const std::vector< std::vector<int> >& rowsPartialSumI,
00731                                         const std::vector<std::vector<int> >& globalIdsPartial, std::vector<std::vector<double> >& denoStrorageInv)
00732   {
00733     //preparing fast access structures
00734     std::map<int,int> procT;
00735     int localProcId=0;
00736     for(std::vector<int>::const_iterator iter1=procsInInteraction.begin();iter1!=procsInInteraction.end();iter1++,localProcId++)
00737       procT[*iter1]=localProcId;
00738     int size=procsInInteraction.size();
00739     std::vector<std::map<int,int> > localToGlobal(size);
00740     for(int i=0;i<size;i++)
00741       {
00742         std::map<int,int>& myLocalToGlobal=localToGlobal[i];
00743         const std::vector<int>& locals=rowsPartialSumI[i];
00744         const std::vector<int>& globals=globalIdsPartial[i];
00745         std::vector<int>::const_iterator iter3=locals.begin();
00746         std::vector<int>::const_iterator iter4=globals.begin();
00747         for(;iter3!=locals.end();iter3++,iter4++)
00748           myLocalToGlobal[*iter3]=*iter4;
00749       }
00750     //
00751     const vector<std::pair<int,int> >& mapping=_mapping.getSendingIds();
00752     std::map<int,double> globalIdVal;
00753     //accumulate for same global id on lazy part.
00754     for(vector<vector<pair<int,double> > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++)
00755       for(vector<pair<int,double> >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00756         {
00757           const std::pair<int,int>& distantLocalLazyId=mapping[(*iter2).first];
00758           int localLazyProcId=procT[distantLocalLazyId.first];
00759           int globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second];
00760           globalIdVal[globalDistantLazyId]+=(*iter2).second;
00761         }
00762     //perform merge
00763     std::vector<std::vector<double> >::iterator iter3=denoStrorageInv.begin();
00764     for(vector<vector<pair<int,double> > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter3++)
00765       {
00766         double val=(*iter3).back();
00767         (*iter3).resize((*iter1).size());
00768         std::vector<double>::iterator iter4=(*iter3).begin();
00769         for(vector<pair<int,double> >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter4++)
00770           {
00771             const std::pair<int,int>& distantLocalLazyId=mapping[(*iter2).first];
00772             int localLazyProcId=procT[distantLocalLazyId.first];
00773             int globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second];
00774             double newVal=globalIdVal[globalDistantLazyId];
00775             if((*iter2).second!=0.)
00776               (*iter4)=val*newVal/(*iter2).second;
00777             else
00778               (*iter4)=std::numeric_limits<double>::max();
00779             (*iter2).second=newVal;
00780           }
00781       }
00782   }
00783 
00784   void InterpolationMatrix::divideByGlobalRowSum(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& resPerProcI,
00785                                                  const std::vector<std::vector<double> >& resPerProcD, std::vector<std::vector<double> >& deno)
00786   {
00787     map<int,double> fastSums;
00788     int procId=0;
00789     for(vector<int>::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++)
00790       {
00791         const std::vector<int>& currentProcI=resPerProcI[procId];
00792         const std::vector<double>& currentProcD=resPerProcD[procId];
00793         vector<double>::const_iterator iter3=currentProcD.begin();
00794         for(vector<int>::const_iterator iter2=currentProcI.begin();iter2!=currentProcI.end();iter2++,iter3++)
00795           fastSums[_col_offsets[std::make_pair(*iter1,*iter2)]]=*iter3;
00796       }
00797     deno.resize(_coeffs.size());
00798     vector<vector<double> >::iterator iter6=deno.begin();
00799     for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++)
00800       {
00801         (*iter6).resize((*iter4).size());
00802         vector<double>::iterator iter7=(*iter6).begin();
00803         for(vector<pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++,iter7++)
00804           *iter7=fastSums[(*iter5).first];
00805       }
00806   }
00807 
00808   void InterpolationMatrix::computeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
00809   {
00810     denoStrorage.resize(_coeffs.size());
00811     vector<vector<double> >::iterator iter2=denoStrorage.begin();
00812     for(vector<vector<pair<int,double> > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++)
00813       {
00814         (*iter2).resize((*iter1).size());
00815         double sumOfCurrentRow=0.;
00816         for(vector<pair<int,double> >::const_iterator iter3=(*iter1).begin();iter3!=(*iter1).end();iter3++)
00817           sumOfCurrentRow+=(*iter3).second;
00818         std::fill((*iter2).begin(),(*iter2).end(),sumOfCurrentRow);
00819       }
00820   }
00821 
00822   void InterpolationMatrix::resizeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
00823   {
00824     vector<vector<double> >::iterator iter2=denoStrorage.begin();
00825     for(vector<vector<pair<int,double> > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++)
00826       {
00827         double val=(*iter2).back();
00828         (*iter2).resize((*iter1).size());
00829         std::fill((*iter2).begin(),(*iter2).end(),val);
00830       }
00831   }
00832 
00833   // ==================================================================
00834   // The call to this method updates the arrays on the target side
00835   //   so that they know which amount of data from which processor they 
00836   //   should expect. 
00837   //   That call makes actual interpolations via multiply method 
00838   //   available.
00839   // ==================================================================
00840 
00841   void InterpolationMatrix::prepare()
00842   {
00843     int nbelems = _source_field->getField()->getNumberOfTuples();
00844     for (int ielem=0; ielem < nbelems; ielem++)
00845       {
00846         _row_offsets[ielem+1]+=_row_offsets[ielem];
00847       }
00848     _mapping.prepareSendRecv();
00849   }
00850 
00851 
00852   //   =======================================================================
00853   //   brief performs t=Ws, where t is the target field, s is the source field
00854 
00855   //   The call to this method must be called both on the working side 
00856   //   and on the idle side. On the working side, the vector  T=VT^(-1).(W.S)
00857   //   is computed and sent. On the idle side, no computation is done, but the 
00858   //   result from the working side is received and the field is updated.
00859 
00860   //   param field source field on processors involved on the source side,
00861   //   target field on processors on the target side
00862   //   =======================================================================
00863 
00864   void InterpolationMatrix::multiply(MEDCouplingFieldDouble& field) const
00865   {
00866     int nbcomp = field.getArray()->getNumberOfComponents();
00867     vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
00868 
00869     //computing the matrix multiply on source side
00870     if (_source_group.containsMyRank())
00871       {
00872         int nbrows = _coeffs.size();
00873 
00874         // performing W.S
00875         // W is the intersection matrix
00876         // S is the source vector
00877 
00878         for (int irow=0; irow<nbrows; irow++)
00879           {
00880             for (int icomp=0; icomp< nbcomp; icomp++)
00881               {
00882                 double coeff_row = field.getIJ(irow,icomp);
00883                 for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
00884                   {
00885                     int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
00886                     double value = _coeffs[irow][icol-_row_offsets[irow]].second;
00887                     double deno = _deno_multiply[irow][icol-_row_offsets[irow]];
00888                     target_value[colid*nbcomp+icomp]+=value*coeff_row/deno;
00889                   }
00890               }
00891           }
00892       }
00893 
00894     if (_target_group.containsMyRank())
00895       {
00896         int nbelems = field.getArray()->getNumberOfTuples() ;
00897         double* value = const_cast<double*> (field.getArray()->getPointer());
00898         for (int i=0; i<nbelems*nbcomp; i++)
00899           {
00900             value[i]=0.0;
00901           }
00902       }
00903 
00904     //on source side : sending  T=VT^(-1).(W.S)
00905     //on target side :: receiving T and storing it in field
00906     _mapping.sendRecv(&target_value[0],field);
00907   }
00908   
00909 
00910   // =========================================================================
00911   // brief performs s=WTt, where t is the target field, s is the source field,
00912   // WT is the transpose matrix from W
00913 
00914   //   The call to this method must be called both on the working side 
00915   //   and on the idle side. On the working side, the target vector T is
00916   //   received and the vector  S=VS^(-1).(WT.T) is computed to update
00917   //   the field. 
00918   //   On the idle side, no computation is done, but the field is sent.
00919 
00920   //   param field source field on processors involved on the source side,
00921   //   target field on processors on the target side
00922   // =========================================================================
00923 
00924   void InterpolationMatrix::transposeMultiply(MEDCouplingFieldDouble& field) const
00925   {
00926     int nbcomp = field.getArray()->getNumberOfComponents();
00927     vector<double> source_value(_col_offsets.size()* nbcomp,0.0);
00928     _mapping.reverseSendRecv(&source_value[0],field);
00929 
00930     //treatment of the transpose matrix multiply on the source side
00931     if (_source_group.containsMyRank())
00932       {
00933         int nbrows    = _coeffs.size();
00934         double *array = field.getArray()->getPointer() ;
00935 
00936         // Initialization
00937         std::fill(array, array+nbrows*nbcomp, 0.0) ;
00938 
00939         //performing WT.T
00940         //WT is W transpose
00941         //T is the target vector
00942         for (int irow = 0; irow < nbrows; irow++)
00943           {
00944             for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
00945               {
00946                 int colid    = _coeffs[irow][icol-_row_offsets[irow]].first;
00947                 double value = _coeffs[irow][icol-_row_offsets[irow]].second;
00948                 double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
00949                 for (int icomp=0; icomp<nbcomp; icomp++)
00950                   {
00951                     double coeff_row = source_value[colid*nbcomp+icomp];
00952                     array[irow*nbcomp+icomp] += value*coeff_row/deno;
00953                   }
00954               }
00955           }
00956       }
00957   }
00958 
00959   bool InterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const
00960   {
00961     return method=="P0";
00962   }
00963 }