1 // Copyright (C) 2007-2022 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"
30 #include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
32 #include "Interpolation1D.txx"
33 #include "Interpolation2DCurve.hxx"
34 #include "Interpolation2D.txx"
35 #include "Interpolation3D.txx"
36 #include "Interpolation3DSurf.hxx"
37 #include "Interpolation2D1D.txx"
38 #include "Interpolation2D3D.txx"
39 #include "Interpolation3D1D.txx"
40 #include "Interpolation1D0D.txx"
41 #include "InterpolationCU.txx"
42 #include "InterpolationCC.txx"
44 using namespace MEDCoupling;
46 MEDCouplingRemapper::MEDCouplingRemapper():_src_ft(0),_target_ft(0),_interp_matrix_pol(IK_ONLY_PREFERED),_nature_of_deno(NoNature),_time_deno_update(0)
50 MEDCouplingRemapper::~MEDCouplingRemapper()
56 * This method is the second step of the remapping process. The remapping process works in three phases :
58 * - Set remapping options appropriately
59 * - The computation of remapping matrix
60 * - Apply the matrix vector multiply to obtain the result of the remapping
62 * 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.
63 * 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).
64 * So this method is less precise but a more user friendly way to compute a remapping matrix.
66 * \param [in] srcMesh the source mesh
67 * \param [in] targetMesh the target mesh
68 * \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.
69 * 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).
71 * \sa MEDCouplingRemapper::prepareEx
73 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
75 MCAuto<MEDCouplingFieldTemplate> src,target;
76 BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
77 return prepareEx(src,target);
81 * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
82 * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
84 * \param [in] src is the field template source side.
85 * \param [in] target is the field template target side.
87 * \sa MEDCouplingRemapper::prepare
89 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
91 restartUsing(src,target);
92 if(isInterpKernelOnlyOrNotOnly())
93 return prepareInterpKernelOnly();
95 return prepareNotInterpKernelOnly();
98 void MEDCouplingRemapper::setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector<std::map<mcIdType,double> >& m)
100 MCAuto<MEDCouplingFieldTemplate> src,target;
101 BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
102 setCrudeMatrixEx(src,target,m);
105 void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector<std::map<mcIdType,double> >& m)
107 restartUsing(src,target);
108 if(ToIdType(m.size())!=target->getNumberOfTuplesExpected())
110 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : input matrix has " << m.size() << " rows whereas there are " << target->getNumberOfTuplesExpected() << " expected !";
111 throw INTERP_KERNEL::Exception(oss.str());
113 auto srcNbElem(src->getNumberOfTuplesExpected());
118 auto idToTest(it2.first);
119 if(idToTest<0 || idToTest>=srcNbElem)
121 std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : presence of elt #" << idToTest << " ! not in [0," << srcNbElem << ") !";
122 throw INTERP_KERNEL::Exception(oss.str());
127 synchronizeSizeOfSideMatricesAfterMatrixComputation(srcNbElem);
130 int MEDCouplingRemapper::prepareInterpKernelOnly()
132 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
140 // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
141 // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
143 // } MEDCouplingMeshType;
145 switch(meshInterpType)
147 case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
148 case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
149 case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
150 case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
151 case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
152 case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
153 case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
154 case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
155 case 85: // UNSTRUCTURED - UNSTRUCTURED
156 return prepareInterpKernelOnlyUU();
157 case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
158 case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
159 case 87: // UNSTRUCTURED - CARTESIAN
160 return prepareInterpKernelOnlyUC();
161 case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
162 case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
163 case 117: // CARTESIAN - UNSTRUCTURED
164 return prepareInterpKernelOnlyCU();
165 case 119: // CARTESIAN - CARTESIAN
166 return prepareInterpKernelOnlyCC();
167 case 136: // EXTRUDED - EXTRUDED
168 return prepareInterpKernelOnlyEE();
170 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
174 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
176 std::string srcm,trgm,method;
177 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
178 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
181 return prepareNotInterpKernelOnlyGaussGauss();
183 return prepareNotInterpKernelOnlyFEFE();
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(trgNbOfCompo != srcField->getNumberOfComponents())
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==0 && srcSpaceDim==1 )
432 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
433 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 1D space ! Select PointLocator as intersection type !");
434 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
435 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
436 INTERP_KERNEL::Interpolation1D interpolation(*this);
437 nbCols=interpolation.interpolateMeshes0D(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
439 else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==2 )
441 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
442 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 2D space ! Select PointLocator as intersection type !");
443 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
444 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
445 INTERP_KERNEL::Interpolation1D interpolation(*this);
446 nbCols=interpolation.interpolateMeshes0D(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
448 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
450 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
451 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
452 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
453 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
455 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
457 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
458 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
459 INTERP_KERNEL::Interpolation2D interpolation(*this);
460 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
462 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
464 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
465 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
466 INTERP_KERNEL::Interpolation3D interpolation(*this);
467 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
469 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
471 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
472 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
473 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
474 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
476 else if(srcMeshDim==3 && trgMeshDim==1 && 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 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
485 else if(srcMeshDim==3 && trgMeshDim==0 && srcSpaceDim==3)
487 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
488 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 0D ! 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);//Not a bug : 3D1D deal with 3D0D
492 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
494 else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
496 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
497 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! Select PointLocator as intersection type !");
498 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
499 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
500 INTERP_KERNEL::Interpolation1D0D interpolation(*this);
501 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
503 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
505 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
506 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
507 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
508 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
509 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
510 std::vector<std::map<mcIdType,double> > matrixTmp;
511 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
512 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
513 ReverseMatrix(matrixTmp,nbCols,_matrix);
514 nbCols=ToIdType(matrixTmp.size());
516 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
518 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
520 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
521 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
522 INTERP_KERNEL::Interpolation2D interpolation(*this);
523 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
527 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
528 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
529 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
530 std::vector<std::map<mcIdType,double> > matrixTmp;
531 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
532 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
533 ReverseMatrix(matrixTmp,nbCols,_matrix);
534 nbCols=ToIdType(matrixTmp.size());
535 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
536 if(!duplicateFaces.empty())
538 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
539 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
541 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
542 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
548 else if(srcMeshDim==2 && trgMeshDim==0 && srcSpaceDim==2)
550 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
551 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 2D and 0D ! Select PointLocator as intersection type !");
552 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
553 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
554 INTERP_KERNEL::Interpolation2D interpolation(*this);
555 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
557 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
559 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
561 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
562 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
563 INTERP_KERNEL::Interpolation2D interpolation(*this);
564 std::vector<std::map<mcIdType,double> > matrixTmp;
565 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
566 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
567 ReverseMatrix(matrixTmp,nbCols,_matrix);
568 nbCols=ToIdType(matrixTmp.size());
572 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
573 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
574 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
575 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
576 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
577 if(!duplicateFaces.empty())
579 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
580 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
582 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
583 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
589 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
591 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
592 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!!");
595 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
596 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
597 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
598 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
599 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
600 if(!duplicateFaces.empty())
602 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
603 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
605 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
606 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
612 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
614 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
616 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
617 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
618 INTERP_KERNEL::Interpolation3D interpolation(*this);
619 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
623 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
624 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
625 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
626 std::vector<std::map<mcIdType,double> > matrixTmp;
627 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
628 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
629 ReverseMatrix(matrixTmp,nbCols,_matrix);
630 nbCols=ToIdType(matrixTmp.size());
631 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
632 if(!duplicateFaces.empty())
634 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
635 for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
637 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
638 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
644 else if(trgMeshDim==-1)
646 if(srcMeshDim==2 && srcSpaceDim==2)
648 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
649 INTERP_KERNEL::Interpolation2D interpolation(*this);
650 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
652 else if(srcMeshDim==3 && srcSpaceDim==3)
654 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
655 INTERP_KERNEL::Interpolation3D interpolation(*this);
656 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
658 else if(srcMeshDim==2 && srcSpaceDim==3)
660 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
661 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
662 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
665 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
667 else if(srcMeshDim==-1)
669 if(trgMeshDim==2 && trgSpaceDim==2)
671 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
672 INTERP_KERNEL::Interpolation2D interpolation(*this);
673 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
675 else if(trgMeshDim==3 && trgSpaceDim==3)
677 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
678 INTERP_KERNEL::Interpolation3D interpolation(*this);
679 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
681 else if(trgMeshDim==2 && trgSpaceDim==3)
683 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
684 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
685 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
688 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
691 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
692 synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols);
696 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
698 std::string srcMeth,trgMeth;
699 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
700 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
701 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
703 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
704 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
705 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
706 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
707 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
708 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
709 std::vector<std::map<mcIdType,double> > matrix2D;
710 mcIdType nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
711 MEDCouplingUMesh *s1D,*t1D;
713 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
714 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
715 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
716 std::vector<std::map<mcIdType,double> > matrix1D;
717 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
718 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
719 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
720 mcIdType nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
723 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
724 target_mesh->getMesh3DIds()->getConstPointer());
726 synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols2D*nbCols1D);
730 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
732 std::string srcMeth,trgMeth;
733 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
734 if(methodCpp!="P0P0")
735 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only P0P0 interpolation supported for the moment !");
736 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
737 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
738 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
739 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
740 const int srcMeshDim=src_mesh->getMeshDimension();
741 const int srcSpceDim=src_mesh->getSpaceDimension();
742 const int trgMeshDim=target_mesh->getMeshDimension();
743 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
744 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!");
745 std::vector<std::map<mcIdType,double> > res;
750 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
751 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
752 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
753 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
758 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
759 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
760 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
761 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
766 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
767 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
768 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
769 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
773 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
775 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
776 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
778 synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
782 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
784 std::string srcMeth,trgMeth;
785 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
786 if(methodCpp!="P0P0")
787 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
788 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
789 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
790 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
791 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
792 const int srcMeshDim=src_mesh->getMeshDimension();
793 const int trgMeshDim=target_mesh->getMeshDimension();
794 const int trgSpceDim=target_mesh->getSpaceDimension();
795 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
796 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!");
801 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
802 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
803 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
804 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
809 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
810 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
811 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
812 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
817 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
818 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
819 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
820 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
824 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
826 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
828 synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
832 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
834 std::string srcMeth,trgMeth;
835 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
836 if(methodCpp!="P0P0")
837 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
838 if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
839 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
840 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
841 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
842 const int srcMeshDim=src_mesh->getMeshDimension();
843 const int trgMeshDim=target_mesh->getMeshDimension();
844 if(trgMeshDim!=srcMeshDim)
845 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
850 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
851 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
852 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
853 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
858 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
859 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
860 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
861 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
866 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
867 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
868 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
869 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
873 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
875 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
877 synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
881 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
883 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
884 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 !");
885 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
886 const double *trgLocPtr=trgLoc->begin();
887 mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents());
888 MCAuto<DataArrayIdType> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
889 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
891 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
892 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
893 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
894 throw INTERP_KERNEL::Exception(oss.str().c_str());
896 const mcIdType *srcOffsetArrPtr=srcOffsetArr->begin();
897 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
898 const double *srcLocPtr=srcLoc->begin();
899 MCAuto<DataArrayIdType> eltsArr,eltsIndexArr;
900 mcIdType trgNbOfGaussPts=trgLoc->getNumberOfTuples();
901 _matrix.resize(trgNbOfGaussPts);
902 _src_ft->getMesh()->getCellsContainingPointsLinearPartOnlyOnNonDynType(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
903 const mcIdType *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
904 MCAuto<DataArrayIdType> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
905 MCAuto<DataArrayIdType> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
906 for(const mcIdType *trgId=ids0->begin();trgId!=ids0->end();trgId++)
908 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
909 mcIdType srcCellId=elts[eltsIndex[*trgId]];
910 double dist=std::numeric_limits<double>::max();
911 mcIdType srcEntry=-1;
912 for(mcIdType srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
914 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
916 for(mcIdType i=0;i<trgSpaceDim;i++)
917 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
919 { dist=tmp; srcEntry=srcId; }
921 _matrix[*trgId][srcEntry]=1.;
923 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
925 MCAuto<DataArrayIdType> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
926 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
927 MCAuto<DataArrayIdType> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
928 const mcIdType *srcIdPerTrgPtr=srcIdPerTrg->begin();
929 for(const mcIdType *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
930 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
932 synchronizeSizeOfSideMatricesAfterMatrixComputation( srcLoc->getNumberOfTuples() );
936 int MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE()
938 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
939 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE : The intersection type is not supported ! Only PointLocator is supported for FE->FE interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !");
940 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
941 mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents());
943 THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for target !")
944 if(_src_ft->getMesh()->getSpaceDimension() != 3)
945 THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for source !")
946 const MEDCouplingUMesh *srcUMesh( dynamic_cast<const MEDCouplingUMesh *>(_src_ft->getMesh()) );
947 const MEDCouplingPointSet *trgMesh( dynamic_cast<const MEDCouplingPointSet *>(_target_ft->getMesh()) );
949 THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D UMesh supported as source !");
951 THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D PointSet mesh supported as target !");
954 _matrix.resize(trgMesh->getNumberOfNodes());
957 auto matrixFeeder = [this,&rowId](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
959 auto& row = this->_matrix[rowId++];
960 MCAuto<DataArrayDouble> resVector( gl.getShapeFunctionValues() );
961 for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt)
963 row[ conn[iPt] ] = resVector->getIJ(0,iPt);
967 MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(srcUMesh,trgMesh->getCoords()->begin(),trgMesh->getNumberOfNodes(),matrixFeeder);
968 synchronizeSizeOfSideMatricesAfterMatrixComputation( srcUMesh->getNumberOfNodes() );
973 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
974 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
976 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
978 if(method=="GAUSSGAUSS")
982 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
983 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
984 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS FEFE !";
985 throw INTERP_KERNEL::Exception(oss.str().c_str());
989 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
990 * 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
991 * only INTERP_KERNEL method should be applied.
993 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
995 std::string srcm,trgm,method;
996 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
997 switch(_interp_matrix_pol)
999 case IK_ONLY_PREFERED:
1003 std::string tmp1,tmp2;
1004 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
1007 catch(INTERP_KERNEL::Exception& /*e*/)
1012 case NOT_IK_ONLY_PREFERED:
1016 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
1019 catch(INTERP_KERNEL::Exception& /*e*/)
1024 case IK_ONLY_FORCED:
1026 case NOT_IK_ONLY_FORCED:
1029 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
1033 void MEDCouplingRemapper::updateTime() const
1037 void MEDCouplingRemapper::checkPrepare() const
1039 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1041 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
1042 if(!s->getMesh() || !t->getMesh())
1043 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
1046 void MEDCouplingRemapper::synchronizeSizeOfSideMatricesAfterMatrixComputation(mcIdType nbOfColsInMatrix)
1048 _deno_multiply.clear();
1049 _deno_multiply.resize(_matrix.size());
1050 _deno_reverse_multiply.clear();
1051 _deno_reverse_multiply.resize(nbOfColsInMatrix);
1056 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
1057 * This method returns 3 information (2 in output parameters and 1 in return).
1059 * \param [out] srcMeth the string code of the discretization of source field template
1060 * \param [out] trgMeth the string code of the discretization of target field template
1061 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
1063 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
1065 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1067 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1068 if(!s->getMesh() || !t->getMesh())
1069 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1070 srcMeth=_src_ft->getDiscretization()->getRepr();
1071 trgMeth=_target_ft->getDiscretization()->getRepr();
1072 return BuildMethodFrom(srcMeth,trgMeth);
1075 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1077 std::string method(meth1); method+=meth2;
1081 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1083 if(!srcMesh || !targetMesh)
1084 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1085 std::string srcMethod,targetMethod;
1086 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1087 src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1088 src->setMesh(srcMesh);
1089 target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1090 target->setMesh(targetMesh);
1093 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1097 if(matrixSuppression)
1100 _deno_multiply.clear();
1101 _deno_reverse_multiply.clear();
1105 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1108 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1109 if(!src->getMesh() || !target->getMesh())
1110 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1112 _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1113 _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1116 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1118 if(!srcField || !targetField)
1119 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1120 srcField->checkConsistencyLight();
1122 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1123 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1124 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1125 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1126 if(srcField->getNature()!=targetField->getNature())
1127 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1128 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1130 std::ostringstream oss;
1131 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();
1132 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1133 throw INTERP_KERNEL::Exception(oss.str().c_str());
1135 DataArrayDouble *array(targetField->getArray());
1136 std::size_t srcNbOfCompo(srcField->getNumberOfComponents());
1139 targetField->checkConsistencyLight();
1140 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1141 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1146 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1147 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1148 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1149 targetField->setArray(tmp);
1151 computeDeno(srcField->getNature(),srcField,targetField);
1152 double *resPointer(targetField->getArray()->getPointer());
1153 const double *inputPointer(srcField->getArray()->getConstPointer());
1154 computeProduct(inputPointer,(int)srcNbOfCompo,isDftVal,dftValue,resPointer);
1157 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1160 return computeDenoFromScratch(nat,srcField,trgField);
1161 else if(nat!=_nature_of_deno)
1162 return computeDenoFromScratch(nat,srcField,trgField);
1163 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1164 return computeDenoFromScratch(nat,srcField,trgField);
1167 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1169 _nature_of_deno=nat;
1170 switch(_nature_of_deno)
1172 case IntensiveMaximum:
1174 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1177 case ExtensiveMaximum:
1179 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1180 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1181 const double *denoPtr=deno->getArray()->getConstPointer();
1182 const double *denoRPtr=denoR->getArray()->getConstPointer();
1183 if(trgField->getMesh()->getMeshDimension()==-1)
1185 double *denoRPtr2=denoR->getArray()->getPointer();
1186 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1188 if(srcField->getMesh()->getMeshDimension()==-1)
1190 double *denoPtr2=deno->getArray()->getPointer();
1191 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1194 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1195 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1197 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1198 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1204 case ExtensiveConservation:
1206 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1209 case IntensiveConservation:
1211 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1212 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1213 const double *denoPtr=deno->getArray()->getConstPointer();
1214 const double *denoRPtr=denoR->getArray()->getConstPointer();
1215 if(trgField->getMesh()->getMeshDimension()==-1)
1217 double *denoRPtr2=denoR->getArray()->getPointer();
1218 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1220 if(srcField->getMesh()->getMeshDimension()==-1)
1222 double *denoPtr2=deno->getArray()->getPointer();
1223 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1226 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1227 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1229 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1230 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1237 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1241 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1244 double *tmp=new double[inputNbOfCompo];
1245 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1247 if((*iter1).empty())
1250 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1254 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1255 std::map<mcIdType,double>::const_iterator iter3=_deno_multiply[idx].begin();
1256 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1258 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,(*iter2).second/(*iter3).second));
1259 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1265 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1267 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1269 double *tmp=new double[inputNbOfCompo];
1270 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1271 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1273 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1275 isReached[(*iter2).first]=true;
1276 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1277 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1282 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1284 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1287 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<mcIdType,double> >& matIn, mcIdType nbColsMatIn, std::vector<std::map<mcIdType,double> >& matOut)
1289 matOut.resize(nbColsMatIn);
1291 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1292 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1293 matOut[(*iter2).first][id]=(*iter2).second;
1296 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1297 std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1299 std::map<mcIdType,double> values;
1301 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1304 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1306 sum+=(*iter2).second;
1307 values[(*iter2).first]+=(*iter2).second;
1309 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1310 deno[idx][(*iter2).first]=sum;
1313 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1315 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1316 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1320 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1321 std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1323 std::map<mcIdType,double> values;
1325 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1328 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1330 sum+=(*iter2).second;
1331 values[(*iter2).first]+=(*iter2).second;
1333 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1334 denoReverse[(*iter2).first][idx]=sum;
1337 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1339 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1340 deno[idx][(*iter2).first]=values[(*iter2).first];
1344 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<mcIdType,double> >& m1D,
1345 const std::vector< std::map<mcIdType,double> >& m2D,
1346 const mcIdType *corrCellIdSrc, mcIdType nbOf2DCellsSrc, mcIdType nbOf1DCellsSrc,
1347 const mcIdType *corrCellIdTrg)
1349 mcIdType nbOf2DCellsTrg=ToIdType(m2D.size());
1350 mcIdType nbOf1DCellsTrg=ToIdType(m1D.size());
1351 mcIdType nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1352 _matrix.resize(nbOf3DCellsTrg);
1354 for(std::vector< std::map<mcIdType,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1356 for(std::map<mcIdType,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1359 for(std::vector< std::map<mcIdType,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1361 for(std::map<mcIdType,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1363 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1370 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<mcIdType,double> >& m)
1373 for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1375 std::cout << "Target Cell # " << id << " : ";
1376 for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1377 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1378 std::cout << std::endl;
1382 const std::vector<std::map<mcIdType,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1388 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1390 mcIdType MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1392 return ToIdType(_deno_reverse_multiply.size());
1396 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1397 * If not the behaviour is unpredictable.
1398 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1399 * 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.
1400 * 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
1403 * \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.
1404 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1405 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1407 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1410 std::vector<std::map<mcIdType,double> > matrixNew(_matrix.size());
1412 for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1414 std::map<mcIdType,double>& rowNew=matrixNew[i];
1415 for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1417 if(fabs((*it2).second)>maxValAbs)
1418 rowNew[(*it2).first]=(*it2).second;
1429 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1430 * If not the behaviour is unpredictable.
1431 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1432 * 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.
1433 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1434 * 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
1437 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1438 * \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
1439 * that all coefficients are null.
1440 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1442 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1444 double maxVal=getMaxValueInCrudeMatrix();
1447 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1451 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1452 * If not the behaviour is unpredictable.
1453 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1454 * The returned value is positive.
1456 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1459 for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1460 for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1461 if(fabs((*it2).second)>ret)
1462 ret=fabs((*it2).second);