1 // Copyright (C) 2007-2016 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
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingRemapper.hxx"
22 #include "MEDCouplingMemArray.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24 #include "MEDCouplingFieldTemplate.hxx"
25 #include "MEDCouplingFieldDiscretization.hxx"
26 #include "MEDCouplingMappedExtrudedMesh.hxx"
27 #include "MEDCouplingCMesh.hxx"
28 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
29 #include "MEDCouplingNormalizedCartesianMesh.txx"
31 #include "Interpolation1D.txx"
32 #include "Interpolation2DCurve.hxx"
33 #include "Interpolation2D.txx"
34 #include "Interpolation3D.txx"
35 #include "Interpolation3DSurf.hxx"
36 #include "Interpolation2D1D.txx"
37 #include "Interpolation2D3D.txx"
38 #include "Interpolation3D1D.txx"
39 #include "Interpolation1D0D.txx"
40 #include "InterpolationCU.txx"
41 #include "InterpolationCC.txx"
43 using namespace MEDCoupling;
45 MEDCouplingRemapper::MEDCouplingRemapper():_src_ft(0),_target_ft(0),_interp_matrix_pol(IK_ONLY_PREFERED),_nature_of_deno(NoNature),_time_deno_update(0)
49 MEDCouplingRemapper::~MEDCouplingRemapper()
55 * This method is the second step of the remapping process. The remapping process works in three phases :
57 * - Set remapping options appropriately
58 * - The computation of remapping matrix
59 * - Apply the matrix vector multiply to obtain the result of the remapping
61 * This method performs the second step (computation of remapping matrix) which may be CPU-time consuming. This phase is also the most critical (where the most tricky algorithm) in the remapping process.
62 * Strictly speaking to perform the computation of the remapping matrix the field templates source-side and target-side is required (which is the case of MEDCouplingRemapper::prepareEx).
63 * So this method is less precise but a more user friendly way to compute a remapping matrix.
65 * \param [in] srcMesh the source mesh
66 * \param [in] targetMesh the target mesh
67 * \param [in] method A string obtained by aggregation of the spatial discretisation string representation of source field and target field. The string representation is those returned by MEDCouplingFieldDiscretization::getStringRepr.
68 * Example : "P0" is for cell discretization. "P1" is for node discretization. So "P0P1" for \a method parameter means from a source cell field (lying on \a srcMesh) to a target node field (lying on \a targetMesh).
70 * \sa MEDCouplingRemapper::prepareEx
72 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
74 if(!srcMesh || !targetMesh)
75 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepare : presence of NULL input pointer !");
76 std::string srcMethod,targetMethod;
77 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
78 MCAuto<MEDCouplingFieldTemplate> src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
79 src->setMesh(srcMesh);
80 MCAuto<MEDCouplingFieldTemplate> target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
81 target->setMesh(targetMesh);
82 return prepareEx(src,target);
86 * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
87 * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
89 * \param [in] src is the field template source side.
90 * \param [in] target is the field template target side.
92 * \sa MEDCouplingRemapper::prepare
94 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
97 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL input pointer !");
98 if(!src->getMesh() || !target->getMesh())
99 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
101 _src_ft=const_cast<MEDCouplingFieldTemplate *>(src); _src_ft->incrRef();
102 _target_ft=const_cast<MEDCouplingFieldTemplate *>(target); _target_ft->incrRef();
103 if(isInterpKernelOnlyOrNotOnly())
104 return prepareInterpKernelOnly();
106 return prepareNotInterpKernelOnly();
109 int MEDCouplingRemapper::prepareInterpKernelOnly()
111 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
119 // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
120 // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
122 // } MEDCouplingMeshType;
124 switch(meshInterpType)
126 case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
127 case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
128 case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
129 case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
130 case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
131 case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
132 case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
133 case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
134 case 85: // UNSTRUCTURED - UNSTRUCTURED
135 return prepareInterpKernelOnlyUU();
136 case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
137 case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
138 case 87: // UNSTRUCTURED - CARTESIAN
139 return prepareInterpKernelOnlyUC();
140 case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
141 case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
142 case 117: // CARTESIAN - UNSTRUCTURED
143 return prepareInterpKernelOnlyCU();
144 case 119: // CARTESIAN - CARTESIAN
145 return prepareInterpKernelOnlyCC();
146 case 136: // EXTRUDED - EXTRUDED
147 return prepareInterpKernelOnlyEE();
149 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
153 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
155 std::string srcm,trgm,method;
156 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
157 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
160 return prepareNotInterpKernelOnlyGaussGauss();
163 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
164 throw INTERP_KERNEL::Exception(oss.str().c_str());
170 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
171 * If meshes of \b srcField and \b targetField do not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
173 * \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on MEDCoupling::MEDCouplingRemapper::prepare.
174 * \param [in/out] targetField the destination field with the allocated array in which all tuples will be overwritten.
175 * \param [in] dftValue is the value that will be assigned in the targetField to each entity of target mesh (entity depending on the method selected on prepare invocation) that is not intercepted by any entity of source mesh.
176 * For example in "P0P0" case (cell-cell) if a cell in target mesh is not overlapped by any source cell the \a dftValue value will be attached on that cell in the returned \a targetField. In some cases a target
177 * cell not intercepted by any source cell is a bug so in this case it is advised to set a huge value (1e300 for example) to \a dftValue to quickly point to the problem. But for users doing parallelism a target cell can
178 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
182 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue)
184 if(!srcField || !targetField)
185 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transfer : input field must be both not NULL !");
186 transferUnderground(srcField,targetField,true,dftValue);
190 * This method is equivalent to MEDCoupling::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
191 * If an entity (cell for example) in targetField is not fetched by any entity (cell for example) of \b srcField, the value in targetField is
193 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
195 * \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on MEDCoupling::MEDCouplingRemapper::prepare.
196 * \param [in,out] targetField the destination field with the allocated array in which only tuples whose entities are fetched by interpolation will be overwritten only.
198 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField)
200 if(!srcField || !targetField)
201 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : input field must be both not NULL !");
202 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
205 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue)
207 if(!srcField || !targetField)
208 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::reverseTransfer : input fields must be both not NULL !");
210 targetField->checkConsistencyLight();
211 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
212 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
213 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
214 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
215 if(srcField->getNature()!=targetField->getNature())
216 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
217 if(targetField->getNumberOfTuplesExpected()!=_target_ft->getNumberOfTuplesExpected())
219 std::ostringstream oss;
220 oss << "MEDCouplingRemapper::reverseTransfer : in given source field the number of tuples required is " << _target_ft->getNumberOfTuplesExpected() << " (on prepare) and number of tuples in given target field is " << targetField->getNumberOfTuplesExpected();
221 oss << " ! It appears that the target support is not the same between the prepare and the transfer !";
222 throw INTERP_KERNEL::Exception(oss.str().c_str());
224 DataArrayDouble *array(srcField->getArray());
225 int trgNbOfCompo=targetField->getNumberOfComponents();
228 srcField->checkConsistencyLight();
229 if(trgNbOfCompo!=srcField->getNumberOfTuplesExpected())
230 throw INTERP_KERNEL::Exception("Number of components mismatch !");
234 MCAuto<DataArrayDouble > tmp(DataArrayDouble::New());
235 tmp->alloc(srcField->getNumberOfTuplesExpected(),trgNbOfCompo);
236 srcField->setArray(tmp);
238 computeDeno(srcField->getNature(),srcField,targetField);
239 double *resPointer(srcField->getArray()->getPointer());
240 const double *inputPointer=targetField->getArray()->getConstPointer();
241 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
245 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
246 * If mesh of \b srcField does not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
248 * \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on MEDCoupling::MEDCouplingRemapper::prepare.
249 * \param [in] dftValue is the value that will be assigned in the targetField to each entity of target mesh (entity depending on the method selected on prepare invocation) that is not intercepted by any entity of source mesh.
250 * For example in "P0P0" case (cell-cell) if a cell in target mesh is not overlapped by any source cell the \a dftValue value will be attached on that cell in the returned \a targetField. In some cases a target
251 * cell not intercepted by any source cell is a bug so in this case it is advised to set a huge value (1e300 for example) to \a dftValue to quickly point to the problem. But for users doing parallelism a target cell can
252 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
253 * \return destination field to be deallocated by the caller.
257 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue)
261 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input srcField is NULL !");
262 srcField->checkConsistencyLight();
263 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
264 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
265 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
266 ret->setNature(srcField->getNature());
267 transfer(srcField,ret,dftValue);
268 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
272 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue)
275 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input targetField is NULL !");
276 targetField->checkConsistencyLight();
278 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
279 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
280 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
281 ret->setNature(targetField->getNature());
282 reverseTransfer(ret,targetField,dftValue);
283 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
288 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
289 * is here only for automatic CORBA generators.
291 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
293 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
297 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
298 * is here only for automatic CORBA generators.
300 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
302 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
306 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
307 * is here only for automatic CORBA generators.
309 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
311 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
315 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or preferred.
316 * If interpolation matrix policy is :
318 * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is preferred. That is to say, if it is possible to treat the case
319 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
320 * If not, the \b not only INTERP_KERNEL method will be attempt.
322 * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is preferred. That is to say, if it is possible to treat the case
323 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
324 * If not, the INTERP_KERNEL only method will be attempt.
326 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
328 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
330 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
332 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
334 return _interp_matrix_pol;
338 * This method sets a new interpolation matrix policy. The default one is IK_PREFERED (0). The input is of type \c int to be dealt by standard Salome
339 * CORBA component generators. This method throws an INTERP_KERNEL::Exception if a the input integer is not in the available possibilities, that is to say not in
340 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
342 * If interpolation matrix policy is :
344 * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is preferred. That is to say, if it is possible to treat the case
345 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
346 * If not, the \b not only INTERP_KERNEL method will be attempt.
348 * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is preferred. That is to say, if it is possible to treat the case
349 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
350 * If not, the INTERP_KERNEL only method will be attempt.
352 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
354 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
356 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c MEDCoupling::InterpolationMatrixPolicy
357 * for automatic generation of CORBA component.
359 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
361 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
363 switch(newInterpMatPol)
366 _interp_matrix_pol=IK_ONLY_PREFERED;
369 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
372 _interp_matrix_pol=IK_ONLY_FORCED;
375 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
378 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::setInterpolationMatrixPolicy : invalid input integer value ! Should be in [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)] ! For information, the default is IK_PREFERED=0 !");
382 int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
384 const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
385 const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_target_ft->getMesh());
386 std::string srcMeth,trgMeth;
387 std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth));
388 const int srcMeshDim=src_mesh->getMeshDimension();
391 srcSpaceDim=src_mesh->getSpaceDimension();
392 const int trgMeshDim=target_mesh->getMeshDimension();
395 trgSpaceDim=target_mesh->getSpaceDimension();
396 if(trgSpaceDim!=srcSpaceDim)
397 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
398 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
400 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
402 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
403 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
404 INTERP_KERNEL::Interpolation1D interpolation(*this);
405 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
407 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
409 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
410 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
411 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
412 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
414 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
416 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
417 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
418 INTERP_KERNEL::Interpolation2D interpolation(*this);
419 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
421 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
423 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
424 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
425 INTERP_KERNEL::Interpolation3D interpolation(*this);
426 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
428 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
430 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
431 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
432 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
433 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
435 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
437 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
438 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
439 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
440 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
441 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
442 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
444 else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
446 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
447 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! Select PointLocator as intersection type !");
448 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
449 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
450 INTERP_KERNEL::Interpolation1D0D interpolation(*this);
451 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
453 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
455 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
456 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
457 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
458 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
459 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
460 std::vector<std::map<int,double> > matrixTmp;
461 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
462 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
463 ReverseMatrix(matrixTmp,nbCols,_matrix);
464 nbCols=matrixTmp.size();
466 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
468 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
470 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
471 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
472 INTERP_KERNEL::Interpolation2D interpolation(*this);
473 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
477 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
478 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
479 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
480 std::vector<std::map<int,double> > matrixTmp;
481 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
482 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
483 ReverseMatrix(matrixTmp,nbCols,_matrix);
484 nbCols=matrixTmp.size();
485 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
486 if(!duplicateFaces.empty())
488 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
489 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
491 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
492 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
498 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
500 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
502 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
503 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
504 INTERP_KERNEL::Interpolation2D interpolation(*this);
505 std::vector<std::map<int,double> > matrixTmp;
506 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
507 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
508 ReverseMatrix(matrixTmp,nbCols,_matrix);
509 nbCols=matrixTmp.size();
513 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
514 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
515 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
516 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
517 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
518 if(!duplicateFaces.empty())
520 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
521 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
523 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
524 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
530 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
532 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
533 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
534 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
535 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
536 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
537 if(!duplicateFaces.empty())
539 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
540 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
542 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
543 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
548 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
550 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
552 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
553 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
554 INTERP_KERNEL::Interpolation3D interpolation(*this);
555 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
559 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
560 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
561 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
562 std::vector<std::map<int,double> > matrixTmp;
563 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
564 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
565 ReverseMatrix(matrixTmp,nbCols,_matrix);
566 nbCols=matrixTmp.size();
567 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
568 if(!duplicateFaces.empty())
570 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
571 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
573 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
574 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
580 else if(trgMeshDim==-1)
582 if(srcMeshDim==2 && srcSpaceDim==2)
584 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
585 INTERP_KERNEL::Interpolation2D interpolation(*this);
586 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
588 else if(srcMeshDim==3 && srcSpaceDim==3)
590 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
591 INTERP_KERNEL::Interpolation3D interpolation(*this);
592 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
594 else if(srcMeshDim==2 && srcSpaceDim==3)
596 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
597 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
598 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
601 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
603 else if(srcMeshDim==-1)
605 if(trgMeshDim==2 && trgSpaceDim==2)
607 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
608 INTERP_KERNEL::Interpolation2D interpolation(*this);
609 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
611 else if(trgMeshDim==3 && trgSpaceDim==3)
613 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
614 INTERP_KERNEL::Interpolation3D interpolation(*this);
615 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
617 else if(trgMeshDim==2 && trgSpaceDim==3)
619 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
620 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
621 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
624 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
627 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
628 _deno_multiply.clear();
629 _deno_multiply.resize(_matrix.size());
630 _deno_reverse_multiply.clear();
631 _deno_reverse_multiply.resize(nbCols);
636 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
638 std::string srcMeth,trgMeth;
639 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
640 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
641 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
643 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
644 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
645 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
646 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
647 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
648 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
649 std::vector<std::map<int,double> > matrix2D;
650 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
651 MEDCouplingUMesh *s1D,*t1D;
653 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
654 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
655 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
656 std::vector<std::map<int,double> > matrix1D;
657 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
658 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
659 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
660 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
663 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
664 target_mesh->getMesh3DIds()->getConstPointer());
666 _deno_multiply.clear();
667 _deno_multiply.resize(_matrix.size());
668 _deno_reverse_multiply.clear();
669 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
674 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
676 std::string srcMeth,trgMeth;
677 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
678 if(methodCpp!="P0P0")
679 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
680 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
681 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
682 const int srcMeshDim=src_mesh->getMeshDimension();
683 const int srcSpceDim=src_mesh->getSpaceDimension();
684 const int trgMeshDim=target_mesh->getMeshDimension();
685 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
686 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : space dim of src unstructured should be equal to mesh dim of src unstructured and should be equal also equal to trg cartesian dimension !");
687 std::vector<std::map<int,double> > res;
692 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
693 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
694 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
695 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
700 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
701 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
702 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
703 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
708 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
709 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
710 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
711 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
715 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
717 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
718 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
720 _deno_multiply.clear();
721 _deno_multiply.resize(_matrix.size());
722 _deno_reverse_multiply.clear();
723 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
728 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
730 std::string srcMeth,trgMeth;
731 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
732 if(methodCpp!="P0P0")
733 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
734 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
735 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
736 const int srcMeshDim=src_mesh->getMeshDimension();
737 const int trgMeshDim=target_mesh->getMeshDimension();
738 const int trgSpceDim=target_mesh->getSpaceDimension();
739 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
740 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : space dim of target unstructured should be equal to mesh dim of target unstructured and should be equal also equal to source cartesian dimension !");
745 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
746 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
747 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
748 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
753 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
754 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
755 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
756 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
761 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
762 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
763 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
764 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
768 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
770 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
772 _deno_multiply.clear();
773 _deno_multiply.resize(_matrix.size());
774 _deno_reverse_multiply.clear();
775 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
780 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
782 std::string srcMeth,trgMeth;
783 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
784 if(methodCpp!="P0P0")
785 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
786 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
787 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
788 const int srcMeshDim=src_mesh->getMeshDimension();
789 const int trgMeshDim=target_mesh->getMeshDimension();
790 if(trgMeshDim!=srcMeshDim)
791 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
796 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
797 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
798 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
799 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
804 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
805 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
806 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
807 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
812 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
813 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
814 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
815 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
819 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
821 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
823 _deno_multiply.clear();
824 _deno_multiply.resize(_matrix.size());
825 _deno_reverse_multiply.clear();
826 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
831 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
833 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
834 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : The intersection type is not supported ! Only PointLocator is supported for Gauss->Gauss interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !");
835 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
836 const double *trgLocPtr=trgLoc->begin();
837 int trgSpaceDim=trgLoc->getNumberOfComponents();
838 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
839 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
841 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
842 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
843 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
844 throw INTERP_KERNEL::Exception(oss.str().c_str());
846 const int *srcOffsetArrPtr=srcOffsetArr->begin();
847 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
848 const double *srcLocPtr=srcLoc->begin();
849 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
850 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
851 _matrix.resize(trgNbOfGaussPts);
852 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
853 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
854 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
855 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
856 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
858 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
859 int srcCellId=elts[eltsIndex[*trgId]];
860 double dist=std::numeric_limits<double>::max();
862 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
864 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
866 for(int i=0;i<trgSpaceDim;i++)
867 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
869 { dist=tmp; srcEntry=srcId; }
871 _matrix[*trgId][srcEntry]=1.;
873 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
875 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
876 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
877 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
878 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
879 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
880 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
882 _deno_multiply.clear();
883 _deno_multiply.resize(_matrix.size());
884 _deno_reverse_multiply.clear();
885 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
891 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
892 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
894 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
896 if(method=="GAUSSGAUSS")
898 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
899 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
900 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
901 throw INTERP_KERNEL::Exception(oss.str().c_str());
905 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
906 * to IK_ONLY_PREFERED, which method will be applied. If \c true is returned the INTERP_KERNEL only method should be applied to \c false the \b not
907 * only INTERP_KERNEL method should be applied.
909 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
911 std::string srcm,trgm,method;
912 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
913 switch(_interp_matrix_pol)
915 case IK_ONLY_PREFERED:
919 std::string tmp1,tmp2;
920 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
923 catch(INTERP_KERNEL::Exception& /*e*/)
928 case NOT_IK_ONLY_PREFERED:
932 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
935 catch(INTERP_KERNEL::Exception& /*e*/)
942 case NOT_IK_ONLY_FORCED:
945 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
949 void MEDCouplingRemapper::updateTime() const
953 void MEDCouplingRemapper::checkPrepare() const
955 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
957 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
958 if(!s->getMesh() || !t->getMesh())
959 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
963 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
964 * This method returns 3 information (2 in output parameters and 1 in return).
966 * \param [out] srcMeth the string code of the discretization of source field template
967 * \param [out] trgMeth the string code of the discretization of target field template
968 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
970 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
972 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
974 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
975 if(!s->getMesh() || !t->getMesh())
976 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
977 srcMeth=_src_ft->getDiscretization()->getRepr();
978 trgMeth=_target_ft->getDiscretization()->getRepr();
979 return BuildMethodFrom(srcMeth,trgMeth);
982 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
984 std::string method(meth1); method+=meth2;
988 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
992 if(matrixSuppression)
995 _deno_multiply.clear();
996 _deno_reverse_multiply.clear();
1000 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1002 if(!srcField || !targetField)
1003 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1004 srcField->checkConsistencyLight();
1006 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1007 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1008 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1009 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1010 if(srcField->getNature()!=targetField->getNature())
1011 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1012 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1014 std::ostringstream oss;
1015 oss << "MEDCouplingRemapper::transferUnderground : in given source field the number of tuples required is " << _src_ft->getNumberOfTuplesExpected() << " (on prepare) and number of tuples in given source field is " << srcField->getNumberOfTuplesExpected();
1016 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1017 throw INTERP_KERNEL::Exception(oss.str().c_str());
1019 DataArrayDouble *array(targetField->getArray());
1020 int srcNbOfCompo(srcField->getNumberOfComponents());
1023 targetField->checkConsistencyLight();
1024 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1025 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1030 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1031 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1032 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1033 targetField->setArray(tmp);
1035 computeDeno(srcField->getNature(),srcField,targetField);
1036 double *resPointer(targetField->getArray()->getPointer());
1037 const double *inputPointer(srcField->getArray()->getConstPointer());
1038 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1041 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1044 return computeDenoFromScratch(nat,srcField,trgField);
1045 else if(nat!=_nature_of_deno)
1046 return computeDenoFromScratch(nat,srcField,trgField);
1047 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1048 return computeDenoFromScratch(nat,srcField,trgField);
1051 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1053 _nature_of_deno=nat;
1054 _time_deno_update=getTimeOfThis();
1055 switch(_nature_of_deno)
1057 case IntensiveMaximum:
1059 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1062 case ExtensiveMaximum:
1064 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1065 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1066 const double *denoPtr=deno->getArray()->getConstPointer();
1067 const double *denoRPtr=denoR->getArray()->getConstPointer();
1068 if(trgField->getMesh()->getMeshDimension()==-1)
1070 double *denoRPtr2=denoR->getArray()->getPointer();
1071 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1073 if(srcField->getMesh()->getMeshDimension()==-1)
1075 double *denoPtr2=deno->getArray()->getPointer();
1076 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1079 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1080 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1082 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1083 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1089 case ExtensiveConservation:
1091 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1094 case IntensiveConservation:
1096 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1097 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1098 const double *denoPtr=deno->getArray()->getConstPointer();
1099 const double *denoRPtr=denoR->getArray()->getConstPointer();
1100 if(trgField->getMesh()->getMeshDimension()==-1)
1102 double *denoRPtr2=denoR->getArray()->getPointer();
1103 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1105 if(srcField->getMesh()->getMeshDimension()==-1)
1107 double *denoPtr2=deno->getArray()->getPointer();
1108 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1111 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1112 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1114 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1115 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1122 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1126 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1129 double *tmp=new double[inputNbOfCompo];
1130 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1132 if((*iter1).empty())
1135 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1139 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1140 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1141 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1143 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1144 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1150 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1152 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1154 double *tmp=new double[inputNbOfCompo];
1155 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1156 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1158 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1160 isReached[(*iter2).first]=true;
1161 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1162 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1167 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1169 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1172 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1174 matOut.resize(nbColsMatIn);
1176 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1177 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1178 matOut[(*iter2).first][id]=(*iter2).second;
1181 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1182 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1184 std::map<int,double> values;
1186 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1189 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1191 sum+=(*iter2).second;
1192 values[(*iter2).first]+=(*iter2).second;
1194 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1195 deno[idx][(*iter2).first]=sum;
1198 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1200 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1201 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1205 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1206 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1208 std::map<int,double> values;
1210 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1213 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1215 sum+=(*iter2).second;
1216 values[(*iter2).first]+=(*iter2).second;
1218 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1219 denoReverse[(*iter2).first][idx]=sum;
1222 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1224 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1225 deno[idx][(*iter2).first]=values[(*iter2).first];
1229 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1230 const std::vector< std::map<int,double> >& m2D,
1231 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1232 const int *corrCellIdTrg)
1234 int nbOf2DCellsTrg=m2D.size();
1235 int nbOf1DCellsTrg=m1D.size();
1236 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1237 _matrix.resize(nbOf3DCellsTrg);
1239 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1241 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1244 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1246 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1248 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1255 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1258 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1260 std::cout << "Target Cell # " << id << " : ";
1261 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1262 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1263 std::cout << std::endl;
1267 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1273 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1275 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1277 return (int)_deno_reverse_multiply.size();
1281 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1282 * If not the behaviour is unpredictable.
1283 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1284 * set to 0. That is to say that its entry disappear from the map storing the corresponding row in the data storage of sparse crude matrix.
1285 * This method is useful to correct at a high level some problems linked to precision. Indeed, with some \ref NatureOfField "natures of field" some threshold effect
1288 * \param [in] maxValAbs is a limit behind which a coefficient is set to 0. \a maxValAbs is expected to be positive, if not this method do nothing.
1289 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1290 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1292 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1295 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1297 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1299 std::map<int,double>& rowNew=matrixNew[i];
1300 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1302 if(fabs((*it2).second)>maxValAbs)
1303 rowNew[(*it2).first]=(*it2).second;
1314 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1315 * If not the behaviour is unpredictable.
1316 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1317 * set to 0. That is to say that its entry disappear from the map storing the corresponding row in the data storage of sparse crude matrix.
1318 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1319 * This method is useful to correct at a high level some problems linked to precision. Indeed, with some \ref NatureOfField "natures of field" some threshold effect
1322 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1323 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged. If -1 is returned it means
1324 * that all coefficients are null.
1325 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1327 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1329 double maxVal=getMaxValueInCrudeMatrix();
1332 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1336 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1337 * If not the behaviour is unpredictable.
1338 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1339 * The returned value is positive.
1341 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1344 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1345 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1346 if(fabs((*it2).second)>ret)
1347 ret=fabs((*it2).second);