1 // Copyright (C) 2007-2012 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
20 #include "OverlapInterpolationMatrix.hxx"
21 #include "ParaMESH.hxx"
22 #include "ParaFIELD.hxx"
23 #include "ProcessorGroup.hxx"
24 #include "TranslationRotationMatrix.hxx"
25 #include "Interpolation.hxx"
26 #include "Interpolation1D.txx"
27 #include "Interpolation2DCurve.hxx"
28 #include "Interpolation2D.txx"
29 #include "Interpolation3DSurf.hxx"
30 #include "Interpolation3D.txx"
31 #include "Interpolation3D2D.txx"
32 #include "Interpolation2D1D.txx"
33 #include "MEDCouplingUMesh.hxx"
34 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
35 #include "InterpolationOptions.hxx"
36 #include "NormalizedUnstructuredMesh.hxx"
37 #include "ElementLocator.hxx"
38 #include "InterpKernelAutoPtr.hxx"
46 OverlapInterpolationMatrix::OverlapInterpolationMatrix(ParaFIELD *source_field,
47 ParaFIELD *target_field,
48 const ProcessorGroup& group,
49 const DECOptions& dec_options,
50 const INTERP_KERNEL::InterpolationOptions& i_opt):
51 INTERP_KERNEL::InterpolationOptions(i_opt),
52 DECOptions(dec_options),
53 _source_field(source_field),
54 _target_field(target_field),
55 _source_support(source_field->getSupport()->getCellMesh()),
56 _target_support(target_field->getSupport()->getCellMesh()),
60 int nbelems = source_field->getField()->getNumberOfTuples();
61 _row_offsets.resize(nbelems+1);
62 _coeffs.resize(nbelems);
63 _target_volume.resize(nbelems);
66 void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayInt *ids)
68 _mapping.keepTracksOfSourceIds(procId,ids);
71 void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayInt *ids)
73 _mapping.keepTracksOfTargetIds(procId,ids);
76 OverlapInterpolationMatrix::~OverlapInterpolationMatrix()
80 void OverlapInterpolationMatrix::addContribution(const MEDCouplingPointSet *src, const DataArrayInt *srcIds, const std::string& srcMeth, int srcProcId,
81 const MEDCouplingPointSet *trg, const DataArrayInt *trgIds, const std::string& trgMeth, int trgProcId)
83 std::string interpMethod(srcMeth);
84 interpMethod+=trgMeth;
85 //creating the interpolator structure
86 vector<map<int,double> > surfaces;
88 //computation of the intersection volumes between source and target elements
89 const MEDCouplingUMesh *trgC=dynamic_cast<const MEDCouplingUMesh *>(trg);
90 const MEDCouplingUMesh *srcC=dynamic_cast<const MEDCouplingUMesh *>(src);
91 if ( src->getMeshDimension() == -1 )
93 if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==2)
95 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trgC);
96 INTERP_KERNEL::Interpolation2D interpolation(*this);
97 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
99 else if(trgC->getMeshDimension()==3 && trgC->getSpaceDimension()==3)
101 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(trgC);
102 INTERP_KERNEL::Interpolation3D interpolation(*this);
103 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
105 else if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==3)
107 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(trgC);
108 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
109 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
112 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
114 else if ( trg->getMeshDimension() == -1 )
116 if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==2)
118 MEDCouplingNormalizedUnstructuredMesh<2,2> local_mesh_wrapper(srcC);
119 INTERP_KERNEL::Interpolation2D interpolation(*this);
120 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
122 else if(srcC->getMeshDimension()==3 && srcC->getSpaceDimension()==3)
124 MEDCouplingNormalizedUnstructuredMesh<3,3> local_mesh_wrapper(srcC);
125 INTERP_KERNEL::Interpolation3D interpolation(*this);
126 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
128 else if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==3)
130 MEDCouplingNormalizedUnstructuredMesh<3,2> local_mesh_wrapper(srcC);
131 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
132 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
135 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
137 else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 3
138 && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
140 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
141 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
143 INTERP_KERNEL::Interpolation3D2D interpolator (*this);
144 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
145 target_wrapper.releaseTempArrays();
146 source_wrapper.releaseTempArrays();
148 else if ( src->getMeshDimension() == 3 && trg->getMeshDimension() == 2
149 && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
151 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
152 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
154 INTERP_KERNEL::Interpolation3D2D interpolator (*this);
155 vector<map<int,double> > surfacesTranspose;
156 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());//not a bug target in source.
157 TransposeMatrix(surfacesTranspose,colSize,surfaces);
158 colSize=surfacesTranspose.size();
159 target_wrapper.releaseTempArrays();
160 source_wrapper.releaseTempArrays();
162 else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2
163 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
165 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
166 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
168 INTERP_KERNEL::Interpolation2D1D interpolator (*this);
169 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
170 target_wrapper.releaseTempArrays();
171 source_wrapper.releaseTempArrays();
173 else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 1
174 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
176 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
177 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
179 INTERP_KERNEL::Interpolation2D1D interpolator (*this);
180 vector<map<int,double> > surfacesTranspose;
181 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfacesTranspose,interpMethod.c_str());//not a bug target in source.
182 TransposeMatrix(surfacesTranspose,colSize,surfaces);
183 colSize=surfacesTranspose.size();
184 target_wrapper.releaseTempArrays();
185 source_wrapper.releaseTempArrays();
187 else if (trg->getMeshDimension() != _source_support->getMeshDimension())
189 throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
191 else if( src->getMeshDimension() == 1
192 && src->getSpaceDimension() == 1 )
194 MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(trgC);
195 MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(srcC);
197 INTERP_KERNEL::Interpolation1D interpolation(*this);
198 colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
199 target_wrapper.releaseTempArrays();
200 source_wrapper.releaseTempArrays();
202 else if( trg->getMeshDimension() == 1
203 && trg->getSpaceDimension() == 2 )
205 MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(trgC);
206 MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(srcC);
208 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
209 colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
210 target_wrapper.releaseTempArrays();
211 source_wrapper.releaseTempArrays();
213 else if ( trg->getMeshDimension() == 2
214 && trg->getSpaceDimension() == 3 )
216 MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(trgC);
217 MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(srcC);
219 INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
220 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
221 target_wrapper.releaseTempArrays();
222 source_wrapper.releaseTempArrays();
224 else if ( trg->getMeshDimension() == 2
225 && trg->getSpaceDimension() == 2)
227 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
228 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
230 INTERP_KERNEL::Interpolation2D interpolator (*this);
231 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
232 target_wrapper.releaseTempArrays();
233 source_wrapper.releaseTempArrays();
235 else if ( trg->getMeshDimension() == 3
236 && trg->getSpaceDimension() == 3 )
238 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
239 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
241 INTERP_KERNEL::Interpolation3D interpolator (*this);
242 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
243 target_wrapper.releaseTempArrays();
244 source_wrapper.releaseTempArrays();
248 throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
250 bool needSourceSurf=isSurfaceComputationNeeded(srcMeth);
251 MEDCouplingFieldDouble *source_triangle_surf=0;
253 source_triangle_surf=src->getMeasureField(getMeasureAbsStatus());
255 fillDistributedMatrix(surfaces,srcIds,srcProcId,trgIds,trgProcId);
258 source_triangle_surf->decrRef();
262 * \b res rows refers to target and column (first param of map) to source.
264 void OverlapInterpolationMatrix::fillDistributedMatrix(const std::vector< std::map<int,double> >& res,
265 const DataArrayInt *srcIds, int srcProc,
266 const DataArrayInt *trgIds, int trgProc)
268 _mapping.addContributionST(res,srcIds,srcProc,trgIds,trgProc);
272 * 'procsInInteraction' gives the global view of interaction between procs.
273 * In 'procsInInteraction' for a proc with id i, is in interaction with procs listed in procsInInteraction[i]
275 void OverlapInterpolationMatrix::prepare(const std::vector< std::vector<int> >& procsInInteraction)
278 _mapping.prepare(procsInInteraction,_target_field->getField()->getNumberOfTuplesExpected());
280 _mapping.prepare(procsInInteraction,0);
283 void OverlapInterpolationMatrix::computeDeno()
285 if(_target_field->getField()->getNature()==ConservativeVolumic)
286 _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected());
288 throw INTERP_KERNEL::Exception("Policy Not implemented yet : only ConservativeVolumic defined !");
291 void OverlapInterpolationMatrix::multiply()
293 _mapping.multiply(_source_field->getField(),_target_field->getField());
296 void OverlapInterpolationMatrix::transposeMultiply()
298 _mapping.transposeMultiply(_target_field->getField(),_source_field->getField());
301 bool OverlapInterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const
306 void OverlapInterpolationMatrix::TransposeMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
308 matOut.resize(nbColsMatIn);
310 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
311 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
312 matOut[(*iter2).first][id]=(*iter2).second;