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<mcIdType,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<mcIdType,double> >& m)
106 restartUsing(src,target);
107 if(ToIdType(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 std::size_t trgNbOfCompo=targetField->getNumberOfComponents();
251 srcField->checkConsistencyLight();
252 if(ToIdType(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,(int)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==3 && trgMeshDim==0 && srcSpaceDim==3)
469 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
470 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 0D ! 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::Interpolation3D1D interpolation(*this);//Not a bug : 3D1D deal with 3D0D
474 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
476 else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
478 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
479 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! 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::Interpolation1D0D interpolation(*this);
483 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
485 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
487 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
488 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
489 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
490 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
491 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
492 std::vector<std::map<mcIdType,double> > matrixTmp;
493 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
494 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
495 ReverseMatrix(matrixTmp,nbCols,_matrix);
496 nbCols=ToIdType(matrixTmp.size());
498 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
500 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
502 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
503 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
504 INTERP_KERNEL::Interpolation2D interpolation(*this);
505 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
509 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
510 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
511 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
512 std::vector<std::map<mcIdType,double> > matrixTmp;
513 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
514 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
515 ReverseMatrix(matrixTmp,nbCols,_matrix);
516 nbCols=ToIdType(matrixTmp.size());
517 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
518 if(!duplicateFaces.empty())
520 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
521 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
523 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
524 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
530 else if(srcMeshDim==2 && trgMeshDim==0 && srcSpaceDim==2)
532 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
533 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 2D and 0D ! Select PointLocator as intersection type !");
534 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
535 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
536 INTERP_KERNEL::Interpolation2D interpolation(*this);
537 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
539 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
541 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
543 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
544 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
545 INTERP_KERNEL::Interpolation2D interpolation(*this);
546 std::vector<std::map<mcIdType,double> > matrixTmp;
547 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
548 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
549 ReverseMatrix(matrixTmp,nbCols,_matrix);
550 nbCols=ToIdType(matrixTmp.size());
554 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
555 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
556 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
557 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
558 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
559 if(!duplicateFaces.empty())
561 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
562 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
564 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
565 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
571 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
573 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
574 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!!");
577 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
578 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
579 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
580 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
581 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
582 if(!duplicateFaces.empty())
584 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
585 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
587 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
588 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
594 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
596 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
598 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
599 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
600 INTERP_KERNEL::Interpolation3D interpolation(*this);
601 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
605 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
606 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
607 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
608 std::vector<std::map<mcIdType,double> > matrixTmp;
609 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
610 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
611 ReverseMatrix(matrixTmp,nbCols,_matrix);
612 nbCols=ToIdType(matrixTmp.size());
613 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
614 if(!duplicateFaces.empty())
616 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
617 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
619 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
620 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
626 else if(trgMeshDim==-1)
628 if(srcMeshDim==2 && srcSpaceDim==2)
630 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
631 INTERP_KERNEL::Interpolation2D interpolation(*this);
632 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
634 else if(srcMeshDim==3 && srcSpaceDim==3)
636 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
637 INTERP_KERNEL::Interpolation3D interpolation(*this);
638 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
640 else if(srcMeshDim==2 && srcSpaceDim==3)
642 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
643 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
644 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
647 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
649 else if(srcMeshDim==-1)
651 if(trgMeshDim==2 && trgSpaceDim==2)
653 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
654 INTERP_KERNEL::Interpolation2D interpolation(*this);
655 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
657 else if(trgMeshDim==3 && trgSpaceDim==3)
659 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
660 INTERP_KERNEL::Interpolation3D interpolation(*this);
661 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
663 else if(trgMeshDim==2 && trgSpaceDim==3)
665 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
666 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
667 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
670 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
673 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
674 _deno_multiply.clear();
675 _deno_multiply.resize(_matrix.size());
676 _deno_reverse_multiply.clear();
677 _deno_reverse_multiply.resize(nbCols);
682 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
684 std::string srcMeth,trgMeth;
685 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
686 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
687 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
689 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
690 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
691 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
692 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
693 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
694 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
695 std::vector<std::map<mcIdType,double> > matrix2D;
696 mcIdType nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
697 MEDCouplingUMesh *s1D,*t1D;
699 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
700 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
701 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
702 std::vector<std::map<mcIdType,double> > matrix1D;
703 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
704 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
705 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
706 mcIdType nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
709 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
710 target_mesh->getMesh3DIds()->getConstPointer());
712 _deno_multiply.clear();
713 _deno_multiply.resize(_matrix.size());
714 _deno_reverse_multiply.clear();
715 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
720 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
722 std::string srcMeth,trgMeth;
723 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
724 if(methodCpp!="P0P0")
725 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only P0P0 interpolation supported for the moment !");
726 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
727 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
728 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
729 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
730 const int srcMeshDim=src_mesh->getMeshDimension();
731 const int srcSpceDim=src_mesh->getSpaceDimension();
732 const int trgMeshDim=target_mesh->getMeshDimension();
733 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
734 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!");
735 std::vector<std::map<mcIdType,double> > res;
740 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
741 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
742 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
743 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
748 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
749 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
750 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
751 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
756 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
757 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
758 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
759 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
763 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
765 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
766 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
768 _deno_multiply.clear();
769 _deno_multiply.resize(_matrix.size());
770 _deno_reverse_multiply.clear();
771 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
776 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
778 std::string srcMeth,trgMeth;
779 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
780 if(methodCpp!="P0P0")
781 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
782 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
783 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
784 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
785 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
786 const int srcMeshDim=src_mesh->getMeshDimension();
787 const int trgMeshDim=target_mesh->getMeshDimension();
788 const int trgSpceDim=target_mesh->getSpaceDimension();
789 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
790 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!");
795 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
796 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
797 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
798 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
803 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
804 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
805 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
806 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
811 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
812 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
813 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
814 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
818 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
820 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
822 _deno_multiply.clear();
823 _deno_multiply.resize(_matrix.size());
824 _deno_reverse_multiply.clear();
825 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
830 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
832 std::string srcMeth,trgMeth;
833 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
834 if(methodCpp!="P0P0")
835 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
836 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
837 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
838 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
839 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
840 const int srcMeshDim=src_mesh->getMeshDimension();
841 const int trgMeshDim=target_mesh->getMeshDimension();
842 if(trgMeshDim!=srcMeshDim)
843 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
848 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
849 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
850 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
851 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
856 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
857 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
858 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
859 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
864 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
865 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
866 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
867 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
871 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
873 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
875 _deno_multiply.clear();
876 _deno_multiply.resize(_matrix.size());
877 _deno_reverse_multiply.clear();
878 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
883 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
885 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
886 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 !");
887 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
888 const double *trgLocPtr=trgLoc->begin();
889 mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents());
890 MCAuto<DataArrayIdType> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
891 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
893 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
894 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
895 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
896 throw INTERP_KERNEL::Exception(oss.str().c_str());
898 const mcIdType *srcOffsetArrPtr=srcOffsetArr->begin();
899 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
900 const double *srcLocPtr=srcLoc->begin();
901 MCAuto<DataArrayIdType> eltsArr,eltsIndexArr;
902 mcIdType trgNbOfGaussPts=trgLoc->getNumberOfTuples();
903 _matrix.resize(trgNbOfGaussPts);
904 _src_ft->getMesh()->getCellsContainingPointsLinearPartOnlyOnNonDynType(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
905 const mcIdType *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
906 MCAuto<DataArrayIdType> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
907 MCAuto<DataArrayIdType> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
908 for(const mcIdType *trgId=ids0->begin();trgId!=ids0->end();trgId++)
910 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
911 mcIdType srcCellId=elts[eltsIndex[*trgId]];
912 double dist=std::numeric_limits<double>::max();
913 mcIdType srcEntry=-1;
914 for(mcIdType srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
916 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
918 for(mcIdType i=0;i<trgSpaceDim;i++)
919 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
921 { dist=tmp; srcEntry=srcId; }
923 _matrix[*trgId][srcEntry]=1.;
925 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
927 MCAuto<DataArrayIdType> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
928 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
929 MCAuto<DataArrayIdType> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
930 const mcIdType *srcIdPerTrgPtr=srcIdPerTrg->begin();
931 for(const mcIdType *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
932 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
934 _deno_multiply.clear();
935 _deno_multiply.resize(_matrix.size());
936 _deno_reverse_multiply.clear();
937 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
943 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
944 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
946 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
948 if(method=="GAUSSGAUSS")
950 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
951 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
952 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
953 throw INTERP_KERNEL::Exception(oss.str().c_str());
957 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
958 * 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
959 * only INTERP_KERNEL method should be applied.
961 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
963 std::string srcm,trgm,method;
964 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
965 switch(_interp_matrix_pol)
967 case IK_ONLY_PREFERED:
971 std::string tmp1,tmp2;
972 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
975 catch(INTERP_KERNEL::Exception& /*e*/)
980 case NOT_IK_ONLY_PREFERED:
984 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
987 catch(INTERP_KERNEL::Exception& /*e*/)
994 case NOT_IK_ONLY_FORCED:
997 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
1001 void MEDCouplingRemapper::updateTime() const
1005 void MEDCouplingRemapper::checkPrepare() const
1007 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1009 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
1010 if(!s->getMesh() || !t->getMesh())
1011 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
1015 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
1016 * This method returns 3 information (2 in output parameters and 1 in return).
1018 * \param [out] srcMeth the string code of the discretization of source field template
1019 * \param [out] trgMeth the string code of the discretization of target field template
1020 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
1022 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
1024 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1026 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1027 if(!s->getMesh() || !t->getMesh())
1028 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1029 srcMeth=_src_ft->getDiscretization()->getRepr();
1030 trgMeth=_target_ft->getDiscretization()->getRepr();
1031 return BuildMethodFrom(srcMeth,trgMeth);
1034 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1036 std::string method(meth1); method+=meth2;
1040 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1042 if(!srcMesh || !targetMesh)
1043 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1044 std::string srcMethod,targetMethod;
1045 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1046 src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1047 src->setMesh(srcMesh);
1048 target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1049 target->setMesh(targetMesh);
1052 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1056 if(matrixSuppression)
1059 _deno_multiply.clear();
1060 _deno_reverse_multiply.clear();
1064 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1067 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1068 if(!src->getMesh() || !target->getMesh())
1069 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1071 _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1072 _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1075 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1077 if(!srcField || !targetField)
1078 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1079 srcField->checkConsistencyLight();
1081 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1082 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1083 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1084 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1085 if(srcField->getNature()!=targetField->getNature())
1086 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1087 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1089 std::ostringstream oss;
1090 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();
1091 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1092 throw INTERP_KERNEL::Exception(oss.str().c_str());
1094 DataArrayDouble *array(targetField->getArray());
1095 std::size_t srcNbOfCompo(srcField->getNumberOfComponents());
1098 targetField->checkConsistencyLight();
1099 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1100 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1105 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1106 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1107 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1108 targetField->setArray(tmp);
1110 computeDeno(srcField->getNature(),srcField,targetField);
1111 double *resPointer(targetField->getArray()->getPointer());
1112 const double *inputPointer(srcField->getArray()->getConstPointer());
1113 computeProduct(inputPointer,(int)srcNbOfCompo,isDftVal,dftValue,resPointer);
1116 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1119 return computeDenoFromScratch(nat,srcField,trgField);
1120 else if(nat!=_nature_of_deno)
1121 return computeDenoFromScratch(nat,srcField,trgField);
1122 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1123 return computeDenoFromScratch(nat,srcField,trgField);
1126 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1128 _nature_of_deno=nat;
1129 std::size_t _time_deno_update=getTimeOfThis();
1130 switch(_nature_of_deno)
1132 case IntensiveMaximum:
1134 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1137 case ExtensiveMaximum:
1139 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1140 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1141 const double *denoPtr=deno->getArray()->getConstPointer();
1142 const double *denoRPtr=denoR->getArray()->getConstPointer();
1143 if(trgField->getMesh()->getMeshDimension()==-1)
1145 double *denoRPtr2=denoR->getArray()->getPointer();
1146 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1148 if(srcField->getMesh()->getMeshDimension()==-1)
1150 double *denoPtr2=deno->getArray()->getPointer();
1151 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1154 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1155 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1157 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1158 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1164 case ExtensiveConservation:
1166 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1169 case IntensiveConservation:
1171 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1172 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1173 const double *denoPtr=deno->getArray()->getConstPointer();
1174 const double *denoRPtr=denoR->getArray()->getConstPointer();
1175 if(trgField->getMesh()->getMeshDimension()==-1)
1177 double *denoRPtr2=denoR->getArray()->getPointer();
1178 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1180 if(srcField->getMesh()->getMeshDimension()==-1)
1182 double *denoPtr2=deno->getArray()->getPointer();
1183 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1186 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1187 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1189 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1190 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1197 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1201 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1204 double *tmp=new double[inputNbOfCompo];
1205 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1207 if((*iter1).empty())
1210 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1214 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1215 std::map<mcIdType,double>::const_iterator iter3=_deno_multiply[idx].begin();
1216 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1218 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1219 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1225 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1227 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1229 double *tmp=new double[inputNbOfCompo];
1230 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1231 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1233 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1235 isReached[(*iter2).first]=true;
1236 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1237 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1242 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1244 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1247 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<mcIdType,double> >& matIn, mcIdType nbColsMatIn, std::vector<std::map<mcIdType,double> >& matOut)
1249 matOut.resize(nbColsMatIn);
1251 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1252 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1253 matOut[(*iter2).first][id]=(*iter2).second;
1256 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1257 std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1259 std::map<mcIdType,double> values;
1261 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1264 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1266 sum+=(*iter2).second;
1267 values[(*iter2).first]+=(*iter2).second;
1269 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1270 deno[idx][(*iter2).first]=sum;
1273 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1275 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1276 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1280 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1281 std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1283 std::map<mcIdType,double> values;
1285 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1288 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1290 sum+=(*iter2).second;
1291 values[(*iter2).first]+=(*iter2).second;
1293 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1294 denoReverse[(*iter2).first][idx]=sum;
1297 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1299 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1300 deno[idx][(*iter2).first]=values[(*iter2).first];
1304 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<mcIdType,double> >& m1D,
1305 const std::vector< std::map<mcIdType,double> >& m2D,
1306 const mcIdType *corrCellIdSrc, mcIdType nbOf2DCellsSrc, mcIdType nbOf1DCellsSrc,
1307 const mcIdType *corrCellIdTrg)
1309 mcIdType nbOf2DCellsTrg=ToIdType(m2D.size());
1310 mcIdType nbOf1DCellsTrg=ToIdType(m1D.size());
1311 mcIdType nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1312 _matrix.resize(nbOf3DCellsTrg);
1314 for(std::vector< std::map<mcIdType,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1316 for(std::map<mcIdType,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1319 for(std::vector< std::map<mcIdType,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1321 for(std::map<mcIdType,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1323 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1330 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<mcIdType,double> >& m)
1333 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1335 std::cout << "Target Cell # " << id << " : ";
1336 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1337 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1338 std::cout << std::endl;
1342 const std::vector<std::map<mcIdType,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1348 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1350 mcIdType MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1352 return ToIdType(_deno_reverse_multiply.size());
1356 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1357 * If not the behaviour is unpredictable.
1358 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1359 * 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.
1360 * 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
1363 * \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.
1364 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1365 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1367 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1370 std::vector<std::map<mcIdType,double> > matrixNew(_matrix.size());
1372 for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1374 std::map<mcIdType,double>& rowNew=matrixNew[i];
1375 for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1377 if(fabs((*it2).second)>maxValAbs)
1378 rowNew[(*it2).first]=(*it2).second;
1389 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1390 * If not the behaviour is unpredictable.
1391 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1392 * 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.
1393 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1394 * 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
1397 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1398 * \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
1399 * that all coefficients are null.
1400 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1402 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1404 double maxVal=getMaxValueInCrudeMatrix();
1407 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1411 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1412 * If not the behaviour is unpredictable.
1413 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1414 * The returned value is positive.
1416 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1419 for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1420 for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1421 if(fabs((*it2).second)>ret)
1422 ret=fabs((*it2).second);