1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
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"
48 Creates an empty matrix structure linking two distributed supports.
49 The method must be called by all processors belonging to source
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
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)
69 int nbelems = source_field->getField()->getNumberOfTuples();
70 _row_offsets.resize(nbelems+1);
71 _coeffs.resize(nbelems);
72 _target_volume.resize(nbelems);
75 InterpolationMatrix::~InterpolationMatrix()
81 \brief Adds the contribution of a distant subdomain to the*
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.
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
95 void InterpolationMatrix::addContribution ( MEDCouplingPointSet& distant_support,
97 const int* distant_elems,
98 const std::string& srcMeth,
99 const std::string& targetMeth)
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 )
110 if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==2)
112 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(source_supportC);
113 INTERP_KERNEL::Interpolation2D interpolation(*this);
114 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
116 else if(source_supportC->getMeshDimension()==3 && source_supportC->getSpaceDimension()==3)
118 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(source_supportC);
119 INTERP_KERNEL::Interpolation3D interpolation(*this);
120 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
122 else if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==3)
124 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(source_supportC);
125 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
126 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
129 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
131 else if ( source_supportC->getMeshDimension() == -1 )
133 if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==2)
135 MEDCouplingNormalizedUnstructuredMesh<2,2> distant_mesh_wrapper(distant_supportC);
136 INTERP_KERNEL::Interpolation2D interpolation(*this);
137 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
139 else if(distant_supportC->getMeshDimension()==3 && distant_supportC->getSpaceDimension()==3)
141 MEDCouplingNormalizedUnstructuredMesh<3,3> distant_mesh_wrapper(distant_supportC);
142 INTERP_KERNEL::Interpolation3D interpolation(*this);
143 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
145 else if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==3)
147 MEDCouplingNormalizedUnstructuredMesh<3,2> distant_mesh_wrapper(distant_supportC);
148 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
149 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
152 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
154 else if ( distant_support.getMeshDimension() == 2
155 && _source_support->getMeshDimension() == 3
156 && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
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();
165 else if ( distant_support.getMeshDimension() == 1
166 && _source_support->getMeshDimension() == 2
167 && distant_support.getSpaceDimension() == 2 && _source_support->getSpaceDimension() == 2)
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();
176 else if ( distant_support.getMeshDimension() == 3
177 && _source_support->getMeshDimension() == 1
178 && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
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();
187 else if (distant_support.getMeshDimension() != _source_support->getMeshDimension())
189 throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
191 else if( distant_support.getMeshDimension() == 1
192 && distant_support.getSpaceDimension() == 1 )
194 MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(distant_supportC);
195 MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(source_supportC);
197 INTERP_KERNEL::Interpolation1D interpolation(*this);
198 interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
199 target_wrapper.releaseTempArrays();
200 source_wrapper.releaseTempArrays();
202 else if( distant_support.getMeshDimension() == 1
203 && distant_support.getSpaceDimension() == 2 )
205 MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(distant_supportC);
206 MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(source_supportC);
208 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
209 interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
210 target_wrapper.releaseTempArrays();
211 source_wrapper.releaseTempArrays();
213 else if ( distant_support.getMeshDimension() == 2
214 && distant_support.getSpaceDimension() == 3 )
216 MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC);
217 MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC);
219 INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
220 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
221 target_wrapper.releaseTempArrays();
222 source_wrapper.releaseTempArrays();
224 else if ( distant_support.getMeshDimension() == 2
225 && distant_support.getSpaceDimension() == 2)
227 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
228 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
230 INTERP_KERNEL::Interpolation2D interpolator (*this);
231 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
232 target_wrapper.releaseTempArrays();
233 source_wrapper.releaseTempArrays();
235 else if ( distant_support.getMeshDimension() == 3
236 && distant_support.getSpaceDimension() == 3 )
238 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
239 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
241 INTERP_KERNEL::Interpolation3D interpolator (*this);
242 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
243 target_wrapper.releaseTempArrays();
244 source_wrapper.releaseTempArrays();
248 throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
250 bool needTargetSurf=isSurfaceComputationNeeded(targetMeth);
252 MEDCouplingFieldDouble *target_triangle_surf=0;
254 target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus());
255 fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf);
258 target_triangle_surf->decrRef();
261 void InterpolationMatrix::fillDSFromVM(int iproc_distant, const int* distant_elems, const std::vector< std::map<int,double> >& values, MEDCouplingFieldDouble *surf)
263 //loop over the elements to build the interpolation
265 int source_size=values.size();
266 for (int ielem=0; ielem < source_size; ielem++)
268 _row_offsets[ielem+1] += values[ielem].size();
269 for(map<int,double>::const_iterator iter=values[ielem].begin();iter!=values[ielem].end();iter++)
273 localId=distant_elems[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));
280 if (iter2 == _col_offsets.end())
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);
290 col_id = iter2->second;
292 //the non zero coefficient is stored
294 //col_id is the number of the column
295 //iter->second is the value of the coefficient
298 double surface = surf->getIJ(iter->first,0);
299 _target_volume[ielem].push_back(surface);
301 _coeffs[ielem].push_back(make_pair(col_id,iter->second));
306 void InterpolationMatrix::serializeMe(std::vector< std::vector< std::map<int,double> > >& data1, std::vector<int>& data2) const
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;
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);
325 for(std::vector< std::vector< std::pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,id++)
327 for(std::vector< std::pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++)
329 const std::pair<int,int>& elt=sendingIds[(*iter5).first];
330 data1[fastProcAcc[elt.first]][id][elt.second]=(*iter5).second;
335 void InterpolationMatrix::initialize()
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();
345 void InterpolationMatrix::finishContributionW(ElementLocator& elementLocator)
347 NatureOfField nature=elementLocator.getLocalNature();
350 case IntensiveMaximum:
351 computeConservVolDenoW(elementLocator);
353 case ExtensiveMaximum:
355 if(!elementLocator.isM1DCorr())
356 computeIntegralDenoW(elementLocator);
358 computeGlobConstraintDenoW(elementLocator);
361 case ExtensiveConservation:
362 computeGlobConstraintDenoW(elementLocator);
364 case IntensiveConservation:
366 if(!elementLocator.isM1DCorr())
367 computeRevIntegralDenoW(elementLocator);
369 computeConservVolDenoW(elementLocator);
373 throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
378 void InterpolationMatrix::finishContributionL(ElementLocator& elementLocator)
380 NatureOfField nature=elementLocator.getLocalNature();
383 case IntensiveMaximum:
384 computeConservVolDenoL(elementLocator);
386 case ExtensiveMaximum:
388 if(!elementLocator.isM1DCorr())
389 computeIntegralDenoL(elementLocator);
391 computeConservVolDenoL(elementLocator);
394 case ExtensiveConservation:
395 //this is not a bug doing like IntensiveMaximum
396 computeConservVolDenoL(elementLocator);
398 case IntensiveConservation:
400 if(!elementLocator.isM1DCorr())
401 computeRevIntegralDenoL(elementLocator);
403 computeConservVolDenoL(elementLocator);
407 throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
412 void InterpolationMatrix::computeConservVolDenoW(ElementLocator& elementLocator)
414 computeGlobalColSum(_deno_reverse_multiply);
415 computeGlobalRowSum(elementLocator,_deno_multiply,_deno_reverse_multiply);
418 void InterpolationMatrix::computeConservVolDenoL(ElementLocator& elementLocator)
420 int pol1=elementLocator.sendPolicyToWorkingSideL();
421 if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
423 elementLocator.recvFromWorkingSideL();
424 elementLocator.sendToWorkingSideL();
426 else if(ElementLocator::CUMULATIVE_POLICY)
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();
441 throw INTERP_KERNEL::Exception("Not managed policy detected on lazy side : not implemented !");
444 void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator)
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++)
452 (*iter6).resize((*iter4).size());
453 std::fill((*iter6).begin(),(*iter6).end(),*values);
455 source_triangle_surf->decrRef();
456 _deno_reverse_multiply=_target_volume;
459 void InterpolationMatrix::computeRevIntegralDenoW(ElementLocator& elementLocator)
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++)
468 (*iter6).resize((*iter4).size());
469 std::fill((*iter6).begin(),(*iter6).end(),*values);
471 source_triangle_surf->decrRef();
475 * Nothing to do because surface computation is on working side.
477 void InterpolationMatrix::computeIntegralDenoL(ElementLocator& elementLocator)
482 * Nothing to do because surface computation is on working side.
484 void InterpolationMatrix::computeRevIntegralDenoL(ElementLocator& elementLocator)
489 void InterpolationMatrix::computeGlobConstraintDenoW(ElementLocator& elementLocator)
491 computeGlobalColSum(_deno_multiply);
492 computeGlobalRowSum(elementLocator,_deno_reverse_multiply,_deno_multiply);
495 void InterpolationMatrix::computeGlobalRowSum(ElementLocator& elementLocator, std::vector<std::vector<double> >& denoStrorage, std::vector<std::vector<double> >& denoStrorageInv)
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)
507 computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
508 elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
509 elementLocator.recvSumFromLazySideW(rowsPartialSumD);
511 else if(pol1==ElementLocator::CUMULATIVE_POLICY)
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);
529 elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
530 elementLocator.recvSumFromLazySideW(rowsPartialSumD);
531 mergeRowSum3(globalIdsPartial,rowsPartialSumD);
532 mergeCoeffs(elementLocator.getDistantProcIds(),rowsPartialSumI,globalIdsPartial,denoStrorageInv);
535 throw INTERP_KERNEL::Exception("Not managed policy detected : not implemented !");
536 divideByGlobalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD,denoStrorage);
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.
545 void InterpolationMatrix::computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<int> >& resPerProcI,
546 std::vector<std::vector<double> >& resPerProcD) const
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;
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++)
559 std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first);
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()]);
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.
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)
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());
581 elementLocator.sendCandidatesForAddElementsW(tmp);
582 elementLocator.recvAddElementsFromLazyProcsW(elementsToAdd);
585 void InterpolationMatrix::addGhostElements(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& elementsToAdd)
587 std::vector< std::vector< std::map<int,double> > > data1;
588 std::vector<int> data2;
589 serializeMe(data1,data2);
591 int nbOfDistProcs=distantProcs.size();
592 for(int i=0;i<nbOfDistProcs;i++)
594 int procId=distantProcs[i];
595 const std::vector<int>& eltsForThisProc=elementsToAdd[i];
596 if(!eltsForThisProc.empty())
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())
602 int rank=iter1-data2.begin();
603 toFeed=&(data1[rank].back());
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());
614 for(std::vector<int>::const_iterator iter2=eltsForThisProc.begin();iter2!=eltsForThisProc.end();iter2++)
615 (*toFeed)[*iter2]=0.;
619 nbOfDistProcs=data2.size();
620 for(int j=0;j<nbOfDistProcs;j++)
621 fillDSFromVM(data2[j],0,data1[j],0);
624 int InterpolationMatrix::mergePolicies(const std::vector<int>& policyPartial)
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())
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());
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'
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)
650 std::map<int,double> sumToReturn;
651 int nbLazyProcs=rowsPartialSumD.size();
652 for(int i=0;i<nbLazyProcs;i++)
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;
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++)
669 *iter3=(*iter5).first;
670 *iter4=(*iter5).second;
675 * This method simply reorganize the result contained in 'sumCorresponding' computed by lazy side into 'rowsPartialSumD' with help of 'globalIdsPartial' and 'globalIdsLazySideInteraction'
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'
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)
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++)
691 int nbLazyProcs=globalIdsPartial.size();
692 for(int i=0;i<nbLazyProcs;i++)
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++)
703 void InterpolationMatrix::mergeRowSum3(const std::vector< std::vector<int> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD)
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++)
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++)
715 iter2=rowsPartialSumD.begin();
716 for(iter1=globalIdsPartial.begin();iter1!=globalIdsPartial.end();iter1++,iter2++)
718 std::vector<int>::const_iterator iter3=(*iter1).begin();
719 std::vector<double>::iterator iter4=(*iter2).begin();
720 for(;iter3!=(*iter1).end();iter3++,iter4++)
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
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)
735 //preparing fast access structures
736 std::map<int,int> procT;
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++)
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;
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++)
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;
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++)
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++)
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;
780 (*iter4)=std::numeric_limits<double>::max();
781 (*iter2).second=newVal;
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)
789 map<int,double> fastSums;
791 for(vector<int>::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++)
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;
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++)
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];
810 void InterpolationMatrix::computeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
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++)
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);
824 void InterpolationMatrix::resizeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
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++)
829 double val=(*iter2).back();
830 (*iter2).resize((*iter1).size());
831 std::fill((*iter2).begin(),(*iter2).end(),val);
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
839 That call makes actual interpolations via multiply method
842 void InterpolationMatrix::prepare()
844 int nbelems = _source_field->getField()->getNumberOfTuples();
845 for (int ielem=0; ielem < nbelems; ielem++)
847 _row_offsets[ielem+1]+=_row_offsets[ielem];
849 _mapping.prepareSendRecv();
855 \brief performs t=Ws, where t is the target field, s is the source field
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.
862 \param field source field on processors involved on the source side,
863 target field on processors on the target side
865 void InterpolationMatrix::multiply(MEDCouplingFieldDouble& field) const
867 int nbcomp = field.getArray()->getNumberOfComponents();
868 vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
870 //computing the matrix multiply on source side
871 if (_source_group.containsMyRank())
873 int nbrows = _coeffs.size();
876 // W is the intersection matrix
877 // S is the source vector
879 for (int irow=0; irow<nbrows; irow++)
881 for (int icomp=0; icomp< nbcomp; icomp++)
883 double coeff_row = field.getIJ(irow,icomp);
884 for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
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;
895 if (_target_group.containsMyRank())
897 int nbelems = field.getArray()->getNumberOfTuples() ;
898 double* value = const_cast<double*> (field.getArray()->getPointer());
899 for (int i=0; i<nbelems*nbcomp; i++)
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);
912 \brief performs s=WTt, where t is the target field, s is the source field,
913 WT is the transpose matrix from W
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
919 On the idle side, no computation is done, but the field is sent.
921 param field source field on processors involved on the source side,
922 target field on processors on the target side
924 void InterpolationMatrix::transposeMultiply(MEDCouplingFieldDouble& field) const
926 int nbcomp = field.getArray()->getNumberOfComponents();
927 vector<double> source_value(_col_offsets.size()* nbcomp,0.0);
928 _mapping.reverseSendRecv(&source_value[0],field);
930 //treatment of the transpose matrix multiply on the source side
931 if (_source_group.containsMyRank())
933 int nbrows = _coeffs.size();
934 double *array = field.getArray()->getPointer() ;
937 std::fill(array, array+nbrows*nbcomp, 0.0) ;
941 //T is the target vector
942 for (int irow = 0; irow < nbrows; irow++)
944 for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
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++)
951 double coeff_row = source_value[colid*nbcomp+icomp];
952 array[irow*nbcomp+icomp] += value*coeff_row/deno;
959 bool InterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const