1 // Copyright (C) 2007-2014 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 "Interpolation3D2D.txx"
33 #include "Interpolation2D1D.txx"
34 #include "MEDCouplingUMesh.hxx"
35 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
36 #include "InterpolationOptions.hxx"
37 #include "NormalizedUnstructuredMesh.hxx"
38 #include "ElementLocator.hxx"
42 // class InterpolationMatrix
43 // This class enables the storage of an interpolation matrix Wij mapping
44 // source field Sj to target field Ti via Ti=Vi^(-1).Wij.Sj.
45 // The matrix is built and stored on the processors belonging to the source
53 // ====================================================================
54 // Creates an empty matrix structure linking two distributed supports.
55 // The method must be called by all processors belonging to source
57 // param source_support local support
58 // param source_group processor group containing the local processors
59 // param target_group processor group containing the distant processors
60 // param method interpolation method
61 // ====================================================================
63 InterpolationMatrix::InterpolationMatrix(const ParaMEDMEM::ParaFIELD *source_field,
64 const ProcessorGroup& source_group,
65 const ProcessorGroup& target_group,
66 const DECOptions& dec_options,
67 const INTERP_KERNEL::InterpolationOptions& interp_options):
68 INTERP_KERNEL::InterpolationOptions(interp_options),
69 DECOptions(dec_options),
70 _source_field(source_field),
71 _source_support(source_field->getSupport()->getCellMesh()),
72 _mapping(source_group, target_group, dec_options),
73 _source_group(source_group),
74 _target_group(target_group)
76 int nbelems = source_field->getField()->getNumberOfTuples();
77 _row_offsets.resize(nbelems+1);
78 _coeffs.resize(nbelems);
79 _target_volume.resize(nbelems);
82 InterpolationMatrix::~InterpolationMatrix()
87 // ======================================================================
88 // \brief Adds the contribution of a distant subdomain to the*
89 // interpolation matrix.
90 // The method adds contribution to the interpolation matrix.
91 // For each row of the matrix, elements are addded as
92 // a (column, coeff) pair in the _coeffs array. This column number refers
93 // to an element on the target side via the _col_offsets array.
94 // It is made of a series of (iproc, ielem) pairs.
95 // The number of elements per row is stored in the row_offsets array.
97 // param distant_support local representation of the distant subdomain
98 // param iproc_distant id of the distant subdomain (in the distant group)
99 // param distant_elems mapping between the local representation of
100 // the subdomain and the actual elem ids on the distant subdomain
101 // ======================================================================
103 void InterpolationMatrix::addContribution ( MEDCouplingPointSet& distant_support,
105 const int* distant_elems,
106 const std::string& srcMeth,
107 const std::string& targetMeth)
109 std::string interpMethod(targetMeth);
110 interpMethod+=srcMeth;
111 //creating the interpolator structure
112 vector<map<int,double> > surfaces;
113 //computation of the intersection volumes between source and target elements
114 MEDCouplingUMesh *distant_supportC=dynamic_cast<MEDCouplingUMesh *>(&distant_support);
115 MEDCouplingUMesh *source_supportC=dynamic_cast<MEDCouplingUMesh *>(_source_support);
116 if ( distant_support.getMeshDimension() == -1 )
118 if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==2)
120 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(source_supportC);
121 INTERP_KERNEL::Interpolation2D interpolation(*this);
122 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
124 else if(source_supportC->getMeshDimension()==3 && source_supportC->getSpaceDimension()==3)
126 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(source_supportC);
127 INTERP_KERNEL::Interpolation3D interpolation(*this);
128 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
130 else if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==3)
132 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(source_supportC);
133 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
134 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
137 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
139 else if ( source_supportC->getMeshDimension() == -1 )
141 if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==2)
143 MEDCouplingNormalizedUnstructuredMesh<2,2> distant_mesh_wrapper(distant_supportC);
144 INTERP_KERNEL::Interpolation2D interpolation(*this);
145 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
147 else if(distant_supportC->getMeshDimension()==3 && distant_supportC->getSpaceDimension()==3)
149 MEDCouplingNormalizedUnstructuredMesh<3,3> distant_mesh_wrapper(distant_supportC);
150 INTERP_KERNEL::Interpolation3D interpolation(*this);
151 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
153 else if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==3)
155 MEDCouplingNormalizedUnstructuredMesh<3,2> distant_mesh_wrapper(distant_supportC);
156 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
157 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
160 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
162 else if ( distant_support.getMeshDimension() == 2
163 && _source_support->getMeshDimension() == 3
164 && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
166 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
167 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
168 INTERP_KERNEL::Interpolation3D2D interpolator (*this);
169 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
170 target_wrapper.releaseTempArrays();
171 source_wrapper.releaseTempArrays();
173 else if ( distant_support.getMeshDimension() == 1
174 && _source_support->getMeshDimension() == 2
175 && distant_support.getSpaceDimension() == 2 && _source_support->getSpaceDimension() == 2)
177 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
178 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
179 INTERP_KERNEL::Interpolation2D1D interpolator (*this);
180 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
181 target_wrapper.releaseTempArrays();
182 source_wrapper.releaseTempArrays();
184 else if ( distant_support.getMeshDimension() == 3
185 && _source_support->getMeshDimension() == 1
186 && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
188 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
189 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
190 INTERP_KERNEL::Interpolation3D interpolator (*this);
191 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
192 target_wrapper.releaseTempArrays();
193 source_wrapper.releaseTempArrays();
195 else if (distant_support.getMeshDimension() != _source_support->getMeshDimension())
197 throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
199 else if( distant_support.getMeshDimension() == 1
200 && distant_support.getSpaceDimension() == 1 )
202 MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(distant_supportC);
203 MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(source_supportC);
205 INTERP_KERNEL::Interpolation1D interpolation(*this);
206 interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
207 target_wrapper.releaseTempArrays();
208 source_wrapper.releaseTempArrays();
210 else if( distant_support.getMeshDimension() == 1
211 && distant_support.getSpaceDimension() == 2 )
213 MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(distant_supportC);
214 MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(source_supportC);
216 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
217 interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
218 target_wrapper.releaseTempArrays();
219 source_wrapper.releaseTempArrays();
221 else if ( distant_support.getMeshDimension() == 2
222 && distant_support.getSpaceDimension() == 3 )
224 MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC);
225 MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC);
227 INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
228 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
229 target_wrapper.releaseTempArrays();
230 source_wrapper.releaseTempArrays();
232 else if ( distant_support.getMeshDimension() == 2
233 && distant_support.getSpaceDimension() == 2)
235 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
236 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
238 INTERP_KERNEL::Interpolation2D interpolator (*this);
239 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
240 target_wrapper.releaseTempArrays();
241 source_wrapper.releaseTempArrays();
243 else if ( distant_support.getMeshDimension() == 3
244 && distant_support.getSpaceDimension() == 3 )
246 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
247 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
249 INTERP_KERNEL::Interpolation3D interpolator (*this);
250 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
251 target_wrapper.releaseTempArrays();
252 source_wrapper.releaseTempArrays();
256 throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
258 bool needTargetSurf=isSurfaceComputationNeeded(targetMeth);
260 MEDCouplingFieldDouble *target_triangle_surf=0;
262 target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus());
263 fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf);
266 target_triangle_surf->decrRef();
269 void InterpolationMatrix::fillDSFromVM(int iproc_distant, const int* distant_elems, const std::vector< std::map<int,double> >& values, MEDCouplingFieldDouble *surf)
271 //loop over the elements to build the interpolation
273 int source_size=values.size();
274 for (int ielem=0; ielem < source_size; ielem++)
276 _row_offsets[ielem+1] += values[ielem].size();
277 for(map<int,double>::const_iterator iter=values[ielem].begin();iter!=values[ielem].end();iter++)
281 localId=distant_elems[iter->first];
284 //locating the (iproc, itriangle) pair in the list of columns
285 map<pair<int,int>,int >::iterator iter2 = _col_offsets.find(make_pair(iproc_distant,localId));
288 if (iter2 == _col_offsets.end())
290 //(iproc, itriangle) is not registered in the list
291 //of distant elements
292 col_id =_col_offsets.size();
293 _col_offsets.insert(make_pair(make_pair(iproc_distant,localId),col_id));
294 _mapping.addElementFromSource(iproc_distant,localId);
298 col_id = iter2->second;
300 //the non zero coefficient is stored
302 //col_id is the number of the column
303 //iter->second is the value of the coefficient
306 double surface = surf->getIJ(iter->first,0);
307 _target_volume[ielem].push_back(surface);
309 _coeffs[ielem].push_back(make_pair(col_id,iter->second));
314 void InterpolationMatrix::serializeMe(std::vector< std::vector< std::map<int,double> > >& data1, std::vector<int>& data2) const
318 const std::vector<std::pair<int,int> >& sendingIds=_mapping.getSendingIds();
319 std::set<int> procsS;
320 for(std::vector<std::pair<int,int> >::const_iterator iter1=sendingIds.begin();iter1!=sendingIds.end();iter1++)
321 procsS.insert((*iter1).first);
322 data1.resize(procsS.size());
323 data2.resize(procsS.size());
324 std::copy(procsS.begin(),procsS.end(),data2.begin());
325 std::map<int,int> fastProcAcc;
327 for(std::set<int>::const_iterator iter2=procsS.begin();iter2!=procsS.end();iter2++,id++)
328 fastProcAcc[*iter2]=id;
329 int nbOfSrcElt=_coeffs.size();
330 for(std::vector< std::vector< std::map<int,double> > >::iterator iter3=data1.begin();iter3!=data1.end();iter3++)
331 (*iter3).resize(nbOfSrcElt);
333 for(std::vector< std::vector< std::pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,id++)
335 for(std::vector< std::pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++)
337 const std::pair<int,int>& elt=sendingIds[(*iter5).first];
338 data1[fastProcAcc[elt.first]][id][elt.second]=(*iter5).second;
343 void InterpolationMatrix::initialize()
345 int lgth=_coeffs.size();
346 _row_offsets.clear(); _row_offsets.resize(lgth+1);
347 _coeffs.clear(); _coeffs.resize(lgth);
348 _target_volume.clear(); _target_volume.resize(lgth);
349 _col_offsets.clear();
350 _mapping.initialize();
353 void InterpolationMatrix::finishContributionW(ElementLocator& elementLocator)
355 NatureOfField nature=elementLocator.getLocalNature();
358 case ConservativeVolumic:
359 computeConservVolDenoW(elementLocator);
363 if(!elementLocator.isM1DCorr())
364 computeIntegralDenoW(elementLocator);
366 computeGlobConstraintDenoW(elementLocator);
369 case IntegralGlobConstraint:
370 computeGlobConstraintDenoW(elementLocator);
374 if(!elementLocator.isM1DCorr())
375 computeRevIntegralDenoW(elementLocator);
377 computeConservVolDenoW(elementLocator);
381 throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
386 void InterpolationMatrix::finishContributionL(ElementLocator& elementLocator)
388 NatureOfField nature=elementLocator.getLocalNature();
391 case ConservativeVolumic:
392 computeConservVolDenoL(elementLocator);
396 if(!elementLocator.isM1DCorr())
397 computeIntegralDenoL(elementLocator);
399 computeConservVolDenoL(elementLocator);
402 case IntegralGlobConstraint:
403 //this is not a bug doing like ConservativeVolumic
404 computeConservVolDenoL(elementLocator);
408 if(!elementLocator.isM1DCorr())
409 computeRevIntegralDenoL(elementLocator);
411 computeConservVolDenoL(elementLocator);
415 throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
420 void InterpolationMatrix::computeConservVolDenoW(ElementLocator& elementLocator)
422 computeGlobalColSum(_deno_reverse_multiply);
423 computeGlobalRowSum(elementLocator,_deno_multiply,_deno_reverse_multiply);
426 void InterpolationMatrix::computeConservVolDenoL(ElementLocator& elementLocator)
428 int pol1=elementLocator.sendPolicyToWorkingSideL();
429 if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
431 elementLocator.recvFromWorkingSideL();
432 elementLocator.sendToWorkingSideL();
434 else if(ElementLocator::CUMULATIVE_POLICY)
436 //ask for lazy side to deduce ids eventually missing on working side and to send it back.
437 elementLocator.recvLocalIdsFromWorkingSideL();
438 elementLocator.sendCandidatesGlobalIdsToWorkingSideL();
439 elementLocator.recvCandidatesForAddElementsL();
440 elementLocator.sendAddElementsToWorkingSideL();
441 //Working side has updated its eventually missing ids updates its global ids with lazy side procs contribution
442 elementLocator.recvLocalIdsFromWorkingSideL();
443 elementLocator.sendGlobalIdsToWorkingSideL();
444 //like no post treatment
445 elementLocator.recvFromWorkingSideL();
446 elementLocator.sendToWorkingSideL();
449 throw INTERP_KERNEL::Exception("Not managed policy detected on lazy side : not implemented !");
452 void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator)
454 MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
455 _deno_multiply.resize(_coeffs.size());
456 vector<vector<double> >::iterator iter6=_deno_multiply.begin();
457 const double *values=source_triangle_surf->getArray()->getConstPointer();
458 for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
460 (*iter6).resize((*iter4).size());
461 std::fill((*iter6).begin(),(*iter6).end(),*values);
463 source_triangle_surf->decrRef();
464 _deno_reverse_multiply=_target_volume;
467 void InterpolationMatrix::computeRevIntegralDenoW(ElementLocator& elementLocator)
469 _deno_multiply=_target_volume;
470 MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
471 _deno_reverse_multiply.resize(_coeffs.size());
472 vector<vector<double> >::iterator iter6=_deno_reverse_multiply.begin();
473 const double *values=source_triangle_surf->getArray()->getConstPointer();
474 for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
476 (*iter6).resize((*iter4).size());
477 std::fill((*iter6).begin(),(*iter6).end(),*values);
479 source_triangle_surf->decrRef();
483 * Nothing to do because surface computation is on working side.
485 void InterpolationMatrix::computeIntegralDenoL(ElementLocator& elementLocator)
490 * Nothing to do because surface computation is on working side.
492 void InterpolationMatrix::computeRevIntegralDenoL(ElementLocator& elementLocator)
497 void InterpolationMatrix::computeGlobConstraintDenoW(ElementLocator& elementLocator)
499 computeGlobalColSum(_deno_multiply);
500 computeGlobalRowSum(elementLocator,_deno_reverse_multiply,_deno_multiply);
503 void InterpolationMatrix::computeGlobalRowSum(ElementLocator& elementLocator, std::vector<std::vector<double> >& denoStrorage, std::vector<std::vector<double> >& denoStrorageInv)
505 //stores id in distant procs sorted by lazy procs connected with
506 vector< vector<int> > rowsPartialSumI;
507 //stores for each lazy procs connected with, if global info is available and if it's the case the policy
508 vector<int> policyPartial;
509 //stores the corresponding values.
510 vector< vector<double> > rowsPartialSumD;
511 elementLocator.recvPolicyFromLazySideW(policyPartial);
512 int pol1=mergePolicies(policyPartial);
513 if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
515 computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
516 elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
517 elementLocator.recvSumFromLazySideW(rowsPartialSumD);
519 else if(pol1==ElementLocator::CUMULATIVE_POLICY)
521 //updateWithNewAdditionnalElements(addingElements);
522 //stores for each lazy procs connected with, the ids in global mode if it exists (regarding policyPartial). This array has exactly the size of rowsPartialSumI,
523 //if policyPartial has CUMALATIVE_POLICY in each.
524 vector< vector<int> > globalIdsPartial;
525 computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
526 elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
527 elementLocator.recvCandidatesGlobalIdsFromLazyProcsW(globalIdsPartial);
528 std::vector< std::vector<int> > addingElements;
529 findAdditionnalElements(elementLocator,addingElements,rowsPartialSumI,globalIdsPartial);
530 addGhostElements(elementLocator.getDistantProcIds(),addingElements);
531 rowsPartialSumI.clear();
532 globalIdsPartial.clear();
533 computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
534 elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
535 elementLocator.recvGlobalIdsFromLazyProcsW(rowsPartialSumI,globalIdsPartial);
537 elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
538 elementLocator.recvSumFromLazySideW(rowsPartialSumD);
539 mergeRowSum3(globalIdsPartial,rowsPartialSumD);
540 mergeCoeffs(elementLocator.getDistantProcIds(),rowsPartialSumI,globalIdsPartial,denoStrorageInv);
543 throw INTERP_KERNEL::Exception("Not managed policy detected : not implemented !");
544 divideByGlobalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD,denoStrorage);
548 * @param distantProcs in parameter that indicates which lazy procs are concerned.
549 * @param resPerProcI out parameter that must be cleared before calling this method. The size of 1st dimension is equal to the size of 'distantProcs'.
550 * It contains the element ids (2nd dimension) of the corresponding lazy proc.
551 * @param resPerProcD out parameter with the same format than 'resPerProcI'. It contains corresponding sum values.
553 void InterpolationMatrix::computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<int> >& resPerProcI,
554 std::vector<std::vector<double> >& resPerProcD) const
556 resPerProcI.resize(distantProcs.size());
557 resPerProcD.resize(distantProcs.size());
558 std::vector<double> res(_col_offsets.size());
559 for(vector<vector<pair<int,double> > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++)
560 for(vector<pair<int,double> >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++)
561 res[(*iter3).first]+=(*iter3).second;
564 const vector<std::pair<int,int> >& mapping=_mapping.getSendingIds();
565 for(vector<std::pair<int,int> >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++)
567 std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first);
569 id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin();
570 resPerProcI[id].push_back((*iter2).second);
571 resPerProcD[id].push_back(res[iter2-mapping.begin()]);
576 * This method is only usable when CUMULATIVE_POLICY detected. This method finds elements ids (typically nodes) lazy side that
577 * are not present in columns of 'this' and that should regarding cumulative merge of elements regarding their global ids.
579 void InterpolationMatrix::findAdditionnalElements(ElementLocator& elementLocator, std::vector<std::vector<int> >& elementsToAdd,
580 const std::vector<std::vector<int> >& resPerProcI, const std::vector<std::vector<int> >& globalIdsPartial)
582 std::set<int> globalIds;
583 int nbLazyProcs=globalIdsPartial.size();
584 for(int i=0;i<nbLazyProcs;i++)
585 globalIds.insert(globalIdsPartial[i].begin(),globalIdsPartial[i].end());
586 std::vector<int> tmp(globalIds.size());
587 std::copy(globalIds.begin(),globalIds.end(),tmp.begin());
589 elementLocator.sendCandidatesForAddElementsW(tmp);
590 elementLocator.recvAddElementsFromLazyProcsW(elementsToAdd);
593 void InterpolationMatrix::addGhostElements(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& elementsToAdd)
595 std::vector< std::vector< std::map<int,double> > > data1;
596 std::vector<int> data2;
597 serializeMe(data1,data2);
599 int nbOfDistProcs=distantProcs.size();
600 for(int i=0;i<nbOfDistProcs;i++)
602 int procId=distantProcs[i];
603 const std::vector<int>& eltsForThisProc=elementsToAdd[i];
604 if(!eltsForThisProc.empty())
606 std::vector<int>::iterator iter1=std::find(data2.begin(),data2.end(),procId);
607 std::map<int,double> *toFeed=0;
608 if(iter1!=data2.end())
610 int rank=iter1-data2.begin();
611 toFeed=&(data1[rank].back());
615 iter1=std::lower_bound(data2.begin(),data2.end(),procId);
616 int rank=iter1-data2.begin();
617 data2.insert(iter1,procId);
618 std::vector< std::map<int,double> > tmp(data1.front().size());
619 data1.insert(data1.begin()+rank,tmp);
620 toFeed=&(data1[rank].back());
622 for(std::vector<int>::const_iterator iter2=eltsForThisProc.begin();iter2!=eltsForThisProc.end();iter2++)
623 (*toFeed)[*iter2]=0.;
627 nbOfDistProcs=data2.size();
628 for(int j=0;j<nbOfDistProcs;j++)
629 fillDSFromVM(data2[j],0,data1[j],0);
632 int InterpolationMatrix::mergePolicies(const std::vector<int>& policyPartial)
634 if(policyPartial.empty())
635 return ElementLocator::NO_POST_TREATMENT_POLICY;
636 int ref=policyPartial[0];
637 std::vector<int>::const_iterator iter1=std::find_if(policyPartial.begin(),policyPartial.end(),std::bind2nd(std::not_equal_to<int>(),ref));
638 if(iter1!=policyPartial.end())
640 std::ostringstream msg; msg << "Incompatible policies between lazy procs each other : proc # " << iter1-policyPartial.begin();
641 throw INTERP_KERNEL::Exception(msg.str().c_str());
647 * This method introduce global ids aspects in computed 'rowsPartialSumD'.
648 * As precondition rowsPartialSumD.size()==policyPartial.size()==globalIdsPartial.size(). Foreach i in [0;rowsPartialSumD.size() ) rowsPartialSumD[i].size()==globalIdsPartial[i].size()
649 * @param rowsPartialSumD : in parameter, Partial row sum computed for each lazy procs connected with.
650 * @param rowsPartialSumI : in parameter, Corresponding local ids for each lazy procs connected with.
651 * @param globalIdsPartial : in parameter, the global numbering, of elements connected with.
652 * @param globalIdsLazySideInteraction : out parameter, constituted from all global ids of lazy procs connected with.
653 * @para sumCorresponding : out parameter, relative to 'globalIdsLazySideInteraction'
655 void InterpolationMatrix::mergeRowSum(const std::vector< std::vector<double> >& rowsPartialSumD, const std::vector< std::vector<int> >& globalIdsPartial,
656 std::vector<int>& globalIdsLazySideInteraction, std::vector<double>& sumCorresponding)
658 std::map<int,double> sumToReturn;
659 int nbLazyProcs=rowsPartialSumD.size();
660 for(int i=0;i<nbLazyProcs;i++)
662 const std::vector<double>& rowSumOfP=rowsPartialSumD[i];
663 const std::vector<int>& globalIdsOfP=globalIdsPartial[i];
664 std::vector<double>::const_iterator iter1=rowSumOfP.begin();
665 std::vector<int>::const_iterator iter2=globalIdsOfP.begin();
666 for(;iter1!=rowSumOfP.end();iter1++,iter2++)
667 sumToReturn[*iter2]+=*iter1;
670 int lgth=sumToReturn.size();
671 globalIdsLazySideInteraction.resize(lgth);
672 sumCorresponding.resize(lgth);
673 std::vector<int>::iterator iter3=globalIdsLazySideInteraction.begin();
674 std::vector<double>::iterator iter4=sumCorresponding.begin();
675 for(std::map<int,double>::const_iterator iter5=sumToReturn.begin();iter5!=sumToReturn.end();iter5++,iter3++,iter4++)
677 *iter3=(*iter5).first;
678 *iter4=(*iter5).second;
683 * This method simply reorganize the result contained in 'sumCorresponding' computed by lazy side into 'rowsPartialSumD' with help of 'globalIdsPartial' and 'globalIdsLazySideInteraction'
685 * @param globalIdsPartial : in parameter, global ids sorted by lazy procs
686 * @param rowsPartialSumD : in/out parameter, with exactly the same size as 'globalIdsPartial'
687 * @param globalIdsLazySideInteraction : in parameter that represents ALL the global ids of every lazy procs in interaction
688 * @param sumCorresponding : in parameter with same size as 'globalIdsLazySideInteraction' that stores the corresponding sum of 'globalIdsLazySideInteraction'
690 void InterpolationMatrix::mergeRowSum2(const std::vector< std::vector<int> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD,
691 const std::vector<int>& globalIdsLazySideInteraction, const std::vector<double>& sumCorresponding)
693 std::map<int,double> acc;
694 std::vector<int>::const_iterator iter1=globalIdsLazySideInteraction.begin();
695 std::vector<double>::const_iterator iter2=sumCorresponding.begin();
696 for(;iter1!=globalIdsLazySideInteraction.end();iter1++,iter2++)
699 int nbLazyProcs=globalIdsPartial.size();
700 for(int i=0;i<nbLazyProcs;i++)
702 const std::vector<int>& tmp1=globalIdsPartial[i];
703 std::vector<double>& tmp2=rowsPartialSumD[i];
704 std::vector<int>::const_iterator iter3=tmp1.begin();
705 std::vector<double>::iterator iter4=tmp2.begin();
706 for(;iter3!=tmp1.end();iter3++,iter4++)
711 void InterpolationMatrix::mergeRowSum3(const std::vector< std::vector<int> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD)
713 std::map<int,double> sum;
714 std::vector< std::vector<int> >::const_iterator iter1=globalIdsPartial.begin();
715 std::vector< std::vector<double> >::iterator iter2=rowsPartialSumD.begin();
716 for(;iter1!=globalIdsPartial.end();iter1++,iter2++)
718 std::vector<int>::const_iterator iter3=(*iter1).begin();
719 std::vector<double>::const_iterator iter4=(*iter2).begin();
720 for(;iter3!=(*iter1).end();iter3++,iter4++)
723 iter2=rowsPartialSumD.begin();
724 for(iter1=globalIdsPartial.begin();iter1!=globalIdsPartial.end();iter1++,iter2++)
726 std::vector<int>::const_iterator iter3=(*iter1).begin();
727 std::vector<double>::iterator iter4=(*iter2).begin();
728 for(;iter3!=(*iter1).end();iter3++,iter4++)
734 * This method updates this->_coeffs attribute in order to take into account hidden (because having the same global number) similar nodes in _coeffs array.
735 * If in this->_coeffs two distant element id have the same global id their values will be replaced for each by the sum of the two.
736 * @param procsInInteraction input parameter : specifies the procId in absolute of distant lazy procs in interaction with
737 * @param rowsPartialSumI input parameter : local ids of distant lazy procs elements in interaction with
738 * @param globalIdsPartial input parameter : global ids of distant lazy procs elements in interaction with
740 void InterpolationMatrix::mergeCoeffs(const std::vector<int>& procsInInteraction, const std::vector< std::vector<int> >& rowsPartialSumI,
741 const std::vector<std::vector<int> >& globalIdsPartial, std::vector<std::vector<double> >& denoStrorageInv)
743 //preparing fast access structures
744 std::map<int,int> procT;
746 for(std::vector<int>::const_iterator iter1=procsInInteraction.begin();iter1!=procsInInteraction.end();iter1++,localProcId++)
747 procT[*iter1]=localProcId;
748 int size=procsInInteraction.size();
749 std::vector<std::map<int,int> > localToGlobal(size);
750 for(int i=0;i<size;i++)
752 std::map<int,int>& myLocalToGlobal=localToGlobal[i];
753 const std::vector<int>& locals=rowsPartialSumI[i];
754 const std::vector<int>& globals=globalIdsPartial[i];
755 std::vector<int>::const_iterator iter3=locals.begin();
756 std::vector<int>::const_iterator iter4=globals.begin();
757 for(;iter3!=locals.end();iter3++,iter4++)
758 myLocalToGlobal[*iter3]=*iter4;
761 const vector<std::pair<int,int> >& mapping=_mapping.getSendingIds();
762 std::map<int,double> globalIdVal;
763 //accumulate for same global id on lazy part.
764 for(vector<vector<pair<int,double> > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++)
765 for(vector<pair<int,double> >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
767 const std::pair<int,int>& distantLocalLazyId=mapping[(*iter2).first];
768 int localLazyProcId=procT[distantLocalLazyId.first];
769 int globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second];
770 globalIdVal[globalDistantLazyId]+=(*iter2).second;
773 std::vector<std::vector<double> >::iterator iter3=denoStrorageInv.begin();
774 for(vector<vector<pair<int,double> > >::iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter3++)
776 double val=(*iter3).back();
777 (*iter3).resize((*iter1).size());
778 std::vector<double>::iterator iter4=(*iter3).begin();
779 for(vector<pair<int,double> >::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter4++)
781 const std::pair<int,int>& distantLocalLazyId=mapping[(*iter2).first];
782 int localLazyProcId=procT[distantLocalLazyId.first];
783 int globalDistantLazyId=localToGlobal[localLazyProcId][distantLocalLazyId.second];
784 double newVal=globalIdVal[globalDistantLazyId];
785 if((*iter2).second!=0.)
786 (*iter4)=val*newVal/(*iter2).second;
788 (*iter4)=std::numeric_limits<double>::max();
789 (*iter2).second=newVal;
794 void InterpolationMatrix::divideByGlobalRowSum(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& resPerProcI,
795 const std::vector<std::vector<double> >& resPerProcD, std::vector<std::vector<double> >& deno)
797 map<int,double> fastSums;
799 for(vector<int>::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++)
801 const std::vector<int>& currentProcI=resPerProcI[procId];
802 const std::vector<double>& currentProcD=resPerProcD[procId];
803 vector<double>::const_iterator iter3=currentProcD.begin();
804 for(vector<int>::const_iterator iter2=currentProcI.begin();iter2!=currentProcI.end();iter2++,iter3++)
805 fastSums[_col_offsets[std::make_pair(*iter1,*iter2)]]=*iter3;
807 deno.resize(_coeffs.size());
808 vector<vector<double> >::iterator iter6=deno.begin();
809 for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++)
811 (*iter6).resize((*iter4).size());
812 vector<double>::iterator iter7=(*iter6).begin();
813 for(vector<pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++,iter7++)
814 *iter7=fastSums[(*iter5).first];
818 void InterpolationMatrix::computeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
820 denoStrorage.resize(_coeffs.size());
821 vector<vector<double> >::iterator iter2=denoStrorage.begin();
822 for(vector<vector<pair<int,double> > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++)
824 (*iter2).resize((*iter1).size());
825 double sumOfCurrentRow=0.;
826 for(vector<pair<int,double> >::const_iterator iter3=(*iter1).begin();iter3!=(*iter1).end();iter3++)
827 sumOfCurrentRow+=(*iter3).second;
828 std::fill((*iter2).begin(),(*iter2).end(),sumOfCurrentRow);
832 void InterpolationMatrix::resizeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
834 vector<vector<double> >::iterator iter2=denoStrorage.begin();
835 for(vector<vector<pair<int,double> > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++)
837 double val=(*iter2).back();
838 (*iter2).resize((*iter1).size());
839 std::fill((*iter2).begin(),(*iter2).end(),val);
843 // ==================================================================
844 // The call to this method updates the arrays on the target side
845 // so that they know which amount of data from which processor they
847 // That call makes actual interpolations via multiply method
849 // ==================================================================
851 void InterpolationMatrix::prepare()
853 int nbelems = _source_field->getField()->getNumberOfTuples();
854 for (int ielem=0; ielem < nbelems; ielem++)
856 _row_offsets[ielem+1]+=_row_offsets[ielem];
858 _mapping.prepareSendRecv();
862 // =======================================================================
863 // brief performs t=Ws, where t is the target field, s is the source field
865 // The call to this method must be called both on the working side
866 // and on the idle side. On the working side, the vector T=VT^(-1).(W.S)
867 // is computed and sent. On the idle side, no computation is done, but the
868 // result from the working side is received and the field is updated.
870 // param field source field on processors involved on the source side,
871 // target field on processors on the target side
872 // =======================================================================
874 void InterpolationMatrix::multiply(MEDCouplingFieldDouble& field) const
876 int nbcomp = field.getArray()->getNumberOfComponents();
877 vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
879 //computing the matrix multiply on source side
880 if (_source_group.containsMyRank())
882 int nbrows = _coeffs.size();
885 // W is the intersection matrix
886 // S is the source vector
888 for (int irow=0; irow<nbrows; irow++)
890 for (int icomp=0; icomp< nbcomp; icomp++)
892 double coeff_row = field.getIJ(irow,icomp);
893 for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
895 int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
896 double value = _coeffs[irow][icol-_row_offsets[irow]].second;
897 double deno = _deno_multiply[irow][icol-_row_offsets[irow]];
898 target_value[colid*nbcomp+icomp]+=value*coeff_row/deno;
904 if (_target_group.containsMyRank())
906 int nbelems = field.getArray()->getNumberOfTuples() ;
907 double* value = const_cast<double*> (field.getArray()->getPointer());
908 for (int i=0; i<nbelems*nbcomp; i++)
914 //on source side : sending T=VT^(-1).(W.S)
915 //on target side :: receiving T and storing it in field
916 _mapping.sendRecv(&target_value[0],field);
920 // =========================================================================
921 // brief performs s=WTt, where t is the target field, s is the source field,
922 // WT is the transpose matrix from W
924 // The call to this method must be called both on the working side
925 // and on the idle side. On the working side, the target vector T is
926 // received and the vector S=VS^(-1).(WT.T) is computed to update
928 // On the idle side, no computation is done, but the field is sent.
930 // param field source field on processors involved on the source side,
931 // target field on processors on the target side
932 // =========================================================================
934 void InterpolationMatrix::transposeMultiply(MEDCouplingFieldDouble& field) const
936 int nbcomp = field.getArray()->getNumberOfComponents();
937 vector<double> source_value(_col_offsets.size()* nbcomp,0.0);
938 _mapping.reverseSendRecv(&source_value[0],field);
940 //treatment of the transpose matrix multiply on the source side
941 if (_source_group.containsMyRank())
943 int nbrows = _coeffs.size();
944 double *array = field.getArray()->getPointer() ;
947 std::fill(array, array+nbrows*nbcomp, 0.0) ;
951 //T is the target vector
952 for (int irow = 0; irow < nbrows; irow++)
954 for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
956 int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
957 double value = _coeffs[irow][icol-_row_offsets[irow]].second;
958 double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
959 for (int icomp=0; icomp<nbcomp; icomp++)
961 double coeff_row = source_value[colid*nbcomp+icomp];
962 array[irow*nbcomp+icomp] += value*coeff_row/deno;
969 bool InterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const