1 // Copyright (C) 2007-2008 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.
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
19 #include "ParaMESH.hxx"
20 #include "ProcessorGroup.hxx"
21 #include "MxN_Mapping.hxx"
22 #include "InterpolationMatrix.hxx"
23 #include "TranslationRotationMatrix.hxx"
24 #include "Interpolation.hxx"
25 #include "Interpolation2D.txx"
26 #include "Interpolation3DSurf.txx"
27 #include "Interpolation3D.txx"
28 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
29 #include "InterpolationOptions.hxx"
30 #include "VolSurfFormulae.hxx"
31 #include "NormalizedUnstructuredMesh.hxx"
33 // class InterpolationMatrix
34 // This class enables the storage of an interpolation matrix Wij mapping
35 // source field Sj to target field Ti via Ti=Vi^(-1).Wij.Sj.
36 // The matrix is built and stored on the processors belonging to the source
44 // ====================================================================
45 // Creates an empty matrix structure linking two distributed supports.
46 // The method must be called by all processors belonging to source
48 // param source_support local support
49 // param source_group processor group containing the local processors
50 // param target_group processor group containing the distant processors
51 // param method interpolation method
52 // ====================================================================
54 InterpolationMatrix::InterpolationMatrix(ParaMEDMEM::ParaMESH *source_support,
55 const ProcessorGroup& source_group,
56 const ProcessorGroup& target_group,
57 const DECOptions& dec_options,
58 const INTERP_KERNEL::InterpolationOptions& interp_options):
59 _source_support(source_support->getCellMesh()),
60 _mapping(source_group, target_group, dec_options),
61 _source_group(source_group),
62 _target_group(target_group),
64 DECOptions(dec_options),
65 INTERP_KERNEL::InterpolationOptions(interp_options)
67 int nbelems = _source_support->getNumberOfCells();
69 _row_offsets.resize(nbelems+1);
70 for (int i=0; i<nbelems+1; i++)
75 _coeffs.resize(nbelems);
78 InterpolationMatrix::~InterpolationMatrix()
83 // ======================================================================
84 // \brief Adds the contribution of a distant subdomain to the*
85 // interpolation matrix.
86 // The method adds contribution to the interpolation matrix.
87 // For each row of the matrix, elements are addded as
88 // a (column, coeff) pair in the _coeffs array. This column number refers
89 // to an element on the target side via the _col_offsets array.
90 // It is made of a series of (iproc, ielem) pairs.
91 // The number of elements per row is stored in the row_offsets array.
93 // param distant_support local representation of the distant subdomain
94 // param iproc_distant id of the distant subdomain (in the distant group)
95 // param distant_elems mapping between the local representation of
96 // the subdomain and the actual elem ids on the distant subdomain
97 // ======================================================================
99 void InterpolationMatrix::addContribution ( MEDCouplingUMesh& distant_support,
102 const std::string& srcMeth,
103 const std::string& targetMeth)
105 if (distant_support.getMeshDimension() != _source_support->getMeshDimension())
107 throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
109 std::string interpMethod(srcMeth);
110 interpMethod+=targetMeth;
111 //creating the interpolator structure
112 vector<map<int,double> > surfaces;
114 //computation of the intersection volumes between source and target elements
116 if ( distant_support.getMeshDimension() == 2
117 && distant_support.getSpaceDimension() == 3 )
119 MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(&distant_support);
120 MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(_source_support);
122 INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
123 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
124 target_wrapper.ReleaseTempArrays();
125 source_wrapper.ReleaseTempArrays();
127 else if ( distant_support.getMeshDimension() == 2
128 && distant_support.getSpaceDimension() == 2)
130 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(&distant_support);
131 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(_source_support);
133 INTERP_KERNEL::Interpolation2D interpolator (*this);
134 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
135 target_wrapper.ReleaseTempArrays();
136 source_wrapper.ReleaseTempArrays();
138 else if ( distant_support.getMeshDimension() == 3
139 && distant_support.getSpaceDimension() ==3 )
141 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(&distant_support);
142 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(_source_support);
144 INTERP_KERNEL::Interpolation3D interpolator (*this);
145 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
146 target_wrapper.ReleaseTempArrays();
147 source_wrapper.ReleaseTempArrays();
151 throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
154 int source_size=surfaces.size();
156 MEDCouplingFieldDouble *target_triangle_surf =
157 getSupportVolumes(&distant_support);
158 MEDCouplingFieldDouble *source_triangle_surf =
159 getSupportVolumes(_source_support) ;
161 // Storing the source volumes
162 _source_volume.resize(source_size);
163 for (int i=0; i<source_size; i++)
165 _source_volume[i] = source_triangle_surf->getIJ(i,0);
168 source_triangle_surf->decrRef();
170 //loop over the elements to build the interpolation
172 for (int ielem=0; ielem < surfaces.size(); ielem++)
174 _row_offsets[ielem+1] += surfaces[ielem].size();
175 //_source_indices.push_back(make_pair(iproc_distant, ielem));
177 for (map<int,double>::const_iterator iter = surfaces[ielem].begin();
178 iter != surfaces[ielem].end();
181 //surface of the target triangle
182 double surf = target_triangle_surf->getIJ(iter->first,0);
184 //locating the (iproc, itriangle) pair in the list of columns
185 vector<pair<int,int> >::iterator iter2 =
186 find(_col_offsets.begin(), _col_offsets.end(),
187 make_pair(iproc_distant,iter->first));
190 if (iter2 == _col_offsets.end())
192 //(iproc, itriangle) is not registered in the list
193 //of distant elements
195 _col_offsets.push_back(make_pair(iproc_distant,iter->first));
196 col_id =_col_offsets.size();
197 _mapping.addElementFromSource(iproc_distant,
198 distant_elems[iter->first]);
199 _target_volume.push_back(surf);
203 col_id = iter2 - _col_offsets.begin() + 1;
206 //the non zero coefficient is stored
208 //col_id is the number of the column
209 //iter->second is the value of the coefficient
211 _coeffs[ielem].push_back(make_pair(col_id,iter->second));
214 target_triangle_surf->decrRef();
218 // ==================================================================
219 // The call to this method updates the arrays on the target side
220 // so that they know which amount of data from which processor they
222 // That call makes actual interpolations via multiply method
224 // ==================================================================
226 void InterpolationMatrix::prepare()
228 int nbelems = _source_support->getNumberOfCells();
229 for (int ielem=0; ielem < nbelems; ielem++)
231 _row_offsets[ielem+1]+=_row_offsets[ielem];
233 _mapping.prepareSendRecv();
237 // =======================================================================
238 // brief performs t=Ws, where t is the target field, s is the source field
240 // The call to this method must be called both on the working side
241 // and on the idle side. On the working side, the vector T=VT^(-1).(W.S)
242 // is computed and sent. On the idle side, no computation is done, but the
243 // result from the working side is received and the field is updated.
245 // param field source field on processors involved on the source side,
246 // target field on processors on the target side
247 // =======================================================================
249 void InterpolationMatrix::multiply(MEDCouplingFieldDouble& field) const
251 int nbcomp = field.getArray()->getNumberOfComponents();
252 vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
254 //computing the matrix multiply on source side
255 if (_source_group.containsMyRank())
257 int nbrows = _source_support->getNumberOfCells();
260 // W is the intersection matrix
261 // S is the source vector
263 for (int irow=0; irow<nbrows; irow++)
265 for (int icomp=0; icomp< nbcomp; icomp++)
267 double coeff_row = field.getIJ(irow,icomp);
268 for (int icol=_row_offsets[irow]; icol< _row_offsets[irow+1];icol++)
270 int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
271 double value = _coeffs[irow][icol-_row_offsets[irow]].second;
272 target_value[(colid-1)*nbcomp+icomp]+=value*coeff_row;
277 // performing VT^(-1).(W.S)
278 // where VT^(-1) is the inverse of the diagonal matrix containing
279 // the volumes of target cells
281 for (int i=0; i<_col_offsets.size();i++)
283 for (int icomp=0; icomp<nbcomp; icomp++)
285 target_value[i*nbcomp+icomp] /= _target_volume[i];
291 if (_target_group.containsMyRank())
293 int nbelems = field.getArray()->getNumberOfTuples() ;
294 double* value = const_cast<double*> (field.getArray()->getPointer());
295 for (int i=0; i<nbelems*nbcomp; i++)
301 //on source side : sending T=VT^(-1).(W.S)
302 //on target side :: receiving T and storing it in field
303 _mapping.sendRecv(&target_value[0],field);
307 // =========================================================================
308 // brief performs s=WTt, where t is the target field, s is the source field,
309 // WT is the transpose matrix from W
311 // The call to this method must be called both on the working side
312 // and on the idle side. On the working side, the target vector T is
313 // received and the vector S=VS^(-1).(WT.T) is computed to update
315 // On the idle side, no computation is done, but the field is sent.
317 // param field source field on processors involved on the source side,
318 // target field on processors on the target side
319 // =========================================================================
321 void InterpolationMatrix::transposeMultiply(MEDCouplingFieldDouble& field) const
323 // int nbcomp = field.getNumberOfComponents();
324 int nbcomp = field.getArray()->getNumberOfComponents();
325 vector<double> source_value(_col_offsets.size()* nbcomp,0.0);
326 _mapping.reverseSendRecv(&source_value[0],field);
328 //treatment of the transpose matrix multiply on the source side
329 if (_source_group.containsMyRank())
331 int nbrows = _source_support->getNumberOfCells();
332 double *array = field.getArray()->getPointer() ;
335 std::fill(array, array+nbrows*nbcomp, 0.0) ;
339 //T is the target vector
340 for (int irow = 0; irow < nbrows; irow++)
342 for (int icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
344 int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
345 double value = _coeffs[irow][icol-_row_offsets[irow]].second;
347 for (int icomp=0; icomp<nbcomp; icomp++)
349 double coeff_row = source_value[(colid-1)*nbcomp+icomp];
350 array[irow*nbcomp+icomp] += value*coeff_row;
355 //performing VS^(-1).(WT.T)
356 //VS^(-1) is the inverse of the diagonal matrix storing
357 //volumes of the source cells
359 for (int irow=0; irow<nbrows; irow++)
361 for (int icomp=0; icomp<nbcomp; icomp++)
363 array[irow*nbcomp+icomp] /= _source_volume[irow];
370 MEDCouplingFieldDouble* InterpolationMatrix::getSupportVolumes(MEDCouplingMesh * mesh)
372 if(!mesh->isStructured())
373 return getSupportUnstructuredVolumes((MEDCouplingUMesh *)mesh);
375 throw INTERP_KERNEL::Exception("Not implemented yet !!!");
378 // ====================================================================
379 // brief returns the volumes of the cells underlying the field \a field
381 // For 2D geometries, the returned field contains the areas.
382 // For 3D geometries, the returned field contains the volumes.
384 // param field field on which cells the volumes are required
385 // return field containing the volumes
386 // ====================================================================
388 MEDCouplingFieldDouble* InterpolationMatrix::getSupportUnstructuredVolumes(MEDCouplingUMesh * mesh)
391 int nbelem = mesh->getNumberOfCells() ;
392 int dim_mesh = mesh->getMeshDimension();
393 int dim_space = mesh->getSpaceDimension() ;
394 double *coords = mesh->getCoords()->getPointer() ;
395 int *connec = mesh->getNodalConnectivity()->getPointer() ;
396 int *connec_index = mesh->getNodalConnectivityIndex()->getPointer() ;
399 MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS);
400 DataArrayDouble* array = DataArrayDouble::New() ;
401 array->alloc(nbelem, 1) ;
402 double *area_vol = array->getPointer() ;
406 case 2: // getting the areas
407 for ( int iel=0 ; iel<nbelem ; iel++ )
409 ipt = connec_index[iel] ;
414 case INTERP_KERNEL::NORM_TRI3 :
415 case INTERP_KERNEL::NORM_TRI6 :
417 int N1 = connec[ipt+1];
418 int N2 = connec[ipt+2];
419 int N3 = connec[ipt+3];
421 area_vol[iel]=INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1),
422 coords+(dim_space*N2),
423 coords+(dim_space*N3),
428 case INTERP_KERNEL::NORM_QUAD4 :
429 case INTERP_KERNEL::NORM_QUAD8 :
431 int N1 = connec[ipt+1];
432 int N2 = connec[ipt+2];
433 int N3 = connec[ipt+3];
434 int N4 = connec[ipt+4];
436 area_vol[iel]=INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1,
444 case INTERP_KERNEL::NORM_POLYGON :
446 // We must remember that the first item is the type. That's
447 // why we substract 1 to get the number of nodes of this polygon
448 int size = connec_index[iel+1] - connec_index[iel] - 1 ;
450 double **pts = new double *[size] ;
452 for ( int inod=0 ; inod<size ; inod++ )
454 // Remember the first item is the type
455 pts[inod] = coords+dim_space*connec[ipt+inod+1] ;
458 area_vol[iel]=INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,
465 throw INTERP_KERNEL::Exception("Bad Support to get Areas on it !");
469 } // End of the loop over the cells
471 case 3: // getting the volumes
472 for ( int iel=0 ; iel<nbelem ; iel++ )
474 ipt = connec_index[iel] ;
479 case INTERP_KERNEL::NORM_TETRA4 :
480 case INTERP_KERNEL::NORM_TETRA10 :
482 int N1 = connec[ipt+1];
483 int N2 = connec[ipt+2];
484 int N3 = connec[ipt+3];
485 int N4 = connec[ipt+4];
487 area_vol[iel]=INTERP_KERNEL::calculateVolumeForTetra(coords+dim_space*N1,
490 coords+dim_space*N4) ;
494 case INTERP_KERNEL::NORM_PYRA5 :
495 case INTERP_KERNEL::NORM_PYRA13 :
497 int N1 = connec[ipt+1];
498 int N2 = connec[ipt+2];
499 int N3 = connec[ipt+3];
500 int N4 = connec[ipt+4];
501 int N5 = connec[ipt+5];
503 area_vol[iel]=INTERP_KERNEL::calculateVolumeForPyra(coords+dim_space*N1,
507 coords+dim_space*N5) ;
511 case INTERP_KERNEL::NORM_PENTA6 :
512 case INTERP_KERNEL::NORM_PENTA15 :
514 int N1 = connec[ipt+1];
515 int N2 = connec[ipt+2];
516 int N3 = connec[ipt+3];
517 int N4 = connec[ipt+4];
518 int N5 = connec[ipt+5];
519 int N6 = connec[ipt+6];
521 area_vol[iel]=INTERP_KERNEL::calculateVolumeForPenta(coords+dim_space*N1,
526 coords+dim_space*N6) ;
530 case INTERP_KERNEL::NORM_HEXA8 :
531 case INTERP_KERNEL::NORM_HEXA20 :
533 int N1 = connec[ipt+1];
534 int N2 = connec[ipt+2];
535 int N3 = connec[ipt+3];
536 int N4 = connec[ipt+4];
537 int N5 = connec[ipt+5];
538 int N6 = connec[ipt+6];
539 int N7 = connec[ipt+7];
540 int N8 = connec[ipt+8];
542 area_vol[iel]=INTERP_KERNEL::calculateVolumeForHexa(coords+dim_space*N1,
549 coords+dim_space*N8) ;
553 case INTERP_KERNEL::NORM_POLYHED :
555 throw INTERP_KERNEL::Exception("Not yet implemented !");
560 throw INTERP_KERNEL::Exception("Bad Support to get Volume on it !");
565 throw INTERP_KERNEL::Exception("interpolation is not available for this dimension");
568 field->setArray(array) ;
570 field->setMesh(mesh) ;