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 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
542 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
543 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
544 INTERP_KERNEL::Interpolation3D interpolation(*this);
545 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
549 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
550 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
551 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
552 std::vector<std::map<int,double> > matrixTmp;
553 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
554 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
555 ReverseMatrix(matrixTmp,nbCols,_matrix);
556 nbCols=matrixTmp.size();
557 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
558 if(!duplicateFaces.empty())
560 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
561 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
563 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
564 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
570 else if(trgMeshDim==-1)
572 if(srcMeshDim==2 && srcSpaceDim==2)
574 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
575 INTERP_KERNEL::Interpolation2D interpolation(*this);
576 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
578 else if(srcMeshDim==3 && srcSpaceDim==3)
580 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
581 INTERP_KERNEL::Interpolation3D interpolation(*this);
582 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
584 else if(srcMeshDim==2 && srcSpaceDim==3)
586 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
587 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
588 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
591 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
593 else if(srcMeshDim==-1)
595 if(trgMeshDim==2 && trgSpaceDim==2)
597 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
598 INTERP_KERNEL::Interpolation2D interpolation(*this);
599 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
601 else if(trgMeshDim==3 && trgSpaceDim==3)
603 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
604 INTERP_KERNEL::Interpolation3D interpolation(*this);
605 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
607 else if(trgMeshDim==2 && trgSpaceDim==3)
609 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
610 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
611 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
614 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
617 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
618 _deno_multiply.clear();
619 _deno_multiply.resize(_matrix.size());
620 _deno_reverse_multiply.clear();
621 _deno_reverse_multiply.resize(nbCols);
626 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
628 std::string srcMeth,trgMeth;
629 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
630 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
631 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
633 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
634 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
635 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
636 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
637 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
638 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
639 std::vector<std::map<int,double> > matrix2D;
640 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
641 MEDCouplingUMesh *s1D,*t1D;
643 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
644 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
645 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
646 std::vector<std::map<int,double> > matrix1D;
647 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
648 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
649 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
650 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
653 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
654 target_mesh->getMesh3DIds()->getConstPointer());
656 _deno_multiply.clear();
657 _deno_multiply.resize(_matrix.size());
658 _deno_reverse_multiply.clear();
659 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
664 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
666 std::string srcMeth,trgMeth;
667 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
668 if(methodCpp!="P0P0")
669 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
670 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
671 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
672 const int srcMeshDim=src_mesh->getMeshDimension();
673 const int srcSpceDim=src_mesh->getSpaceDimension();
674 const int trgMeshDim=target_mesh->getMeshDimension();
675 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
676 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 !");
677 std::vector<std::map<int,double> > res;
682 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
683 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
684 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
685 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
690 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
691 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
692 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
693 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
698 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
699 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
700 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
701 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
705 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
707 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
708 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
710 _deno_multiply.clear();
711 _deno_multiply.resize(_matrix.size());
712 _deno_reverse_multiply.clear();
713 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
718 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
720 std::string srcMeth,trgMeth;
721 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
722 if(methodCpp!="P0P0")
723 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
724 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
725 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
726 const int srcMeshDim=src_mesh->getMeshDimension();
727 const int trgMeshDim=target_mesh->getMeshDimension();
728 const int trgSpceDim=target_mesh->getSpaceDimension();
729 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
730 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 !");
735 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
736 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
737 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
738 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
743 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
744 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
745 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
746 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
751 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
752 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
753 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
754 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
758 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
760 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
762 _deno_multiply.clear();
763 _deno_multiply.resize(_matrix.size());
764 _deno_reverse_multiply.clear();
765 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
770 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
772 std::string srcMeth,trgMeth;
773 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
774 if(methodCpp!="P0P0")
775 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
776 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
777 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
778 const int srcMeshDim=src_mesh->getMeshDimension();
779 const int trgMeshDim=target_mesh->getMeshDimension();
780 if(trgMeshDim!=srcMeshDim)
781 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
786 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
787 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
788 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
789 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
794 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
795 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
796 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
797 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
802 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
803 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
804 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
805 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
809 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
811 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
813 _deno_multiply.clear();
814 _deno_multiply.resize(_matrix.size());
815 _deno_reverse_multiply.clear();
816 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
821 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
823 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
824 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 !");
825 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
826 const double *trgLocPtr=trgLoc->begin();
827 int trgSpaceDim=trgLoc->getNumberOfComponents();
828 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
829 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
831 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
832 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
833 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
834 throw INTERP_KERNEL::Exception(oss.str().c_str());
836 const int *srcOffsetArrPtr=srcOffsetArr->begin();
837 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
838 const double *srcLocPtr=srcLoc->begin();
839 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
840 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
841 _matrix.resize(trgNbOfGaussPts);
842 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
843 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
844 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
845 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
846 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
848 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
849 int srcCellId=elts[eltsIndex[*trgId]];
850 double dist=std::numeric_limits<double>::max();
852 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
854 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
856 for(int i=0;i<trgSpaceDim;i++)
857 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
859 { dist=tmp; srcEntry=srcId; }
861 _matrix[*trgId][srcEntry]=1.;
863 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
865 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
866 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
867 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
868 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
869 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
870 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
872 _deno_multiply.clear();
873 _deno_multiply.resize(_matrix.size());
874 _deno_reverse_multiply.clear();
875 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
881 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
882 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
884 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
886 if(method=="GAUSSGAUSS")
888 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
889 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
890 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
891 throw INTERP_KERNEL::Exception(oss.str().c_str());
895 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
896 * 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
897 * only INTERP_KERNEL method should be applied.
899 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
901 std::string srcm,trgm,method;
902 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
903 switch(_interp_matrix_pol)
905 case IK_ONLY_PREFERED:
909 std::string tmp1,tmp2;
910 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
913 catch(INTERP_KERNEL::Exception& /*e*/)
918 case NOT_IK_ONLY_PREFERED:
922 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
925 catch(INTERP_KERNEL::Exception& /*e*/)
932 case NOT_IK_ONLY_FORCED:
935 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
939 void MEDCouplingRemapper::updateTime() const
943 void MEDCouplingRemapper::checkPrepare() const
945 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
947 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
948 if(!s->getMesh() || !t->getMesh())
949 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
953 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
954 * This method returns 3 informations (2 in ouput parameters and 1 in return).
956 * \param [out] srcMeth the string code of the discretization of source field template
957 * \param [out] trgMeth the string code of the discretization of target field template
958 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
960 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
962 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
964 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
965 if(!s->getMesh() || !t->getMesh())
966 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
967 srcMeth=_src_ft->getDiscretization()->getRepr();
968 trgMeth=_target_ft->getDiscretization()->getRepr();
969 return BuildMethodFrom(srcMeth,trgMeth);
972 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
974 std::string method(meth1); method+=meth2;
978 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
982 if(matrixSuppression)
985 _deno_multiply.clear();
986 _deno_reverse_multiply.clear();
990 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
992 if(!srcField || !targetField)
993 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
994 srcField->checkConsistencyLight();
996 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
997 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
998 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
999 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1000 if(srcField->getNature()!=targetField->getNature())
1001 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1002 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1004 std::ostringstream oss;
1005 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();
1006 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1007 throw INTERP_KERNEL::Exception(oss.str().c_str());
1009 DataArrayDouble *array(targetField->getArray());
1010 int srcNbOfCompo(srcField->getNumberOfComponents());
1013 targetField->checkConsistencyLight();
1014 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1015 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1020 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1021 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1022 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1023 targetField->setArray(tmp);
1025 computeDeno(srcField->getNature(),srcField,targetField);
1026 double *resPointer(targetField->getArray()->getPointer());
1027 const double *inputPointer(srcField->getArray()->getConstPointer());
1028 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1031 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1034 return computeDenoFromScratch(nat,srcField,trgField);
1035 else if(nat!=_nature_of_deno)
1036 return computeDenoFromScratch(nat,srcField,trgField);
1037 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1038 return computeDenoFromScratch(nat,srcField,trgField);
1041 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1043 _nature_of_deno=nat;
1044 _time_deno_update=getTimeOfThis();
1045 switch(_nature_of_deno)
1047 case IntensiveMaximum:
1049 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1052 case ExtensiveMaximum:
1054 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1055 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1056 const double *denoPtr=deno->getArray()->getConstPointer();
1057 const double *denoRPtr=denoR->getArray()->getConstPointer();
1058 if(trgField->getMesh()->getMeshDimension()==-1)
1060 double *denoRPtr2=denoR->getArray()->getPointer();
1061 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1063 if(srcField->getMesh()->getMeshDimension()==-1)
1065 double *denoPtr2=deno->getArray()->getPointer();
1066 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1069 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1070 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1072 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1073 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1079 case ExtensiveConservation:
1081 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1084 case IntensiveConservation:
1086 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1087 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1088 const double *denoPtr=deno->getArray()->getConstPointer();
1089 const double *denoRPtr=denoR->getArray()->getConstPointer();
1090 if(trgField->getMesh()->getMeshDimension()==-1)
1092 double *denoRPtr2=denoR->getArray()->getPointer();
1093 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1095 if(srcField->getMesh()->getMeshDimension()==-1)
1097 double *denoPtr2=deno->getArray()->getPointer();
1098 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1101 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1102 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1104 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1105 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1112 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1116 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1119 double *tmp=new double[inputNbOfCompo];
1120 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1122 if((*iter1).empty())
1125 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1129 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1130 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1131 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1133 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1134 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1140 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1142 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1144 double *tmp=new double[inputNbOfCompo];
1145 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1146 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1148 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1150 isReached[(*iter2).first]=true;
1151 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1152 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1157 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1159 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1162 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1164 matOut.resize(nbColsMatIn);
1166 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1167 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1168 matOut[(*iter2).first][id]=(*iter2).second;
1171 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1172 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1174 std::map<int,double> values;
1176 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1179 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1181 sum+=(*iter2).second;
1182 values[(*iter2).first]+=(*iter2).second;
1184 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1185 deno[idx][(*iter2).first]=sum;
1188 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1190 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1191 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1195 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1196 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1198 std::map<int,double> values;
1200 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1203 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1205 sum+=(*iter2).second;
1206 values[(*iter2).first]+=(*iter2).second;
1208 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1209 denoReverse[(*iter2).first][idx]=sum;
1212 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1214 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1215 deno[idx][(*iter2).first]=values[(*iter2).first];
1219 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1220 const std::vector< std::map<int,double> >& m2D,
1221 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1222 const int *corrCellIdTrg)
1224 int nbOf2DCellsTrg=m2D.size();
1225 int nbOf1DCellsTrg=m1D.size();
1226 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1227 _matrix.resize(nbOf3DCellsTrg);
1229 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1231 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1234 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1236 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1238 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1245 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1248 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1250 std::cout << "Target Cell # " << id << " : ";
1251 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1252 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1253 std::cout << std::endl;
1257 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1263 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1265 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1267 return (int)_deno_reverse_multiply.size();
1271 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1272 * If not the behaviour is unpredictable.
1273 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1274 * 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.
1275 * 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
1278 * \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.
1279 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1280 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1282 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1285 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1287 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1289 std::map<int,double>& rowNew=matrixNew[i];
1290 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1292 if(fabs((*it2).second)>maxValAbs)
1293 rowNew[(*it2).first]=(*it2).second;
1304 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1305 * If not the behaviour is unpredictable.
1306 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1307 * 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.
1308 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1309 * 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
1312 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1313 * \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
1314 * that all coefficients are null.
1315 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1317 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1319 double maxVal=getMaxValueInCrudeMatrix();
1322 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1326 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1327 * If not the behaviour is unpredictable.
1328 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1329 * The returned value is positive.
1331 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1334 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1335 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1336 if(fabs((*it2).second)>ret)
1337 ret=fabs((*it2).second);