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 MCAuto<MEDCouplingFieldTemplate> src,target;
75 BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
76 return prepareEx(src,target);
80 * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
81 * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
83 * \param [in] src is the field template source side.
84 * \param [in] target is the field template target side.
86 * \sa MEDCouplingRemapper::prepare
88 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
90 restartUsing(src,target);
91 if(isInterpKernelOnlyOrNotOnly())
92 return prepareInterpKernelOnly();
94 return prepareNotInterpKernelOnly();
97 void MEDCouplingRemapper::setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector<std::map<int,double> >& m)
99 MCAuto<MEDCouplingFieldTemplate> src,target;
100 BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
101 setCrudeMatrixEx(src,target,m);
104 void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector<std::map<int,double> >& m)
106 restartUsing(src,target);
107 if(m.size()!=target->getNumberOfTuplesExpected())
109 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : input matrix has " << m.size() << " rows whereas there are " << target->getNumberOfTuplesExpected() << " expected !";
110 throw INTERP_KERNEL::Exception(oss.str());
112 auto srcNbElem(src->getNumberOfTuplesExpected());
117 auto idToTest(it2.first);
118 if(idToTest<0 || idToTest>=srcNbElem)
120 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : presence of elt #" << idToTest << " ! not in [0," << srcNbElem << ") !";
121 throw INTERP_KERNEL::Exception(oss.str());
126 _deno_multiply.clear();
127 _deno_multiply.resize(_matrix.size());
128 _deno_reverse_multiply.clear();
129 _deno_reverse_multiply.resize(srcNbElem);
132 int MEDCouplingRemapper::prepareInterpKernelOnly()
134 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
142 // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
143 // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
145 // } MEDCouplingMeshType;
147 switch(meshInterpType)
149 case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
150 case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
151 case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
152 case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
153 case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
154 case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
155 case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
156 case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
157 case 85: // UNSTRUCTURED - UNSTRUCTURED
158 return prepareInterpKernelOnlyUU();
159 case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
160 case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
161 case 87: // UNSTRUCTURED - CARTESIAN
162 return prepareInterpKernelOnlyUC();
163 case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
164 case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
165 case 117: // CARTESIAN - UNSTRUCTURED
166 return prepareInterpKernelOnlyCU();
167 case 119: // CARTESIAN - CARTESIAN
168 return prepareInterpKernelOnlyCC();
169 case 136: // EXTRUDED - EXTRUDED
170 return prepareInterpKernelOnlyEE();
172 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
176 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
178 std::string srcm,trgm,method;
179 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
180 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
183 return prepareNotInterpKernelOnlyGaussGauss();
186 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
187 throw INTERP_KERNEL::Exception(oss.str().c_str());
193 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
194 * 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.
196 * \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.
197 * \param [in/out] targetField the destination field with the allocated array in which all tuples will be overwritten.
198 * \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.
199 * 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
200 * 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
201 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
205 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue)
207 if(!srcField || !targetField)
208 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transfer : input field must be both not NULL !");
209 transferUnderground(srcField,targetField,true,dftValue);
213 * This method is equivalent to MEDCoupling::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
214 * 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
216 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
218 * \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.
219 * \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.
221 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField)
223 if(!srcField || !targetField)
224 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : input field must be both not NULL !");
225 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
228 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue)
230 if(!srcField || !targetField)
231 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::reverseTransfer : input fields must be both not NULL !");
233 targetField->checkConsistencyLight();
234 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
235 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
236 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
237 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
238 if(srcField->getNature()!=targetField->getNature())
239 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
240 if(targetField->getNumberOfTuplesExpected()!=_target_ft->getNumberOfTuplesExpected())
242 std::ostringstream oss;
243 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();
244 oss << " ! It appears that the target support is not the same between the prepare and the transfer !";
245 throw INTERP_KERNEL::Exception(oss.str().c_str());
247 DataArrayDouble *array(srcField->getArray());
248 int trgNbOfCompo=targetField->getNumberOfComponents();
251 srcField->checkConsistencyLight();
252 if(trgNbOfCompo!=srcField->getNumberOfTuplesExpected())
253 throw INTERP_KERNEL::Exception("Number of components mismatch !");
257 MCAuto<DataArrayDouble > tmp(DataArrayDouble::New());
258 tmp->alloc(srcField->getNumberOfTuplesExpected(),trgNbOfCompo);
259 srcField->setArray(tmp);
261 computeDeno(srcField->getNature(),srcField,targetField);
262 double *resPointer(srcField->getArray()->getPointer());
263 const double *inputPointer=targetField->getArray()->getConstPointer();
264 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
268 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
269 * If mesh of \b srcField does not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
271 * \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.
272 * \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.
273 * 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
274 * 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
275 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
276 * \return destination field to be deallocated by the caller.
280 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue)
284 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input srcField is NULL !");
285 srcField->checkConsistencyLight();
286 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
287 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
288 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
289 ret->setNature(srcField->getNature());
290 transfer(srcField,ret,dftValue);
291 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
295 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue)
298 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input targetField is NULL !");
299 targetField->checkConsistencyLight();
301 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
302 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
303 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
304 ret->setNature(targetField->getNature());
305 reverseTransfer(ret,targetField,dftValue);
306 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
311 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
312 * is here only for automatic CORBA generators.
314 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
316 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
320 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
321 * is here only for automatic CORBA generators.
323 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
325 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
329 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
330 * is here only for automatic CORBA generators.
332 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
334 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
338 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or preferred.
339 * If interpolation matrix policy is :
341 * - 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
342 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
343 * If not, the \b not only INTERP_KERNEL method will be attempt.
345 * - 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
346 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
347 * If not, the INTERP_KERNEL only method will be attempt.
349 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
351 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
353 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
355 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
357 return _interp_matrix_pol;
361 * 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
362 * 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
363 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
365 * If interpolation matrix policy is :
367 * - 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
368 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
369 * If not, the \b not only INTERP_KERNEL method will be attempt.
371 * - 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
372 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
373 * If not, the INTERP_KERNEL only method will be attempt.
375 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
377 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
379 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c MEDCoupling::InterpolationMatrixPolicy
380 * for automatic generation of CORBA component.
382 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
384 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
386 switch(newInterpMatPol)
389 _interp_matrix_pol=IK_ONLY_PREFERED;
392 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
395 _interp_matrix_pol=IK_ONLY_FORCED;
398 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
401 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 !");
405 int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
407 const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
408 const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_target_ft->getMesh());
409 std::string srcMeth,trgMeth;
410 std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth));
411 const int srcMeshDim=src_mesh->getMeshDimension();
414 srcSpaceDim=src_mesh->getSpaceDimension();
415 const int trgMeshDim=target_mesh->getMeshDimension();
418 trgSpaceDim=target_mesh->getSpaceDimension();
419 if(trgSpaceDim!=srcSpaceDim)
420 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
421 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
423 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
425 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
426 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
427 INTERP_KERNEL::Interpolation1D interpolation(*this);
428 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
430 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
432 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
433 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
434 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
435 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
437 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
439 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
440 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
441 INTERP_KERNEL::Interpolation2D interpolation(*this);
442 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
444 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
446 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
447 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
448 INTERP_KERNEL::Interpolation3D interpolation(*this);
449 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
451 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
453 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
454 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
455 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
456 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
458 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
460 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
461 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
462 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
463 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
464 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
465 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
467 else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
469 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
470 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! Select PointLocator as intersection type !");
471 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
472 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
473 INTERP_KERNEL::Interpolation1D0D interpolation(*this);
474 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
476 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
478 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
479 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
480 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
481 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
482 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
483 std::vector<std::map<int,double> > matrixTmp;
484 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
485 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
486 ReverseMatrix(matrixTmp,nbCols,_matrix);
487 nbCols=matrixTmp.size();
489 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
491 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
493 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
494 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
495 INTERP_KERNEL::Interpolation2D interpolation(*this);
496 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
500 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
501 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
502 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
503 std::vector<std::map<int,double> > matrixTmp;
504 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
505 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
506 ReverseMatrix(matrixTmp,nbCols,_matrix);
507 nbCols=matrixTmp.size();
508 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
509 if(!duplicateFaces.empty())
511 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
512 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
514 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
515 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
521 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
523 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
525 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
526 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
527 INTERP_KERNEL::Interpolation2D interpolation(*this);
528 std::vector<std::map<int,double> > matrixTmp;
529 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
530 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
531 ReverseMatrix(matrixTmp,nbCols,_matrix);
532 nbCols=matrixTmp.size();
536 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
537 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
538 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
539 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
540 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
541 if(!duplicateFaces.empty())
543 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
544 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
546 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
547 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
553 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
555 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
556 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
557 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
558 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
559 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
560 if(!duplicateFaces.empty())
562 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
563 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
565 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
566 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
571 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
573 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
575 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
576 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
577 INTERP_KERNEL::Interpolation3D interpolation(*this);
578 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
582 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
583 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
584 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
585 std::vector<std::map<int,double> > matrixTmp;
586 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
587 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
588 ReverseMatrix(matrixTmp,nbCols,_matrix);
589 nbCols=matrixTmp.size();
590 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
591 if(!duplicateFaces.empty())
593 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
594 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
596 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
597 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
603 else if(trgMeshDim==-1)
605 if(srcMeshDim==2 && srcSpaceDim==2)
607 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
608 INTERP_KERNEL::Interpolation2D interpolation(*this);
609 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
611 else if(srcMeshDim==3 && srcSpaceDim==3)
613 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
614 INTERP_KERNEL::Interpolation3D interpolation(*this);
615 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
617 else if(srcMeshDim==2 && srcSpaceDim==3)
619 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
620 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
621 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
624 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
626 else if(srcMeshDim==-1)
628 if(trgMeshDim==2 && trgSpaceDim==2)
630 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
631 INTERP_KERNEL::Interpolation2D interpolation(*this);
632 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
634 else if(trgMeshDim==3 && trgSpaceDim==3)
636 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
637 INTERP_KERNEL::Interpolation3D interpolation(*this);
638 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
640 else if(trgMeshDim==2 && trgSpaceDim==3)
642 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
643 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
644 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
647 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
650 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
651 _deno_multiply.clear();
652 _deno_multiply.resize(_matrix.size());
653 _deno_reverse_multiply.clear();
654 _deno_reverse_multiply.resize(nbCols);
659 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
661 std::string srcMeth,trgMeth;
662 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
663 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
664 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
666 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
667 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
668 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
669 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
670 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
671 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
672 std::vector<std::map<int,double> > matrix2D;
673 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
674 MEDCouplingUMesh *s1D,*t1D;
676 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
677 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
678 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
679 std::vector<std::map<int,double> > matrix1D;
680 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
681 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
682 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
683 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
686 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
687 target_mesh->getMesh3DIds()->getConstPointer());
689 _deno_multiply.clear();
690 _deno_multiply.resize(_matrix.size());
691 _deno_reverse_multiply.clear();
692 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
697 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
699 std::string srcMeth,trgMeth;
700 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
701 if(methodCpp!="P0P0")
702 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only P0P0 interpolation supported for the moment !");
703 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
704 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
705 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
706 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
707 const int srcMeshDim=src_mesh->getMeshDimension();
708 const int srcSpceDim=src_mesh->getSpaceDimension();
709 const int trgMeshDim=target_mesh->getMeshDimension();
710 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
711 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: space dimension of unstructured source mesh should be equal to mesh dimension of unstructured source mesh, and should also be equal to target cartesian dimension!");
712 std::vector<std::map<int,double> > res;
717 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
718 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
719 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
720 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
725 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
726 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
727 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
728 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
733 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
734 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
735 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
736 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
740 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
742 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
743 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
745 _deno_multiply.clear();
746 _deno_multiply.resize(_matrix.size());
747 _deno_reverse_multiply.clear();
748 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
753 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
755 std::string srcMeth,trgMeth;
756 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
757 if(methodCpp!="P0P0")
758 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
759 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
760 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
761 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
762 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
763 const int srcMeshDim=src_mesh->getMeshDimension();
764 const int trgMeshDim=target_mesh->getMeshDimension();
765 const int trgSpceDim=target_mesh->getSpaceDimension();
766 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
767 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: space dimension of unstructured target mesh should be equal to mesh dimension of unstructured target mesh, and should also be equal to source cartesian dimension!");
772 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
773 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
774 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
775 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
780 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
781 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
782 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
783 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
788 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
789 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
790 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
791 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
795 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
797 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
799 _deno_multiply.clear();
800 _deno_multiply.resize(_matrix.size());
801 _deno_reverse_multiply.clear();
802 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
807 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
809 std::string srcMeth,trgMeth;
810 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
811 if(methodCpp!="P0P0")
812 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
813 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
814 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
815 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
816 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
817 const int srcMeshDim=src_mesh->getMeshDimension();
818 const int trgMeshDim=target_mesh->getMeshDimension();
819 if(trgMeshDim!=srcMeshDim)
820 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
825 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
826 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
827 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
828 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
833 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
834 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
835 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
836 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
841 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
842 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
843 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
844 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
848 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
850 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
852 _deno_multiply.clear();
853 _deno_multiply.resize(_matrix.size());
854 _deno_reverse_multiply.clear();
855 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
860 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
862 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
863 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 !");
864 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
865 const double *trgLocPtr=trgLoc->begin();
866 int trgSpaceDim=trgLoc->getNumberOfComponents();
867 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
868 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
870 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
871 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
872 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
873 throw INTERP_KERNEL::Exception(oss.str().c_str());
875 const int *srcOffsetArrPtr=srcOffsetArr->begin();
876 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
877 const double *srcLocPtr=srcLoc->begin();
878 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
879 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
880 _matrix.resize(trgNbOfGaussPts);
881 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
882 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
883 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
884 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
885 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
887 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
888 int srcCellId=elts[eltsIndex[*trgId]];
889 double dist=std::numeric_limits<double>::max();
891 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
893 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
895 for(int i=0;i<trgSpaceDim;i++)
896 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
898 { dist=tmp; srcEntry=srcId; }
900 _matrix[*trgId][srcEntry]=1.;
902 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
904 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
905 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
906 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
907 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
908 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
909 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
911 _deno_multiply.clear();
912 _deno_multiply.resize(_matrix.size());
913 _deno_reverse_multiply.clear();
914 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
920 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
921 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
923 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
925 if(method=="GAUSSGAUSS")
927 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
928 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
929 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
930 throw INTERP_KERNEL::Exception(oss.str().c_str());
934 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
935 * 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
936 * only INTERP_KERNEL method should be applied.
938 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
940 std::string srcm,trgm,method;
941 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
942 switch(_interp_matrix_pol)
944 case IK_ONLY_PREFERED:
948 std::string tmp1,tmp2;
949 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
952 catch(INTERP_KERNEL::Exception& /*e*/)
957 case NOT_IK_ONLY_PREFERED:
961 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
964 catch(INTERP_KERNEL::Exception& /*e*/)
971 case NOT_IK_ONLY_FORCED:
974 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
978 void MEDCouplingRemapper::updateTime() const
982 void MEDCouplingRemapper::checkPrepare() const
984 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
986 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
987 if(!s->getMesh() || !t->getMesh())
988 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
992 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
993 * This method returns 3 information (2 in output parameters and 1 in return).
995 * \param [out] srcMeth the string code of the discretization of source field template
996 * \param [out] trgMeth the string code of the discretization of target field template
997 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
999 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
1001 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1003 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1004 if(!s->getMesh() || !t->getMesh())
1005 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1006 srcMeth=_src_ft->getDiscretization()->getRepr();
1007 trgMeth=_target_ft->getDiscretization()->getRepr();
1008 return BuildMethodFrom(srcMeth,trgMeth);
1011 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1013 std::string method(meth1); method+=meth2;
1017 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1019 if(!srcMesh || !targetMesh)
1020 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1021 std::string srcMethod,targetMethod;
1022 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1023 src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1024 src->setMesh(srcMesh);
1025 target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1026 target->setMesh(targetMesh);
1029 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1033 if(matrixSuppression)
1036 _deno_multiply.clear();
1037 _deno_reverse_multiply.clear();
1041 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1044 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1045 if(!src->getMesh() || !target->getMesh())
1046 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1048 _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1049 _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1052 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1054 if(!srcField || !targetField)
1055 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1056 srcField->checkConsistencyLight();
1058 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1059 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1060 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1061 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1062 if(srcField->getNature()!=targetField->getNature())
1063 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1064 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1066 std::ostringstream oss;
1067 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();
1068 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1069 throw INTERP_KERNEL::Exception(oss.str().c_str());
1071 DataArrayDouble *array(targetField->getArray());
1072 int srcNbOfCompo(srcField->getNumberOfComponents());
1075 targetField->checkConsistencyLight();
1076 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1077 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1082 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1083 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1084 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1085 targetField->setArray(tmp);
1087 computeDeno(srcField->getNature(),srcField,targetField);
1088 double *resPointer(targetField->getArray()->getPointer());
1089 const double *inputPointer(srcField->getArray()->getConstPointer());
1090 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1093 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1096 return computeDenoFromScratch(nat,srcField,trgField);
1097 else if(nat!=_nature_of_deno)
1098 return computeDenoFromScratch(nat,srcField,trgField);
1099 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1100 return computeDenoFromScratch(nat,srcField,trgField);
1103 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1105 _nature_of_deno=nat;
1106 _time_deno_update=getTimeOfThis();
1107 switch(_nature_of_deno)
1109 case IntensiveMaximum:
1111 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1114 case ExtensiveMaximum:
1116 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1117 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1118 const double *denoPtr=deno->getArray()->getConstPointer();
1119 const double *denoRPtr=denoR->getArray()->getConstPointer();
1120 if(trgField->getMesh()->getMeshDimension()==-1)
1122 double *denoRPtr2=denoR->getArray()->getPointer();
1123 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1125 if(srcField->getMesh()->getMeshDimension()==-1)
1127 double *denoPtr2=deno->getArray()->getPointer();
1128 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1131 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1132 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1134 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1135 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1141 case ExtensiveConservation:
1143 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1146 case IntensiveConservation:
1148 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1149 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1150 const double *denoPtr=deno->getArray()->getConstPointer();
1151 const double *denoRPtr=denoR->getArray()->getConstPointer();
1152 if(trgField->getMesh()->getMeshDimension()==-1)
1154 double *denoRPtr2=denoR->getArray()->getPointer();
1155 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1157 if(srcField->getMesh()->getMeshDimension()==-1)
1159 double *denoPtr2=deno->getArray()->getPointer();
1160 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1163 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1164 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1166 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1167 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1174 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1178 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1181 double *tmp=new double[inputNbOfCompo];
1182 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1184 if((*iter1).empty())
1187 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1191 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1192 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1193 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1195 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1196 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1202 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1204 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1206 double *tmp=new double[inputNbOfCompo];
1207 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1208 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1210 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1212 isReached[(*iter2).first]=true;
1213 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1214 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1219 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1221 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1224 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1226 matOut.resize(nbColsMatIn);
1228 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1229 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1230 matOut[(*iter2).first][id]=(*iter2).second;
1233 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1234 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1236 std::map<int,double> values;
1238 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1241 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1243 sum+=(*iter2).second;
1244 values[(*iter2).first]+=(*iter2).second;
1246 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1247 deno[idx][(*iter2).first]=sum;
1250 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1252 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1253 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1257 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1258 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1260 std::map<int,double> values;
1262 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1265 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1267 sum+=(*iter2).second;
1268 values[(*iter2).first]+=(*iter2).second;
1270 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1271 denoReverse[(*iter2).first][idx]=sum;
1274 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1276 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1277 deno[idx][(*iter2).first]=values[(*iter2).first];
1281 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1282 const std::vector< std::map<int,double> >& m2D,
1283 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1284 const int *corrCellIdTrg)
1286 int nbOf2DCellsTrg=m2D.size();
1287 int nbOf1DCellsTrg=m1D.size();
1288 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1289 _matrix.resize(nbOf3DCellsTrg);
1291 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1293 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1296 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1298 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1300 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1307 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1310 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1312 std::cout << "Target Cell # " << id << " : ";
1313 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1314 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1315 std::cout << std::endl;
1319 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1325 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1327 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1329 return (int)_deno_reverse_multiply.size();
1333 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1334 * If not the behaviour is unpredictable.
1335 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1336 * 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.
1337 * 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
1340 * \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.
1341 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1342 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1344 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1347 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1349 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1351 std::map<int,double>& rowNew=matrixNew[i];
1352 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1354 if(fabs((*it2).second)>maxValAbs)
1355 rowNew[(*it2).first]=(*it2).second;
1366 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1367 * If not the behaviour is unpredictable.
1368 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1369 * 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.
1370 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1371 * 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
1374 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1375 * \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
1376 * that all coefficients are null.
1377 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1379 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1381 double maxVal=getMaxValueInCrudeMatrix();
1384 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1388 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1389 * If not the behaviour is unpredictable.
1390 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1391 * The returned value is positive.
1393 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1396 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1397 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1398 if(fabs((*it2).second)>ret)
1399 ret=fabs((*it2).second);