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 #if __cplusplus >= 201103L
107 restartUsing(src,target);
108 if(m.size()!=target->getNumberOfTuplesExpected())
110 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : input matrix has " << m.size() << " rows whereas there are " << target->getNumberOfTuplesExpected() << " expected !";
111 throw INTERP_KERNEL::Exception(oss.str());
113 auto srcNbElem(src->getNumberOfTuplesExpected());
118 auto idToTest(it2.first);
119 if(idToTest<0 || idToTest>=srcNbElem)
121 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : presence of elt #" << idToTest << " ! not in [0," << srcNbElem << ") !";
122 throw INTERP_KERNEL::Exception(oss.str());
127 _deno_multiply.clear();
128 _deno_multiply.resize(_matrix.size());
129 _deno_reverse_multiply.clear();
130 _deno_reverse_multiply.resize(srcNbElem);
132 throw INTERP_KERNEL::Exception("Breaking news : 10% off for C++11 compiler :)");
136 int MEDCouplingRemapper::prepareInterpKernelOnly()
138 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
146 // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
147 // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
149 // } MEDCouplingMeshType;
151 switch(meshInterpType)
153 case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
154 case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
155 case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
156 case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
157 case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
158 case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
159 case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
160 case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
161 case 85: // UNSTRUCTURED - UNSTRUCTURED
162 return prepareInterpKernelOnlyUU();
163 case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
164 case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
165 case 87: // UNSTRUCTURED - CARTESIAN
166 return prepareInterpKernelOnlyUC();
167 case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
168 case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
169 case 117: // CARTESIAN - UNSTRUCTURED
170 return prepareInterpKernelOnlyCU();
171 case 119: // CARTESIAN - CARTESIAN
172 return prepareInterpKernelOnlyCC();
173 case 136: // EXTRUDED - EXTRUDED
174 return prepareInterpKernelOnlyEE();
176 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
180 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
182 std::string srcm,trgm,method;
183 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
184 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
187 return prepareNotInterpKernelOnlyGaussGauss();
190 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
191 throw INTERP_KERNEL::Exception(oss.str().c_str());
197 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
198 * 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.
200 * \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.
201 * \param [in/out] targetField the destination field with the allocated array in which all tuples will be overwritten.
202 * \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.
203 * 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
204 * 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
205 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
209 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue)
211 if(!srcField || !targetField)
212 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transfer : input field must be both not NULL !");
213 transferUnderground(srcField,targetField,true,dftValue);
217 * This method is equivalent to MEDCoupling::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
218 * 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
220 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
222 * \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.
223 * \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.
225 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField)
227 if(!srcField || !targetField)
228 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : input field must be both not NULL !");
229 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
232 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue)
234 if(!srcField || !targetField)
235 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::reverseTransfer : input fields must be both not NULL !");
237 targetField->checkConsistencyLight();
238 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
239 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
240 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
241 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
242 if(srcField->getNature()!=targetField->getNature())
243 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
244 if(targetField->getNumberOfTuplesExpected()!=_target_ft->getNumberOfTuplesExpected())
246 std::ostringstream oss;
247 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();
248 oss << " ! It appears that the target support is not the same between the prepare and the transfer !";
249 throw INTERP_KERNEL::Exception(oss.str().c_str());
251 DataArrayDouble *array(srcField->getArray());
252 int trgNbOfCompo=targetField->getNumberOfComponents();
255 srcField->checkConsistencyLight();
256 if(trgNbOfCompo!=srcField->getNumberOfTuplesExpected())
257 throw INTERP_KERNEL::Exception("Number of components mismatch !");
261 MCAuto<DataArrayDouble > tmp(DataArrayDouble::New());
262 tmp->alloc(srcField->getNumberOfTuplesExpected(),trgNbOfCompo);
263 srcField->setArray(tmp);
265 computeDeno(srcField->getNature(),srcField,targetField);
266 double *resPointer(srcField->getArray()->getPointer());
267 const double *inputPointer=targetField->getArray()->getConstPointer();
268 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
272 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
273 * If mesh of \b srcField does not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
275 * \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.
276 * \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.
277 * 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
278 * 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
279 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
280 * \return destination field to be deallocated by the caller.
284 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue)
288 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input srcField is NULL !");
289 srcField->checkConsistencyLight();
290 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
291 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
292 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
293 ret->setNature(srcField->getNature());
294 transfer(srcField,ret,dftValue);
295 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
299 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue)
302 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input targetField is NULL !");
303 targetField->checkConsistencyLight();
305 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
306 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
307 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
308 ret->setNature(targetField->getNature());
309 reverseTransfer(ret,targetField,dftValue);
310 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
315 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
316 * is here only for automatic CORBA generators.
318 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
320 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
324 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
325 * is here only for automatic CORBA generators.
327 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
329 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
333 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
334 * is here only for automatic CORBA generators.
336 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
338 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
342 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or preferred.
343 * If interpolation matrix policy is :
345 * - 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
346 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
347 * If not, the \b not only INTERP_KERNEL method will be attempt.
349 * - 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
350 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
351 * If not, the INTERP_KERNEL only method will be attempt.
353 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
355 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
357 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
359 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
361 return _interp_matrix_pol;
365 * 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
366 * 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
367 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
369 * If interpolation matrix policy is :
371 * - 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
372 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
373 * If not, the \b not only INTERP_KERNEL method will be attempt.
375 * - 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
376 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
377 * If not, the INTERP_KERNEL only method will be attempt.
379 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
381 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
383 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c MEDCoupling::InterpolationMatrixPolicy
384 * for automatic generation of CORBA component.
386 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
388 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
390 switch(newInterpMatPol)
393 _interp_matrix_pol=IK_ONLY_PREFERED;
396 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
399 _interp_matrix_pol=IK_ONLY_FORCED;
402 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
405 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 !");
409 int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
411 const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
412 const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_target_ft->getMesh());
413 std::string srcMeth,trgMeth;
414 std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth));
415 const int srcMeshDim=src_mesh->getMeshDimension();
418 srcSpaceDim=src_mesh->getSpaceDimension();
419 const int trgMeshDim=target_mesh->getMeshDimension();
422 trgSpaceDim=target_mesh->getSpaceDimension();
423 if(trgSpaceDim!=srcSpaceDim)
424 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
425 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
427 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
429 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
430 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
431 INTERP_KERNEL::Interpolation1D interpolation(*this);
432 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
434 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
436 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
437 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
438 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
439 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
441 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
443 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
444 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
445 INTERP_KERNEL::Interpolation2D interpolation(*this);
446 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
448 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
450 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
451 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
452 INTERP_KERNEL::Interpolation3D interpolation(*this);
453 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
455 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
457 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
458 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
459 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
460 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
462 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
464 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
465 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
466 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
467 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
468 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
469 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
471 else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
473 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
474 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! Select PointLocator as intersection type !");
475 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
476 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
477 INTERP_KERNEL::Interpolation1D0D interpolation(*this);
478 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
480 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
482 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
483 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
484 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
485 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
486 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
487 std::vector<std::map<int,double> > matrixTmp;
488 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
489 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
490 ReverseMatrix(matrixTmp,nbCols,_matrix);
491 nbCols=matrixTmp.size();
493 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
495 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
497 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
498 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
499 INTERP_KERNEL::Interpolation2D interpolation(*this);
500 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
504 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
505 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
506 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
507 std::vector<std::map<int,double> > matrixTmp;
508 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
509 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
510 ReverseMatrix(matrixTmp,nbCols,_matrix);
511 nbCols=matrixTmp.size();
512 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
513 if(!duplicateFaces.empty())
515 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
516 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
518 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
519 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
525 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
527 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
529 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
530 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
531 INTERP_KERNEL::Interpolation2D interpolation(*this);
532 std::vector<std::map<int,double> > matrixTmp;
533 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
534 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
535 ReverseMatrix(matrixTmp,nbCols,_matrix);
536 nbCols=matrixTmp.size();
540 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
541 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
542 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
543 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
544 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
545 if(!duplicateFaces.empty())
547 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
548 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
550 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
551 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
557 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
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 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
563 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
564 if(!duplicateFaces.empty())
566 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
567 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
569 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
570 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
575 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
577 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
579 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
580 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
581 INTERP_KERNEL::Interpolation3D interpolation(*this);
582 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
586 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
587 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
588 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
589 std::vector<std::map<int,double> > matrixTmp;
590 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
591 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
592 ReverseMatrix(matrixTmp,nbCols,_matrix);
593 nbCols=matrixTmp.size();
594 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
595 if(!duplicateFaces.empty())
597 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
598 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
600 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
601 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
607 else if(trgMeshDim==-1)
609 if(srcMeshDim==2 && srcSpaceDim==2)
611 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
612 INTERP_KERNEL::Interpolation2D interpolation(*this);
613 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
615 else if(srcMeshDim==3 && srcSpaceDim==3)
617 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
618 INTERP_KERNEL::Interpolation3D interpolation(*this);
619 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
621 else if(srcMeshDim==2 && srcSpaceDim==3)
623 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
624 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
625 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
628 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
630 else if(srcMeshDim==-1)
632 if(trgMeshDim==2 && trgSpaceDim==2)
634 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
635 INTERP_KERNEL::Interpolation2D interpolation(*this);
636 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
638 else if(trgMeshDim==3 && trgSpaceDim==3)
640 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
641 INTERP_KERNEL::Interpolation3D interpolation(*this);
642 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
644 else if(trgMeshDim==2 && trgSpaceDim==3)
646 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
647 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
648 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
651 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
654 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
655 _deno_multiply.clear();
656 _deno_multiply.resize(_matrix.size());
657 _deno_reverse_multiply.clear();
658 _deno_reverse_multiply.resize(nbCols);
663 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
665 std::string srcMeth,trgMeth;
666 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
667 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
668 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
670 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
671 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
672 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
673 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
674 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
675 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
676 std::vector<std::map<int,double> > matrix2D;
677 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
678 MEDCouplingUMesh *s1D,*t1D;
680 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
681 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
682 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
683 std::vector<std::map<int,double> > matrix1D;
684 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
685 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
686 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
687 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
690 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
691 target_mesh->getMesh3DIds()->getConstPointer());
693 _deno_multiply.clear();
694 _deno_multiply.resize(_matrix.size());
695 _deno_reverse_multiply.clear();
696 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
701 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
703 std::string srcMeth,trgMeth;
704 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
705 if(methodCpp!="P0P0")
706 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
707 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
708 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
709 const int srcMeshDim=src_mesh->getMeshDimension();
710 const int srcSpceDim=src_mesh->getSpaceDimension();
711 const int trgMeshDim=target_mesh->getMeshDimension();
712 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
713 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 !");
714 std::vector<std::map<int,double> > res;
719 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
720 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
721 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
722 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
727 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
728 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
729 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
730 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
735 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
736 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
737 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
738 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
742 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
744 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
745 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
747 _deno_multiply.clear();
748 _deno_multiply.resize(_matrix.size());
749 _deno_reverse_multiply.clear();
750 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
755 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
757 std::string srcMeth,trgMeth;
758 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
759 if(methodCpp!="P0P0")
760 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
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::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 !");
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 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
814 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
815 const int srcMeshDim=src_mesh->getMeshDimension();
816 const int trgMeshDim=target_mesh->getMeshDimension();
817 if(trgMeshDim!=srcMeshDim)
818 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
823 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
824 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
825 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
826 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
831 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
832 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
833 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
834 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
839 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
840 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
841 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
842 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
846 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
848 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
850 _deno_multiply.clear();
851 _deno_multiply.resize(_matrix.size());
852 _deno_reverse_multiply.clear();
853 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
858 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
860 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
861 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 !");
862 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
863 const double *trgLocPtr=trgLoc->begin();
864 int trgSpaceDim=trgLoc->getNumberOfComponents();
865 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
866 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
868 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
869 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
870 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
871 throw INTERP_KERNEL::Exception(oss.str().c_str());
873 const int *srcOffsetArrPtr=srcOffsetArr->begin();
874 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
875 const double *srcLocPtr=srcLoc->begin();
876 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
877 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
878 _matrix.resize(trgNbOfGaussPts);
879 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
880 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
881 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
882 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
883 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
885 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
886 int srcCellId=elts[eltsIndex[*trgId]];
887 double dist=std::numeric_limits<double>::max();
889 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
891 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
893 for(int i=0;i<trgSpaceDim;i++)
894 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
896 { dist=tmp; srcEntry=srcId; }
898 _matrix[*trgId][srcEntry]=1.;
900 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
902 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
903 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
904 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
905 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
906 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
907 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
909 _deno_multiply.clear();
910 _deno_multiply.resize(_matrix.size());
911 _deno_reverse_multiply.clear();
912 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
918 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
919 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
921 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
923 if(method=="GAUSSGAUSS")
925 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
926 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
927 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
928 throw INTERP_KERNEL::Exception(oss.str().c_str());
932 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
933 * 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
934 * only INTERP_KERNEL method should be applied.
936 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
938 std::string srcm,trgm,method;
939 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
940 switch(_interp_matrix_pol)
942 case IK_ONLY_PREFERED:
946 std::string tmp1,tmp2;
947 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
950 catch(INTERP_KERNEL::Exception& /*e*/)
955 case NOT_IK_ONLY_PREFERED:
959 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
962 catch(INTERP_KERNEL::Exception& /*e*/)
969 case NOT_IK_ONLY_FORCED:
972 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
976 void MEDCouplingRemapper::updateTime() const
980 void MEDCouplingRemapper::checkPrepare() const
982 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
984 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
985 if(!s->getMesh() || !t->getMesh())
986 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
990 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
991 * This method returns 3 information (2 in output parameters and 1 in return).
993 * \param [out] srcMeth the string code of the discretization of source field template
994 * \param [out] trgMeth the string code of the discretization of target field template
995 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
997 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
999 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1001 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1002 if(!s->getMesh() || !t->getMesh())
1003 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1004 srcMeth=_src_ft->getDiscretization()->getRepr();
1005 trgMeth=_target_ft->getDiscretization()->getRepr();
1006 return BuildMethodFrom(srcMeth,trgMeth);
1009 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1011 std::string method(meth1); method+=meth2;
1015 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1017 if(!srcMesh || !targetMesh)
1018 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1019 std::string srcMethod,targetMethod;
1020 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1021 src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1022 src->setMesh(srcMesh);
1023 target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1024 target->setMesh(targetMesh);
1027 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1031 if(matrixSuppression)
1034 _deno_multiply.clear();
1035 _deno_reverse_multiply.clear();
1039 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1042 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1043 if(!src->getMesh() || !target->getMesh())
1044 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1046 _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1047 _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1050 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1052 if(!srcField || !targetField)
1053 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1054 srcField->checkConsistencyLight();
1056 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1057 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1058 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1059 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1060 if(srcField->getNature()!=targetField->getNature())
1061 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1062 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1064 std::ostringstream oss;
1065 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();
1066 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1067 throw INTERP_KERNEL::Exception(oss.str().c_str());
1069 DataArrayDouble *array(targetField->getArray());
1070 int srcNbOfCompo(srcField->getNumberOfComponents());
1073 targetField->checkConsistencyLight();
1074 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1075 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1080 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1081 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1082 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1083 targetField->setArray(tmp);
1085 computeDeno(srcField->getNature(),srcField,targetField);
1086 double *resPointer(targetField->getArray()->getPointer());
1087 const double *inputPointer(srcField->getArray()->getConstPointer());
1088 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1091 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1094 return computeDenoFromScratch(nat,srcField,trgField);
1095 else if(nat!=_nature_of_deno)
1096 return computeDenoFromScratch(nat,srcField,trgField);
1097 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1098 return computeDenoFromScratch(nat,srcField,trgField);
1101 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1103 _nature_of_deno=nat;
1104 _time_deno_update=getTimeOfThis();
1105 switch(_nature_of_deno)
1107 case IntensiveMaximum:
1109 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1112 case ExtensiveMaximum:
1114 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1115 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1116 const double *denoPtr=deno->getArray()->getConstPointer();
1117 const double *denoRPtr=denoR->getArray()->getConstPointer();
1118 if(trgField->getMesh()->getMeshDimension()==-1)
1120 double *denoRPtr2=denoR->getArray()->getPointer();
1121 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1123 if(srcField->getMesh()->getMeshDimension()==-1)
1125 double *denoPtr2=deno->getArray()->getPointer();
1126 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1129 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1130 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1132 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1133 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1139 case ExtensiveConservation:
1141 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1144 case IntensiveConservation:
1146 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1147 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1148 const double *denoPtr=deno->getArray()->getConstPointer();
1149 const double *denoRPtr=denoR->getArray()->getConstPointer();
1150 if(trgField->getMesh()->getMeshDimension()==-1)
1152 double *denoRPtr2=denoR->getArray()->getPointer();
1153 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1155 if(srcField->getMesh()->getMeshDimension()==-1)
1157 double *denoPtr2=deno->getArray()->getPointer();
1158 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1161 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1162 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1164 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1165 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1172 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1176 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1179 double *tmp=new double[inputNbOfCompo];
1180 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1182 if((*iter1).empty())
1185 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1189 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1190 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1191 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1193 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1194 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1200 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1202 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1204 double *tmp=new double[inputNbOfCompo];
1205 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1206 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1208 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1210 isReached[(*iter2).first]=true;
1211 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1212 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1217 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1219 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1222 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1224 matOut.resize(nbColsMatIn);
1226 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1227 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1228 matOut[(*iter2).first][id]=(*iter2).second;
1231 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1232 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1234 std::map<int,double> values;
1236 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1239 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1241 sum+=(*iter2).second;
1242 values[(*iter2).first]+=(*iter2).second;
1244 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1245 deno[idx][(*iter2).first]=sum;
1248 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1250 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1251 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1255 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1256 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1258 std::map<int,double> values;
1260 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1263 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1265 sum+=(*iter2).second;
1266 values[(*iter2).first]+=(*iter2).second;
1268 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1269 denoReverse[(*iter2).first][idx]=sum;
1272 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1274 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1275 deno[idx][(*iter2).first]=values[(*iter2).first];
1279 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1280 const std::vector< std::map<int,double> >& m2D,
1281 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1282 const int *corrCellIdTrg)
1284 int nbOf2DCellsTrg=m2D.size();
1285 int nbOf1DCellsTrg=m1D.size();
1286 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1287 _matrix.resize(nbOf3DCellsTrg);
1289 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1291 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1294 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1296 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1298 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1305 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1308 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1310 std::cout << "Target Cell # " << id << " : ";
1311 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1312 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1313 std::cout << std::endl;
1317 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1323 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1325 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1327 return (int)_deno_reverse_multiply.size();
1331 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1332 * If not the behaviour is unpredictable.
1333 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1334 * 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.
1335 * 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
1338 * \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.
1339 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1340 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1342 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1345 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1347 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1349 std::map<int,double>& rowNew=matrixNew[i];
1350 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1352 if(fabs((*it2).second)>maxValAbs)
1353 rowNew[(*it2).first]=(*it2).second;
1364 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1365 * If not the behaviour is unpredictable.
1366 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1367 * 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.
1368 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1369 * 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
1372 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1373 * \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
1374 * that all coefficients are null.
1375 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1377 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1379 double maxVal=getMaxValueInCrudeMatrix();
1382 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1386 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1387 * If not the behaviour is unpredictable.
1388 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1389 * The returned value is positive.
1391 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1394 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1395 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1396 if(fabs((*it2).second)>ret)
1397 ret=fabs((*it2).second);