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 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
708 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
709 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
710 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
711 const int srcMeshDim=src_mesh->getMeshDimension();
712 const int srcSpceDim=src_mesh->getSpaceDimension();
713 const int trgMeshDim=target_mesh->getMeshDimension();
714 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
715 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!");
716 std::vector<std::map<int,double> > res;
721 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
722 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
723 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
724 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
729 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
730 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
731 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
732 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
737 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
738 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
739 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
740 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
744 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
746 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
747 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
749 _deno_multiply.clear();
750 _deno_multiply.resize(_matrix.size());
751 _deno_reverse_multiply.clear();
752 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
757 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
759 std::string srcMeth,trgMeth;
760 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
761 if(methodCpp!="P0P0")
762 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
763 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
764 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
765 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
766 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
767 const int srcMeshDim=src_mesh->getMeshDimension();
768 const int trgMeshDim=target_mesh->getMeshDimension();
769 const int trgSpceDim=target_mesh->getSpaceDimension();
770 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
771 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!");
776 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
777 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
778 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
779 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
784 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
785 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
786 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
787 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
792 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
793 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
794 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
795 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
799 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
801 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
803 _deno_multiply.clear();
804 _deno_multiply.resize(_matrix.size());
805 _deno_reverse_multiply.clear();
806 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
811 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
813 std::string srcMeth,trgMeth;
814 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
815 if(methodCpp!="P0P0")
816 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
817 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
818 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
819 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
820 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
821 const int srcMeshDim=src_mesh->getMeshDimension();
822 const int trgMeshDim=target_mesh->getMeshDimension();
823 if(trgMeshDim!=srcMeshDim)
824 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
829 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
830 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
831 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
832 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
837 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
838 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
839 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
840 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
845 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
846 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
847 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
848 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
852 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
854 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
856 _deno_multiply.clear();
857 _deno_multiply.resize(_matrix.size());
858 _deno_reverse_multiply.clear();
859 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
864 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
866 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
867 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 !");
868 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
869 const double *trgLocPtr=trgLoc->begin();
870 int trgSpaceDim=trgLoc->getNumberOfComponents();
871 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
872 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
874 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
875 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
876 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
877 throw INTERP_KERNEL::Exception(oss.str().c_str());
879 const int *srcOffsetArrPtr=srcOffsetArr->begin();
880 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
881 const double *srcLocPtr=srcLoc->begin();
882 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
883 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
884 _matrix.resize(trgNbOfGaussPts);
885 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
886 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
887 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
888 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
889 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
891 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
892 int srcCellId=elts[eltsIndex[*trgId]];
893 double dist=std::numeric_limits<double>::max();
895 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
897 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
899 for(int i=0;i<trgSpaceDim;i++)
900 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
902 { dist=tmp; srcEntry=srcId; }
904 _matrix[*trgId][srcEntry]=1.;
906 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
908 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
909 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
910 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
911 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
912 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
913 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
915 _deno_multiply.clear();
916 _deno_multiply.resize(_matrix.size());
917 _deno_reverse_multiply.clear();
918 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
924 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
925 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
927 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
929 if(method=="GAUSSGAUSS")
931 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
932 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
933 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
934 throw INTERP_KERNEL::Exception(oss.str().c_str());
938 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
939 * 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
940 * only INTERP_KERNEL method should be applied.
942 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
944 std::string srcm,trgm,method;
945 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
946 switch(_interp_matrix_pol)
948 case IK_ONLY_PREFERED:
952 std::string tmp1,tmp2;
953 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
956 catch(INTERP_KERNEL::Exception& /*e*/)
961 case NOT_IK_ONLY_PREFERED:
965 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
968 catch(INTERP_KERNEL::Exception& /*e*/)
975 case NOT_IK_ONLY_FORCED:
978 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
982 void MEDCouplingRemapper::updateTime() const
986 void MEDCouplingRemapper::checkPrepare() const
988 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
990 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
991 if(!s->getMesh() || !t->getMesh())
992 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
996 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
997 * This method returns 3 information (2 in output parameters and 1 in return).
999 * \param [out] srcMeth the string code of the discretization of source field template
1000 * \param [out] trgMeth the string code of the discretization of target field template
1001 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
1003 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
1005 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1007 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1008 if(!s->getMesh() || !t->getMesh())
1009 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1010 srcMeth=_src_ft->getDiscretization()->getRepr();
1011 trgMeth=_target_ft->getDiscretization()->getRepr();
1012 return BuildMethodFrom(srcMeth,trgMeth);
1015 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1017 std::string method(meth1); method+=meth2;
1021 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1023 if(!srcMesh || !targetMesh)
1024 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1025 std::string srcMethod,targetMethod;
1026 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1027 src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1028 src->setMesh(srcMesh);
1029 target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1030 target->setMesh(targetMesh);
1033 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1037 if(matrixSuppression)
1040 _deno_multiply.clear();
1041 _deno_reverse_multiply.clear();
1045 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1048 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1049 if(!src->getMesh() || !target->getMesh())
1050 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1052 _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1053 _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1056 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1058 if(!srcField || !targetField)
1059 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1060 srcField->checkConsistencyLight();
1062 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1063 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1064 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1065 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1066 if(srcField->getNature()!=targetField->getNature())
1067 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1068 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1070 std::ostringstream oss;
1071 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();
1072 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1073 throw INTERP_KERNEL::Exception(oss.str().c_str());
1075 DataArrayDouble *array(targetField->getArray());
1076 int srcNbOfCompo(srcField->getNumberOfComponents());
1079 targetField->checkConsistencyLight();
1080 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1081 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1086 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1087 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1088 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1089 targetField->setArray(tmp);
1091 computeDeno(srcField->getNature(),srcField,targetField);
1092 double *resPointer(targetField->getArray()->getPointer());
1093 const double *inputPointer(srcField->getArray()->getConstPointer());
1094 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1097 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1100 return computeDenoFromScratch(nat,srcField,trgField);
1101 else if(nat!=_nature_of_deno)
1102 return computeDenoFromScratch(nat,srcField,trgField);
1103 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1104 return computeDenoFromScratch(nat,srcField,trgField);
1107 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1109 _nature_of_deno=nat;
1110 _time_deno_update=getTimeOfThis();
1111 switch(_nature_of_deno)
1113 case IntensiveMaximum:
1115 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1118 case ExtensiveMaximum:
1120 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1121 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1122 const double *denoPtr=deno->getArray()->getConstPointer();
1123 const double *denoRPtr=denoR->getArray()->getConstPointer();
1124 if(trgField->getMesh()->getMeshDimension()==-1)
1126 double *denoRPtr2=denoR->getArray()->getPointer();
1127 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1129 if(srcField->getMesh()->getMeshDimension()==-1)
1131 double *denoPtr2=deno->getArray()->getPointer();
1132 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1135 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1136 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1138 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1139 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1145 case ExtensiveConservation:
1147 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1150 case IntensiveConservation:
1152 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1153 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1154 const double *denoPtr=deno->getArray()->getConstPointer();
1155 const double *denoRPtr=denoR->getArray()->getConstPointer();
1156 if(trgField->getMesh()->getMeshDimension()==-1)
1158 double *denoRPtr2=denoR->getArray()->getPointer();
1159 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1161 if(srcField->getMesh()->getMeshDimension()==-1)
1163 double *denoPtr2=deno->getArray()->getPointer();
1164 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1167 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1168 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1170 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1171 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1178 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1182 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1185 double *tmp=new double[inputNbOfCompo];
1186 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1188 if((*iter1).empty())
1191 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1195 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1196 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1197 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1199 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1200 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1206 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1208 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1210 double *tmp=new double[inputNbOfCompo];
1211 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1212 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1214 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1216 isReached[(*iter2).first]=true;
1217 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1218 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1223 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1225 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1228 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1230 matOut.resize(nbColsMatIn);
1232 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1233 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1234 matOut[(*iter2).first][id]=(*iter2).second;
1237 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1238 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1240 std::map<int,double> values;
1242 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1245 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1247 sum+=(*iter2).second;
1248 values[(*iter2).first]+=(*iter2).second;
1250 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1251 deno[idx][(*iter2).first]=sum;
1254 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1256 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1257 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1261 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1262 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1264 std::map<int,double> values;
1266 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1269 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1271 sum+=(*iter2).second;
1272 values[(*iter2).first]+=(*iter2).second;
1274 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1275 denoReverse[(*iter2).first][idx]=sum;
1278 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1280 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1281 deno[idx][(*iter2).first]=values[(*iter2).first];
1285 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1286 const std::vector< std::map<int,double> >& m2D,
1287 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1288 const int *corrCellIdTrg)
1290 int nbOf2DCellsTrg=m2D.size();
1291 int nbOf1DCellsTrg=m1D.size();
1292 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1293 _matrix.resize(nbOf3DCellsTrg);
1295 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1297 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1300 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1302 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1304 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1311 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1314 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1316 std::cout << "Target Cell # " << id << " : ";
1317 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1318 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1319 std::cout << std::endl;
1323 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1329 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1331 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1333 return (int)_deno_reverse_multiply.size();
1337 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1338 * If not the behaviour is unpredictable.
1339 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1340 * 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.
1341 * 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
1344 * \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.
1345 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1346 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1348 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1351 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1353 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1355 std::map<int,double>& rowNew=matrixNew[i];
1356 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1358 if(fabs((*it2).second)>maxValAbs)
1359 rowNew[(*it2).first]=(*it2).second;
1370 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1371 * If not the behaviour is unpredictable.
1372 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1373 * 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.
1374 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1375 * 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
1378 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1379 * \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
1380 * that all coefficients are null.
1381 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1383 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1385 double maxVal=getMaxValueInCrudeMatrix();
1388 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1392 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1393 * If not the behaviour is unpredictable.
1394 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1395 * The returned value is positive.
1397 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1400 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1401 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1402 if(fabs((*it2).second)>ret)
1403 ret=fabs((*it2).second);