1 // Copyright (C) 2007-2013 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 // Author : Anthony Geay (CEA/DEN)
21 #include "OverlapInterpolationMatrix.hxx"
22 #include "ParaMESH.hxx"
23 #include "ParaFIELD.hxx"
24 #include "ProcessorGroup.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"
39 #include "InterpKernelAutoPtr.hxx"
47 OverlapInterpolationMatrix::OverlapInterpolationMatrix(ParaFIELD *source_field,
48 ParaFIELD *target_field,
49 const ProcessorGroup& group,
50 const DECOptions& dec_options,
51 const INTERP_KERNEL::InterpolationOptions& i_opt):
52 INTERP_KERNEL::InterpolationOptions(i_opt),
53 DECOptions(dec_options),
54 _source_field(source_field),
55 _target_field(target_field),
56 _source_support(source_field->getSupport()->getCellMesh()),
57 _target_support(target_field->getSupport()->getCellMesh()),
61 int nbelems = source_field->getField()->getNumberOfTuples();
62 _row_offsets.resize(nbelems+1);
63 _coeffs.resize(nbelems);
64 _target_volume.resize(nbelems);
67 void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayInt *ids)
69 _mapping.keepTracksOfSourceIds(procId,ids);
72 void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayInt *ids)
74 _mapping.keepTracksOfTargetIds(procId,ids);
77 OverlapInterpolationMatrix::~OverlapInterpolationMatrix()
81 void OverlapInterpolationMatrix::addContribution(const MEDCouplingPointSet *src, const DataArrayInt *srcIds, const std::string& srcMeth, int srcProcId,
82 const MEDCouplingPointSet *trg, const DataArrayInt *trgIds, const std::string& trgMeth, int trgProcId)
84 std::string interpMethod(srcMeth);
85 interpMethod+=trgMeth;
86 //creating the interpolator structure
87 vector<map<int,double> > surfaces;
89 //computation of the intersection volumes between source and target elements
90 const MEDCouplingUMesh *trgC=dynamic_cast<const MEDCouplingUMesh *>(trg);
91 const MEDCouplingUMesh *srcC=dynamic_cast<const MEDCouplingUMesh *>(src);
92 if ( src->getMeshDimension() == -1 )
94 if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==2)
96 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trgC);
97 INTERP_KERNEL::Interpolation2D interpolation(*this);
98 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
100 else if(trgC->getMeshDimension()==3 && trgC->getSpaceDimension()==3)
102 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(trgC);
103 INTERP_KERNEL::Interpolation3D interpolation(*this);
104 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
106 else if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==3)
108 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(trgC);
109 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
110 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,surfaces,trgMeth.c_str());
113 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
115 else if ( trg->getMeshDimension() == -1 )
117 if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==2)
119 MEDCouplingNormalizedUnstructuredMesh<2,2> local_mesh_wrapper(srcC);
120 INTERP_KERNEL::Interpolation2D interpolation(*this);
121 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
123 else if(srcC->getMeshDimension()==3 && srcC->getSpaceDimension()==3)
125 MEDCouplingNormalizedUnstructuredMesh<3,3> local_mesh_wrapper(srcC);
126 INTERP_KERNEL::Interpolation3D interpolation(*this);
127 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
129 else if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==3)
131 MEDCouplingNormalizedUnstructuredMesh<3,2> local_mesh_wrapper(srcC);
132 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
133 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,surfaces,srcMeth.c_str());
136 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
138 else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 3
139 && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
141 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
142 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
144 INTERP_KERNEL::Interpolation3D2D interpolator (*this);
145 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
146 target_wrapper.releaseTempArrays();
147 source_wrapper.releaseTempArrays();
149 else if ( src->getMeshDimension() == 3 && trg->getMeshDimension() == 2
150 && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
152 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
153 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
155 INTERP_KERNEL::Interpolation3D2D interpolator (*this);
156 vector<map<int,double> > surfacesTranspose;
157 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());//not a bug target in source.
158 TransposeMatrix(surfacesTranspose,colSize,surfaces);
159 colSize=surfacesTranspose.size();
160 target_wrapper.releaseTempArrays();
161 source_wrapper.releaseTempArrays();
163 else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2
164 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
166 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
167 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
169 INTERP_KERNEL::Interpolation2D1D interpolator (*this);
170 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
171 target_wrapper.releaseTempArrays();
172 source_wrapper.releaseTempArrays();
174 else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 1
175 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
177 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
178 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
180 INTERP_KERNEL::Interpolation2D1D interpolator (*this);
181 vector<map<int,double> > surfacesTranspose;
182 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfacesTranspose,interpMethod.c_str());//not a bug target in source.
183 TransposeMatrix(surfacesTranspose,colSize,surfaces);
184 colSize=surfacesTranspose.size();
185 target_wrapper.releaseTempArrays();
186 source_wrapper.releaseTempArrays();
188 else if (trg->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( src->getMeshDimension() == 1
193 && src->getSpaceDimension() == 1 )
195 MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(trgC);
196 MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(srcC);
198 INTERP_KERNEL::Interpolation1D interpolation(*this);
199 colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
200 target_wrapper.releaseTempArrays();
201 source_wrapper.releaseTempArrays();
203 else if( trg->getMeshDimension() == 1
204 && trg->getSpaceDimension() == 2 )
206 MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(trgC);
207 MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(srcC);
209 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
210 colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
211 target_wrapper.releaseTempArrays();
212 source_wrapper.releaseTempArrays();
214 else if ( trg->getMeshDimension() == 2
215 && trg->getSpaceDimension() == 3 )
217 MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(trgC);
218 MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(srcC);
220 INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
221 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
222 target_wrapper.releaseTempArrays();
223 source_wrapper.releaseTempArrays();
225 else if ( trg->getMeshDimension() == 2
226 && trg->getSpaceDimension() == 2)
228 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
229 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
231 INTERP_KERNEL::Interpolation2D interpolator (*this);
232 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
233 target_wrapper.releaseTempArrays();
234 source_wrapper.releaseTempArrays();
236 else if ( trg->getMeshDimension() == 3
237 && trg->getSpaceDimension() == 3 )
239 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
240 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
242 INTERP_KERNEL::Interpolation3D interpolator (*this);
243 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,surfaces,interpMethod.c_str());
244 target_wrapper.releaseTempArrays();
245 source_wrapper.releaseTempArrays();
249 throw INTERP_KERNEL::Exception("no interpolator exists for these mesh and space dimensions ");
251 bool needSourceSurf=isSurfaceComputationNeeded(srcMeth);
252 MEDCouplingFieldDouble *source_triangle_surf=0;
254 source_triangle_surf=src->getMeasureField(getMeasureAbsStatus());
256 fillDistributedMatrix(surfaces,srcIds,srcProcId,trgIds,trgProcId);
259 source_triangle_surf->decrRef();
263 * \b res rows refers to target and column (first param of map) to source.
265 void OverlapInterpolationMatrix::fillDistributedMatrix(const std::vector< std::map<int,double> >& res,
266 const DataArrayInt *srcIds, int srcProc,
267 const DataArrayInt *trgIds, int trgProc)
269 _mapping.addContributionST(res,srcIds,srcProc,trgIds,trgProc);
273 * 'procsInInteraction' gives the global view of interaction between procs.
274 * In 'procsInInteraction' for a proc with id i, is in interaction with procs listed in procsInInteraction[i]
276 void OverlapInterpolationMatrix::prepare(const std::vector< std::vector<int> >& procsInInteraction)
279 _mapping.prepare(procsInInteraction,_target_field->getField()->getNumberOfTuplesExpected());
281 _mapping.prepare(procsInInteraction,0);
284 void OverlapInterpolationMatrix::computeDeno()
286 if(_target_field->getField()->getNature()==ConservativeVolumic)
287 _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected());
289 throw INTERP_KERNEL::Exception("Policy Not implemented yet : only ConservativeVolumic defined !");
292 void OverlapInterpolationMatrix::multiply()
294 _mapping.multiply(_source_field->getField(),_target_field->getField());
297 void OverlapInterpolationMatrix::transposeMultiply()
299 _mapping.transposeMultiply(_target_field->getField(),_source_field->getField());
302 bool OverlapInterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const
307 void OverlapInterpolationMatrix::TransposeMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
309 matOut.resize(nbColsMatIn);
311 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
312 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
313 matOut[(*iter2).first][id]=(*iter2).second;