+ }
+
+ void InterpolationMatrix::serializeMe(std::vector< std::vector< std::map<int,double> > >& data1, std::vector<int>& data2) const
+ {
+ data1.clear();
+ data2.clear();
+ const std::vector<std::pair<int,int> >& sendingIds=_mapping.getSendingIds();
+ std::set<int> procsS;
+ for(std::vector<std::pair<int,int> >::const_iterator iter1=sendingIds.begin();iter1!=sendingIds.end();iter1++)
+ procsS.insert((*iter1).first);
+ data1.resize(procsS.size());
+ data2.resize(procsS.size());
+ std::copy(procsS.begin(),procsS.end(),data2.begin());
+ std::map<int,int> fastProcAcc;
+ int id=0;
+ for(std::set<int>::const_iterator iter2=procsS.begin();iter2!=procsS.end();iter2++,id++)
+ fastProcAcc[*iter2]=id;
+ int nbOfSrcElt=_coeffs.size();
+ for(std::vector< std::vector< std::map<int,double> > >::iterator iter3=data1.begin();iter3!=data1.end();iter3++)
+ (*iter3).resize(nbOfSrcElt);
+ id=0;
+ for(std::vector< std::vector< std::pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,id++)
+ {
+ for(std::vector< std::pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++)
+ {
+ const std::pair<int,int>& elt=sendingIds[(*iter5).first];
+ data1[fastProcAcc[elt.first]][id][elt.second]=(*iter5).second;
+ }
+ }
+ }
+
+ void InterpolationMatrix::initialize()
+ {
+ int lgth=_coeffs.size();
+ _row_offsets.clear(); _row_offsets.resize(lgth+1);
+ _coeffs.clear(); _coeffs.resize(lgth);
+ _target_volume.clear(); _target_volume.resize(lgth);
+ _col_offsets.clear();
+ _mapping.initialize();
+ }
+
+ void InterpolationMatrix::finishContributionW(ElementLocator& elementLocator)
+ {
+ NatureOfField nature=elementLocator.getLocalNature();
+ switch(nature)
+ {
+ case ConservativeVolumic:
+ computeConservVolDenoW(elementLocator);
+ break;
+ case Integral:
+ {
+ if(!elementLocator.isM1DCorr())
+ computeIntegralDenoW(elementLocator);
+ else
+ computeGlobConstraintDenoW(elementLocator);
+ break;
+ }
+ case IntegralGlobConstraint:
+ computeGlobConstraintDenoW(elementLocator);
+ break;
+ case RevIntegral:
+ {
+ if(!elementLocator.isM1DCorr())
+ computeRevIntegralDenoW(elementLocator);
+ else
+ computeConservVolDenoW(elementLocator);
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
+ break;
+ }
+ }
+
+ void InterpolationMatrix::finishContributionL(ElementLocator& elementLocator)
+ {
+ NatureOfField nature=elementLocator.getLocalNature();
+ switch(nature)
+ {
+ case ConservativeVolumic:
+ computeConservVolDenoL(elementLocator);
+ break;
+ case Integral:
+ {
+ if(!elementLocator.isM1DCorr())
+ computeIntegralDenoL(elementLocator);
+ else
+ computeConservVolDenoL(elementLocator);
+ break;
+ }
+ case IntegralGlobConstraint:
+ //this is not a bug doing like ConservativeVolumic
+ computeConservVolDenoL(elementLocator);
+ break;
+ case RevIntegral:
+ {
+ if(!elementLocator.isM1DCorr())
+ computeRevIntegralDenoL(elementLocator);
+ else
+ computeConservVolDenoL(elementLocator);
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
+ break;
+ }
+ }
+
+ void InterpolationMatrix::computeConservVolDenoW(ElementLocator& elementLocator)
+ {
+ computeGlobalColSum(_deno_reverse_multiply);
+ computeGlobalRowSum(elementLocator,_deno_multiply,_deno_reverse_multiply);
+ }
+
+ void InterpolationMatrix::computeConservVolDenoL(ElementLocator& elementLocator)
+ {
+ int pol1=elementLocator.sendPolicyToWorkingSideL();
+ if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
+ {
+ elementLocator.recvFromWorkingSideL();
+ elementLocator.sendToWorkingSideL();
+ }
+ else if(ElementLocator::CUMULATIVE_POLICY)
+ {
+ //ask for lazy side to deduce ids eventually missing on working side and to send it back.
+ elementLocator.recvLocalIdsFromWorkingSideL();
+ elementLocator.sendCandidatesGlobalIdsToWorkingSideL();
+ elementLocator.recvCandidatesForAddElementsL();
+ elementLocator.sendAddElementsToWorkingSideL();
+ //Working side has updated its eventually missing ids updates its global ids with lazy side procs contribution
+ elementLocator.recvLocalIdsFromWorkingSideL();
+ elementLocator.sendGlobalIdsToWorkingSideL();
+ //like no post treatment
+ elementLocator.recvFromWorkingSideL();
+ elementLocator.sendToWorkingSideL();
+ }
+ else
+ throw INTERP_KERNEL::Exception("Not managed policy detected on lazy side : not implemented !");
+ }
+
+ void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator)
+ {
+ MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
+ _deno_multiply.resize(_coeffs.size());
+ vector<vector<double> >::iterator iter6=_deno_multiply.begin();
+ const double *values=source_triangle_surf->getArray()->getConstPointer();
+ for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
+ {
+ (*iter6).resize((*iter4).size());
+ std::fill((*iter6).begin(),(*iter6).end(),*values);
+ }
+ source_triangle_surf->decrRef();
+ _deno_reverse_multiply=_target_volume;
+ }
+
+ void InterpolationMatrix::computeRevIntegralDenoW(ElementLocator& elementLocator)
+ {
+ _deno_multiply=_target_volume;
+ MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
+ _deno_reverse_multiply.resize(_coeffs.size());
+ vector<vector<double> >::iterator iter6=_deno_reverse_multiply.begin();
+ const double *values=source_triangle_surf->getArray()->getConstPointer();
+ for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
+ {
+ (*iter6).resize((*iter4).size());
+ std::fill((*iter6).begin(),(*iter6).end(),*values);
+ }
+ source_triangle_surf->decrRef();
+ }
+
+ /*!
+ * Nothing to do because surface computation is on working side.
+ */
+ void InterpolationMatrix::computeIntegralDenoL(ElementLocator& elementLocator)
+ {
+ }
+
+ /*!
+ * Nothing to do because surface computation is on working side.
+ */
+ void InterpolationMatrix::computeRevIntegralDenoL(ElementLocator& elementLocator)
+ {
+ }
+
+
+ void InterpolationMatrix::computeGlobConstraintDenoW(ElementLocator& elementLocator)
+ {
+ computeGlobalColSum(_deno_multiply);
+ computeGlobalRowSum(elementLocator,_deno_reverse_multiply,_deno_multiply);
+ }
+
+ void InterpolationMatrix::computeGlobalRowSum(ElementLocator& elementLocator, std::vector<std::vector<double> >& denoStrorage, std::vector<std::vector<double> >& denoStrorageInv)
+ {
+ //stores id in distant procs sorted by lazy procs connected with
+ vector< vector<int> > rowsPartialSumI;
+ //stores for each lazy procs connected with, if global info is available and if it's the case the policy
+ vector<int> policyPartial;
+ //stores the corresponding values.
+ vector< vector<double> > rowsPartialSumD;
+ elementLocator.recvPolicyFromLazySideW(policyPartial);
+ int pol1=mergePolicies(policyPartial);
+ if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
+ {
+ computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
+ elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
+ elementLocator.recvSumFromLazySideW(rowsPartialSumD);
+ }
+ else if(pol1==ElementLocator::CUMULATIVE_POLICY)
+ {
+ //updateWithNewAdditionnalElements(addingElements);
+ //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,
+ //if policyPartial has CUMALATIVE_POLICY in each.
+ vector< vector<int> > globalIdsPartial;
+ computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
+ elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
+ elementLocator.recvCandidatesGlobalIdsFromLazyProcsW(globalIdsPartial);
+ std::vector< std::vector<int> > addingElements;
+ findAdditionnalElements(elementLocator,addingElements,rowsPartialSumI,globalIdsPartial);
+ addGhostElements(elementLocator.getDistantProcIds(),addingElements);
+ rowsPartialSumI.clear();
+ globalIdsPartial.clear();
+ computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
+ elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
+ elementLocator.recvGlobalIdsFromLazyProcsW(rowsPartialSumI,globalIdsPartial);
+ //
+ elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
+ elementLocator.recvSumFromLazySideW(rowsPartialSumD);
+ mergeRowSum3(globalIdsPartial,rowsPartialSumD);
+ mergeCoeffs(elementLocator.getDistantProcIds(),rowsPartialSumI,globalIdsPartial,denoStrorageInv);
+ }
+ else
+ throw INTERP_KERNEL::Exception("Not managed policy detected : not implemented !");
+ divideByGlobalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD,denoStrorage);
+ }
+
+ /*!
+ * @param distantProcs in parameter that indicates which lazy procs are concerned.
+ * @param resPerProcI out parameter that must be cleared before calling this method. The size of 1st dimension is equal to the size of 'distantProcs'.
+ * It contains the element ids (2nd dimension) of the corresponding lazy proc.
+ * @param resPerProcD out parameter with the same format than 'resPerProcI'. It contains corresponding sum values.
+ */
+ void InterpolationMatrix::computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<int> >& resPerProcI,
+ std::vector<std::vector<double> >& resPerProcD) const
+ {
+ resPerProcI.resize(distantProcs.size());
+ resPerProcD.resize(distantProcs.size());
+ std::vector<double> res(_col_offsets.size());
+ for(vector<vector<pair<int,double> > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++)
+ for(vector<pair<int,double> >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++)
+ res[(*iter3).first]+=(*iter3).second;
+ set<int> procsSet;
+ int id=-1;
+ const vector<std::pair<int,int> >& mapping=_mapping.getSendingIds();
+ for(vector<std::pair<int,int> >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++)
+ {
+ std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first);
+ if(isIns.second)
+ id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin();
+ resPerProcI[id].push_back((*iter2).second);
+ resPerProcD[id].push_back(res[iter2-mapping.begin()]);
+ }
+ }
+
+ /*!
+ * This method is only usable when CUMULATIVE_POLICY detected. This method finds elements ids (typically nodes) lazy side that
+ * are not present in columns of 'this' and that should regarding cumulative merge of elements regarding their global ids.
+ */
+ void InterpolationMatrix::findAdditionnalElements(ElementLocator& elementLocator, std::vector<std::vector<int> >& elementsToAdd,
+ const std::vector<std::vector<int> >& resPerProcI, const std::vector<std::vector<int> >& globalIdsPartial)
+ {
+ std::set<int> globalIds;
+ int nbLazyProcs=globalIdsPartial.size();
+ for(int i=0;i<nbLazyProcs;i++)
+ globalIds.insert(globalIdsPartial[i].begin(),globalIdsPartial[i].end());
+ std::vector<int> tmp(globalIds.size());
+ std::copy(globalIds.begin(),globalIds.end(),tmp.begin());
+ globalIds.clear();
+ elementLocator.sendCandidatesForAddElementsW(tmp);
+ elementLocator.recvAddElementsFromLazyProcsW(elementsToAdd);
+ }
+
+ void InterpolationMatrix::addGhostElements(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& elementsToAdd)
+ {
+ std::vector< std::vector< std::map<int,double> > > data1;
+ std::vector<int> data2;
+ serializeMe(data1,data2);
+ initialize();
+ int nbOfDistProcs=distantProcs.size();
+ for(int i=0;i<nbOfDistProcs;i++)
+ {
+ int procId=distantProcs[i];
+ const std::vector<int>& eltsForThisProc=elementsToAdd[i];
+ if(!eltsForThisProc.empty())
+ {
+ std::vector<int>::iterator iter1=std::find(data2.begin(),data2.end(),procId);
+ std::map<int,double> *toFeed=0;
+ if(iter1!=data2.end())
+ {//to test
+ int rank=iter1-data2.begin();
+ toFeed=&(data1[rank].back());
+ }
+ else
+ {
+ iter1=std::lower_bound(data2.begin(),data2.end(),procId);
+ int rank=iter1-data2.begin();
+ data2.insert(iter1,procId);
+ std::vector< std::map<int,double> > tmp(data1.front().size());
+ data1.insert(data1.begin()+rank,tmp);
+ toFeed=&(data1[rank].back());
+ }
+ for(std::vector<int>::const_iterator iter2=eltsForThisProc.begin();iter2!=eltsForThisProc.end();iter2++)
+ (*toFeed)[*iter2]=0.;
+ }
+ }
+ //
+ nbOfDistProcs=data2.size();
+ for(int j=0;j<nbOfDistProcs;j++)
+ fillDSFromVM(data2[j],0,data1[j],0);
+ }
+
+ int InterpolationMatrix::mergePolicies(const std::vector<int>& policyPartial)
+ {
+ if(policyPartial.empty())
+ return ElementLocator::NO_POST_TREATMENT_POLICY;
+ int ref=policyPartial[0];
+ std::vector<int>::const_iterator iter1=std::find_if(policyPartial.begin(),policyPartial.end(),std::bind2nd(std::not_equal_to<int>(),ref));
+ if(iter1!=policyPartial.end())
+ {
+ std::ostringstream msg; msg << "Incompatible policies between lazy procs each other : proc # " << iter1-policyPartial.begin();
+ throw INTERP_KERNEL::Exception(msg.str().c_str());
+ }
+ return ref;
+ }
+
+ /*!
+ * This method introduce global ids aspects in computed 'rowsPartialSumD'.
+ * As precondition rowsPartialSumD.size()==policyPartial.size()==globalIdsPartial.size(). Foreach i in [0;rowsPartialSumD.size() ) rowsPartialSumD[i].size()==globalIdsPartial[i].size()
+ * @param rowsPartialSumD : in parameter, Partial row sum computed for each lazy procs connected with.
+ * @param rowsPartialSumI : in parameter, Corresponding local ids for each lazy procs connected with.
+ * @param globalIdsPartial : in parameter, the global numbering, of elements connected with.
+ * @param globalIdsLazySideInteraction : out parameter, constituted from all global ids of lazy procs connected with.
+ * @para sumCorresponding : out parameter, relative to 'globalIdsLazySideInteraction'
+ */
+ void InterpolationMatrix::mergeRowSum(const std::vector< std::vector<double> >& rowsPartialSumD, const std::vector< std::vector<int> >& globalIdsPartial,
+ std::vector<int>& globalIdsLazySideInteraction, std::vector<double>& sumCorresponding)
+ {
+ std::map<int,double> sumToReturn;
+ int nbLazyProcs=rowsPartialSumD.size();
+ for(int i=0;i<nbLazyProcs;i++)
+ {
+ const std::vector<double>& rowSumOfP=rowsPartialSumD[i];
+ const std::vector<int>& globalIdsOfP=globalIdsPartial[i];
+ std::vector<double>::const_iterator iter1=rowSumOfP.begin();
+ std::vector<int>::const_iterator iter2=globalIdsOfP.begin();
+ for(;iter1!=rowSumOfP.end();iter1++,iter2++)
+ sumToReturn[*iter2]+=*iter1;
+ }
+ //
+ int lgth=sumToReturn.size();
+ globalIdsLazySideInteraction.resize(lgth);
+ sumCorresponding.resize(lgth);
+ std::vector<int>::iterator iter3=globalIdsLazySideInteraction.begin();
+ std::vector<double>::iterator iter4=sumCorresponding.begin();
+ for(std::map<int,double>::const_iterator iter5=sumToReturn.begin();iter5!=sumToReturn.end();iter5++,iter3++,iter4++)
+ {
+ *iter3=(*iter5).first;
+ *iter4=(*iter5).second;
+ }
+ }
+
+ /*!
+ * This method simply reorganize the result contained in 'sumCorresponding' computed by lazy side into 'rowsPartialSumD' with help of 'globalIdsPartial' and 'globalIdsLazySideInteraction'
+ *
+ * @param globalIdsPartial : in parameter, global ids sorted by lazy procs
+ * @param rowsPartialSumD : in/out parameter, with exactly the same size as 'globalIdsPartial'
+ * @param globalIdsLazySideInteraction : in parameter that represents ALL the global ids of every lazy procs in interaction
+ * @param sumCorresponding : in parameter with same size as 'globalIdsLazySideInteraction' that stores the corresponding sum of 'globalIdsLazySideInteraction'
+ */
+ void InterpolationMatrix::mergeRowSum2(const std::vector< std::vector<int> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD,
+ const std::vector<int>& globalIdsLazySideInteraction, const std::vector<double>& sumCorresponding)
+ {
+ std::map<int,double> acc;
+ std::vector<int>::const_iterator iter1=globalIdsLazySideInteraction.begin();
+ std::vector<double>::const_iterator iter2=sumCorresponding.begin();
+ for(;iter1!=globalIdsLazySideInteraction.end();iter1++,iter2++)
+ acc[*iter1]=*iter2;
+ //
+ int nbLazyProcs=globalIdsPartial.size();
+ for(int i=0;i<nbLazyProcs;i++)
+ {
+ const std::vector<int>& tmp1=globalIdsPartial[i];
+ std::vector<double>& tmp2=rowsPartialSumD[i];
+ std::vector<int>::const_iterator iter3=tmp1.begin();
+ std::vector<double>::iterator iter4=tmp2.begin();
+ for(;iter3!=tmp1.end();iter3++,iter4++)
+ *iter4=acc[*iter3];
+ }