1 // Copyright (C) 2007-2020 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_field local source field
52 \param source_group processor group containing the local processors
53 \param target_group processor group containing the distant processors
54 \param dec_options DEC options (including projection method P0, P1)
55 \param interp_options interpolation options
57 InterpolationMatrix::InterpolationMatrix(const MEDCoupling::ParaFIELD *source_field,
58 const ProcessorGroup& source_group,
59 const ProcessorGroup& target_group,
60 const DECOptions& dec_options,
61 const INTERP_KERNEL::InterpolationOptions& interp_options):
62 INTERP_KERNEL::InterpolationOptions(interp_options),
63 DECOptions(dec_options),
64 _source_field(source_field),
65 _source_support(source_field->getSupport()->getCellMesh()),
66 _mapping(source_group, target_group, dec_options),
67 _source_group(source_group),
68 _target_group(target_group)
70 mcIdType nbelems = source_field->getField()->getNumberOfTuples();
71 _row_offsets.resize(nbelems+1);
72 _coeffs.resize(nbelems);
73 _target_volume.resize(nbelems);
76 InterpolationMatrix::~InterpolationMatrix()
82 \brief Adds the contribution of a distant subdomain to the*
84 The method adds contribution to the interpolation matrix.
85 For each row of the matrix, elements are added as
86 a (column, coeff) pair in the _coeffs array. This column number refers
87 to an element on the target side via the _col_offsets array.
88 It is made of a series of (iproc, ielem) pairs.
89 The number of elements per row is stored in the row_offsets array.
91 param distant_support local representation of the distant subdomain
92 param iproc_distant id of the distant subdomain (in the distant group)
93 param distant_elems mapping between the local representation of
94 the subdomain and the actual elem ids on the distant subdomain
96 void InterpolationMatrix::addContribution ( MEDCouplingPointSet& distant_support,
98 const mcIdType* distant_elems,
99 const std::string& srcMeth,
100 const std::string& targetMeth)
102 std::string interpMethod(targetMeth);
103 interpMethod+=srcMeth;
104 //creating the interpolator structure
105 vector<map<mcIdType,double> > surfaces;
106 //computation of the intersection volumes between source and target elements
107 MEDCouplingUMesh *distant_supportC=dynamic_cast<MEDCouplingUMesh *>(&distant_support);
108 MEDCouplingUMesh *source_supportC=dynamic_cast<MEDCouplingUMesh *>(_source_support);
109 if ( distant_support.getMeshDimension() == -1 )
111 if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==2)
113 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(source_supportC);
114 INTERP_KERNEL::Interpolation2D interpolation(*this);
115 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
117 else if(source_supportC->getMeshDimension()==3 && source_supportC->getSpaceDimension()==3)
119 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(source_supportC);
120 INTERP_KERNEL::Interpolation3D interpolation(*this);
121 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
123 else if(source_supportC->getMeshDimension()==2 && source_supportC->getSpaceDimension()==3)
125 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(source_supportC);
126 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
127 interpolation.fromIntegralUniform(source_mesh_wrapper,surfaces,srcMeth);
130 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
132 else if ( source_supportC->getMeshDimension() == -1 )
134 if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==2)
136 MEDCouplingNormalizedUnstructuredMesh<2,2> distant_mesh_wrapper(distant_supportC);
137 INTERP_KERNEL::Interpolation2D interpolation(*this);
138 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
140 else if(distant_supportC->getMeshDimension()==3 && distant_supportC->getSpaceDimension()==3)
142 MEDCouplingNormalizedUnstructuredMesh<3,3> distant_mesh_wrapper(distant_supportC);
143 INTERP_KERNEL::Interpolation3D interpolation(*this);
144 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
146 else if(distant_supportC->getMeshDimension()==2 && distant_supportC->getSpaceDimension()==3)
148 MEDCouplingNormalizedUnstructuredMesh<3,2> distant_mesh_wrapper(distant_supportC);
149 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
150 interpolation.toIntegralUniform(distant_mesh_wrapper,surfaces,srcMeth);
153 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
155 else if ( distant_support.getMeshDimension() == 2
156 && _source_support->getMeshDimension() == 3
157 && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
159 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
160 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
161 INTERP_KERNEL::Interpolation2D3D interpolator (*this);
162 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
163 target_wrapper.releaseTempArrays();
164 source_wrapper.releaseTempArrays();
166 else if ( distant_support.getMeshDimension() == 1
167 && _source_support->getMeshDimension() == 2
168 && distant_support.getSpaceDimension() == 2 && _source_support->getSpaceDimension() == 2)
170 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
171 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
172 INTERP_KERNEL::Interpolation2D1D interpolator (*this);
173 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
174 target_wrapper.releaseTempArrays();
175 source_wrapper.releaseTempArrays();
177 else if ( distant_support.getMeshDimension() == 3
178 && _source_support->getMeshDimension() == 1
179 && distant_support.getSpaceDimension() == 3 && _source_support->getSpaceDimension() == 3)
181 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
182 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
183 INTERP_KERNEL::Interpolation3D interpolator (*this);
184 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
185 target_wrapper.releaseTempArrays();
186 source_wrapper.releaseTempArrays();
188 else if (distant_support.getMeshDimension() != _source_support->getMeshDimension())
190 throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
192 else if( distant_support.getMeshDimension() == 1
193 && distant_support.getSpaceDimension() == 1 )
195 MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(distant_supportC);
196 MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(source_supportC);
198 INTERP_KERNEL::Interpolation1D interpolation(*this);
199 interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
200 target_wrapper.releaseTempArrays();
201 source_wrapper.releaseTempArrays();
203 else if( distant_support.getMeshDimension() == 1
204 && distant_support.getSpaceDimension() == 2 )
206 MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(distant_supportC);
207 MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(source_supportC);
209 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
210 interpolation.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
211 target_wrapper.releaseTempArrays();
212 source_wrapper.releaseTempArrays();
214 else if ( distant_support.getMeshDimension() == 2
215 && distant_support.getSpaceDimension() == 3 )
217 MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC);
218 MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC);
220 INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
221 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
222 target_wrapper.releaseTempArrays();
223 source_wrapper.releaseTempArrays();
225 else if ( distant_support.getMeshDimension() == 2
226 && distant_support.getSpaceDimension() == 2)
228 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
229 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
231 INTERP_KERNEL::Interpolation2D interpolator (*this);
232 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
233 target_wrapper.releaseTempArrays();
234 source_wrapper.releaseTempArrays();
236 else if ( distant_support.getMeshDimension() == 3
237 && distant_support.getSpaceDimension() == 3 )
239 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
240 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
242 INTERP_KERNEL::Interpolation3D interpolator (*this);
243 interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod);
244 target_wrapper.releaseTempArrays();
245 source_wrapper.releaseTempArrays();
249 throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
251 bool needTargetSurf=isSurfaceComputationNeeded(targetMeth);
253 MEDCouplingFieldDouble *target_triangle_surf=0;
255 target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus());
256 fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf);
259 target_triangle_surf->decrRef();
262 void InterpolationMatrix::fillDSFromVM(int iproc_distant, const mcIdType* distant_elems, const std::vector< std::map<mcIdType,double> >& values, MEDCouplingFieldDouble *surf)
264 //loop over the elements to build the interpolation
266 std::size_t source_size=values.size();
267 for (std::size_t ielem=0; ielem < source_size; ielem++)
269 _row_offsets[ielem+1] += ToIdType(values[ielem].size());
270 for(map<mcIdType,double>::const_iterator iter=values[ielem].begin();iter!=values[ielem].end();iter++)
274 localId=distant_elems[iter->first];
277 //locating the (iproc, itriangle) pair in the list of columns
278 map<pair<int,mcIdType>,mcIdType >::iterator iter2 = _col_offsets.find(make_pair(iproc_distant,localId));
281 if (iter2 == _col_offsets.end())
283 //(iproc, itriangle) is not registered in the list
284 //of distant elements
285 col_id =ToIdType(_col_offsets.size());
286 _col_offsets.insert(make_pair(make_pair(iproc_distant,localId),col_id));
287 _mapping.addElementFromSource(iproc_distant,localId);
291 col_id = iter2->second;
293 //the non zero coefficient is stored
295 //col_id is the number of the column
296 //iter->second is the value of the coefficient
299 double surface = surf->getIJ(iter->first,0);
300 _target_volume[ielem].push_back(surface);
302 _coeffs[ielem].push_back(make_pair(col_id,iter->second));
307 void InterpolationMatrix::serializeMe(std::vector< std::vector< std::map<mcIdType,double> > >& data1, std::vector<int>& data2) const
311 const std::vector<std::pair<int,mcIdType> >& sendingIds=_mapping.getSendingIds();
312 std::set<int> procsS;
313 for(std::vector<std::pair<int,mcIdType> >::const_iterator iter1=sendingIds.begin();iter1!=sendingIds.end();iter1++)
314 procsS.insert((*iter1).first);
315 data1.resize(procsS.size());
316 data2.resize(procsS.size());
317 std::copy(procsS.begin(),procsS.end(),data2.begin());
318 std::map<int,int> fastProcAcc;
320 for(std::set<int>::const_iterator iter2=procsS.begin();iter2!=procsS.end();iter2++,id++)
321 fastProcAcc[*iter2]=id;
322 mcIdType nbOfSrcElt=ToIdType(_coeffs.size());
323 for(std::vector< std::vector< std::map<mcIdType,double> > >::iterator iter3=data1.begin();iter3!=data1.end();iter3++)
324 (*iter3).resize(nbOfSrcElt);
326 for(std::vector< std::vector< std::pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,id++)
328 for(std::vector< std::pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++)
330 const std::pair<int,mcIdType>& elt=sendingIds[(*iter5).first];
331 data1[fastProcAcc[elt.first]][id][elt.second]=(*iter5).second;
336 void InterpolationMatrix::initialize()
338 mcIdType lgth=ToIdType(_coeffs.size());
339 _row_offsets.clear(); _row_offsets.resize(lgth+1);
340 _coeffs.clear(); _coeffs.resize(lgth);
341 _target_volume.clear(); _target_volume.resize(lgth);
342 _col_offsets.clear();
343 _mapping.initialize();
346 void InterpolationMatrix::finishContributionW(ElementLocator& elementLocator)
348 NatureOfField nature=elementLocator.getLocalNature();
351 case IntensiveMaximum:
352 computeConservVolDenoW(elementLocator);
354 case ExtensiveMaximum:
356 if(!elementLocator.isM1DCorr())
357 computeIntegralDenoW(elementLocator);
359 computeGlobConstraintDenoW(elementLocator);
362 case ExtensiveConservation:
363 computeGlobConstraintDenoW(elementLocator);
365 case IntensiveConservation:
367 if(!elementLocator.isM1DCorr())
368 computeRevIntegralDenoW(elementLocator);
370 computeConservVolDenoW(elementLocator);
374 throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
379 void InterpolationMatrix::finishContributionL(ElementLocator& elementLocator)
381 NatureOfField nature=elementLocator.getLocalNature();
384 case IntensiveMaximum:
385 computeConservVolDenoL(elementLocator);
387 case ExtensiveMaximum:
389 if(!elementLocator.isM1DCorr())
390 computeIntegralDenoL(elementLocator);
392 computeConservVolDenoL(elementLocator);
395 case ExtensiveConservation:
396 //this is not a bug doing like IntensiveMaximum
397 computeConservVolDenoL(elementLocator);
399 case IntensiveConservation:
401 if(!elementLocator.isM1DCorr())
402 computeRevIntegralDenoL(elementLocator);
404 computeConservVolDenoL(elementLocator);
408 throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
413 void InterpolationMatrix::computeConservVolDenoW(ElementLocator& elementLocator)
415 computeGlobalColSum(_deno_reverse_multiply);
416 computeGlobalRowSum(elementLocator,_deno_multiply,_deno_reverse_multiply);
419 void InterpolationMatrix::computeConservVolDenoL(ElementLocator& elementLocator)
421 int pol1=elementLocator.sendPolicyToWorkingSideL();
422 if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
424 elementLocator.recvFromWorkingSideL();
425 elementLocator.sendToWorkingSideL();
427 else if(ElementLocator::CUMULATIVE_POLICY)
429 //ask for lazy side to deduce ids eventually missing on working side and to send it back.
430 elementLocator.recvLocalIdsFromWorkingSideL();
431 elementLocator.sendCandidatesGlobalIdsToWorkingSideL();
432 elementLocator.recvCandidatesForAddElementsL();
433 elementLocator.sendAddElementsToWorkingSideL();
434 //Working side has updated its eventually missing ids updates its global ids with lazy side procs contribution
435 elementLocator.recvLocalIdsFromWorkingSideL();
436 elementLocator.sendGlobalIdsToWorkingSideL();
437 //like no post treatment
438 elementLocator.recvFromWorkingSideL();
439 elementLocator.sendToWorkingSideL();
442 throw INTERP_KERNEL::Exception("Not managed policy detected on lazy side : not implemented !");
445 void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator)
447 MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
448 _deno_multiply.resize(_coeffs.size());
449 vector<vector<double> >::iterator iter6=_deno_multiply.begin();
450 const double *values=source_triangle_surf->getArray()->getConstPointer();
451 for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
453 (*iter6).resize((*iter4).size());
454 std::fill((*iter6).begin(),(*iter6).end(),*values);
456 source_triangle_surf->decrRef();
457 _deno_reverse_multiply=_target_volume;
460 void InterpolationMatrix::computeRevIntegralDenoW(ElementLocator& elementLocator)
462 _deno_multiply=_target_volume;
463 MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
464 _deno_reverse_multiply.resize(_coeffs.size());
465 vector<vector<double> >::iterator iter6=_deno_reverse_multiply.begin();
466 const double *values=source_triangle_surf->getArray()->getConstPointer();
467 for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
469 (*iter6).resize((*iter4).size());
470 std::fill((*iter6).begin(),(*iter6).end(),*values);
472 source_triangle_surf->decrRef();
476 * Nothing to do because surface computation is on working side.
478 void InterpolationMatrix::computeIntegralDenoL(ElementLocator& elementLocator)
483 * Nothing to do because surface computation is on working side.
485 void InterpolationMatrix::computeRevIntegralDenoL(ElementLocator& elementLocator)
490 void InterpolationMatrix::computeGlobConstraintDenoW(ElementLocator& elementLocator)
492 computeGlobalColSum(_deno_multiply);
493 computeGlobalRowSum(elementLocator,_deno_reverse_multiply,_deno_multiply);
496 void InterpolationMatrix::computeGlobalRowSum(ElementLocator& elementLocator, std::vector<std::vector<double> >& denoStrorage, std::vector<std::vector<double> >& denoStrorageInv)
498 //stores id in distant procs sorted by lazy procs connected with
499 vector< vector<mcIdType> > rowsPartialSumI;
500 //stores for each lazy procs connected with, if global info is available and if it's the case the policy
501 vector<int> policyPartial;
502 //stores the corresponding values.
503 vector< vector<double> > rowsPartialSumD;
504 elementLocator.recvPolicyFromLazySideW(policyPartial);
505 int pol1=mergePolicies(policyPartial);
506 if(pol1==ElementLocator::NO_POST_TREATMENT_POLICY)
508 computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
509 elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
510 elementLocator.recvSumFromLazySideW(rowsPartialSumD);
512 else if(pol1==ElementLocator::CUMULATIVE_POLICY)
514 //updateWithNewAdditionnalElements(addingElements);
515 //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,
516 //if policyPartial has CUMALATIVE_POLICY in each.
517 vector< vector<mcIdType> > globalIdsPartial;
518 computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
519 elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
520 elementLocator.recvCandidatesGlobalIdsFromLazyProcsW(globalIdsPartial);
521 std::vector< std::vector<mcIdType> > addingElements;
522 findAdditionnalElements(elementLocator,addingElements,rowsPartialSumI,globalIdsPartial);
523 addGhostElements(elementLocator.getDistantProcIds(),addingElements);
524 rowsPartialSumI.clear();
525 globalIdsPartial.clear();
526 computeLocalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD);
527 elementLocator.sendLocalIdsToLazyProcsW(rowsPartialSumI);
528 elementLocator.recvGlobalIdsFromLazyProcsW(rowsPartialSumI,globalIdsPartial);
530 elementLocator.sendSumToLazySideW(rowsPartialSumI,rowsPartialSumD);
531 elementLocator.recvSumFromLazySideW(rowsPartialSumD);
532 mergeRowSum3(globalIdsPartial,rowsPartialSumD);
533 mergeCoeffs(elementLocator.getDistantProcIds(),rowsPartialSumI,globalIdsPartial,denoStrorageInv);
536 throw INTERP_KERNEL::Exception("Not managed policy detected : not implemented !");
537 divideByGlobalRowSum(elementLocator.getDistantProcIds(),rowsPartialSumI,rowsPartialSumD,denoStrorage);
541 * @param distantProcs in parameter that indicates which lazy procs are concerned.
542 * @param resPerProcI out parameter that must be cleared before calling this method. The size of 1st dimension is equal to the size of 'distantProcs'.
543 * It contains the element ids (2nd dimension) of the corresponding lazy proc.
544 * @param resPerProcD out parameter with the same format than 'resPerProcI'. It contains corresponding sum values.
546 void InterpolationMatrix::computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<mcIdType> >& resPerProcI,
547 std::vector<std::vector<double> >& resPerProcD) const
549 resPerProcI.resize(distantProcs.size());
550 resPerProcD.resize(distantProcs.size());
551 std::vector<double> res(_col_offsets.size());
552 for(vector<vector<pair<int,double> > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++)
553 for(vector<pair<int,double> >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++)
554 res[(*iter3).first]+=(*iter3).second;
557 const vector<std::pair<int,mcIdType> >& mapping=_mapping.getSendingIds();
558 for(vector<std::pair<int,mcIdType> >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++)
560 std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first);
562 id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin();
563 resPerProcI[id].push_back((*iter2).second);
564 resPerProcD[id].push_back(res[iter2-mapping.begin()]);
569 * This method is only usable when CUMULATIVE_POLICY detected. This method finds elements ids (typically nodes) lazy side that
570 * are not present in columns of 'this' and that should regarding cumulative merge of elements regarding their global ids.
572 void InterpolationMatrix::findAdditionnalElements(ElementLocator& elementLocator, std::vector<std::vector<mcIdType> >& elementsToAdd,
573 const std::vector<std::vector<mcIdType> >& resPerProcI, const std::vector<std::vector<mcIdType> >& globalIdsPartial)
575 std::set<mcIdType> globalIds;
576 std::size_t nbLazyProcs=globalIdsPartial.size();
577 for(std::size_t i=0;i<nbLazyProcs;i++)
578 globalIds.insert(globalIdsPartial[i].begin(),globalIdsPartial[i].end());
579 std::vector<mcIdType> tmp(globalIds.size());
580 std::copy(globalIds.begin(),globalIds.end(),tmp.begin());
582 elementLocator.sendCandidatesForAddElementsW(tmp);
583 elementLocator.recvAddElementsFromLazyProcsW(elementsToAdd);
586 void InterpolationMatrix::addGhostElements(const std::vector<int>& distantProcs, const std::vector<std::vector<mcIdType> >& elementsToAdd)
588 std::vector< std::vector< std::map<mcIdType,double> > > data1;
589 std::vector<int> data2;
590 serializeMe(data1,data2);
592 std::size_t nbOfDistProcs=distantProcs.size();
593 for(std::size_t i=0;i<nbOfDistProcs;i++)
595 int procId=distantProcs[i];
596 const std::vector<mcIdType>& eltsForThisProc=elementsToAdd[i];
597 if(!eltsForThisProc.empty())
599 std::vector<int>::iterator iter1=std::find(data2.begin(),data2.end(),procId);
600 std::map<mcIdType,double> *toFeed=0;
601 if(iter1!=data2.end())
603 std::size_t rank=iter1-data2.begin();
604 toFeed=&(data1[rank].back());
608 iter1=std::lower_bound(data2.begin(),data2.end(),procId);
609 std::size_t rank=iter1-data2.begin();
610 data2.insert(iter1,procId);
611 std::vector< std::map<mcIdType,double> > tmp(data1.front().size());
612 data1.insert(data1.begin()+rank,tmp);
613 toFeed=&(data1[rank].back());
615 for(std::vector<mcIdType>::const_iterator iter2=eltsForThisProc.begin();iter2!=eltsForThisProc.end();iter2++)
616 (*toFeed)[*iter2]=0.;
620 nbOfDistProcs=data2.size();
621 for(std::size_t j=0;j<nbOfDistProcs;j++)
622 fillDSFromVM(data2[j],0,data1[j],0);
625 int InterpolationMatrix::mergePolicies(const std::vector<int>& policyPartial)
627 if(policyPartial.empty())
628 return ElementLocator::NO_POST_TREATMENT_POLICY;
629 int ref=policyPartial[0];
630 std::vector<int>::const_iterator iter1=std::find_if(policyPartial.begin(),policyPartial.end(),std::bind2nd(std::not_equal_to<int>(),ref));
631 if(iter1!=policyPartial.end())
633 std::ostringstream msg; msg << "Incompatible policies between lazy procs each other : proc # " << iter1-policyPartial.begin();
634 throw INTERP_KERNEL::Exception(msg.str().c_str());
640 * This method introduce global ids aspects in computed 'rowsPartialSumD'.
641 * As precondition rowsPartialSumD.size()==policyPartial.size()==globalIdsPartial.size(). Foreach i in [0;rowsPartialSumD.size() ) rowsPartialSumD[i].size()==globalIdsPartial[i].size()
642 * @param rowsPartialSumD : in parameter, Partial row sum computed 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<mcIdType> >& globalIdsPartial,
648 std::vector<mcIdType>& globalIdsLazySideInteraction, std::vector<double>& sumCorresponding)
650 std::map<mcIdType,double> sumToReturn;
651 std::size_t nbLazyProcs=rowsPartialSumD.size();
652 for(std::size_t i=0;i<nbLazyProcs;i++)
654 const std::vector<double>& rowSumOfP=rowsPartialSumD[i];
655 const std::vector<mcIdType>& globalIdsOfP=globalIdsPartial[i];
656 std::vector<double>::const_iterator iter1=rowSumOfP.begin();
657 std::vector<mcIdType>::const_iterator iter2=globalIdsOfP.begin();
658 for(;iter1!=rowSumOfP.end();iter1++,iter2++)
659 sumToReturn[*iter2]+=*iter1;
662 std::size_t lgth=sumToReturn.size();
663 globalIdsLazySideInteraction.resize(lgth);
664 sumCorresponding.resize(lgth);
665 std::vector<mcIdType>::iterator iter3=globalIdsLazySideInteraction.begin();
666 std::vector<double>::iterator iter4=sumCorresponding.begin();
667 for(std::map<mcIdType,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<mcIdType> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD,
683 const std::vector<mcIdType>& globalIdsLazySideInteraction, const std::vector<double>& sumCorresponding)
685 std::map<mcIdType,double> acc;
686 std::vector<mcIdType>::const_iterator iter1=globalIdsLazySideInteraction.begin();
687 std::vector<double>::const_iterator iter2=sumCorresponding.begin();
688 for(;iter1!=globalIdsLazySideInteraction.end();iter1++,iter2++)
691 std::size_t nbLazyProcs=globalIdsPartial.size();
692 for(std::size_t i=0;i<nbLazyProcs;i++)
694 const std::vector<mcIdType>& tmp1=globalIdsPartial[i];
695 std::vector<double>& tmp2=rowsPartialSumD[i];
696 std::vector<mcIdType>::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<mcIdType> >& globalIdsPartial, std::vector< std::vector<double> >& rowsPartialSumD)
705 std::map<mcIdType,double> sum;
706 std::vector< std::vector<mcIdType> >::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<mcIdType>::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<mcIdType>::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<mcIdType> >& rowsPartialSumI,
733 const std::vector<std::vector<mcIdType> >& 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 std::size_t size=procsInInteraction.size();
741 std::vector<std::map<mcIdType,mcIdType> > localToGlobal(size);
742 for(std::size_t i=0;i<size;i++)
744 std::map<mcIdType,mcIdType>& myLocalToGlobal=localToGlobal[i];
745 const std::vector<mcIdType>& locals=rowsPartialSumI[i];
746 const std::vector<mcIdType>& globals=globalIdsPartial[i];
747 std::vector<mcIdType>::const_iterator iter3=locals.begin();
748 std::vector<mcIdType>::const_iterator iter4=globals.begin();
749 for(;iter3!=locals.end();iter3++,iter4++)
750 myLocalToGlobal[*iter3]=*iter4;
753 const vector<std::pair<int,mcIdType> >& mapping=_mapping.getSendingIds();
754 std::map<mcIdType,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,mcIdType>& distantLocalLazyId=mapping[(*iter2).first];
760 int localLazyProcId=procT[distantLocalLazyId.first];
761 mcIdType 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,mcIdType>& distantLocalLazyId=mapping[(*iter2).first];
774 int localLazyProcId=procT[distantLocalLazyId.first];
775 mcIdType 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<mcIdType> >& resPerProcI,
787 const std::vector<std::vector<double> >& resPerProcD, std::vector<std::vector<double> >& deno)
789 map<mcIdType,double> fastSums;
791 for(vector<int>::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++)
793 const std::vector<mcIdType>& currentProcI=resPerProcI[procId];
794 const std::vector<double>& currentProcD=resPerProcD[procId];
795 vector<double>::const_iterator iter3=currentProcD.begin();
796 for(vector<mcIdType>::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 mcIdType nbelems = _source_field->getField()->getNumberOfTuples();
845 for (mcIdType 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 mcIdType nbcomp = ToIdType(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 mcIdType nbrows = ToIdType(_coeffs.size());
876 // W is the intersection matrix
877 // S is the source vector
879 for (mcIdType irow=0; irow<nbrows; irow++)
881 for (mcIdType icomp=0; icomp< nbcomp; icomp++)
883 double coeff_row = field.getIJ(irow,icomp);
884 for (mcIdType 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 mcIdType nbelems = field.getArray()->getNumberOfTuples() ;
898 double* value = const_cast<double*> (field.getArray()->getPointer());
899 for (mcIdType 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 std::size_t 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 mcIdType nbrows = ToIdType( _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 (mcIdType irow = 0; irow < nbrows; irow++)
944 for (mcIdType 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 (std::size_t 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