1 // Copyright (C) 2007-2019 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingRemapper.hxx"
22 #include "MEDCouplingMemArray.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24 #include "MEDCouplingFieldTemplate.hxx"
25 #include "MEDCouplingFieldDiscretization.hxx"
26 #include "MEDCouplingMappedExtrudedMesh.hxx"
27 #include "MEDCouplingCMesh.hxx"
28 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
29 #include "MEDCouplingNormalizedCartesianMesh.txx"
31 #include "Interpolation1D.txx"
32 #include "Interpolation2DCurve.hxx"
33 #include "Interpolation2D.txx"
34 #include "Interpolation3D.txx"
35 #include "Interpolation3DSurf.hxx"
36 #include "Interpolation2D1D.txx"
37 #include "Interpolation2D3D.txx"
38 #include "Interpolation3D1D.txx"
39 #include "Interpolation1D0D.txx"
40 #include "InterpolationCU.txx"
41 #include "InterpolationCC.txx"
43 using namespace MEDCoupling;
45 MEDCouplingRemapper::MEDCouplingRemapper():_src_ft(0),_target_ft(0),_interp_matrix_pol(IK_ONLY_PREFERED),_nature_of_deno(NoNature),_time_deno_update(0)
49 MEDCouplingRemapper::~MEDCouplingRemapper()
55 * This method is the second step of the remapping process. The remapping process works in three phases :
57 * - Set remapping options appropriately
58 * - The computation of remapping matrix
59 * - Apply the matrix vector multiply to obtain the result of the remapping
61 * This method performs the second step (computation of remapping matrix) which may be CPU-time consuming. This phase is also the most critical (where the most tricky algorithm) in the remapping process.
62 * Strictly speaking to perform the computation of the remapping matrix the field templates source-side and target-side is required (which is the case of MEDCouplingRemapper::prepareEx).
63 * So this method is less precise but a more user friendly way to compute a remapping matrix.
65 * \param [in] srcMesh the source mesh
66 * \param [in] targetMesh the target mesh
67 * \param [in] method A string obtained by aggregation of the spatial discretisation string representation of source field and target field. The string representation is those returned by MEDCouplingFieldDiscretization::getStringRepr.
68 * Example : "P0" is for cell discretization. "P1" is for node discretization. So "P0P1" for \a method parameter means from a source cell field (lying on \a srcMesh) to a target node field (lying on \a targetMesh).
70 * \sa MEDCouplingRemapper::prepareEx
72 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
74 MCAuto<MEDCouplingFieldTemplate> src,target;
75 BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
76 return prepareEx(src,target);
80 * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
81 * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
83 * \param [in] src is the field template source side.
84 * \param [in] target is the field template target side.
86 * \sa MEDCouplingRemapper::prepare
88 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
90 restartUsing(src,target);
91 if(isInterpKernelOnlyOrNotOnly())
92 return prepareInterpKernelOnly();
94 return prepareNotInterpKernelOnly();
97 void MEDCouplingRemapper::setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector<std::map<int,double> >& m)
99 MCAuto<MEDCouplingFieldTemplate> src,target;
100 BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
101 setCrudeMatrixEx(src,target,m);
104 void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector<std::map<int,double> >& m)
106 restartUsing(src,target);
107 if(m.size()!=target->getNumberOfTuplesExpected())
109 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : input matrix has " << m.size() << " rows whereas there are " << target->getNumberOfTuplesExpected() << " expected !";
110 throw INTERP_KERNEL::Exception(oss.str());
112 auto srcNbElem(src->getNumberOfTuplesExpected());
117 auto idToTest(it2.first);
118 if(idToTest<0 || idToTest>=srcNbElem)
120 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : presence of elt #" << idToTest << " ! not in [0," << srcNbElem << ") !";
121 throw INTERP_KERNEL::Exception(oss.str());
126 _deno_multiply.clear();
127 _deno_multiply.resize(_matrix.size());
128 _deno_reverse_multiply.clear();
129 _deno_reverse_multiply.resize(srcNbElem);
132 int MEDCouplingRemapper::prepareInterpKernelOnly()
134 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
142 // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
143 // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
145 // } MEDCouplingMeshType;
147 switch(meshInterpType)
149 case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
150 case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
151 case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
152 case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
153 case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
154 case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
155 case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
156 case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
157 case 85: // UNSTRUCTURED - UNSTRUCTURED
158 return prepareInterpKernelOnlyUU();
159 case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
160 case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
161 case 87: // UNSTRUCTURED - CARTESIAN
162 return prepareInterpKernelOnlyUC();
163 case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
164 case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
165 case 117: // CARTESIAN - UNSTRUCTURED
166 return prepareInterpKernelOnlyCU();
167 case 119: // CARTESIAN - CARTESIAN
168 return prepareInterpKernelOnlyCC();
169 case 136: // EXTRUDED - EXTRUDED
170 return prepareInterpKernelOnlyEE();
172 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
176 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
178 std::string srcm,trgm,method;
179 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
180 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
183 return prepareNotInterpKernelOnlyGaussGauss();
186 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
187 throw INTERP_KERNEL::Exception(oss.str().c_str());
193 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
194 * If meshes of \b srcField and \b targetField do not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
196 * \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on MEDCoupling::MEDCouplingRemapper::prepare.
197 * \param [in/out] targetField the destination field with the allocated array in which all tuples will be overwritten.
198 * \param [in] dftValue is the value that will be assigned in the targetField to each entity of target mesh (entity depending on the method selected on prepare invocation) that is not intercepted by any entity of source mesh.
199 * For example in "P0P0" case (cell-cell) if a cell in target mesh is not overlapped by any source cell the \a dftValue value will be attached on that cell in the returned \a targetField. In some cases a target
200 * cell not intercepted by any source cell is a bug so in this case it is advised to set a huge value (1e300 for example) to \a dftValue to quickly point to the problem. But for users doing parallelism a target cell can
201 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
205 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue)
207 if(!srcField || !targetField)
208 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transfer : input field must be both not NULL !");
209 transferUnderground(srcField,targetField,true,dftValue);
213 * This method is equivalent to MEDCoupling::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
214 * If an entity (cell for example) in targetField is not fetched by any entity (cell for example) of \b srcField, the value in targetField is
216 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
218 * \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on MEDCoupling::MEDCouplingRemapper::prepare.
219 * \param [in,out] targetField the destination field with the allocated array in which only tuples whose entities are fetched by interpolation will be overwritten only.
221 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField)
223 if(!srcField || !targetField)
224 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : input field must be both not NULL !");
225 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
228 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue)
230 if(!srcField || !targetField)
231 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::reverseTransfer : input fields must be both not NULL !");
233 targetField->checkConsistencyLight();
234 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
235 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
236 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
237 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
238 if(srcField->getNature()!=targetField->getNature())
239 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
240 if(targetField->getNumberOfTuplesExpected()!=_target_ft->getNumberOfTuplesExpected())
242 std::ostringstream oss;
243 oss << "MEDCouplingRemapper::reverseTransfer : in given source field the number of tuples required is " << _target_ft->getNumberOfTuplesExpected() << " (on prepare) and number of tuples in given target field is " << targetField->getNumberOfTuplesExpected();
244 oss << " ! It appears that the target support is not the same between the prepare and the transfer !";
245 throw INTERP_KERNEL::Exception(oss.str().c_str());
247 DataArrayDouble *array(srcField->getArray());
248 int trgNbOfCompo=targetField->getNumberOfComponents();
251 srcField->checkConsistencyLight();
252 if(trgNbOfCompo!=srcField->getNumberOfTuplesExpected())
253 throw INTERP_KERNEL::Exception("Number of components mismatch !");
257 MCAuto<DataArrayDouble > tmp(DataArrayDouble::New());
258 tmp->alloc(srcField->getNumberOfTuplesExpected(),trgNbOfCompo);
259 srcField->setArray(tmp);
261 computeDeno(srcField->getNature(),srcField,targetField);
262 double *resPointer(srcField->getArray()->getPointer());
263 const double *inputPointer=targetField->getArray()->getConstPointer();
264 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
268 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
269 * If mesh of \b srcField does not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
271 * \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on MEDCoupling::MEDCouplingRemapper::prepare.
272 * \param [in] dftValue is the value that will be assigned in the targetField to each entity of target mesh (entity depending on the method selected on prepare invocation) that is not intercepted by any entity of source mesh.
273 * For example in "P0P0" case (cell-cell) if a cell in target mesh is not overlapped by any source cell the \a dftValue value will be attached on that cell in the returned \a targetField. In some cases a target
274 * cell not intercepted by any source cell is a bug so in this case it is advised to set a huge value (1e300 for example) to \a dftValue to quickly point to the problem. But for users doing parallelism a target cell can
275 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
276 * \return destination field to be deallocated by the caller.
280 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue)
284 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input srcField is NULL !");
285 srcField->checkConsistencyLight();
286 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
287 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
288 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
289 ret->setNature(srcField->getNature());
290 transfer(srcField,ret,dftValue);
291 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
295 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue)
298 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input targetField is NULL !");
299 targetField->checkConsistencyLight();
301 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
302 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
303 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
304 ret->setNature(targetField->getNature());
305 reverseTransfer(ret,targetField,dftValue);
306 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
311 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
312 * is here only for automatic CORBA generators.
314 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
316 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
320 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
321 * is here only for automatic CORBA generators.
323 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
325 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
329 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
330 * is here only for automatic CORBA generators.
332 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
334 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
338 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or preferred.
339 * If interpolation matrix policy is :
341 * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is preferred. That is to say, if it is possible to treat the case
342 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
343 * If not, the \b not only INTERP_KERNEL method will be attempt.
345 * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is preferred. That is to say, if it is possible to treat the case
346 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
347 * If not, the INTERP_KERNEL only method will be attempt.
349 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
351 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
353 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
355 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
357 return _interp_matrix_pol;
361 * This method sets a new interpolation matrix policy. The default one is IK_PREFERED (0). The input is of type \c int to be dealt by standard Salome
362 * CORBA component generators. This method throws an INTERP_KERNEL::Exception if a the input integer is not in the available possibilities, that is to say not in
363 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
365 * If interpolation matrix policy is :
367 * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is preferred. That is to say, if it is possible to treat the case
368 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
369 * If not, the \b not only INTERP_KERNEL method will be attempt.
371 * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is preferred. That is to say, if it is possible to treat the case
372 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
373 * If not, the INTERP_KERNEL only method will be attempt.
375 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
377 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
379 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c MEDCoupling::InterpolationMatrixPolicy
380 * for automatic generation of CORBA component.
382 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
384 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
386 switch(newInterpMatPol)
389 _interp_matrix_pol=IK_ONLY_PREFERED;
392 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
395 _interp_matrix_pol=IK_ONLY_FORCED;
398 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
401 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::setInterpolationMatrixPolicy : invalid input integer value ! Should be in [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)] ! For information, the default is IK_PREFERED=0 !");
405 int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
407 const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
408 const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_target_ft->getMesh());
409 std::string srcMeth,trgMeth;
410 std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth));
411 const int srcMeshDim=src_mesh->getMeshDimension();
414 srcSpaceDim=src_mesh->getSpaceDimension();
415 const int trgMeshDim=target_mesh->getMeshDimension();
418 trgSpaceDim=target_mesh->getSpaceDimension();
419 if(trgSpaceDim!=srcSpaceDim)
420 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
421 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
423 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
425 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
426 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
427 INTERP_KERNEL::Interpolation1D interpolation(*this);
428 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
430 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
432 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
433 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
434 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
435 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
437 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
439 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
440 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
441 INTERP_KERNEL::Interpolation2D interpolation(*this);
442 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
444 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
446 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
447 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
448 INTERP_KERNEL::Interpolation3D interpolation(*this);
449 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
451 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
453 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
454 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
455 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
456 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
458 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
460 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
461 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
462 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
463 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
464 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
465 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
467 else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
469 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
470 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! Select PointLocator as intersection type !");
471 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
472 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
473 INTERP_KERNEL::Interpolation1D0D interpolation(*this);
474 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
476 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
478 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
479 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
480 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
481 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
482 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
483 std::vector<std::map<int,double> > matrixTmp;
484 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
485 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
486 ReverseMatrix(matrixTmp,nbCols,_matrix);
487 nbCols=matrixTmp.size();
489 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
491 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
493 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
494 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
495 INTERP_KERNEL::Interpolation2D interpolation(*this);
496 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
500 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
501 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
502 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
503 std::vector<std::map<int,double> > matrixTmp;
504 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
505 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
506 ReverseMatrix(matrixTmp,nbCols,_matrix);
507 nbCols=matrixTmp.size();
508 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
509 if(!duplicateFaces.empty())
511 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
512 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
514 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
515 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
521 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
523 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
525 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
526 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
527 INTERP_KERNEL::Interpolation2D interpolation(*this);
528 std::vector<std::map<int,double> > matrixTmp;
529 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
530 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
531 ReverseMatrix(matrixTmp,nbCols,_matrix);
532 nbCols=matrixTmp.size();
536 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
537 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
538 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
539 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
540 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
541 if(!duplicateFaces.empty())
543 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
544 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
546 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
547 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
553 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
555 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
556 throw INTERP_KERNEL::Exception("Using point locator to transfer a mesh of dim 2 to a mesh of dim 3 does not make sense: 3D centers of mass can not be localized in a mesh having mesh-dim=2, space-dim=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," "));
576 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
578 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
580 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
581 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
582 INTERP_KERNEL::Interpolation3D interpolation(*this);
583 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
587 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
588 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
589 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
590 std::vector<std::map<int,double> > matrixTmp;
591 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
592 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
593 ReverseMatrix(matrixTmp,nbCols,_matrix);
594 nbCols=matrixTmp.size();
595 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
596 if(!duplicateFaces.empty())
598 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
599 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
601 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
602 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
608 else if(trgMeshDim==-1)
610 if(srcMeshDim==2 && srcSpaceDim==2)
612 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
613 INTERP_KERNEL::Interpolation2D interpolation(*this);
614 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
616 else if(srcMeshDim==3 && srcSpaceDim==3)
618 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
619 INTERP_KERNEL::Interpolation3D interpolation(*this);
620 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
622 else if(srcMeshDim==2 && srcSpaceDim==3)
624 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
625 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
626 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
629 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
631 else if(srcMeshDim==-1)
633 if(trgMeshDim==2 && trgSpaceDim==2)
635 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
636 INTERP_KERNEL::Interpolation2D interpolation(*this);
637 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
639 else if(trgMeshDim==3 && trgSpaceDim==3)
641 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
642 INTERP_KERNEL::Interpolation3D interpolation(*this);
643 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
645 else if(trgMeshDim==2 && trgSpaceDim==3)
647 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
648 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
649 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
652 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
655 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
656 _deno_multiply.clear();
657 _deno_multiply.resize(_matrix.size());
658 _deno_reverse_multiply.clear();
659 _deno_reverse_multiply.resize(nbCols);
664 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
666 std::string srcMeth,trgMeth;
667 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
668 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
669 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
671 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
672 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
673 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
674 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
675 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
676 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
677 std::vector<std::map<int,double> > matrix2D;
678 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
679 MEDCouplingUMesh *s1D,*t1D;
681 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
682 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
683 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
684 std::vector<std::map<int,double> > matrix1D;
685 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
686 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
687 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
688 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
691 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
692 target_mesh->getMesh3DIds()->getConstPointer());
694 _deno_multiply.clear();
695 _deno_multiply.resize(_matrix.size());
696 _deno_reverse_multiply.clear();
697 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
702 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
704 std::string srcMeth,trgMeth;
705 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
706 if(methodCpp!="P0P0")
707 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only P0P0 interpolation supported for the moment !");
708 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
709 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
710 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
711 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
712 const int srcMeshDim=src_mesh->getMeshDimension();
713 const int srcSpceDim=src_mesh->getSpaceDimension();
714 const int trgMeshDim=target_mesh->getMeshDimension();
715 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
716 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!");
717 std::vector<std::map<int,double> > res;
722 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
723 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
724 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
725 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
730 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
731 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
732 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
733 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
738 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
739 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
740 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
741 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
745 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
747 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
748 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
750 _deno_multiply.clear();
751 _deno_multiply.resize(_matrix.size());
752 _deno_reverse_multiply.clear();
753 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
758 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
760 std::string srcMeth,trgMeth;
761 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
762 if(methodCpp!="P0P0")
763 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
764 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
765 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
766 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
767 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
768 const int srcMeshDim=src_mesh->getMeshDimension();
769 const int trgMeshDim=target_mesh->getMeshDimension();
770 const int trgSpceDim=target_mesh->getSpaceDimension();
771 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
772 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!");
777 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
778 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
779 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
780 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
785 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
786 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
787 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
788 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
793 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
794 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
795 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
796 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
800 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
802 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
804 _deno_multiply.clear();
805 _deno_multiply.resize(_matrix.size());
806 _deno_reverse_multiply.clear();
807 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
812 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
814 std::string srcMeth,trgMeth;
815 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
816 if(methodCpp!="P0P0")
817 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
818 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
819 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
820 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
821 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
822 const int srcMeshDim=src_mesh->getMeshDimension();
823 const int trgMeshDim=target_mesh->getMeshDimension();
824 if(trgMeshDim!=srcMeshDim)
825 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
830 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
831 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
832 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
833 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
838 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
839 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
840 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
841 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
846 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
847 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
848 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
849 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
853 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
855 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
857 _deno_multiply.clear();
858 _deno_multiply.resize(_matrix.size());
859 _deno_reverse_multiply.clear();
860 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
865 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
867 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
868 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 !");
869 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
870 const double *trgLocPtr=trgLoc->begin();
871 int trgSpaceDim=trgLoc->getNumberOfComponents();
872 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
873 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
875 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
876 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
877 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
878 throw INTERP_KERNEL::Exception(oss.str().c_str());
880 const int *srcOffsetArrPtr=srcOffsetArr->begin();
881 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
882 const double *srcLocPtr=srcLoc->begin();
883 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
884 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
885 _matrix.resize(trgNbOfGaussPts);
886 _src_ft->getMesh()->getCellsContainingPointsLinearPartOnlyOnNonDynType(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
887 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
888 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
889 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
890 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
892 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
893 int srcCellId=elts[eltsIndex[*trgId]];
894 double dist=std::numeric_limits<double>::max();
896 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
898 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
900 for(int i=0;i<trgSpaceDim;i++)
901 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
903 { dist=tmp; srcEntry=srcId; }
905 _matrix[*trgId][srcEntry]=1.;
907 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
909 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
910 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
911 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
912 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
913 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
914 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
916 _deno_multiply.clear();
917 _deno_multiply.resize(_matrix.size());
918 _deno_reverse_multiply.clear();
919 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
925 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
926 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
928 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
930 if(method=="GAUSSGAUSS")
932 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
933 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
934 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
935 throw INTERP_KERNEL::Exception(oss.str().c_str());
939 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
940 * 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
941 * only INTERP_KERNEL method should be applied.
943 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
945 std::string srcm,trgm,method;
946 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
947 switch(_interp_matrix_pol)
949 case IK_ONLY_PREFERED:
953 std::string tmp1,tmp2;
954 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
957 catch(INTERP_KERNEL::Exception& /*e*/)
962 case NOT_IK_ONLY_PREFERED:
966 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
969 catch(INTERP_KERNEL::Exception& /*e*/)
976 case NOT_IK_ONLY_FORCED:
979 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
983 void MEDCouplingRemapper::updateTime() const
987 void MEDCouplingRemapper::checkPrepare() const
989 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
991 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
992 if(!s->getMesh() || !t->getMesh())
993 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
997 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
998 * This method returns 3 information (2 in output parameters and 1 in return).
1000 * \param [out] srcMeth the string code of the discretization of source field template
1001 * \param [out] trgMeth the string code of the discretization of target field template
1002 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
1004 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
1006 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1008 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1009 if(!s->getMesh() || !t->getMesh())
1010 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1011 srcMeth=_src_ft->getDiscretization()->getRepr();
1012 trgMeth=_target_ft->getDiscretization()->getRepr();
1013 return BuildMethodFrom(srcMeth,trgMeth);
1016 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1018 std::string method(meth1); method+=meth2;
1022 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1024 if(!srcMesh || !targetMesh)
1025 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1026 std::string srcMethod,targetMethod;
1027 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1028 src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1029 src->setMesh(srcMesh);
1030 target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1031 target->setMesh(targetMesh);
1034 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1038 if(matrixSuppression)
1041 _deno_multiply.clear();
1042 _deno_reverse_multiply.clear();
1046 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1049 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1050 if(!src->getMesh() || !target->getMesh())
1051 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1053 _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1054 _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1057 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1059 if(!srcField || !targetField)
1060 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1061 srcField->checkConsistencyLight();
1063 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1064 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1065 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1066 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1067 if(srcField->getNature()!=targetField->getNature())
1068 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1069 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1071 std::ostringstream oss;
1072 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();
1073 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1074 throw INTERP_KERNEL::Exception(oss.str().c_str());
1076 DataArrayDouble *array(targetField->getArray());
1077 int srcNbOfCompo(srcField->getNumberOfComponents());
1080 targetField->checkConsistencyLight();
1081 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1082 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1087 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1088 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1089 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1090 targetField->setArray(tmp);
1092 computeDeno(srcField->getNature(),srcField,targetField);
1093 double *resPointer(targetField->getArray()->getPointer());
1094 const double *inputPointer(srcField->getArray()->getConstPointer());
1095 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1098 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1101 return computeDenoFromScratch(nat,srcField,trgField);
1102 else if(nat!=_nature_of_deno)
1103 return computeDenoFromScratch(nat,srcField,trgField);
1104 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1105 return computeDenoFromScratch(nat,srcField,trgField);
1108 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1110 _nature_of_deno=nat;
1111 _time_deno_update=getTimeOfThis();
1112 switch(_nature_of_deno)
1114 case IntensiveMaximum:
1116 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1119 case ExtensiveMaximum:
1121 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1122 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1123 const double *denoPtr=deno->getArray()->getConstPointer();
1124 const double *denoRPtr=denoR->getArray()->getConstPointer();
1125 if(trgField->getMesh()->getMeshDimension()==-1)
1127 double *denoRPtr2=denoR->getArray()->getPointer();
1128 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1130 if(srcField->getMesh()->getMeshDimension()==-1)
1132 double *denoPtr2=deno->getArray()->getPointer();
1133 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1136 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1137 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1139 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1140 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1146 case ExtensiveConservation:
1148 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1151 case IntensiveConservation:
1153 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1154 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1155 const double *denoPtr=deno->getArray()->getConstPointer();
1156 const double *denoRPtr=denoR->getArray()->getConstPointer();
1157 if(trgField->getMesh()->getMeshDimension()==-1)
1159 double *denoRPtr2=denoR->getArray()->getPointer();
1160 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1162 if(srcField->getMesh()->getMeshDimension()==-1)
1164 double *denoPtr2=deno->getArray()->getPointer();
1165 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1168 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1169 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1171 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1172 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1179 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1183 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1186 double *tmp=new double[inputNbOfCompo];
1187 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1189 if((*iter1).empty())
1192 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1196 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1197 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1198 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1200 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1201 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1207 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1209 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1211 double *tmp=new double[inputNbOfCompo];
1212 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1213 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1215 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1217 isReached[(*iter2).first]=true;
1218 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1219 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1224 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1226 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1229 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1231 matOut.resize(nbColsMatIn);
1233 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1234 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1235 matOut[(*iter2).first][id]=(*iter2).second;
1238 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1239 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1241 std::map<int,double> values;
1243 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1246 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1248 sum+=(*iter2).second;
1249 values[(*iter2).first]+=(*iter2).second;
1251 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1252 deno[idx][(*iter2).first]=sum;
1255 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1257 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1258 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1262 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1263 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1265 std::map<int,double> values;
1267 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1270 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1272 sum+=(*iter2).second;
1273 values[(*iter2).first]+=(*iter2).second;
1275 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1276 denoReverse[(*iter2).first][idx]=sum;
1279 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1281 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1282 deno[idx][(*iter2).first]=values[(*iter2).first];
1286 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1287 const std::vector< std::map<int,double> >& m2D,
1288 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1289 const int *corrCellIdTrg)
1291 int nbOf2DCellsTrg=m2D.size();
1292 int nbOf1DCellsTrg=m1D.size();
1293 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1294 _matrix.resize(nbOf3DCellsTrg);
1296 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1298 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1301 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1303 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1305 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1312 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1315 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1317 std::cout << "Target Cell # " << id << " : ";
1318 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1319 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1320 std::cout << std::endl;
1324 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1330 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1332 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1334 return (int)_deno_reverse_multiply.size();
1338 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1339 * If not the behaviour is unpredictable.
1340 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1341 * 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.
1342 * 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
1345 * \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.
1346 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1347 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1349 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1352 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1354 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1356 std::map<int,double>& rowNew=matrixNew[i];
1357 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1359 if(fabs((*it2).second)>maxValAbs)
1360 rowNew[(*it2).first]=(*it2).second;
1371 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1372 * If not the behaviour is unpredictable.
1373 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1374 * 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.
1375 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1376 * 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
1379 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1380 * \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
1381 * that all coefficients are null.
1382 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1384 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1386 double maxVal=getMaxValueInCrudeMatrix();
1389 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1393 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1394 * If not the behaviour is unpredictable.
1395 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1396 * The returned value is positive.
1398 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1401 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1402 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1403 if(fabs((*it2).second)>ret)
1404 ret=fabs((*it2).second);