1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingRemapper.hxx"
22 #include "MEDCouplingMemArray.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24 #include "MEDCouplingFieldTemplate.hxx"
25 #include "MEDCouplingFieldDiscretization.hxx"
26 #include "MEDCouplingMappedExtrudedMesh.hxx"
27 #include "MEDCouplingCMesh.hxx"
28 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
29 #include "MEDCouplingNormalizedCartesianMesh.txx"
31 #include "Interpolation1D.txx"
32 #include "Interpolation2DCurve.hxx"
33 #include "Interpolation2D.txx"
34 #include "Interpolation3D.txx"
35 #include "Interpolation3DSurf.hxx"
36 #include "Interpolation2D1D.txx"
37 #include "Interpolation2D3D.txx"
38 #include "Interpolation3D1D.txx"
39 #include "InterpolationCU.txx"
40 #include "InterpolationCC.txx"
42 using namespace MEDCoupling;
44 MEDCouplingRemapper::MEDCouplingRemapper():_src_ft(0),_target_ft(0),_interp_matrix_pol(IK_ONLY_PREFERED),_nature_of_deno(NoNature),_time_deno_update(0)
48 MEDCouplingRemapper::~MEDCouplingRemapper()
54 * This method is the second step of the remapping process. The remapping process works in three phases :
56 * - Set remapping options appropriately
57 * - The computation of remapping matrix
58 * - Apply the matrix vector multiply to obtain the result of the remapping
60 * 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.
61 * 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).
62 * So this method is less precise but a more user friendly way to compute a remapping matrix.
64 * \param [in] srcMesh the source mesh
65 * \param [in] targetMesh the target mesh
66 * \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.
67 * 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).
69 * \sa MEDCouplingRemapper::prepareEx
71 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
73 if(!srcMesh || !targetMesh)
74 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepare : presence of NULL input pointer !");
75 std::string srcMethod,targetMethod;
76 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
77 MCAuto<MEDCouplingFieldTemplate> src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
78 src->setMesh(srcMesh);
79 MCAuto<MEDCouplingFieldTemplate> target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
80 target->setMesh(targetMesh);
81 return prepareEx(src,target);
85 * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
86 * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
88 * \param [in] src is the field template source side.
89 * \param [in] target is the field template target side.
91 * \sa MEDCouplingRemapper::prepare
93 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
96 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL input pointer !");
97 if(!src->getMesh() || !target->getMesh())
98 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
100 _src_ft=const_cast<MEDCouplingFieldTemplate *>(src); _src_ft->incrRef();
101 _target_ft=const_cast<MEDCouplingFieldTemplate *>(target); _target_ft->incrRef();
102 if(isInterpKernelOnlyOrNotOnly())
103 return prepareInterpKernelOnly();
105 return prepareNotInterpKernelOnly();
108 int MEDCouplingRemapper::prepareInterpKernelOnly()
110 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
118 // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
119 // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
121 // } MEDCouplingMeshType;
123 switch(meshInterpType)
125 case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
126 case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
127 case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
128 case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
129 case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
130 case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
131 case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
132 case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
133 case 85: // UNSTRUCTURED - UNSTRUCTURED
134 return prepareInterpKernelOnlyUU();
135 case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
136 case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
137 case 87: // UNSTRUCTURED - CARTESIAN
138 return prepareInterpKernelOnlyUC();
139 case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
140 case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
141 case 117: // CARTESIAN - UNSTRUCTURED
142 return prepareInterpKernelOnlyCU();
143 case 119: // CARTESIAN - CARTESIAN
144 return prepareInterpKernelOnlyCC();
145 case 136: // EXTRUDED - EXTRUDED
146 return prepareInterpKernelOnlyEE();
148 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
152 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
154 std::string srcm,trgm,method;
155 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
156 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
159 return prepareNotInterpKernelOnlyGaussGauss();
162 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
163 throw INTERP_KERNEL::Exception(oss.str().c_str());
169 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
170 * 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.
172 * \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.
173 * \param [in/out] targetField the destination field with the allocated array in which all tuples will be overwritten.
174 * \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.
175 * 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
176 * 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
177 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
181 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue)
183 if(!srcField || !targetField)
184 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transfer : input field must be both not NULL !");
185 transferUnderground(srcField,targetField,true,dftValue);
189 * This method is equivalent to MEDCoupling::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
190 * 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
192 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
194 * \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.
195 * \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.
197 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField)
199 if(!srcField || !targetField)
200 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : input field must be both not NULL !");
201 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
204 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue)
206 if(!srcField || !targetField)
207 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::reverseTransfer : input fields must be both not NULL !");
209 targetField->checkConsistencyLight();
210 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
211 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
212 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
213 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
214 if(srcField->getNature()!=targetField->getNature())
215 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
216 if(targetField->getNumberOfTuplesExpected()!=_target_ft->getNumberOfTuplesExpected())
218 std::ostringstream oss;
219 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();
220 oss << " ! It appears that the target support is not the same between the prepare and the transfer !";
221 throw INTERP_KERNEL::Exception(oss.str().c_str());
223 DataArrayDouble *array(srcField->getArray());
224 int trgNbOfCompo=targetField->getNumberOfComponents();
227 srcField->checkConsistencyLight();
228 if(trgNbOfCompo!=srcField->getNumberOfTuplesExpected())
229 throw INTERP_KERNEL::Exception("Number of components mismatch !");
233 MCAuto<DataArrayDouble > tmp(DataArrayDouble::New());
234 tmp->alloc(srcField->getNumberOfTuplesExpected(),trgNbOfCompo);
235 srcField->setArray(tmp);
237 computeDeno(srcField->getNature(),srcField,targetField);
238 double *resPointer(srcField->getArray()->getPointer());
239 const double *inputPointer=targetField->getArray()->getConstPointer();
240 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
244 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
245 * If mesh of \b srcField does not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
247 * \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.
248 * \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.
249 * 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
250 * 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
251 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
252 * \return destination field to be deallocated by the caller.
256 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue)
260 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input srcField is NULL !");
261 srcField->checkConsistencyLight();
262 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
263 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
264 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
265 ret->setNature(srcField->getNature());
266 transfer(srcField,ret,dftValue);
267 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
271 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue)
274 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input targetField is NULL !");
275 targetField->checkConsistencyLight();
277 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
278 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
279 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
280 ret->setNature(targetField->getNature());
281 reverseTransfer(ret,targetField,dftValue);
282 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
287 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
288 * is here only for automatic CORBA generators.
290 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
292 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
296 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
297 * is here only for automatic CORBA generators.
299 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
301 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
305 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
306 * is here only for automatic CORBA generators.
308 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
310 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
314 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or prefered.
315 * If interpolation matrix policy is :
317 * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is prefered. That is to say, if it is possible to treat the case
318 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
319 * If not, the \b not only INTERP_KERNEL method will be attempt.
321 * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is prefered. That is to say, if it is possible to treat the case
322 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
323 * If not, the INTERP_KERNEL only method will be attempt.
325 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
327 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
329 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
331 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
333 return _interp_matrix_pol;
337 * 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
338 * 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
339 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
341 * If interpolation matrix policy is :
343 * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is prefered. That is to say, if it is possible to treat the case
344 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
345 * If not, the \b not only INTERP_KERNEL method will be attempt.
347 * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is prefered. That is to say, if it is possible to treat the case
348 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
349 * If not, the INTERP_KERNEL only method will be attempt.
351 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
353 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
355 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c MEDCoupling::InterpolationMatrixPolicy
356 * for automatic generation of CORBA component.
358 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
360 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
362 switch(newInterpMatPol)
365 _interp_matrix_pol=IK_ONLY_PREFERED;
368 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
371 _interp_matrix_pol=IK_ONLY_FORCED;
374 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
377 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 !");
381 int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
383 const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
384 const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_target_ft->getMesh());
385 std::string srcMeth,trgMeth;
386 std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth));
387 const int srcMeshDim=src_mesh->getMeshDimension();
390 srcSpaceDim=src_mesh->getSpaceDimension();
391 const int trgMeshDim=target_mesh->getMeshDimension();
394 trgSpaceDim=target_mesh->getSpaceDimension();
395 if(trgSpaceDim!=srcSpaceDim)
396 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
397 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
399 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
401 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
402 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
403 INTERP_KERNEL::Interpolation1D interpolation(*this);
404 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
406 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
408 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
409 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
410 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
411 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
413 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
415 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
416 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
417 INTERP_KERNEL::Interpolation2D interpolation(*this);
418 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
420 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
422 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
423 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
424 INTERP_KERNEL::Interpolation3D interpolation(*this);
425 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
427 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
429 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
430 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
431 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
432 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
434 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
436 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
437 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
438 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
439 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
440 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
441 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
443 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
445 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
446 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
447 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
448 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
449 INTERP_KERNEL::Interpolation3D1D interpolation(*this);
450 std::vector<std::map<int,double> > matrixTmp;
451 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
452 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
453 ReverseMatrix(matrixTmp,nbCols,_matrix);
454 nbCols=matrixTmp.size();
456 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
458 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
460 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
461 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
462 INTERP_KERNEL::Interpolation2D interpolation(*this);
463 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
467 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
468 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
469 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
470 std::vector<std::map<int,double> > matrixTmp;
471 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
472 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
473 ReverseMatrix(matrixTmp,nbCols,_matrix);
474 nbCols=matrixTmp.size();
475 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
476 if(!duplicateFaces.empty())
478 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
479 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
481 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
482 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
488 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
490 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
492 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
493 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
494 INTERP_KERNEL::Interpolation2D interpolation(*this);
495 std::vector<std::map<int,double> > matrixTmp;
496 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
497 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
498 ReverseMatrix(matrixTmp,nbCols,_matrix);
499 nbCols=matrixTmp.size();
503 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
504 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
505 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
506 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
507 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
508 if(!duplicateFaces.empty())
510 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
511 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
513 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
514 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
520 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
522 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
523 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
524 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
525 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
526 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
527 if(!duplicateFaces.empty())
529 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
530 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
532 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
533 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
538 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
540 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
541 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
542 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
543 std::vector<std::map<int,double> > matrixTmp;
544 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
545 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
546 ReverseMatrix(matrixTmp,nbCols,_matrix);
547 nbCols=matrixTmp.size();
548 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
549 if(!duplicateFaces.empty())
551 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
552 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
554 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
555 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
560 else if(trgMeshDim==-1)
562 if(srcMeshDim==2 && srcSpaceDim==2)
564 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
565 INTERP_KERNEL::Interpolation2D interpolation(*this);
566 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
568 else if(srcMeshDim==3 && srcSpaceDim==3)
570 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
571 INTERP_KERNEL::Interpolation3D interpolation(*this);
572 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
574 else if(srcMeshDim==2 && srcSpaceDim==3)
576 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
577 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
578 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
581 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
583 else if(srcMeshDim==-1)
585 if(trgMeshDim==2 && trgSpaceDim==2)
587 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
588 INTERP_KERNEL::Interpolation2D interpolation(*this);
589 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
591 else if(trgMeshDim==3 && trgSpaceDim==3)
593 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
594 INTERP_KERNEL::Interpolation3D interpolation(*this);
595 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
597 else if(trgMeshDim==2 && trgSpaceDim==3)
599 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
600 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
601 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
604 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
607 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
608 _deno_multiply.clear();
609 _deno_multiply.resize(_matrix.size());
610 _deno_reverse_multiply.clear();
611 _deno_reverse_multiply.resize(nbCols);
616 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
618 std::string srcMeth,trgMeth;
619 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
620 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
621 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
623 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
624 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
625 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
626 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
627 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
628 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
629 std::vector<std::map<int,double> > matrix2D;
630 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
631 MEDCouplingUMesh *s1D,*t1D;
633 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
634 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
635 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
636 std::vector<std::map<int,double> > matrix1D;
637 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
638 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
639 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
640 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
643 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
644 target_mesh->getMesh3DIds()->getConstPointer());
646 _deno_multiply.clear();
647 _deno_multiply.resize(_matrix.size());
648 _deno_reverse_multiply.clear();
649 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
654 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
656 std::string srcMeth,trgMeth;
657 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
658 if(methodCpp!="P0P0")
659 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
660 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
661 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
662 const int srcMeshDim=src_mesh->getMeshDimension();
663 const int srcSpceDim=src_mesh->getSpaceDimension();
664 const int trgMeshDim=target_mesh->getMeshDimension();
665 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
666 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : space dim of src unstructured should be equal to mesh dim of src unstructured and should be equal also equal to trg cartesian dimension !");
667 std::vector<std::map<int,double> > res;
672 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
673 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
674 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
675 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
680 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
681 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
682 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
683 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
688 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
689 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
690 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
691 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
695 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
697 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
698 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
700 _deno_multiply.clear();
701 _deno_multiply.resize(_matrix.size());
702 _deno_reverse_multiply.clear();
703 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
708 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
710 std::string srcMeth,trgMeth;
711 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
712 if(methodCpp!="P0P0")
713 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
714 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
715 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
716 const int srcMeshDim=src_mesh->getMeshDimension();
717 const int trgMeshDim=target_mesh->getMeshDimension();
718 const int trgSpceDim=target_mesh->getSpaceDimension();
719 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
720 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : space dim of target unstructured should be equal to mesh dim of target unstructured and should be equal also equal to source cartesian dimension !");
725 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
726 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
727 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
728 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
733 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
734 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
735 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
736 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
741 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
742 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
743 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
744 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
748 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
750 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
752 _deno_multiply.clear();
753 _deno_multiply.resize(_matrix.size());
754 _deno_reverse_multiply.clear();
755 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
760 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
762 std::string srcMeth,trgMeth;
763 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
764 if(methodCpp!="P0P0")
765 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
766 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
767 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
768 const int srcMeshDim=src_mesh->getMeshDimension();
769 const int trgMeshDim=target_mesh->getMeshDimension();
770 if(trgMeshDim!=srcMeshDim)
771 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
776 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
777 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
778 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
779 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
784 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
785 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
786 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
787 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
792 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
793 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
794 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
795 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
799 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
801 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
803 _deno_multiply.clear();
804 _deno_multiply.resize(_matrix.size());
805 _deno_reverse_multiply.clear();
806 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
811 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
813 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
814 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 !");
815 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
816 const double *trgLocPtr=trgLoc->begin();
817 int trgSpaceDim=trgLoc->getNumberOfComponents();
818 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
819 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
821 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
822 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
823 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
824 throw INTERP_KERNEL::Exception(oss.str().c_str());
826 const int *srcOffsetArrPtr=srcOffsetArr->begin();
827 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
828 const double *srcLocPtr=srcLoc->begin();
829 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
830 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
831 _matrix.resize(trgNbOfGaussPts);
832 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
833 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
834 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
835 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
836 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
838 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
839 int srcCellId=elts[eltsIndex[*trgId]];
840 double dist=std::numeric_limits<double>::max();
842 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
844 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
846 for(int i=0;i<trgSpaceDim;i++)
847 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
849 { dist=tmp; srcEntry=srcId; }
851 _matrix[*trgId][srcEntry]=1.;
853 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
855 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
856 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
857 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
858 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
859 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
860 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
862 _deno_multiply.clear();
863 _deno_multiply.resize(_matrix.size());
864 _deno_reverse_multiply.clear();
865 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
871 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
872 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
874 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
876 if(method=="GAUSSGAUSS")
878 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
879 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
880 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
881 throw INTERP_KERNEL::Exception(oss.str().c_str());
885 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
886 * 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
887 * only INTERP_KERNEL method should be applied.
889 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
891 std::string srcm,trgm,method;
892 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
893 switch(_interp_matrix_pol)
895 case IK_ONLY_PREFERED:
899 std::string tmp1,tmp2;
900 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
903 catch(INTERP_KERNEL::Exception& /*e*/)
908 case NOT_IK_ONLY_PREFERED:
912 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
915 catch(INTERP_KERNEL::Exception& /*e*/)
922 case NOT_IK_ONLY_FORCED:
925 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
929 void MEDCouplingRemapper::updateTime() const
933 void MEDCouplingRemapper::checkPrepare() const
935 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
937 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
938 if(!s->getMesh() || !t->getMesh())
939 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
943 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
944 * This method returns 3 informations (2 in ouput parameters and 1 in return).
946 * \param [out] srcMeth the string code of the discretization of source field template
947 * \param [out] trgMeth the string code of the discretization of target field template
948 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
950 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
952 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
954 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
955 if(!s->getMesh() || !t->getMesh())
956 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
957 srcMeth=_src_ft->getDiscretization()->getRepr();
958 trgMeth=_target_ft->getDiscretization()->getRepr();
959 return BuildMethodFrom(srcMeth,trgMeth);
962 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
964 std::string method(meth1); method+=meth2;
968 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
972 if(matrixSuppression)
975 _deno_multiply.clear();
976 _deno_reverse_multiply.clear();
980 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
982 if(!srcField || !targetField)
983 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
984 srcField->checkConsistencyLight();
986 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
987 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
988 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
989 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
990 if(srcField->getNature()!=targetField->getNature())
991 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
992 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
994 std::ostringstream oss;
995 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();
996 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
997 throw INTERP_KERNEL::Exception(oss.str().c_str());
999 DataArrayDouble *array(targetField->getArray());
1000 int srcNbOfCompo(srcField->getNumberOfComponents());
1003 targetField->checkConsistencyLight();
1004 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1005 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1010 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1011 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1012 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1013 targetField->setArray(tmp);
1015 computeDeno(srcField->getNature(),srcField,targetField);
1016 double *resPointer(targetField->getArray()->getPointer());
1017 const double *inputPointer(srcField->getArray()->getConstPointer());
1018 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1021 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1024 return computeDenoFromScratch(nat,srcField,trgField);
1025 else if(nat!=_nature_of_deno)
1026 return computeDenoFromScratch(nat,srcField,trgField);
1027 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1028 return computeDenoFromScratch(nat,srcField,trgField);
1031 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1033 _nature_of_deno=nat;
1034 _time_deno_update=getTimeOfThis();
1035 switch(_nature_of_deno)
1037 case IntensiveMaximum:
1039 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1042 case ExtensiveMaximum:
1044 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1045 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1046 const double *denoPtr=deno->getArray()->getConstPointer();
1047 const double *denoRPtr=denoR->getArray()->getConstPointer();
1048 if(trgField->getMesh()->getMeshDimension()==-1)
1050 double *denoRPtr2=denoR->getArray()->getPointer();
1051 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1053 if(srcField->getMesh()->getMeshDimension()==-1)
1055 double *denoPtr2=deno->getArray()->getPointer();
1056 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1059 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1060 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1062 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1063 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1069 case ExtensiveConservation:
1071 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1074 case IntensiveConservation:
1076 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1077 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1078 const double *denoPtr=deno->getArray()->getConstPointer();
1079 const double *denoRPtr=denoR->getArray()->getConstPointer();
1080 if(trgField->getMesh()->getMeshDimension()==-1)
1082 double *denoRPtr2=denoR->getArray()->getPointer();
1083 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1085 if(srcField->getMesh()->getMeshDimension()==-1)
1087 double *denoPtr2=deno->getArray()->getPointer();
1088 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1091 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1092 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1094 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1095 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1102 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1106 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1109 double *tmp=new double[inputNbOfCompo];
1110 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1112 if((*iter1).empty())
1115 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1119 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1120 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1121 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1123 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1124 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1130 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1132 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1134 double *tmp=new double[inputNbOfCompo];
1135 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1136 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1138 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1140 isReached[(*iter2).first]=true;
1141 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1142 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1147 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1149 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1152 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1154 matOut.resize(nbColsMatIn);
1156 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1157 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1158 matOut[(*iter2).first][id]=(*iter2).second;
1161 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1162 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1164 std::map<int,double> values;
1166 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1169 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1171 sum+=(*iter2).second;
1172 values[(*iter2).first]+=(*iter2).second;
1174 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1175 deno[idx][(*iter2).first]=sum;
1178 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1180 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1181 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1185 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1186 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1188 std::map<int,double> values;
1190 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1193 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1195 sum+=(*iter2).second;
1196 values[(*iter2).first]+=(*iter2).second;
1198 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1199 denoReverse[(*iter2).first][idx]=sum;
1202 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1204 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1205 deno[idx][(*iter2).first]=values[(*iter2).first];
1209 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1210 const std::vector< std::map<int,double> >& m2D,
1211 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1212 const int *corrCellIdTrg)
1214 int nbOf2DCellsTrg=m2D.size();
1215 int nbOf1DCellsTrg=m1D.size();
1216 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1217 _matrix.resize(nbOf3DCellsTrg);
1219 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1221 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1224 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1226 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1228 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1235 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1238 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1240 std::cout << "Target Cell # " << id << " : ";
1241 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1242 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1243 std::cout << std::endl;
1247 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1253 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1255 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1257 return (int)_deno_reverse_multiply.size();
1261 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1262 * If not the behaviour is unpredictable.
1263 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1264 * 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.
1265 * 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
1268 * \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.
1269 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1270 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1272 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1275 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1277 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1279 std::map<int,double>& rowNew=matrixNew[i];
1280 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1282 if(fabs((*it2).second)>maxValAbs)
1283 rowNew[(*it2).first]=(*it2).second;
1294 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1295 * If not the behaviour is unpredictable.
1296 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1297 * 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.
1298 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1299 * 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
1302 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1303 * \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
1304 * that all coefficients are null.
1305 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1307 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1309 double maxVal=getMaxValueInCrudeMatrix();
1312 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1316 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1317 * If not the behaviour is unpredictable.
1318 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1319 * The returned value is positive.
1321 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1324 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1325 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1326 if(fabs((*it2).second)>ret)
1327 ret=fabs((*it2).second);