1 // Copyright (C) 2007-2024 CEA, EDF
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
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 "Interpolation2D3D.txx"
33 #include "Interpolation2D1D.txx"
34 #include "MEDCouplingUMesh.hxx"
35 #include "MEDCouplingFieldDouble.hxx"
36 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
37 #include "InterpolationOptions.hxx"
38 #include "NormalizedUnstructuredMesh.hxx"
39 #include "ElementLocator.hxx"
40 #include "InterpKernelAutoPtr.hxx"
48 OverlapInterpolationMatrix::OverlapInterpolationMatrix(ParaFIELD *source_field,
49 ParaFIELD *target_field,
50 const ProcessorGroup& group,
51 const DECOptions& dec_options,
52 const INTERP_KERNEL::InterpolationOptions& i_opt,
53 const OverlapElementLocator & locator):
54 INTERP_KERNEL::InterpolationOptions(i_opt),
55 DECOptions(dec_options),
56 _source_field(source_field),
57 _target_field(target_field),
58 _source_support(source_field->getSupport()->getCellMesh()),
59 _target_support(target_field->getSupport()->getCellMesh()),
60 _mapping(group, locator)
64 void OverlapInterpolationMatrix::keepTracksOfSourceIds(int procId, DataArrayIdType *ids)
66 _mapping.keepTracksOfSourceIds(procId,ids);
69 void OverlapInterpolationMatrix::keepTracksOfTargetIds(int procId, DataArrayIdType *ids)
71 _mapping.keepTracksOfTargetIds(procId,ids);
74 OverlapInterpolationMatrix::~OverlapInterpolationMatrix()
78 // TODO? Merge with MEDCouplingRemapper::prepareInterpKernelOnlyUU() ?
80 * Local run (on this proc) of the sequential interpolation algorithm.
82 * @param srcIds is null if the source mesh is on the local proc
83 * @param trgIds is null if the source mesh is on the local proc
85 * One of the 2 is necessarily null (the two can be null together)
87 void OverlapInterpolationMatrix::computeLocalIntersection(const MEDCouplingPointSet *src, const DataArrayIdType *srcIds, const std::string& srcMeth, int srcProcId,
88 const MEDCouplingPointSet *trg, const DataArrayIdType *trgIds, const std::string& trgMeth, int trgProcId)
90 std::string interpMethod(srcMeth);
91 interpMethod+=trgMeth;
92 //creating the interpolator structure
93 vector<SparseDoubleVec > sparse_matrix_part;
95 //computation of the intersection volumes between source and target elements
96 const MEDCouplingUMesh *trgC=dynamic_cast<const MEDCouplingUMesh *>(trg);
97 const MEDCouplingUMesh *srcC=dynamic_cast<const MEDCouplingUMesh *>(src);
98 if ( src->getMeshDimension() == -1 )
100 if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==2)
102 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trgC);
103 INTERP_KERNEL::Interpolation2D interpolation(*this);
104 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth);
106 else if(trgC->getMeshDimension()==3 && trgC->getSpaceDimension()==3)
108 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(trgC);
109 INTERP_KERNEL::Interpolation3D interpolation(*this);
110 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth);
112 else if(trgC->getMeshDimension()==2 && trgC->getSpaceDimension()==3)
114 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(trgC);
115 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
116 colSize=interpolation.fromIntegralUniform(target_mesh_wrapper,sparse_matrix_part,trgMeth);
119 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
121 else if ( trg->getMeshDimension() == -1 )
123 if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==2)
125 MEDCouplingNormalizedUnstructuredMesh<2,2> local_mesh_wrapper(srcC);
126 INTERP_KERNEL::Interpolation2D interpolation(*this);
127 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth);
129 else if(srcC->getMeshDimension()==3 && srcC->getSpaceDimension()==3)
131 MEDCouplingNormalizedUnstructuredMesh<3,3> local_mesh_wrapper(srcC);
132 INTERP_KERNEL::Interpolation3D interpolation(*this);
133 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth);
135 else if(srcC->getMeshDimension()==2 && srcC->getSpaceDimension()==3)
137 MEDCouplingNormalizedUnstructuredMesh<3,2> local_mesh_wrapper(srcC);
138 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
139 colSize=interpolation.toIntegralUniform(local_mesh_wrapper,sparse_matrix_part,srcMeth);
142 throw INTERP_KERNEL::Exception("No para interpolation available for the given mesh and space dimension of distant mesh to -1D sourceMesh");
144 else if ( src->getMeshDimension() == 2 && trg->getMeshDimension() == 3
145 && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
147 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
148 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
150 INTERP_KERNEL::Interpolation2D3D interpolator (*this);
151 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
153 else if ( src->getMeshDimension() == 3 && trg->getMeshDimension() == 2
154 && trg->getSpaceDimension() == 3 && src->getSpaceDimension() == 3 )
156 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
157 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
159 INTERP_KERNEL::Interpolation2D3D interpolator (*this);
160 vector<SparseDoubleVec > matrixTranspose;
161 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,sparse_matrix_part,interpMethod);//not a bug target in source.
162 TransposeMatrix(matrixTranspose,colSize,sparse_matrix_part);
163 colSize=ToIdType(matrixTranspose.size());
165 else if ( src->getMeshDimension() == 1 && trg->getMeshDimension() == 2
166 && trg->getSpaceDimension() == 2 && src->getSpaceDimension() == 2 )
168 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
169 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
171 INTERP_KERNEL::Interpolation2D1D interpolator (*this);
172 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
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<SparseDoubleVec > matrixTranspose;
182 colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,matrixTranspose,interpMethod);//not a bug target in source.
183 TransposeMatrix(matrixTranspose,colSize,sparse_matrix_part);
184 colSize=ToIdType(matrixTranspose.size());
186 else if (trg->getMeshDimension() != _source_support->getMeshDimension())
188 throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
190 else if( src->getMeshDimension() == 1
191 && src->getSpaceDimension() == 1 )
193 MEDCouplingNormalizedUnstructuredMesh<1,1> target_wrapper(trgC);
194 MEDCouplingNormalizedUnstructuredMesh<1,1> source_wrapper(srcC);
196 INTERP_KERNEL::Interpolation1D interpolation(*this);
197 colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
199 else if( trg->getMeshDimension() == 1
200 && trg->getSpaceDimension() == 2 )
202 MEDCouplingNormalizedUnstructuredMesh<2,1> target_wrapper(trgC);
203 MEDCouplingNormalizedUnstructuredMesh<2,1> source_wrapper(srcC);
205 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
206 colSize=interpolation.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
208 else if ( trg->getMeshDimension() == 2
209 && trg->getSpaceDimension() == 3 )
211 MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(trgC);
212 MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(srcC);
214 INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
215 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
217 else if ( trg->getMeshDimension() == 2
218 && trg->getSpaceDimension() == 2)
220 MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(trgC);
221 MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(srcC);
223 INTERP_KERNEL::Interpolation2D interpolator (*this);
224 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
226 else if ( trg->getMeshDimension() == 3
227 && trg->getSpaceDimension() == 3 )
229 MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(trgC);
230 MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(srcC);
232 INTERP_KERNEL::Interpolation3D interpolator (*this);
233 colSize=interpolator.interpolateMeshes(source_wrapper,target_wrapper,sparse_matrix_part,interpMethod);
237 throw INTERP_KERNEL::Exception("No interpolator exists for these mesh and space dimensions!");
239 /* Fill distributed matrix:
240 In sparse_matrix_part rows refer to target, and columns (=first param of map in SparseDoubleVec)
243 _mapping.addContributionST(sparse_matrix_part,srcIds,srcProcId,trgIds,trgProcId);
247 * 'procsToSendField' gives the list of procs field data has to be sent to.
248 * See OverlapElementLocator::computeBoundingBoxesAndTodoList()
250 void OverlapInterpolationMatrix::prepare(const std::vector< int >& procsToSendField)
253 _mapping.prepare(procsToSendField,_target_field->getField()->getNumberOfTuplesExpected());
255 _mapping.prepare(procsToSendField,0);
258 void OverlapInterpolationMatrix::computeSurfacesAndDeno()
260 if(_target_field->getField()->getNature()==IntensiveMaximum)
261 _mapping.computeDenoConservativeVolumic(_target_field->getField()->getNumberOfTuplesExpected());
263 throw INTERP_KERNEL::Exception("OverlapDEC: Policy not set (did you call setNature()?) or not implemented yet: only IntensiveMaximum!");
265 // if(_target_field->getField()->getNature()==IntensiveConservation)
267 // MCAuto<MEDCouplingFieldDouble> f;
268 // int orient = getOrientation(); // From InterpolationOptions inheritance
269 // if(orient == 2) // absolute areas
270 // f = _target_support->getMeasureField(true);
272 // if(orient == 0) // relative areas
273 // f = _target_support->getMeasureField(false);
275 // throw INTERP_KERNEL::Exception("OverlapDEC: orientation policy not impl. yet!");
276 // _mapping.computeDenoRevIntegral(*f->getArray());
279 // throw INTERP_KERNEL::Exception("OverlapDEC: Policy not implemented yet: only IntensiveMaximum and IntensiveConservation defined!");
283 void OverlapInterpolationMatrix::multiply(double default_val)
285 _mapping.multiply(_source_field->getField(),_target_field->getField(), default_val);
288 void OverlapInterpolationMatrix::transposeMultiply()
290 _mapping.transposeMultiply(_target_field->getField(),_source_field->getField());
293 void OverlapInterpolationMatrix::TransposeMatrix(const std::vector<SparseDoubleVec >& matIn,
294 mcIdType nbColsMatIn, std::vector<SparseDoubleVec >& matOut)
296 matOut.resize(nbColsMatIn);
298 for(std::vector<SparseDoubleVec >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
299 for(SparseDoubleVec::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
300 matOut[(*iter2).first][id]=(*iter2).second;